Finite Time Blowup in a Realistic Food-Chain Model

Abstract

We investigate a realistic three-species food-chain model, with generalist top predator. The model based on a modified version of the Leslie-Gower scheme incorporates mutual interference in all the three populations and generalizes several other known models in the ecological literature. We show that the model exhibits finite time blowup in certain parameter range and for large enough initial data. This result implies that finite time blowup is possible in a large class of such three-species food-chain models. We propose a modification to the model and prove that the modified model has globally existing classical solutions, as well as a global attractor. We reconstruct the attractor using nonlinear time series analysis and show that it pssesses rich dynamics, including chaos in certain parameter regime, whilst avoiding blowup in any parameter regime. We also provide estimates on its fractal dimension as well as provide numerical simulations to visualise the spatiotemporal chaos.

1. Introduction

Interaction networks in natural ecosystems can be visualized as consisting of simple units known as food chains or food webs that are made up of a number of species linked by trophic interaction [1]. A food chain model essentially comprises of the predator-prey relationship between interacting species in a given ecosystem [2]. In their seminal work [3], Hastings and Powell for the first time demonstrated that the evolution of species participating in a tritrophic relationship might be chaotic [3]. This led to a great deal of research activity in analyzing the dynamical behavior of food chain models. Upadhyay and Rai [4] provided a new example of a chaotic population system in a simple three-species food chain with Holling type II functional response. This model is different from the Hastings and Powell model, in that it considers a generalist top predator, one that can switch its food source, in the absence of its favorite prey. Letellier and Aziz-Alaoui [5] and Aziz-Alaoui [6] revisited the Upadhyay and Rai model and found that the chaotic dynamics is observed via sequences of period-doubling bifurcation of limit cycles which however suddenly break down and reverse giving rise to a sequence of period-halving bifurcation leading to limit cycles. Upadhyay in [7] next proposed a three-species food-chain model, by incorporating mutual interference, in the original model [4], thus generalizing the models in [3, 4, 8]. Parshad and Upadhyay [9] considered the diffusive form of the model proposed in [7]. Under certain restrictions on the parameter space, they proved global existence of solutions as well as the existence of a finite dimensional global attractor, for the diffusive model. The original three-species model for a generalist top predator [4] and its variants are very rich dynamically and have led to a number of works in the literature [4–9]. They continue to be recently investigated, for rich dynamical behavior, such as turing instability [10]. Despite these rich dynamics, the Parshad and Upadhyay model possesses a drawback, in that there is the possibility of finite time blowup in the model. Recall in general we say that the ‖⋅‖_𝑋 norm of a solution blows up in finite time if lim_{𝑡→𝑇∗<∞}‖𝑢‖_𝑋=+∞, (1) where 𝑋 is a suitable function space where one looks for the solution 𝑢, to the PDE or ODE at hand. Blow-up phenomenon in physical and biological problems has a rich history and has been studied in the context of gas dynamics, chemotaxis, cooperative, and competitive lotka-Volterra type systems [11–20]. Understanding blowup in a model can help to determine how robust a system is, as well as providing insight on the characteristics of the corresponding physical or biological phenomena.

As the Parshad and Upadhyay model generalises many of the aforementioned three-species models, our result actually implies finite time blowup in all of those models as well. This, however, has not been reported in the literature. It is well documented in the literature [5] that in the absence of the middle predator the top predator becomes extinct if 𝑐𝐷₃<𝑤₃. This essentially means that if the middle predator is absent the top predator has to switch its food source. However, it is probably not able to predate efficiently enough, and thus the population cannot be sustained, unless the rate of sexual reproduction is above a certain threshold. If this is not so, that is, 𝑐𝐷₃<𝑤₃, then the species goes extinct. Also, it grows unboundedly if the inequality reverses. We remark that the situation can be much worse than unbounded growth and that finite time blowup is actually possible.

In the current paper our primary contributions are the following. (1)We show that the diffusive model of Parshad and Upadhyay [9] can blow up in finite time, for certain parameter regime and large enough initial data, via Theorem 2. This immediately shows the blowup of the diffusive form of the model of Upadhyay and Rai [7]. The temporal model proposed in [4] is also shown to blow up in finite time, in certain parameter regime. (2)We propose a modification to the diffusive model and show global existence of classical solutions and global attractor, for this model, via Theorems 6 and 14. The modified model is well posed and does not blow up in finite time, for any data or in any parameter regime. (3)We reconstruct the attractor of the modified model, via nonlinear time series analysis, and show that it also exhibits rich dynamics, including chaos in certain parameter regime, just like the old model, whilst avoiding finite time blowup in any parameter regime. Numerical simulations are also provided to show the blowup in the classical models, as well as spatiotemporal chaos, in the modified model.

This paper is organized as follows . In Section 2, we recount the diffusive three-species food-chain model [9], that generalises several other known models in the literature. Finite time blowup results are summarized in Section 3. We propose a modification to the general model in Section 4 and show global existence of classical solutions and global attractor for this model. Various estimates made here concerning the middle predator and prey species are very similar to [21]. We reconstruct the chaotic attractor of the modified model, using nonlinear time series analysis, and estimate its fractal dimension. This, as well as some numerical simulations, is presented in Section 5. Some important conclusions are discussed in Section 6.

2. Food Chain Model System

