Long time dynamics of a three-species food chain model with Allee effect in the top predator

Abstract

The Allee effect is an important phenomenon in population biology characterized by positive density dependence, that is a positive correlation between population density and individual fitness. However, the effect is not well studied in multi-level trophic food chains. We consider a ratio dependent spatially explicit three species food chain model, where the top predator is subjected to a strong Allee effect. We show the existence of a global attractor for the the model, that is upper semicontinuous in the Allee threshold parameter $m$ . To the best of our knowledge this is the first robustness result, for a spatially explicit three species food chain model with an Allee effect. Next, we numerically investigate the decay rate to a target attractor, that is when $m=0$ , in terms of $m$ . We find decay estimates that are $\mathcal{O}(m^{\gamma})$ , where $\gamma$ is found explicitly. Furthermore, we prove various overexploitation theorems for the food chain model, showing that overexploitation has to be driven by the middle predator. In particular overexploitation is not possible without an Allee effect in place. We also uncover a rich class of Turing patterns in the model which depend significantly on the Allee threshold parameter $m$ . Our results have potential applications to trophic cascade control, conservation efforts in food chains, as well as Allee mediated biological control.

Rana D. Parshad, Emmanuel Quansah and Kelly Black

Department of Mathematics,

Clarkson University,

Potsdam, New York 13699, USA.

Ranjit K. Upadhyay and S.K.Tiwari

Department of Applied Mathematics,

Indian School of Mines,

Dhanbad 826004, Jharkhand, India.

Nitu Kumari

School of Basic Sciences,

Indian Institute of Technology Mandi,

Mandi, Himachal Pradesh 175 001, India.

Keywords:

three species reaction diffusion food chain model, global existence, global attractor, upper semi-continuity, Allee effect, Turing instability .

1 Introduction

Interactions of predator and prey species are ubiquitous in spatial ecology. Therein a predator or a “hunting” organism, hunts down and attempts to kill a prey, in order to feed. Food webs, which comprise all of the predator-prey interactions in a given ecosystem, are inherently more complex. The food chains or linear links of these food webs in real ecosystems, have multiple levels of predator-prey interaction, across various trophic levels. To better understand natural ecosystems, with multiple levels of trophic interaction, where predators and prey alike, disperse in space in search of food, mates and refuge, a natural starting point is to deviate from the classical predator-prey two species models, and investigate spatially explicit three species food chains. Such models are appropriate to model populations of generalist predators, predating on a specialist predator, which in turn predates on a prey species. Also, they can model two specialist predators competing for a single prey [15, 30]. The spatial component of ecological interaction has been identified as an important factor in how ecological communities are shaped, and thus the effects of space and spatially dispersing populations, via partial differential equations (PDE)/spatially explicit models of three and more interacting species, have been very well studied [33, 15, 16, 23, 19]. In most of these models, the prey is regulated or “inhibited” from growing to carrying capacity due to predation by the predator. Whereas loss in the predator is due to death or intraspecific competition terms. However, there are various other natural self regulating mechanisms in a population of predators (or prey). For instance, much research in two species models, has focused on one such mechanism: the so called Allee effect. However, less attention has been paid to this mechanism in the three species case. This effect, named after the ecologist Walter Clyde Allee, can occur whenever fitness of an individual in a small or sparse population, decreases as the population size or density does [8]. Since the pioneering work of Allee [2, 3], Allee dynamics has been regarded as one of the central issues in the population and community ecology.

The effect can be best understood by the following equation for a single species $u$ ,

\frac{du}{dt}=u(u-m)(1-u/K),

(1)

essentially a modification to the logistic equation. Here $u(t)$ the state variable, represents the numbers of a certain species at a given time $t$ , $K$ is the carrying capacity of the environment that $u$ resides in, and $m$ is the Allee threshold, with $m<K$ . A strong Allee effect occurs if $m>0$ , [31]. This essentially means that if $u$ falls below the threshold population $m$ , its growth rate is negative, and the species will go extinct [17]. If $m<0$ , then there is a weak Allee effect, and we have a compensatory growth function [13, 5, 7]. For our purposes we assume a strong Allee effect is in place, that is $m>0$ . Note, $u^{*}=0,K$ are stable fixed points for (1), and $u^{*}=m$ is unstable. Thus dynamically speaking the global attractor, which is the repository of all the long time dynamics of for (1) is $[0,K]$ . An interesting question can now be posed: what happens to this attractor as $m\rightarrow 0$ ? Ecologically speaking, this is asking: what happens to the species $u$ as the Allee threshold $m$ is decreased? When $m=0$ the only stable fixed point is $u^{*}=K$ , with $u^{*}=0$ being half stable. Thus the global attractor is now reduced to a single point $K$ . What we observe is that the difference between having a slight Allee effect (say $0<m<<1$ ) and having no Allee effect ( $m=0$ ), can change completely what the global attractor of (1) is. That is to say, the global attractor changes from an entire set $[0,K]$ , to a single fixed point $K$ . From an ecological point of view, this says as long as there is a slight Allee effect there is an extinction risk for $u$ , but without one there is none, and $u$ will always grow to carrying capacity. Proving that an attractor $\mathcal{A}_{m}$ approaches a target attractor $\mathcal{A}_{0}$ , in a continuous sense, in the case of spatially explicit/PDE models requires a fair amount of work. This involves making several estimates of functional norms, independent of the parameter $m$ [22]. If this can be proven however, the attractor $\mathcal{A}_{m}$ is said to be robust at $m=0$ . However, to the best of our knowledge there are no robustness results in three species food chain models, where the parameter of interest is the Allee threshold $m$ . Note, unless a systems dynamics are robust, there is no possibility to capture the same in a laboratory experiment or natural setting.

The critical importance of the Allee effect has widely been realized in the conservation biology: That is: it is most likely to increase the extinction risk of low-density populations. This is also known as critical depensation in fisheries sciences. Essentially individual species population must surpass the threshold $m$ to grow. Hence a small introduced population under strong Allee effect can only succeed if it is faced with favorable ecological conditions, or experiences rapid adaptive evolution, or simply has good luck, so that it can surpass this critical threshold. There are several factors that might cause an Allee effect. These include difficulty of finding mates at low population densities, inbreeding depression and environmental conditioning[28]. The Allee effect can be regarded not only as a suite of problems associated with rarity, but also as the basis of animal sociality [26]. Petrovskii et al. [21] found that a deterministic system with Allee effect, can induce patch invasion. Morozov et al. [14] found that the temporal population oscillations can exhibit chaotic dynamics even when the distribution of the species in the space was regular. Also, Sharma and Samanta [24] developed a ratio-dependent predator-prey model with disease in prey, as well as an Allee effect in prey. Furthermore, Invasion biologists have recently attempted to consider Allee effects as a benefit in limiting establishment of an invading species[28], which without the effect would possibly run amuck, an cause excessive damages to native species because of predation and competition.

On the same lines, a phenomenon that occurs in large food chains, is the excessive harvesting/predation of certain species in the chain, via the predators in the trophic level above them. In many aquatic food chains, this can lead to trophic cascades [6], and is among the major activities that is threatening global biodiversity. Formally, overexploitation refers to the phenomenon where excessive harvesting of a species can result in its extinction. Mathematically, for a two species predator-prey system, one can prove an overexploitation type theorem, if one shows for a large enough initial density of the predator, the prey will be predated on till extinction. In this case, If the predator does not have an alternate food source, it will also subsequently go extinct. Although there is a fair amount of literature on this in two species models [32], it is less studied in three species models. In particular, the Allee effect itself on overexploitation, in multi-trophic level food chains, has not been sufficiently explored.

Our primary goal in the current manuscript is to propose and analyze a reaction-diffusion three species food chain model, with an Allee effect in the top predator. In particular, we aim to investigate the link between the Allee threshold parameter $m$ , and the longtime dynamics of the three species food chain, in terms of

•

Robustness of global attractors in terms of the Allee threshold parameter $m$ .
•

Overexploitation phenomenon in the food chain model, as it is effected by the Allee threshold parameter $m$ .
•

Pattern formation and the effect on patterns that form in the food chain model due to the Allee threshold parameter $m$ .

Our hope is that this model could be used as a feasible toy model, to better understand the following in multi-trophic level food chains,

•

Top-down pressure in multi-trophic level food chains leading to overexploitation phenomenon.
•

Trophic cascades and cascade control via an Allee effect.
•

Conservation efforts in food webs mediated via an Allee effect.
•

Biological control of an invasive top predator in a food chain via Allee effect

Thus, we consider a situation where a prey species $u$ serves as the only food for a specialist predator $v$ which is itself predated by a generalist top predator $r$ .
This is a typical situation often seen in nature in various food chains [20]. The governing equations for populations $u$ and $v$ follow ratio dependent functional responses, and are modeled by the Volterra scheme i.e., the middle predator population dies out exponentially in the absence of its prey. There is recently much debate about functional responses used in ecology, and the ratio dependent response is considered to be more realistic than its Holling type counterparts [1]. The above situation is described via the following system of PDE. This is already nondimensionalised, see appendix 7.1 for the details of the nondimensionalisation.

	$\displaystyle\frac{\partial u}{\partial t}=d_{1}\Delta u+u-u^{2}-w_{1}\frac{uv}{u+v},$		(2)
	$\displaystyle\frac{\partial v}{\partial t}=d_{2}\Delta v-a_{2}v+w_{2}\frac{uv}{u+v}-w_{3}\left(\frac{vr}{v+r}\right),$		(3)
	$\displaystyle\frac{\partial r}{\partial t}=d_{3}\Delta r+r\big{(}{{r-m}}\big{)}\bigg{(}c-\frac{w_{4}r}{v+D_{3}}\bigg{)}.,$		(4)

the spatial domain for the above is $\Omega\subset\mathbb{R}^{n}$ , $n=1,2,3$ . $\Omega$ is assumed bounded, and we prescribe Neuman boundary conditions $\nabla u\cdot\textbf{n}=\nabla v\cdot\textbf{n}=\nabla r\cdot\textbf{n}=0$ . In the above $a_{2}$ is the intrinsic death rate of the specialist predator $v$ in the absence of its only food $u$ , $c$ measures the rate of self-reproduction of the generalist predator $r$ . $w_{i}^{\prime}s$ are the maximum values which per capita growth rate can attain. $D_{3}$ shows that the top predator $r$ is a generalist and can switch its food source in the absence of $v$ . Also we model the top predator via a Leslie Gower scheme [30], and assume that it is subject to the Allee effect. The parameter $m$ is the Allee threshold or minimum of viable population level. In this work, we limit ourselves to the case where $m>0$ , that is we assume a strong Allee effect.

In particular we ask: How does the Allee threshold $m$ affect the various dynamical aspects of model system (2)-(4)? To this end we list our primary findings.

•

We show that there is a $(L^{2}(\Omega),H^{1}(\Omega))$ global attractor for (2)-(4), via theorem 3.2.
•

We show that this attractor is upper semi-continuous w.r.t the Allee threshold $m$ , via theorem 3.3. That is the attractor $\mathcal{A}_{m}$ is robust at $m=0$ . This result requires making estimates of various functional norms, independent of the parameter $m$ , and then proving and using a chain of theorems that facilitate passing to the limit as $m\rightarrow 0$ . The calculations are detailed, so we confine them to an appendix section 7.
•

The decay rate of the global attractor $\mathcal{A}_{m},m>0$ , for the system (2)-(4), to a target attractor $\mathcal{A}_{0}$ (that is when $m=0$ ) is estimated in terms of $m$ numerically. We find a decay rate of the order $\mathcal{O}(m^{\gamma})$ , where $\gamma$ is close to 1.
•

We investigate overexploitation phenomenon in the model system (2)-(4), showing that the Allee threshold $m$ can effect overexploitation in the system via theorems 4.1, 4.3 and lemma 4.5. In particular overexploitation is not possible without an Allee effect in place.
•

Turing instability exists in the model system (2)-(4), via theorem 5.3. Furthermore, the Allee threshold $m$ has a significant impact on the type of Turing patterns that form, see figures 3, 4, 5, 6, 7, 9.

2 Preliminary Estimates

2.1 Preliminaries

We now present various notations and definitions that will be used frequently. The usual norms in the spaces $\mathbb{L}^{p}(\Omega)$ , $\mathbb{L}^{\infty}(\Omega)$ and $\mathbb{C}\left(\overline{\Omega}\right)$ are respectively denoted by

\left\|u\right\|_{p}^{p}\text{=}\frac{1}{\left|\Omega\right|}\int\limits_{\Omega}\left|u(x)\right|^{p}dx,

\left\|u\right\|_{\infty}\text{=}\underset{x\in\Omega}{max}\left|u(x)\right|.

We define the following phase spaces

H=L^{2}(\Omega)\times L^{2}(\Omega)\times L^{2}(\Omega),

E=H^{1}(\Omega)\times H^{1}(\Omega)\times H^{1}(\Omega),

X=H^{2}(\Omega)\times H^{2}(\Omega)\times H^{2}(\Omega).

We now recall the following lemma.

Lemma 2.1 (Uniform Gronwall Lemma)

Let $\beta,\zeta$ , and $h$ be nonnegative functions in $L^{1}_{loc}[0,\infty;\mathbb{R})$ . Assume that $\beta$ is absolutely continuous on $(0,\infty)$ and the following differential inequality is satisfied:

\frac{d\beta}{dt}\leq\zeta\beta+h,\ \mbox{for}\ t>0.

(5)

If there exists a finite time $t_{1}>0$ and some $q>0$ such that

\int^{t+q}_{t}\zeta(\tau)d\tau\leq A,\ \int^{t+q}_{t}\beta(\tau)d\tau\leq B,\ \mbox{and}\ \int^{t+q}_{t}h(\tau)d\tau\leq C,

(6)

for any $t>t_{1}$ , where $A,B$ , and $C$ are some positive constants, then

\beta(t)\leq\left(\frac{B}{q}+C\right)e^{A},\ \mbox{for \ any}\ t>t_{1}+q.

(7)

2.2 Uniform $L^{2}(\Omega)$ and $H^{1}(\Omega)$ estimates

In all estimates made hence forth the constants $C,C_{1},C_{2},C_{3},C_{\infty}$ are generic constants, that can change in value from line to line, and sometimes within the same line if so required. Using positivity of the solution, and comparison arguments, one is easily able to prove global existence of classical solutions to system (2)-(4). We state the following result.

Proposition 2.2

All positive solutions of the system (2)-(4), with initial data in $\mathbb{L}^{\infty}(\Omega)$ are classical and global.

See appendix 7.2 for the proof. Although we have global existence of classical via proposition 2.2, we can actually derive uniform $L^{\infty}(\Omega)$ bounds, for initial data only in $L^{2}(\Omega)$ . This will facilitate estimates required to show existence of bounded absorbing sets and the compactness of trajectories, which will lead to the existence of a global attractor in the product space $L^{2}(\Omega)$ . Since the construction of global attractors require a Hilbert space setting, where the appropriate function spaces are the Hilbert spaces $L^{2}(\Omega)$ and $H^{s}(\Omega)$ , $s=1,2$ , we want to proceed by showing there is a weak solution in these function spaces. We then make estimates on this class of solutions, to proceed with the requisite analysis for the global attractors. We state the following theorem:

