# Analytical detection of stationary and dynamic patterns in a prey-predator
model with reproductive Allee effect in prey growth

Subrata Dey, S Ghorai, Malay Banerjee Corresponding author

###### Abstract

Allee effect in population dynamics has a major impact in suppressing the
paradox of enrichment through global bifurcation, and it can generate highly
complex dynamics. The influence of the reproductive Allee effect, incorporated
in the prey’s growth rate of a prey-predator model with Beddington-DeAngelis
functional response, is investigated here. Preliminary local and global
bifurcations are identified of the temporal model. Existence and non-existence
of heterogeneous steady-state solutions of the spatio-temporal system are
established for suitable ranges of parameter values. The spatio-temporal model
satisfies Turing instability conditions, but numerical investigation reveals
that the heterogeneous patterns corresponding to unstable Turing eigen modes
acts as a transitory pattern. Inclusion of the reproductive Allee effect in
the prey population has a destabilising effect on the coexistence equilibrium.
For a range of parameter values, various branches of stationary solutions
including mode-dependent Turing solutions and localized pattern solutions are
identified using numerical bifurcation technique. The model is also capable to
produce some complex dynamic patterns such as travelling wave, moving pulse
solution, and spatio-temporal chaos for certain range of parameters and
diffusivity along with appropriate choice of initial conditions Judicious
choices of parametrization for the Beddington-DeAngelis functional response
help us to infer about the resulting patterns for similar prey-predator models
with Holling type-II functional response and ratio-dependent functional
response.

Department of Mathematics and Statistics, Indian Institute of Technology
Kanpur, Kanpur 208016, India

E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

##  1 Introduction

Understanding self-organized spatio-temporal patterns and mechanisms of
spatial dispersal in interacting populations have always been important areas
of mathematical biology and theoretical ecology. Since the pioneering work of
Lotka [1] and Volterra [2], a wide variety of mathematical models have been
developed to investigate various complex ecological phenomena, e.g., Allee
effect [3], group defense [4], parental care [5], hunting cooperation [6],
intra-guild predation [7], etc. In addition, persistence, stability of steady-
state, local and global bifurcation, etc., of these models play significant
roles in dynamical system theory. The study of such phenomena under spatial
distribution makes it more interesting and realistic in spatial ecology. The
spatio-temporal distribution of species has been extensively studied to
explain pattern formation observed in fish skin [8], prey-predator
interactions [9], terrestrial vegetation [10], mussel bed [11], etc.

The classical explanation of spatial pattern formation comes from Turing’s
pioneering work in reaction-diffusion (RD) system for chemical morphogenesis
[12]. The framework of the diffusion-driven instability in Turing’s work needs
two basic criteria: an activator-inhibitor interaction structure of at least
two interacting species, and an order of magnitude difference in their
dispersal coefficient. Segel and Jackson [13] first used Turing’s model
framework to explain the spatial pattern formation for ecological systems.
Turing instability produces various stationary patterns such as spots,
stripes, and a mixture of both [14, 15]. For one-dimensional habitats, Turing
patterns lead to stationary spatially periodic patterns [16, 17]. As
parameters are varied towards certain thresholds, the Turing pattern solution
becomes a localized stationary pattern solution due to another instability,
the so-called Belyakov-Devaney (BD) transition [18, 19]. The bounds and the
region of existence of the localized pattern solution, and the homoclinic
snaking mechanism [20, 21] of the patterns can be described by semi-strong
asymptotic analysis [22]. These spatial instabilities depend on the reaction
kinetics of the system. For example, Rosenzweig-MacArthur(RM) prey-predator
system, one of the most popular models in ecology, does not follow Turing’s
activator-inhibitor structure and it does not produce any stationary patterns
[23, 24]. Some additional nonlinear effects such as predator-dependent
functional response, generalist type predator, nonlinear death rate of the
predator, Allee effect, etc., [15, 25, 17, 26] are required to produce a
stationary Turing pattern.

The Allee effect, introduced by renowned ecologist W. C. Allee [27], has
important consequences on population dynamics and their persistence. It refers
to the correlation between the population density and per capita growth rate
of the population at low densities [3]. It plays an important role in
ecological conservation and wildlife management since it is closely related to
population extinction. Generally speaking, there are two kinds of Allee
effects [28]: component Allee effect and demographic Allee effect. Component
Allee effect focuses on components of individual fitness, whereas demographic
Allee effect is at the level of the overall mean individual fitness of the
whole population. Allee effect can be induced by a variety of factors that
include difficulties in social interaction, cooperative breeding, reproductive
facilitation, anti-predator behavior, mate finding, and environmental
conditioning, among others [28, 29, 30]. The demographic Allee effect is a
consequence of one or more components Allee effect. The reproductive Allee
effect is due to many different mechanisms. Due to the difficulty in finding a
suitable mate during their receptive period, individuals in small or sparse
populations face the mate-finding Allee effect and have less reproductive
success. Mate-finding Allee effect is observed in many species such as
Glanville fritillary butterfly [31], sheep ticks [32], polar bears [33], queen
conch [34], and many others [see [28, 35] for further examples]. Other
mechanisms of the reproductive Allee effect include broadcast spawning in
corals, pollen limitation in field gentian, reproductive facilitation in
whiptail lizard, and sperm limitation in blue crab, among others [28, 35].
This component or demographic Allee effect can be strong or weak. If there is
a threshold population density, called the Allee threshold, below which the
per capita growth rate of the population is negative, then it is called the
strong Allee effect. For the weak Allee effect, there is no such Allee
threshold. If the population density rises above this threshold, it will
persist, whereas it will become extinct if it falls below [36, 28].

Let $N(T)$ be the density of a single species population (prey) at time $T$.
The single species logistic growth is described by the equation
$\displaystyle{\frac{dN}{dT}=rN\left(1-\frac{N}{K}\right)}$; $r$ and $K$
respectively denote the intrinsic growth rate and carrying capacity of the
population. This description has a limitation in that the per capita growth
rate is maximum at low population density. Single species growth with the
Allee effect assumes that the per capita growth rate is an increasing function
of the population at low population density. There exist several
parametrizations of the growth of a single species population subject to the
Allee effect. In terms of mathematical formulations, those can be broadly
classified into two types: additive Allee effect and multiplicative Allee
effect. In the case of additive Allee effect, the logistic growth equation is
modified as [3, 37]

$\frac{dN}{dT}\,=\,rN\left(1-\frac{N}{K}-\frac{s}{N+\phi}\right),$

where $s$ and $\phi$ indicate the severity of the additive Allee effect. On
the other hand, some well-known mathematical formulations for the
multiplicative Allee effect include the parametrization
$\displaystyle{rN\left(1-\frac{N}{K}\right)\left(\frac{N}{K}-\frac{\theta}{K}\right)}$
[38], $\displaystyle{rN\left(1-\frac{N}{K}\right)\frac{N}{N+\theta}}$ [39],
and
$\displaystyle{rN\left(1-\frac{N}{K}\right)\left(1-\frac{N+C}{N+\theta}\right)}$
[29, 28]. In these parametrizations, $\theta$ is the strength of the Allee
effect, and $C$ characterizes the shape of the per capita growth rate curve.
The most widely used form of the multiple Allee effect is

$\frac{dN}{dT}=rN\left(1-\frac{N}{K}\right)\left(\frac{N}{K}-\frac{\theta}{K}\right),$
(1)

where $0<\theta<K$ corresponds to strong Allee effect and $-K<\theta<0$ stands
for weak Allee effect. This formulation has a drawback as the intra-specific
cooperation and competition appear jointly in the formulation in a
multiplicative way. To overcome this situation, Petrovskii et al. [40]
proposed a different parametrization to model the single species population
growth subject to Allee effect as follows:

$\frac{dN}{dT}=N\big{(}B(N)-c-I(N)\big{)},$ (2)

where $B(N)$ and $I(N)$ represent population multiplication due to
reproduction and an additional density-dependent mortality rate due to intra-
species competition respectively and $c$ denotes the intrinsic mortality rate.
Taking the parametrization of sexual reproduction of the population and
density-dependent mortality rate in the form $B(N)=abN$ and $I(N)=aN^{2}$
(following the approach in [41, 42]), we can write the single species growth
equation as

$\frac{dN}{dT}=aN^{2}(b-N)-cN,$ (3)

where the positive constants $a$ and $b$ can be interpreted as intrinsic
growth rate and threshold for positive growth rate respectively. This
description can be considered as the growth equation of a two-sex population
model [43]. To illustrate this claim, let $N_{m}$ and $N_{f}$ denote the
densities of male and female populations of a species ($N_{m}+N_{f}=N$), and
let $R$ be the available resources. Then the reproduction rate is proportional
to $N_{m}N_{f}R,$ and the available resources $R$ can be expressed as
$R=K-S(N_{m}+N_{f})$, where $K$ is the carrying capacity and $S$ denotes the
per capita resource consumption rate. Assuming an equal sex ratio, i.e.,
$N_{m}=N_{f}=N/2$, we find the growth rate due to reproduction to be
proportional to $N^{2}(K-SN)/4$. Therefore, the growth equation can be written
as

$\frac{dN}{dT}=R_{1}N^{2}(1-N/K_{1})-cN,$ (4)

where $R_{1}$ and $K_{1}$ are positive parameters. We recover equation (3)
from (4) by taking $R_{1}/K_{1}=a$ and $K_{1}=b$.

The per capita growth rate, $G(N)=aN(b-N)-c$, in (3) is negative below the
threshold
$\displaystyle{N_{1}=\frac{b}{2}-\sqrt{\frac{b^{2}}{4}-\frac{c}{a}}}$ and
above the threshold
$\displaystyle{N_{2}=\frac{b}{2}+\sqrt{\frac{b^{2}}{4}-\frac{c}{a}}}$, when
$ab^{2}>4c$. Thus, equation (3) can be compared with the single species
population growth with strong Allee effect (1) with the revised
parametrization $r=aN_{2}^{2}$, $\theta=N_{1}$, and $K=N_{2}$. In summary, we
can consider equation (3) as the description of sex-structured single species
population growth which implicitly includes the mate-finding Allee effect with
the Allee threshold $N_{1}$. In this work, we consider equation (3) to
describe the growth rate of prey in the absence of specialist predator.

Predator-independent functional responses are widely used in prey-predator
systems, which is not always realistic in ecology [44]. If predator density
$P$ increases, then the per capita growth rate of predators may decrease. This
phenomenon is called predator interference [45] that can be incorporated into
system dynamics through a predator-dependent functional response. A well-known
model incorporating predator interference comes from the RM system by
replacing Holling type-II functional response with Beddington–DeAngelis
functional response

$F(N,P)=\frac{mN}{p+N+qP},$ (5)

where parameters $m,p$ and $q$ denote maximum predation rate, self-saturation
constant and predator mutual interference respectively [46, 47]. Inclusion of
mutual interference in predators induces Turing’s activator-inhibitor
structure in the system dynamics [47]. Depending on the parameter values $p$
and $q$, the Beddington–DeAngelis functional response can be converted into a
ratio-dependent functional response or predator-independent Holling type-II
functional response. Here, we incorporate prey-predator interaction through
the Beddington–DeAngelis functional response and compare the system dynamics
with the ratio-dependent functional response and Holling type-II functional
response.

Homogeneous and non-homogeneous stationary solutions of prey and predator
distributions correspond to one of the possible feeding strategies of
predators within a permanent habitat and hunting in a nearby area. To avoid
predators, some prey species migrate to nearby places [48, 49]. To be
successful in capturing prey, predator species need to follow the prey species
with some strategy which leads to the formation of various complex dynamic
patterns that include travelling wave, spatio-temporal chaos, and moving pulse
solution [50]. Predator invasion into the prey habitat leads to the formation
of a travelling wave that connects a predator-free steady state to a
coexisting steady state [51, 50]. The existence and non-existence of
travelling wave solutions and invasion speed are also interesting topics in
dynamical system theory [52, 53]. The impact of spatio-temporal chaos on
population dynamics is an active area of research in ecology. Appearance of
spatio-temporal chaos is widespread and it has also been observed in systems
that do not produce stationary pattern [23]. The stationary and moving pulse
solutions are studied in a prey-predator model following the prey growth rate
equation (3) and Holling type-II functional response [54].

Here, we consider a prey-predator model with prey’s growth rate as given in
(3) and Beddington-DeAngelis functional response (5) in the prey-predator
interaction. First, we perform a bifurcation analysis of the temporal model
and illustrate some representative dynamics through bifurcation diagrams. The
global existence of the solutions of the diffusive system and their global
asymptotic behavior are explored in various scenarios. Using energy estimates,
we obtain a priori bounds of global steady-state solutions and identify the
range of diffusion parameter for the non-existence of spatial patterns. The
temporal model exhibits bistable dynamics but the corresponding spatio-
temporal model possesses many non-constant spatio-temporal patterns that
include localized patterns, various Turing mode patterns, and biological
invasion, among others. We use stability analysis and Leray–Schauder degree
theory to show the existence of such non-constant steady states. A variety of
stable and unstable non-constant stationary solutions are obtained through
numerical continuation. Using numerical continuation, we illustrate the
bifurcation scenario of the spatial patterns that include Turing and localized
patterns. We also investigate the parametric region for the existence and non-
existence of the travelling wave solution and its profile inside the wedge-
shaped region of invasion.

The contents of this paper are as follows. Section 2 contains the description
of the temporal model along with its equilibria and their bifurcations. In
section 3, we extend it to spatio-temporal model with no-flux boundary
condition and investigate the existence and non-existence of constant and non-
constant solutions with their prior bounds. Using extensive numerical
continuation, we present various stationary mode-specific Turing pattern
solutions and localized pattern solutions together with their bifurcations in
section 4. In section 5, we discuss some complex dynamic patterns that include
travelling wave, moving pulse solution, and spatio-temporal chaos. Finally,
conclusions and discussions are presented in section 6.

##  2 Temporal Model

Let $N(T)$ and $P(T)$ respectively be the prey and predator densities at time
$T$. Suppose that the growth rate of the prey population is subjected to
reproductive Allee effect (3) and Beddington–DeAngelis functional response (5)
is chosen to represent the prey-predator interaction. Then, the governing
prey-predator system, subject to non-negative initial conditions, is

$\displaystyle\frac{dN}{dT}$
$\displaystyle=aN^{2}(b-N)-cN-\frac{mNP}{p+N+qP},$ (6a)
$\displaystyle\frac{dP}{dT}$ $\displaystyle=\frac{emNP}{p+N+qP}-dP,$ (6b)

where $d$ and $e$ respectively denote the per capita natural death rate of the
predator population and the conversation coefficient. All other parameters
have been described before. Introducing dimensionless variables
$u=\frac{N}{b}$, $v=\frac{mP}{bd}$ and $t=dT$, we obtain dimensionless version
of (6):