In this section, we recount the model from [9]. Consider a situation where a prey population 𝑥₁ is predated by individuals of population 𝑥₂. The population 𝑥₂ in turn serves as a favourite food for individuals of population 𝑥₃. This interaction is represented by the following system of a simple prey-specialist predator-generalist predator interaction with the inclusion of spatial spread:𝜕𝑥₁𝜕𝑡=𝑑₁Δ𝑥₁+𝑎₁𝑥₁−𝑏₁𝑥²₁−𝑤₀𝑥^𝑚2₂𝑥₁𝑥₁+𝐷₀^𝑚1=𝑔₁𝑥₁,𝑥₂,𝑥₃+𝑑₁Δ𝑥₁,𝜕𝑥₂𝜕𝑡=𝑑₂Δ𝑥₂−𝑎₂𝑥₂+𝑤₁𝑥^𝑚2₂𝑥₁𝑥₁+𝐷₁^𝑚1−𝑤₂𝑥^𝑚2₃𝑥₂𝑥₂+𝐷₂^𝑚2=𝑔₂𝑥₁,𝑥₂,𝑥₃+𝑑₂Δ𝑥₂,𝜕𝑥₃𝜕𝑡=𝑑₃Δ𝑥₃+𝑐𝑥^𝑚3₃−𝑤₃𝑓₁𝑥₂𝑥^𝑚3₃=𝑔₃𝑥₁,𝑥₂,𝑥₃+𝑑₃Δ𝑥₃, (2)(3)(4) where 𝑎₁, 𝑎₂, 𝑏₁, 𝑐, 𝑤₀, 𝑤₁, 𝑤₂, 𝑤₃, 𝐷₀, 𝐷₁, 𝐷₂, 𝐷₃≥0, 𝑓₁(𝑥₂)=1/(𝑥₂+𝐷₃)∈𝐶(𝑅₊). The parameters 𝑚_𝑖 for 𝑖=1,2,3 are mutual interference parameters that model the intraspecific competition among predators when hunting for prey [22–27]. The problem is posed on Ω⊂ℝ².  Ω is bounded, and 𝜕Ω is assumed to be smooth. We consider the Dirichlet boundary conditions 𝑥₁=𝑥₂=𝑥₃=0on𝜕Ω. (5)

We also assume suitable initial conditions. In this model, prey population of size 𝑥₁ serves as the only food for the specialist predator population of size 𝑥₂. This predator population, in turn, serves as a favorite food for the generalist predator population of size 𝑥₃. The equations for rate of change of population size for prey and specialist predator have been written following the Volterra scheme; that is, predator population dies out exponentially in the absence of its lone prey. The interaction between this predator 𝑥₂ and the generalist predator 𝑥₃ is modeled by the modified version of the Leslie-Gower scheme where the loss in a predator population is proportional to the reciprocal of per capita availability of its most favorite food. 𝑎₁ is the intrinsic growth rate of the prey population 𝑥₁, 𝑎₂ is the intrinsic death rate of the predator population 𝑥₂ in the absence of the only food 𝑥₁, 𝑐 measures the rate of self-reproduction of generalist predator 𝑥₃, and 𝑤₀, 𝑤₁, 𝑤₂, and 𝑤₃ are the maximum values which per capita growth rate can attain. 𝑏₁ measures the strength of intraspecific competition among the individuals of the prey species 𝑥₁. 𝐷₀ and 𝐷₁ quantify the extent to which environment provides protection to the prey 𝑥₁ and may be thought of as a refuge or a measure of the effectiveness of the prey in evading a predator’s attack. 𝐷₂ is the value of 𝑥₂ at which per capita removal rate of 𝑥₂ becomes 𝑤₂/2. 𝐷₃ represents the residual loss in 𝑥₃ population due to severe scarcity of its favourite food 𝑥₂. For 𝑚₁=𝑚₂=1 the coefficient 𝑤₀/(𝑥₁+𝐷₀) of the third term on the right-hand side of (2) is obtained by considering the probable effect of the density of the prey’s population on predator’s attack rate. If this coefficient is multiplied by 𝑥₁ (the prey population at any instant of time), it gives the attack rate on the prey per predator. Denote 𝑝(𝑥₁)=(𝑤₀𝑥₁)/(𝑥₁+𝐷₀), when 𝑥₁→∞, 𝑝(𝑥₁)→𝑤₀, which is the maximum that it can reach. The third term (𝑤₂𝑥₂𝑥₃)/(𝑥₂+𝐷₂) on the right-hand side of (3) represents the per capita functional response of the invertebrate predator 𝑥₃ and was first introduced by Holling [28] in the ecological literature. He carried out specially designed experiments to determine the nature of these functional responses. The ecological role of per capita functional response has been well described by May [29]. The interaction terms appearing in the rate equation restore to some extent the symmetry which characterizes the Lotka-Volterra model. The generalist predator 𝑥₃, in (4), is sexually reproducing species. We assume that males and females are equal in number and every individual has got equal opportunity to meet an individual of opposite sex. The first term of (4) represents growth rate of the sexually reproducing species in well-mixed conditions. 𝑤₃ measures the limitation on the growth of the generalist predator 𝑥₃ by its dependence on per capita availability of its most favorite prey 𝑥₂.

The study of ecological problems render’s one to regard communities of individuals as subpopulations affecting the survival of the individuals of other populations. One aspect of the dynamics of community interactions would be mutual interference among the interacting subpopulations, which in our model system are represented by the parameters 𝑚_𝑖,𝑖=1,2,3. It has been shown in [24–26] that mutual interference is generally a “stabilizing” process. The parameters 𝑑₁, 𝑑₂, and 𝑑₃ are diffusion coefficients of the populations.

3. Finite Time Blowup

In this section, we show that (2)–(4) blow up in finite time, for certain parameter regime and large enough initial data. In particular, the equation for 𝑥₃ as described via (4) blows up. To this end we employ the first eigenvalue method for blow up,  [30]. First recall Jensen’s inequality with a probability measure.

Lemma 1. Let 𝑝(𝑥) be a probability measure on ℝ, and let 𝑢 be a nonnegative function. Then the following holds whenever 𝑓 is a convex function: 𝑓_𝑅𝑢(𝑥)𝑝(𝑥)≤_𝑅𝑓(𝑢(𝑥))𝑝(𝑥). (6)

Note, this is easily generalised to higher dimensions. We next present our blow-up result.

Theorem 2. Consider the three-species food-chain model as described via (2)–(4). For 𝑐>(𝑤₃/𝐷₃) and initial data 𝑥₃₀(𝑥) large enough, that is _Ω𝑥₃₀(𝑥)𝑊₁(𝑥)𝑑𝑥>𝜆₁𝑑₃1𝑐−𝑤₃/𝐷₃^{1/(𝑚3−1)}, (7) there exists a finite time 𝑇^∗>0, such, that for any 𝑚₃>1, lim_{𝑡→𝑇∗Ω}𝑥₃(𝑥,𝑡)𝑊₁(𝑥)𝑑𝑥=∞. (8)
Here, 𝑊₁(𝑥) is a function such that −Δ𝑊₁(𝑥)=𝜆₁𝑊₁(𝑥) on Ω, 𝑊₁(𝑥)=0 on 𝜕Ω. Furthermore, 𝑊₁(𝑥)>0 and ∫_Ω𝑊₁(𝑥)𝑑𝑥=1.