Theorem 2.3

Consider the three species food chain model described via (2)-(4). For any initial data $(u_{0},v_{0},r_{0})$ in $L^{2}(\Omega)$ , and spatial dimension $n=1,2,3$ , there exists a global weak solution $(u,v,r)$ to the system, which becomes a strong solution and then a classical solution.

The proofs of proposition 2.2 and theorem 2.3 follow via remark 2.7 and the estimates referred to therein.

We state the following result next,

Lemma 2.4

Consider $(u,v,r)$ that are solutions to the diffusive three species food chain model described via (2)-(4), for any $(u_{0},v_{0},r_{0})\in L^{2}(\Omega)$ , and spatial dimension $n=1,2,3$ , there exists a time $t^{*}(||u_{0}||_{2},||v_{0}||_{2})$ , and a constant $C$ independent of time, initial data and the allee threshold $m$ , and dependent only on the other parameters in (2)-(4), such that for any $t>t^{*}$ the following uniform estimates hold:

||u||^{2}_{2}\leq C,||v||^{2}_{2}\leq C,\int^{t+1}_{t}||\nabla u||^{2}_{2}ds\leq C,\int^{t+1}_{t}||\nabla v||^{2}_{2}ds\leq C,||\nabla u||^{2}_{2}\leq C,||\nabla v||^{2}_{2}\leq C.

(8)

This follows easily via the estimates and methods of [20, 18, 19]. Note what differs here from [20, 18, 19] is the equation for $r$ , particularly due to a Allee effect now being modeled.

We now proceed with the estimate on $r$ . We multiply (4) by $r$ and integrate by parts to obtain

\frac{1}{2}\frac{d}{dt}||r||^{2}_{2}+d_{3}||\nabla r||^{2}_{2}+\frac{w_{3}}{||v||_{\infty}+D_{3}}||r||^{4}_{4}\leq\left(c+\frac{mw_{3}}{v+D_{3}}\right)||r||^{3}_{3}.

(9)

Thus we obtain

\frac{1}{2}\frac{d}{dt}||r||^{2}_{2}+d_{3}||\nabla r||^{2}_{2}+cm||r||^{2}_{2}+\frac{w_{3}}{||v||_{\infty}+D_{3}}||r||^{4}_{4}\leq\left(c+\frac{mw_{3}}{v+D_{3}}\right)||r||^{3}_{3}.

(10)

We then use Hölder’s inequality followed by Young’s inequality to obtain

	$\displaystyle\frac{1}{2}\frac{d}{dt}\|\|r\|\|^{2}_{2}+d_{3}\|\|\nabla r\|\|^{2}_{2}+cm\|\|r\|\|^{2}_{2}+\frac{w_{3}}{\|\|v\|\|_{\infty}+D_{3}}\|\|r\|\|^{4}_{4},$
	$\displaystyle\leq\frac{w_{3}}{\|\|v\|\|_{\infty}+D_{3}}\|\|r\|\|^{4}_{4}+C_{1}(\|\|v\|\|_{\infty}+D_{3}+m^{4})\|\Omega\|.$

Thus we obtain

\frac{1}{2}\frac{d}{dt}||r||^{2}_{2}+cm||r||^{2}_{2}\leq C_{1}(||v||_{\infty}+D_{3}+m^{4})|\Omega|,

(12)

which implies

||r||^{2}_{2}\leq e^{-cmt}||r_{0}||^{2}_{2}+\frac{C_{1}(||v||_{\infty}+D_{3}+m^{4})|\Omega|}{cm}.

(13)

We now use the estimate for $||v||_{\infty}$ via (27) (which does not depend on $r$ see remark 2.8) to obtain the following estimates for $r$ ,

||r||^{2}_{2}\leq 1+\frac{C_{1}(C+D_{3}+m^{4})|\Omega|}{cm},\ \mbox{for}\ t>\max\left(t^{*}+1,\frac{\ln(||r_{0}||^{2}_{2})}{cm}\right).

(14)

Integrating (7.4) in the time interval $[t,t+1]$ we obtain

\int^{t+1}_{t}||\nabla r||^{2}_{2}ds\leq 1+\frac{C_{1}(C+D_{3}+m^{4})|\Omega|}{cm}+2|\Omega|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}},

(15)

for $t>\max\left(t^{*}+1,\frac{\ln(||r_{0}||^{2}_{2})}{cm}\right)$ .

Here, $C,C_{1}$ do not depend on $m$ .

Remark 2.5

What we notice is that the estimate via (14), (15) depend singularly on $m$ . That is they yield no information, if we try to pass to the limit as $m\rightarrow 0$

We now estimate the gradient of $r$ . See the appendix 7.3 for details. What the standard analysis in appendix 7.3 shows, is that we have an estimate of the form

\mathop{\limsup}_{t\rightarrow\infty}||\nabla r||^{2}_{2}\leq C,

(16)

Remark 2.6

The constant $C$ in (16) is independent of time and initial conditions, but depend singularly on $m$ . In our estimates we apply the uniform Gronwall lemma which requires us to use the estimate via (15). This also depends singularly on $m$ , thus the estimate (16) is uniform with respect to time and initial data but depend singularly on $m$ .

Remark 2.7

The uniform $H^{1}(\Omega)$ estimates via lemma 2.4 and (95) give us uniform $L^{6}(\Omega)$ bounds in $\mathbb{R}^{3}$ . We can now prove the reaction terms are in $L^{p}(\Omega)$ for $p>\frac{3}{2}$ . It suffices to show that

||u-u^{2}-w_{1}\frac{uv}{u+v}||_{2}\leq C||u||^{2}_{4}\leq C,

(17)

||-a_{2}v+w_{2}\frac{uv}{u+v}-w_{3}\left(\frac{vr}{v+r}\right)||_{2}\leq C||v||_{2}\leq C,

(18)

||r\big{(}{{r-m}}\big{)}\bigg{(}c-\frac{w_{4}r}{v+D_{3}}\bigg{)}||_{2}\leq C||r||^{3}_{6}\leq C.

(19)

But the above follows via the uniform $L^{6}(\Omega)$ bounds on $u,v,r$ . This proves proposition 2.2. Furthermore the estimate via (16) and lemma 2.4 allow us to derive the appropriate $L^{2}(\Omega)$ and $H^{1}(\Omega)$ bounds on a Galerkin truncation of (2)-(4), extract appropriate subsequebces and pass to the limit as is standard [27, 22], to prove theorem 2.3.

2.3 Uniform in the parameter $m$ $L^{2}(\Omega)$ and $H^{1}(\Omega)$ estimates for $r$

Note the estimates via (14), (16) are singular in $m$ . This causes extensive difficulties if we try to pass to the limit as $m\rightarrow 0$ , hence making it difficult to prove upper semicontinuity. Our next goal is to derive uniform in $m$ estimate on $r$ . These are derived in detail in the appendix 7.4. However, the key estimate derived therein is

\displaystyle||\nabla r||^{2}_{2}\leq C_{K}e^{KC_{K}},\ \mbox{for}\ t>t_{1}=t_{0}+1.

(20)

Here $C_{K},K$ are independent of $m$ . Thus the $H^{1}(\Omega)$ estimate for $r$ can be made independent of $m$ .

2.4 Uniform $H^{2}(\Omega)$ estimates

We will now estimate the $H^{2}(\Omega)$ norms of the solution. We rewrite (2) as

u_{t}-d_{1}\Delta u=u-u^{2}-w_{1}\left(\frac{uv}{u+v}\right).

We square both sides of the equation and integrate by parts over $\Omega$ to obtain

			$\displaystyle\|\|u_{t}\|\|^{2}_{2}+d_{1}\|\|\Delta u\|\|^{2}_{2}+\frac{d}{dt}\|\|\nabla u\|\|^{2}_{2},$
		$\displaystyle=$	$\displaystyle\left\\|\left(u-u^{2}-w_{1}\left(\frac{uv}{u+v}\right)\right)\right\\|_{2}^{2},$
		$\displaystyle\leq$	$\displaystyle C\left(\|\|u\|\|^{2}_{2}+(w_{1})^{2}\|\|u\|\|^{2}_{2}+\|\|u\|\|^{4}_{4}\right)\leq C.$

This result follows by the embedding of $H^{1}(\Omega)\hookrightarrow L^{4}(\Omega)\hookrightarrow L^{2}(\Omega)$ . Therefore, we obtain

||u_{t}||^{2}_{2}+d_{1}||\Delta u||^{2}_{2}+\frac{d}{dt}||\nabla u||^{2}_{2}\leq C\left(||u||^{2}_{2}+(w_{1})^{2}||u||^{2}_{2}+||u||^{4}_{4}\right)\leq C.

We now make uniform in time estimates of the higher order terms. Integrating the estimates of (2.4) in the time interval $[T,T+1]$ , for $T>t^{*}$

\int^{T+1}_{T}||u_{t}||^{2}_{2}dt\leq C_{1},\ \int^{T+1}_{T}||\Delta u||^{2}_{2}dt\leq C_{2}.

Similarly we adopt the same procedure for $v$ to obtain

\int^{T+1}_{T}||v_{t}||^{2}_{2}dt\leq C_{1},\ \int^{T+1}_{T}||\Delta v||^{2}_{2}dt\leq C_{2}.

Next, consider the gradient of (2). Following the same technique as in deriving (2.4) we obtain for the left hand side

	$\displaystyle\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2}+\int_{\partial\Omega}\Delta u\nabla u_{t}\cdot\textbf{n}dS,$
	$\displaystyle=\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2}+\int_{\partial\Omega}\Delta u\frac{\partial}{\partial t}(\nabla u\cdot\textbf{n})dS,$
	$\displaystyle=\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2}.$		(21)

This follows via the boundary condition. Thus we have

	$\displaystyle\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2},$
	$\displaystyle=\left(\nabla(u-u^{2}-w_{1}\left(\frac{uv}{u+v}\right))\right)^{2},$
	$\displaystyle\leq C(\|\|u\|\|^{4}_{4}+\|\|v\|\|^{4}_{4}+\|\|\Delta u\|\|^{4}_{2}+\|\|\nabla v\|\|^{2}_{2}+\|\|\nabla u\|\|^{2}_{2}).$		(22)

This follows via Young’s inequality with epsilon, as well as the embedding of $H^{2}(\Omega)\hookrightarrow W^{1,4}(\Omega)\hookrightarrow H^{1}(\Omega)$ . This implies that

\frac{d}{dt}||\Delta u||^{2}_{2}\leq C||\Delta u||^{4}_{2}+C(||u||^{4}_{4}+||v||^{4}_{4}+||\nabla v||^{2}_{2}+||\nabla u||^{2}_{2}).

We now use the uniform Gronwall lemma with

\beta(t)=||\Delta u||^{2}_{2},\ \zeta(t)=||\Delta u||^{2}_{2},\ h(t)=C(||u||^{4}_{4}+||v||^{4}_{4}+||\nabla v||^{2}_{2}+||\nabla u||^{2}_{2}),~~q=1,

to obtain

||\Delta u||_{2}\leq C,\ ||\Delta v||_{2}\leq C,\ \mbox{for}\ t>t^{*}+1.

The estimate for $v$ is derived similarly.

Remark 2.8

It is critical to notice that the estimate for $||\Delta v||_{2}$ does not involve $r$ . This is because in order to derive it, we adopt same procedure as in (2.4), instead for the $v$ equation, and square both sides to proceed. Therein notice that the term involving $r$ , $\frac{vr}{v+r}$ can be estimated as $||\frac{vr}{v+r}||_{2}\leq C||\frac{r}{v+r}||_{\infty}||v||_{2}\leq C_{1}||v||_{2}$ . So the right hand side has no dependence on $r$ .

In order to estimate the $H^{2}(\Omega)$ norm of $r$ we proceed similarly to obtain

$\displaystyle\|\|r_{t}\|\|^{2}_{2}+d_{3}\|\|\Delta r\|\|^{2}_{2}+\frac{d}{dt}\|\|\nabla r\|\|^{2}_{2}$	$\displaystyle=$	$\displaystyle\left\\|\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}-mcr\right)\right\\|_{2}^{2},$
	$\displaystyle\leq$	$\displaystyle C\left(\|\|r\|\|^{2}_{2}+\|\|r\|\|^{4}_{4}\right)+\|\|r\|\|^{6}_{6},$
	$\displaystyle\leq$	$\displaystyle C\left(\|\|\nabla r\|\|^{2}_{2}\right)^{3}.$

We now make uniform in time estimates of the higher order terms. Integrating the estimates of (2.4) in the time interval $[T,T+1]$ , for $T>t_{1}$ , and using the estimates via (7.4) yields

\int^{T+1}_{T}||r_{t}||^{2}_{2}dt\leq C_{1},\ \int^{T+1}_{T}||\Delta r||^{2}_{2}dt\leq C_{2}.

Next, consider the gradient of (4). Following the same technique as in deriving (2.4) we obtain

	$\displaystyle\|\|\nabla r_{t}\|\|^{2}_{2}+d_{3}\|\|\nabla(\Delta r)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta r\|\|^{2}_{2}$		$\displaystyle=\left(\nabla\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}-mcr\right)\right)^{2},$
			$\displaystyle\leq C\left(\|\|\nabla r\|\|^{2}_{2}\right)^{2}.$

This yields

\frac{d}{dt}||\Delta r||^{2}_{2}\leq C\left(||\nabla r||^{2}_{2}\right)^{2}\left(||\Delta r||^{2}_{2}\right)+C_{2}\left(||r||^{2}_{2}\right)^{2}.

(25)

We can now use the estimates via (99), (7.4) in conjunction with the uniform Gronwall lemma to yield a uniform bound on $||\Delta r||_{2}$ . Thus, via elliptic regularity,

||u||_{H^{2}(\Omega)}\leq C,\ ||v||_{H^{2}(\Omega)}\leq C,\ ||r||_{H^{2}(\Omega)}\ \leq C\ \mbox{for}\ t>\max(t_{1}+1,t^{*}+1).

(26)

Since $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ , the following estimate is valid in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$

||u||_{\infty}\leq C,\ ||v||_{\infty}\leq C,\ ||r||_{\infty}\leq C\ \mbox{for}\ t>\max(t_{1}+1,t^{*}+1).

(27)

3 Long Time Dynamics

3.1 Existence of global attractor

In this section we prove the existence of a $(H,H)$ global attractor for system (2)-(4), which is subsequently demonstrated to be a $(H,E)$ attractor. We now state the following theorem,

Theorem 3.1

Consider the reaction diffusion equation described via (2)-(4) where $\Omega$ is of spatial dimension $n=1,2,3$ . There exists a $(H,H)$ global attractor $\mathcal{A}$ for the system. This is compact and invariant in $H$ , and it attracts all bounded subsets of $H$ in the $H$ metric.

proof We have shown that the system is well posed via theorems 2.3. Thus there exists a well defined semigroup $\left\{S(t)\right\}_{t\geq 0}:H\rightarrow H$ . The estimates derived in Lemma 2.4 demonstrate the existence of bounded absorbing sets in $H$ and $E$ . Thus given a sequence $\left\{u_{0,n}\right\}^{\infty}_{n=1}$ , that is bounded in $L^{2}(\Omega)$ , we know that for $t>t^{*}$ ,