$\displaystyle\frac{du}{dt}$ $\displaystyle=\big{(}\sigma u^{2}(1-u)-\eta
u\big{)}-\frac{uv}{\alpha+u+\beta v}\equiv F_{1}(u,v),$ (7a)
$\displaystyle\frac{dv}{dt}$ $\displaystyle=\frac{\gamma uv}{\alpha+u+\beta
v}-v\equiv F_{2}(u,v),$ (7b)

where $\alpha=\frac{p}{b},\;\beta=\frac{qd}{m},\;\gamma=\frac{em}{d},$
$\sigma=\frac{ab}{d}$ and $\eta=\frac{c}{d}$ are dimensionless parameters.

###  2.1 Existence and stability of equilibria

The system (7) has trivial equilibrium point $E_{0}(0,0)$ irrespective of
parameter values. It has none, one or two axial equilibria depending on
parameter values. We denote the axial equilibria (whenever they exist) by
$E_{j}(u_{j},0)$, $j=1,2$. The $u$ components of the axial equilibria are
roots of the equation

$P(u)\equiv\sigma u(1-u)-\eta=0.$

For $\sigma>4\eta$, the system (7) has two axial equilibria $E_{1}(u_{1},0)$
and $E_{2}(u_{2},0)$, where

$u_{1}=\frac{\sigma+\sqrt{\sigma^{2}-4\sigma\eta}}{2\sigma}\text{ and
}u_{2}=\frac{\sigma-\sqrt{\sigma^{2}-4\sigma\eta}}{2\sigma}.$ (8)

On the other hand, it has only one axial equilibrium point $E_{1}(1/2,0)$ for
$\sigma=4\eta$ and no axial equilibrium point for $\sigma<4\eta$.

A interior equilibrium $E_{*}(u_{*},v_{*})$ is a point of intersection of the
nontrivial prey and predator nullclines $f_{1}(u,v)=0$ and $f_{2}(u,v)=0$,
where

$f_{1}(u,v)=\sigma u(1-u)-\eta-\frac{v}{\alpha+u+\beta v}\text{ and
}f_{2}(u,v)=\frac{\gamma u}{\alpha+u+\beta v}-1.$ (9)

The component $u_{*}$ of the interior equilibrium $E_{*}$ is a root of the
cubic equation

$Q(u)\equiv\sigma\gamma\beta{u}^{3}-\sigma\gamma\beta{u}^{2}+\left(\beta\eta\gamma+\gamma-1\right)u-\alpha=0.$

From (9), we get

$v_{*}={\frac{\gamma u_{*}-\alpha-u_{*}}{\beta}}.$

We must have $u_{*}>\frac{\alpha}{\gamma-1}\geq 0$ for the feasibility of
$E_{*}$. The system has only the trivial equilibrium point for $\sigma<4\eta$
since $f_{1}(u,v)<0$ in the positive quadrant of $u$-$v$ plane in this case.
For $\sigma>4\eta$, the coexisting equilibria $E_{*}$ is feasible if $u_{*}$
satisfies

$\text{max}\left\\{\frac{\alpha}{\gamma-1},u_{2}\right\\}<u_{*}<u_{1}.$

The Jacobian matrix of the system (7) evaluated at an equilibrium point
$E(u,v)$ is

$J(E)=\begin{bmatrix}\sigma u\left(2-3u\right)-\eta-{\frac{v}{\alpha+u+\beta
v}}+{\frac{uv}{\left(\alpha+u+\beta
v\right)^{2}}}&-{\frac{u\left(\alpha+u\right)}{\left(\alpha+u+\beta
v\right)^{2}}}\vspace{0.06in}\\\ {\frac{\gamma v\left(\alpha+\beta
v\right)}{\left(\alpha+u+\beta v\right)^{2}}}&{\frac{\gamma u}{\alpha+u+\beta
v}}-{\frac{\gamma uv\beta}{\left(\alpha+u+\beta
v\right)^{2}}}-1\end{bmatrix}.$ (10)

Clearly, trivial equilibrium point $E_{0}(0,0)$ is asymptotically stable since
both the eigenvalues of $J(E_{0})$ are negative. The eigenvalue of Jacobian
matrix evaluated at an axial equilibrium point $E_{j}(u_{j},0)$ ($j=1$ or $2$)
are $\sigma u_{j}(1-2{u_{j}})$ and $-1+{\frac{\gamma u_{j}}{\alpha+u_{j}}}.$
Hence, $E_{j}$ is asymptotically stable if
$\frac{1}{2}<u_{j}<\frac{\alpha}{\gamma-1}$, unstable if
$\frac{\alpha}{\gamma-1}<u_{j}<\frac{1}{2}$, and a saddle point if
$u_{j}>\text{max}\\{\frac{\alpha}{\gamma-1},\frac{1}{2}\\}$ or
$u_{j}<\text{min}\\{\frac{\alpha}{\gamma-1},\frac{1}{2}\\}.$

The Jacobian matrix evaluated at a coexisting equilibrium point
$E_{*}(u_{*},v_{*})$ is given by

$J(E_{*})=\begin{bmatrix}-u_{*}\left(2u_{*}-1\right)\sigma+{\frac{v_{*}}{\left(\beta
v_{*}+\alpha+u_{*}\right)^{2}}}&-{\frac{u_{*}\left(\alpha+u_{*}\right)}{\left(\beta
v_{*}+\alpha+u_{*}\right)^{2}}}\vspace{0.06in}\\\ {\frac{\gamma
v_{*}\left(\beta v_{*}+\alpha\right)}{\left(\beta
v_{*}+\alpha+u_{*}\right)^{2}}}&-{\frac{\beta v_{*}}{\beta
v_{*}+\alpha+u_{*}}}\end{bmatrix}\equiv\begin{bmatrix}a_{10}&a_{01}\vspace{0.06in}\\\
b_{10}&b_{01}\end{bmatrix}.$ (11)

Using Routh-Hurwitz stability criteria [55], $E_{*}$ is asymptotically stable
if

$\mathrm{T}:=\text{trace}(J(E_{*}))<0\text{ and
}\mathrm{D}:=\text{det}((J(E_{*}))>0.$

Due to unavailability of explicit expression of a interior equilibrium point
$E_{*}$, it is difficult to determine the analytic condition for the stability
of coexisting equilibrium point. However, we discuss the existence and
stability of of $E_{*}$ through various bifurcations that occur due to the
variation of control parameters.

###  2.2 Local bifurcation analysis

Here we discuss some preliminary bifurcation results that are required to
study of spatio-temporal pattern formation. In this subsection, we follow the
same notations as given in [55].

####  2.2.1 Saddle-node bifurcation

######  Theorem 1.

For $\gamma>2\alpha+1$, the temporal system (7) undergoes a saddle-node
bifurcation when the quadratic polynomial $P(u)$ has a double root.

###### Proof.

Suppose that $P(u)$ has a double root $u=u_{s}$ at $\sigma=\sigma_{\mbox{\tiny
SN}},$ i.e., $P(u_{s})=P^{\prime}(u_{s})=0$. Now, $P^{\prime}(u_{s})=0$ gives
$u_{s}=\frac{1}{2}$ and $P(\frac{1}{2})=0$ implies $\sigma_{\mbox{\tiny
SN}}=4\eta.$ From (8), we observe that $u_{1}>\frac{1}{2}$ and
$u_{2}<\frac{1}{2}$ for $\sigma>\sigma_{\mbox{\tiny SN}}$. Since
$\frac{\alpha}{\gamma-1}<\frac{1}{2}$, therefore $E_{1}$ is always a saddle
point and $E_{2}$ is an unstable node or a saddle point depending on the
parameter $\sigma$. These two axial equilibria merge into a single axial
equilibrium point $E_{\mbox{\tiny SN}}(\frac{1}{2},0)$ at
$\sigma=\sigma_{\mbox{\tiny SN}}$ and disappear for
$\sigma<\sigma_{\mbox{\tiny SN}}$. The Jacobian matrix evaluated at
$E_{\mbox{\tiny SN}}$ for $\sigma=\sigma_{\mbox{\tiny SN}}$ is

$J(E_{\mbox{\tiny
SN}})=\begin{bmatrix}0&-\frac{1}{(2\alpha+1)}\vspace{0.06in}\\\
0&\frac{2\gamma}{2\alpha+1}-1\end{bmatrix},$

which has a zero eigenvalue. Let $\phi=[0,\;1]^{T}$ and
$\psi=[2(\gamma-\alpha)-1,\;1]^{T}$ respectively be two eigenvectors of
$J(E_{\mbox{\tiny SN}})$ and $J(E_{\mbox{\tiny SN}})^{T}$ corresponding to
zero eigenvalue. Now, we verify the following transversality conditions:

$\displaystyle\psi^{T}{F}_{\sigma}({E_{\mbox{\tiny
SN}}};\sigma=\sigma_{\mbox{\tiny SN}})$ $\displaystyle=$
$\displaystyle\frac{2(\gamma-\alpha)-1}{8}>0,$
$\displaystyle\psi^{T}D^{2}{F}({E_{\mbox{\tiny
SN}}};\sigma=\sigma_{\mbox{\tiny SN}})(\phi,\phi)$ $\displaystyle=$
$\displaystyle-{\frac{4\beta\left(-\gamma+2\alpha+1\right)}{\left(2\alpha+1\right)^{2}}}\neq
0,$

where ${F}=\left[F_{1}(u,v),F_{2}(u,v)\right]^{T}$. Hence, the system (7)
undergoes a saddle-node bifurcation at $\sigma=4\eta.$ The pictorial
representation of the saddle-node bifurcation is shown in $\sigma$-$u$ plane
in Fig. 1(a). ∎

####  2.2.2 Transcritical bifurcation

######  Theorem 2.

For $\gamma>2\alpha+1$, the system (7) undergoes a transcritical bifurcation
around the axial equilibrium point $E_{\mbox{\tiny TC}}(u_{\mbox{\tiny
TC}},0)$ at $\sigma\equiv\sigma_{\mbox{\tiny
TC}}={\frac{\eta\left(\gamma-1\right)^{2}}{\alpha\left(\gamma-\alpha-1\right)}}$,
where $u_{{}_{\mbox{\tiny TC}}}=({\sigma_{\mbox{\tiny
TC}}-\sqrt{\sigma_{\mbox{\tiny TC}}^{2}-4\eta\sigma_{\mbox{\tiny
TC}}}})/{2\sigma_{\mbox{\tiny TC}}}$.

###### Proof.

For $\sigma<\sigma_{\mbox{\tiny TC}},$ the axial equilibrium point $E_{2}$ is
unstable which becomes saddle for $\sigma>\sigma_{\mbox{\tiny TC}}$ (see Fig.
1(a)). The unstable coexisting equilibrium point $E_{*}$ becomes feasible for
$\sigma>\sigma_{\mbox{\tiny TC}}.$

The Jacobian matrix evaluated at $E_{\mbox{\tiny TC}}(u_{{\mbox{\tiny
TC}}},0)$ for $\sigma=\sigma_{\mbox{\tiny TC}}$ is

$J(E_{\mbox{\tiny TC}})=\begin{bmatrix}\sigma(2u_{{\mbox{\tiny
TC}}}-3u_{{\mbox{\tiny TC}}}^{2})-\eta&-{u_{{\mbox{\tiny
TC}}}}/{(\alpha+u_{{\mbox{\tiny TC}}})}\vspace{0.06in}\\\ 0&0\end{bmatrix},$

which has a zero eigenvalue. Let $\phi=[p,\;1]^{T}$ and $\psi=[0,\;1]^{T}$ be
two eigenvectors corresponding to the zero eigenvalue of $J(E_{\mbox{\tiny
TC}})$ and $J(E_{\mbox{\tiny TC}})^{T}$ respectively, where

$p=\frac{u_{{\mbox{\tiny TC}}}}{(\alpha+u_{{\mbox{\tiny
TC}}})\big{(}\sigma(2u_{{\mbox{\tiny TC}}}-3u_{{\mbox{\tiny
TC}}}^{2})-\eta\big{)}}.$

Now, we verify the following transversality conditions:

$\displaystyle\psi^{T}{F}_{\sigma}({E_{\mbox{\tiny
TC}}};\sigma=\sigma_{{\mbox{\tiny TC}}})$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\psi^{T}D{F}_{\sigma}({E_{\mbox{\tiny
TC}}};\sigma=\sigma_{{\mbox{\tiny TC}}})\phi$ $\displaystyle=$ $\displaystyle
0,$ $\displaystyle\psi^{T}D^{2}{F}({E_{\mbox{\tiny
TC}}};\sigma=\sigma_{{\mbox{\tiny TC}}})(\phi,\phi)$ $\displaystyle=$
$\displaystyle{\frac{2\gamma u_{{\mbox{\tiny
TC}}}\beta}{\left(\alpha+u_{{\mbox{\tiny
TC}}}\right)^{2}}}\left(\frac{\alpha}{(\alpha+u_{{\mbox{\tiny TC}}})\sigma
u_{{\mbox{\tiny TC}}}(1-2u_{{\mbox{\tiny TC}}})}+\beta\right)\neq 0.$

However, the second transversality condition should be nonzero for non-
degenerate transcritical bifurcation [55]. Hence, the system (7) undergoes a
degenerate transcritical bifurcation [56] around $E_{\mbox{\tiny TC}}$ at
$\sigma=\sigma_{\mbox{\tiny TC}}.$

∎

(a)

(b)

Figure 1: (a) Bifurcation diagram in the $\sigma$-$u$ plane. (b) A zoom
version of Fig. (a) for $\sigma\in[1.7,1.9]$. Points marked TC and SN denote
the transcritical and saddle-node bifurcation thresholds respectively. Here,
solid green, red dashed and black dashed curves represent stable, unstable and
saddle branches of equilibria respectively. Further, solid red color curve
represents the maximum and minimum of $u$ of the stable limit cycle. Parameter
values are $\alpha=0.07,\beta=0.2,\gamma=1.2$ and $\eta=0.1$.

####  2.2.3 Hopf bifurcation

An unstable coexisting equilibrium point $E_{*}$ can become stable via a Hopf
Bifurcation when the trace of the Jacobian matrix evaluated at $E_{*}$ changes
from positive to negative value due to variation of a model parameter. Here,
we choose $\sigma$ as a Hopf bifurcation parameter.

Solving

$T_{\sigma_{{}_{H}}}\equiv tr[J(E_{*};\sigma=\sigma_{{}_{H}})]=0,$

we find

$\sigma_{{}_{H}}={\frac{\left(1-v_{*}{\beta}^{2}-\left(\alpha+u_{*}\right)\beta\right)v_{*}}{\gamma^{2}u_{*}^{3}\left(2u_{*}-1\right)}}.$

This is an implicit expression for $\sigma$ since $E_{*}(u_{*},v_{*})$ depends
on $\sigma$. The system (7) undergoes a Hopf bifurcation at
$\sigma=\sigma_{{}_{H}}$ if the following non-hyperbolicity and transversality
conditions are satisfied:

$\displaystyle C1:$
$\displaystyle\;\;D_{\sigma_{{}_{H}}}\equiv\mathrm{det}[J(E_{*};\sigma=\sigma_{{}_{H}})]>0,$
$\displaystyle C2:$
$\displaystyle\;\;\frac{d}{d\sigma}(T_{\sigma})|_{\sigma=\sigma_{{}_{H}}}\neq
0.$

For $\sigma>\sigma_{{}_{H}}$, $E_{*}$ is always asymptotically stable. A limit
cycle surrounding $E_{*}$ is generated through a Hopf bifurcation at
$\sigma=\sigma_{{}_{H}}.$ The stability of the limit cycle is determined by
the first Lyapunov coefficient $l_{1}$ [55]. The Hopf bifurcation is
supercritical when $l_{1}<0$ and results in a stable limit cycle. In contrary,
the subcritical Hopf bifurcation results in an unstable limit cycle for
$l_{1}>0$. Due to unavailability of explicit expression of the interior
equilibrium $E_{*}$, it is impossible to obtain the sign of $l_{1}$
analytically. However, we obtain Hopf bifurcation and its stability
numerically. We find that Hopf bifurcation occurs around
${E_{*}}=(0.5986,0.2486)$ at $\sigma_{{}_{H}}=1.8566$ for parameter values
$\alpha=0.07,\beta=0.2,\gamma=1.2$ and $\eta=0.1$. The corresponding first
Lyapunov coefficient $l_{1}=-22.7488\pi<0$, indicating that the Hopf
bifurcation is supercritical. The corresponding temporal bifurcation diagram
is plotted in the Fig. 1(b). The Hopf generating stable limit cycle vanishes
through a heteroclinic bifurcation at $\sigma_{Het}=1.789<\sigma_{H}.$ For
$\sigma<\sigma_{Het}$, trivial equilibrium point $E_{0}$ is globally
asymptotically stable for the temporal system (7). A decrease in prey growth
drives the system from stable coexistence to oscillatory coexistence and
further decrease results in system collapse through a global bifurcation.

##  3 Spatio-temporal Model

Random movement of the population species is taken into account by
incorporating diffusion term in the temporal model. For simplicity, we
consider the spatio-temporal system in one-dimensional spatial domain
$\Omega:=(0,L)$. The dimensionless spatio-temporal system is

$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial
x^{2}}+F_{1}(u,v),$ (12a) $\frac{\partial v}{\partial
t}=d\frac{\partial^{2}u}{\partial x^{2}}+F_{2}(u,v),$ (12b)