Proof. Consider the equation for 𝑥₃ rewritten as 𝜕𝑥₃𝜕𝑡+𝑤₃𝑥^𝑚3₃𝑥₂+𝐷₃=𝑑₃Δ𝑥₃+𝑐𝑥^𝑚3₃. (9) Since 𝑥₂ is assumed positive, we obtain 𝜕𝑥₃𝜕𝑡+𝑤₃𝑥^𝑚3₃𝐷₃≥𝑑₃Δ𝑥₃+𝑐𝑥^𝑚3₃. (10) Thus, 𝜕𝑥₃𝜕𝑡≥𝑑₃Δ𝑥₃+𝑤𝑐−³𝐷₃𝑥^𝑚3₃. (11)
Now we multiply (11) by 𝑊₁, consider 𝐹(𝑡)=∫_Ω𝑥₃(𝑥,𝑡)𝑊₁(𝑥)𝑑𝑥, and integrate by parts to obtain 𝐹^′(𝑡)=_Ω𝜕𝑥₃𝑊𝜕𝑡¹(𝑥)≥𝑑𝑥_Ω𝑑₃Δ𝑥₃+𝑤𝑐−³𝐷₃𝑥^𝑚3₃𝑊₁(𝑥)𝑑𝑥=−𝑑₃𝜆₁𝐹(𝑡)+𝑤𝑐−³𝐷_3Ω𝑥^𝑚3₃𝑊₁(𝑥)𝑑𝑥. (12)
Via assumptions on the eigenfunction 𝑊₁, it is seen to be a probability measure. That is, we can define a probability measure 𝑝(𝑥)=𝑊₁ on Ω and 𝑝(𝑥)=0onℝ²/Ω. Now we consider the convex function 𝑓(𝑥₃)=𝑥^𝑚3₃ (note that any 𝑚₃>1 works as the 𝑥₃ in question is always positive). Thus, the application of Jensen’s inequality yields 𝑤𝑐−³𝐷_3Ω𝑥^𝑚3₃𝑊₁(𝑥)≥𝑤𝑑𝑥𝑐−³𝐷_3Ω𝑥₃(𝑥,𝑡)𝑊₁(𝑥)𝑑𝑥^𝑚3=𝑤𝑐−³𝐷₃(𝐹(𝑡))^𝑚3. (13)
Thus, substituting the previous equation into (12) yields 𝐹^′(𝑡)≥−𝑑₃𝜆₁𝐹(𝑡)+𝑤𝑐−³𝐷₃(𝐹(𝑡))^𝑚3. (14)
Simple integration in time of the previous equation yields that 𝐹(𝑡) blows up in finite time; that is lim_{𝑡→𝑇∗}𝐹(𝑡)=lim_{𝑡→𝑇∗Ω}𝑥₃(𝑥,𝑡)𝑊₁(𝑥)𝑑𝑥=∞, (15) where 𝑇^∗∗=1𝑑₃𝜆₁𝑤log𝑐−₃/𝐷₃1−𝑚₃𝑑₃𝜆₁−(𝐹(0))^1−𝑚31−𝑚₃ (16) and 𝑇^∗≤𝑇^∗∗.
After some algebra, we obtain that for blowup the data has to satisfy 𝐹(0)=_Ω𝑥₃₀(𝑥)𝑊₁(𝑥)𝑑𝑥>𝜆₁𝑑₃1𝑐−𝑤₃/𝐷₃^{1/(𝑚3−1)}. (17) This proves the theorem.

Remark 3. Note that the blowup of 𝐹(𝑡) immediately implies that the 𝐿^∞ norm of 𝑟 blows up as follows: lim_{𝑡→𝑇∗}𝐹(𝑡)≤lim_{𝑡→𝑇∗}‖𝑟(𝑥,𝑡)‖_∞Ω𝑊₁(𝑥)𝑑𝑥=lim_{𝑡→𝑇∗}‖𝑟(𝑥,𝑡)‖_∞=∞. (18)

Remark 4. We also show blow-up in the Upadhyay-Rai model proposed in [4], that is, the temporal case with 𝑚₁=1,𝑚₂=1, and𝑚₃=2. Also, for the model in [4], the blow-up time is given by 𝑇^∗∗=1𝑑₃𝜆₁log𝑐−𝑤₃/𝐷₃𝐹(0)𝑐−𝑤₃/𝐷₃𝐹(0)−𝜆₁𝑑₃. (19)

Remark 5. We note that finite time blowup is also possible in the same weighted 𝐿¹(Ω) norm, for Neumann boundary conditions. In this case, the eigenfunction 𝑊₁ will have to be differently defined, that is, satisfying Neumann boundary conditions itself.

4. Modification to the Old Model

The aim of the current section is to introduce an improvement to (2)–(4), which is well posed. The key issue in amending the old model, (2)–(4), is to correct (4). Upadhyay and Rai introduced sexual reproduction in (4), by considering 𝑐𝑥²₃. Under this assumption, the number of males and females in the population is assumed to be equal, [5]. Thus, mating takes place as 𝑐×(1/2)𝑥₃×(1/2)𝑥₃, and consequently the population grows at rate 𝑐. We point out that this assumption may fail in various structured populations, where due to specific demographics, the numbers of females and males are rarely equal. We modify this by considering the mating term to be 𝑐𝑥₃^2−𝜖=𝑐𝑥₃×𝑥₃^(1−𝜖), where we restrict 0<𝜖<1. This modification may be interpreted as modeling a population where the number of males is different from the number of females (even if very slightly so). If the numbers of males and females are in fact equal or close to, one can choose 0<𝜖≪1. Next, the mutual interference parameter is incorporated, by raising the reaction terms to the 𝑚₃ power. The primary role of the previously stated modification is that now the reaction term in (4) is 𝑐𝑥^𝑚3₃−𝑤₃𝑓₁(𝑥₂)𝑥^𝑚4₃, and because we enforce 𝑚₃<𝑚₄, standard parabolic theory [31] guarantees global existence. Based on the above, we introduce the following three-species food-chain model, with a generalist top predator, 𝜕𝑥₁𝜕𝑡=𝑑₁Δ𝑥₁+𝑎₁𝑥₁−𝑏₁𝑥²₁−𝑤₀𝑥^𝑚2₂𝑥₁𝑥₁+𝐷₀^𝑚1,𝜕𝑥₂𝜕𝑡=𝑑₂Δ𝑥₂−𝑎₂𝑥₂+𝑤₁𝑥^𝑚2₂𝑥₁𝑥₁+𝐷₀^𝑚1−𝑤₂𝑥₂𝑥₂+𝐷₂^𝑚2𝑥^𝑚2₃,𝜕𝑥₃𝜕𝑡=𝑑₃Δ𝑥₃+𝑐𝑥^𝑚3₃−𝑤₃𝑓₁𝑥₂𝑥^𝑚4₃. (20)(21)(22)