S(t)(u_{0,n})\subset B\subset H^{1}(\Omega).

(28)

Here $B$ is the bounded absorbing set in $H^{1}(\Omega)$ from lemma 2.4. Now for n large enough $t_{n}>t^{*}$ , thus for such $t_{n}$ we have

S(t_{n})(u_{0,n})\subset B\subset H^{1}(\Omega).

(29)

This implies that we have the following uniform bound,

||S(t_{n})(u_{0,n})||_{H^{1}(\Omega)}\leq C,

(30)

This implies via standard functional analysis theory, see [27], the existence of a subsequence still labeled $S(t_{n})(u_{0,n})$ such that

S(t_{n})(u_{0,n})\rightharpoonup u\ \mbox{in}\ H^{1}(\Omega),

(31)

Which implies via the compact Sobolev embedding of

E\hookrightarrow H,

(32)

that

S(t_{n})(u_{0,n})\rightarrow u\ \mbox{in}\ L^{2}(\Omega).

(33)

This yields the asymptotic compactness of the semigroup $\left\{S(t)\right\}_{t\geq 0}$ in $H$ . The convergences for the $v,r$ components follow similarly. The theorem is now proved. $\square$

Next we can state the following theorem

Theorem 3.2

Consider the reaction diffusion equation described via (2)-(4) where $\Omega$ is of spatial dimension $n=1,2,3$ . There exists a $(H,E)$ global attractor $\mathcal{A}$ for the system. This is compact and invariant in $H$ , and it attracts all bounded subsets of $H$ in the $E$ metric.

proof The proof essentially follows verbatim theorem 3.1. The existence of a bounded absorbing set in $E$ follows from lemma 2.4. In order to prove the asymptotic compactness in $E$ we use the uniform $H^{2}(\Omega)$ estimates in (26), and the compact Sobolev embedding of $H^{2}(\Omega)\hookrightarrow H^{1}(\Omega)$ . $\square$

3.2 Upper-semicontinuity of Global attractor with respect to the Allee parameter $m$

Here we state and prove the main semicontinuity result. This result follows via a series of estimates, and theorems, all of which we derive for the reader in appendix 7.5

Theorem 3.3

Consider the reaction diffusion equation described via (2)-(4) where $\Omega$ is of spatial dimension $n=1,2,3$ . Given a set of positive parameters excluding $m$ , the family of global attractors is upper semicontinuous with respect to the Allee threshold $m\geq 0$ as it converges to zero. That is

dist_{E}(\mathcal{A}_{m},\mathcal{A}_{0})\rightarrow 0,\ \mbox{as}\ m\rightarrow 0^{+}.

(34)

Proof 3.4

We know that for $\epsilon>0$ , the global attractor $\mathcal{A}_{0}$ for the reaction diffusion equation described via (2)-(4) attracts $\mathcal{U}$ . This follows via theorems 7.6, 7.7 and 7.9. Thus there exists a finite time $t_{\epsilon}>0$ , such that

S_{0}(t_{\epsilon})\mathcal{U}\subset\mathcal{N}_{E}(\mathcal{A}_{0},\frac{\epsilon}{2}).

(35)

Here, $\mathcal{N}_{E}(\mathcal{A}_{0},\frac{\epsilon}{2})$ is the $\frac{\epsilon}{2}$ ball of $\mathcal{A}_{0}$ in the space $E$ . Now we also know by the uniform convergence theorem, theorem 7.10 that there exists a $m_{\epsilon}\in(0,1]$ such that

\sup_{g_{0}\in\mathcal{U}}||S_{m}(t_{\epsilon})g_{0}-S_{0}(t_{\epsilon})g_{0}||_{E}\leq\frac{\epsilon}{2},\ \mbox{for}\ \mbox{any}\ m\in m_{\epsilon}.

(36)

Since we know $\mathcal{A}_{m}$ is invariant we have

\mathcal{A}_{m}=S_{m}(t_{\epsilon})\mathcal{A}_{m}\subset S_{m}(t_{\epsilon})\mathcal{U}\subset\mathcal{N}_{E}(S_{0}(t_{\epsilon})\mathcal{U},\frac{\epsilon}{2})\subset\mathcal{N}_{E}(\mathcal{A}_{0},\epsilon).

(37)

Thus one obtains the upper semicontinuity,

dist_{E}(\mathcal{A}_{m},\mathcal{A}_{0})\rightarrow 0,\ \mbox{as}\ m\rightarrow 0^{+}.

(38)

This proves the theorem.

3.3 Computing the decay rate with respect to the Allee threshold parameter $m$

In the previous subsection 3.2, the role of the Allee threshold parameter $m$ on the long time dynamics of system (2)-(4) is explored. Theorem 3.3 shows that the global attractors $\mathcal{A}_{m}$ converge to $\mathcal{A}_{0}$ in the $H^{1}(\Omega)$ norm, as $m$ decreases to zero. It is also seen how changing $m$ effects the Turing dynamics in subsection 5. Theorem 3.3 is a purely theoretical result. It gives no information as to the exact decay rate in terms of the parameter $m$ . In general, there are no results in the literature to estimate this decay rate. We now aim to examine the role of $m$ with respect to convergence, as $m$ decreases to zero as discussed earlier. We aim to explicitly find the decay rate of $\mathcal{A}_{m}$ to a target attractor $\mathcal{A}_{0}$ (that is when $m=0$ ). We choose Turing patterns as the state of interest. In particular the numerical scheme discussed in 5 was employed, and the Turing patterns at an extended time are examined. The goal is to gather numerical evidence as to the dependence on $m$ with respect to the convergence rate,

\displaystyle\|U_{m}-U_{0}\|_{H^{1}(\Omega)}

\displaystyle<

\displaystyle Cm^{\gamma},

(39)

where $U_{0}=(u_{,}v_{0},r_{0})$ is the solution for $m=0$ , $U_{m}=(u_{m},v_{m},r_{m})$ is the solution for any given value of $m$ , and both $C$ and $\gamma$ are constants.

The numerical approximation for the case $m=0$ was established at the scaled time, $t=10,000$ . The approximations, $u_{m}(x,t)$ , $v_{m}(x,t)$ , and $r_{m}(x,t)$ , were then established for a range of values of $m$ with 121 equally spaced values from 0 to 0.0035. The approximation for each value of $m$ was established up to $t=10,000$ , and it was determined that for each value of $m$ in this range a fixed spatial pattern was found. The $H^{1}(\Omega)$ errors for the numerical approximation were then calculated with respect to the $m=0$ case,

	$\displaystyle\\|u_{m}(\cdot,10,000)-u_{0}(\cdot,10,000)\\|_{H^{1}(\Omega)},$
	$\displaystyle\\|v_{m}(\cdot,10,000)-v_{0}(\cdot,10,000)\\|_{H^{1}(\Omega)},$
	$\displaystyle\\|r_{m}(\cdot,10,000)-r_{0}(\cdot,10,000)\\|_{H^{1}(\Omega)}.$

Refer to caption — Figure 1: The $H^{1}(\Omega)$ errors for the differences $\|u_{m}(\cdot,10,000)-u_{0}(\cdot,10,000)\|$ , $\|v_{m}(\cdot,10,000)-v_{0}(\cdot,10,000)\|$ , and $\|r_{m}(\cdot,10,000)-r_{0}(\cdot,10,000)\|$ are shown for a range of values of $m$ . The solid line is the linear least squares regression line.

The errors are shown in Figure 1 and Figure 2. The first figure, Figure 1, shows the raw errors, and the second figure, Figure 2, shows the errors as a log-log plot. A linear least squares fit for the data was also calculated, and the best fit straight line is shown on both plots. In both cases the close linear fit for larger values of $m$ contribute to strong, positive correlations in the data. However, for smaller values $m$ the log-log data moves away from a straight line. For the smaller values of $m$ the errors in the numerical approximation are of the same order as the error in the $H^{1}(\Omega)$ norm, and the data no longer follows the same trend.

Examining the errors for $u_{m}$ , the correlation for the raw data is 0.99995, and the data reveals a strong, positive linear relationship. In this scenario, it is assumed that the value of $\gamma$ in equation (39) is one. The approximation for the slope of the raw data gives an estimate of 20.6 with a 95% confidence interval between 20.6215 and 20.6206. The 95% confidence interval for the intercept is between -5.1E-5 and -9.0E-4. Due to the large number of samples, the error in the confidence interval is very small, but the estimate for the parameter is close to zero.

Examining the errors for the log of the errors correlation is 0.993 for the log-log data. The approximation for the slope of the log-log data gives an estimate of 1.003 with a 95% confidence interval between 0.72 and 1.28. The 95% confidence interval for the intercept is between 2.7 and 3.3. Note that the value of one is within the confidence interval making it difficult to claim that the relationship is not linear. The estimate for the slope, in fact, is very close to one indicating that the value of $\gamma$ in equation (39) is close to one.

4 Overexploitation Phenomenon

Overexploitation refers to the phenomenon where excessive harvesting of a species can result in its extinction. This finds applications in conservation, biodiversity, cascade control and so on. Although there is a fair amount of literature on this in two species models [32], it is much less studied in three species models. To begin, we state the following theorem

Theorem 4.1 (Middle predator and Allee effect mediated overexploitation)

Consider $(u,v,r)$ that are solutions to the diffusive three species food chain model described via (2)-(4). If $w_{1}>a_{2}+1+w_{3}$ , then for any given initial prey density $u_{0}$ there exists a threshold $M_{1}=\frac{1}{\alpha}||u_{0}||_{p}$ , s.t if $M_{1}<||v_{0}||_{p}$ and $||r_{0}||_{\infty}<\min(m,\frac{cD_{3}}{w_{4}})$ , then $(u,v,r)\rightarrow(0,0,0)$ uniformly for $x\in\bar{\Omega}$ as $t\rightarrow\infty$ .

proof: We begin by looking at the equation for $u$ , and following the ideas in [12]

\frac{\partial u}{\partial t}=d_{1}\Delta u+u-u^{2}-w_{1}\frac{uv}{u+v}=d_{1}\Delta u+u-u^{2}-w_{1}\frac{u}{\frac{u}{v}+1}.

(40)

We will claim that for all $t$ , $||\frac{u}{v}||_{\infty}<\alpha$ , and $\lim_{t\rightarrow\infty}||u||_{p}\rightarrow 0$ , $\forall p$ , if $\frac{||u_{0}||_{p}}{||v_{0}||_{p}}<\alpha$ . If not there exists a first time $t_{1}$ s.t., $\frac{||u||_{p}}{||v||_{p}}<\alpha$ on $t\in[0,t_{1}]$ and $\frac{||u(t_{1})||_{p}}{||v(t_{1})||_{p}}=\alpha$ .

Due to $w_{1}>a_{2}+1+w_{3}$ , there exists an $\alpha$ s.t $\frac{w_{1}}{1+\alpha}=1+a_{2}+w_{3}$ , thus

$\displaystyle\frac{\partial u}{\partial t}=d_{1}\Delta u+u-u^{2}-w_{1}\frac{u}{\frac{u}{v}+1}$	$\displaystyle\leq d_{1}\Delta u+u-w_{1}\frac{u}{\alpha+1},$	(41)
	$\displaystyle=d_{1}\Delta u+u-u-a_{2}u-w_{3}u,$
	$\displaystyle=d_{1}\Delta u-(a_{2}+w_{3})u.$

Thus $u$ satisfies

||u||_{p}\leq||u_{0}||_{p}e^{-(a_{2}+w_{3})t}.

(42)

This easily follows via multiplying (2) by $|u|^{p-1}$ , and integrating by parts. We also notice that the $\tilde{v}$ solving the following equation, with $\tilde{v}_{0}=\sup v_{0}(x)$ , is a super solution to (3)

\frac{\partial\tilde{v}}{\partial t}=-(a_{2}+w_{3})\tilde{v}.

(43)

Multiplying through by $|v|^{p-1}$ , and integrating by parts yields

||v||_{p}\geq||v_{0}||_{p}e^{-(a_{2}+w_{3})t}.

(44)

This implies

\frac{||u||_{p}}{||v||_{p}}\leq\frac{||u_{0}||_{p}}{||v_{0}||_{p}}<\alpha.

(45)

This implies a contradiction. Thus $\lim_{t\rightarrow\infty}||u||_{p}\rightarrow 0$ , which implies uniform convergence of a subsequence, say $u_{n_{j}}$ to 0 (where $u_{n}$ is a Galerkin truncation of $u$ and standard theory is adopted [22]), and by uniqueness of solutions, $u\rightarrow 0$ uniformly, as $t\rightarrow\infty$ . This reduces the $v$ equation to

\frac{\partial v}{\partial t}=d_{2}\Delta v-a_{2}v-w_{3}\left(\frac{vr}{v+r}\right)\leq d_{2}\Delta v-a_{2}v,

(46)

and it is easy to see that $v\rightarrow 0$ uniformly, as $t\rightarrow\infty$ .

Now the $r$ equation reduces to

\frac{\partial r}{\partial t}=d_{3}\Delta r+r\big{(}{{r-m}}\big{)}\bigg{(}c-\frac{w_{4}r}{D_{3}}\bigg{)},

(47)

and if $||r_{0}||_{\infty}<\min(m,\frac{cD_{3}}{w_{4}})$ standard analysis yields

$r\rightarrow 0$ uniformly.

This proves the theorem.

$\square$

Remark 4.2

Note, in the three species food chain model (2)-(4), overexploitation has to be middle predator mediated, and also depends on the allele threshold. Clearly if $||r_{0}||_{p}>\max(m,\frac{cD_{3}}{w_{4}})$ , or $\frac{cD_{3}}{w_{4}}<||r_{0}||_{p}<m$ , or $\frac{cD_{3}}{w_{4}}>||r_{0}||_{p}>m$ , $r$ cannot go extinct.

Also, if the attack rate of the top predator $r$ is large enough, it could cause $v$ to go extinct, before $u$ does. Once $v$ goes extinct $u\rightarrow 1$ , again prohibiting overexploitation.

Theorem 4.3 (Top predator mediated overexploitation avoidance)

Consider $(u,v,r)$ that are solutions to the diffusive three species food chain model described via (2)-(4). If $w_{3}>w_{2}+1+w_{1}$ , then for any given initial prey density $u_{0}$ there exists a threshold $M_{2}=\frac{1}{\alpha_{1}}||r_{0}||_{p}$ , s.t if $\max\left(M_{2},m,\frac{cD_{3}}{w_{4}}\right)<||r_{0}||_{p}$ then $(u,v,r)\rightarrow(1,0,r^{*})$ uniformly as $t\rightarrow\infty$ .

proof: We first derive a lower bound on $u$ . Trivially we have $||u||_{\infty}\leq 1$ . Now

\frac{\partial u}{\partial t}\geq-u^{2}-w_{1}\frac{vu}{v+u}\geq-u-w_{1}u.

(48)

This yields

||u||_{p}\geq||u_{0}||_{p}e^{-(1+w_{1})t}.

(49)