where $d$ is is ratio of the diffusion coefficients of predator and prey
populations. The above system is subjected to initial conditions
$u_{0}(x)\equiv u(x,0)\geq 0$, $v_{0}(x)\equiv v(x,0)\geq 0$ for $x\in\Omega$
and no-flux boundary conditions $\frac{\partial u}{\partial n}=\frac{\partial
v}{\partial n}=0$ for $x\in\partial\Omega:=\\{0,L\\}$.

###  3.1 Existence and boundedness of global solution

In order to ensure well-posedness of the model (12b), we establish the global
existence of the solutions of the system (12b) and their priori bounds under
specified initial and boundary condition. Let
$\mathscrsfs{D}=\Omega\times(0,\infty)$,
$\bar{\mathscrsfs{D}}={\bar{\Omega}}\times(0,\infty)$ and
$\mathscrsfs{B}=\partial\Omega\times(0,\infty)$.

######  Theorem 3.

For the system (12b), the following hold whenever $\sigma>4\eta$:

  * (a)

If $u_{0}(x)\geq 0$ and $v_{0}(x)\geq 0$, then $u(x,t)>0$ and $v(x,t)>0$ in
$\bar{\mathscrsfs{D}}$ unless $\big{(}u(x,t),v(x,t)\big{)}\equiv(0,0)$ in
$\bar{\mathscrsfs{D}}$.

  * (b)

If $\sup_{\bar{\Omega}}u_{0}(x)<u_{2}$, then
$\big{(}u(x,t),v(x,t)\big{)}\rightarrow(0,0)$ uniformly as
$t\rightarrow\infty.$

  * (c)

If $\gamma<1$, then $\big{(}u(x,t),v(x,t)\big{)}$ $\rightarrow$
$(u_{A}(x),0),$ where $u_{A}(x)$ is the solution of

$\begin{cases}\nabla^{2}u+\big{(}\sigma u^{2}(1-u)-\eta u\big{)}=0,\text{ for
}x\in\Omega,\\\ \frac{\partial u}{\partial n}=0,\text{ for
}x\in\partial\Omega.\end{cases}$

  * (d)

For any solution $\big{(}u(x,t),v(x,t)\big{)}$ of (12b),

$\limsup_{t\rightarrow\infty}\int_{\Omega}u(x,t)\,dx\leq u_{1}|\Omega|\mbox{
and
}\limsup_{t\rightarrow\infty}\int_{\Omega}v(x,t)\,dx\leq\gamma\left(1+\frac{\sigma}{4}-\eta\right)u_{1}|\Omega|.$

###### Proof.

(a) Since

$\frac{\partial F_{1}(u,v)}{\partial v}\leq 0\;\,\mbox{ and
}\;\,\frac{\partial F_{2}(u,v)}{\partial u}\geq
0\quad\mbox{for}\;\,(u,v)\in\mathbb{R}^{2}_{\geq 0}=\\{u\geq 0,v\geq 0\\},$

the system (12b) is a mixed quasi-monotone system [57, 58]. Let
$\big{(}\hat{u}(t),\hat{v}(t)\big{)}$ be the solution of the differential
equations

$\begin{cases}\frac{du}{dt}=\sigma u^{2}(1-u)-\eta u\equiv\sigma
u(u-u_{1})(u_{2}-u),\\\\[5.0pt] \frac{dv}{dt}=\displaystyle{\frac{\gamma
uv}{\alpha+u+\beta v}-v},\\\\[5.0pt] u(0)=\sup_{\bar{\Omega}}u_{0}(x)\text{
and }v(0)=\sup_{\bar{\Omega}}v_{0}(x),\end{cases}$ (13)

where $u_{1}$ and $u_{2}$ are given in (8). Clearly $\hat{u}(t)<\rho_{1}$ and
$\hat{v}(t)<\rho_{2}$ [17], where

$\rho_{1}=\text{max}\left\\{\hat{u}(0),u_{1}\right\\}\;\text{and}\;\rho_{2}=\text{max}\left\\{\hat{v}(0),u_{1}\left(1+\frac{\gamma}{\beta}+\frac{\sigma}{4}-\eta\right)\right\\}.$

Using the definition 8.1.2 in [57], we find that
$\big{(}\underline{u}(x,t),\underline{v}(x,t)\big{)}=(0,0)$ and
$\big{(}\overline{u}(x,t),\overline{v}(x,t)\big{)}=\big{(}\hat{u}(t),\hat{v}(t)\big{)}$
are lower and upper solutions of system (12b) since

$\frac{\partial\overline{u}(x,t)}{\partial
t}-\nabla^{2}\overline{u}(x,t)-F_{1}(\overline{u}(x,t),\underline{v}(x,t))=0\geq
0=\frac{\partial\underline{u}(x,t)}{\partial
t}-\nabla^{2}\underline{u}(x,t)-F_{1}(\underline{u}(x,t),\overline{v}(x,t)),$
$\frac{\partial\overline{v}(x,t)}{\partial
t}-\nabla^{2}\overline{v}(x,t)-F_{2}(\overline{u}(x,t),\overline{v}(x,t))=0\geq
0=\frac{\partial\underline{v}(x,t)}{\partial
t}-\nabla^{2}\underline{v}(x,t)-F_{2}(\underline{u}(x,t),\underline{v}(x,t)),$

for $(x,t)\in{\mathscrsfs{D}}$ with boundary conditions

$\frac{\partial\overline{u}(x,t)}{\partial n}\geq
0\geq\frac{\partial\underline{u}(x,t)}{\partial
n},\quad\frac{\partial\overline{v}(x,t)}{\partial n}\geq
0\geq\frac{\partial\underline{v}(x,t)}{\partial n}\;\mbox{ for
}(x,t)\in\mathscrsfs{B},$

and initial conditions

$\underline{u}(x,0)\leq u_{0}(x)\leq\overline{u}(x,0)\text{ and
}\underline{v}(x,0)\leq v_{0}(x)\leq\overline{v}(x,0)\;\text{ for
}x\in\Omega.$

Now, using theorem 8.3.3 of [57], we conclude that the system (12b) admits
unique globally defined solution $\big{(}u(x,t),v(x,t)\big{)}$ which satisfies

$0\leq u(x,t)\leq\hat{u}(t)\text{ and }0\leq v(x,t)\leq\hat{v}(t).$

From the expressions of $\rho_{1}$ and $\rho_{2}$ given above, we observe that

$F_{1}(\rho_{1},0)\leq 0\leq F_{1}(0,\rho_{2})\text{ and
}F_{2}(\rho_{1},\rho_{2})\leq 0\leq F_{2}(0,0).$

Hence, if we define

$S=\left\\{(u,v)\in C^{2}(\bar{\mathscrsfs{D}}):\;0\leq u\leq\rho_{1},0\leq
v\leq\rho_{2}\right\\},$

and $\big{(}u(x,0),v(x,0)\big{)}\in S$, then $u(x,t)>0,v(x,t)>0$ or
$\big{(}u(x,t),v(x,t)\big{)}\equiv(0,0)$ for $(x,t)\in\bar{\mathscrsfs{D}}$
(see theorem 8.3.3 of [57]).

(b) Since $\hat{u}(0)=\sup_{\bar{\Omega}}u_{0}(x)$, we have
$\hat{u}(0)<u_{2}$. Now, from the first equation of (13) we find that
$\hat{u}(t)\rightarrow 0$ as $t\rightarrow\infty$ and consequently
$\hat{v}(t)\rightarrow 0$ as $t\rightarrow\infty$. Thus,
$\big{(}u(x,t),v(x,t)\big{)}\rightarrow(0,0)$ uniformly as
$t\rightarrow\infty.$

(c) If $\gamma<1,$ then

$\frac{d\hat{v}}{dt}\leq(\gamma-1)\hat{v}<0,$

which gives $\hat{v}(t)\rightarrow 0$ as $t\rightarrow\infty$. Since, $0\leq
v(x,t)\leq\hat{v}(t)$, we conclude that ${v}(x,t)\rightarrow 0$ uniformly as
$t\rightarrow\infty$ for $x\in\bar{\Omega}.$ Therefore the limiting behavior
of $u(x,t)$ is determined by the semi-flow generated by the parabolic equation

$\begin{cases}u_{t}=\nabla^{2}u+\big{(}\sigma u^{2}(1-u)-\eta u\big{)},\text{
for }x\in\Omega,\,t>0,\\\ \frac{\partial u}{\partial n}=0,\text{ for
}x\in\partial\Omega.\end{cases}$ (14)

Now, every orbit of the gradient system (14) converges to a steady state
$u_{A}(x)$ [59]. Then, using the theory of asymptotically autonomous dynamical
system [59], we conclude that the solution
$(u(x,t),v(x,t))\rightarrow(u_{A}(x),0)$ as $t\rightarrow\infty$.

(d) If $\hat{u}(0)>u_{2},$ then $\hat{u}(t)\rightarrow u_{1}$ as
$t\rightarrow\infty$. Thus for any $\epsilon>0,$ there exist $T_{0}>0$ such
that $u(x,t)\leq u_{1}+\epsilon$ for
$(x,t)\in\bar{\Omega}\times(T_{0},\infty)$. Thus

$\limsup_{t\rightarrow\infty}\int_{\Omega}u(x,t)\,dx\leq u_{1}|\Omega|.$

To estimate the bound of $v(x,t)$, we introduce $P(t)=\int_{\Omega}u(x,t)\,dx$
and $Q(t)=\int_{\Omega}v(x,t)\,dx$. Then

$\displaystyle\frac{dP}{dt}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nabla^{2}u\;dx+\int_{\Omega}\left[\sigma
u^{2}(1-u)-\eta u-\frac{uv}{\alpha+u+\beta v}\right]dx,$
$\displaystyle\frac{dQ}{dt}$ $\displaystyle=$ $\displaystyle
d\int_{\Omega}\nabla^{2}v\;dx+\int_{\Omega}\left[\frac{\gamma
uv}{\alpha+u+\beta v}-v\right]dx.$

Using Neumann boundary condition, we obtain

$\begin{split}\frac{d}{dt}(\gamma P+Q)&=\int_{\Omega}\gamma\big{(}\sigma
u^{2}(1-u)-\eta u\big{)}\,dx-Q\,\leq\,-(\gamma
P+Q)+\gamma(1+\frac{\sigma}{4}-\eta)P\\\ &\leq-(\gamma
P+Q)+\gamma(1+\frac{\sigma}{4}-\eta)(u_{1}+\epsilon)|\Omega|\;\text{ for
}t>T_{0}.\end{split}$

Integrating the above gives

$\gamma
P(t)+Q(t)\leq\gamma(1+\frac{\sigma}{4}-\eta)(u_{1}+\epsilon)|\Omega|+\epsilon\;\,\text{
for }t>T_{1},$

where $T_{1}>T_{0}$. Since $P(t)\geq 0$, we find

$\limsup_{t\rightarrow\infty}Q(t)\leq\gamma(1+\frac{\sigma}{4}-\eta)u_{1}|\Omega|.$

∎

###  3.2  Homogeneous steady state analysis

The equilibrium points of the temporal model (7) correspond to homogeneous
steady-states of spatio-temporal model. Therefore, the system (12b) has
homogeneous steady-states corresponding to the trivial steady state $E_{0}$,
semi-trivial steady states $E_{1}$ and $E_{2}$, and the coexisting steady
state $E_{*}.$ Here we discuss the stability of these homogeneous steady state
solutions.

######  Theorem 4.

Suppose $\sigma>\sigma_{TC}$ and $\gamma>1$, then the following conditions
hold:

  * (a)

Trivial steady state $E_{0}(0,0)$ is always asymptotically stable.

  * (b)

Semi-trivial steady state $E_{1}(u_{1},0)$ and $E_{2}(u_{2},0)$ are always
unstable.

  * (c)

Coexisting steady state $E_{*}(u_{*},v_{*})$ is asymptotically stable if

$\sigma>\frac{v_{*}}{\gamma^{2}u_{*}^{3}(2u_{*}-1)}\equiv\sigma_{S},$

and unstable if $\sigma<\sigma_{{}_{H}}.$

###### Proof.

The linearized system of (12b) about a homogeneous steady state
$E(u_{e},v_{e})$ can be expressed as

$\frac{\partial W}{\partial t}=\mathcal{L}(W):=\mathcal{D}\nabla^{2}W+J(E)W,$