The problem is posed on Ω⊂ℝ². Ω is bounded, and 𝜕Ω is assumed to be smooth. We consider the Dirichlet boundary conditions 𝑥₁=𝑥₂=𝑥₃=0on𝜕Ω. (23)

Note, the proofs of global existence to follow will work the with Neumann boundary conditions as well. The forms of (2)-(3) are kept essentially the same; however, we assume that 𝐷₀=𝐷₁, as commonly done in numerical simulations, [5]. We also assume 0<𝑚₂≤1. This facilitates the estimates to follow. Restricting 𝑚₃ to the range 1<𝑚₃<𝑚₄ has distinct advantages over the old model (2)–(4). Firstly, finite time blowup is prevented, and this is independent of the sign of 𝑐−𝑤₃/𝐷₃. Secondly, in the absence of 𝑥₂, if 𝑐−𝑤₃/𝐷₃<0, the dynamics of (2)–(4) is quite boring, and 𝑥₃ decays to 0, [5, 6]. However, if 1<𝑚₃<𝑚₄, even if 𝑥₂=0 and 𝑐−𝑤₃/𝐷₃<0, 𝑥₃ can still persist. Thus, this formulation plausibly models a true generalist 𝑥₃. One that does not go extinct, even if its favorite food source does.

4.1. Preliminaries

The system (20)–(22) is an example of a good semilinear parabolic system, and the standard theory for parabolic systems [31] ensures the existence of a nonnegative local solution. Next, comparison arguments can be used to derive global existence for initial data in 𝐿^∞(Ω). We assume that 𝑥₁, 𝑥₂, and 𝑥₃ are nonnegative because they are densities. This assumption ensures that 𝑥₁+𝐷₀>0,  𝑥₁+𝐷₁>0,  𝑥₂+𝐷₂>0,  and 𝑥₂+𝐷₃>0. The nonlinearities are thus well defined. Standard parabolic theory guarantees the uniform boundedness of 𝑥₃. We state the following theorem.

Theorem 6. Consider the three-species food-chain model described by (20)–(22). For any positive initial data (𝑥₁₀,𝑥₂₀,𝑥₃₀)∈𝐿²(Ω) and 1<𝑚₃<𝑚₄, 0<𝑚₂≤1, there exists a global classical solution (𝑥₁,𝑥₂,𝑥₃) to the system.

Remark 7. In this section, we demonstrate a stronger version of global existence. That is, we derive uniform 𝐿^∞(Ω) bounds on the solution with the initial data only in 𝐿²(Ω).

The opposite signed nonlinearities here cause problems in making direct 𝐻¹ estimates and proving the existence of bounded absorbing sets, for the 𝑥₂ variable. To circumvent this difficulty, we resort to grouping estimates via appropriate addition of the equations at hand. This technique is very similar to what we adopted in [21]. This condition plays a key role in the dissipation process. If it is satisfied, there is often a global attractor. We will also show the existence of a global attractor for the reaction diffusion system (20)–(22).

Remark 8. Note, in [9] in order to proceed, we need to assume 𝑐≤𝑤₃/(ℎ+𝑣). Thus, our results in [9] are only true in the parameter regime, such that 𝑣<(𝑤₃/𝑐)−𝐷₃. Furthermore in [9], strong 𝐿^∞ boundedness requirements on the solutions are assumed. Here, we enforce no such requirement, and the methods of analysis are quite different from [9].

Definition 9. Consider a semigroup, 𝑆(𝑡), acting on a reflexive Banach space, 𝐻. Then, the global attractor, 𝒜⊂𝐻, for this semigroup is an object that has the following properties: (i)𝒜 is compact in 𝐻,(ii)𝒜 is invariant; that is, 𝑆(𝑡)𝒜=𝒜,𝑡≥0,(iii)if 𝐵 is bounded in 𝑀, then dist_𝑀(𝑆(𝑡)𝐵,𝒜)→0, 𝑡→∞.
Next, various preliminaries are presented, detailing the phase spaces of interest and recalling certain standard theories. Let us define our phase spaces of interest as 𝐻=𝐿²(Ω)³,𝑉=𝐻¹(Ω)³. (24) To prove the existence of a global attractor, we are required to show [32](i)the existence of a bounded absorbing set in the phase space; (ii)the asymptotic compactness property of the semigroup in question. These are defined next.

Definition 10 (bounded absorbing set). A bounded set, ℬ, in a reflexive Banach space, 𝐻, is called a bounded absorbing set if, for each bounded subset 𝑈 of 𝐻, there is a time, 𝑇=𝑇(𝑈), such that 𝑆(𝑡)𝑈⊂ℬ for all 𝑡>𝑇. The number 𝑇=𝑇(𝑈) is referred to as the compactification time for 𝑆(𝑡)𝑈. This is essentially the time after which the semigroup compactifies.

Definition 11 (asymptotic compactness). The semigroup {𝑆(𝑡)}_𝑡≥0∶𝐻→𝐻 associated with a dynamical system is said to be asymptotically compact in 𝐻 if, for any {𝑢_𝑛(0)}^∞_𝑛=1 bounded in 𝐻 and a sequence of times {𝑡_𝑛→∞},   𝑆(𝑡_𝑛)𝑢_𝑛(0) possesses a convergent subsequence in 𝐻.

4.2. Global Existence and Existence of Absorbing Sets in 𝐿²(Ω)