Next we will claim that for all $t$ , $||\frac{v}{r}||_{\infty}<\alpha_{1}$ , and $\lim_{t\rightarrow\infty}||v||_{p}\rightarrow 0$ , $\forall p$ , if $\frac{||v_{0}||_{p}}{||r_{0}||_{p}}<\alpha_{1}$ . If not there exists a first time $t_{2}$ s.t., $\frac{||v||_{p}}{||r||_{p}}<\alpha_{1}$ on $t\in[0,t_{2}]$ and $\frac{||u(t_{2})||_{p}}{||v(t_{2})||_{p}}=\alpha_{1}$ .

Due to the choice of $w_{3}>w_{2}+1+w_{1}$ , there exists an $\alpha_{1}$ s.t $\frac{w_{3}}{1+\alpha_{1}}=1+w_{2}+w_{1}$ , thus

	$\displaystyle\frac{\partial v}{\partial t}$		$\displaystyle=d_{2}\Delta v-a_{2}v+w_{2}\frac{vu}{v+u}-w_{3}\frac{v}{\frac{v}{r}+1},$
			$\displaystyle\leq d_{2}\Delta v-a_{2}v+w_{2}v-\frac{w_{3}}{1+\alpha_{1}}v,$
			$\displaystyle=d_{2}\Delta v-a_{2}v+w_{2}v-(1+w_{2}+w_{1})v.$

Thus $v$ satisfies

||v||_{p}\leq||v_{0}||_{p}e^{-(a_{2}+1+w_{3})t}.

(51)

By the restriction on $r_{0}$ , we see that trivially $r\rightarrow m$ or $r\rightarrow\frac{cD_{3}}{w_{3}}$ . This implies

\frac{||v||_{p}}{||r||_{p}}\leq\frac{||v_{0}||_{p}}{||r_{0}||_{p}}<\alpha_{1}.

(52)

This implies a contradiction. Thus $\lim_{t\rightarrow\infty}||v||_{p}\rightarrow 0$ , which implies uniform convergence to 0 of a subsequence, say $v_{n_{j}}$ (where $v_{n}$ is a Galerkin truncation of $v$ and standard theory is adopted [22]), and by uniqueness $v\rightarrow 0$ uniformly as $t\rightarrow\infty$ . Looking at the decay rate for $v$ , we see it converges to 0 faster than $u$ , so $u>0$ , when $v\rightarrow 0$ uniformly, and so then the $u$ equation reduces to

\frac{\partial u}{\partial t}=d_{1}\Delta u+u-u^{2},

(53)

and we know via standard theory that now $u\rightarrow 1$ uniformly.

This proves the theorem. $\square$

Remark 4.4

This also happens if $\frac{cD_{3}}{w_{4}}<||r_{0}||_{\infty}<m$ , or $\frac{cD_{3}}{w_{4}}>||r_{0}||_{\infty}>m$ .

The following lemma describes persistence in the top predator.

Lemma 4.5 (Persistence in top predator)

Consider $(u,v,r)$ that are solutions to the diffusive three species food chain model described via (2)-(4). If $\min(m,\frac{cD_{3}}{w_{4}})<||r_{0}||_{\infty}$ , then $r\rightarrow\max(m,\frac{cD_{3}}{w_{4}})$ uniformly as $t\rightarrow\infty$ , hence persists.

proof: Note $c-\frac{w_{4}r}{v+D_{3}}>c-\frac{w_{4}r}{D_{3}}$ , due to positivity of $r,v$ . Thus we can compare the solution of (4) to the solution of the ODE

\frac{d\tilde{r}}{dt}=\tilde{r}(\tilde{r}-m)\bigg{(}c-\frac{w_{4}\tilde{r}}{v+D_{3}}\bigg{)},

(54)

with $\tilde{r}_{0}=\min(r_{0}(x))>\min(m,\frac{cD_{3}}{w_{4}})$ . Clearly $\tilde{r}\rightarrow\max(m,\frac{cD_{3}}{w_{4}})$ as $t\rightarrow\infty$ , and since $r$ solving (4) is a supersolution to $\tilde{r}$ solving (54), we have $\liminf_{t\rightarrow\infty}\min_{\bar{\Omega}}(r(x,t)\geq\max(m,\frac{cD_{3}}{w_{4}})$ .

$\square$

Remark 4.6

In the event that there is no Allee effect, or $m=0$ , $r\rightarrow\frac{cD_{3}}{w_{4}}$ uniformly as $t\rightarrow\infty$ , hence always persists, independent of all other parameters or initial data $(u_{0},v_{0},r_{0})$ . Thus in a three species system where the top predator is a generalist, and can switch its favorite food source, true overexploitation can only take place if an Allee effect is in place, and not otherwise.

Remark 4.7

All of the over exploitation theorems were verified numerically, for various ranges of parameter values.

5 Turing Instability

In this section we establish the condition needed to ensure a positive interior steady state for model system (2-4). We focus on deriving the condition necessary and sufficient for Turing instability to occur as a result of the introduction of diffusion. This phenomena is referred to as diffusion driven instability, which was first introduced by Alan Turing[29]. We first establish the steady state for model system (55-57)

5.1 Steady States

Before we can begin the analysis of Turing instability, we have to establish the conditions needed to ensure a positive interior steady state for the ODE version of model system (2)-(4), which is given as

	$\displaystyle\frac{du}{dt}=u-u^{2}-w_{1}\frac{uv}{u+v},$		(55)
	$\displaystyle\frac{dv}{dt}=-a_{2}v+w_{2}\frac{uv}{u+v}-w_{3}\left(\frac{vr}{v+r}\right),$		(56)
	$\displaystyle\frac{dr}{dt}=r\big{(}{{r-m}}\big{)}\bigg{(}c-\frac{w_{4}r}{v+D_{3}}\bigg{)}.$		(57)

Also all parameters associated with model system (55-57) are assumed to be positive constants and have the usual biological meaning.

The model system (55)-(57) has the following non-negative steady states

1.

Total Extinction of the three species: $E_{0}(0,0,0)$ .
2.

Predator free steady state: $E_{1}(1,0,0)$ .
3.

No primary food source for predator steady state: $E_{2}(0,0,m)$ and $E_{3}(0,0,\frac{cD_{3}}{w_{4}})$ .
4.

Top Predator free steady state: $E_{4}(\tilde{u},\tilde{v},0)$ where $\tilde{u}=\frac{w_{2}-w_{1}w_{2}+a_{2}w_{1}}{w_{2}},\tilde{v}=\frac{(w_{2}-a_{2})\tilde{u}}{a_{2}}$ if $w_{2}>w_{1}(w_{2}-a_{2}),w_{2}>a_{2}$ .
5.

Middle Predator free steady state: $E_{5}(1,0,\frac{cD_{3}}{w_{4}})$ and $E_{6}(1,0,m)$ .
6.

Coexistence of Species: $E_{7}(\tilde{\tilde{u}},\tilde{\tilde{v}},m)$ where $\tilde{\tilde{v}}=\frac{(1-\tilde{\tilde{u}})\tilde{\tilde{u}}}{(w_{1}+\tilde{\tilde{u}}-1)},(1-{w_{1}})<\tilde{\tilde{u}}<1,w_{1}<1$ and $\tilde{\tilde{u}}$ is the real positive root of following equation:

$\nu_{1}u^{3}+\nu_{2}u^{2}+\nu_{3}u+\nu_{4}=0,$ (58)

where $\nu_{1}=w_{2},\,\nu_{2}=-[a_{2}w_{1}+w_{2}(-w_{1}+(2+m))],\nu_{3}=[a_{2}w_{1}-w_{1}w_{2}+w_{2}+mw_{1}w_{3}-m(-a_{2}w_{1}+2w_{2}(w_{1}-1))],\,\nu_{4}=m(w_{1}-1)[w_{2}+w_{1}(w_{3}-(w_{2}-a_{2}))].$

Coexistence of Species with interior equilibrium population: $E_{8}(u^{*},v^{*},r^{*})$ is the solution of following equations:


	$\displaystyle 1-u-\frac{w_{1}v}{u+v}=0,$		(59a)
	$\displaystyle-a_{2}+\frac{w_{2}u}{u+v}-\frac{w_{3}r}{r+v}=0,$		(59b)
	$\displaystyle c-\frac{w_{4}r}{v+D_{3}}=0.$		(59c)

From (59) we have,

	$\displaystyle v=\frac{(1-u)u}{w_{1}+u-1},$		(60)
	$\displaystyle r=\frac{c(v+D_{3})}{w_{4}}.$		(61)

Putting the values of $v$ and $r$ from (60) and (61) into (59b), we get

\displaystyle\alpha_{1}u^{3}-\alpha_{2}u^{2}-\alpha_{3}u-\alpha_{4}=0,

(62)

where
$\alpha_{1}=w_{2}(w_{4}+c)$ ,
$\alpha_{2}=[w_{4}(a_{2}w_{1}-w_{1}w_{2}+2w_{2})+c(w_{2}(2+D_{3})+w_{1}(w_{3}+a_{2}-w_{2}))]$ ,
$\alpha_{3}=[-w_{2}(w_{4}+c)-cD_{3}w_{1}(w_{3}+(a_{2}-2w_{2}))-(w_{1}(a_{2}-w_{2})w_{4}+c(2D_{3}w_{2}+w_{1}(w_{3}+a_{2}-w_{2})))]$ ,
$\alpha_{4}=-cD_{3}(w_{1}-1)[w_{2}+w_{1}(w_{3}+a_{2}-w_{2})]$ .
$u^{*}$ is the real positive root of (62). Hence interior equilibrium point $E_{8}(u^{*},v^{*},r^{*})$ exists if $({1}-{w_{1}})<u<{1},w_{1}<1$ . Knowing the value of $u^{*}$ we can find the values of $v^{*}$ and $r^{*}$ from (60) and (61) respectively.

Remark 5.1

From a realistic biological point of view, we are only interested in the dynamical behavior of model system (2)-(4) around the positive interior equilibrium point $E_{8}(u^{*},v^{*},r^{*})$ . This ensures the coexistence of all three species.

Remark 5.2

Equilibrium points $E_{0},E_{1},E_{2},E_{3}$ exist, even though an indeterminate form appears in (56), (57) due to the ratio dependent functional response. This has been rigorosly proven in [12].

5.2 Turing conditions

In this subsection we derive sufficient conditions for Turing instability to occur in (2)-(4) where the positive interior equilibrium point $E_{8}(u^{*},v^{*},r^{*})$ is stable in the absence of diffusion and unstable due to the addition of diffusion, under a small perturbation to $E_{8}(u^{*},v^{*},r^{*})$ . This is achieved by first linearizing model (2)-(4) about the homogenous steady state, by introducing both space and time-dependent fluctuations around $E_{8}(u^{*},v^{*},r^{*})$ . This is given as


	$\displaystyle u=u^{*}+\hat{u}(\xi,t),$		(63a)
	$\displaystyle v=v^{*}+\hat{v}(\xi,t),$		(63b)
	$\displaystyle r=r^{*}+\hat{r}(\xi,t),$		(63c)

where $|\hat{u}(\xi,t)|\ll u^{*}$ , $|\hat{v}(\xi,t)|\ll v^{*}$ and $|\hat{r}(\xi,t)|\ll r^{*}$ . Conventionally we choose

\left[{\begin{array}[]{cc}\hat{u}(\xi,t)\\ \hat{v}(\xi,t)\\ \hat{r}(\xi,t)\end{array}}\right]=\left[{\begin{array}[]{cc}\epsilon_{1}\\ \epsilon_{2}\\ \epsilon_{3}\end{array}}\right]e^{\lambda t+ik\xi},

where $\epsilon_{i}$ for $i=1,2,3$ are the corresponding amplitudes, $k$ is the wavenumber, $\lambda$ is the growth rate of perturbation in time $t$ and $\xi$ is the spatial coordinate. Substituing (63) into (2)-(4) and ignoring higher order terms including nonlinear terms, we obtain the characteristic equation which is given as

\displaystyle({\bf J}-\lambda{\bf I}-k^{2}{\bf D})\left[{\begin{array}[]{cc}\epsilon_{1}\\ \epsilon_{2}\\ \epsilon_{3}\end{array}}\right]=0,

(67)

where

\quad\bf{D}=\left[{\begin{array}[]{ccc}d_{1}&0&0\\ 0&d_{2}&0\\ 0&0&d_{3}\end{array}}\right],

$\bf{J}=\begin{bmatrix}u^{*}\left(-1+\frac{w_{1}v^{*}}{(u^{*}+v^{*})^{2}}\right)&-\frac{w_{1}u^{*^{2}}}{(u^{*}+v^{*})^{2}}&0\\ \frac{w_{2}v^{*^{2}}}{(u^{*}+v^{*})^{2}}&v^{*}\left(-\frac{w_{2}u^{*}}{(u^{*}+v^{*})^{2}}+\frac{w_{3}r^{*}}{(v^{*}+r^{*})^{2}}\right)&-\frac{w_{3}v^{*^{2}}}{(v^{*}+r^{*})^{2}}\\ 0&\frac{r^{*^{2}}(r^{*}-m)w_{4}}{(v^{*}+D_{3})^{2}}&-\frac{r^{*}(r^{*}-m)w_{4}}{(v^{*}+D_{3})}\\ \end{bmatrix}=\begin{bmatrix}J_{11}&J_{12}&J_{13}\\ J_{21}&J_{22}&J_{23}\\ J_{31}&J_{32}&J_{33}\\ \end{bmatrix},$

and $\bf{I}$ is a $3\times 3$ identity matrix.
For the non-trivial solution of (67), we require that

\left|\begin{array}[]{ccc}J_{11}-\lambda-k^{2}d_{1}&J_{12}&J_{13}\\ J_{21}&J_{22}-\lambda-k^{2}d_{2}&J_{23}\\ J_{31}&J_{32}&J_{33}-\lambda-k^{2}d_{3}\\ \end{array}\right|=0,

which gives a dispersion relation corresponding to $E_{8}(u^{*},v^{*},r^{*})$ . To determine the stability domain associated with $E_{8}(u^{*},v^{*},r^{*})$ , we rewrite the dispersion relation as a cubic polynomial function given as

\displaystyle P(\lambda(k^{2}))=\lambda^{3}+\boldsymbol{\mu_{2}}(k^{2})\lambda^{2}+\boldsymbol{\mu_{1}}(k^{2})\lambda+\boldsymbol{\mu_{0}}(k^{2}),

(68)

with coefficients

	$\displaystyle\boldsymbol{\mu_{2}}(k^{2})$	$\displaystyle=(d_{1}+d_{2}+d_{3})k^{2}-(J_{11}+J_{22}+J_{33}),$
	$\displaystyle\boldsymbol{\mu_{1}}(k^{2})$	$\displaystyle=J_{11}J_{33}+J_{11}J_{22}+J_{22}J_{33}-J_{32}J_{23}-J_{12}J_{21}-k^{2}\big{(}(d_{3}+d_{1})J_{22}+(d_{2}+d_{1})J_{33}+(d_{2}+d_{3})J_{11}\big{)}$
		$\displaystyle+k^{4}(d_{2}d_{3}+d_{2}d_{1}+d_{1}d_{3}),$
	$\displaystyle\boldsymbol{\mu_{0}}(k^{2})$	$\displaystyle=J_{11}J_{32}J_{23}-J_{11}J_{22}J_{33}+J_{12}J_{21}J_{33}+k^{2}\big{(}d_{1}(J_{22}J_{33}-J_{32}J_{23})+d_{2}J_{11}J_{33}+d_{3}(J_{22}J_{11}-J_{12}J_{21})\big{)}$
		$\displaystyle-k^{4}\big{(}d_{2}d_{1}J_{33}+d_{1}d_{3}J_{22}+d_{2}d_{3}J_{11}\big{)}+k^{6}d_{1}d_{2}d_{3}.$

According to Routh-Hurwitz criterion for stability, $\mathbb{R}e(\lambda(k^{2}))<0$ in model (2)-(4) around equilibrium point $E_{8}(u^{*},v^{*},r^{*})$ (i.e $E_{8}$ is stable) if and only if these conditions holds:

\displaystyle\boldsymbol{\mu_{2}}(k^{2})>0,\,\boldsymbol{\mu_{1}}(k^{2})>0,\,\boldsymbol{\mu_{0}}(k^{2})>0\quad\text{and}\quad[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2})>0.