where $J(E)$ is defined in (10), $\mathcal{D}=$diag$(1,d)$ and $W\in
C^{2}(\bar{\mathscrsfs{D}})\times C^{2}(\bar{\mathscrsfs{D}}).$

Let $0=k_{0}<k_{1}<\cdots<k_{j}<\cdots$ be the eigenvalues and $E(k_{j})$ be
the eigenfunction space corresponding to $k_{j}$ for the eigenvalue problem

$\displaystyle-\nabla^{2}w=kw\quad\mbox{in}\;\Omega,$
$\displaystyle\frac{\partial w}{\partial
n}=0\qquad\quad\;\mbox{on}\;\partial\Omega.$

Further, suppose that
$\\{\psi_{i,j}:i=1,\cdots,\text{dim}\big{(}E(k_{j})\big{)}\\}$ be an
orthogonal basis set of $E(k_{j})$ and
$\mathcal{W}_{ij}=\\{c\psi_{i,j}:c=(c_{1},c_{2})^{T}\\}$. Let
$\mathcal{W}_{j}=\bigoplus_{i=1}^{\mbox{dim}(E(k_{j}))}\mathcal{W}_{ij}$ be
the direct sum of $\mathcal{W}_{ij}$. It can be shown that

$\mathcal{W}\equiv\left\\{\left(\phi,\psi\right)^{T}\in
C^{2}(\bar{\Omega})\times C^{2}(\bar{\Omega}):\frac{\partial\phi}{\partial
n}=\frac{\partial\psi}{\partial n}=0\;\text{ for
}x\in\partial\Omega\right\\}=\bigoplus_{j=1}^{\infty}\mathcal{W}_{j},$ (15)

and $\mathcal{W}_{j}$ is invariant under the operator $\mathcal{L}$. Now,
$\lambda$ is an eigenvalue of $\mathcal{L}$ if and only if $\lambda$ is an
eigenvalue of the matrix $\mathcal{L}_{j}=-k_{j}\mathcal{D}+J(E)$ for some
$j\geq 0.$ The characteristic equation of $\mathcal{L}_{j}$ is given by

$\text{det}(\lambda
I-\mathcal{L}_{j})=\lambda^{2}-\text{trace}(\mathcal{L}_{j})\lambda+\text{det}(\mathcal{L}_{j}),$

where
$\text{trace}(\mathcal{L}_{j})=a_{10}+b_{01}-(1+d)k_{j},\;\text{det}(\mathcal{L}_{j})=dk_{j}^{2}-(da_{10}+b_{01})k_{j}+(a_{10}b_{01}-a_{01}b_{10}).$

(a) At $E_{0}$, we find $J(E_{0})=\begin{pmatrix}-\eta&0\\\
0&-1\end{pmatrix},$ and
$\text{trace}(\mathcal{L}_{j})=-1-\eta-(1+d)k_{j}<0,\;\text{det}(\mathcal{L}_{j})=dk_{j}^{2}+(d\eta+1)k_{j}+\eta>0.$
Hence, $E_{0}$ is always asymptotically stable.

(b) For $E(u,0)$, we find $J(E)=\begin{pmatrix}\sigma
u(1-2{u})&\frac{u}{\alpha+u}\\\ 0&-1+{\frac{\gamma
u}{\alpha+u}}\end{pmatrix}.$ For $j=0$, $\mathcal{L}_{j}$ has one positive
eigenvalue $-1+\frac{\gamma u}{\alpha+u}>0$ since $\gamma>1$ and
$\sigma>\sigma_{TC}$. Thus, both $E_{1}$ and $E_{2}$ are unstable.

(c) From the Jacobian matrix $J(E_{*})$ given in (11), we find that
$a_{01}<0$, $b_{10}>0$ and $b_{01}<0.$ Now, $E_{*}$ is asymptotically stable
if and only if

$\text{trace}(\mathcal{L}_{j})<0\text{ and
}\text{det}(\mathcal{L}_{j})>0\text{ for all }j.$ (16)

If $\sigma>\sigma_{S}$, then $a_{10}<0$ and both the conditions in (16) are
satisfied. Thus, the steady state solution $E_{*}$ is stable when
$\sigma>\sigma_{S}$. On the other hand, $a_{10}+b_{01}>0$ for
$\sigma<\sigma_{H}$ which leads to $\text{trace}(\mathcal{L}_{j})>0$ for
$j=0$. Therefore, one of the eigenvalues of $\mathcal{L}_{0}$ must have
positive real part and the steady state $E_{*}$ becomes unstable. ∎

######  Remark 1.

From Theorem 5, we observe that $E_{*}$ is asymptotically stable for
$\sigma>\sigma_{S}$ and unstable for $\sigma<\sigma_{H}$ independent of the
diffusion parameter value $d$. Hence, the system (12b) shows bistability
between the trivial steady state $E_{0}$ and the coexisting steady state
$E_{*}$ for $\sigma>\sigma_{S}$. The system may show various non-constant
stationary solutions for $\sigma<\sigma_{S}$, which we discuss in section 4.

######  Remark 2.

Since the semi-trivial steady states are always saddle points, there is a
possibility of appearance of travelling wave solution which we discuss in
section 5.1.

###  3.3 Steady state solution

The spatio-temporal system (12b) admits time independent solutions
$(u(x),v(x))$ which satisfy

$\frac{d^{2}u}{dx^{2}}\,+\,F_{1}(u,v)\;=\;0,\;\;d\frac{d^{2}v}{dx^{2}}\,+\,F_{2}(u,v)\;=\;0,\quad\text{
for }x\in\Omega,$ (17)

along with the boundary conditions $\displaystyle{\frac{\partial u}{\partial
n}=\frac{\partial v}{\partial n}=0\text{ on }\partial\Omega.}$ Here we
determine the conditions which establish the existence and non-existence of
such solutions.

####  3.3.1 Nonexistence of non-constant steady state solution

Here we show the non-existence of non-constant solution. Before proving the
main theorem, we need the following proposition.

######  Proposition 1.

Suppose that $(u(x),v(x))$ is a non-negative nontrivial solution of (17). Then
$(u(x),v(x))$ is of the form of semi-trivial solution $(u(x),0)$ which satisfy
(14) or coexisting solution $(u(x),v(x))$ which satisfy

$0<u(x)<u_{1}\text{ and }0<v(x)<M^{*}\equiv\gamma
u_{1}\left(\frac{\sigma}{4}-\eta\right)\text{ for }x\in\Omega.$

###### Proof.

If $v(x_{0})=0$ for some $x\in\Omega$, then strong maximum principle implies
$v(x)\equiv 0$ and $u(x)$ satisfies (14). On the other hand, if there exists
$x_{0}\in\Omega$ such that $u(x_{0})=0$, then strong maximum principle gives
$u(x)\equiv 0$, which also implies that $v(x)\equiv 0$. Hence, the other
possibility is $u(x)>0$ and $v(x)>0$ for $x\in\Omega$.

From theorem 3(d), we have $u(x)\leq u_{1}$ for $x\in\Omega$ as
$t\rightarrow\infty$, and thus the strong maximum principle implies
$u(x)<u_{1}\text{ for }x\in\Omega.$ Now from (17), we find

$\begin{split}-(\gamma\nabla^{2}u+d\nabla^{2}v)=\gamma\big{(}\sigma
u^{2}(1-u)-\eta u\big{)}-v&=-\frac{1}{d}(\gamma u+dv)+\gamma u\left(\sigma
u(1-u)-\eta+\frac{1}{d}\right)\\\ &\leq-\frac{1}{d}(\gamma u+dv)+\gamma
u_{1}\left(\frac{\sigma}{4}-\eta+\frac{1}{d}\right).\end{split}$

Then using maximum principle theorem [60], we find

$(\gamma u+dv)\leq\gamma
du_{1}\left(\frac{\sigma}{4}-\eta+\frac{1}{d}\right)\;\text{ for }x\in\Omega,$

which leads to the desired estimate

$0<v(x)<M^{*}\;\text{ for }x\in\Omega.$

∎

######  Theorem 5.

If $\sigma>4\eta$ and $\gamma>1$, then there exists
$d^{*}=d^{*}(\alpha,\beta,\gamma,\sigma,\eta,|\Omega|)$ such that the system
(12b) does not admit a non-constant steady state solution for $d^{*}<d<1$.

###### Proof.

Let $\big{(}u(x),v(x)\big{)}$ be a non-constant solution of (17). Then

$\int_{\Omega}(u-\bar{u})dx=\int_{\Omega}(v-\bar{v})dx=0,\,\,\textrm{where}\,\,\bar{u}=\frac{1}{\Omega}\int_{\Omega}u(x)dx,\text{
and }\quad\bar{v}=\frac{1}{\Omega}\int_{\Omega}v(x)dx.$

Multiplying the first equation of (17) by $u-\bar{u}$ and integrating over
$\Omega$ using the no-flux boundary condition, we find

$\int_{\Omega}|\nabla(u-\bar{u})|^{2}dx=\int_{\Omega}(u-\bar{u})F_{1}(u,v)dx=\int_{\Omega}(u-\bar{u})\big{(}F_{1}(u,v)-F_{1}(\bar{u},\bar{v})\big{)}dx\equiv
I_{1}+I_{2},$

where.

$\displaystyle I_{1}$
$\displaystyle=\int_{\Omega}(u-\bar{u})\big{(}u(u-u_{1})(u_{2}-u)-\bar{u}(\bar{u}-u_{1})(u_{2}-\bar{u})\big{)}dx,\;$
$\displaystyle I_{2}$
$\displaystyle=\int_{\Omega}(\bar{u}-u)\left(\frac{uv}{\alpha+u+\beta
v}-\frac{\bar{u}\bar{v}}{\alpha+\bar{u}+\beta\bar{v}}\right)dx.$

Similarly, multiplying the second equation of (17) by $v-\bar{v}$ and
integrating over $\Omega$ using no-flux boundary condition, we find

$\begin{split}d\int_{\Omega}|\nabla(v-\bar{v})|^{2}dx=&\int_{\Omega}(v-\bar{v})F_{2}(u,v)dx\\\
=&\int_{\Omega}(v-\bar{v})^{2}\left(\frac{\gamma u}{\alpha+u+\beta
v}-1\right)dx+I_{3},\\\ \end{split}$

where

$\begin{split}I_{3}=&\int_{\Omega}(v-\bar{v})\bar{v}\left(\frac{\gamma
u}{\alpha+u+\beta v}-1\right)dx.\end{split}$

Using upper bounds of $I_{1}$, $I_{2}$ and $I_{3}$ (see Appendix A), we get

$\int_{\Omega}|\nabla(u-\bar{u})|^{2}dx\leq\left(\frac{u_{1}}{2\alpha}+u_{1}(2-u_{2})\right)\int_{\Omega}(u-\bar{u})^{2}dx+\frac{u_{1}}{2\alpha}\int_{\Omega}(v-\bar{v})^{2}dx.$
(18)

$d\int_{\Omega}|\nabla(v-\bar{v})|^{2}dx\leq\frac{\gamma}{2\beta}\int_{\Omega}(u-\bar{u})^{2}dx+\left(\frac{\gamma}{2\beta}+(\gamma-1)\right)\int_{\Omega}(v-\bar{v})^{2}dx.$
(19)

Adding (18), (19) and using Poincare inequality, we find

$\int_{\Omega}|\nabla(u-\bar{u})|^{2}dx+d\int_{\Omega}|\nabla(v-\bar{v})|^{2}dx\leq\frac{1}{k_{1}}\left(A\int_{\Omega}|\nabla(u-\bar{u})|^{2}dx+B\int_{\Omega}|\nabla(v-\bar{v})|^{2}dx\right),$
(20)

where

$A=\frac{\gamma}{2\beta}+\frac{u_{1}}{2\alpha}+u_{1}(2-u_{2})\;\text{ and
}B=\frac{\gamma}{2\beta}+(\gamma-1)+\frac{u_{1}}{2\alpha}.$

If

$\text{min}\\{1,d\\}\geq\frac{1}{k_{1}}\text{max}\\{A,B\\}\equiv d^{*},$

then

$\nabla(u-\bar{u})=\nabla(v-\bar{v})=0,$

which leads to constant solution. ∎

######  Remark 3.

Theorem 5 is not applicable for any parameter set when $d^{*}>1$. For
$d^{*}<1,$ the first wavenumber $k_{1}$ must be large which leads to a small
domain size. To verify theorem 5 with numerical simulation, consider a one-
dimensional domain of length $|\Omega|=1$ with parameters values $\eta=0.1,$
$\alpha=0.2,$ $\beta=2.4,$ $\gamma=1.3,$ and $\sigma=1.5.$ Then we obtain
$d^{*}=0.4439$ with these parameter values and the system always shows a
homogeneous steady-state corresponding to either the trivial steady state
$E_{0}$ or the coexistence steady state $E_{*}$ for $d^{*}<d<1.$

###  3.4  Existence of non-constant steady state solution

Here we establish the existence of positive non-constant steady solution of
(12b) using Leray–Schauder degree theory. We define
$\displaystyle{W=(u,v)^{T},\;W_{*}=(u_{*},v_{*})^{T},\;\mathcal{W}^{+}=\\{W\in\mathcal{W}\;:u,v>0\,\text{
on }\bar{\Omega}\\},}$ and
$\displaystyle{\mathcal{F}(W)=(F_{1},d^{-1}F_{2})^{T},}$ where $\mathcal{W}$
has been defined in (15). We write the system (17) as

$\begin{cases}-\nabla^{2}W=\mathcal{F}(W),&x\in\Omega,\\\ \frac{\partial
W}{\partial n}=0,&x\in\partial\Omega.\end{cases}$ (21)

Now the system (21) has a positive solution $W$ if and only if

$\mathcal{T}(W):=W-(\mathcal{I}-\nabla^{2})^{-1}\\{\mathcal{F}(W)+W\\}=0\;\text{
for }W\in\mathcal{W}^{+},$

where $\mathcal{I}$ is the identity operator, $(\mathcal{I}-\nabla^{2})^{-1}$
is the inverse of $(\mathcal{I}-\nabla^{2})$ under Neumann boundary condition.
We observe that

$D_{W}\mathcal{F}(W_{*})=\mathcal{A},\;\text{and
}D_{W}\mathcal{T}(W_{*})=\mathcal{I}-(\mathcal{I}-\nabla^{2})^{-1}(\mathcal{I}+\mathcal{A}),$

where

$\mathcal{A}=\begin{pmatrix}a_{10}&a_{01}\\\
d^{-1}b_{10}&d^{-1}b_{01}\end{pmatrix}.$

Now, $\lambda$ is an eigenvalue of $D_{W}\mathcal{T}(W_{*})$ if and only if
$\lambda(1+k_{j})$ is an eigenvalue of the matrix
$(k_{j}\mathcal{I}-\mathcal{A}).$ Hence, $D_{W}\mathcal{T}(W_{*})$ is
invertible if and only if