We begin by multiplying (20) by 𝑥₁ and integrating by parts over Ω. This yields 12𝑑𝑥𝑑𝑡¹²²=−𝑑₁∇𝑥₁²₂+𝑎₁𝑥₁²₂−𝑏₁𝑥₁³₃−𝑤_0Ω𝑥^𝑚1₁⁺¹𝑥^𝑚2₂𝑥₁+𝐷₀^𝑚1𝑑𝑥. (25)

We then use Hölder’s inequality followed by Young’s inequality to obtain 12𝑑𝑥𝑑𝑡¹²²≤−𝑑₁∇𝑥₁²₂+𝑏₁𝑥₁³₃+𝐶−𝑏₁𝑥₁³₃−𝑤_0Ω𝑥^𝑚1₁⁺¹𝑥^𝑚2₂𝑥₁+𝐷₀^𝑚1𝑑𝑥. (26) Next, using the compact Sobolev embedding, 𝐻¹(Ω)↪𝐿²(Ω), and the positivity of 𝑥₁ and 𝑥₂, we find that 𝑑𝑥𝑑𝑡¹²²+𝐶₁𝑥₁²₂≤𝐶, (27) where 𝐶 depends explicitly on 𝑏₁ and 𝑎₁.

Thus, application of Gronwall’s Lemma gives us the following estimate: 𝑥₁²₂≤𝑒^−𝐶1𝑡𝑥₁(0)²₂+𝐶𝐶₁. (28)

Therefore, there exists time 𝑡₁ given explicitly by 𝑡₁𝑥=max0,ln₁(0)²₂𝐶₁, (29) such that, for all 𝑡≥𝑡₁, the following estimate holds uniformly: 𝑥₁²₂𝐶≤1+²𝐶₁. (30)

We now make a local in time estimate for ∇𝑥₁. By integrating (26) in the time interval [𝑡₁,𝑡₁+1], we obtain ^𝑡1_𝑡⁺¹₁∇𝑥₁²₂𝑑𝑡≤𝑥₁𝑡₁²₂+^𝑡1_𝑡⁺¹₁𝐶₂𝑑𝑡≤𝐶. (31) Thus, via a mean value theorem for integrals, there exists a time 𝑡₂∈[𝑡₁,𝑡₁+1] such that the following estimate holds: ∇𝑥₁𝑡₂²₂≤𝐶. (32)

We now move on to showing the existence of an absorbing set for 𝑥₂ in 𝐿²(Ω). By multiplying (20) by 𝑤₁ and (21) by 𝑤₀, adding the two together, and setting 𝑤=𝑤₁𝑥₁+𝑤₀𝑥₂, we obtain 𝜕𝑤𝜕𝑡=𝑑₂Δ𝑤+𝑑₁−𝑑₂𝑤₁Δ𝑥₁+𝑤₁𝑎₁𝑥₁−𝑤₁𝑏₁𝑥²₁−𝑤₀𝑎₂𝑥₂−𝑤₀𝑤₂𝑥₂𝑥₁+𝐷₂^𝑚2𝑥^𝑚2₃. (33) By multiplying (33) by 𝑤, integrating by parts over Ω, and using Hölder’s inequality followed by Young’s inequality, keeping in mind the positivity of the solutions, we find that 12𝑑‖𝑤‖𝑑𝑡²²+𝑑₂‖∇𝑤‖²₂≤𝑤₁𝑎_1Ω𝑥₁𝑤𝑑𝑥+_Ω𝑑₂−𝑑₁𝑤₁∇𝑥₁‖𝑤‖⋅∇𝑤𝑑𝑥≤𝐶²₂+𝐶₂𝑥₁²₂+𝐶₁∇𝑥₁²₂+𝑑₂2‖∇𝑤‖²₂. (34) Thus, by using the compact Sobolev embedding, 𝐻¹(Ω)↪𝐿²(Ω), we obtain 𝑑‖𝑤‖𝑑𝑡²²+𝐶‖𝑤‖²₂≤𝐶₁∇𝑥₁²₂+𝐶₂𝑥₁²₂. (35) Next, we multiply the previous equation by 𝑒^𝐶𝑡 and integrate in time from 0 to 𝑡 to find ‖𝑤‖²₂≤𝑒^−𝐶𝑡‖𝑤(0)‖²₂+^𝑡₀𝑒^{−(𝑡−𝜏)}∇𝑥₁²₂𝑑𝜏+𝐶₁𝑥₁₀,𝑡, (36) where 𝐶₁(𝑥₁₀,𝑡)=𝑒^−𝐶𝑡∫^𝑡₀𝑒^−𝐶1𝜏‖𝑥₁(0)‖²₂𝑑𝜏+(𝐶/𝐶₁). This follows via (30).

We now make the following estimate via time integration of (26), after multiplying through by 𝑒^𝐶𝑡: ^𝑡₀𝑒^𝐶𝜏∇𝑥₁²₂𝑑𝜏≤𝐶𝑒^𝐶𝑡−12^𝑡₀𝑒^𝜏𝑑𝑥𝑑𝜏¹²²=𝐶𝑒^𝐶𝑡−12𝑒^𝐶𝑡𝑥₁(𝑡)²₂−𝑥₁₀²₂−^𝑡₀𝑒^𝐶𝜏𝑥₁(𝜏)²₂𝑑𝜏≤𝐶𝑒^𝐶𝑡+𝑥₁₀²₂+^𝑡₀𝑒^𝐶𝜏𝑥₁(𝜏)²₂𝑑𝜏≤𝐶𝑒^𝐶𝑡+𝑥₁₀²₂. (37) This follows via (30).

We now substitute the above into (36) and use the form of 𝐶₁(𝑥₁₀,𝑡) to obtain ‖𝑤‖²₂≤𝑒^−𝐶𝑡‖𝑤(0)‖²₂+𝐶₁𝑥₁₀,𝑡+𝐶𝑒^−𝐶𝑡𝐶𝑒^𝐶𝑡+𝑥₁₀²₂. (38) Therefore, there exists time 𝑡₃, given explicitly by 𝑡₃𝑤=max0,ln₁𝑥₁(0)²₂+𝑤₀𝑥₁(0)²₂𝐶, (39) such that, for all 𝑡≥𝑡₃, the following estimate holds uniformly: 𝑥₂²₂≤1𝑤₀‖𝑤‖²₂𝐶≤1+𝐶+𝐶₁. (40) Therefore, limsup_𝑡→∞𝑥₂²₂≤1𝑤₀limsup_𝑡→∞‖𝑤‖²₂≤𝐶, (41) where 𝐶 is a constant that is independent of the time and the initial data.