(69)

As violating either of the above conditions implies instability (i.e $\mathbb{R}e(\lambda(k^{2}))>0$ ).
We now require conditions where an homogenous steady state ( $E_{8}(u^{*},v^{*},r^{*})$ ) will be stable to small perturbation in the absence of diffusion and unstable in the present of diffusion with certain $k$ vaules. Meaning, we require that around the homogenous steady state $E_{8}(u^{*},v^{*},r^{*})$

\displaystyle\mathbb{R}e(\lambda(k^{2}>0))>0,\,\text{for some}\,k\,\text{and}\,\mathbb{R}e(\lambda(k^{2}=0))<0,

where we consider $k$ to be real and positive even though $k$ can be complex. This behavior is called diffusion driven instability. Models that exhibits this sort of behavior in $2$ and $3$ species has been extensively studied in [33], where several different patterns was observed.
In order for homogenous steady state $E_{8}(u^{*},v^{*},r^{*})$ to be stable( in the absence of diffusion) we need

\displaystyle\boldsymbol{\mu_{2}}(k^{2}=0)>0,\,\boldsymbol{\mu_{1}}(k^{2}=0)>0,\,\boldsymbol{\mu_{0}}(k^{2}=0)>0\quad\text{and}\quad[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2}=0)>0,

where as with diffusion ( $k^{2}>0$ ) we look for conditions where we can drive the homogenous steady state to be unstable, this can be achieved by studying the coefficient of (68). In order to achieve this we need to reverse at least one of the signs in (69). By which we first study $\boldsymbol{\mu_{2}}(k^{2})$ . Irrespective of the value of $k^{2}$ , $\boldsymbol{\mu_{2}}(k^{2})$ will be positive since $J_{11}+J_{22}+J_{33}$ is always less than zero. Therefore we cannot depend on $\boldsymbol{\mu_{2}}(k^{2})$ for diffusion driven instability to occur. Hence for diffusion driven instability to occur in our case, we only depend on the $2$ conditions which are

\displaystyle\boldsymbol{\mu_{0}}(k^{2})\quad\text{and}\quad[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2}).

(70)

Both functions are cubic functions of $k^{2}$ , which are generally of the form

\displaystyle G(k^{2})=H_{H}+k^{2}D_{D}+(k^{2})^{2}C_{C}+(k^{2})^{3}B_{B},\,\text{with}\,B_{B}>0,\,\text{and}\,H_{H}>0.

The coefficient of $G(k^{2})$ are given in Table 1.

Table 1: Coefficients of cubic functions

\boldsymbol{\mu_{0}}(k^{2})

and

[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2})

used in determining conditions for Turing instability

Coefficient of $G(k^{2})$	$\boldsymbol{\mu_{0}}(k^{2})$	$[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2})$
$H_{H}$	$J_{11}J_{32}J_{23}+J_{12}J_{21}J_{33}$	$J_{11}J_{22}J_{33}-(J_{11}+J_{22}+J_{33})(J_{11}J_{22}-J_{12}J_{21}+J_{11}J_{33}+J_{22}J_{33}-J_{23}J_{32})$
$H_{H}$	$-J_{11}J_{22}J_{33}$	$-J_{11}J_{23}J_{32}-J_{12}J_{21}J_{33}$
$D_{D}$
	$d_{1}(J_{22}J_{33}-J_{32}J_{23})+d_{2}J_{11}J_{33}$	$d_{1}(2J_{11}J_{33}+2J_{11}J_{22}+2J_{22}J_{33}+J_{33}J_{33}+J_{22}J_{22}-J_{12}J_{21})$
	$+d_{3}(J_{11}J_{22}-J_{12}J_{21})$	$+d_{2}(2J_{22}J_{11}+2J_{22}J_{33}+2J_{33}J_{11}+J_{11}J_{11}+J_{33}J_{33}-J_{21}J_{12}-J_{23}J_{32})$
		$d_{3}(2J_{22}J_{11}+2J_{22}J_{33}+2J_{33}J_{11}+J_{11}J_{11}+J_{22}J_{22}-J_{23}J_{32})$
$C_{C}$
	$-d_{1}d_{2}J_{33}-d_{1}d_{3}J_{22}-d_{2}d_{3}J_{11}$	$-J_{11}(d_{2}+d_{3})(2d_{1}+d_{2}+d_{3})-J_{22}(d_{1}+d_{3})(d_{1}+2d_{2}+d_{3})$
		$-J_{33}(d_{1}+d_{2})(d_{1}+d_{2}+2d_{3})$
$B_{B}$
	$d_{1}d_{2}d_{3}$	$(d_{2}+d_{3})(d_{1}d_{1}+d_{2}d_{3}+d_{1}d_{2}$
		$+d_{1}d_{3})$

To drive either $\boldsymbol{\mu_{0}}(k^{2})$ or $[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2})$ to negative for some $k$ , we need to find the minimum $k^{2}$ referred to as the minimum turing point ( $k^{2}_{T}$ ) such that $G(k^{2}=k^{2}_{T})<0$ . This minimum turing point occurs when

{\partial G/\partial(k^{2})}=0,

this when solved for $k^{2}$ yields

\displaystyle k^{2}=k^{2}_{T}={-C_{C}+\sqrt{C_{C}^{2}-3B_{B}D_{D}}\over 3B_{B}}.

To ensures $k^{2}$ is real and positive such that ${\partial^{2}G/\partial^{2}(k^{2})}>0$ , we require either

\displaystyle D_{D}<0\quad\text{or}\quad C_{C}<0,

(71)

which ensures that

\displaystyle C_{C}^{2}-3B_{B}D_{D}>0.

Therefore $G(k^{2})<0$ , if at $k^{2}=k^{2}_{T}$

\displaystyle G_{min}(k^{2})=2C_{C}^{3}-9D_{D}C_{C}B_{B}-2(C_{C}^{2}-3D_{D}B_{B})^{3/2}+27B_{B}H_{H}^{2}<0.

(72)

Hence (71)-(72) are necessary and sufficient conditions for $E_{8}(u^{*},v^{*},r^{*})$ to produce diffusion driven instability, leading to the emergence of patterns. Also to first establish stability when $k=0$ , $H_{H}$ in each case has to be positive. We state this via the following theorem

Theorem 5.3

Consider the three species food chain model described via (2)-(4), with equilibrium point $E_{8}(u^{*},v^{*},r^{*})$ . Supposing we have a set of parameters such that the following conditions are satisfied,

1.

$H_{H}>0$ when $k=0$ ,
2.

$D_{D}<0$ or $C_{C}<0$ ( and $C_{C}^{2}>B_{B}D_{D}$ ),
3.

$G_{min}(k^{2})=2C_{C}^{3}-9D_{D}C_{C}B_{B}-2(C_{C}^{2}-3D_{D}B_{B})^{3/2}+27B_{B}H_{H}^{2}<0.$

Then Turing instability will always occur in the model.

See [33] for the connection between the signs of the coefficients in table 1 and the associated patterns.

5.3 Effect of Allee threshold parameter $m$ on Turing Instability

In this section we investigate the effects of Allee threshold $m$ on Turing patterns. We focus on whether $m$ can induce Turing instabilities. We perform numerical simulation of model system (2)-(4) and the use of coefficient of dispersion relation, given in Table 1.
A Chebychev collocation scheme is employed for the approximation of the time varying, one-dimensional problem[11, 9].

Our one-dimensional numerical simulation employs the zero-flux boundary condition with a spatial domain of size $(256\times 256)$ . We define the spatial domain to be $X\in(0,L)$ , in the nondimensionalistaion in appendix 7.1, where $\quad x={X\pi\over L}$ .

The initial condition as defined is a small perturbation around the positive homogenous steady state given as

	$\displaystyle u=u^{*}+\epsilon_{1}cos^{2}(nx)(x>0)(x<\pi),$
	$\displaystyle v=v^{*}+\epsilon_{2}cos^{2}(nx)(x>0)(x<\pi),$
	$\displaystyle r=r^{*}+\epsilon_{3}cos^{2}(nx)(x>0)(x<\pi),$

where $\epsilon_{i}=0.05$ $\forall i$ .

In the two-dimensional case, we numerically solve model system (2)-(4) using a finite difference method. A forward difference scheme is used for the reaction part. A standard five-point explicit finite difference scheme is used for the two-dimensional diffusion terms. Model system (2)-(4) is numerically solved in two-dimension over $200\times 200$ mesh points with spatial resolution $\Delta x=\Delta y=0.1$ and time step $\Delta t=0.1$ .
In our numerical simulation both in the one-dimensional and two-dimensional case, with specific parameter set as show in figure (3)-figure (7), different types of dynamics are observed as a result of the value of the Allee threshold $m$ .
From figure (3)-figure (4), we numerically observe that change in $m$ can result in disappearance of patterns or the appearance of patterns. Specifically from figure (3)-figure (4), patterns where observed when $m=0.01$ but these patterns change as $m$ gets larger and finally at $m=0.5$ , these patterns completely disappear. The critical Allee threshold $m$ , for patterns to disappear, depends on the parameter set being used. In some cases, patterns where observed when $m=0$ but as $m\to 0.2$ and beyond, these patterns disappeared.

It is also conclusive from figure (5)-figure (7) that Allee threshold has an effect on the type of patterns that do form particularly, as in the case of figure (5) and (7), changing $m$ , can change spatio-temporal patterns to a fixed spatial patterns. This behavior can be observed using the dispersion relation as shown in figure (9). In choosing $\mu_{0}(k^{2})<0$ for some $k$ , which will mostly result in fixed spatial patterns, changing the value of $m$ does not result in a change in the type of patterns formed, but parameters move out of the Turing space. Therefore in our numerical experiments we do not find parameters such that fixed spatial patterns can change to spatial-temporal patterns, as $m$ increases from zero.

Table (2), also displays the fact that no patterns were observed as $m$ increases from $m=0.01$ to $m=0.5$ from figure (3) and figure (4) because $m$ can destabilizes an already stable equilibrium point in the ode model as in when $k=0$ . Hence as in figure (3) and figure (4) where $m=0.01$ produces patterns since $E_{8}(u^{*},v^{*},r^{*})$ is stable when $k=0$ , this stable ( $+$ ) equilibrium point becomes unstable ( $-$ ) even when $k=0$ .

$m$	$H_{H}=\mu_{0}(k^{2}=0)$	$H_{H}=[\mu_{2}\mu_{1}-\mu_{0}](k^{2}=0)$	Pattern
$0$	$+$	$+$	Patterns may occur
$0.1$	$+$	$+$	Patterns may occur
$0.2$	$-$	$-$	No Patterns
$0.3$	$-$	$-$	No Patterns
$0.4$	$-$	$-$	No Patterns
$0.5$	$-$	$-$	No Patterns

Table 2: The influence of

m

on the sign of

\boldsymbol{\mu_{0}}(k^{2})

and

[\boldsymbol{\mu_{2}}\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{0}}](k^{2})

. Describing the stability of

E_{8}(u^{*},v^{*},r^{*})

when

k=0

and how its might lead to pattern formation. The parameters are:

w_{1}=0.96,w_{2}=0.52,w_{3}=1.06,w_{4}=0.37,a_{2}=0.014,D_{3}=0.1,c=0.1,d_{1}=10^{-3},d_{2}=10^{-5},d_{3}=10^{-5}

6 Discussion

In this manuscript we have considered a spatially explicit ratio dependent three species food chain model, with a strong Allee effect in the top predator. The model represents dynamics of three interacting species in “well-mixed” conditions and is applicable to ecological problems in terrestrial, as well as aquatic systems. Our goal is to understand the effect of the Allee threshold on the system dynamics, as there is not much work on Allee effects in multi-trophic level food chains, and we hope to address this here.

We first prove the existence of a global attractor for the model. In many systems an interesting question is to consider the effect on the global attractor, under perturbation of the physical parameters in the system. Most often continuity cannot be proven, but upper semi-continuity can. We show that the global attractor is upper semi-continuous w.r.t the Allee threshold parameter $m$ , that is $dist_{E}(\mathcal{A}_{m},\mathcal{A}_{0})\rightarrow 0,\ \mbox{as}\ m\rightarrow 0^{+}$ , via theorem 3.3. To the best of our knowledge this is the first robustness result for a three-species model with strong Allee effect. Unless a systems dynamics are robust, there is no possibility to capture the same ecological behavior in a laboratory experiment, or natural setting.

The next question of interest, after one has proved such a result, is estimating the rate of decay, in terms of the physical parameter in question. Since there are no theoretical results to estimate the decay rate to the target attractor (that is when $m=0$ ), we choose to investigate this rate numerically. We find decay estimates of the order of $\mathcal{O}(m^{\gamma})$ , where $\gamma$ is found explicitly to be slightly less, but very close to 1.

We also find that $m$ effects the spatiotemporal dynamics of the system in two distinct ways.
(1) It changes the Turing patterns that occur in the system, in both one and two spatial dimensions. That is varying $m$ , has a distinct effect on the sorts of Turing patterns that form. This can have interesting effects on the patchiness of spatially dispersing species. Our results say that if one can effect the Allee threshold in the top predator, one can effectively change the spatial concentrations of the species involved. This is best visualized in figure 8. This could have many potential applications in Allee mediated biological control, and conservation, as recently biologists have begun to consider Allee effects as beneficial in limiting establishment of an invading species [28].
(2) It facilitates overexploitation phenomenon. That is without $m$ , or when $m=0$ , one does not see extinction in top predator, as the origin becomes unstable. However, introducing it stabilises the origin, and low initial concentrations of $r$ , may yield extinction, via theorems 4.1. This has many possible applications to top down trophic cascades and cascade control. This also tells us that in a three species system, where the top predator can switch its food source, only having an Allee effect in place can possibly drive it to extinction, thus bring about true overexploitation. This could have potential applications in Allee mediated biological control, if the top predator is an invasive species.

