Oscillations and pattern formation in

a slow-fast prey-predator system

Pranali Roy Chowdhury1, Sergei Petrovskii2,3, Malay Banerjee1

1 Indian Institute of Technology Kanpur, Kanpur - 208016, India

2 School of Mathematics & Actuarial Science, University of Leicester,
Leicester LE1 7RH, UK

3 Peoples Friendship University of Russia (RUDN University), 6 Miklukho-
Maklaya St,

Moscow 117198, Russian Federation

###### Abstract

We consider the properties of a slow-fast prey-predator system in time and
space. We first argue that the simplicity of prey-predator system is apparent
rather than real and there are still many of its hidden properties that have
been poorly studied or overlooked altogether. We further focus on the case
where, in the slow-fast system, the prey growth is affected by a weak Allee
effect. We first consider this system in the non-spatial case and make its
comprehensive study using a variety of mathematical techniques. In particular,
we show that the interplay between the Allee effect and the existence of
multiple timescales may lead to a regime shift where small-amplitude
oscillations in the population abundances abruptly change to large-amplitude
oscillations. We then consider the spatially explicit slow-fast prey-predator
system and reveal the effect of different time scales on the pattern
formation. We show that a decrease in the timescale ratio may lead to another
regime shift where the spatiotemporal pattern becomes spatially correlated
leading to large-amplitude oscillations in spatially average population
densities and potential species extinction.

Keywords: Slow-fast time scale; relaxation oscillation; canard cycle; spatial
pattern; regime shift

##  1 Introduction

In the natural environment, interactions in a population community are usually
quite complex [5, 25, 77]. This ubiquitous complexity have several different
sources such as the complexity of the wood web, nonlinearity of species
feedbacks, multiplicity of temporal and spatial scales, etc. It is extremely
difficult, in fact hardly possible at all to capture the entire complexity of
ecological interactions in a single mathematical model or framework. Instead,
the usual means of analysis tend to focus on a particular aspect or feature of
the ecological system. For instance, while the food web theory endeavours to
link the properties of a realistic population community to the complexity of
the corresponding food web, in particular by analysing the web connectivity
and revealing the bottlenecks [1, 60], a lot of attention focuses on the
properties of simpler ‘building blocks’ from which the web is made [24]. A
variety of blocks of intermediate complexity have been considered, a few
examples are given by the three-species competition system [17], intraguild
predation [19] and a three-species resource-consumer-predator food chain [15].

Arguably, the most basic block is the prey-predator system. It has been a
focus of research for almost a century [64, 80] and there is a tendency to
think about it as a fully studied, textbook material [36]. However, this is
far from true. The apparent mathematical simplicity of the prey-predator
system (usually associated with the classical Rosenzweig-McCarthur model as a
paradigm [64, 77]) is superficial rather than real, and there has recently
been a surge of interest and an increase in mathematical modellng literature
dealing with its ‘hidden’, overlooked properties, with more than a hundred of
papers published in the first quarter of 2021 alone111Data are taken from the
Web of Science. New properties readily arise as soon as one introduces
relatively small (i.e. preserving the defining structure of the model),
biologically motivated changes into the paradigmatic system, e.g. adding
explicit heterogeneous space [83], changing the specialist predator to a
generalist one [62, 69], changing the properties of predator’s functional
response [2, 20], considering different types of density dependence in the
population growth or mortality [10, 23], or taking into account the fact that
the intraspecific dynamics of prey and predator often occur on a very
different time scale [26, 59].

While one of the generic properties of a prey-predator system is its intrinsic
capability to produce sustained population cycles (due to the emergence of a
stable limit cycle in a certain parameter range [37]), with many fundamental
implications for the population dynamics, another equally important property
is its capacity to exhibit pattern formation, in particualar due to the Turing
instability [68, 78]. The latter has been a focus of many groundbreaking
studies that linked the patterns observed in various biological and ecological
systems to the dissipative instability in a prey-predator (or, more
generically, activator-inhibitor) system, e.g. see [13, 14, 22, 38, 40, 44,
45, 46, 47, 48, 49, 63], also [50] for an exhaustive review of earlier
research. Other studies also discovered and considered in detail a possibility
of non-Turing pattern formation, in particular due to the interplay between
the Hopf bifurcation and diffusion [51, 52, 53, 57, 67, 73] as well as pattern
formation resulting from the Turing-Hopf bifurcation [3].

Interestingly, in spite of the large number of modelling papers concerned with
the prey-predator system, there are still a number of issues poorly
investigated. One such issue is the interaction between different types of
density dependence and the existence of different time scales, either in a
spatial or nonspatial system. Indeed, one context where the prey-predation
framework has been particularly successful to provide a new insight into the
mechanisms of ecological interactions are is large-magnitude nearly-periodical
fluctuations in population size that has been observed in many species and
ecosystems. In such a case, typically, a large outbreak in population
abundance is followed by a population decline, often to a small population
size or density. For instance, the fluctuations in populations of snowshoe
hares and Canadian lynx in the Canadian Boreal forest was modeled with the
help of a tri-trophic food web model [75] where the population explosion of
the lynx was observed every 9-11 years followed by a rapid decline in the
population of hares. Also in case of plankton ecosystem within lake, the
seasonal abundance of zooplankton (particularly Daphnia) is frequently
observed, which completely grazes down the algal biomass thus resulting in
clear-water phases in lakes [67]. The exact causes of these fluctuations are
still a debatable issue among various researchers. However, one of the common
trait observed in the above examples is that the bottom level of a multi-
trophic system or the basal prey has faster growth and decay compared to their
consumers. On the other hand, the population of the budworm can increase
several hundred fold within a span of few years whereas the leaves of adult
trees do not grow at a comparable rate. This resulted in the outbreak of
spruce budworm which destroyed the balsam forest of eastern Northern America
[34]. To capture this type of rapid growth/decay for the interacting species,
researchers introduced the mathematical models with slow-fast time scale. A
small time scale parameter is introduced either in the prey growth equation or
in the predator growth depending upon the species under consideration.

In a rather general case, the prey-predator interaction in a nonspatial system
can be modeled by a system of coupled ordinary differential equations

$\displaystyle u^{\prime}$ $\displaystyle=uf(u)-vg(u,v),$ (1) $\displaystyle
v^{\prime}$ $\displaystyle=evg(u,v)-m(v)v,$

where $u$ and $v$ denote the prey and predator densities, respectively, at
time $t$. (A spatially explicit approach includes diffusion terms, hence
turning the ODEs to PDEs, see Section 5). Here both the species are assumed to
be distributed homogeneously within their habitat. The function $f(u)$
represents the per capita growth rate of prey, $g(u,v)$ describes the prey-
predator interaction and $\mu$ describes the natural death rate of predators
in absence of prey, $e$ is known as the conversion efficiency. Assuming the
rate of growth of prey population much faster than its predator, a time-scale
parameter $0<\varepsilon\ll 1$ is introduced which transform the original
model to a slow-fast model as follows

$\displaystyle\varepsilon u^{\prime}$ $\displaystyle=uf(u)-vg(u,v),$ (2)
$\displaystyle v^{\prime}$ $\displaystyle=evg(u,v)-m(v)v.$

This type of slow-fast prey predator models were first studied by Rinaldi and
Muratori [61], so far as our knowledge goes, where the cyclic coexistence of
the slow-fast limit cycle was discussed. They also analyzed the cyclic
fluctuation in population densities of three species model in a slow-fast
setting with one and two multiple time scale parameters.

In mathematical literature, the slow-fast systems are considered as singularly
perturbed ordinary differential equation, where $\varepsilon$ is the singular
perturbation parameter. The standard stability and bifurcation analysis
performed for the prey-predator models was not enough to analyze the complete
dynamics exhibited by the slow-fast systems. Many mathematical techniques were
developed to study this class of systems. In the late 1970s, Neil Fenichel
[11] introduced a geometric approach based on the invariant manifold theory to
study the singularly perturbed coupled systems, known as Geometric Singular
Perturbation Theory (GSPT). Using this theory the dynamics of the full slow-
fast system are studied by reducing it to sub-systems of lower dimension and
thereby studying the complete dynamics of the subsystems. The application of
Fenichel’s theory in the context of biology was well explained by G. Hek [16].
But this theory fails to approximate the dynamics near the non-hyperbolic
equilibrium points where the system encounters a singularity. Later in 2001,
Krupa and Szmolyan [27, 28] extended Fenichel’s theory to overcome the
difficulty around non-hyperbolic points using the blow-up technique. This was
based on the pioneering work of Dumortier [7, 8, 9]. The main idea behind this
was to blow up the non-hyperbolic equilibrium points of the system by a
4-dimensional unit sphere $S^{3}$ and the trajectories of the blow-up system
is mapped on and around the sphere. In case, the blow-up space still has non-
hyperbolic points, sequence of blow-up maps can be used to desingularize the
system.

Before the development of mathematical tools to study this class of systems, a
Dutch Physicist Van der Pol [9, 79] observed large amplitude periodic
oscillation consisting of slow and fast dynamics, which he named relaxation
oscillation. These are periodic solutions consisting of slow curvilinear
motion and sudden fast jumps. These types of slow-fast limit cycles were later
observed in many chemical and biological systems [26, 42, 61, 81, 82]. Another
type of periodic solution observed in singularly perturbed systems are canard
solutions. This was first investigated by E. Benoît et. al [4] while studying
the Van der Pol Oscillator. Dumortier and Roussarie, in their seminal work
[9], analyzed this phenomenon through a geometric approach, using blow-up
technique and with the help of invariant manifold theory. A canard is a
solution of a singularly perturbed system which follows an attracting slow
manifold, closely passing through the bifurcation point of the critical
manifold, and then following a repelling slow manifold for $\mathcal{O}(1)$
time. It was observed that for the existence of canard solution Hopf
bifurcation is necessary [28]. The fast transition from small stable limit
cycles appearing through Hopf bifurcation to large amplitude relaxation
oscillation via a sequence of canard cycles within an exponentially small
range of the parameter is known as canard explosion. In real-world ecosystems
this phenomenon can be related to sudden outbreak or decline of a particular
species [34, 66, 67, 74, 75].

Over the last few years, several works have been done on prey-predator systems
with slow-fast time scale. In [42, 43, 61], the authors have analyzed the
periodic bursting of high and low-frequency oscillations in interacting
population models with two and three-trophic level with slow-fast time scale.
A novel 1-fast-3-slow dynamical system have been developed in [58] to consider
the adaptive change of diet of a predator population that switches its feeding
between two prey populations. The classical Rosenzweig–MacArthur (RM) model
and the Mass Balance chemostat model in the slow-fast setting is studied in
[26], where the authors have shown that the RM model exhibits canard explosion
in the oscillatory regime of the parameter space whereas the later model does
not exhibit such phenomenon. They have used the asymptotic expansion technique
to determine the canard explosion point. In [59], the authors have used the
blow-up technique to obtain an analytical expression of the bifurcation
thresholds for which maximal canard solution occurs in the RM-model. The
existence and uniqueness of the relaxation oscillation cycle have been studied
for the Leslie-Gower model with the help of entry-exit function and GSPT in
[82]. The rich and complex slow-fast dynamics of the predator-prey model with
Beddington-DeAngelis functional response is studied in [65]. To the best of
our knowledge, there is no work so far in literature considering the Allee
effect in the slow-fast prey-predator model. In this paper, first we consider
the slow-fast dynamics of the classical Rosenzweig–MacArthur (RM) model with
multiplicative weak Allee effect in prey growth equation using GSPT and blow-
up technique.

In population ecology, the Allee effect is a widely observed phenomenon
especially at low population density, which describes a positive relationship
between species population and per capita population growth rate of species.
The main causes of Allee effect include difficulties in mate finding,
inbreeding depression, cooperative defense mechanism etc. Mostly, we are
concerned about the demographic Allee effect which can be classified as:
strong Allee and weak Allee effect. For strong Allee effect, the per capita
growth rate is negative below some critical population density (Allee
threshold), and the growth rate becomes positive above that threshold.
Whereas, in case of weak Allee effect, per capita growth rate is small and
remains positive even at low population densities. But with the introduction
of the time scale parameter the per capita growth rate becomes much higher
even at low density. Thus the species can recover itself from the endemic
level, and extinction is prevented. In this paper, we have incorporated the
weak Allee effect in prey’s growth in order to capture the true essence of the
slow-fast cycle.

The main objective of this paper is to provide a detailed slow-fast analysis
of the temporal model based on the various mathematical approach discussed
above, and numerically investigating the corresponding spatially extended
slow-fast model. In prey-predator models, the oscillatory dynamics of the
system arises from the Hopf bifurcation but in the slow-fast setting other
than Hopf bifurcating limit cycle, the system exhibit various other
interesting periodic solution namely canard and relaxation oscillation. Here
we are interested to explore these solutions analytically and numerically with
the help of sophisticated slow-fast techniques as discussed above.

Over the last few decades, significant work has been done to study the
mechanism of spatial dispersal of species, with the help of reaction-diffusion
systems. In this regard, the study of invasion of the exotic species emerged
to be of particular interest for many theoretical or field ecologists.
Biological invasion is a complex phenomenon that starts with a local
introduction of exotic species and once it gets established in a particular
region, they start spreading and occupying new areas [33]. The study of rate
and pattern of spread is of primary importance as they have a huge
environmental impact and also can be used as a control measure for other
species. The invasion of exotic species takes place via propagation of
continuous wave fronts, as well as via irregular movement of separate
population patches [32]. In [41, 56], the authors have shown that patchy
invasion is possible in deterministic models as a result of the Allee effect.
To the best of our knowledge, there exist hardly any work in literature which
such complete analysis of the slow-fast prey-predator model as well as on the
effect of explicit time-scale parameters on pattern formation. Here, first we
perform exhaustive numerical simulations to examine the pattern spread and
patchy invasion of the species. And then, we examine how the invasion of the
species are affected by the time-scale parameter.

This paper is divided into two parts, in the first part we will provide
mathematical analysis of the temporal slow-fast model and the second part is
supported by exhaustive numerical simulations to investigate the effect of
varying time-scale in biological invasion. In section 2, we introduce the non-
dimensionalised temporal model and standard stability analysis is performed.
Then in section 3 we introduced the slow-fast system. In section 4 we
discussed GSPT and blow-up technique for a detailed mathematical analysis of
slow-fast systems. The existence and uniqueness of the relaxation oscillation
is studied here followed by the phenomenon of canard explosion. In section 5,
we consider the corresponding slow-fast spatio-temporal model to examine how
the spread of invasive species is affected by time-scale parameters. Finally,
we draw the conclusion of our work in section 6.

##  2 Temporal Model and its linear stability analysis

We consider the classical Rosenzweig-MacArthur prey-predator model with the
multiplicative weak Allee effect in prey growth [6, 50, 70]. Let $u$ and $v$
be the prey and its specialist predator densities, respectively, at time $t$.
In appropriately chosen dimensionless variables and parameter (see [41] for
details), the model is given by the following equations:

$\displaystyle\dfrac{du}{dt}$ $\displaystyle=$ $\displaystyle f(u,v):=\gamma
u(1-u)(u+\beta)-v\Big{(}\frac{u}{1+\alpha u}\Big{)},$ (3a)
$\displaystyle\dfrac{dv}{dt}$ $\displaystyle=$ $\displaystyle
g(u,v):=v\Big{(}\frac{u}{1+\alpha u}-\delta\Big{)}.$ (3b)