We next derive local in time estimates for ∇𝑥₂. We proceed by multiplying (21) by 𝑥₂ and integrating by parts over Ω to obtain 12𝑑𝑥𝑑𝑡²²²=−𝑑₂∇𝑥₂²₂−𝑎₂𝑥₂²₂+𝑤_1Ω𝑥₁𝑥₁+𝐷₀^𝑚1𝑥^𝑚2₂⁺¹𝑑𝑥−𝑤_2Ω𝑥^𝑚2₂⁺¹𝑥₂+𝐷₂^𝑚2𝑥^𝑚2₃𝑑𝑥. (42)

We now use Hölder’s inequality, positivity of solutions, the embedding 𝐿²↪𝐿^𝑚2+1 (via assumption on 𝑚₂), and the earlier estimate via (40) to yield 𝑑𝑥𝑑𝑡²²²+𝑑₂∇𝑥₂²₂≤𝑥₁𝑥₁+𝐷₀^𝑚1_∞𝑥₂²₂≤𝐶₁. (43) By integrating this inequality in the time interval [𝑡,𝑡+1], for 𝑡≥𝑡₃, we obtain _𝑡^𝑡+1∇𝑥₂²₂𝑑𝑠≤𝑥₂(𝑡)²₂+^𝑡1_𝑡⁺¹₁𝐶₁|Ω|𝑑𝑠≤𝐶. (44) The existence of an absorbing set for 𝑥₃ in 𝐿²(Ω) follows easily via the form of (22).

Standard estimates imply the existence of a time 𝑡₅ given explicitly by 𝑡₅𝑥=max0,ln₃(0)²₂𝐶, (45) such that, for all 𝑡≥𝑡₅, the following estimates hold uniformly: 𝑥₃²₂𝐶≤1+²𝐶_𝑡𝑡+1,∇𝑥₃²₂𝑑𝑠≤𝐶. (46)

We can now state the following lemma.

Lemma 12. Consider (𝑥₁,𝑥₂,𝑥₃) which are solutions to the diffusive three-species food-chain model described via (20)–(22), with (𝑥₁₀,𝑥₂₀,𝑥₃₀)∈𝐿²(Ω), 1<𝑚₃<𝑚₄, and 0<𝑚₂≤1. There exists a time 𝑡^∗=max(𝑡₁,𝑡₃,𝑡₅) and a constant C independent of time and the initial data and dependent only on the parameters in (20)–(22), such that, for any 𝑡>𝑡^∗, the following uniform estimates hold: 𝑥₁²₂≤𝐶,𝑥₂²₂≤𝐶,𝑥₃²₂≤𝐶. (47)

We now make the 𝐻¹ estimate. We take the inner product of (20) with −Δ𝑥₁ and integrate by parts over Ω. Thus, we obtain 12𝑑𝑑𝑡∇𝑥₁²₂+𝑑₁Δ𝑥₁²₂+2𝑏_1Ω𝑥₁∇𝑥₁²𝑑𝑥=𝑎₁∇𝑥₁²₂+𝑤_0Ω𝑥₁𝑥₁+𝐷₀^𝑚1−Δ𝑥₁𝑥₂^𝑚2𝑑𝑥. (48) We apply the Cauchy-Schwartz inequality, the Cauchy with epsilon, and Young’s inequalities on the last term, along with the embedding 𝐿²↪𝐿^2𝑚2 (via assumption on 𝑚₂) to obtain 12𝑑𝑑𝑡∇𝑥₁²₂+𝑑₁Δ𝑥₁²₂+2𝑏_1Ω𝑥₁∇𝑥₁²𝑑𝑥≤𝑎₁∇𝑥₁²₂+𝑑₁2Δ𝑥₁²₂+𝐶₁𝑥₁𝑥₁+𝐷₀^𝑚1_∞Ω𝑥₂²𝑑𝑥. (49) This yields 12𝑑𝑑𝑡∇𝑥₁²₂≤𝐶∇𝑥₁²₂+𝐶₁𝑥₂²₂. (50)

We now recall the following lemma.

Lemma 13 (the uniform Gronwall lemma). Let 𝛽, 𝜁, and ℎ be nonnegative functions in 𝐿¹_𝑙𝑜𝑐[0,∞;𝑅). Assume that 𝛽 is absolutely continuous on (0,∞) and the following differential inequality is satisfied: 𝑑𝛽𝑑𝑡≤𝜁𝛽+ℎ,for𝑡>0. (51) If there exists a finite time 𝑡₁>0 and some 𝑞>0 such that _𝑡^𝑡+𝑞𝜁(𝜏)𝑑𝜏≤𝐴,_𝑡^𝑡+𝑞𝛽(𝜏)𝑑𝜏≤𝐵,_{𝑡𝑡+𝑞}ℎ(𝜏)𝑑𝜏≤𝐶, (52) for any 𝑡>𝑡₁, where 𝐴,𝐵, and 𝐶 are some positive constants, then 𝛽(𝑡)≤𝐵𝑞+𝐶𝑒^𝐴,forany𝑡>𝑡₁+𝑞. (53)

We apply the previous lemma to (50) with 𝜁(𝜏)=𝐶,𝛽(𝜏)=∇𝑥₁²₂,ℎ(𝜏)=𝐶₁𝑥₂²₂, (54) and 𝑡₁=𝑡^∗, 𝑞=1, and we use the estimates via (31) and Lemma 12 to obtain ∇𝑥₁(𝑡+1)²₂≤𝐶,𝑡>𝑡^∗. (55) The 𝐻¹ estimates on the other components, 𝑥₂ and 𝑥₃, are obtained similarly.

Recall that in ℝ² we have via the Sobolev embedding that 𝐻¹(Ω)↪𝐿^𝑝(Ω),  for all 𝑝, limsup_𝑡→∞𝑥_1𝑝≤𝐶, (56) and similarly, for 𝑥₂ and 𝑥₃, we have limsup_𝑡→∞𝑥_2𝑝≤𝐶,limsup_𝑡→∞𝑥_3𝑝≤𝐶. (57)