As future investigations it would be interesting to try to model weak Allee effects in top predator, and/or Allee effects in the other species as well, be they weak or strong. It has been stipulated with evidence, that two or more Allee effects can occur simultaneously in the same population[4]. Thus an extremely interesting question would be to consider different type of Allee effects in the different species, and investigate the interplay of the various Allee thresholds, as they affect the dynamics of the system. All in all our results would be of interest to both the mathematical and ecological communities, and in particular to groups that are interested in conservation efforts in food chain systems, cascade control, such as in many aquatic systems, and even Allee mediated biological control.

References

[1] Peter A Abrams and Lev R Ginzburg. The nature of predation: prey dependent, ratio dependent or neither? Trends in Ecology & Evolution, 15(8):337–341, 2000.
[2] Warder Clyde Allee. Co-operation among animals. American Journal of Sociology, pages 386–398, 1931.
[3] Warder Clyde Allee, Orlando Park, Alfred Edwards Emerson, Thomas Park, Karl P Schmidt, et al. Principles of animal ecology. Number Edn 1. WB Saundere Co. Ltd., 1949.
[4] Luděk Berec, Elena Angulo, and Franck Courchamp. Multiple allee effects and population management. Trends in Ecology & Evolution, 22(4):185–191, 2007.
[5] Mathematical Bioeconomics. The optimal management of renewable resources, 1990.
[6] Stephen R Carpenter, James F Kitchell, and James R Hodgson. Cascading trophic interactions and lake productivity. BioScience, pages 634–639, 1985.
[7] Colin W Clark. The worldwide crisis in fisheries: economic models and human behavior. Cambridge University Press, 2006.
[8] JM Drake and AM Kramer. Allee effects. Nat. Educ. Knowl, 3(10):2, 2011.
[9] David Gottlieb and Liviu Lustman. The spectrum of the chebyshev collocation operator for the heat equation. SIAM journal on numerical analysis, 20(5):909–921, 1983.
[10] Dan Henry. Geometric theory of semilinear parabolic equations, volume 840. Springer-Verlag New York, 1981.
[11] Jan S Hesthaven, Sigal Gottlieb, and David Gottlieb. Spectral methods for time-dependent problems, volume 21. Cambridge University Press, 2007.
[12] Yang Kuang and Edoardo Beretta. Global qualitative analysis of a ratio-dependent predator–prey system. Journal of Mathematical Biology, 36(4):389–406, 1998.
[13] Martin Liermann and Ray Hilborn. Depensation: evidence, models and implications. Fish and Fisheries, 2(1):33–58, 2001.
[14] Andrew Morozov, Sergei Petrovskii, and Bai-Lian Li. Bifurcations and chaos in a predator-prey system with the allee effect. Proceedings of the Royal Society of London. Series B: Biological Sciences, 271(1546):1407–1414, 2004.
[15] James D Murray. Mathematical biology i: An introduction, vol. 17 of interdisciplinary applied mathematics, 2002.
[16] Akira Okubo and Smon A Levin. Diffusion and ecological problems: modern perspectives, volume 14. Springer Science & Business Media, 2013.
[17] Pallav Jyoti Pal and Prashanta Kumar Mandal. Bifurcation analysis of a modified leslie–gower predator–prey model with beddington–deangelis functional response and strong allee effect. Mathematics and Computers in Simulation, 97:123–146, 2014.
[18] Rana D Parshad, Hamid Ait Abderrahmane, Ranjit Kumar Upadhyay, and Nitu Kumari. Finite time blowup in a realistic food-chain model. ISRN Biomathematics, 2013, 2013.
[19] Rana D Parshad, Nitu Kumari, Aslan R Kasimov, and Hamid Ait Abderrahmane. Turing patterns and long-time behavior in a three-species food-chain model. Mathematical biosciences, 254:83–102, 2014.
[20] Rana D Parshad and Ranjit K Upadhyay. Investigation of long time dynamics of a diffusive three species aquatic model. Dynamics of Partial Differential Equations, 7(3):217–244, 2010.
[21] Sergei V Petrovskii, Andrew Y Morozov, and Ezio Venturino. Allee effect makes possible patchy invasion in a predator–prey system. Ecology Letters, 5(3):345–352, 2002.
[22] George R Sell and Yuncheng You. Dynamics of evolutionary equations, volume 143. Springer Science & Business Media, 2013.
[23] Moitri Sen, Malay Banerjee, and Andrew Morozov. Bifurcation analysis of a ratio-dependent prey–predator model with the allee effect. Ecological Complexity, 11:12–27, 2012.
[24] Swarnali Sharma and GP Samanta. A ratio-dependent predator-prey model with allee effect and disease in prey. Journal of Applied Mathematics and Computing, 47(1-2):345–364, 2014.
[25] Joel Smoller. Shock waves and reaction?diffusion equations, volume 258. Springer Science & Business Media, 1994.
[26] Philip A Stephens and William J Sutherland. Consequences of the allee effect for behaviour, ecology and conservation. Trends in Ecology & Evolution, 14(10):401–405, 1999.
[27] Roger Temam. Infinite-dimensional dynamical systems in mechanics and physics, volume 68. Springer Science & Business Media, 2012.
[28] Patrick C Tobin, Luděk Berec, and Andrew M Liebhold. Exploiting allee effects for managing biological invasions. Ecology letters, 14(6):615–624, 2011.
[29] Alan Mathison Turing. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 237(641):37–72, 1952.
[30] Ranjit Kumar Upadhyay and Vikas Rai. Why chaos is rarely observed in natural populations. Chaos, Solitons & Fractals, 8(12):1933–1939, 1997.
[31] George AK Van Voorn, Lia Hemerik, Martin P Boer, and Bob W Kooi. Heteroclinic orbits indicate overexploitation in predator–prey systems with a strong allee effect. Mathematical biosciences, 209(2):451–469, 2007.
[32] Jinfeng Wang, Junping Shi, and Junjie Wei. Dynamics and pattern formation in a diffusive predator–prey system with strong allee effect in prey. Journal of Differential Equations, 251(4):1276–1304, 2011.
[33] KAJ White and CA Gilligan. Spatial heterogeneity in three species, plant–parasite–hyperparasite, systems. Philosophical Transactions of the Royal Society B: Biological Sciences, 353(1368):543–557, 1998.

7 Appendix

7.1 Nondimensionalisation

The model system in (2)-(4) is a nondimensionalised version of the following system


$\displaystyle{\partial U\over{\partial T}}$	$\displaystyle=D_{U}\Delta U+A_{1}U-B_{1}U^{2}-\frac{W_{1}UV}{\beta_{1}V+U},$	(73a)
$\displaystyle{\partial V\over{\partial T}}$	$\displaystyle=D_{V}\Delta V-A_{2}V+\frac{W_{2}UV}{\beta_{1}V+U}-\frac{W_{3}VR}{\beta_{3}R+V},$	(73b)
$\displaystyle{\partial R\over{\partial T}}$	$\displaystyle=D_{R}\Delta R+R\big{(}{{R-M}}\big{)}\bigg{(}C-\frac{W_{4}R}{V+A_{3}}\bigg{)}.$	(73c)

$A_{1}$ is the self-growth rate of the prey population $U$ . $A_{2}$ the intrinsic death rate of the specialist predator $V$ in the absence of its only food $U$ , $C$ measures the rate of self-reproduction of the generalist predator $R$ . $B_{1}$ measures the strength of competition among individuals of the prey species $U$ . $W_{i}^{\prime}s$ are the maximum values which per capita growth rate can attain. $A_{3}$ represents the residual loss in $R$ population due to severe scarcity of its favorite food $V$ . $\beta_{1}$ is the parameters that describes the handling time of the prey $U$ by predator $V$ , and $\beta_{3}$ is the parameter that describes the handling time of the prey $V$ by predator $R$ .

We will now go over the details of the nondimensionalisation. We consider (73), and aim to introduce a change of variables and time scaling which reduces the number of parameters of model system (73)‎. We take

	$\displaystyle u$	$\displaystyle={UB_{1}\over A_{1}},\quad v={VB_{1}\beta_{1}\over A_{1}},\quad r={Rb_{1}\beta_{1}\beta_{3}\over A_{1}},$
	$\displaystyle t$	$\displaystyle={T\over{A_{1}}},\quad w_{1}={W_{1}\over\beta_{1}A_{1}},\quad a_{2}={A_{2}\over A_{1}},$
	$\displaystyle w_{2}$	$\displaystyle={W_{2}\over A_{1}},\quad w_{3}={W_{3}\over\beta_{3}A_{1}},\quad c={CA_{1}\over{A_{1}B_{1}\beta_{1}\beta_{3}}},$
	$\displaystyle m$	$\displaystyle={MB_{1}\beta_{1}\beta_{3}\over A_{1}},\quad D_{3}={A_{3}B_{1}\beta_{1}\over A_{1}},\quad w_{4}={W_{4}B_{1}\over{B_{1}^{2}\beta_{1}\beta_{3}^{2}}},$
	$\displaystyle d_{1}$	$\displaystyle={D_{U}\pi^{2}\over B_{1}L^{2}},\quad d_{2}={D_{V}\pi^{2}\over B_{1}\beta_{1}L^{2}},\quad d_{3}={D_{R}\pi^{2}\over{B_{1}\beta_{1}\beta_{3}L^{2}}}.$

Considering Neuman boundary conditions, model system (73) reduces to

	$\displaystyle\frac{\partial u}{\partial t}=d_{1}\Delta u+u-u^{2}-w_{1}\frac{uv}{u+v},$		(74)
	$\displaystyle\frac{\partial v}{\partial t}=d_{2}\Delta v-a_{2}v+w_{2}\frac{uv}{u+v}-w_{3}\left(\frac{vr}{v+r}\right),$		(75)
	$\displaystyle\frac{\partial r}{\partial t}=d_{3}\Delta r+r\big{(}{{r-m}}\big{)}\bigg{(}c-\frac{w_{4}r}{v+D_{3}}\bigg{)}.$		(76)

Also all parameters associated with model system (74-76) are assumed to be positive constants and have the usual biological meaning.

7.2 Global existence

Here we prove proposition 2.2. The positivity of solutions follow trivially from the form of the reaction terms, which are quasi positive. By the positivity of the first component $u(t,.)$ of the solution to (2) on $[0,T_{\max})$ $\times\Omega$ , we get from equation (2)

\frac{\partial u}{\partial t}-d_{1}\Delta u\leq a_{1}u,\text{ \ on }[0,T_{\max})\times\Omega,

(77)

then we use a comparison argument [25]. That is, we can compare the solution $u$ of (77) to the solution $u^{*}$ , of the linear heat equation

\frac{\partial u^{*}}{\partial t}-d_{1}\Delta u^{*}=a_{1}u^{*},\text{ \ on }[0,T_{\max})\times\Omega,

(78)

where $u^{*}$ satisfies the same initial and boundary conditions as $u$ . Clearly $u\leq u^{*}$ , and since $u^{*}$ (being the solution of a linear heat equation) is bounded, so is $u$ , and we deduce

u(t,.)\leq C_{1},

(79)

where $C_{1}$ is a constant depending only on $T_{\max}$ . Then equation (3) gives

\frac{\partial v}{\partial t}-d_{2}\Delta v\leq w_{1}C_{1}v,\text{ \ on }[0,T_{\max})\times\Omega,

(80)

which implies by the same arguments

v(t,.)\leq C_{2},\text{ \ on }[0,T_{\max})\times\Omega,

(81)

where $C_{2}$ is a constant depending only on $T_{\max}$ . To get a bound on the component $r$ we apply Young’s inequality to yield

	$\displaystyle\frac{\partial r}{\partial t}-d_{3}\Delta r+Mcr$
	$\displaystyle=\left(c+\frac{Mw_{3}}{v+D_{3}}\right)r^{2}-\left(\frac{w_{3}}{v+D_{3}}\right)r^{3}$
	$\displaystyle\leq\left(c+\frac{Mw_{3}}{D_{3}}\right)r^{2}-\left(\frac{w_{3}}{v+D_{3}}\right)r^{3}$
	$\displaystyle=\left(c+\frac{Mw_{3}}{D_{3}}\right)\left(\frac{v+D_{3}}{w_{3}}\right)^{\frac{2}{3}}\left(\frac{v+D_{3}}{w_{3}}\right)^{-\frac{2}{3}}r^{2}-\left(\frac{w_{3}}{v+D_{3}}\right)r^{3}$
	$\displaystyle\leq\left(c+\frac{Mw_{3}}{D_{3}}\right)^{3}\left(\frac{v+D_{3}}{w_{3}}\right)^{2}+\left(\frac{w_{3}}{v+D_{3}}\right)r^{3}-\left(\frac{w_{3}}{v+D_{3}}\right)r^{3}$
	$\displaystyle\leq\left(c+\frac{Mw_{3}}{D_{3}}\right)^{3}\left(\frac{C_{2}+D_{3}}{w_{3}}\right)^{2}$		(82)

Here we use the bound on $v$ via (81). Thus there exists a positive constant $C_{3}$ such that

r(t,.)\leq C_{3}=\left(c+\frac{Mw_{3}}{D_{3}}\right)^{3}\left(\frac{C_{2}+D_{3}}{w_{3}}\right)^{2},\text{ \ on }[0,T_{\max})\times\Omega.

(83)

At this step via (79), (81), (83) standard theory [10] is applicable and the solution is global ( $T_{\max}=+\infty$ ). This proves the proposition.

7.3 Apriori estimates

In all estimates made hence forth the constants $C,C_{1},C_{2},C_{3},C_{\infty}$ are generic constants, that can change in value from line to line, and sometimes within the same line if so required. We estimate the gradient of $r$ by multiplying (4) by $-\Delta r$ , and integrating by parts and using standard methods to obtain

\frac{1}{2}\frac{d}{dt}||\nabla r||^{2}_{2}+d_{3}||\Delta r||^{2}_{2}=\int_{\Omega}r(r-m)(c-w_{4}\frac{r}{v+D_{3}})(-\Delta r)dx.

(84)

We have to estimate $\int_{\Omega}r(r-m)(c-w_{4}\frac{r}{v+D_{3}})(-\Delta r)dx$ . Note

	$\displaystyle\int_{\Omega}r(r-m)(c-\frac{r}{v+D_{3}})(-\Delta r)dx$	$\displaystyle=\int_{\Omega}\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}-mcr\right)(-\Delta r)dx,$
		$\displaystyle=\int_{\Omega}\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}\right)(-\Delta r)dx-cm\|\|\nabla r\|\|^{2}_{2}.$

Thus we obtain

\frac{1}{2}\frac{d}{dt}||\nabla r||^{2}_{2}+d_{3}||\Delta r||^{2}_{2}+cm||\nabla r||^{2}_{2}=\int_{\Omega}\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}\right)(-\Delta r)dx.

(85)

We now proceed in 2 steps. We first assume $\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}\right)\geq 0$ .