Here and below, the sign “:=” means “is defined”. We focus on the case where
the growth rate of the prey population is damped by the weak Allee effect, so
that $0<\beta<1$. For $\beta<0$, the Allee effect becomes strong (in this
case, the prey population has another [usntable] equilibrium at $u=\beta$);
for $\beta\geq 1$, the Allee effect is absent [31]. The per capita growth rate
$f(u,v)/u$ is increasing for $0<u<\dfrac{1-\beta}{2}$ and decreasing for
$\dfrac{1-\beta}{2}<u<1$. The predator is a specialist predator as they do not
have any alternative food source to survive apart from $u$. The prey-dependent
functional response is taken to be Holling type II [18]. The system contains
four positive dimensionless parameters where $\beta$ quantifies the weak Allee
parameter, $\gamma$ is the coefficient proportional to the maximum per capita
growth rate, called characteristic growth rate [21]. The parameter $\alpha$
characterizes the inverse saturation level of the functional response and
$\delta$ is the natural mortality rate of the predator. Throughout this paper
we will consider $\delta$ as the bifurcation parameter to determine the
stability conditions of the coexisting steady-state for the model 3.

Depending on the species traits, the prey population often grows much faster
than its predator; one well known example is given by hare and lynx where
hares reproduce much faster than lynx [75]. This motivated researchers to
introduce a small time-scale parameter $\varepsilon,\ 0<\varepsilon<1$ in the
basic model (3). The parameter $\varepsilon$ is interpreted as the ratio
between the linear death rate of the predator and the linear growth rate of
the prey [16, 61]. And the assumption $\varepsilon<1$ implies that one
generation of predator can encounter several generations of prey [18, 29].
Therefore considering the difference in the time scale, the slow-fast version
of the dimensionless model (3) can be written as

$\displaystyle\dfrac{du}{dt}$ $\displaystyle=$ $\displaystyle f(u,v)=\gamma
u(1-u)(u+\beta)-\frac{uv}{1+\alpha u},$ (4a) $\displaystyle\dfrac{dv}{dt}$
$\displaystyle=$ $\displaystyle\varepsilon g(u,v)=\varepsilon
v\Big{(}\frac{u}{1+\alpha u}-\delta\Big{)},$ (4b)