Using the previous equation it is easy to show that the reaction terms 𝑔₁, 𝑔₂, and 𝑔₃ are in 𝐿^𝑝, for 𝑝>1. In general, 𝑝>𝑁/2 is required [33], but in the current scenario 𝑁=2; thus the condition reduces to 𝑝>1. The estimates trivially follow from Hölder’s inequality, applied on the reaction terms along with (56) and (57). This proves the global existence result.

4.3. Existence of Global Attractor in 𝐿²(Ω)

We can now use the standard theory [32] to prove that the semigroup for (20)–(22) possesses a global attractor. Thus, we can state the following result.

Theorem 14. Consider the three-species food-chain model described by (20)–(22), 1<𝑚₃<𝑚₄, and 0<𝑚₂≤1. There exists a global attractor 𝒜 in 𝐻, for the solution semigroup {𝑆(𝑡)}_𝑡≥0 generated by this model.

Proof. We have shown that the system is well posed via Theorem 6. Thus, there exists a well-defined semigroup {𝑆(𝑡)}_𝑡≥0∶𝐻→𝐻. Lemma 12 establishes the existence of bounded absorbing sets in 𝐻. Thus, given a sequence {𝑥_1𝑛(0)}^∞_𝑛=1 that is bounded in 𝐿²(Ω), we know that, for 𝑡>𝑡^∗, 𝑆(𝑡)𝑥_1𝑛(0)⊂𝐵⊂𝐻¹(Ω). (58) Here, 𝐵 is the bounded absorbing set in 𝐻¹(Ω). Now, 𝑡_𝑛>𝑡^∗, for any large enough 𝑛. Therefore, for 𝑡_𝑛, we have 𝑆𝑡_𝑛𝑥_1𝑛(0)⊂𝐵⊂𝐻¹(Ω). (59) This implies that we have the following uniform bound, via (57): 𝑆𝑡_𝑛𝑥_1𝑛(0)_𝐻1(Ω)≤𝐶₁. (60) By standard functional analysis theory [32], a subsequence, still labelled as 𝑆(𝑡_𝑛)(𝑥_1𝑛(0)), exists such that 𝑆𝑡_𝑛𝑥_1𝑛(0)⇀𝑥₁in𝐻¹(Ω), (61) which implies, via the compact Sobolev embedding of 𝑉↪𝐻, that 𝑆𝑡_𝑛𝑥_1𝑛(0)⟶𝑥₁in𝐿²(Ω). (62) This yields the asymptotic compactness of the semigroup {𝑆(𝑡)}_𝑡≥0 in 𝐻. Similar analysis is possible for components 𝑥₂ and 𝑥₃. We have the existence of an absorbing set via Lemma 12 and the asymptotic compactness via (62). These in conjunction yield the existence of the global attractor, via standard methods [32].

5. Numerical Simulation and Attractor Reconstruction

We now carry out numerical simulations of (2)–(4). Our goal is to firstly show spatiotemporal chaos, even if 𝑚₃<𝑚₄=2. This is indeed the case in certain parameter regime; see Figure 3. We also want to numerically validate the finite time blow-up results in 1D as well as in 2D. For 1D simulation, The system is simulated in MATLAB (R2010) via the PDEPE solver for partial differential equation. The algorithm essentially uses a central difference in space and an implicit time stepping method. All the simulations performed have been refined several times on spatial grids with 300, 600, 900, and 1200 points on the domain of length 𝐿=𝜋. These refinements lead to the same general shape and structure of the figures. Table 1 lists the parameter values used in the simulations.

Table 1 

Parameters used in simulations.

In order to explore the dynamics of the model in 2D, we use a finite difference method. A forward difference scheme is used for the reaction terms. For the diffusion terms, a standard five-point explicit finite difference scheme is used. The numerical simulation is carried out at different time levels for two dimensional spatial model system. The system of equations is numerically solved over 200 × 200 mesh points, on a domain of size 𝐿_𝑥×𝐿_𝑦, with spatial resolution Δ𝑥 = Δ𝑦 = 1 and time step Δ𝑡=0.1. The initial condition used is a small perturbation about (2.1,2.9,1.9), and the boundary conditions used are the no flux Neumann conditions. Values for 𝐿_𝑥 and 𝐿_𝑦 are taken to be 200 each. Note this is fine, as the blow-up results hold for Neumann boundary conditions as well. See Figure 1 to visualise the blowup in 1D, and Figure 2 to visualise the blowup in 2D.

(a)

(b)

Figure 1 

Blowup is numerically demonstrated for the system (2)–(4), in these plots, when the mutual interference parameters are 𝑚1=1,𝑚2=1,and𝑚3=2. This is exactly the model proposed in [4]. We see blowup in the ODE model in the left, at 𝑡=24, and PDE model in the right, at 𝑡=41. Note here we take 𝑐=0.07>(𝑤3/𝐷3)=0.06.

Figure 2 

Here, finite time blowup for the model system in two dimensional spatial domain is seen. This occurs at time 𝑡=29.4.

(a)

(b)

Figure 3 

The densities of the 𝑥1and𝑥2 species are shown as contour plots in the 𝑥−𝑡-plane. The corresponding parameters are listed in Table 1. The long-time dynamics (we run up to time 𝑡=3000) is seen to be that of spatiotemporal chaos. This shows that one can have interesting dynamics even in the case that 𝑚3<𝑚4, whilst avoiding finite time blowup.

Our next goal is to quantify the chaotic dynamics present in the system. The dynamics of the three-species model described via (2)–(4) is depicted through the reconstruction of the trajectories in phase space with the aide of delay-time method. This method, proposed by Takens and Mane [34] ensures equivalence between the topological properties of actual and reconstructed attractors. Numerical simulations were conducted for three values of the mutual interference parameter 𝑚₃. That is 𝑚₃=2, 𝑚₃=1.9, and 𝑚₃=1.2. In all cases, 𝑚₄ was set equal to 2. The system of PDE (2)–(4) was solved again via the MATLAB routine PDEPE. We then fixed certain points in space, and recorded the trajectory of the solution in time, yielding a time series. Very interesting behavior was observed; this consists of suppression of chaos. In fact for a value of 𝑚₃=2, chaos is found when 𝑐=0.03 in Table 1 in Upadhyay et al. [35]; see Figure 4(a).