Young’s Inequality with epsilon then gives us the following estimate

	$\displaystyle\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}\right)$	$\displaystyle\leq\left(\frac{v+D_{3}}{w_{4}}\right)^{\frac{1}{3}}\left(c+\frac{mw_{4}}{v+D_{3}}\right)r\left(\frac{w_{4}}{v+D_{3}}\right)^{\frac{1}{3}}r-w_{4}\frac{r^{3}}{v+D_{3}},$
		$\displaystyle\leq\left(\frac{v+D_{3}}{w_{4}}\right)^{\frac{4}{9}}\left(c+\frac{mw_{4}}{v+D_{3}}\right)^{\frac{4}{3}}r^{\frac{4}{3}}+w_{4}\frac{r^{3}}{v+D_{3}}-w_{4}\frac{r^{3}}{v+D_{3}},$
		$\displaystyle\leq\left(\frac{\|\|v+D_{3}\|\|_{\infty}}{w_{4}}\right)\left(c+\frac{mw_{4}}{v+D_{3}}\right)^{\frac{4}{3}}r^{\frac{4}{3}}+w_{4}\frac{r^{3}}{v+D_{3}}-w_{4}\frac{r^{3}}{v+D_{3}},$
		$\displaystyle\leq\left(\frac{C+D_{3}}{w_{4}}\right)\left(c+\frac{mw_{4}}{D_{3}}\right)^{\frac{4}{3}}r^{\frac{4}{3}}+w_{4}\frac{r^{3}}{v+D_{3}}-w_{4}\frac{r^{3}}{v+D_{3}},$
		$\displaystyle=\left(\frac{C+D_{3}}{w_{4}}\right)\left(c+\frac{mw_{4}}{D_{3}}\right)^{\frac{4}{3}}r^{\frac{4}{3}}.$

Then we obtain

	$\displaystyle\frac{1}{2}\frac{d}{dt}\|\|\nabla r\|\|^{2}_{2}+d_{3}\|\|\Delta r\|\|^{2}_{2}+2cm\|\|\nabla r\|\|^{2}_{2}$	$\displaystyle=$	$\displaystyle\int_{\Omega}\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}\right)(-\Delta r)dx,$
		$\displaystyle\leq$	$\displaystyle\left(\frac{C+D_{3}}{w_{4}}\right)\left(c+\frac{mw_{4}}{D_{3}}\right)^{\frac{4}{3}}\int_{\Omega}r^{\frac{4}{3}}\|-\Delta r\|dx$

This follows via (7.3). Thus we obtain

\frac{1}{2}\frac{d}{dt}||\nabla r||^{2}_{2}+d_{3}||\Delta r||^{2}_{2}\leq C_{1}||r||^{\frac{8}{3}}_{\frac{8}{3}}+\frac{d_{3}}{2}||\Delta r||^{2}_{2}.

(88)

Here $C_{1}=\left(\frac{C+D_{3}}{w_{4}}\right)\left(c+\frac{w_{4}}{D_{3}}\right)^{\frac{4}{3}}$ . Now we can use the sobolev embedding of $H^{1}(\Omega)\hookrightarrow L^{\frac{8}{3}}(\Omega)$ , where $C_{3}$ is the embedding constant, to obtain

\frac{d}{dt}||\nabla r||^{2}_{2}+d_{3}||\Delta r||^{2}_{2}\leq C_{3}C_{1}\left(||\nabla r||^{2}_{2}\right)\left(||\nabla r||^{2}_{2}\right).

(89)

Thus

\frac{d}{dt}||\nabla r||^{2}_{2}\leq C_{3}C_{1}(||\nabla r||^{2}_{2})(||\nabla r||^{2}_{2}).

(90)

Now we invoke the estimate via (15), and the uniform Gronwall lemma to obtain

\mathop{\limsup}_{t\rightarrow\infty}||\nabla r||^{2}_{2}\leq C.

(91)

Now assume $\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}\right)<0$ . Then we obtain

			$\displaystyle\frac{1}{2}\frac{d}{dt}\|\|\nabla r\|\|^{2}_{2}+d_{3}\|\|\Delta r\|\|^{2}_{2}+cm\|\|\nabla r\|\|^{2}_{2},$
		$\displaystyle\leq$	$\displaystyle\int_{\Omega}\left(c+\frac{mw_{4}}{D_{3}}\right)r^{2}(-\Delta r)dx,$
		$\displaystyle\leq$	$\displaystyle\frac{2}{d_{3}}\left(c+\frac{mw_{4}}{D_{3}}\right)\|\|r\|\|^{4}_{4}+\frac{d_{3}}{2}\|\|\Delta r\|\|^{2}_{2}.$

We now use the Sobolev embedding of $H^{1}(\Omega)\hookrightarrow L^{4}(\Omega)$ to obtain

\frac{d}{dt}||\nabla r||^{2}_{2}+d_{3}||\Delta r||^{2}_{2}+cm||\nabla r||^{2}_{2}\leq\frac{2}{d_{3}}\left(c+\frac{w_{4}}{D_{3}}\right)||\nabla r||^{4}_{2},

(93)

dropping the $cm||\nabla r||^{2}_{2}$ term yields

\frac{d}{dt}||\nabla r||^{2}_{2}\leq\frac{2}{d_{3}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(||\nabla r||^{2}_{2}\right)\left(||\nabla r||^{2}_{2}\right),

(94)

\mathop{\limsup}_{t\rightarrow\infty}||\nabla r||^{2}_{2}\leq C,

(95)

7.4 Apriori estimates independent of the parameter $m$

Recall the following $L^{2}$ estimate

\frac{1}{2}\frac{d}{dt}||r||^{2}_{2}+d_{3}||\nabla r||^{2}_{2}+cm||r||^{2}_{2}+\int_{\Omega}\frac{w_{4}}{v+D_{3}}r^{4}dx\leq\int_{\Omega}\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{3}dx.

We then use Hölder’s inequality followed by Young’s inequality to obtain

	$\displaystyle\frac{1}{2}\frac{d}{dt}\|\|r\|\|^{2}_{2}+d_{3}\|\|\nabla r\|\|^{2}_{2}+cm\|\|r\|\|^{2}_{2}+\int_{\Omega}\frac{w_{4}}{v+D_{3}}r^{4}dx,$
	$\displaystyle\leq\frac{3}{4}\int_{\Omega}\frac{w_{4}}{v+D_{3}}r^{4}dx+\|\Omega\|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}.$

Now we drop the $cm||r||^{2}_{2}$ term from the left hand side, to avoid the singularity caused by $\frac{1}{m}$ , the $d_{3}||\nabla r||^{2}_{2}$ , and use the estimate on $||v||_{\infty}$ via (27), and the embedding of $L^{4}(\Omega)\hookrightarrow L^{2}(\Omega)$ to obtain

\frac{d}{dt}||r||^{2}_{2}+\frac{Cw_{4}}{C+D_{3}}||r||^{2}_{2}\leq 2|\Omega|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}.

(97)

Here, we assume $||r||_{2}\geq 1$ , else we already have a bound for $||r||_{2}$ . Next we use Gronwall’s inequality to obtain

||r||^{2}_{2}\leq e^{-\left(\frac{Cw_{4}}{C+D_{3}}\right)t}||r_{0}||^{2}_{2}+\frac{2|\Omega|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}.\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}}.

(98)

Thus for $t>t_{0}=\max(t^{*}+1,\frac{\ln(||r_{0}||^{2}_{2})}{\left(\frac{Cw_{4}}{C+D_{3}}\right)})$ , we obtain

||r||^{2}_{2}\leq 1+\frac{2|\Omega|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}},

(99)

which is a uniform $L^{2}(\Omega)$ bound in $r$ that is independent of the Allee parameter $m$ , time and initial data.

We next focus on making $H^{1}(\Omega)$ estimates, that are independent of the Allee parameter $m$ .

We choose the estimate derived in (99) for $||r||_{2}$ and insert this in (15) to obtain

	$\displaystyle\int^{t+1}_{t}\|\|\nabla r\|\|^{2}_{2}ds$
	$\displaystyle\leq 1+\frac{2\|\Omega\|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}}+2\|\Omega\|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}},$
	$\displaystyle\ \mbox{for}\ t>t_{0}.$		(100)

Now we choose

K_{m}=\max\left(\frac{2}{d_{3}}\left(c+\frac{mw_{4}}{D_{3}}\right),C_{3}\left(\frac{C+D_{3}}{w_{4}}\right)\left(c+\frac{mw_{4}}{D_{3}}\right)^{\frac{4}{3}}\right).

(101)

Then, we can majorise the right hand side of the above by plugging in $m=1$ . Since the above quantities are not singular in $m$ , this is possible. Thus we obtain

K=\max\left(\frac{2}{d_{3}}\left(c+\frac{w_{4}}{D_{3}}\right),C_{3}\left(\frac{C+D_{3}}{w_{4}}\right)\left(c+\frac{w_{4}}{D_{3}}\right)^{\frac{4}{3}}\right),

(102)

going back to the $H^{1}(\Omega)$ estimates in (91), and (95) one obtains,

\frac{d}{dt}||\nabla r||^{2}_{2}\leq K\left(||\nabla r||^{2}_{2}\right)\left(||\nabla r||^{2}_{2}\right).

(103)

Also note the integral in time estimate via (7.4) is now independent of $m$ . This allows us to apply the uniform Gronwall lemma on (103) with $\beta(t)=||\nabla r||^{2}_{2},\ \zeta(t)=K||\nabla r||^{2}_{2},\ h(t)=0,q=1$ , and the estimate via (7.4) to yield the following bound

	$\displaystyle\|\|\nabla r\|\|^{2}_{2}$
	$\displaystyle\leq C_{K}e^{KC_{K}},\ \mbox{for}\ t>t_{1}=t_{0}+1.$		(104)

Thus the $H^{1}(\Omega)$ estimate for $r$ can be made independent of $m$ .

7.5 Upper semi-continuity of global attractors

In this section we prove upper semi continuity of global attractors. We first introduce certain concepts pertinent to the upper semicontinuity of attractors.

Definition 7.1 (Uniform dissipativity)

Suppose there is a family of semiflows $\{\{S_{\lambda}(t)\}_{t\geq 0}\}_{\lambda\in\Lambda}$ on a Banach space $H$ , where $\Lambda$ is an open set in a Euclidean space of parameter $\lambda$ , is called uniformly dissipative at $\lambda_{0}\in\Lambda$ , if there is a neighborhood $\mathcal{N}$ of $\lambda_{0}$ in $\Lambda$ and there is a bounded set $\mathcal{B}\subset H$ such that $\mathcal{B}$ is an absorbing set for each semiflow $S_{\lambda}(t)$ , $\lambda\in\mathcal{N}$ , in common.

Definition 7.2 (Upper semicontinuity)

A family of semiflows $\{\{S_{\lambda}(t)\}_{t\geq 0}\}_{\lambda\in\Lambda}$ on a Banach space $H$ , where $\Lambda$ is an open set in a Euclidean space of parameter $\lambda$ , and that there exists a global attractor $\mathcal{A}_{\lambda}$ in $\mathcal{H}$ for each semiflow $\{\{S_{\lambda}(t)\}_{t\geq 0}\}_{\lambda\in\Lambda}$ . If $\lambda_{0}\in\Lambda$ and

dist_{\mathcal{H}}(\mathcal{A}_{\lambda},\mathcal{A}_{\lambda_{0}})\rightarrow 0,\ \mbox{as}\ \lambda\rightarrow\lambda_{0}\ \mbox{in}\ \Lambda,

with respect to the Hausdorff semidistance, then we say that the family of global attractors $\{\mathcal{A}_{\lambda}\}_{\lambda\in\Lambda}$ is upper semicontinuous at $\lambda_{0}$ , or that $\mathcal{A}_{\lambda}$ is robust at $\lambda_{0}$ .

We next recall the following lemma

Lemma 7.3 (Gronwall-Henry Inequality)

Let $\psi(t)$ be a nonnegative function in $L^{\infty}_{loc}[0,T;\mathbb{R})$ and $\zeta(.)$ $\in$ $L^{1}_{loc}[0,T;\mathbb{R})$ , such that the following inequality is satisfied:

\psi(t)\leq\zeta(t)+\mu\int^{t}_{0}(t-s)^{r-1}\psi(s)ds,\ t\in(0,T),

(105)

where $0<T\leq\infty$ , and $\mu,r$ are positive constants. Then it holds that

\psi(t)\leq\zeta(t)+\kappa\int^{t}_{0}\Phi(\kappa(t-s))\psi(s)ds,\ t\in(0,T),

(106)

where $\kappa=(\mu\Gamma(r))^{\frac{1}{r}}$ , $\Gamma(.)$ is the gamma function, and the function $\phi(t)$ is given by

\phi(t)=\sum^{\infty}_{n=1}\frac{nr}{\Gamma(nr+1)}t^{nr-1}=\sum^{\infty}_{n=1}\frac{nr}{\Gamma(nr)}t^{nr-1}.

(107)

We state the following theorem

Theorem 7.4

Consider the reaction diffusion equation described via (2)-(4) where $\Omega$ is of spatial dimension $n=1,2,3$ . There exists a universal constant $K_{\infty}>0$ , independent of the Allee parameter $m$ that bounds uniformly in $L^{\infty}$ the family of global attractors $\mathcal{A}_{m}$ . That is,

\bigcup_{0\leq m\leq 1}\mathcal{A}_{m}\subset B_{L^{\infty}}(0,K_{\infty}).

(108)

Where $B_{L^{\infty}}(0,K_{\infty})$ is the closed ball of radius $K_{\infty}$ in the space $L^{\infty}(\Omega)$ , with $K_{\infty}=3C$ , where the $C$ is from (27).

This follows from the estimates via (27).

We next focus on the uniform dissipativity in $L^{2}(\Omega)$ . The uniform in parameter $m$ $L^{2}(\Omega)$ estimate via (99) enables us to state the following theorem

Theorem 7.5

Consider the reaction diffusion system described via (2)-(4) where $\Omega$ is of spatial dimension $n=1,2,3$ . The family of semiflows $\{\{S_{m}(t)\}_{t\geq 0}\}_{m\geq 0}$ for this system on $H$ , is uniformly dissipative at $m=0$ . Specifically there exists a constant $K_{H}>0$ such that the ball $B_{H}(0,K_{H})$ , which is the closed ball of radius $K_{H}$ in the space $H$ is a common absorbing set for the semiflows $\{\{S_{m}(t)\}_{t\geq 0}\}$ for all $m\in[0,1]$ . Here $K_{H}=1+\frac{2|\Omega|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}}+2C$ , where the $2C$ comes from lemma 2.4.

Next, let us choose $(u_{0},v_{0},r_{0})$ $\in$ $\mathcal{U}$ . We know $\mathcal{U}\subset B_{H}(0,K_{H})$ , the common absorbing ball in $H$ for the semiflow $S_{m}(t)$ , $m\in[0,1]$ , from theorem 7.5. We use the earlier $H^{1}(\Omega)$ and $L^{2}(\Omega)$ estimates on $u,v$ from lemma 2.4, noticing that they do not depend on $m$ , we also use (99) and (7.4) to obtain

	$\displaystyle\|\|S_{m}(t)((u_{0},v_{0},r_{0}))\|\|^{2}_{E}$
	$\displaystyle=\|\|(\nabla u,\nabla v,\nabla r)\|\|^{2}+\|\|(u,v,r)\|\|^{2}$
	$\displaystyle\leq 4C+C_{K}e^{KC_{K}}+1+\frac{2\|\Omega\|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}}.$		(109)