$\mathcal{H}(k_{j},d):=\text{det}(k_{j}\mathcal{I}-\mathcal{A})=d^{-1}\big{(}dk_{j}^{2}-k_{j}(da_{10}+b_{01})+\mathrm{D}\big{)}\neq
0\text{ for all }j\geq 0.$

By the Leray–Schauder theorem, we know that if $D_{W}\mathcal{T}(W_{*})$ is
invertible then

$\text{index}(\mathcal{T(.),W_{*}})=(-1)^{\xi},$

where $\xi$ is the sum of the algebraic multiplicities of all the negative
eigenvalues of $D_{W}\mathcal{T}(W_{*}).$ This leads to the following lemma
[61]:

######  Lemma 6.

If $\mathcal{H}(k_{j},d)\neq 0$ for all $j\geq 1$, then

$\mathrm{index}(\mathcal{T(.),W_{*}})=(-1)^{\xi},$

where

$\xi:=\sum_{j\geq 1,\;\mathcal{H}(k_{j},d)<0}\mathrm{dim}\,E(k_{j}).$

We now describe a theorem which guarantees the existence of at least one
positive non-constant solution. For this we assume that

$a_{10}>0\;\mathrm{and}\;a_{10}d+b_{01}>2\sqrt{d\mathrm{D}}>0.$ (22)

Then $\mathcal{H}(k,d)=0$ has two real roots $k_{\pm}$, where

$k_{\pm}=\frac{(da_{10}+b_{01})\pm\sqrt{(da_{10}+b_{01})^{2}-4d\mathrm{D}}}{2d}.$

######  Theorem 7.

Suppose that (22) holds. If there exist some integers $0\leq n_{1}<n_{2}$ such
that

$k_{-}\in(k_{n_{1}},k_{n_{1}+1}),\;k_{+}\in(k_{n_{2}},k_{n_{2}+1}),$ (23)

and $\sum_{j=n_{1}+1}^{n_{2}}\text{dim}E(k_{j})$ is odd, then model (12b) has
at least one positive non-constant solution.

###### Proof.

We prove this theorem using topological degree with homotopy invariance.
Theorem 5 implies the existence of $\bar{d}<d$ such that

  1. (i)

system (17) with diffusion coefficient $\bar{d}$ has no non-constant
solutions,

  2. (ii)

$\mathcal{H}(k_{j},\bar{d})>0$ for all $j\geq 0.$

For $\tau\in[0,1]$, we define

$\overline{\mathcal{F}}(\tau,W)=\left(F_{1},\big{(}\tau
d+(1-\tau)\bar{d})^{-1}\big{)}F_{2}\right)^{T}$

and consider the problem

$\begin{cases}-\nabla^{2}W=\overline{\mathcal{F}}(\tau,W)&x\in\Omega,\\\
\frac{\partial W}{\partial n}=0&x\in\partial\Omega.\end{cases}$ (24)

Then the system (24) has positive solution $W$ if and only if

$\overline{\mathcal{T}}(\tau,W)=W-(\mathcal{I}-\nabla^{2})^{-1}\big{(}\overline{\mathcal{F}}(\tau,W)+W\big{)}=0\text{
for }W\in\mathcal{W}^{+}.$

A straightforward calculation yields

$D_{W}\overline{\mathcal{T}}(\tau,W_{*})=\mathcal{I}-(\mathcal{I}-\nabla)^{-1}(\mathcal{I}+\overline{\mathcal{A}}),$

where

$\overline{\mathcal{A}}=\begin{pmatrix}a_{10}&a_{01}\\\ (\tau
d+(1-\tau)\bar{d})^{-1})b_{10}&(\tau
d+(1-\tau)\bar{d})^{-1})b_{01}\end{pmatrix}.$

Thus, $W$ is a positive non-constant solution of (17) if and only if it is a
solution of (24) when $\tau=1.$

Observe that $\overline{\mathcal{F}}(1,W)=\mathcal{F}(W),\text{ and
}D_{W}\overline{\mathcal{T}}(1,W_{*})=D_{W}{\mathcal{T}}(W_{*}).$ Now from
(23), we find

$\begin{cases}\mathcal{H}(k_{j},d)<0&\text{ for }n_{1}+1\leq j\leq n_{2},\\\
\mathcal{H}(k_{j},d)>0&\text{ otherwise}.\end{cases}$ (25)

Further, using Lemma 1, we get

$\xi=\sum_{j\geq
0,\;{H}(k_{j},d)<0}\text{dim}E(k_{j})=\sum_{j=n_{1}+1}^{n_{2}}\text{dim}E(k_{j}),$

which is an odd number. Therefore