with initial conditions $u(0)\geq 0,\ v(0)\geq 0$ . Since the prey population
grows faster compared to the predator, $u$ and $v$ are referred to as fast and
slow variables, respectively and time $t$ is called fast time. The equilibrium
points for the system are independent of $\varepsilon$, thus system (3) and
(4) has same equilibrium points. The extinction equilibrium point and prey
only equilibrium point of system (3) (as well as for (4) are given by
$E_{0}=(0,0)$ and $E_{1}=(1,0)$ respectively. The interior equilibrium point
$E_{*}(u_{*},v_{*})$ of the system is the point where the non-trivial prey
nullcline intersect with non-trivial predator nullcline in the interior of the
positive quadrant, and we have,

$u_{*}=\dfrac{\delta}{1-\alpha\delta},\ \
v_{*}=\gamma(1-u_{*})(u_{*}+\beta)(1+\alpha u_{*}).$

$E_{*}$ is feasible if the parametric restriction $\delta(\alpha+1)<1$ holds.
With the help of linear stability analysis, we find $E_{0}$ is always a saddle
point. $E_{1}$ is stable for $\delta>\dfrac{1}{1+\alpha}$ and saddle point for
$\delta<\dfrac{1}{1+\alpha}$. $E_{*}$ bifurcates from predator free
equilibrium point $E_{1}$ through transcritical bifurcation at
$\delta=\delta_{T}\equiv\dfrac{1}{1+\alpha}.$
Now evaluating the Jacobian matrix for the system (4) at the interior
equilibrium point $E_{*}(u_{*},v_{*})$ we have

$J_{*}=\begin{pmatrix}\gamma(u_{*}(2-3u-2\beta)+\beta)-\dfrac{v_{*}}{(1+u_{*}\alpha)^{2}}&-\dfrac{u_{*}}{1+u_{*}\alpha}\\\
\dfrac{\varepsilon
v_{*}}{(1+u_{*}\alpha)^{2}}&\varepsilon\Big{(}\dfrac{u_{*}}{1+u_{*}\alpha}-\delta\Big{)}\end{pmatrix}.$

From the feasibility condition of $E_{*}$ we always have Det$(J_{*})>0$. The
interior equilibrium point is stable if Tr$(J_{*})<0$, and it loses its
stability via super-critical Hopf bifurcation when Tr$(J_{*})=0$ and is
unstable for Tr$(J_{*})>0$. The Hopf threshold $\delta=\delta_{H}$ can be
obtained by solving Tr$(J_{*})=0$ which on simplification gives

$\delta_{H}=\dfrac{1+\alpha^{2}\beta-\sqrt{1+\alpha+\alpha^{2}-\alpha\beta+\alpha^{2}\beta+\alpha^{2}\beta^{2}}}{\alpha(-1-\alpha+\alpha\beta+\alpha^{2}\beta)}.$
(5)

Transversality condition for Hopf bifurcation is satisfied at
$\delta=\delta_{H}.$ The coexistence steady state $E_{*}(u_{*},v_{*})$ is
stable for $\delta>\delta_{H}$ and it destabilizes for $\delta<\delta_{H}$,
surrounded by a stable limit cycle. The bifurcation diagrams of the system (3)
with $\delta$ as bifurcation parameter and for two different values of $\beta$
are plotted in Fig. 1. It is readily seen that, in case $\beta<1/\alpha$, with
an increase in the Allee threshold $\beta$ (hence making the Allee effect
weaker) the Hopf bifurcation point shifts to the left. Correspondingly, the
steady species coexistence occurs for a broader parameter range of $\delta$
(see Fig. 1). Furthermore, with the increase of $\beta$ the size of the limit
cycle is reduced. Increase in the strength of the weak Allee effect not only
enhance the stable coexistence rather reduces the amplitude of stable
oscillatory coexistence.

(a) $\beta=0.2$

(b) $\beta=0.8$

Figure 1: The bifurcation diagram of system (3) using $\delta$ as the
bifurcation parameter shown for two different values of $\beta$ and other
parameters as $\alpha=0.5,\ \gamma=3.$ Here the blue solid curve shows either
the stable steady state (for $\delta>\delta_{H}$) or the size of the stable
limit cycle (for $\delta<\delta_{H}$), the solid red line shows the
semitrivial ‘prey-only’ state in the range where it is a saddle (hence
unstable), and the dashed red line shows the unstable coexistence steady
state.

Interestingly, the linear stability results remains unaltered in the presence
of slow-fast time scale as the analytical conditions are independent of
$\varepsilon$. The linear stability analysis fails to capture the complete
dynamics of the slow-fast system (4) for $0<\varepsilon\ll 1$. The system (4)
exhibit catastrophic transition which cannot be captured by standard stability
analysis, rather the model may sometimes overestimate ecological resilience
[74]. Therefore, to study the complete dynamics of the system we take help of
geometric singular perturbation theory and blow-up technique which will be
discussed in next sections.

##  3 Slow-fast system

In this section, we shall describe the dynamics of the slow-fast system (4)
when $0<\varepsilon\ll 1$. To understand the dynamics of the system (4) for
sufficiently small $\varepsilon\ (>0)$ we need to consider the behaviour of
two subsystems corresponding to (4), which can be obtained for
$\varepsilon=0$. The system in its singular limit, $\varepsilon=0$ is obtained
as follows

$\displaystyle\dfrac{du}{dt}$ $\displaystyle=$ $\displaystyle f(u,v)=\gamma
u(1-u)(u+\beta)-\frac{uv}{1+\alpha u},$ (6a) $\displaystyle\dfrac{dv}{dt}$
$\displaystyle=$ $\displaystyle 0.$ (6b)

The above system is known as fast subsystem or layer system corresponding to
the slow-fast system (4). The fast flow consists with constant predator
density determined by the initial condition $v(0)=c$, and by integrating the
differential equation

$\displaystyle\dfrac{du}{dt}=\gamma u(1-u)(u+\beta)-\frac{uc}{1+\alpha u},$
(7)

with initial condition $u(0)>0$. The direction of the fast flow depends on the
choice of initial conditions $u(0)$, $v(0)$ and other parameter values. Green
horizontal lines are solution trajectories with appropriate direction as shown
in Fig. 2a. Now writing system (4) in terms of the slow time
$\tau:=\varepsilon t$, we get the equivalent system in terms of slow time
derivatives,

$\displaystyle\varepsilon\dfrac{du}{d\tau}$ $\displaystyle=$ $\displaystyle
f(u,v)=\gamma u(1-u)(u+\beta)-\frac{uv}{1+\alpha u},$ (8a)
$\displaystyle\dfrac{dv}{d\tau}$ $\displaystyle=$ $\displaystyle
g(u,v)=v\Big{(}\frac{u}{1+\alpha u}-\delta\Big{)}.$ (8b)

Substituting $\varepsilon=0$ in the above system to find the following
differential algebraic equation (DAE),

$\displaystyle 0$ $\displaystyle=$ $\displaystyle f(u,v)=\gamma
u(1-u)(u+\beta)-\frac{uv}{1+\alpha u},$ (9a) $\displaystyle\dfrac{dv}{d\tau}$
$\displaystyle=$ $\displaystyle g(u,v)=v\Big{(}\frac{u}{1+\alpha
u}-\delta\Big{)},$ (9b)

which is known as the slow subsystem corresponding to the slow-fast system
(8). The solution of the above system is constrained to the set
$\\{(u,v)\in\mathbb{R}^{2}_{+}:f(u,v)=0\\}$ and is known as critical manifold
$C_{0}$. This set has one-one correspondence with the set of equilibrium of
the system (7). This critical manifold consists of two of different manifolds

$\displaystyle C_{0}^{0}$ $\displaystyle=$
$\displaystyle\\{(u,v)\in\mathbb{R}^{2}_{+}:u=0,v\geq 0\\},$ $\displaystyle
C_{0}^{1}$ $\displaystyle=$
$\displaystyle\Big{\\{}(u,v)\in\mathbb{R}^{2}_{+}:v=q(u):=\gamma(1-u)(u+\beta)(1+\alpha
u),\ 0<u<1,\ v>0\Big{\\}},$

such that $C_{0}=C^{0}_{0}\cup C^{1}_{0}$ where $C_{0}^{0}$ is the positive
$v$-axis and $C_{0}^{1}$ is a portion of a parabola as shown in Fig. 2a,
marked with black colour. The slow flow on the critical manifold is given by

$\displaystyle\dfrac{du}{d\tau}$ $\displaystyle=$
$\displaystyle\dfrac{g(u,q(u))}{\dot{q}(u)},$ (10)

where $`.`$ refers to the differentiation w.r.t $u$. The solution of the
system (4) for sufficiently small $\varepsilon>0$ cannot be approximated from
its limiting solution at $\varepsilon=0$. Therefore, $\varepsilon=0$ is the
singular limit of the system (4). The solution of the full system is obtained
by combining the solution of the system in its singular limits. And depending
on the region in the phase space we use either of the subsystems.

For $\alpha=0.5$, $\beta=0.22$, $\gamma=3$, and $\delta=0.3$ the coexistence
steady-state is unstable for $0<\varepsilon\leq 1$ and is surrounded by a
stable limit cycle. Interestingly the size and shape of stable limit cycle
change with the variation of $\varepsilon$ which is shown in Fig. 2(b). The
size and shape of closed curve attractor (blue), obtained for
$\varepsilon=0.001$ is quite different from stable limit cycle (magenta) which
is obtained for $\varepsilon=1$. This shape does not change if we further
decrease $\varepsilon$, keeping other parameters fixed. This observation is
based upon the numerical simulation and we need detailed analysis to
understand the possible shape of the trajectories in singular limit
$\epsilon\rightarrow 0$. For $\varepsilon=0.001$, the closed curve attractor
(blue) consists of two horizontal segments on which flow is fast and one
curvilinear and vertical segment where the flow is slow. They are obtained by
concatenating the solution of the layer and the reduced subsystem
respectively. The two horizontal segments of the attractor (blue) are the
perturbed trajectories corresponding to the layer system. This signifies the
fast growth or decay of the prey species while predator density remains
unaltered. The vertical portion close to $v$-axis and curvilinear part are
close to the critical manifolds $C^{0}_{0}$ and $C^{1}_{0}$. Change in shape
and size of the attractor does not solely depend upon the magnitude of
$\varepsilon$ rather determined by the magnitude of the parameters involved
with the reaction kinetics and time scale parameter.

Now we fix $\varepsilon=0.01$ and other parameters as mentioned above, except
$\delta$. Small variation in $\delta$ just below the $\delta_{H}$ results in
rapid change in size and shape of the periodic attractor (see Fig. 2(c)). A
small limit cycle (cyan) appears for $\delta=0.3762$ known as canard cycle
without head. This cycle encounter a change in curvature when
$\delta=0.376165$ and the resulting cycle is known as canard cycle with head
(blue). Further decreasing $\delta$ to $0.36$ the system settles down to a
closed cycle known as relaxation oscillation. Further decrease in $\delta$
does not alter the size and shape of the closed attractor and the trajectories
converge to the stable relaxation oscillation cycle even for $\varepsilon$
sufficiently small. In the next section, we will derive the analytical
conditions for the existence of canard cycle and relaxation oscillation. The
analytical results will help us to identify the domains in the prarametric
plane where we can find these different types of closed curve attractors.

Figure 2: (a) Dynamics of the slow-fast system (4) where single arrow represnt
slow flow and double arrow represent fast flow, (b) limit cycles for different
values of time-scale parameter $\varepsilon=1$ (magenta), $\varepsilon=0.1$
(green), $\varepsilon=0.001$ (blue), (c) canard cycles for different values of
the bifurcation parameter $\delta$ for $\varepsilon=0.01$, canard cycle
without head for $\delta=0.3762$ (cyan), canard cycle with head for
$\delta=0.376165$ (blue), relaxation oscillation for $\delta=0.36$ (magenta)
and other parameters are mentioned in the text.

##  4 Analysis of slow-fast system

The critical manifold $C_{0}^{1}$ can be divided into two parts, one part
consists of the attractors of the fast sub-system and another part is
repelling in nature. The attracting and repelling part of the manifold is
separated by a non-degenerate fold point $P$. The fold point $P(u_{f},v_{f})$
is characterized by the following conditions [27].

$\frac{\partial{f}}{\partial u}(u_{f},v_{f})=0,\ \frac{\partial f}{\partial
v}(u_{f},v_{f})\neq 0,\ \frac{\partial^{2}f}{\partial u^{2}}(u_{f},v_{f})\neq
0,\ \text{and}\ g(u_{f},v_{f})\neq 0.$

The components of the fold point are given by

$u_{f}=\dfrac{(\alpha-\alpha\beta-1)+(1+\alpha+\alpha^{2}-\alpha\beta+\alpha^{2}\beta+\alpha^{2}\beta^{2})^{\frac{1}{2}}}{3\alpha},$
$v_{f}=\gamma(1-u_{f})(u_{f}+\beta)(1+\alpha u_{f}),$

which is the maximum point on the critical manifold. The fold point divides
the critical manifold into normally hyperbolic attracting ($C^{1,a}_{0}$) and
repelling ($C^{1,r}_{0}$) submanifolds given by

$\displaystyle C_{0}^{1,a}$
$\displaystyle=\Big{\\{}(u,v)\in\mathbb{R}^{2}_{+}:v=q(u),u_{f}<u\leq
1\Big{\\}}$ $\displaystyle C_{0}^{1,r}$
$\displaystyle=\Big{\\{}(u,v)\in\mathbb{R}^{2}_{+}:v=q(u),0\leq
u<u_{f}\Big{\\}}.$

The point of intersection of $C^{1}_{0}$ with the vertical $v$-axis is
$T_{C}(0,\beta\gamma)$, which is the transcritical bifurcation point for the
fast subsystem. It follows from Fenichel’s theorem [11, 29], there exist
locally invariant slow sub-manifolds $C^{1}_{\varepsilon}$ and
$C_{\varepsilon}^{0}$ which are diffeomorphic to the respective critical
manifolds $C^{1}_{0}$ and $C^{0}_{0}$, except at the non-hyperbolic points $P$
and $T$. $C^{1}_{0}$ can be written explicitly as $v=q(u)$, we assume the
invariant manifold $C_{\varepsilon}^{1}$ can be written as a perturbation of
$v=q(u)$ as follows, with $\varepsilon$ as perturbation parameter,

$C_{\varepsilon}^{1}=\Big{\\{}(u,v)\in\mathbb{R}^{2}_{+}:v=q(u,\varepsilon),\
0<u<1,\ v>0\Big{\\}},$ (11)

where $q(u,\varepsilon)=q_{0}(u)+\varepsilon
q_{1}(u)+\varepsilon^{2}q_{2}(u)+\cdots$, and
$C^{0}_{\varepsilon}=\\{(u,v)\in\mathbb{R}^{2}_{+}:u=0,v\geq 0\\}$. Using the
invariance condition and the asymptotic expansion of $q(u,\varepsilon)$, we
can find the perturbed invariant manifold approximated up to the desired
order. The approximation of $q(u,\varepsilon)$ up to second order is provided
in Appendix A with explicit expressions for $q_{0}$, $q_{1}$, $q_{2}$. The
approximations of invariant manifold for different values of $\varepsilon$ are
shown in Fig. 3. This approximation has two non-removable discontinuities in
the vicinity of the non-hyperbolic points $P$ and $T_{C}$.

The critical manifold $C^{1}_{0}$ is normally hyperbolic except at the points
$P(u_{f},v_{f})$ and $T_{C}(0,\beta\gamma)$ and so is $C^{1}_{\varepsilon}$.
Thus any trajectory starting near the attracting (repelling) submanifold
$C^{1,a}_{0}$ $(C^{1,r}_{0})$ cannot cross the fold point $P$ (transcritical
point $T_{C}$). We can see from Fig. 2 that for sufficiently small values of
$\varepsilon$ the trajectories pass enough close to the attracting manifold
$C^{1,a}_{0}$ and cross the point $P$. Fenichel’s theory is not adequate to
determine the analytical expression for perturbed sub-manifolds close to
$C_{0}^{1}$ and is continuous in the vicinity of the non-hyperbolic points.

Figure 3:  Second order approximation of perturbed manifold from GSPT with
$\varepsilon=1$ (magenta), $\varepsilon=0.1$(green), $\varepsilon=0.01$(red),
$\varepsilon=0.001$(blue); for the parameter values $\alpha=0.5,\ \beta=0.22,\
\gamma=3,\ \delta=0.3$ .

Therefore, to construct a trajectory passing through the vicinity of the point
$P$ we must remove the singularity at this point. Depending on parameter
$\delta,$ the predator nullcline intersects either $C^{1,a}_{0}$ or
$C_{0}^{1,r}$ or passes through the point $P$. Thus the coexistence
equilibrium point $E_{*}$ lies either on $C^{1,a}_{0}$ or $C_{0}^{1,r}$ or
coincides with $P$. When $E_{*}$ lie on $C_{0}^{1,a}$ it is globally
asymptotically stable and every trajectory converges to $E_{*}$. When $E_{*}$
coincides with the fold point $P$ then $f_{u}(E_{*})=0$, $f_{v}(E_{*})\neq 0$,
$f_{uu}(E_{*})\neq 0$, and $g(E_{*})=0$, this point is called the canard
point. For the model under consideration (4), the Hopf point coincides with
the canard point. The solution passing through this point is known as the
canard solution. For $\delta<\delta_{H}$, the coexistence equilibrium point
$E_{*}$ lies on the repelling sub-manifold $C_{0}^{1,r}$, it is unstable and
we obtain a special kind of periodic solution consisting of two fast flow
(almost horizontal) and two slow flow (passing close to $C^{1,a}_{0}$ and
$C^{0}_{0}$) called relaxation oscillation. We will discuss the existence of
such solutions in consequent subsections.

To remove the singularity at the fold point, we use blow-up transformation at
the non-hyperbolic fold point which will extend the system over a 3-sphere in
$\mathbb{R}^{4}$, denoted by $S^{3}=\\{x\in\mathbb{R}^{4}:||x||=1\\}$. Using
the blow-up technique we remove the singularity from the system and determine
the canard solution passing through this point.

To apply the blow-up technique, first, we transform the slow-fast system (4)
into its desired slow-fast normal form.

###  4.1 Slow-fast normal form

Here we consider a topologically equivalent form of the system (4) by re-
scaling the time with the help of a transformation $t\rightarrow(1+\alpha
u)t$, where $(1+\alpha u)>0$ [30]. The transformed system is

$\displaystyle\dfrac{du}{dt}$ $\displaystyle=\gamma u(1-u)(u+\beta)(1+\alpha
u)-uv\equiv F(u,v,\delta),$ (12) $\displaystyle\dfrac{dv}{dt}$
$\displaystyle=\varepsilon\left(uv-\delta v(1+\alpha
u)\right)\equiv\varepsilon G(u,v,\delta).$

The fold point $P$ coincides with the coexistence equilibrium point at
$\delta=\delta_{*}$. As a consequence, the following conditions hold

$\displaystyle F(u_{*},v_{*},\delta_{*})=0,\ F_{u}(u_{*},v_{*},\delta_{*})=0,\
F_{v}(u_{*},v_{*},\delta_{*})\neq 0,\ F_{uu}(u_{*},v_{*},\delta_{*})\neq 0,$
(13) $\displaystyle\ G_{u}(u_{*},v_{*},\delta_{*})\neq 0,\
G_{\delta}(u_{*},v_{*},\delta_{*})\neq 0\ \text{and}\
G(u_{*},v_{*},\delta_{*})=0.$

Using the transformation $U=u-u_{*},\ V=v-v_{*},\ \lambda=\delta-\delta_{*}$,
we translate the fold point to the origin, and together with the conditions
(13) the system reduces to the slow-fast normal form near $(0,0)$ as follows

$\displaystyle\dfrac{dU}{dt}$
$\displaystyle=-Vh_{1}(U,V)+U^{2}h_{2}(U,V)+\varepsilon h_{3}(U,V),$ (14)
$\displaystyle\dfrac{dV}{dt}$
$\displaystyle=\varepsilon\left(Uh_{4}(U,V)-\lambda
h_{5}(U,V)+Vh_{6}(U,V)\right),$

where $h_{i}^{\prime}s,\ i=1,2,3\cdots 6$ are given in Appendix B.

Here $\lambda$ measures the perturbation of $\delta$ from $\delta_{*}$ and is
considered as bifurcation parameter for the system (14). The bifurcation
parameter $\lambda$ and time-scale parameter $\varepsilon$ are assumed to be
independent of time. We now extend the above system to $\mathbb{R}^{4}$ by
augmenting the equations $\dfrac{d\lambda}{dt}=0$ and
$\dfrac{d\varepsilon}{dt}=0$ to system (14) and study the dynamics of the
system in the vicinity of $(0,0,0,0)$.

###  4.2 Blow-up transformation

The fold point $P$ of the system (4) and the equilibrium point $E_{*}$
coincides at the Hopf bifurcation threshold. Hence $P$ is now a Canard point.
Let us consider the blow-up space
$S^{3}=\\{(\bar{U},\bar{V},\bar{\lambda},\bar{\varepsilon})\in\mathbb{R}^{4}:\bar{U}^{2}+\bar{V}^{2}+\bar{\lambda}^{2}+\bar{\varepsilon}^{2}=1\\}$
and an interval $\mathcal{I}:=[0,\rho]$ where $\rho>0$ is a small constant. We
define a manifold $\mathcal{M}:=S^{3}\times\mathcal{I}$ and the blow-up map
$\Phi$, $\Phi:\mathcal{M}\rightarrow\mathbb{R}^{4}$ where

$\Phi(\bar{U},\bar{V},\bar{\lambda},\bar{\varepsilon},\bar{r})=(\bar{r}\bar{U},\bar{r}^{2}\bar{V},\bar{r}\bar{\lambda},\bar{r}^{2}\bar{\varepsilon}):=(U,V,\lambda,\varepsilon).$
(15)

Using the above map we can write the transformed system as follows

$\displaystyle\dfrac{d\bar{U}}{dt}$
$\displaystyle=\dfrac{1}{\bar{r}}\left(\dfrac{dU}{dt}-\bar{U}\dfrac{d\bar{r}}{dt}\right),\
\dfrac{d\bar{V}}{dt}=\dfrac{1}{\bar{r}^{2}}\left(\dfrac{dV}{dt}-2\bar{r}\bar{V}\dfrac{d\bar{r}}{dt}\right),$
(16) $\displaystyle\dfrac{d\bar{\lambda}}{dt}$
$\displaystyle=\dfrac{1}{\bar{r}}\left(\dfrac{d\lambda}{dt}-\bar{\lambda}\dfrac{d\bar{r}}{dt}\right),\
\dfrac{d\bar{\varepsilon}}{dt}=\dfrac{1}{\bar{r}^{2}}\left(\dfrac{d\varepsilon}{dt}-2\bar{r}\bar{\varepsilon}\dfrac{d\bar{r}}{dt}\right),$

where
$\dfrac{d\bar{U}}{dt},\dfrac{d\bar{V}}{dt},\dfrac{d\bar{\lambda}}{dt},\dfrac{d\bar{\varepsilon}}{dt}$
are given in (14). To study the dynamics of the transformed system on and
around the hemisphere $S^{3}_{\varepsilon\geq 0}$ we will introduce the charts
with direction blow-up maps [29, 76]. Along each direction of the coordinate
axis we define the charts $K_{1}$, $K_{2}$, $K_{3}$ and $K_{4}$ by setting
$\bar{V}=1$, $\bar{\varepsilon}=1$, $\bar{U}=1$ and $\bar{\lambda}=1$
respectively, in (16). The charts $K_{1}$ and $K_{3}$ describe the dynamics in
the neighborhood of the equator of $S^{3}$ and $K_{2}$ describes the dynamics
in a neighborhood of the positive hemisphere. Here, we mainly focus on chart
$K_{2}$ to prove the existence of a periodic solution for
$0<\bar{\varepsilon}\ll 1$. Re-scaling the time with the transformation
$\bar{t}:=\bar{r}t$, we desingularize the system (16) so that the
multiplicative factor $\bar{r}$ disappears. The transformed version of the
system (16) can be written in chart $K_{2}$ as follows

$\displaystyle\dfrac{d\bar{U}}{dt}$
$\displaystyle=-\bar{V}b_{1}+\bar{U}^{2}b_{2}+\bar{r}\left(a_{1}\bar{U}-a_{2}\bar{U}\bar{V}+a_{3}\bar{U}^{3}\right)+O(\bar{r}(\bar{\lambda}+\bar{r})),$
(17) $\displaystyle\dfrac{d\bar{V}}{dt}$
$\displaystyle=\bar{U}b_{3}-\bar{\lambda}b_{4}+\bar{r}\left(a_{4}\bar{U}^{2}+a_{5}\bar{V}\right)+O(\bar{r}(\bar{\lambda}+\bar{r}),$
$\displaystyle\dfrac{d\bar{\lambda}}{dt}$ $\displaystyle=0,$
$\displaystyle\dfrac{d\bar{\varepsilon}}{dt}$ $\displaystyle=0,$

where $b_{j}$’s are given in Appendix B,

$\ a_{1}=a_{5}=0,\ \ a_{2}=1,\ \
a_{3}=-\left(\gamma+\alpha\gamma(4u_{*}+\beta-1)\right),\ \
a_{5}=u_{*}-(1+u_{*}\alpha)\delta_{*}.$ (18)

The condition for the destabilization of the coexistence equilibrium point
through singular Hopf-bifurcation is summarized in the following theorem. The
Hopf bifurcation of the system (14) occurs at $\lambda=0$. At Hopf
bifurcation, the purely imaginary eigenvalues of the corresponding Jacobian
matrix are a function of $\varepsilon$, which tends to zero as
$\varepsilon\rightarrow 0$. Also, in slow-time derivative, the eigenvalues are
a function of $1/\varepsilon.$ Thus on both the time scales the Hopf
bifurcation is singular as $\varepsilon\rightarrow 0.$

######  Theorem 4.1.

Let $(U,V)=(0,0)$ is the canard point of the transformed system (14) at
$\lambda=0$ such that $(0,0)$ is a folded singularity and $G(0,0,0)=0$. Then
for sufficiently small $\varepsilon$ there exist a singular Hopf bifurcation
curve $\lambda=\mathcal{\lambda_{H}}(\sqrt{\varepsilon})$ such that the
equilibrium point $p$ of the system (14) is stable for
$\lambda>\mathcal{\lambda_{H}}(\sqrt{\varepsilon})$ and

$\mathcal{\lambda_{H}}(\sqrt{\varepsilon})=-\dfrac{b_{3}(a_{1}+a_{5})}{2b_{2}b_{4}}\varepsilon+O(\varepsilon^{3/2}).$
(19)

###### Proof.

The proof of the theorem is given in Appendix B.∎

The singular Hopf bifurcation curve for the system (12) is thus given by

$\displaystyle\mathcal{\delta_{H}}(\sqrt{\varepsilon})=$
$\displaystyle\dfrac{1+\alpha^{2}\beta-\sqrt{1+\alpha+\alpha^{2}-\alpha\beta+\alpha^{2}\beta+\alpha^{2}\beta^{2}}}{\alpha(-1-\alpha+\alpha\beta+\alpha^{2}\beta)}-$
(20)
$\displaystyle\quad\dfrac{b_{3}(a_{1}+a_{5})}{2b_{2}b_{4}}\varepsilon+O(\varepsilon^{3/2}).$

In Fig. 4 the singular Hopf bifurcation curve (red) is plotted in
$\delta-\varepsilon$ parametric plane, it clearly explains how the singular
Hopf-bifurcation threshold changes with the variation in $\varepsilon$. Once
the coexistence equilibrium loses stability through Hopf bifurcation, at the
Canard point, we find a closed orbit as attractor surrounding the unstable
equilibrium point. From this point small amplitude stable canard cycle
originates enclosing the point $P$ and then forms canard cycle with head
depending on the parameter values. The following theorem provides an
anlaytical expression of the maximal canard curve in $(\lambda-\varepsilon)$
plane.

######  Theorem 4.2.

Let $(U,V)=(0,0)$ is the canard point of the slow-fast normal form (14) at
$\lambda=0$ such that $(0,0)$ is a folded singularity and $G(0,0,0)=0$. Then
for $\varepsilon>0$ sufficiently small there exists maximal canard curve
$\lambda=\lambda_{c}(\sqrt{\varepsilon})$ such that the slow flow on the
normally hyperbolic invariant submanifolds $\mathcal{M}_{\varepsilon}^{1,a}$
connects with $\mathcal{M}_{\varepsilon}^{1,r}$ in the blow-up space. And
$\lambda_{c}(\sqrt{\varepsilon})$ is given by

$\displaystyle\lambda_{c}(\sqrt{\varepsilon})$
$\displaystyle=-\dfrac{1}{A_{5}}\Big{(}\dfrac{3A_{1}}{4A_{4}^{2}}+\dfrac{A_{2}}{2A_{4}}+A_{3}\Big{)}\varepsilon+O(\varepsilon^{3/2})$
(21)

###### Proof.

The proof of this theorem is given in Appendix C. ∎

The maximal canard curve, along which the canard cycle with head appears for
the system (12) is given by

$\displaystyle\delta_{c}(\sqrt{\varepsilon})=$
$\displaystyle\dfrac{1+\alpha^{2}\beta-\sqrt{1+\alpha+\alpha^{2}-\alpha\beta+\alpha^{2}\beta+\alpha^{2}\beta^{2}}}{\alpha(-1-\alpha+\alpha\beta+\alpha^{2}\beta)}-$
(22)
$\displaystyle\quad\dfrac{1}{A_{5}}\Big{(}\dfrac{3A_{1}}{4A_{4}^{2}}+\dfrac{A_{2}}{2A_{4}}+A_{3}\Big{)}\varepsilon+O(\varepsilon^{3/2}).$

Keeping $\alpha$, $\beta$ and $\varepsilon\ (>0)$ fixed,
$\delta_{c}(\sqrt{\varepsilon})$ gives the threshold for the existence of
canard cycle with head. A schematic diagram of the threshold curves in
$\delta-\varepsilon$-plane is illustrated in Fig. 4 and it divides the
$\delta-\varepsilon$ parametric plane into four domains.

Figure 4: Schematic diagram showing singular Hopf bifurcation curve
$\mathcal{\delta_{H}}$ (red), maximal canard curve $\delta_{c}$ (blue), and
$\delta_{ro}$ (dashed black) relaxation oscillation cycle.

In domain I, when $\delta>\delta_{\mathcal{H}}$, the coexistence equilibrium
point is stable. For a fixed $\varepsilon>0$, as we decrease $\delta$ from
domain I to domain II small amplitude canard cycles appear after crossing the
Hopf bifurcation threshold $\delta=\delta_{\mathcal{H}}.$ In domain II, that
is when $\delta_{c}<\delta<\delta_{\mathcal{H}},$ the system experiences a
transition from canard cycle with head to canard cycle without head. The size
of the canard cycle increases on decreasing $\delta$ and the shape of the
cycle changes to canard with a head at $\delta=\delta_{c}.$ The canard cycle
whith head persists in a narrow domain III, where
$\delta_{ro}<\delta<\delta_{c}$. On further decreasing $\delta,$ that is, when
$\delta\leq\delta_{ro}$ the unstable equilibrium point is surrounded by a
stable periodic attractor called relaxation oscillation. This periodic
attractor consists of two concatenated slow (close to the critical manifold)
and fast (almost horizontal and away from the critical manifold) flow. We can
see that for sufficiently small $\varepsilon$, this transition, from small
canard cycle to relaxation oscillation through canard cycle with head, takes
place within a narrow interval of the parameter $\delta$ and the phenomenon is
known as canard explosion. This mechanism is further illustrated with the help
of numerical example in the sub-section 4.4.

###  4.3 Entry-Exit Function

The canard cycle and relaxation oscillation pass through the fold point $P$
when $\varepsilon=0$. These closed attractors pass through the vicinity of $P$
for $\varepsilon\ll 1$. We now prove the existence of a trajectory that jumps
from the fold point to the other attracting slow manifold $C^{0}_{0}$ through
fast horizontal flow and continuing there for a constant time the trajectory
leaves $C^{0}_{0}$ at a certain point. This is determined by the entry-exit
function and we can find the coordinates of the exit point from the slow
manifold $C^{0}_{0}$. To do this, first, rewrite the system (4) as a
Kolmogorov system [12] as follows

$\displaystyle\dfrac{du}{dt}$ $\displaystyle=$ $\displaystyle uf_{1}(u,v)$
$\displaystyle=u\Big{(}\gamma(1-u)(u+\beta)-\frac{v}{1+\alpha u}\Big{)},$ (23)
$\displaystyle\dfrac{dv}{dt}$ $\displaystyle=$ $\displaystyle\varepsilon
vg_{1}(u,v)$ $\displaystyle=\varepsilon v\Big{(}\frac{u}{1+\alpha
u}-\delta\Big{)}.$

We can verify $f_{1}(0,v)=\gamma\beta-v$, $g_{1}(0,v)=-\delta<0$, which
implies that $f_{1}(0,v)<0$ if $v>\gamma\beta$, $f_{1}(0,v)>0$ if
$v<\gamma\beta$. $T\ (0,\beta\gamma)$ on the vertical axis is the
transcritical bifurcation point and we can divide the slow manifold
$C^{0}_{0}$ into two parts $V^{+}:=\\{(u,v):u=0,v>\gamma\beta\\}$ and
$V^{-}:=\\{(u,v):u=0,\ v<\gamma\beta\\}$. Clearly $V^{+}$ is attracting and
$V^{-}$ is repelling.

Let us fix $\varepsilon>0$ and let $u_{max}$ be the point of maximum of the
critical manifold $C^{1}_{0}$ obtained from the extremum condition
$\dot{q_{0}}(u)=0,$ where $q_{0}(u)$ is given in Appendix A. Solving for
$u_{max}$ we find

$u_{max}=\dfrac{(\alpha-\alpha\beta-1)+\sqrt{1+\alpha+\alpha^{2}-\alpha\beta+\alpha^{2}\beta+\alpha^{2}\beta^{2}}}{3\alpha}$
(24)

and from the expression of $C_{0}^{1}$ we have

$\displaystyle v_{max}=q(u_{max})\ =\
\gamma(1-u_{max})(u_{max}+\beta)(1+\alpha u_{max}).$ (25)

Now we consider a trajectory starting from a point, say $(u_{1},v_{1})$, where
$u_{1}<u_{max}$ and $v_{1}=v_{max}$. The trajectory gets attracted toward the
attracting manifold $V_{+}$ and starts moving downward maintaining proximity
to $V^{+}.$ It was expected that the trajectory would leave the vertical axis
at the bifurcation point $T$ where it loses its stability [42]. The trajectory
crosses the point $T$ and continues to move vertically downward remaining
close to the repelling part $V^{-}$, for a certain time, until a minimum
predator population $p(v_{1})$ is attained s.t. $0<p(v_{1})<\gamma\beta$.
After leaving the slow manifold near the point $p(v_{1})$, the trajectory
starts moving along a fast horizontal segment and gets attracted to attracting
slow manifold $C_{0}^{1,a}$. This point of exit is determined by an implicit
function $p(v_{1})$, called entry-exit function, which is defined implicitly
as

$\int_{p(v_{1})}^{v_{1}}\dfrac{f_{1}(0,v)}{vg_{1}(0,v)}dv=0.$

For simplicity we define $v_{0}:=p(v_{1})$, then we have

$\displaystyle\int_{v_{0}}^{v_{1}}\dfrac{v-\gamma\beta}{v\delta}dv=0\,\,\implies(v_{1}-v_{0})-\gamma\beta\ln\Big{(}\frac{v_{1}}{v_{0}}\Big{)}=0$
(26)

Substituting (25) into equation (26) we obtain a transcendental equation in
$v_{0}$ which we solve numerically to obtain the exit point.
For the parameter values $\alpha=0.5$, $\beta=0.2$, $\gamma=3$, $\delta=0.3$
we obtain $u_{max}=0.472$, $v_{1}=1.316$, and solving the transcendental
equation (26) we get $p(v_{1})=0.207509$, which is the exit point from the
manifold $C_{0}^{0}$.

######  Theorem 4.3.

Let $P$ be the fold point on the critical manifold $C^{1}_{0}$ where the slow
flow on the attracting manifold $C_{0}^{1,a}$ is given by (10). Also assume
that the coexistence equilibrium point lies on the normally hyperbolic
repelling critical submanifold under the parametric restriction,

$\dfrac{\delta}{1-\alpha\delta}<\dfrac{(\alpha-\alpha\beta-1)+\sqrt{1+\alpha+\alpha^{2}-\alpha\beta+\alpha^{2}\beta+\alpha^{2}\beta^{2}}}{3\alpha}$

and let $U$ denotes a small neighborhood of a singular trajectory $\gamma_{0}$
consisting of alternate slow and fast trajectories. Then for sufficiently
small $\varepsilon$ there exist a unique attracting limit cycle
$\gamma_{\varepsilon}\subset U$ such that
$\gamma_{\varepsilon}\rightarrow\gamma_{0}$ as $\varepsilon\rightarrow 0$.

###### Proof.

The proof is given in Appendix D. ∎

These cycles are shown with the help of a numerical example in Appendix D.
When $\varepsilon\rightarrow 0$, all the trajectories asymptotically converge
to this stable limit cycle consisting of alternate slow and fast transitions
of prey and predator densities. This cycle can be interpreted as: when the
predator population reaches some high density there is a rapid decline in the
prey population due to excessive consumption by the specialist predator and
the prey reaches a considerably low level. As a consequence, the predator
population declines slowly until it reaches a low threshold density at which
the prey population again starts growing. Consequently, the prey regenerates
within a very short time while predator density remains more or less fixed.
Finally, the predator population starts growing slowly due to the abundance of
resources. Finally, when the predator density reaches its maximum level, the
slow-fast cycle completes and this dynamics continues with time.

###  4.4 Canard Explosion

In the previous sub-sections, we have observed the periodic dynamics of the
slow-fast system near the canard point, where the predator nullcline
intersects the non-trivial prey nullcline at the fold point. This occurs at a
certain threshold of the parameter $\delta$. At this point, the coexistence
equilibrium point loses stability through singular Hopf bifurcation and a
small amplitude stable limit cycle is observed. Due to the decrease of the
parameter $\delta$, the Hopf bifurcating stable cycle grows in size and
settles down to relaxation oscillation. The fast transition in the size of the
limit cycle from small canard cycles to relaxation oscillation occurs in an
exponentially small range of the parameter $\delta$. This phenomenon is known
as the canard explosion.

The family of canard cycles are already shown in Fig. 2(c) for fixed
$\varepsilon$ and three values of $\delta$ close to the singular Hopf
bifurcation threshold $\delta_{H}$. The coexistence equilibrium is stable for
$\delta>\delta_{H}$ and the trajectory converges to the stable steady state
for any initial condition as it is the global attractor. We can see that for
$\delta$ just below $\delta_{H}$, a stable limit cycle grows in size and a new
periodic solution emerges known as the canard cycle without a head Fig. 2(c)
(cyan color). This marks the onset of the canard explosion. Further decreasing
$\delta$ slightly we obtain another canard cycle known as canard with head
Fig. 2(c) (blue color). This cycle is special in the sense that from the
vicinity of the fold point it follows the repelling slow manifold
$C^{1,r}_{0}$ for $O(1)$ time, before jumping to another attracting manifold.
A maximal canard is obtained at $\delta=\delta_{c}$. After crossing the
maximal canard threshold, the system settles down to a large stable periodic
solution called relaxation oscillation, which marks the end of canard
explosion. This orbit is characterized by the fact that the slow flow on
reaching the vicinity of the fold point directly jumps to another attracting
slow manifold, as studied in the previous section. The strength of the Allee
effect has a significant influence on the amplitude of stable oscillatory
coexistence of both the species. The change in the amplitude of the limit
cycle corresponding to canard explosion is shown in Fig. 5.

Figure 5: The bifurcation diagram showing the change in the amplitude of the
canard cycles is plotted against $\delta$ for $\alpha=0.5,\ \gamma=3,\
\varepsilon=0.01$ $\beta=0.22$

For smaller values of $\beta$, the size of the canard cycle is very large and
the canard explosion occurs in an exponentially small interval. However, on
increasing the value of $\beta$ the size of the limit cycle shrinks and
instead of a sudden change in the size of the cycle, we observe a gradual
increase in the amplitude of the periodic solution Fig. 6. Though the
transition from canard cycle to relaxation oscillation takes place in a much
wider parametric interval, in this case, it is difficult to distinctively
identify the different periodic solutions. For smaller values of $\beta$, when
the prey population is almost absent, the predator population also slowly
declines to an almost endemic level. But on increasing the strength of the
Allee effect the predator density never collapses but stays at the system at a
much higher density. After which the system experience a sudden outbreak in
the prey population within a very short interval of time. Again because of the
abundance of resources, the predator population grows slowly and reaches the
maximum capacity. Due to the exploitation of the resources, there is a fast
decline in the prey density and this cycle continues.

Figure 6: (a) The canard cycles for $\alpha=0.5,\ \beta=0.8,\ \gamma=3,\
\varepsilon=0.01$ and for different values of $\delta$ i.e $\delta=0.234$
(green), $\delta=0.233$ (blue), $\delta=0.231$ (magenta), (b) the bifurcation
diagram showing the change in the size of the cycles.

##  5 Spatio-temporal model

We now consider the spatio-temporal model corresponding to the model (4) with
slow-fast time scale. Here, we consider that the prey and predator densities
are functions of time and space, $u(t,{\bf x})$ and $v(t,{\bf x})$ denote prey
and predator densities, respectively, at time $t$ and at spatial location
${\bf x}$. In case of one dimensional (1D) space ${\bf x}=x\in\mathbb{R}$ and
for two dimensional (2D) space ${\bf x}=(x,y)\in\mathbb{R}^{2}$. For
simplicity we assume that ${\bf x}$ belongs to a bounded domain
$D\subset\mathbb{R}$ and $D\subset\mathbb{R}^{2}$ respectively. The spatio-
temporal dynamics of the prey-predator interaction is described by the
following reaction-diffusion equation

$\displaystyle u_{t}$ $\displaystyle=$ $\displaystyle\gamma
u(1-u)(u+\beta)-\frac{uv}{1+\alpha u}+\nabla^{2}u,$ (27a) $\displaystyle
v_{t}$ $\displaystyle=$ $\displaystyle\varepsilon\left(\frac{uv}{1+\alpha
u}-\delta v\right)+d\nabla^{2}v,$ (27b)

where $d$ is the ratio of diffusivity coefficients of predator to prey and
$\nabla^{2}$ is the Laplacian operator. The above spatio-temporal model is
subject to no-flux boundary condition and non-negative initial condition. The
model (27) can not produce any stationary Turing pattern and it can be proved
that the Turing instability condition is not satisfied. However, instability
of the coexistence steady-state due to Hopf bifurcation combined with the
diffusivity of two species leads to some dynamic pattern due to the formation
of traveling wave, wave of invasion and spatio-temporal chaos. The mechanisms
responsible for such kind of pattern formation are described in [33].

Analytical condition for the existence of traveling wave leads to successful
invasion by specialist predator is derived in the next subsection. For
simplicity of mathematical calculation, we restrict ourselves to one
dimensional space to explain the existence of traveling wave. In case of two
dimensional spatial domain, the analogous patterns are presented separately.

###  5.1 Existence of Traveling wave

To study the successful invasion by the predator we consider the system (27)
in one dimensional space, we re-write the above system as

$\displaystyle\frac{\partial u(t,x)}{\partial t}$ $\displaystyle=$
$\displaystyle f(u,v)+\frac{\partial^{2}u}{\partial x^{2}},$ (28a)
$\displaystyle\frac{\partial v(t,x)}{\partial t}$ $\displaystyle=$
$\displaystyle\varepsilon g(u,v)+d\frac{\partial^{2}v}{\partial x^{2}},$ (28b)

where $f(u,v)=\gamma u(1-u)(u+\beta)-\frac{uv}{1+\alpha u},\
g(u,v)=\frac{uv}{1+\alpha u}-\delta v.$ The predator is introduced in a small
domain where the prey density is at its carrying capacity. The successful
invasion of the predator is characterized by the existence of a traveling wave
joining the predator free steady-state with the coexistence steady-state.
Depending upon the stability of the coexistence steady-state, we can find
monotone traveling wave, non-monotonic traveling wave, and periodic traveling
wave as explained below with the help of numerical examples.

We begin with deriving the minimum speed of the traveling wave which result in
the successful invasion of the specialist predator into the space already
inhabited by its prey. For this, we first consider a single-species model with
the linear growth:

$\frac{\partial v(t,x)}{\partial t}=\alpha v+D\frac{\partial^{2}v}{\partial
x^{2}},$

where $\alpha$, $D\,>\,0$ are parameters with obvious meaning. Stricktly
speaking, the above equation does not possess a traveling wave solution.
However, for a compact initial condition, it is known that the tail of the
profile propagates with the constant speed given by $c_{min}=2\sqrt{D\alpha}$,
see [33], sometimes referred to as the Fisher spreading speed.

For the invasion of predator into space inhabited by its prey, we consider the
tail of the profile where $u\approx 1$ and $v\approx 0$ and linearize (28b)
around $(1,0)$:

$\frac{\partial v(t,x)}{\partial
t}=\varepsilon\Big{(}\frac{1}{\alpha+1}-\delta\Big{)}v+d\frac{\partial^{2}v}{\partial
x^{2}}.$ (29)

Clearly, the speed of the traveling wave, at the onset of successful invasion,
is given by

$c_{v}=2\Big{(}\varepsilon
d\Big{[}\frac{1}{\alpha+1}-\delta\Big{]}\Big{)}^{1/2}.$ (30)

The feasibility condition is $\delta(\alpha+1)<1$. The expression for $c_{v}$
indicates that the speed of traveling wave reduces in the order of
$\sqrt{\varepsilon}$. We consider the existence of traveling wave starting
from the predator free steady-state, which leads to the successful
establishment of the predators, this requires the consideration of the system
(28) with the following conditions

$u(t,x)=1,\ \text{and \ }v(t,x)=0,\ \text{as\ }x\rightarrow-\infty,\
\forall\,t,$ $u(t,x)=u_{*},\ \text{and \ }v(t,x)=v_{*},\ \text{as\
}x\rightarrow\infty,\ \forall\,t.$

We consider the traveling wave solution of the system (28) in the form
$u(t,x)=\phi(\xi)$, $v(t,x)=\psi(\xi)$ where $\xi=x-ct$ and $c$ is the wave
speed. The functions $\phi(\xi)$ and $\psi(\xi)$ thus satisfy the equations

$\displaystyle\dfrac{d^{2}\phi}{d\xi^{2}}+c\dfrac{d\phi}{d\xi}+f(\phi,\psi)$
$\displaystyle=0,$ (31) $\displaystyle
d\dfrac{d^{2}\psi}{d\xi^{2}}+c\dfrac{d\psi}{d\xi}+\varepsilon g(\phi,\psi)$
$\displaystyle=0.$

Substituting $p(\xi)=-\dfrac{d\phi}{d\xi}$ and $q(\xi)=-\dfrac{d\psi}{d\xi}$,
from (31) we can derive four coupled ordinary differential equations as
follows

$\displaystyle\dfrac{d\phi}{d\xi}$ $\displaystyle=-p,$ (32)
$\displaystyle\dfrac{dp}{d\xi}$ $\displaystyle=-cp+f(\phi,\psi),$
$\displaystyle\dfrac{d\psi}{d\xi}$ $\displaystyle=-q,$
$\displaystyle\dfrac{dq}{d\xi}$ $\displaystyle=\dfrac{1}{d}(-cq+\varepsilon
g(\phi,\psi)).$

Three homogeneous steady states ($E_{0},\ E_{1},\ E_{*}$) of the spatio-
temporal model (28) corresponds to three steady states of system (32) are
$Q_{0}(0,0,0,0)$, $Q_{1}(1,0,0,0)$ and $Q_{*}(u_{*},0,v_{*},0)$. To ensure the
successful invasion of the predator, we focus on the dynamics of the system
(32) around $Q_{1}$ and $Q_{*}$. The Jacobian matrix of the system (32)
evaluated at $Q_{1}$ is

$J_{Q_{1}}=\begin{pmatrix}0&-1&0&0\\\
-\gamma(1+\beta)&-c&-\frac{1}{1+\alpha}&0\\\ 0&0&0&1\\\
0&0&\frac{\varepsilon}{d}(\frac{1}{1+\alpha}-\delta)&-\frac{c}{d}\end{pmatrix}$
(33)

The eigenvalues of the matrix $J_{Q_{1}}$ are
$\lambda_{1,2}=\frac{c}{2}\pm\frac{\sqrt{c^{2}+4g(1+\beta)}}{2}$ and
$\lambda_{3,4}=-\frac{c(1+\alpha)}{2}\pm\frac{\sqrt{\Gamma}}{d(1+\alpha)}$
where $\Gamma=(1+\alpha)^{2}c^{2}-4\varepsilon
d(1+\alpha)(1-\delta-\alpha\delta)$. First two eigenvalues are real, whereas
$\lambda_{3,4}$ are real for $c^{2}\geq\frac{4\varepsilon
d(1-\delta-\alpha\delta)}{1+\alpha}$. The traveling wave exist if all the
eigenvalues are real, otherwise the trajectories will spiral around $Q_{1}$
and leads to negative population density. Hence, the minimum speed of
traveling wave originating from predator free steady-state is

$c_{min}=\Big{[}\frac{4\varepsilon
d(1-\delta-\alpha\delta)}{1+\alpha}\Big{]}^{1/2}.$ (34)

Note that $c_{min}$ depends on $\varepsilon$, thus for $\varepsilon<1$, the
wave speed decreases. The minimum wave speed derived here is the same as it
was derived earlier with the help of the linearized equation. The expressions
for the eigenvalues of the Jacobian matrix $J_{Q_{*}}$ are quite complicated
and hence we avoid writing them here explicitly for the sake of brevity. For
$c\geq c_{min}$, the eigenvalues of $J_{Q_{*}}$ can be real or complex
depending on which we find monotone traveling wave and non-monotone as well as
periodic traveling wave originating from the steady states $E_{1}$.

(a) $t=300$

(b) $t=240$ Figure 7: Monotone traveling wave, non-monotone traveling wave and
periodic traveling wave obtained for the model (28) with $\alpha=0.5$,
$\beta=0.22$, $\gamma=3$, $\varepsilon=1$ and (a) $\delta=0.6$, (b)
$\delta=0.38$ and (c) $\delta=0.3$ at different time as mentioned with the
figures.

(c) $t=220$

To illustrate the existence of various type of traveling waves, we fix the
parameter values $\alpha=0.5$, $\beta=0.22$, $\gamma=3$, $\varepsilon=1$ and
consider $\delta$ as variable parameter. From (34) we find that traveling wave
exists for $\delta<2/3$. The eigenvalues of $J_{Q_{*}}$ are real for
$0.52<\delta<0.667$ and the eigenvalues are complex conjugate for $\delta\leq
0.52$. Complex conjugate eigenvalues with negative real parts correspond to
non-monotone traveling wave and with positive real part correspond to periodic
traveling wave. Three types of traveling waves are shown in Fig. 7(b) for
three different values of $\delta$. For $\delta\,=\,0.38$, the minimum wave
speed $c_{min}\approx 1.07$, and the complex conjugate eigenvalues with
negative real part are $-1.22\pm 0.423i$, whereas for $\delta=0.3,\
c_{min}\approx 1.211$ and the complex eigenvalues with positive real part are
$0.034\pm 0.405i$. The initial condition used in Fig. 7(b) is given by

$u(0,x)\,=\,\left\\{\begin{array}[]{ll}u_{*},&0\leq x\leq 3\\\ 1,&3<x\leq
300\\\
\end{array}\right.,\,\,\,v(0,x)\,=\,\left\\{\begin{array}[]{ll}v_{*},&0\leq
x\leq 3\\\ 0,&3<x\leq 300\\\ \end{array}\right..$

The traveling wave emerging from the above initial conditions connects the
prey-only steady state $E_{1}$ to the coexistence state $E_{*}$. Its shape
depends on parameter values; for $\delta=0.6$, its profile is monotone (see
Fig. 7(b)a). A non-monotone traveling wave emerges for $\delta=0.38$ (see Fig.
7(b)b). The range of spatio-temporal oscillation increases for values of
$\delta$ close to the temporal Hopf bifurcation threshold. For $\delta=0.3$,
we find periodic traveling wave, with a plateau behind the oscillatory front
corresponding to steady state $E_{*}$ which is, for these parameter values,
unstable: a phenomenon known as the dynamical stabilization [35, 54, 55, 71].

Figure 8: Periodic traveling waves for the model (28) with $\alpha=0.5$,
$\beta=0.22$, $\gamma=3$, $\varepsilon=1$ and (a) $\delta=0.3$, (b)
$\delta=0.2$ and (c) $\delta=0.1$ at the same instant of time $t=160$.

We mention here that the properties of emerging traveling wave are rather
robust with regard to the choice of initial conditions. For instance, for a
different initial condition as given by

$u(0,x)\,=\,\left\\{\begin{array}[]{ll}1,&0\leq x\leq 3\\\ 0,&3<x\leq 300\\\
\end{array}\right.,\,\,\,v(0,x)\,=\,\left\\{\begin{array}[]{ll}0.2,&0\leq
x\leq 2\\\ 0,&2<x\leq 300\\\ \end{array}\right.,$

the emerging periodic traveling waves shown in Fig. 7(b)c and Fig. 8a are
qualitatively similar. Numerical simulation for smaller values of $\delta$
shows an increase in the magnitude and period of the oscillating front as
shown in Fig. 8b-c. For numerical simulations we have chosen a spatial domain
of size $[0,300]$, further increase in domain size does not effect the
qualitative property of the traveling waves.

To understand the effect of slow-fast time scale on the resulting pattern
formation, here we consider the change in periodic traveling wave shown in
Fig. 8a for $\varepsilon<1$. The model (28) is simulated for three different
values of $\varepsilon$ as mentioned at the caption of Fig. 9. With the
decrease in $\varepsilon$, we observe that the oscillatory wake of the
invading predator front separating the predator-free area and the onset of
spatiotemporal shrinks and eventually disappears for smaller values of
$\varepsilon$, so that dynamical stabilization does not occur. Also the size
of predator-free area where prey exists at its carrying capacity increases
with the decrease in $\varepsilon$. The size of predator-free patch increases
further for $\varepsilon<0.25$.

Figure 9: Periodic traveling waves for the model (28) with $\alpha=0.5$,
$\beta=0.22$, $\gamma=3$, $\delta=0.3$ and (a) $\varepsilon=0.75$, (b)
$\varepsilon=0.5$ and (c) $\varepsilon=0.25$ at the same instant of time
$t=160$.

###  5.2 Patterns in two dimension

Finally, we consider the spatio-temporal pattern formation over two
dimensional spatial domain. The nonlinear reaction-diffusion system (27) is
solved numerically using five-point finite difference scheme for the Laplacian
operator and forward Euler scheme for the temporal part with the initial
conditions (35) and (36). Equal diffusivities are considered throughout i.e,
$d=1$. Two types of initial conditions are used to study the successful
invasion and establishment of both species. In [41] Morozov et.al used the
following initial condition which explains that a small amount of population
of both the species are introduced within a small elliptic domain

$u(0,x,y)\,=\,\left\\{\begin{array}[]{ll}u_{0},&\frac{(x-x_{1})^{2}}{\Delta_{11}}+\frac{(y-y_{1})^{2}}{\Delta_{12}}\,\leq\,1\\\
0,&\text{otherwise}\\\ \end{array}\right.\text{,}\ \
v(0,x,y)\,=\,\left\\{\begin{array}[]{ll}v_{0},&\frac{(x-x_{2})^{2}}{\Delta_{21}}+\frac{(y-y_{2})^{2}}{\Delta_{22}}\,\leq\,1\\\
0,&\text{otherwise}\\\ \end{array}\right.$ (35)

where $u_{0}$ and $v_{0}$ measure the initial densities of prey (native) and
predator (invasive) species respectively. The other initial condition we
consider here, is a small amplitude heterogeneous perturbation from the
homogeneous steady states (see [39] for details), is

$\displaystyle u(x,y,0)$ $\displaystyle=$ $\displaystyle
u_{*}-e_{1}(x-0.1y-225)(x-0.1y-675),$ (36) $\displaystyle v(x,y,0)$
$\displaystyle=$ $\displaystyle v_{*}-e_{2}(x-450)-e_{3}(y-450),$

where $(u_{*},v_{*})$ is the homogeneous steady state and for numerical
simulation we choose $e_{1}=2\times 10^{-7}$, $e_{2}=3\times 10^{-5}$, and
$e_{3}=2\times 10^{-4}$.

First we simulate the spatio-temporal model with the initial condition (35) in
a square domain $L\times L$ with $L=300$, with grid spacing $\Delta x=\Delta
y=1$ and time step $\Delta t=0.01$. The simulation results are verified with
other choices of $\Delta x$ and $\Delta t$ to ensure that the obtained results
are free from numerical artifact. We consider a small elliptic domain within
which the prey is at its carrying capacity ($u_{0}=1$) and a small amount of
population ($v_{0}=0.2$) is introduced. Other parameter values are
$x_{1}=153.5$, $y_{1}=145$, $x_{2}=150$, $y_{2}=150$, $\Delta_{11}=12.5$,
$\Delta_{12}=12.5$, $\Delta_{21}=5$, $\Delta_{22}=10$. We simulate the model
for a sufficiently long time so that the invading waves can cover the whole
domain and hits the domain boundary. Parameter values of $\alpha$, $\beta$,
$\gamma$ and $\varepsilon$ are same as mentioned in the previous subsection.
We will first check the change in resulting pattern by varying the parameter
$\delta$. In Fig 10, two dimensional spatial distribution of prey population
density at two different time points are shown for $\delta=0.3$. We have
omitted the spatial distribution of the predator population as they exhibit
similar patterns as exhibited by the prey species. Choosing $\delta=0.3$, just
below the temporal Hopf bifurcation threshold ($\delta_{H}=0.3768$), we
observe concentric circular rings as the initial transient pattern which
eventually settle down to an interacting spiral pattern once the transients
are over. With the advancement of time, the expanding circular rings hit the
domain boundary and break into irregular patches. These irregular spiral
patches cover the whole domain, and the system dynamics can be identified as
interacting spiral chaos (see Fig. 10b). Note that the initial invading waves
are the periodic traveling waves (Fig. 10a) but in large-time we observe
irregular spatio-temporal oscillations (Fig. 10b). This type of chaotic
dynamics persists in the vicinity of temporal Hopf bifurcation threshold
($0.2<\delta<\ \delta_{H}=0.3768$).

(a) $t=80$

(b) $t=500$

Figure 10: Spatio-temporal pattern of prey with initial condition 35 obtained
for $\alpha=0.5,\ \beta=0.22,\ \gamma=3,\ \delta=0.3,\ \varepsilon=1$ at
different intervals of time.

Keeping other parameters fixed, we further decrease $\delta\ (\leq 0.2)$ and
find propagating circular rings which are periodic traveling waves. The number
of rings, that is the number of population patches within the fixed domain,
decreases with the decrease in magnitude of $\delta$. These periodic traveling
waves do not break after hitting the boundary and the spatio-temporal dynamics
remains unaltered. This result is in agreement with Sherrat et. al [72] that
the system exhibits oscillatory dynamics as a successful invasion. Periodic
traveling fronts for $\delta\leq 0.2.$ are shown in Fig. 11(b).

(a) $\delta=0.2$

(b) $\delta=0.1$ Figure 11: Spatio-temporal pattern of prey with initial
condition 35 for $\alpha=0.5,\ \beta=0.22,\ \gamma=3,\ \varepsilon=1$ at
$t=200$.

(c) $\delta=0.05$

Now we consider the effect of the slow-fast time scale on the resulting
patterns. In the previous subsection, we have explained the reduction of
traveling wave speed with the decrease in the magnitude of $\varepsilon$. As a
result, the time taken by the predators to invade over the entire domain for
$\varepsilon\ll 1$ is much longer as compared to $\varepsilon=1$. For
$\varepsilon<1$ we find two kinds of distinctive changes in the resulting
patterns. The width of the population patches increase (see Fig 12a), and the
spatio-temporal chaotic dynamics changes to periodic temporal oscillation of
nearly homogeneous distribution of prey and predator densities (see Fig 12d).
The time evolution of the spatial average of the prey and predator population
is analogous to the temporal Canard cycle.

Figure 12: Spatio-temporal pattern of prey and plots of average spatial of
prey and predator population is given w.r.t time for $\alpha=0.5,\
\beta=0.22,\ \gamma=3,\ \delta=0.3$ for $\varepsilon=0.1$(left) and
$\varepsilon=0.01$(right). Upper panel shows the pattern at (a) $t=10000$, (b)
$t=5000$; lower panel shows the phase trajectory of spatially averaged
densities (c) $t\in[2000,10000]$, (d) $t\in[5000,10000].$

The choice of the initial condition and the domain size plays an important
role in pattern formation. We simulate the system (27) with the second kind of
initial condition (36) and over a square domain $L\times L$ with $L=900$. The
initial condition indicates a small heterogeneous perturbation from the
homogeneous steady state $(u_{*},v_{*})$. Choosing the same parameter set as
Fig. 10, initially, we find two spirals rotating about their fixed centers.
The regular spirals are destroyed with the advancement of time and the
interacting spiral pattern engulfs the whole domain (see Fig. 13b-c). These
patches move, break and form new patches, but qualitatively, the dynamics of
the system does not alter with time.

(a) $t=360$

(b) $t=900$

(c) $t=1500$

Figure 13: Spatio-temporal pattern of prey density for $\alpha=0.5,\
\beta=0.22,\ \gamma=3,\ \delta=0.3,\ \varepsilon=1,$ at different time
intervals.

Exhaustive numerical simulation indicates that for $\delta$ close to the
temporal Hopf bifurcation threshold, the system always exhibits spatio-
temporal chaos. The duration and type of transient patterns depend upon the
initial condition and the size of the domain. Now considering $\varepsilon<1$,
we found the persistent interacting spirals with thick arms (see Fig. 14a for
$\varepsilon=0.1$). The regular spiral grows in size and doesn’t breakdown
even after hitting the boundary when $\varepsilon$ is significantly small, say
$\varepsilon=0.01$. The irregularity of the temporal evolution of spatial
averages of both the population decreases and moves towards periodic or quasi-
periodic oscillation with the decrease in the magnitude of $\varepsilon$. This
claim is justified from the time evolution of spatial averages as presented in
the lower panel of Fig. 14.

(a) $\varepsilon=1$

(b) $\varepsilon=0.1$

(c) $\varepsilon=0.01$

(d) $\varepsilon=1$

(e) $\varepsilon=0.1$

(f) $\varepsilon=0.01$

Figure 14: Upper panel (a,b,c) represents spatial distribution of prey for
$\alpha=0.5,\ \beta=0.22,\ \gamma=3,\ \delta=0.3$ for different values of
$\varepsilon$ (a) $t=12000$, (b) $t=55000$, (c) $t=10000$ ; lower panel
(d,e,f) represents phase trajectory of spatially averaged densities obtained
after removing initial transients (d) $t\in[3000,12000],$ (e)
$t\in[3000,55000]$, (f) $t\in[5000,10000]$.

##  6 Discussion and Conclusions

Understanding the effects that the existence of multiple time scales may have
on the population dynamics of corresponding interacting species, in particular
by promoting or hampering their persistence, has been attracting an increasing
attention over the last two decades. In particular, some preliminary yet
significant work has been done to understand changes in the oscillatory
coexistence in the presence of slow-fast time scales [16, 26, 42, 61].
However, the effects of the slow-fast dynamics in the spatially explicit
systems, e.g. as given by the corresponding reaction–diffusion equations,
remains poorly investigated. This paper aims to bridge this gap, at least
partially. As our baseline system, we consider the classical Rosenzweig-
McCarthur prey-predator model with the multiplicative weak Allee effect in
prey’s growth. We pay particular attention to the interplay between the
strength of the weak Allee effect (quantified by parameter $0<\beta<1$) and
the difference in the time scales for prey and predator (quantified by
$\epsilon\leq 1$).

We first provide a detailed slow-fast analysis for the corresponding
nonspatial system. In doing that, we have obtained the following results:

  * •

in the presence of slow-fast dynamics ($\epsilon\ll 1$) and a weak Allee
effect, a decrease in the predator mortality may lead to a regime shift where
small-amplitudfe oscillations in the populaton abundance change to large-
amplitude oscillations (see Fig. 5). This change becomes more abrupt in case
the Allee effect is ‘not too weak’ (i.e. $\beta$ is sufficiently small), cf.
Figs. 5 and 6.

On a more technical side, we have derived an asymptotic expansion in
$\varepsilon$ for the invariant approximated manifolds and have explained the
dynamics of the system near the hyperbolic submanifolds. This theory cannot be
extended at non-hyperbolic points. To unravel the complete geometry of the
manifolds and their intersection as they pass through the non-hyperbolic
points we followed the blow-up technique [9, 27, 28]. We considered the slow-
fast normal form of the model by translating the fold point to the origin. As
the transformed system has a singularity at the origin, it is then blown up to
a sphere $S^{3}$ and the trajectories of the blow-up system are mapped on and
around the sphere. Using the blow-up analysis we have found the analytical
expression for the singular Hopf bifurcation curve
$(\mathcal{\lambda_{H}}(\sqrt{\varepsilon}))$ along which the eigenvalues
become singular as $\varepsilon\rightarrow 0$. A particular kind of slow-fast
solution known as canards (with or without head) has been found explicitly
with the help of Melnikov’s distance function in the blow-up space. We have
also calculated an analytical expression for the maximal canard curve
$(\lambda_{c}(\sqrt{\varepsilon}))$. Another type of periodic solution is
obtained which consists of two concatenated slow and fast flow, known as
relaxation oscillation. Analytically, we have proved the existence and
uniqueness of the relaxation oscillation cycle using the entry-exit function
[61, 82] and validated our results numerically.

The difference in the time scale for the growth and decay in prey and predator
species capture some interesting feature of respective populations. As the
prey population growth takes place over a faster time scale, the predator
population remains unchanged during the rapid growth and the decay of prey
population. On the other hand, the change in predator population occurs slowly
compared to the prey population. This type of growth and decay in two
constituent species is observed for steady-state coexistence as well as
oscillatory coexistence. The presence of weak Allee effect in prey growth acts
as a system saver. The size of the limiting relaxation oscillation cycle is
smaller in size when the magnitude of Allee effect is comparatively large (cf.
Fig. 6). It reduces the chance of extinction (maybe localized) as the periodic
attractor remains away from both the axes.

To understand the change in dynamic behavior, we have chosen the predator
mortality rate $\delta$, as the bifurcation parameter. For predator mortality
rate greater than Hopf threshold the system stabilizes at coexistence steady
state, whereas, the system shows oscillatory dynamics for mortality rate less
than the threshold. For the model under consideration, the Hopf threshold is
independent of the time scale parameter. But the combination of the mortality
rate along with time scale parameter ($0<\varepsilon\ll 1$) has enormous
effect on the nature of the oscillatory coexistence. For a fixed
$\varepsilon>0$, as we decrease $\delta$ below the Hopf threshold, we observe
a fast transition from small amplitude oscillatory coexistence to relaxation
oscillation within an exponentially small range of the parameter $\delta$ via
a family of canard cycles (Fig. 5), known as canard explosion. This type of
dynamics has been observed in an ecosystem where the growth rate of the
interacting species (resource-consumer type) differ on some orders of
magnitude. The reason behind this can be interpreted as: when the predator
density is at a maximum level, due to over-consumption the prey population
collapses rapidly. Since the predator is specialist, due to the lack of food
source the predator population starts decaying slowly and reaches a lower
density. It reduces grazing pressure on the prey, as a result the prey
population revives leading to a sudden outbreak. Again with increasing food
resources, the predator population increases slowly until it reaches a
desirable level which can be supported by the abundance of prey and thus the
cycle continues. Empirical evidence suggests that this type of oscillatory
dynamics is observed in real world, for example in the food web of Canadian
boreal forest [75] where outbreak of hare population follows a cycle of almost
11 years. In the aquatic ecosystem the seasonal cycle of Daphnia and algae
[66, 67] are observed and also in a forest ecosystem where insect pests
defoliate the adult trees [34].

We then considered the effect of multiple time scales in one-dimensional and
two-dimensional spatial extension of our slow-fast system. In the 1D case, the
minimum speed of the traveling wave of predator invading into the space
already occupied by its prey (observed in case of compact initial conditions)
is found analytically, while the patterns emerging in the wake of the front
are investigated by means of numerical simulations. In the 2D case, the effect
of the interplay between the weak Allee effect and the multiple timescales is
studied in simulations. The following result is worth of highlighting:

  * •

in the presence of a weak Allee effect, a decrease in the time scale ratio
(i.e. for $\epsilon\ll 1$) may lead to a regime shift where the pattern
becomes correlated across the whole spatial domain resulting in large-
amplitude oscillations of spatially average population density; see Fig. 14.
Since the corresponding trajectory in the phase plane ($<u>,<v>$) comes close
to the vertical axis, the immediate ecological implication of this is a likely
extinction of prey.

On a more technical side, our main interest was to study how the invasion of
the species is taking place and how it is getting affected with the
introduction of the time scale parameter. For the values of $\delta$ less than
Hopf bifurcation threshold ($\delta<\delta_{H}$), we find spatio-temporal
chaotic patterns. The onset of spatio-temporal chaos and the duration of
transient oscillation is completely influenced by the initial distribution of
the two species. Fig. 10 shows periodic traveling waves as transient dynamics
before spatio-temporal chaos sets in. However, small amplitude heterogeneous
perturbation around the homogeneous steady states reduces the time length for
transient dynamics and the system quickly enters spatio-temporal chaotic
regime. For $\delta$ significantly less than $\delta_{H}$ (see Fig. 11(b)) we
find only periodic traveling waves which indicate that continuous alteration
of population patches mimics the temporal dynamics of large amplitude
oscillations.

Consideration of time scale difference in the growth rates of prey and
predator have some stabilizing effect on the spatio-temporal pattern formation
scenario. On one hand, it increases the size of the coexisting population
patches over the domain, and on other hand, it drives the spatio-temporal
chaotic pattern to periodic or quasi-periodic oscillatory dynamics as shown in
Fig. 14. One prominent feature can be visualized from the numerical simulation
that spatio-temporal chaotic pattern engulf the entire domain of size
$900\times 900$ at $t=1700,$ for $\varepsilon=1,$ whereas it takes quite a
long time for $\varepsilon=0.1.$ The pattern obtained for $\varepsilon=0.1$ in
Fig. 14(b) is interacting spiral but not chaotic, as it does not show any
sensitivity to initial condition. The time evolution of average prey-predator
density changes from chaotic nature to quasi-periodic oscillation with the
decrease in magnitude of $\varepsilon$ as shown in the lower panel of Fig. 14.

In this work, we have considered a prey-predator model with specialist
predator and the consumption of prey by the predator follows prey-dependent
functional response. As a result, the periodic solution arising through Hopf
instability is stable and there is no possibility of global bifurcation
through which one or more species can collapse. The large amplitude
oscillatory coexistence obtained from the temporal model changes to periodic
traveling wave once the individuals are assumed to be distributed over their
habitat heterogeneously. The difference in the time-scale for the growth of
resource and consumer leads to the establishment of the species over larger
patches and the speed of invasive wave reduces with the decrease of
$\varepsilon$. The irregularity of population patches (spatio-temporal chaotic
patterns) as a part of successful invasion is observed for parameter values
close to the Hopf bifurcation threshold and $\varepsilon$ is close to or equal
to 1. Decrease in the magnitude of $\varepsilon$ reduces the irregular
oscillation but the duration of transient oscillations are enhanced. The study
of spatio-temporal pattern formation with a difference in time scales in the
context of ecological systems is quite unexplored in literature. This kind of
study can provide a better insight in establishment of the invasive species.
More realistic phenomena can be captured if we consider long food chain model
with multiple time scales and two species model with generalist predator with
predator dependent functional response, which we will study in our future
works.

## References

  * [1] Allesina, S., Bodini, A.: Who dominates whom in the ecosystem? Energy flow bottlenecks and cascading extinctions. Journal of Theoretical Biology 230, 351-358 (2004)
  * [2] Arditi, R., Ginzburg, L. R.: Coupling in predator-prey dynamics: Ratio-dependence. Journal of Theoretical Biology 139, 311-326 (1989)
  * [3] Baurmann, M ., Gross, T., Feudel, U.: Instabilities in spatially extended predator–prey systems: Spatio-temporal patterns in the neighborhood of Turing–Hopf bifurcations. Journal of Theoretical Biology 245, 220-229 (2007)
  * [4] Benoit, E., Callot, J.F., Diener, F., Diener, M.: Chasse au canard. Collectanea Mathematica 31-32, 37-119 (1981)
  * [5] Brown, J.H.: Complex ecological systems, in: Complexity: metaphors, models, and reality. (Cowan G.A., Pines D., Melzer D., Eds.), Santa Fe Institute studies in the Science of Complexity, proceedings, vol. XVIII. Reading MA: Addison-Wesley, 419-449 (1994)
  * [6] Courchamp, F., Berec, L., Gascoigne, J.: Allee Effects in Ecology and Conservation. Oxford University Press, UK, 2008
  * [7] Dumortier. F.: Singularities of Vector Fields. IMPA, Rio de Janeiro, Brazil, 1978
  * [8] Dumortier, F.: Techniques in the theory of local bifurcations: Blow-up, normal forms, nilpotent bifurcations, singular perturbations. In: Bifurcations and periodic orbits of vector fields, Springer, 19–73 (1993)
  * [9] Dumortier, F., Roussarie, R.: Canard cycles and center manifolds. Memoirs of the American Mathematical Society, 121, 577 (1996)
  * [10] Edwards, A.M., Brindley, J.: Zooplankton mortality and the dynamical behaviour of plankton population models. Bulletin of Mathematical Biology 61, 303-339 (1999)
  * [11] Fenichel, N.: Geometric Singular Perturbation Theory for Ordinary Differential Equations. Journal of Differential Equations 31, 53-98 (1979)
  * [12] Freedman, H.I.: Deterministic Mathematical Models in Population Ecology. Marcel Dekker, New York, 1980
  * [13] Gurney, W.S.C., Veitch, A.R., Cruickshank, I., McGeachin, G.: Circles and spirals: population persistence in a spatially explicit predator-prey model. Ecology 79, 2516-2530 (1998)
  * [14] Hastings, A., Harrison, S., McCann, K.: Unexpected spatial patterns in an insect outbreak match a predator diffusion model. Proceedings of The Royal Society B 264, 1837-1840 (1997)
  * [15] Hastings A, Powell T.: Chaos in a three-species food chain. Ecology 72, 896-903 (1991)
  * [16] Hek, G.: Geometric singular perturbation theory in biological practice. Journal of Mathematical Biology 60, 347–386 (2010)
  * [17] Hofbauer, J., Sigmund, K.: Evolutionary games and population dynamics. Cambridge University Press, UK, 1998
  * [18] Holling, C.S.: The functional response of predators to prey density and its role in mimicry and population regulation. Memoirs of the Entomological Society of Canada 97(S45), 5-60 (1965)
  * [19] Holt, R.D., Polis, G.A.: A theoretical framework for intraguild predation. American Naturalist 149, 745–764 (1997)
  * [20] Huang, J., Ruan, S., Song, J.: Bifurcations in a predator-prey system of Leslie type with generalized Holling type III functional response. Journal of Differential Equations 257, 1721-1752 (2014)
  * [21] Jankovic, M., Petrovskii, S.: Are time delays always destabilizing? Revisiting the role of time delays and the Allee effect. Theoretical Ecology 7, 335–349 (2014)
  * [22] Jansen, V.A.A.: Regulation of predator-prey systems through spatial interactions: a possible solution to the paradox of enrichment. Oikos 74, 384-390 (1995)
  * [23] Jiang, X., She, Z., Ruan, S.: Global dynamics of a predator-prey system with density-dependent mortality and ratio-dependent functional response. Discrete and Continuous Dynamical Systems B 26, 1967-1990 (2021)
  * [24] Jordan, F., Scheuring, I., Vida, G.: Species positions and extinction dynamics in simple food webs. Journal of Theoretical Biology 215, 441-448 (2002)
  * [25] Kareiva, P.M.: Population dynamics in spatially complex environments: theory and data. Philosophical Transactions of the Royal Society B 330, 175-190 (1990)
  * [26] Kooi, B.W., Poggiale, J.C.: Modelling, singular perturbation and bifurcation analyses of bitrophic food chains. Mathematical Biosciences 301, 93–110 (2018)
  * [27] Krupa, M., Szmolyan, P.: Extending Geometric Singular Perturbation Thoery to nonhyperbolic points- folds and canards in two dimension. SIAM Journal of Mathematical Analysis 33(2), 286–314 (2001)
  * [28] Krupa, M., Szmolyan, P.: Relaxation oscillation and Canard explosion. Journal of Differential Equations 174, 312-368 (2001)
  * [29] Kuehn, C.: Multiple Time Scale Dynamics. Springer, New York, 2015
  * [30] Kuznetsov, Y.A.: Elements of Applied Bifurcation Theory. Springer, New York, 2004
  * [31] Lewis, M.A., Kareiva, P.: Allee dynamics and the spread of invading organisms. Theoretical Population Biology 43, 141-158 (1993)
  * [32] Lewis, M.A., Pacala, S.: Modeling and analysis of stochastic invasion processes. Journal of Mathematical Biology 41, 387–429 (2000)
  * [33] Lewis, M.A., Petrovskii, S., Potts, J.R.: The Mathematics Behind Biological Invasions, Springer, New York, 2016
  * [34] Ludwig, D., Jones, D.D., Holling, C.S.: Qualitative analysis of insect outbreak systems: the spruce budworm and forest. Journal of Animal Ecology 47, 315–332 (1978)
  * [35] Malchow, H., Petrovskii, S.V.: Dynamical stabilization of an unstable equilibrium in chemical and biological systems. Mathematics and Computers in Modelling 36, 307-319 (2002)
  * [36] Matts Gyllenberg, personal communication.
  * [37] May, R.M.: Limit cycles in predator–prey communities. Science 177, 900-902 (1972)
  * [38] McCauley, E., Murdoch, W.W.: Predator–prey dynamics in environments rich and poor in nutrients. Nature 343, 455-457 (1990)
  * [39] Medvinsky, A., Petrovskii, S., Tikhonova., I., Malchow., H., Li, B.L.: Spatiotemporal complexity of plankton and fish dynamics. SIAM Review 44(3), 311–370 (2002)
  * [40] Mimura, M., Murray, J.D. On a diffusive prey-predator model which exhibits patchiness. Journal of Theoretical Biology 75, 249-262 (1978)
  * [41] Morozov., A., Petrovskii., S., Li., B.L.: Spatiotemporal complexity of patchy invasion in a predator-prey system with the Allee effect. Journal of Theoretical Biology 238, 18–35 (2006)
  * [42] Muratori, S., Rinaldi, S.: Remarks on competitive coexistence. SIAM Journal of Applied Mathematics 49(5), 1462-1472 (1989)
  * [43] Muratori, S., Rinaldi, S.: Low and High frequency oscillation in three-dimensional food chain systems. SIAM Journal of Applied Mathematics (1992)
  * [44] Murray, J.D.: Singular perturbations of a class of nonlinear hyperbolic and parabolic equations. Journal of Mathematics and Physics 47, 111–133 (1968)
  * [45] Murray, J.D.: Non-existence of wave solutions for a class of reaction diffusion equations given by the Volterra interacting-population equations with diffusion. Journal of Theoretical Biology 52(2), 459–469 (1975)
  * [46] Murray, J.D.: On travelling wave solutions in a model for the Belousov–Zhabotinskii reaction. Journal of Theoretical Biology 56(2), 329–353 (1976)
  * [47] Murray, J.D. A pre-pattern formation mechanism for animal coat marking. Journal of Theoretical Biology 88, 161-199 (1981)
  * [48] Murray, J.D.: Parameter space for Turing instability in reaction-diffusion mechanisms: A comparison of models. Journal of Theoretical Biology 98, 143-163 (1982)
  * [49] Murray, J.D. How the leopard gets its spots. Scientific American 258, 80-87 (1988)
  * [50] Murray, J.D.: Mathematical Biology. Springer, New York, 1989
  * [51] Pascual, M.: Diffusion-induced chaos in a spatial predator–prey system. Proceedings of The Royal Society B 251, 1-7 (1993)
  * [52] Petrovskii, S., Kawasaki, K., Takasu, F., Shigesada, N.: Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species. Japan Journal of Industrial and Applied Mathematics 18, 459-481 (2001)
  * [53] Petrovskii, S.V., Malchow, H.: A minimal model of pattern formation in a prey-predator system. Mathematics and Computers in Modelling 29, 49-63 (1999)
  * [54] Petrovskii, S.V., Malchow, H.: Critical phenomena in plankton communities: KISS model revisited. Nonlinear Analysis Real World Applications 1(1), 37-51 (2000)
  * [55] Petrovskii, S.V., Malchow, H.: Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theoretical Population Biology 59(2), 157-174 (2001)
  * [56] Petrovskii, S., Morozov, A., Venturino, E.: Allee effect makes possible patchy invasion in a predator–prey system. Ecology Letters 5, 345–352 (2002)
  * [57] Petrovskii, S., Vinogradov, M.E., Morozov, A.: Formation of the patchiness in the plankton horizontal distribution due to biological invasion in a two-species model with account for the Allee effect. Oceanology 42, 363–372 (2002)
  * [58] Piltz, S.H., Veerman, F., Maini, P.H., Porter, M.A.: A predator-2 prey fast-slow dynamical system for rapid predator evolution. SIAM Journal of Applied Dynamical Systems 16(1), 54-90 (2017)
  * [59] Poggiale, J.C., Aldebert, C., Girardot, B., Kooi, B.W.: Analysis of a predator–prey model with specific time scales: a geometrical approach proving the occurrence of canard solutions. Journal of Mathematical Biology 80, 39-60 (2020)
  * [60] Polis, G.A., Strong, D.R.: Food web complexity and community dynamics. American Naturalist 147, 813-846 (1996)
  * [61] Rinaldi, S., Muratori, S.: Slow-fast limit cycles in predator-prey models. Ecological Modelling 61, 287-308 (1992)
  * [62] Rodrigues, V.W., Mistro, D.C., Rodrigues, L.A.D.: Pattern Formation and Bistability in a Generalist Predator-Prey Model. Mathematics 8, 20 (2020)
  * [63] Rosenzweig, M.L.: Paradox of enrichment: destabilization of exploitation ecosystem in ecological time. Science 171, 385-387 (1971)
  * [64] Rosenzweig, M. L., MacArthur, R.: Graphical representation and stability conditions of predator-prey interaction. American Naturalist 97, 209-223 (1963)
  * [65] Saha, T., Pal, P.J., Banerjee, M.: Relaxation oscillation and canard explosion in a slow-fast predator-prey model with Beddington-DeAngelis functional response. Nonlinear Dynamics 103, 1195-1217 (2021)
  * [66] Scheffer, M., Rinaldi, S., Kuznetsov, Y.A.: Effects of fish on plankton dynamics: A theoretical analysis, Canadian Journal of Fisheries and Aquatic Sciences 57(6), 1208-1219 (2000)
  * [67] Scheffer, M., Rinaldi, S., Kuznetsov, Y.A., Van Nes E.H.: Seasonal dynamics of Daphnia and algae explained as a periodically forced predator–prey system. Oikos 80, 519–532 (1997)
  * [68] Segel, L.A., Jackson, J.L.: Dissipative structure: an explanation and an ecological example. Journal of Theoretical Biology 37, 545-559 (1972)
  * [69] Sen, D., Petrovskii, S., Ghorai, S., Banerjee, M.: Rich Bifurcation Structure of Prey-Predator Model Induced by the Allee Effect in the Growth of Generalist Predator. International Journal of Bifurcation Chaos 30, 2050084 (2020)
  * [70] Sen, M., Banerjee, M., Morozov, A.: Bifurcation analysis of a ratio-dependent prey-predator model with the Allee effect. Ecological Complexity 11, 12-27 (2011)
  * [71] Sherratt, J.A.: Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D 117(1-4), 145-166 (1998)
  * [72] Sherratt, J.A., Eagan, B.T., Lewis, M.A.: Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality? Philosophical Transactions of the Royal Society B 352, 21-38 (1997)
  * [73] Sherratt, J. A., Lewis, M.A., Fowler, A.: Ecological chaos in the wake of invasion. Proceedings of the National Academy of Sciences of the USA 92, 2524-2528 (1995)
  * [74] Siteur, K., Eppinga, M.B., Doelman, A., Siero, E., Rietkerk, M.: Ecosystems off track: rate-induced critical transitions in ecological models. Oikos 125, 1689–1699 (2016)
  * [75] Stenseth, N.C., Falck, W., Bjornstad, O.N., Krebs, C.J.: Population regulation in snowshoe hare and Canadian lynx: Asymmetric food web configurations between hare and lynx. Proceedings of the National Academy of Sciences of the USA 94(10), 5147-5152 (1997)
  * [76] Tu, L.W.: An Introduction to Manifolds. Springer, New York, 2008
  * [77] Turchin, P.: Complex population dynamics: a theoretical/empirical synthesis. Princeton University Press, Princeton, 2003
  * [78] Turing, A.M.: On the chemical basis of morphogenesis. Philosophical Transactions of the Royal Society B 237, 37-72 (1952)
  * [79] Van der Pol B.: On ”relaxation-oscillations”, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science Ser.7, 2, 978-992 (1926)
  * [80] Volterra, V.: Fluctuations in the abundance of a species considered mathematically. Nature 118, 558-560 (1926)
  * [81] Wang, C., Zhang, X.: Canards, heteroclinic and homoclinic orbits for a slow-fast predator-prey model of generalized Holling type IIIs. Journal of Differential Equations 267, 3397–3441 (2019)
  * [82] Wang, C., Zhang, X.: Relaxation oscillations in a slow–fast modified Leslie–Gower model. Applied Mathematics Letters 87, 147–153 (2019)
  * [83] Zou, R., Guo, S.: Dynamics of a diffusive Leslie-Gower predator-prey model in spatially heterogeneous environment. Discrete and Continuous Dynamical Systems B 25, 4189-4210 (2020)

## Appendix A

Here we will follow geometric singular perturbation technique as given by
Fenichel [11] to find the analytical expression of locally perturbed invariant
manifold $C^{1}_{\varepsilon}$. Since $v=q(u,\varepsilon)$, from the
invariance condition we have

$\dfrac{dv}{dt}=\dfrac{dq(u,\varepsilon)}{du}\dfrac{du}{dt}.$

Using the explicit expression for $\dfrac{du}{dt}$ and $\dfrac{dv}{dt}$ from
(4) we get

$\displaystyle\varepsilon
q(u,\varepsilon)(u(1-\alpha\delta)-\delta)=u\dfrac{dq(u,\varepsilon)}{du}(\gamma(1-u)(u+\beta)(1+\alpha
u)-q(u,\varepsilon)).$ (37)

Substituting the asymptotic expansion of $q(u,\varepsilon)$ from (11 and
assuming $u\neq 0$, $\dot{q_{0}}(u)\neq 0,$ we equate $\varepsilon$ free terms
from both sides to obtain

$\displaystyle q_{0}(u)$ $\displaystyle=\gamma(1-u)(u+\beta)(1+\alpha u),$
(38)

which is exactly the critical manifold. Now equating the coefficients of
$\varepsilon$ from both sides of (37) we get

$\displaystyle
q_{1}(u)=\dfrac{q_{0}(u)(u(1-\alpha\delta)-\delta)}{-u\dot{q}_{0}(u)}.$ (39)

Similarly we obtain $q_{2}(u)$ by equating the coefficients of
$\varepsilon^{2},$

$\displaystyle
q_{2}(u)=\dfrac{q_{1}(u)(u(1-\alpha\delta)-\delta)+uq_{1}\dot{q_{1}}(u)}{-u\dot{q_{0}}(u)}.$
(40)

Proceeding as above we find $q_{r}(u),\ r=3,4,\cdots$ by equating the
coefficients of $\varepsilon^{r}$ from (37). Therefore the second order
approximation of the perturbed invariant manifold is given by

$q(u,\varepsilon)=q_{0}(u)+\varepsilon q_{1}(u)+\varepsilon^{2}q_{2}(u),$

where $q_{0},q_{1},q_{2}$ are given in (38)-(40).

## Appendix B

We apply the blow-up transformation in the slow-fast normal form (14) where

$h_{1}(U,V)=u_{*}+U,\ \ h_{3}(U,V)=0,\ \ h_{5}(U,V)=(v_{*}+V)(1+\alpha
u_{*})+Uv_{*}\alpha,$
$h_{2}(U,V)=-\gamma(-1+6u_{*}^{2}\alpha+3u_{*}(1+\alpha(\beta-1))+\beta-\alpha\beta)-U\gamma(1+\alpha(4u_{*}+\beta-1)),$
$h_{4}(U,V)=(v_{*}+V)(1-\alpha\delta_{*}),\ \
h_{6}(U,V)=u_{*}-(1+u_{*}\alpha)\delta_{*},,$

On chart $K_{2}$, $\bar{\varepsilon}=1$ so the blow-up transformation as
defined in (15) reduces to:

$\bar{r}=\sqrt{\varepsilon},\ U=\sqrt{\varepsilon}\bar{U},\
V=\varepsilon\bar{V},\ \lambda=\sqrt{\varepsilon}\bar{\lambda}.$ (41)

Using the transformation (41) we can write the system (17) by removing the
overbars as

$\displaystyle U_{t}$
$\displaystyle=-b_{1}V+b_{2}U^{2}+\sqrt{\varepsilon}\mathcal{G}_{1}(U,V)+O(\sqrt{\varepsilon}(\lambda+\sqrt{\varepsilon})),$
(42) $\displaystyle V_{t}$
$\displaystyle=b_{3}U-b_{4}\lambda+\sqrt{\varepsilon}\mathcal{G}_{2}(U,V)+O(\sqrt{\varepsilon}(\lambda+\sqrt{\varepsilon})),$

where

$\displaystyle b_{1}=u_{*},\ \
b_{2}=-\gamma(-1+6u_{*}^{2}\alpha+3u_{*}(1+\alpha(\beta-1))+\beta-\alpha\beta),$
(43) $\displaystyle b_{3}=v_{*}(1-\alpha\delta_{*}),\ \ b_{4}=v_{*}(1+\alpha
u_{*}),$

and

$\displaystyle\mathcal{G}_{1}(U,V)=a_{1}U-a_{2}UV+a_{3}U^{3},\ \
\mathcal{G}_{2}(u,V)=a_{4}U^{2}+a_{5}V.$ (44)

Let the equilibrium point of the system (42) is $(U_{e},V_{e})$,
$U_{e}=\dfrac{b_{4}\lambda}{b_{3}}+O(2)$ and $V_{e}=O(2)$ where
$O(2):=O(\lambda^{2},\lambda\sqrt{\varepsilon},\lambda)$. Linearizing the
system about this equilibrium point we have the Jacobian matrix as

$\mathcal{J}:=\begin{pmatrix}2U_{e}b_{2}+a_{1}\sqrt{\varepsilon}+O(2)&-b_{1}+O(2)\\\
b_{3}+O(2)&a_{5}\sqrt{\varepsilon}+O(2)\end{pmatrix}$ (45)

At the Hopf bifurcation we have Trace $\mathcal{J}=0$ which implies

$\dfrac{2b_{2}b_{4}\lambda}{b_{3}}+\sqrt{\varepsilon}(a_{1}+a_{5})+O(2)=0.$
(46)

and applying the blow-down map
$\mathcal{\lambda_{H}}=\lambda\sqrt{\varepsilon}$ we get the singular Hopf
bifurcation curve $\mathcal{\lambda_{H}}(\sqrt{\varepsilon})$ for the slow-
fast normal form (14) as

$\mathcal{\lambda_{H}}(\sqrt{\varepsilon})=-\dfrac{b_{3}(a_{1}+a_{5})}{2b_{2}b_{4}}\varepsilon+O(\varepsilon^{3/2}).$
(47)

## Appendix C

Here we prove the existence of maximal canard curve and will give an
analytical expression for the same. For that we will first prove the following
proposition. In chart $K_{2}$ of the blow up space we consider the
desingularized system (42) as

$\displaystyle U_{t}$
$\displaystyle=-b_{1}V+b_{2}U^{2}+r\mathcal{G}_{1}(U,V)+O(\lambda r,r),$ (48)
$\displaystyle V_{t}$
$\displaystyle=b_{3}U-b_{4}\lambda+r\mathcal{G}_{2}(U,V)+O(\lambda r,r),$
$\displaystyle r_{t}$ $\displaystyle=0,$ $\displaystyle\lambda_{t}$
$\displaystyle=0,$

where $b_{1},\ b_{2},\ b_{3},\ b_{4},\mathcal{G}_{1}$ and $\mathcal{G}_{2}$
are computed above (43), (44). The dynamics of the system on the sphere is
obtained by putting $r=0$ in (48) for different values of $\lambda$ in the
vicinity of $0$. Thus, by taking $r=0,\ \lambda=0$, the above system is
integrable and we have

$\displaystyle U_{t}$ $\displaystyle=-b_{1}V+b_{2}U^{2},$ (49) $\displaystyle
V_{t}$ $\displaystyle=b_{3}U.$

This is a Riccati equation and the solution of this equation helps in proving
our main theorem.
Proposition 1 The solution of the system (49) is given by $H(U,V)=c,$ where

$H(U,V)=e^{-\dfrac{2b_{2}}{b_{3}}V}\Big{(}\dfrac{b_{3}}{2}U^{2}-\dfrac{b_{1}b_{3}^{2}}{4b_{2}^{2}}-\dfrac{b_{1}b_{3}}{2b_{2}}V\Big{)}$

and

$\dfrac{dU}{dt}=-e^{\dfrac{2b_{2}}{b_{3}}V}\dfrac{\partial H}{\partial V},$
$\displaystyle\dfrac{dV}{dt}=e^{\dfrac{2b_{2}}{b_{3}}V}\dfrac{\partial
H}{\partial U}.$ (50)

###### Proof.

We can write the above Riccati system (49) as

$\dfrac{dV}{dU}=\dfrac{b_{3}U}{-b_{1}V+b_{2}U^{2}}$ (51)

where the integrating factor is $e^{-\dfrac{2b_{2}}{b_{3}}V}$. Multiplying
both sides with the I.F and integrating we get

$e^{-\dfrac{2b_{2}}{b_{3}}V}\Big{(}U^{2}-\dfrac{b_{1}}{b_{2}}V-\dfrac{b_{1}b_{3}}{2b_{2}^{2}}\Big{)}=c_{0}.$

Multiplying with $\dfrac{b_{3}}{2}$ we obtain the solution of the system( 49)
as

$e^{-\dfrac{2b_{2}}{b_{3}}V}\Big{(}\dfrac{b_{3}}{2}U^{2}-\dfrac{b_{1}b_{3}^{2}}{4b_{2}^{2}}-\dfrac{b_{1}b_{3}}{2b_{2}}V\Big{)}=c,$

where $c=c_{0}\dfrac{b_{3}}{2}$ is a constant. The solution determined by
$c=0$ is a parabola of the form

$U^{2}=\dfrac{b_{1}b_{3}}{2b_{2}^{2}}+\dfrac{b_{1}}{b_{2}}V.$

∎

Proof of theorem 4.2: We write the solution of the system (49) in the
parametric form

$\eta(t)=(U(t),V(t))=\Big{(}t,\dfrac{b_{2}}{b_{1}}t^{2}-\dfrac{b_{3}}{2b_{2}}\Big{)},\
t\in\mathbb{R}$ (52)

For $\varepsilon=0$ the attracting and repelling submanifolds of the critical
manifold $\mathcal{M}^{1}_{0}$ intersect along the equator of the blow-up
space $S^{3}$. From Fenichel’s theory, for $\varepsilon>0$ there exist
invariant perturbed attracting $(\mathcal{M}_{\varepsilon}^{1,a})$ and
repelling submanifold $(\mathcal{M}_{\varepsilon}^{1,r})$. Along the curve
(52), the attracting $(\mathcal{M}_{\varepsilon}^{1,a})$ and repelling
$(\mathcal{M}_{\varepsilon}^{1,r})$ invariant submanifolds in the blow-up
space intersect and the solution trajectory lying in that intersection is
called maximal canard. We use Melnikov function to calculate the distance
between these invariant manifolds [27], [29] which is given by

$D_{r,\lambda}=d_{r}r+d_{\lambda}\lambda+O(r^{2}),$ (53)

where

$\displaystyle d_{r}=\int_{-\infty}^{\infty}\nabla
H(\eta(t))^{T}\mathcal{G}(\eta(t))dt,$ (54) $\displaystyle
d_{\lambda}=\int_{-\infty}^{\infty}\nabla H(\eta(t))^{T}\begin{pmatrix}0\\\
-b_{4}\end{pmatrix}dt,$

where $\mathcal{G},\ H$ and $b_{4}$ are defined in (44), (50) and (43)
respectively. The distance between the submanifolds
$\mathcal{M}_{\varepsilon}^{1,a}$ and $\mathcal{M}_{\varepsilon}^{1,r}$ is
given by the eq. (53). And since the maximal canard lie in the intersection of
these manifolds, so we must have $D_{r,\lambda}=0$. For that we now calculate
the Melnikov-type integrals $d_{r}$ and $d_{\lambda}$ (52 ) and (54).
Therefore,

$\displaystyle d_{r}$
$\displaystyle=\int_{-\infty}^{\infty}\Big{[}(a_{1}U-a_{2}UV+a_{3}U^{3})\dfrac{\partial
H(\eta(t))}{\partial U}+(a_{4}U^{2}+a_{5}V)\dfrac{\partial
H(\eta(t))}{\partial V}\Big{]}dt$ (55)
$\displaystyle=\int_{-\infty}^{\infty}e^{-\dfrac{2b_{2}}{b_{3}}V}\Big{[}(a_{1}U-a_{2}UV+a_{3}U^{3})b_{3}U+(a_{4}U^{2}+a_{5}V)(b_{1}V-b_{2}U^{2})\Big{]}dt$
$\displaystyle=e\int_{-\infty}^{\infty}e^{-A_{4}t^{2}}\Big{(}A_{1}t^{4}+A_{2}t^{2}+A_{3}\Big{)}dt$

where,

$A_{1}=a_{3}b_{3}-\dfrac{a_{2}b_{2}b_{3}}{b_{1}},\
A_{2}=a_{1}b_{3}+\dfrac{a_{2}b_{3}^{2}}{2b_{2}}-\dfrac{a_{4}b_{1}b_{3}}{2b_{2}}-\dfrac{a_{5}b_{3}}{2},\
A_{3}=\dfrac{a_{5}b_{1}b_{3}^{2}}{4b_{2}^{2}},\
A_{4}=\dfrac{2b_{2}^{2}}{b_{1}b_{3}}.$

Now substituting $z=t^{2}$ and by repeated integration by parts we obtain

$\displaystyle
d_{r}=e\Big{(}\dfrac{3A_{1}}{4A_{4}^{2}}+\dfrac{A_{2}}{2A_{4}}+A_{3}\Big{)}\int_{-\infty}^{\infty}e^{-A_{4}t^{2}}dt,$
(56)

and

$\displaystyle d_{\lambda}$
$\displaystyle=-\int_{-\infty}^{\infty}b_{4}\dfrac{\partial H}{\partial V}dt$
(57)
$\displaystyle=b_{4}\int_{-\infty}^{\infty}e^{-\dfrac{2b_{2}}{b_{3}}V}(-b_{1}V+b_{2}U^{2})dt$
$\displaystyle=eA_{5}\int_{-\infty}^{\infty}e^{-A_{4}t^{2}}dt,$

where $A_{5}=\dfrac{b_{1}b_{3}b_{4}}{2b_{2}}$. Since $d_{\lambda}\neq 0$
therefore using implicit function theorem we can explicitly solve for
$\lambda$ from (53)

$\displaystyle\lambda(r)$
$\displaystyle=-\dfrac{d_{r}}{d_{\lambda}}r+O(r^{2})=-\dfrac{1}{A_{5}}\Big{(}\dfrac{3A_{1}}{4A_{4}^{2}}+\dfrac{A_{2}}{2A_{4}}+A_{3}\Big{)}r+O(r^{2}).$
(58)

Now using blow down map $\lambda_{c}=\lambda\sqrt{\varepsilon}$ we obtain the
maximal canard curve for the slow-fast normal form (14).

$\displaystyle\lambda_{c}(\sqrt{\varepsilon})$
$\displaystyle=-\dfrac{1}{A_{5}}\Big{(}\dfrac{3A_{1}}{4A_{4}^{2}}+\dfrac{A_{2}}{2A_{4}}+A_{3}\Big{)}\varepsilon+O(\varepsilon^{3/2}).$
(59)

## Appendix D

Here we prove the existence of a unique attracting limit cycle called
relaxation oscillation. To study the dynamics of the system (23) we define two
section of the flow as

$\displaystyle\Delta^{in}=\\{(u_{+},v):u_{+}<<u_{max},v\in(v_{1}-\rho,v_{1}+\rho)\\},$
$\displaystyle\Delta^{out}=\\{(u_{+},v):u_{+}<<u_{max},v\in(v_{0}-\rho^{2},v_{0}+\rho^{2})\\},$

where $u_{max},\ v_{1},\ v_{0}$ are defined in subsection 4.3 and $\rho$ is
sufficiently small positive number.

Let us define a return map $\Pi:\Delta^{in}\rightarrow\Delta^{in}$ which is a
composition of two maps

$\Phi:\Delta^{in}\rightarrow\Delta^{out},\ \
\Psi:\Delta^{out}\rightarrow\Delta^{in},$

such that $\Pi=\Psi\circ\Phi$. Let us fix $\varepsilon>0$ and we take a point
$(u_{+},v_{+})$ on the section $\Delta^{in}$. Now we consider a trajectory of
the system (23) starting from the initial point $(u_{+},v_{+})$. From the
analysis of the entry-exit function we can say that this trajectory will be
attracted to $V_{+}$ and will leave $V_{-}$ at point $(0,p(v_{+})),$ where $p$
is the entry-exit function. The trajectory then jumps into the section
$\Delta^{out}$ at the point $(u_{+},p(v_{+})).$ Thus, the map $\Phi$ is
defined with the help of entry-exit function as
$\Phi(u_{+},v_{+})=(u_{+},p(v_{+})).$
Now to study the map $\Psi$ we consider two trajectories
$\gamma_{\varepsilon}^{1},\gamma_{\varepsilon}^{2}$ starting from the section
$\Delta^{out}$. These trajectories get attracted toward
$C_{\varepsilon}^{1,a}$ where the slow flow is given by
$\dfrac{du}{d\tau}=\dfrac{g(u,q(u,\varepsilon))}{\dot{q}(u,\varepsilon)}.$

They follow the slow perturbed manifold until the vicinity of the fold point
where they contract exponentially toward each other [82] and jump into
$\Delta^{in}.$ From Theorem 2.1 of [27] we have that the map $\Pi$ is a
contraction. Using contraction mapping theorem we conclude that $\Pi$ has a
unique fixed point which gives rise to a unique relaxation oscillation cycle
$\gamma_{\varepsilon}$. Further from Fenichel’s theory we infer that
$\gamma_{\varepsilon}$ converges to $\gamma_{0}$ as $\varepsilon\rightarrow
0.$

Now for the parameter values $\alpha=0.5,\ \beta=0.2,\ \delta=0.3,$ the unique
attracting cycle $\gamma_{\varepsilon}$ for $\varepsilon=0.1$, is shown below
which converges to $\gamma_{0}$ as $\varepsilon\rightarrow 0.$

Figure 15: Singular trajectory $\gamma_{0}$ (blue) and unique attracting limit
cycle $\gamma_{\varepsilon}$ for $\varepsilon=0.1$ (green) for $\alpha=0.5,\
\beta=0.2,\ \delta=0.3.$