This is true for any $(u_{0},v_{0},r_{0})\in\mathcal{U}$ and $m\in[0,1]$ . Now by invariance $S_{m}(t)\mathcal{A}_{m}=\mathcal{A}_{m}$ , and so $\mathcal{A}_{m}\subset B_{E}(0,K_{E})$ . Where

K_{E}=4C+C_{K}e^{KC_{K}}+1+\frac{2|\Omega|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}}.

(110)

This yields the following theorem

Theorem 7.6

Consider the reaction diffusion equation described via (2)-(4) where $\Omega$ is of spatial dimension $n=1,2,3$ . There exists a universal constant $K_{E}>0$ , independent of the Allee parameter $m$ that bounds uniformly in $E$ the family of global attractors $\mathcal{A}_{m}$ . That is,

\bigcup_{0\leq m\leq 1}\mathcal{A}_{m}\subset B_{E}(0,K_{E}),

(111)

where $B_{E}(0,K_{E})$ is the closed ball of radius $K_{E}$ explicitly given in (110).

We next show that the trajectories $S_{m}(t)U,t\geq 0$ , are uniformly $E$ -bounded, with respect to $m$ . Since the estimate in (110) does not depend on $m$ , we can take a supremum over $m\in[0,1]$ and still achieve the same bound. Thus we have the following result

Theorem 7.7

Consider the reaction diffusion equation described via (2)-(4). There exists a universal constant $K_{E}>0$ , independent of the Allee parameter $m$ such that

\sup_{0\leq m\leq 1}\sup_{t\geq 0}S_{m}(t)\left(\bigcup_{0\leq m\leq 1}\mathcal{A}_{m}\right)\subset B_{E}(0,K_{E}),

(112)

where, $B_{E}(0,K_{E})$ is the closed ball of radius $K_{E}$ in the space $E$ .

We next set $m=0$ in (4) and note that the corresponding estimates in lemma 2.4 hold for $u,v$ , and do not depend on the parameter $m$ . Thus we carry out the analysis as in (7.4) - (110) to still yield

\displaystyle||\nabla r||^{2}_{2}\leq C_{K}e^{KC_{K}},\ \mbox{for}\ t>t_{1}+1.

(113)

Where $C_{K}e^{KC_{K}}$ is independent of the parameter $m$ . Thus we can state the following result

Theorem 7.8

Consider the reaction diffusion equation described via (2)-(4), when $m=0$ , and where $\Omega$ is of spatial dimension $n=1,2,3$ . There exists a $(H,E)$ global attractor $\mathcal{A}$ for the system. This is compact and invariant in $H$ , and it attracts all bounded subsets of $H$ in the $E$ metric.

The proof of the above theorem follows essentially via theorem 3.2.

Note, since $\left(\bigcup_{0\leq m\leq 1}\mathcal{A}_{m}\right)$ is a bounded set in $H$ , we can state the following result ,

Theorem 7.9

Consider the reaction diffusion equation described via (2)-(4), with $m=0$ . The global attractor for this system $\mathcal{A}_{0}$ , attracts the set $\left(\bigcup_{0\leq m\leq 1}\mathcal{A}_{m}\right)$ , with respect to the $E$ -norm. That is

\lim_{t\rightarrow\infty}dist_{E}\left(S_{0}(t)\left(\bigcup_{0\leq m\leq 1}\mathcal{A}_{m}\right),\mathcal{A}_{0}\right).

(114)

We now use the Gronwall-Henry lemma, to show uniform convergence in the parameter $m$ .

Theorem 7.10

Consider the reaction diffusion equation described via (2)-(4), and associated semi-group $S_{m}(t)$ . Also consider (2)-(4), when $m=0$ , that is with associated semi-group $S_{0}(t)$ . For any given $t\geq 0$ , it holds that

\sup_{g_{0}\in\mathcal{U}}||S_{m}(t)g_{0}-S_{0}(t)g_{0}||_{E}\rightarrow 0,\mbox{as}\ m\rightarrow 0^{+}.

(115)

For any given initial data $g_{0}\in\mathcal{U}$ , we note $S_{m}(t)g_{0}$ and $S_{0}(t)g_{0}$ , are both classical solutions, hence mild solutions. Thus we denote

w(t)=S_{m}(t)g_{0}-S_{0}(t)g_{0},t\geq 0,\ \mbox{with}\ w(0)=0.

(116)

Using the form of mild solutions, we can write down the following,

w(t)=\int^{t}_{0}e^{A(t-\sigma)}[f_{0}(S_{m}(\sigma)g_{0})-f_{0}(S_{0}(\sigma)g_{0})]d\sigma+m\int^{t}_{0}e^{A(t-\sigma)}h(S_{m}(\sigma)g_{0})d\sigma,\ t\geq 0.

(117)

Here $e^{At}$ , $t\geq 0$ is the $C_{0}$ semi-group generated by $A:D(A)\rightarrow H$ . Note

f_{m}(u,v,r)=\left(\begin{array}[]{c}u-u^{2}-w_{1}\frac{uv}{u+v}\\ -a_{2}v+w_{2}\frac{uv}{u+v}-w_{3}\left(\frac{vr}{v+r}\right)\\ -\frac{w_{4}}{v+D_{3}}r^{3}+cr^{2}\\ \end{array}\ \right),h(r)=\left(\begin{array}[]{c}0\\ 0\\ -mcr+\frac{mw_{4}}{v+D_{3}}r^{2}\\ \end{array}\right).

(118)

The above is easily seen to be Lipschitz continuous in $E$ , due to the Sobolev embedding of $H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$ . Here $f_{0}$ , which is gotten by setting $m=0$ in the above, is also Lipschitz continuous, thus there is a Lipschitz constant $L_{f_{0}}(K_{E})>0$ , depending only on $K_{E}$ , where $K_{E}$ is given in (110), s.t

||f_{0}(g_{1})-f_{0}(g_{2})||\leq L_{f_{0}}(K_{E})||g_{1}-g_{2}||_{E},

(119)

for any $g_{1},g_{2}\in B_{E}(0,K_{E})$ . From theorems 7.6, 7.7 we have

	$\displaystyle\|\|w(t)\|\|_{E}$
	$\displaystyle\leq\int^{t}_{0}\|\|e^{A(t-\sigma)}\|\|_{\mathcal{L}(H,E)}L_{f_{0}}(K_{E})\|\|S_{m}(\sigma)g_{0}-S_{0}(\sigma)g_{0}\|\|_{E}d\sigma,$
	$\displaystyle+m\int^{t}_{0}\|\|e^{A(t-\sigma)}\|\|_{\mathcal{L}(H,E)}h(S_{m}(\sigma)g_{0})d\sigma,$
	$\displaystyle\leq Cm(K_{E})^{3}\int^{t}_{0}N(t-\sigma)^{-\frac{1}{2}}d\sigma+NL_{f_{0}}(K_{E})\int^{t}_{0}(t-\sigma)^{-\frac{1}{2}}\|\|w(t)\|\|_{E}d\sigma,$
	$\displaystyle\leq Cm(K_{E})^{3}t^{\frac{1}{2}}+NL_{f_{0}}(K_{E})\int^{t}_{0}(t-\sigma)^{-\frac{1}{2}}\|\|w(t)\|\|_{E}d\sigma.$		(120)

This follows via standard decay estimates on the semi-group [22], theorems 7.6, 7.7 and the Sobolev embedding of $H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$ . We now apply the Gronwall-Henry inequality with $\Psi(t)=||w(t)||_{E},\zeta(t)=Cm(K_{E})^{3}t^{\frac{1}{2}},\mu=NL_{f_{0}}(K_{E}),r=\frac{1}{2}$ to obtain

||S_{m}(t)g_{0}-S_{0}(t)g_{0}||_{E}=||w(t)||_{E}\leq m\left(\Psi(t)+k\int^{t}_{0}\Phi(k(t-s))\Psi(s)ds\right),\ t\geq 0.

(121)

Here

\zeta(t)=C(K_{E})^{3}t^{\frac{1}{2}},\ k=NL_{f_{0}}(K_{E})\left(\Gamma\left(\frac{1}{2}\right)\right)^{\frac{1}{2}},\ \Phi(t)=\sum^{\infty}_{n=1}\frac{1}{\Gamma(\frac{n}{2})}t^{\frac{n}{2}-1}.

(122)

Thus we get the uniform convergence as required, and

\sup_{g_{0}\in\mathcal{U}}||S_{m}(t)g_{0}-S_{0}(t)g_{0}||_{E}=||w(t)||_{E}\rightarrow 0,\mbox{as}\ m\rightarrow 0^{+}.

(123)

where the convergence is uniform in the parameter $m$ .

			$\displaystyle\|\|u_{t}\|\|^{2}_{2}+d_{1}\|\|\Delta u\|\|^{2}_{2}+\frac{d}{dt}\|\|\nabla u\|\|^{2}_{2},$
		$\displaystyle=$	$\displaystyle\left\\|\left(u-u^{2}-w_{1}\left(\frac{uv}{u+v}\right)\right)\right\\|_{2}^{2},$
		$\displaystyle\leq$	$\displaystyle C\left(\|\|u\|\|^{2}_{2}+(w_{1})^{2}\|\|u\|\|^{2}_{2}+\|\|u\|\|^{4}_{4}\right)\leq C.$

	$\displaystyle\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2}+\int_{\partial\Omega}\Delta u\nabla u_{t}\cdot\textbf{n}dS,$
	$\displaystyle=\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2}+\int_{\partial\Omega}\Delta u\frac{\partial}{\partial t}(\nabla u\cdot\textbf{n})dS,$
	$\displaystyle=\|\|\nabla u_{t}\|\|^{2}_{2}+d_{1}\|\|\nabla(\Delta u)\|\|^{2}_{2}+\frac{d}{dt}\|\|\Delta u\|\|^{2}_{2}.$		(21)

$\displaystyle\|\|r_{t}\|\|^{2}_{2}+d_{3}\|\|\Delta r\|\|^{2}_{2}+\frac{d}{dt}\|\|\nabla r\|\|^{2}_{2}$	$\displaystyle=$	$\displaystyle\left\\|\left(\left(c+\frac{mw_{4}}{v+D_{3}}\right)r^{2}-\frac{w_{4}}{v+D_{3}}r^{3}-mcr\right)\right\\|_{2}^{2},$
	$\displaystyle\leq$	$\displaystyle C\left(\|\|r\|\|^{2}_{2}+\|\|r\|\|^{4}_{4}\right)+\|\|r\|\|^{6}_{6},$
	$\displaystyle\leq$	$\displaystyle C\left(\|\|\nabla r\|\|^{2}_{2}\right)^{3}.$

	$\displaystyle\|\|S_{m}(t)((u_{0},v_{0},r_{0}))\|\|^{2}_{E}$
	$\displaystyle=\|\|(\nabla u,\nabla v,\nabla r)\|\|^{2}+\|\|(u,v,r)\|\|^{2}$
	$\displaystyle\leq 4C+C_{K}e^{KC_{K}}+1+\frac{2\|\Omega\|^{\frac{1}{4}}\left(c+\frac{w_{4}}{D_{3}}\right)\left(\frac{w_{4}}{D_{3}}\right)^{\frac{1}{4}}(C+D_{3})}{Cw_{4}}.$		(109)

	$\displaystyle\|\|w(t)\|\|_{E}$
	$\displaystyle\leq\int^{t}_{0}\|\|e^{A(t-\sigma)}\|\|_{\mathcal{L}(H,E)}L_{f_{0}}(K_{E})\|\|S_{m}(\sigma)g_{0}-S_{0}(\sigma)g_{0}\|\|_{E}d\sigma,$
	$\displaystyle+m\int^{t}_{0}\|\|e^{A(t-\sigma)}\|\|_{\mathcal{L}(H,E)}h(S_{m}(\sigma)g_{0})d\sigma,$
	$\displaystyle\leq Cm(K_{E})^{3}\int^{t}_{0}N(t-\sigma)^{-\frac{1}{2}}d\sigma+NL_{f_{0}}(K_{E})\int^{t}_{0}(t-\sigma)^{-\frac{1}{2}}\|\|w(t)\|\|_{E}d\sigma,$
	$\displaystyle\leq Cm(K_{E})^{3}t^{\frac{1}{2}}+NL_{f_{0}}(K_{E})\int^{t}_{0}(t-\sigma)^{-\frac{1}{2}}\|\|w(t)\|\|_{E}d\sigma.$		(120)

Long time dynamics of a three-species food chain model with Allee effect in the top predator

Abstract

1 Introduction

2 Preliminary Estimates

2.1 Preliminaries

Lemma 2.1 (Uniform Gronwall Lemma)

2.2 Uniform L2​(Ω)L^{2}(\Omega) and H1​(Ω)H^{1}(\Omega) estimates

Proposition 2.2

Theorem 2.3

Lemma 2.4

Remark 2.5

Remark 2.6

Remark 2.7

2.3 Uniform in the parameter mm L2​(Ω)L^{2}(\Omega) and H1​(Ω)H^{1}(\Omega) estimates for rr

2.4 Uniform H2​(Ω)H^{2}(\Omega) estimates

Remark 2.8

3 Long Time Dynamics

3.1 Existence of global attractor

Theorem 3.1

Theorem 3.2

3.2 Upper-semicontinuity of Global attractor with respect to the Allee parameter mm

Theorem 3.3

Proof 3.4

3.3 Computing the decay rate with respect to the Allee threshold parameter mm

4 Overexploitation Phenomenon

Theorem 4.1 (Middle predator and Allee effect mediated overexploitation)

Remark 4.2

Theorem 4.3 (Top predator mediated overexploitation avoidance)

Remark 4.4

Lemma 4.5 (Persistence in top predator)

Remark 4.6

Remark 4.7

5 Turing Instability

5.1 Steady States

Remark 5.1

Remark 5.2

5.2 Turing conditions

Theorem 5.3

5.3 Effect of Allee threshold parameter mm on Turing Instability

6 Discussion

References

7 Appendix

7.1 Nondimensionalisation

7.2 Global existence

7.3 Apriori estimates

7.4 Apriori estimates independent of the parameter mm

7.5 Upper semi-continuity of global attractors

Definition 7.1 (Uniform dissipativity)

Definition 7.2 (Upper semicontinuity)

Lemma 7.3 (Gronwall-Henry Inequality)

Theorem 7.4

Theorem 7.5

Theorem 7.6

Theorem 7.7

Theorem 7.8

Theorem 7.9

Theorem 7.10

2.2 Uniform $L^{2}(\Omega)$ and $H^{1}(\Omega)$ estimates

2.3 Uniform in the parameter $m$ $L^{2}(\Omega)$ and $H^{1}(\Omega)$ estimates for $r$

2.4 Uniform $H^{2}(\Omega)$ estimates

3.2 Upper-semicontinuity of Global attractor with respect to the Allee parameter $m$

3.3 Computing the decay rate with respect to the Allee threshold parameter $m$

5.3 Effect of Allee threshold parameter $m$ on Turing Instability

7.4 Apriori estimates independent of the parameter $m$