$\text{index}(\overline{\mathcal{T}}(1,W_{*})=(-1)^{\xi}=-1.$

Using (ii), we obtain

$\text{index}(\overline{\mathcal{T}}(0,W_{*})=(-1)^{0}=1.$

Using proposition 1, we observe the existence $\tilde{M}>0$ such that
$\tilde{M}^{-1}<u,v<\tilde{M}$ for all $\tau\in[0,1]$ and thus
$\overline{\mathcal{T}}(\tau,W)\neq 0.$ Consider the set
$\mathbb{B}(\tilde{M})$ defined by

$\mathbb{B}(\tilde{M})=\left\\{W\in\mathcal{W}\;:\tilde{M}^{-1}<u,v<\tilde{M}\right\\}.$

If the system has unique solution $W_{*}$ in $\mathbb{B}(\tilde{M}),$ then
application of the homotopy invariance of the topological degree gives

$\text{deg}\left(\overline{\mathcal{T}}(0,\cdot),0,\mathbb{B}(\tilde{M})\right)=\text{deg}\left(\overline{\mathcal{T}}(1,\cdot),0,\mathbb{B}(\tilde{M})\right).$
(26)

Also, both the equations $\overline{\mathcal{T}}(0,W)=0$ and
$\overline{\mathcal{T}}(1,W)=0$ have unique positive solution $W_{*}$ in
$\mathbb{B}(\tilde{M})$. Now,

$\text{deg}\left(\overline{\mathcal{T}}(0,\cdot),0,\mathbb{B}(\tilde{M})\right)=\text{index}(\overline{\mathcal{T}}(0,W_{*})=(-1)^{0}=1$

and

$\text{deg}\left(\overline{\mathcal{T}}(1,\cdot),0,\mathbb{B}(\tilde{M})\right)=\text{index}(\overline{\mathcal{T}}(1,W_{*})=(-1)^{\xi}=-1.$

This contradicts (26) and therefore, $\overline{\mathcal{T}}(1,W)$ has at
least one positive solution other than $W_{*}$. Since the system (12b) has
unique homogeneous steady state $W_{*}$, the system (12b) has at least one
positive non-constant steady-sate solution in $\mathbb{B}(\tilde{M})$. ∎

##  4 Continuation of stationary solution

In the previous section, we have discussed the existence and the non-existence
of non-constant solutions. Now, we are interested to find the origin of
various non-constant solutions that include mode-dependent Turing and
localized solutions and their continuation. The stationary solution of the
system (12b) can be analyzed by converting the system (17) into a system of
four first-order ordinary differential equations:

$\frac{du}{dx}=w,\;\frac{dv}{dx}=z,\;\frac{dw}{dx}=-F_{1}(u,v),\;\frac{dz}{dx}=-\frac{1}{d}F_{2}(u,v).$

The corresponding linearized system at $(u_{*},v_{*},0,0)$ can be written as

$\frac{dZ}{dx}=BZ,$

where

$Z=\begin{pmatrix}u\\\ v\\\ w\\\ z\end{pmatrix}\text{ and
}B=\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\ -a_{10}&-a_{01}&0&0\\\
-b_{10}/d&-b_{01}/d&0&0\end{bmatrix}.$

The characteristic equation of the matrix $B$ is given by

$H(\lambda^{2}):=d\lambda^{4}+(da_{10}+b_{01})\lambda^{2}+\mathrm{D}=0.$ (27)

$H(\lambda^{2})$ is a quadratic polynomial in $\lambda^{2}.$ The zeros of (27)
are always symmetric about real axis and imaginary axis in the complex plane.
At the Turing bifurcation and BD transition thresholds, (27) has a pair of
eigenvalues with multiplicity two. These occur when

$\mathrm{D}=\frac{(da_{10}+b_{01})^{2}}{4d}$ (28)

and $H(\lambda^{2})$ reduces to

$H(\lambda^{2})=(\lambda^{2}-K)^{2},\text{ where
}K=\frac{da_{10}+b_{01}}{2d}.$

Note that a Turing bifurcation corresponds to $K<0$, and a BD transition
corresponds to $K>0$ [19, 18, 62]. The threshold values of Turing bifurcation
and BD transition satisfy (28). Now, we verify these results with numerical
simulation and numerical continuation.

We keep the parameter values $\eta=0.1,$ $\alpha=0.07,$ $\beta=0.2$ and
$\gamma=1.2$ fixed in all the calculations. We consider $\sigma$ as a
bifurcation parameter. The corresponding Turing bifurcation threshold and BD
transition threshold are $\sigma_{T}=1.861$ and $\sigma_{BD}=2.098$
respectively. Both the bifurcation thresholds are marked in Fig. 2(b) along
with $||u||_{2}$ of the homogeneous steady state, where

$||u||_{2}=\sqrt{\int_{\Omega}(u^{2})dx}.$

Various localized structured solutions, which are different from Turing
solution, emerge from the BD transition point.

(a)

(b)

Figure 2: (a) Real and imaginary part of eigenvalues of the characteristic
equation (27). (b) Locations of the eigenvalues along with coexisting
homogeneous steady state (HSS) branch. Other parameter values are
$\alpha=0.07,\beta=0.2,\gamma=1.2,\eta=0.1$ and $d=46$.

(a)

(b)

(c)

Figure 3:  (a) Bifurcation diagram of stationary Turing patterns and localized
patterns. Here red color curve represents the homogeneous branch and green
color curves represent localized patterns. Other color curves represent Turing
stationary branches of different modes. Solid color curves represent the
stable branches and light color curves represent unstable branches. (b) A zoom
version of the bifurcation diagram for Turing mode solutions. (c) A schematic
bifurcation diagram of localized pattern solutions. Other parameter values are
$\eta=0.1,\alpha=0.07,\beta=0.2,\gamma=1.2$ and $d=46.$

To track a stationary solution against a bifurcation parameter, we use the
Matlab software package pde2path [63]. We consider the spatial-domain size
$L=200$, and diffusion parameter $d=46$. Numerical continuation of Turing mode
solutions and localized structure solutions are plotted in Fig. 3(a).
Different modes of Turing solutions emerge from different points (BP) [see
Fig. 3(b)]. The BP points can be found by solving for $\sigma$ (implicitly) in
the following equation

$d\left(\frac{n\pi}{L}\right)^{4}+(da_{10}+b_{01})\left(\frac{n\pi}{L}\right)^{2}+\mathrm{D}=0,$
(29)

where $n$ is the corresponding unstable Turing mode. For example, BP1, BP2,
and BP3 correspond to $n=19,20$ and $21$ respectively. We have shown 19-, 20-
and 30- mode Turing solutions and localized solutions in the Appendix B.

In Fig. 3, we observe existence of multiple branches of various Turing mode
solutions and localized solutions. The parameter value $\sigma=1.86$ satisfies
Turing instability criteria. However, the system (12b) does not produce a
stable Turing solution when a numerical simulation is performed with small
amplitude spatial heterogeneity around the coexisting steady state $E_{*}$ as
an initial condition. The evolution of numerical solution shows a long
oscillatory transition between a localized solution and Turing solution which
finally settles down into a steady localized solution (see Fig. 4). Thus we
conclude that the localized solution dominates the Turing solution.

(a)

(b)

(c)

Figure 4: Evolution of solution $u$: (a) spatial average of $u$ against time,
(b) space-time plot corresponding to the red rectangle [see (a)] of the time-
series plot, (c) space-time plot corresponding to the green rectangle [see
(a)] of the time-series plot. Here, parameter values are
$\alpha=0.07,\beta=0.2,\gamma=1.2,\eta=0.1,\sigma=1.86$ and $d=46$.

##  5 Time varying solutions

Here, we study some time-dependent solutions, namely, travelling wave, spatio-
temporal chaos and moving pulse solutions.

###  5.1 Travelling wave solution

Heterogeneous stationary solutions are obtained with small amplitude
heterogeneous perturbation around the unstable coexisting homogeneous steady-
states. However, the introduction of predator over a small domain with the
rest filled with prey population reveals invasive spread of the predators. We
have observed that the coexisting homogeneous steady-state
$E_{*}(u_{*},v_{*})$ may be stable or unstable and the predator free
homogeneous steady-state $E_{1}(u_{1},0)$ is always saddle. Consider the case
in which $E_{*}$ is a stable homogeneous steady-state. Here, we investigate
the stable and unstable regions of the travelling wave solution

$u(x,t)=U(\xi)\quad\mbox{and}\quad v(x,t)=V(\xi)$

connecting the steady states $E_{1}$ and $E_{*},$ where $\xi=x+ct$ and $c$ is
the wave speed. Now, $\big{(}U(\xi),V(\xi)\big{)}$ satisfy

$c\frac{dU}{d\xi}=\frac{d^{2}U}{d\xi^{2}}+F_{1}(U,V),$ (30a)
$c\frac{dV}{d\xi}=d\frac{d^{2}V}{d\xi^{2}}+F_{2}(U,V),$ (30b)

along with

$\big{(}U(-\infty),V(-\infty)\big{)}=E_{1}\text{ and
}\big{(}U(\infty),V(\infty)\big{)}=E_{*}.$

In terms new variables defined by [52]

$\displaystyle
X(t)\,=\,U(ct),\,\,Y(t)\,=\,\frac{1}{c}\left(cU(ct)-U^{\prime}(ct)\right),\,\,$
$\displaystyle
W(t)\,=\,V(ct),\,\,Z(t)\,=\,\frac{1}{c}\left(cV(ct)-dV^{\prime}(ct)\right),$

the system (30b) transforms into a system four first order ordinary
differential equations:

$X^{\prime}=c^{2}(X-Y),\;Y^{\prime}=F_{1}(X,W),\;W^{\prime}=\frac{c^{2}}{d}(W-Z)),\;Z^{\prime}=F_{2}(X,W).$
(31)

The steady states of the system (31) corresponding to $E_{1}$ and $E_{*}$ are
$\mathcal{E}_{1}(u_{1},u_{1},0,0)$ and
$\mathcal{E}_{*}(u_{*},u_{*},v_{*},v_{*}).$ Now, we construct a wedge-shaped
region in $\mathbb{R}^{4}$ such that $\mathcal{E}_{*}$ lies in the interior of
this region, and look for a solution connecting the steady states
$\mathcal{E}_{1}$ and $\mathcal{E}_{*}$ confined within this region. If
$0<X<u_{1}$ and $W>0,$ then

$-W<Z^{\prime}=F_{2}(X,W)<(\gamma-1)W.$

Using the above inequality, our first task is to compare the vector field of
the system (31) with the following planner system

$W^{\prime}=\frac{c^{2}}{d}(W-Z))\text{ and }Z^{\prime}=-W.$ (32)

######  Proposition 2.

The system (32) has a strictly monotone decreasing solution
$W_{1}(t),Z_{1}(t)$ for $t\in(-\infty,\infty)$ which satisfy
$Z_{1}(t)>W_{1}(t)>0$ and

$W_{1}(t)\rightarrow 0,Z_{1}(t)\rightarrow 0\text{ as }t\rightarrow\infty,$
$W_{1}(t)\rightarrow\infty,Z_{1}(t)\rightarrow\infty\text{ as
}t\rightarrow-\infty.$

###### Proof.

The linear system (32) possesses a negative eigenvalue

$\lambda_{1}=\frac{c^{2}-\sqrt{c^{4}+4dc^{2}}}{2d}$

with an eigenvector

$v_{1}=(1,m)^{T},\text{ where }m=\frac{c^{2}+\sqrt{c^{4}+4dc^{2}}}{2c^{2}}>1.$

Thus, the system (32) has a local one-dimensional stable manifold in the
triangular region

$S=\big{\\{}(W,Z):Z>W\geq 0\big{\\}}\cup(0,0).$

If a solution $(W_{1}(t),Z_{1}(t))$ starts from a point $(W_{0},Z_{0})\in S,$
then both $W_{1}(t),\;Z_{1}(t)\rightarrow 0\text{ as }t\rightarrow\infty.$

Further, we can easily observe that the vector field of the system (31) on the
side $\\{Z=W\\}$ of $S$ points vertically downward and the vector field in the
side $W=0$ points to the left horizontally. Hence, the solution
$\big{(}W_{1}(t),Z_{1}(t)\big{)}$ can be extended for all $t<0$ and
$\big{(}W_{1}(t),Z_{1}(t)\big{)}$ remains in the triangular region $S$ for all
$t<0$. Since $W_{1}^{\prime}(t)=\frac{c^{2}}{d}(W_{1}-Z_{1})<0$ and
$Z_{1}^{\prime}(t)<0$, we conclude that both
$W_{1}(t),Z_{1}(t)\rightarrow\infty\text{ as }t\rightarrow-\infty.$ ∎

Since, $W_{1}(t)$ is a strictly decreasing function, $W_{1}$ has an inverse
$W_{1}^{-1}:(0,\infty):\rightarrow\mathbb{R}$ such that $W_{1}(t)=W$ if and
only if $t=W_{1}^{-1}(W)$ for all $W>0.$ Therefore, we can express the stable
manifold of the system (32) inside $S$ as the graph of the function
$\zeta:[0,\infty)\rightarrow[0,\infty)$ by

$\zeta(W)=Z_{1}\big{(}W_{1}^{-1}(W)\big{)}\text{ for }W>0\;\text{ with
}\zeta(0)=0.$

Now we define a wedge-shaped region $\Sigma\in\mathbb{R}^{4}$ [52] as follows:

$\Sigma=\left\\{(X,Y,W,Z):0\leq X\leq u_{1},Y\in\mathbb{R},W\geq
0,\frac{1}{2}W\leq Z\leq\zeta(W)\right\\}.$

The boundary of $\Sigma$ consists of surfaces $B_{1}-B_{4}$ and $C_{3}-C_{5}$
defined by

$\displaystyle B_{1}$
$\displaystyle=\left\\{0<X<u_{1},Y\in\mathbb{R},W>0,Z=\zeta(W)\right\\},$
$\displaystyle B_{2}$
$\displaystyle=\left\\{0<X<u_{1},Y\in\mathbb{R},W>0,Z=\frac{1}{2}W\right\\},$
$\displaystyle B_{3}$
$\displaystyle=\left\\{X=u_{1},Y<u_{1},0<\frac{1}{2}W\leq
Z\leq\zeta(W)\right\\},$ $\displaystyle B_{4}$
$\displaystyle=\left\\{X=0,0<Y,0<\frac{1}{2}W\leq Z\leq\zeta(W)\right\\},$
$\displaystyle C_{3}$ $\displaystyle=\left\\{X=u_{1},Y\geq
u_{1},0<\frac{1}{2}W\leq Z\leq\zeta(W)\right\\},$ $\displaystyle C_{4}$
$\displaystyle=\left\\{X=0,0\leq Y,0<\frac{1}{2}W\leq Z\leq\zeta(W)\right\\},$
$\displaystyle C_{5}$ $\displaystyle=\left\\{0\leq X\leq
u_{1},W=Z=0\right\\}.$

Let $\Phi(t,\textbf{p})=(X(t),Y(t),W(t),Z(t))$ with
$\Phi(0,\textbf{p})=\textbf{p}\in\mathbb{R}^{4}$ be the flow of the system
(31). If $M={\frac{1}{d}\left(\frac{\gamma u_{1}}{\alpha+u_{1}}-1\right)},$
then the flow $\Phi(t,\textbf{p})$ satisfies the following proposition [52]:

######  Proposition 3.

  * (i)

If $c\geq 2d\sqrt{M},$ then the system (31) has a positive bounded travelling
wave $\Phi(t,\textbf{p})=(X(t),Y(t),W(t),Z(t))\in\Sigma$ with
$\Phi(t,\textbf{p})\rightarrow\mathcal{E}_{1}$ as $t\rightarrow-\infty$ and
$\Phi(t,\textbf{p})\rightarrow\mathcal{E}_{*}$ as $t\rightarrow+\infty$.

  * (ii)

If $0<c<2d\sqrt{M}$, then the system (31) does not have a non-negative,
nontrivial travelling wave solution connected to $\mathcal{E}_{1}$.

Next, we verify the heteroclinic connection between $\mathcal{E_{*}}$ and
$\mathcal{E}_{1}$ together with stable and unstable regions of monotonic and
non-monotonic travelling waves using numerical simulations.

####  5.1.1 Numerical Results

The Jacobian of the system (31) evaluated at $\mathcal{E}_{1}$ is

$J(\mathcal{E}_{1})=\begin{pmatrix}c^{2}&-c^{2}&0&0\\\ j_{1}&0&j_{2}&0\\\
0&0&\frac{c^{2}}{d}&-\frac{c^{2}}{d}\\\ 0&0&j_{3}&0\end{pmatrix},$

where

$j_{1}=\sigma
u_{1}(1-2{u_{1}})<0,\;j_{2}=-{\frac{u_{1}}{\alpha+u_{1}}}<0\;\;\mbox{ and
}\;\;j_{3}={\frac{1}{d}\left(\frac{\gamma u_{1}}{\alpha+u_{1}}-1\right)}>0.$

The eigenvalues of $J(\mathcal{E}_{1})$ are

$\lambda_{1,2}=\frac{c^{2}\pm
c\sqrt{c^{2}-4j_{1}}}{2},\;\,\lambda_{3,4}=\frac{c^{2}\pm
c\sqrt{c^{2}-4d^{2}j_{3}}}{2d},$

where $\lambda_{1,2}$ are always real but $\lambda_{3,4}$ may be real or
complex conjugate. The complex conjugate eigenvalues correspond to spiral
solution around $\mathcal{E}_{1}$ which leads to negative population density.
Hence, $\lambda_{3,4}$ must be real for the existence of travelling wave,
which implies $c\geq 2d\sqrt{j_{3}}$. Hence, the minimum wave speed is
$c_{min}=2d\sqrt{j_{3}}$. The Jacobian of the system (31) evaluated at
$\mathcal{E}_{*}$ is given by

$J(\mathcal{E}_{*})=\begin{pmatrix}c^{2}&-c^{2}&0&0\\\ a_{10}&0&a_{01}&0\\\
0&0&\frac{c^{2}}{d}&-\frac{c^{2}}{d}\\\ b_{10}&0&b_{01}&0\end{pmatrix}.$

We find that two of the eigenvalues of $J(\mathcal{E}_{*})$ have negative real
parts. There is a heteroclinic orbit $\Phi(t)$ contained in stable manifold of
$\mathcal{E}_{*}$ and unstable manifold of $\mathcal{E}_{1},$ i.e.,
$\Phi(t)\rightarrow\mathcal{E}_{1}$ as $t\rightarrow-\infty$ and
$\Phi(t)\rightarrow\mathcal{E}_{*}$ as $t\rightarrow+\infty$. We have plotted
that heteroclinic connection $\Phi(t)$ in Fig. 5. We also plot the boundaries
$B_{1}$ and $B_{2}$ in Fig. 5(b) which shows that $\Phi(t)$ is contained in
the wedged-shaped region $\Sigma.$

(a)

(b)

Figure 5:  Heteroclinic connection of the system (31) between
$\mathcal{E}_{1}(0.9615,0.9615,0,0)$ and
$\mathcal{E}_{*}(0.7291,0.7291,0.3791,0.3791)$ for $d=46$, $\sigma=2.7$ and
$c=5.9$. (b) Same heteroclinic connection in between the planes $W=2Z$ and
$Z=\zeta(W)$.

To capture the travelling wave with numerical simulation, we fix the temporal
parameters $\eta=0.1,$ $\alpha=0.07,\;\beta=0.2,\gamma=1.2,$ and consider the
following initial conditions:

$u(x,0)=\begin{cases}u_{*}&0\leq x<200\\\ u_{1}&200\leq x\leq
2000\end{cases},\quad\mbox{and}\quad v(x,0)=\begin{cases}v_{*}&0\leq x<200\\\
0&200\leq x\leq 2000\end{cases}.\;$ (33)

The eigenvalues of $J(\mathcal{E}_{*})$ determine the nature of the travelling
wave. If all the eigenvalues of the $J(\mathcal{E}_{*})$ are real, then the
system (12b) shows monotonic travelling wave. On the other hand, it exhibits
non-monotonic travelling wave if some of the eigenvalues of
$J(\mathcal{E}_{*})$ are complex. These complex eigenvalues are responsible
for the oscillation in the non-monotonic travelling wave-front. Snapshots of
monotonic and non-monotonic travelling waves at different times $t=100,200$
and $300$ are plotted in Fig. 6. The shape of the solution profiles are
similar for both the monotonic and non-monotonic cases with the advancement of
time. The minimum wave propagation speed $c_{min}=2d\sqrt{j_{3}}\approx 4.68$
for $\sigma=2.7$ and $d=46$. Comparing the locations of travelling wave at two
different times in Fig. 6(a), we calculate the speed of the travelling wave to
be $4.71$ approximately. Hence, the numerical value and analytical value of
the speed of propagation are in close agreement.

We have also plotted a diagram (see Fig. 7) in the $\sigma$-$c$ parametric
plane for the system (31), which shows the existence and nature of the
travelling wave solution of the system (12b). The system (31) does not admit
any travelling wave solution $\Phi(t,\textbf{p})$ for $c<2d\sqrt{j_{3}}.$ For
$c\geq 2d\sqrt{j_{3}}$, we obtain non-monotonic travelling wave when
$\sigma>2.68$ and monotonic travelling wave when $\sigma_{H}<\sigma<2.68.$

(a)

(b)

Figure 6: Travelling wave solution: (a) monotonic profile for $\sigma=2.7,$
(b) non-monotonic profile for $\sigma=1.9$. Other parameter values are
$d=46,\eta=0.1,\alpha=0.07,\beta=0.2$ and $\gamma=1.2.$
Figure 7: Existence and nature of travelling wave solution in the $\sigma$-$c$
parametric plane. Black dots indicate region of non-existence of travelling
wave solution. Red asterisks and magenta circles denote non-monotonic and
monotonic travelling waves respectively. Other parameter values are
$\eta=0.1,$ $\alpha=0.07,\;\beta=0.2,\;\gamma=1.2$ and $d=46$.

###  5.2 Spatio-temporal Chaos

The temporal model (7) shows bistability between trivial equilibrium point and
coexisting periodic solution for $\sigma_{Het}<\sigma<\sigma_{H}$. But,
introduction of a small amplitude spatial heterogeneity around the unstable
coexisting steady state leads to non-homogeneous non-stationary solution for
small values of $d$. We take parameter value $d=5$, $\sigma=1.8,$
$|\Omega|=500,$ and the corresponding results are shown in Fig. 8. The plot of
the spatial average of the prey population ($U_{av}$) against time shown in
Fig. 8(a) reveals the chaotic nature of the solution. This is further
supported by Fig. 8(b), which shows that the emerging pattern does not
converge to any stationary steady state. Using XPPAUT [64], we obtain the
largest Lyapunov exponent $\lambda_{max}=0.012>0,$ which confirms the chaotic
nature of the dynamics. Thus, the prey and predator populations exhibit
spatio-temporal chaos for certain parameter value in the Hopf region.

(a)

(b)

Figure 8:  Spatio-temporal chaos: (a) spatial average of the prey population
against time, (b) space-time plot of the prey population distribution after
the initial transient state. Parameter values are
$\alpha=0.07,\beta=0.2,\gamma=1.2,\eta=0.1$, $\sigma=1.8$ and $d=5.$

(a)

(b)

(c)

(d)

Figure 9: Multiple moving pule solution profiles at different times. Here
parameter values are $\sigma=0.6$, $d=20,\eta=0.1,\alpha=0.07,\beta=0.2$ and
$\gamma=1.2.$

###  5.3 Moving Pulse solution

Now we take $\sigma<\sigma_{Het}$ for which the co-existing steady state is
unstable. Suppose that the initial population distribution consists of co-
existing state in a small region around the centre of the domain. Introduction
of a small amplitude spatial heterogeneity leads to moving pulse solution
which shows intriguing complex dynamics. We fix the parameter set to
$\eta=0.1,\alpha=0.07,\gamma=1.2$, $\sigma=0.6$, $L=1000$, $\beta=0.2$ and
$d=20,$ and consider the following initial distribution of the populations:

$u(x,0)=\begin{cases}u_{*}+\xi(x)&495\leq x\leq 505\\\
0&\text{otherwise}\end{cases}\;\mbox{and}\;v(x,0)=\begin{cases}v_{*}+\xi(x)&495\leq
x\leq 505\\\ 0&\text{otherwise}\end{cases},\;$ (34)

where $\xi(x)$ is the Gaussian noise of amplitude $10^{-2}$. Initially, the
system (12b) produces a standing pulse solution consisting of a single island
of nonzero populations surrounded by dead zones. The populations in the dead
zones are almost extinct [see Fig. 9(a)]. As time increases, the width of the
island increases and splits into two islands as shown in Fig. 9(b). These new
islands are also moving towards the boundaries and they are symmetric about
the centre $x=L/2.$ These two islands again split giving rise to four islands
as shown in Fig. 9(c). Through collision between the adjacent islands, the
system either produces new islands or destroys some parent islands. The
maximum possible number of islands in the considered domain is $16$. After the
collision among the maximum number of islands, the number of islands reduces
to two as shown in Fig. 9(d) and the previous cycle repeats. The final spatio-
temporal dynamics of the system is periodic in time with time period $2558$
approximately. The space-time plot of the prey population is shown in Fig. 10,
which is a complex pattern consisting of many triangle shapes of invasion.
Such a pattern of many triangle shapes is often called the Sierpinski gasket
pattern [65, 66].

Figure 10: Space-time plot of the prey population for moving pulse solution.
Here parameter values are $\alpha=0.07,\beta=0.2,\gamma=1.2,\eta=0.1$,
$\sigma=0.6$ and $d=20.$

##  6 Conclusion

Our main objective of this work is to investigate the effect of the
reproductive Allee effect in prey growth on population establishment for prey-
predator type interaction model with specialist predator. Parametrization of
growth function and the choice of functional response are prominent in
determining factors behind the dynamics produced by the interacting
populations under consideration. In the absence of the Allee effect, the
models with logistic prey growth, prey-dependent functional response, and
linear mortality rate of the specialist predator can not provide any
indication of extinction of both species when we consider the temporal
dynamics [28]. The spatio-temporal extension of such models can capture the
population’s dynamic evolution over their habitats and predict localized
extinction. These spatio-temporal models admit travelling waves, waves of
invasion, and spatio-temporal chaos [23, 24]. On the contrary, mutual
interference among the predators can lead to localized establishment of both
species through the generation of stationary spatial patterns depending upon
Turing bifurcation. On the other hand, the temporal models with a strong Allee
effect in prey growth indicate the possibility of extinction of both species
either due to a global bifurcation or due to initial population density.
Uncontrolled or accelerated grazing rate of a specialist predator on the only
food source, which suffers from the Allee effect, drives the prey towards
localized or global extinction and hence system collapses. A spatio-temporal
prey-predator model with Allee effect in prey growth, prey-dependent
functional response, and linear mortality rate for specialist predator can
support the establishment of both species through successful invasion [67].
However, a model with predator-dependent functional response and/or intra-
specific competition among the predators supports stationary pattern formation
as well as oscillatory or chaotic coexistence depending upon the parameter
regime and the rates of diffusivity.

Here we have presented preliminary stability and bifurcation results for the
temporal model with the prey growth rate $\sigma$ as the bifurcation
parameter. As per model formulation, the prey growth rate incorporates the
mating success and it significantly contributes towards stable and oscillatory
coexistence of both species. A decline in mating success is not only harmful
to the specialist predators but also the initial abundance of the prey
population can not save them from extinction. A decline in the value of
$\sigma$ drives the system from stable steady-state coexistence to total
extinction through an oscillatory coexistence (super-critical Hopf
bifurcation) and the system collapses due to global heteroclinic bifurcation.
It is interesting to note that the existence of Hopf-bifurcation within
temporal setup is necessary for the Turing instability but not sufficient
[15]. RM type reaction kinetics supports supercritical Hopf bifurcation but
one can’t find the Turing pattern in a spatio-temporal model with RM reaction
kinetics.

The extinction scenario of populations depends not only on the Allee effect of
the prey population but also on the functional response which characterizes
any prey-predator type interaction. Appropriate parametrization of
Beddington–DeAngelis functional response allows us to compare the population
establishment mechanism with two different types of functional responses.
Beddington–DeAngelis functional response becomes ratio-dependent functional
response for $\alpha=0,$ which shows some qualitative changes in the system
dynamics for same set of other parameter values. Here, the destabilization of
coexisting equilibrium point leads to the eventual extinction of both
populations through a subcritical Hopf bifurcation. Contrary to the
Beddington–DeAngelis functional form, here global bifurcation (homoclinic)
doesn’t indicate the extinction of both populations rather it depends upon the
initial distribution of the populations for their survival. In this case, the
basin of attraction of the coexisting equilibrium point is bounded by a global
bifurcation generated (unstable) limit cycle. The evolution of the prey
population with the variation of prey growth rate is discussed in Appendix C.
However, neglecting predator interference, i.e., in case of Holling type-II
functional response, the mating success does not affect the stability of
coexisting equilibria (see more details in Appendix C).

Our systematic and through investigation of the spatio-temporal model reveals
that the emergence of stationary and non-stationary patterns solely linked
with the strength of species interaction, the rate of diffusivity and the
species’ abundance. By constructing upper and lower solutions, the asymptotic
behaviour of the solutions has been established under suitable parametric
restrictions. The solution of the diffusive system is always bounded
regardless of any parametric restriction, which proves the global existence of
the solution. If the maximum of the initial prey distribution is less than the
reproductive Allee threshold $u_{2},$ then both populations become extinct in
the final state, which is analogous to the temporal dynamics. Existence of
non-constant steady state of the diffusive system has been established using
Leray–Schauder degree theory. With the help of Poincare inequality, we have
also shown the absence of non-constant steady state solution when the ratio of
diffusivity lies in a particular range. However, the existence and non-
existence of heterogeneous solution somehow related to the domain size. Non-
constant steady-state does not exist in small domain size.

There are two prominent mechanisms behind the formation of non-constant steady
state patterns: the Turing instability and the BD transition. We have found
the Turing and BD transition thresholds by converting the steady-state
equations into a system of four coupled ordinary differential equations. A
variety of stable and unstable non-constant stationary solutions are found
through numerical continuation for admissible range of parameter values. These
solutions, which include various Turing mode solutions and localized
solutions, emerge through Turing instablity [17, 16] and BD points [19, 62].
We observe multi-stability between various non-constant solutions along with
the homogeneous steady-state solutions. The final stationary distribution of
the population at large time limit depends on the initial distribution of both
populations. Here we have shown that though the parameter setup satisfies the
Turing instability condition but the numerical simulation reveals the
settlement of population to localized pattern. In ecology, several prey
population has these types of patchy stationary distribution based on niche
separation with a continuous habitat [68, 69]. Prey and predator populations
exhibit identical patterns since patches with high prey density are favored by
predators. There is evidence of this phenomenon in the spatio-temporal
interactions of toxic newts and arrow snakes in western North America [70].

Introducing predator in a small domain leads to specific invasive pattern of
predators distributed over the entire domain which corresponds to travelling
wave solution. Generally, two types of travelling waves exist in prey-predator
models: one corresponds to a heteroclinic trajectory connecting two
homogeneous steady-states and the other corresponds to periodic travelling
wave that involve limit cycle surrounding an unstable homogeneous state. Our
model admits monotone and non-monotone travelling waves which are established
by showing a heteroclinic connection between two equilibria of the
corresponding system of four first-order ordinary differential equations. The
boundedness of the heteroclinic connection is established inside a wedge-
shaped region [52]. The parametric regions for the existence of monotonic and
non-monotonic travelling waves and non-existence of travelling wave have been
obtained with the help of exhaustive numerical simulations. We have also
validated the theoretical results with numerical simulations and discussed the
invasion profile over the parametric region. For a mobility rate of prey
individuals close to predator, the temporal extinction scenario alters to
localized extinction and regeneration of patches at nearby locations,
resulting in spatio-temporal chaotic or moving pulse solution. For lower
values of $d$ in the Hopf region and some restricted prey growth rate
$(\sigma_{Het}<\sigma<\sigma_{H})$, the spatio-temporal system exhibits time
aperiodic and non-homogeneous in space spatio-temporal chaos. Further decrease
in $\sigma$ gives rise to multiple moving pulse solution forming Sierpinski
gasket pattern. In real world, such types of spatio-temporal regular and
irregular oscillations have been reported for some prey-predator type
interactions that include interaction between Daphnia and Bythotrephes [71]
and interaction between tephritid flies and thistle population [72].

The novelty of this work lies in the identification of Turing bifurcation
along with the BD transition which leads to the formation of local patterns
through transient Turing-like patterns. Surprisingly, the Turing patterns
appeared to be transitory patterns and both the species get established at
large stationary patches through BD transition mechanism. This combination of
Turing instability and BD transition is an addition to the list of known
mechanisms behind long transients in spatio-temporal pattern formation [73].
The choice of parametrization for the Beddington-Deangelis functional response
helps us to conclude that this localized pattern formation scenario is
influenced by the reproductive Allee effect and obtained results can be
verified with other two functional responses namely Holling type-II and ratio-
dependent functional response. It is well-known that systems with a prey-
dependent functional response and linear death rate, do not produce stationary
Turing patterns. But if we choose a slightly different parameter setting for
which the temporal system with Holling type-II functional response has a
stable coexisting equilibrium, then the corresponding spatio-temporal model
supports the formation of localized stationary patterns due to BD transition.
A few localized patterns for different values of the parameter $d$ are shown
in Appendix C. Instead of the reproductive Allee effect, if we consider only
the logistic growth of the prey population, then the system fails to produce
any localized pattern. Thus, we conclude that the reproductive Allee effect
plays a central role in the formation of localized patterns.

## Declarations

Conflict of interest: The authors declare that they have no conflict of
interest.
Data availability statement: The authors declare that the manuscript has no
associated data.

##  Appendix A

Here we derive upper bounds for $I_{1}$, $I_{2}$ and $I_{3}$ used in Theorem
5. Using proposition 1, we find

$\displaystyle I_{1}=$
$\displaystyle\int_{\Omega}(u-\bar{u})^{2}\big{(}-(u^{2}+u\bar{u}+\bar{u}^{2})+(u+\bar{u})(u_{1}+u_{2})-u_{1}u_{2}\big{)}dx$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}(u-\bar{u})^{2}\left[(u+\bar{u})(u_{1}+u_{2})-u_{1}u_{2}\right]dx$
$\displaystyle\leq$ $\displaystyle
u_{1}(2-u_{2})\int_{\Omega}(u-\bar{u})^{2}dx,$
$\displaystyle\allowdisplaybreaks I_{2}=$
$\displaystyle\frac{1}{\alpha+\bar{u}+\beta\bar{v}}\int_{\Omega}(u-\bar{u})\left(\frac{\alpha(\bar{u}\bar{v}-uv)+\beta
v\bar{v}(\bar{u}-u)+u\bar{u}(\bar{v}-v)}{\alpha+u+\beta v}\right)dx$
$\displaystyle=$
$\displaystyle\frac{1}{\alpha+\bar{u}+\beta\bar{v}}\int_{\Omega}(u-\bar{u})\left(\frac{\bar{v}(\bar{u}-u)(\alpha+\beta
v)+u(\bar{v}-v)(\alpha+\bar{u})}{\alpha+u+\beta v}\right)dx$
$\displaystyle\leq$
$\displaystyle\frac{(\alpha+\bar{u})}{\alpha+\bar{u}+\beta\bar{v}}\int_{\Omega}(u-\bar{u})\frac{u(\bar{v}-v)}{\alpha+u+\beta
v}dx$ $\displaystyle\leq$
$\displaystyle\frac{u_{1}}{\alpha}\int_{\Omega}|(u-\bar{u})\,(v-\bar{v})|dx$
$\displaystyle\leq$
$\displaystyle\frac{u_{1}}{2\alpha}\int_{\Omega}(u-\bar{u})^{2}dx+\frac{u_{1}}{2\alpha}\int_{\Omega}(v-\bar{v})^{2}dx,$

and

$\displaystyle\allowdisplaybreaks I_{3}=$
$\displaystyle\int_{\Omega}(v-\bar{v})\bar{v}\left(\frac{\gamma
u}{\alpha+u+\beta v}-1\right)dx$ $\displaystyle=$
$\displaystyle\int_{\Omega}(v-\bar{v})\bar{v}\left(\frac{\gamma
u}{\alpha+u+\beta
v}-\frac{\gamma\bar{u}}{\alpha+\bar{u}+\beta\bar{v}}\right)dx$
$\displaystyle=$
$\displaystyle\frac{\gamma\bar{v}}{\alpha+\bar{u}+\beta\bar{v}}\int_{\Omega}(v-\bar{v})\frac{\alpha(u-\bar{u})+\beta(u\bar{v}-\bar{u}v)}{\alpha+u+\beta
v}\,dx$ $\displaystyle=$
$\displaystyle\frac{\gamma\bar{v}}{\alpha+\bar{u}+\beta\bar{v}}\int_{\Omega}(v-\bar{v})\frac{\beta
u(\bar{v}-v)+(u-\bar{u})(\alpha+\beta v)}{\alpha+u+\beta v}\,dx$
$\displaystyle\leq$
$\displaystyle\frac{\gamma\bar{v}}{\alpha+\bar{u}+\beta\bar{v}}\int_{\Omega}(v-\bar{v})\frac{(u-\bar{u})(\alpha+\beta
v)}{\alpha+u+\beta v}\,dx$ $\displaystyle\leq$
$\displaystyle\frac{\gamma}{\beta}\int_{\Omega}(v-\bar{v})(u-\bar{u})\,dx$
$\displaystyle\leq$
$\displaystyle\frac{\gamma}{2\beta}\int_{\Omega}(u-\bar{u})^{2}dx+\frac{\gamma}{2\beta}\int_{\Omega}(v-\bar{v})^{2}dx.$

##  Appendix B

Here we show 19-, 20- and 30- mode Turing solutions in Fig. 11 and few
localized solutions in Fig. 12.

(a)

(b)

(c)

Figure 11: Turing pattern solution of various modes: (a) 19- mode solution for
$\sigma=1.939$, (b) 20- mode solution for $\sigma=1.867$, (c) 21- mode
solution for $\sigma=1.812$. The color of a solution curve corresponds to the
same color diamond point shown in Fig. 3(b). Other parameter values are
$\alpha=0.07,\beta=0.2,\gamma=1.2,\eta=0.1$ and $d=46.$

(a)

(b)

(c)

Figure 12: Localized pattern solutions where the colour of a solution curve
corresponds to the same color point shown in Fig. 3(c). Other parameter values
are $\alpha=0.07,\beta=0.2,\gamma=1.2,\sigma=2.36,\eta=0.1$ and $d=46.$

##  Appendix C

(a)

(b)

Figure 13: (a) Bifurcation diagram in the $\sigma$-$u$ plane for the ratio-
dependent functional response. (b) A zoom version of Fig. (a) for
$\sigma\in[3.9,3.98]$. The point marked SN denotes the saddle-node bifurcation
threshold. Here, solid green, red dashed and black dashed curves represent
stable, unstable and saddle branches of equilibria respectively. Further,
solid red color curve represents the maximum and minimum of $u$ of the
unstable limit cycle. Other parameter values are
$\alpha=0,\beta=0.2,\gamma=1.2$ and $\eta=0.1.$

Here we compare our findings with other functional responses. If we assume
$\alpha=0$, then the resulting functional response is called a ratio-dependent
functional response. The temporal dynamics of the system, keeping all other
parameters the same with $\alpha=0$, changes significantly. Here, two
coexisting equilibria $E_{*}^{1}$ and $E_{*}^{2}$ emerge through a saddle-node
bifurcation, where $E_{*}^{2}$ is always a saddle-node and the stability of
$E_{*}^{1}$ depends on the Hopf bifurcation threshold $\sigma_{H}.$ An
unstable limit cycle is generated due to subcritical Hopf bifurcation at
$\sigma=\sigma_{H}$ which vanishes due to a global homoclinic bifurcation. The
temporal dynamics is summarized in the one-parameter bifurcation diagram shown
in Fig. 13. The corresponding spatio-temporal model shows stationary Turing
and localized patterns. It also exhibits dynamic patterns that include
multiple moving pulse solution, spatio-temporal chaos, and traveling wave.

However, the chosen functional response becomes Holling type-II functional
response for $\beta=0$. Then the system can have at most one coexisting
equilibrium point $E_{*}$, whose feasibility comes from a transcritical
bifurcation. However, the stability of $E_{*}$ is independent of the parameter
$\sigma.$ If we take the parameter $\gamma=1.11,$ then the unique coexisting
equilibrium point $E_{*}$ is asymptotically stable (whenever it exists) for
all value of $\sigma$. Although the system does not produce Turing pattern but
it exhibits localized pattern. We have shown few localized patterns for
different values of diffusion parameter $d$ in Fig. 14. Apart from the
stationary pattern, the system also shows dynamic patterns that include
multiple moving pulse solution, spatio-temporal chaos and travelling wave
solution.

(a)

(b)

(c)

Figure 14: Localized patterns for the system (12b) with Holling type-II
functional response: (a) $d=90$ (b)$d=100$ (c) $d=110$. Other parameter values
are $\sigma=2,\gamma=1.11,\eta=0.1,\alpha=0.07$ and $\beta=0.$

## References

  * [1] AJ Lotka.  Undamped oscillations derived from the law of mass action.  Journal of the American Chemical Society, 42(8):1595–1599, 1920\.
  * [2] V Volterra.  Variazioni e fluttuazioni del numero d’individui in specie animali conviventi.  C. Ferrari, 1926.
  * [3] PA Stephens, WJ Sutherland, and RP Freckleton.  What is the allee effect?  Oikos, pages 185–190, 1999.
  * [4] E Venturino and S Petrovskii.  Spatiotemporal behavior of a prey–predator system with a group defense for prey.  Ecological Complexity, 14:37–47, 2013.
  * [5] W Wang, Y Takeuchi, Y Saito, and S Nakaoka.  Prey–predator system with parental care for predators.  Journal of theoretical biology, 241(3):451–458, 2006.
  * [6] MT Alves and F Hilker.  Hunting cooperation and allee effects in predators.  Journal of theoretical biology, 419:13–22, 2017.
  * [7] Y Kang and L Wedekin.  Dynamics of a intraguild predation model with generalist or specialist predator.  Journal of mathematical biology, 67(5):1227–1259, 2013.
  * [8] S Kondo and R Asai.  A reaction–diffusion wave on the skin of the marine angelfish pomacanthus.  Nature, 376:765–768, 1995.
  * [9] A Ducrots and M Langlais.  A singular reaction–diffusion system modelling prey–predator interactions: Invasion and co-extinction waves.  Journal of Differential Equations, 253(2):502–532, 2012.
  * [10] CA Klausmeier.  Regular and irregular patterns in semiarid vegetation.  Science, 284:1826–1828, 1999.
  * [11] RA Cangelosi, DJ Wollkind, BJ Kealy-Dichone, and I Chaiya.  Nonlinear stability analyses of turing patterns for a mussel-algae model.  Journal of Mathematical Biology, 70(6):1249–1294, 2015.
  * [12] AM Turing.  The chemical basis of morphogenesis.  Phil. Trans. Royal Society, 237:37–72, 1952.
  * [13] LA Segel and JL Jackson.  Dissipative structure: an explanation and an ecological example.  J. Theor. Biol., 37(3):545–59, 1972.
  * [14] MC Cross and PC Hohenberg.  Pattern formation outside of equilibrium.  Reviews of modern physics, 65(3):851, 1993.
  * [15] M Banerjee and S Petrovskii.  Self-organised spatial patterns and chaos in a ratio-dependent predator–prey system.  Theoretical Ecology, pages 37–53, 2011.
  * [16] T Hillen.  A turing model with correlated random walk.  Journal of Mathematical Biology, 35(1):49–72, 1996.
  * [17] S Dey, M Banerjee, and S Ghorai.  Analytical detection of stationary turing pattern in a predator-prey system with generalist predator.  Mathematical Modelling of Natural Phenomena, 17:33, 2022.
  * [18] LA Belyakov, L Yu Glebsky, and LM Lerman.  Abundance of stable stationary localized solutions to the generalized 1d swift-hohenberg equation.  Computers & Mathematics with Applications, 34(2-4):253–266, 1997\.
  * [19] AR Champneys.  Homoclinic orbits in reversible systems and their applications in mechanics, fluids and optics.  Physica D: Nonlinear Phenomena, 112(1-2):158–186, 1998.
  * [20] PD Woods and AR Champneys.  Heteroclinic tangles and homoclinic snaking in the unfolding of a degenerate reversible hamiltonian–hopf bifurcation.  Physica D: Nonlinear Phenomena, 129(3-4):147–170, 1999.
  * [21] J Burke and E Knobloch.  Snakes and ladders: localized states in the swift–hohenberg equation.  Physics Letters A, 360(6):681–688, 2007.
  * [22] F Al Saadi, A Champneys, C Gai, and T Kolokolnikov.  Spikes and localised patterns for a novel schnakenberg model in the semi-strong interaction regime.  European Journal of Applied Mathematics, 33(1):133–152, 2022.
  * [23] SV Petrovskii and H Malchow.  A minimal model of pattern formation in a prey-predator system.  Mathematical and Computer Modelling, 29(8):49–63, 1999.
  * [24] M Banerjee and V Volpert.  Spatio-temporal pattern formation in rosenzweig-macarthur model: effect of nonlocal interactions.  Ecological complexity, 30:2–10, 2017.
  * [25] M Mimura and JD Murray.  On a diffusive prey-predator model which exhibits patchiness.  Journal of Theoretical Biology, 75(3):249–262, 1978.
  * [26] X Wang, Y Cai, and H Ma.  Dynamics of a diffusive predator-prey model with allee effect on predator.  Discrete Dynamics in Nature and Society, 2013, 2013.
  * [27] WC Allee.  Animal aggregations: A study in general sociology.  University of Chicago Press, Chicago, 1931.
  * [28] F Courchamp, L Berec, and J Gascoigne.  Allee effects in ecology and conservation.  OUP Oxford, 2008.
  * [29] DS Boukal and L Berec.  Single-species models of the allee effect: extinction boundaries, sex ratios and mate encounters.  Journal of Theoretical Biology, 218(3):375–394, 2002.
  * [30] JC Gascoigne and RN Lipcius.  Allee effects driven by predation.  Journal of Applied Ecology, 41(5):801–810, 2004.
  * [31] M Kuussaari, I Saccheri, M Camara, and I Hanski.  Allee effect and population dynamics in the glanville fritillary butterfly.  Oikos, pages 384–392, 1998.
  * [32] FJ Rohlf.  The effect of clumped distributions in sparse populations.  Ecology, 50(4):716–721, 1969.
  * [33] PK Molnar, AE Derocher, MA Lewis, and MK Taylor.  Modelling the mating system of polar bears: a mechanistic approach to the allee effect.  Proceedings of the Royal Society B: Biological Sciences, 275(1631):217–226, 2008.
  * [34] AW Stoner, MH Davis, and CJ Booker.  Negative consequences of allee effect are compounded by fishing pressure: comparison of queen conch reproduction in fishing grounds and a marine protected area.  Bulletin of Marine Science, 88(1):89–104, 2012.
  * [35] AM Kramer, B Dennis, AM Liebhold, and JM Drake.  The evidence for allee effects.  Population Ecology, 51(3):341–354, 2009.
  * [36] L Berec, E Angulo, and F Courchamp.  Multiple allee effects and population management.  Trends in Ecology & Evolution, 22(4):185–191, 2007.
  * [37] P Aguirre, E González-Olivares, and E Sáez.  Three limit cycles in a leslie–gower predator-prey model with additive allee effect.  SIAM Journal on Applied Mathematics, 69(5):1244–1262, 2009.
  * [38] MA Lewis and P Kareiva.  Allee dynamics and the spread of invading organisms.  Theoretical Population Biology, 43(2):141–158, 1993.
  * [39] B Dennis.  Allee effects: population growth, critical density, and the chance of extinction.  Natural Resource Modeling, 3(4):481–538, 1989.
  * [40] S Petrovskii, R Blackshaw, and B Li.  Consequences of the allee effect and intraspecific competition on population persistence under adverse environmental conditions.  Bulletin of Mathematical Biology, 70(2):412–437, 2008.
  * [41] M Jankovic and S Petrovskii.  Are time delays always destabilizing? revisiting the role of time delays and the allee effect.  Theoretical ecology, 7(4):335–349, 2014.
  * [42] N Mukherjee and V Volpert.  Bifurcation scenario of turing patterns in prey-predator model with nonlocal consumption in the prey dynamics.  Communications in Nonlinear Science and Numerical Simulation, 96:105677, 2021.
  * [43] AJ Terry.  Predator–prey models with component allee effect for predator reproduction.  Journal of mathematical biology, 71(6):1325–1352, 2015.
  * [44] GT Skalski and JF Gilliam.  Functional responses with predator interference: viable alternatives to the holling type ii model.  Ecology, 82(11):3083–3092, 2001.
  * [45] L Přibylová and L Berec.  Predator interference and stability of predator–prey dynamics.  Journal of mathematical biology, 71(2):301–323, 2015.
  * [46] J R Beddington.  Mutual interference between parasites or predators and its effect on searching efficiency.  The Journal of Animal Ecology, pages 331–340, 1975.
  * [47] D Alonso, F Bartumeus, and J Catalan.  Mutual interference between predators can give rise to turing spatial patterns.  Ecology, 83(1):28–34, 2002.
  * [48] C Skov, BB Chapman, H Baktoft, J Brodersen, C Brönmark, LA Hansson, K Hulthén, and PA Nilsson.  Migration confers survival benefits against avian predators for partially migratory freshwater fish.  Biology letters, 9(2):20121178, 2013.
  * [49] J Fryxell and P Lundberg.  Individual behavior and community dynamics, volume 20.  Springer Science & Business Media, 2012.
  * [50] M Banerjee and V Volpert.  Prey-predator model with a nonlocal consumption of prey.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(8):083120, 2016.
  * [51] SR Dunbar.  Travelling wave solutions of diffusive lotka-volterra equations.  Journal of Mathematical Biology, 17(1):11–32, 1983.
  * [52] W Huang.  A geometric approach in the study of traveling waves for some classes of non-monotone reaction–diffusion systems.  Journal of Differential Equations, 260(3):2190–2224, 2016.
  * [53] S Ai, Y Du, and R Peng.  Traveling waves for a generalized holling–tanner predator–prey model.  Journal of Differential Equations, 263(11):7782–7814, 2017.
  * [54] M Banerjee, N Mukherjee, and V Volpert.  Prey-predator model with nonlocal and global consumption in the prey dynamics.  Discrete & Continuous Dynamical Systems-S, 13(8):2109, 2020.
  * [55] L Perko.  Differential Equations and Dynamical Systems.  Springer-Verlag, New York, 2000.
  * [56] P Liu, J Shi, and Y Wang.  Bifurcation from a degenerate simple eigenvalue.  Journal of Functional Analysis, 264(10):2269–2299, 2013.
  * [57] CV Pao.  Nonlinear parabolic and elliptic equations.  Springer Science & Business Media, 2012.
  * [58] J Smoller.  Shock waves and reaction—diffusion equations, volume 258.  Springer Science & Business Media, 2012.
  * [59] JK Hale.  Asymptotic behavior of dissipative systems.  American Mathematical Soc., 2010.
  * [60] Y Lou and WM Ni.  Diffusion, self-diffusion and cross-diffusion.  Journal of Differential Equations, 131(1):79–131, 1996.
  * [61] PY Pang and M Wang.  Strategy and stationary pattern in a three-species predator–prey model.  Journal of Differential Equations, 200(2):245–273, 2004.
  * [62] F Al Saadi and A Champneys.  Unified framework for localized patterns in reaction–diffusion systems; the gray–scott and gierer–meinhardt cases.  Philosophical Transactions of the Royal Society A, 379(2213):20200277, 2021.
  * [63] H Uecker, D Wetzel, and JD Rademacher.  pde2path-a matlab package for continuation and bifurcation in 2d elliptic systems.  Numerical Mathematics: Theory, Methods and Applications, 7(1):58–106, 2014.
  * [64] B Ermentrout and A Mahajan.  Simulating, analyzing, and animating dynamical systems: a guide to xppaut for researchers and students.  Appl. Mech. Rev., 56(4):B53–B53, 2003.
  * [65] Y Hayase and T Ohta.  Self-replicating pulses and sierpinski gaskets in excitable media.  Physical Review E, 62(5):5998, 2000.
  * [66] VB Kazantsev, VI Nekorkin, S Binczak, and JM Bilbault.  Spiking patterns emerging from wave instabilities in a one-dimensional neural lattice.  Physical Review E, 68(1):017201, 2003.
  * [67] A Morozov, S Petrovskii, and BL Li.  Spatiotemporal complexity of patchy invasion in a predator-prey system with the allee effect.  Journal of theoretical Biology, 238(1):18–35, 2006.
  * [68] RT Holmes and JC Schultz.  Food availability for forest birds: effects of prey distribution and abundance on bird foraging.  Canadian journal of Zoology, 66(3):720–728, 1988.
  * [69] JM Fryxell, ARE Sinclair, and G Caughley.  Wildlife ecology, conservation, and management.  John Wiley & Sons, 2014.
  * [70] Edmund D Brodie J, BJ Ridenhour, and ED III Brodie.  The evolutionary response of predators to dangerous prey: hotspots and coldspots in the geographic mosaic of coevolution between garter snakes and newts.  Evolution, 56(10):2067–2082, 2002.
  * [71] JT Lehman and CE Cáceres.  Food-web responses to species invasion by a predatory invertebrate: Bythotrephes in lake michigan.  Limnology and Oceanography, 38(4):879–891, 1993.
  * [72] F Jeltsch, C Wissel, S Eber, and R Brandl.  Oscillating dispersal patterns of tephritid fly populations.  Ecological modelling, 60(1):63–75, 1992.
  * [73] A Morozov, K Abbott, K Cuddington, T Francis, G Gellner, A Hastings, YC Lai, S Petrovskii, K Scranton, and ML Zeeman.  Long transients in ecology: Theory and applications.  Physics of Life Reviews, 32:1–40, 2020.