(a)

(b)

(c)

Figure 4 

The figures above show the dependence of the dynamics on 𝑚3. For 𝑚3=2, one sees chaotic dynamics with the largest Lyapunov exponent equal to 2.23. Reducing 𝑚3 to 1.9 brings the largest lyapunov exponent equal to 0.26. Further reducing it to 1.2 brings the dynamics to a limit cycle.

We next calculated the dimension of this attractor, using a combination of singular-value decomposition and the Grassberger-Procaccia algorithm, [36]. We essentially want to embed the time series in an 𝑚 dimensional embedding space of embedded vectors, which take the typical form 𝑦_𝑠=𝑥𝑡₀+(𝑠−1)𝑙,𝑥𝑡₀+()𝑥𝑡𝑠−1𝑙+𝜏,…,₀+(𝑠−1)𝑙+(𝑚−1)𝜏^𝑇. (63)

Here, 𝜏 is called the “delay-time," 𝑙 is the sampling interval, and 𝑤=(𝑚−1)𝜏 is the “window length" which represents the time spanned by each embedding vector. This algorithm entails choosing 𝑚 and the delay, so that the window is sufficiently large. We then perform a singular value decomposition and calculate the correlation dimension, from the correlation integral. We continue to refine the embedding dimension till we have maximized the straight line parts in the log-log plot of the correlation integrals. The dimension of this attractor is approximately equal to 2.3 while the value of the largest Lyapunov exponent, estimated with the method proposed by Lai and Chen [37], is equal to 2.23. This clearly indicates that the dynamics is chaotic. Decreasing slightly the value of 𝑚₃ to 1.9, the dynamics remained chaotic and the dimension of the attractor remained almost unchanged, that is, equal to 2.3 while the largest Lyapunov exponent had decreased to 0.26; see Figure 4(b). This tells us that the divergence of close trajectories becomes now very slow. The important point here is that the proposed modification to the Parshad-Upadhyay model and subsequently the classical Upadhyay-Rai models still exhibits chaos in certain parameter regime. That is for 1.9<𝑚₃<𝑚₄=2. Note since 𝑤₀<𝑤₁, this parameter set is not very realistic biologically. However, when we decrease further the value of 𝑚₃ to 𝑚₃=1.2, the dynamics of the system changes radically and becomes oscillatory with single period; this is depicted by a limit cycle in phase space; see Figure 4(c). From this numerical analysis, it seems plausible that the mutual interference parameter 𝑚₃ plays a crucial role in the dynamics of the three-species model. Setting a value of the parameter 𝑚₃ (𝑚₃=1.2), very different from the classical one 𝑚₃=2, suppresses the chaos predicted by the model when 𝑚₃=2. However, keeping the parameters close to 2, that is 1.9<𝑚₃<2, still retains chaotic dynamics, whilst eliminating the possibility of finite time blowup.

6. Discussion

Various nonlinear effects, like structural instability, dissipative structures, dynamical chaos, and finite time blowup, do exist in food chain models. However, the problem of detecting these phenomena in real ecosystems is far from being well understood. Various studies have suggested that the biology of the top generalist predator has an important role to play as far as the dynamics of food-chain models is concerned [10, 38]. The nature of interactions between the populations and communities may be responsible for chaos evading its capture in the wild. In this work, we have shown that a large class of three-species food-chain models [4, 7, 9] can blowup in finite time, for certain parameter regime and large enough initial data. We reiterate that the situation can be worse than unbounded growth if 𝑐>𝑤₃/𝐷₃, as reported earlier [5]. We propose a possible correction based on the analysis of the way mating is modeled in (4). That is, if we restrict 𝑚₃ to the range 1<𝑚₃<𝑚₄, then finite time blowup is prevented, and this is independent of the sign of 𝑐−𝑤₃/𝐷₃. However, if 𝑐−𝑤₃/𝐷₃<0 and 𝑥₂=0, the top predator can still persist, for 1<𝑚₃<𝑚₄. Thus, this formulation plausibly models a true generalist predator 𝑥₃, one that does not go extinct, even if its favorite food source does. Finally, by varying 𝑚₃ in the range 1<𝑚₃<𝑚₄, one can adapt the new model to populations where the numbers of males and females vary. If the numbers are very close, we take 𝑚₃=2−𝜖, 𝜖≪1. One could also model a population where the males far outnumber the females or vice versa. Here one might take 𝑚₃=1.5, saying we have roughly 10 females per 100 males in a population or vice versa. Thus, making the case for the model in structured populations, where such demographics may be dominant. Also, simulations performed by varying the mutual interference parameter 𝑚₃ indicate clearly that there is a bifurcation point 𝑚^∗₃, in the 𝑚₃ parameter space, such that, chaotic dynamics is seen for 𝑚₃>𝑚^∗₃ and oscillatory dynamics for 𝑚₃<𝑚^∗₃. It would be interesting to obtain the actual value of the mating parameter 𝑚₃ from real field data. Comparing this “real value" to 𝑚^∗₃ might then be a means to again answer the long standing question of why chaos has evaded its capture in the wild. Lastly, we remark that the blow-up result obtained via Theorem 2 is sufficient but not necessary. It is not optimal in terms of data. That is, it might still be possible for blowup to occur if the data does not meet the largeness condition _Ω𝑥₃₀(𝑥)𝑊₁(𝑥)𝑑𝑥>𝜆₁𝑑₃1𝑐−𝑤₃/𝐷₃^{1/(𝑚3−1)}. (64)

For the classical diffusive Upadhyay-Rai model [4], we show blowup numerically; see Figure 1. Here, we have precisely all the parameters previously stated, including the first eigenvalue and eigenfunction of the Dirichlet laplacian. We see that, for the parameters chosen to demonstrate blowup in Figure 1, _Ω𝑥₃₀(𝑥)𝑊₁(𝑥)𝑑𝑥=0.8<1.5=𝜆₁𝑑₃1𝑐−𝑤₃/𝐷₃. (65) Thus blowup takes place, even for data smaller than that required by the largeness condition. It would thus be very interesting to try and sharpen this largeness requirement on the data, perhaps via employing certain other techniques available in the blow-up literature, than the version of the first eigenvalue method that we have resorted to currently.