Finite Time Blowup in a Realistic Food-Chain Model

1 Department of Mathematics and Computer Science, Clarkson University, Potsdam, NY 13699, USA 2King Abdullah University of Science and Technology, Applied Mathematics and Computational Science, Thuwal 23955-6900, Saudi Arabia 3 Department of Applied Mathematics, Indian School of Mines, Dhanbad, Jharkhand 826004, India 4 School of Basic Sciences, Indian Institute of Technology Mandi, Mandi, Himachal Pradesh 175 001, India


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 foodchain 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][5][6][7][8][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 ISRN Biomathematics 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 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][12][13][14][15][16][17][18][19][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  3 <  3 .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,  3 <  3 , 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]  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.

Food Chain Model System
In this section, we recount the model from [9].Consider a situation where a prey population  1 is predated by individuals of population  2 .The population  2 in turn serves as a favourite food for individuals of population  3 .This interaction is represented by the following system of a simple prey-specialist predator-generalist predator interaction with the inclusion of spatial spread: where The parameters   for  = 1, 2, 3 are mutual interference parameters that model the intraspecific competition among predators when hunting for prey [22][23][24][25][26][27].The problem is posed on Ω ⊂ R 2 .Ω is bounded, and Ω is assumed to be smooth.We consider the Dirichlet boundary conditions We also assume suitable initial conditions.In this model, prey population of size  1 serves as the only food for the specialist predator population of size  2 .This predator population, in turn, serves as a favorite food for the generalist predator population of size  3 .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  2 and the generalist predator  3 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. 1 is the intrinsic growth rate of the prey population  1 ,  2 is the intrinsic death rate of the predator population  2 in the absence of the only food  1 ,  measures the rate of self-reproduction of generalist predator  3 , and  0 ,  1 ,  2 , and  3 are the maximum values which per capita growth rate can attain. 1 measures the strength of intraspecific competition among the individuals of the prey species  1 . 0 and  1 quantify the extent to which environment provides protection to the prey  1 and may be thought of as a refuge or a measure of the effectiveness of the prey in evading a predator's attack. 2 is the value of  2 at which per capita removal rate of  2 becomes  2 /2. 3 represents the residual loss in  3 population due to severe scarcity of its favourite food  2 .For  1 =  2 = 1 the coefficient  0 /( 1 +  0 ) 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  1 (the prey population at any instant of time), it gives the attack rate on the prey per predator.Denote ( 1 ) = ( 0  1 )/( 1 +  0 ), when  1 → ∞, ( 1 ) →  0 , which is the maximum that it can reach.The third term ( 2  2  3 )/( 2 +  2 ) on the righthand side of (3) represents the per capita functional response of the invertebrate predator  3 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  3 , 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. 3 measures the limitation on the growth of the generalist predator  3 by its dependence on per capita availability of its most favorite prey  2 .
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][25][26] that mutual interference is generally a "stabilizing" process.The parameters  1 ,  2 , and  3 are diffusion coefficients of the populations.

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  3 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 R, and let  be a nonnegative function.Then the following holds whenever  is a convex function: 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  > ( 3 / 3 ) and initial data  30 () large enough, that is there exists a finite time  * > 0, such, that for any  3 > 1, Proof.Consider the equation for  3 rewritten as Since  2 is assumed positive, we obtain Thus, Now we multiply (11) by  1 , consider () = ∫ Ω  3 (, ) 1 (), and integrate by parts to obtain Via assumptions on the eigenfunction  1 , it is seen to be a probability measure.That is, we can define a probability measure () =  1 on Ω and () = 0 on R 2 /Ω.Now we consider the convex function ( 3 ) =   3 3 (note that any  3 > 1 works as the  3 in question is always positive).Thus, the application of Jensen's inequality yields Thus, substituting the previous equation into (12) yields Simple integration in time of the previous equation yields that () blows up in finite time; that is where and  * ≤  * * .
After some algebra, we obtain that for blowup the data has to satisfy .
This proves the theorem.
Remark 3. Note that the blowup of () immediately implies that the  ∞ norm of  blows up as follows: Remark 4. We also show blow-up in the Upadhyay-Rai model proposed in [4], that is, the temporal case with  1 = 1,  2 = 1, and  3 = 2. Also, for the model in [4], the blowup time is given by Remark 5. We note that finite time blowup is also possible in the same weighted  1 (Ω) norm, for Neumann boundary conditions.In this case, the eigenfunction  1 will have to be differently defined, that is, satisfying Neumann boundary conditions itself.

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  2 3 .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) 3 × (1/2) 3 , 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 , 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  3 power.The primary role of the previously stated modification is that now the reaction term in ( 4) is 3 , and because we enforce  3 <  4 , 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, The problem is posed on Ω ⊂ R 2 .Ω is bounded, and Ω is assumed to be smooth.We consider the Dirichlet boundary conditions 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  0 =  1 , as commonly done in numerical simulations, [5].We also assume 0 <  2 ≤ 1.This facilitates the estimates to follow.Restricting  3 to the range 1 <  3 <  4 has distinct advantages over the old model ( 2)-(4).Firstly, finite time blowup is prevented, and this is independent of the sign of  −  3 / 3 .Secondly, in the absence of  2 , if  −  3 / 3 < 0, the dynamics of ( 2)-( 4) is quite boring, and  3 decays to 0, [5,6].However, if 1 <  3 <  4 , even if  2 = 0 and  −  3 / 3 < 0,  3 can still persist.Thus, this formulation plausibly models a true generalist  3 .One that does not go extinct, even if its favorite food source does.
The opposite signed nonlinearities here cause problems in making direct  1 estimates and proving the existence of bounded absorbing sets, for the  2 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  ≤  3 /(ℎ + V).Thus, our results in [9] are only true in the parameter regime, such that V < ( 3 /) −  3 .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, A ⊂ , for this semigroup is an object that has the following properties: 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 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, B, in a reflexive Banach space, , is called a bounded absorbing set if, for each bounded subset  of , there is a time,  = (), such that () ⊂ B for all  > .The number  = () is referred to as the compactification time for ().This is essentially the time after which the semigroup compactifies.

Global Existence and Existence of Absorbing Sets in 𝐿 2 (Ω).
We begin by multiplying (20) by  1 and integrating by parts over Ω.This yields We then use Hölder's inequality followed by Young's inequality to obtain Next, using the compact Sobolev embedding,  1 (Ω) →  2 (Ω), and the positivity of  1 and  2 , we find that where  depends explicitly on  1 and  1 .Thus, application of Gronwall's Lemma gives us the following estimate: Therefore, there exists time  1 given explicitly by such that, for all  ≥  1 , the following estimate holds uniformly: We now make a local in time estimate for ∇ 1 .By integrating (26) in the time interval [ 1 ,  1 + 1], we obtain Thus, via a mean value theorem for integrals, there exists a time  2 ∈ [ 1 ,  1 + 1] such that the following estimate holds: We now move on to showing the existence of an absorbing set for  2 in  2 (Ω).By multiplying (20) by  1 and ( 21) by  0 , adding the two together, and setting  =  1  1 +  0  2 , we obtain 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 Thus, by using the compact Sobolev embedding,  1 (Ω) →  2 (Ω), we obtain Next, we multiply the previous equation by   and integrate in time from 0 to  to find where  1 ( 10 , ) =  − ∫  0  − 1  ‖ 1 (0)‖ 2 2  + (/ 1 ).This follows via (30).
We now make the following estimate via time integration of (26), after multiplying through by   : This follows via (30).We now substitute the above into (36) and use the form of  1 ( 10 , ) to obtain ‖‖ 2  2 ≤  − ‖ (0)‖ Therefore, there exists time  3 , given explicitly by such that, for all  ≥  3 , the following estimate holds uniformly: Therefore, lim sup where  is a constant that is independent of the time and the initial data.We next derive local in time estimates for ∇ 2 .We proceed by multiplying ( 21) by  2 and integrating by parts over Ω to obtain We now use Hölder's inequality, positivity of solutions, the embedding  2 → The existence of an absorbing set for  3 in  2 (Ω) follows easily via the form of (22).Standard estimates imply the existence of a time  5 given explicitly by such that, for all  ≥  5 , the following estimates hold uniformly: 3 We now make the  1 estimate.We take the inner product of (20) with −Δ 1 and integrate by parts over Ω.Thus, we obtain We apply the Cauchy-Schwartz inequality, the Cauchy with epsilon, and Young's inequalities on the last term, along with the embedding  2 →  2 2 (via assumption on  2 ) to obtain This yields We now recall the following lemma.
Lemma 13 (the uniform Gronwall lemma).Let , , and ℎ be nonnegative functions in  1  [0, ∞; ).Assume that  is absolutely continuous on (0, ∞) and the following differential inequality is satisfied: If there exists a finite time  1 > 0 and some  > 0 such that for any  >  1 , where , , and  are some positive constants, then We apply the previous lemma to (50) with and  1 =  * ,  = 1, and we use the estimates via (31) and Lemma 12 to obtain The  1 estimates on the other components,  2 and  3 , are obtained similarly.
Recall that in R 2 we have via the Sobolev embedding that  1 (Ω) →   (Ω), for all , lim sup and similarly, for  2 and  3 , we have lim sup Using the previous equation it is easy to show that the reaction terms  1 ,  2 , and  3 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.
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 Here,  is the bounded absorbing set in  1 (Ω).Now,   >  * , for any large enough .Therefore, for   , we have This implies that we have the following uniform bound, via (57): By standard functional analysis theory [32], a subsequence, still labelled as (  )( 1 (0)), exists such that which implies, via the compact Sobolev embedding of  → , that This yields the asymptotic compactness of the semigroup {()} ≥0 in .Similar analysis is possible for components  2 and  3 .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].

Numerical Simulation and Attractor Reconstruction
We now carry out numerical simulations of ( 2)-(4).Our goal is to firstly show spatiotemporal chaos, even if  3 <  4 = 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.
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.
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  3 .That is  3 = 2,  3 = 1.9, and  3 = 1.2.In all cases,  4 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  3 = 2, chaos is found when  = 0.03 in Table 1 in Upadhyay et al. [35]; see Figure 4(a).
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   = ( ( 0 + ( − 1) ) ,  ( 0 + ( − 1)  + ) , . . .,  ( 0 + ( − 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  3 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 <  3 <  4 = 2.Note since  0 <  1 , this parameter set is not very realistic biologically.However, when we decrease further the value of  3 to  3 = 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  3 plays a crucial role in the dynamics of the three-species model.Setting a value of the parameter  3 ( 3 = 1.2), very different from the classical one  3 = 2, suppresses the chaos predicted by the model when  3 = 2.However, keeping the parameters close to 2, that is 1.9 <  3 < 2, still retains chaotic dynamics, whilst eliminating the possibility of finite time blowup.

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  >  3 / 3 , 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  3 to the range 1 <  3 <  4 , then finite time blowup is prevented, and this is independent of the sign of  −  3 / 3 .However, if  −  3 / 3 < 0 and  2 = 0, the top predator can still persist, for 1 <  3 <  4 .Thus, this formulation plausibly models a true generalist predator  3 , one that does not go extinct, even if its favorite food source does.Finally, by varying  3 in the range 1 <  3 <  4 , one can adapt the new model to populations where the numbers of males and females vary.If the numbers are very close, we take  3 = 2−,  ≪ 1.One could also model a population where the males far outnumber the females or vice versa.Here one might take  3 = 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  3 indicate clearly that there is a bifurcation point  * 3 , in the  3 parameter space, such that, chaotic dynamics is seen for  3 >  * 3 and oscillatory dynamics for  3 <  * 3 .It would be interesting to obtain the actual value of the mating parameter  3 from real field data.Comparing this "real value" to  * 3 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 ISRN Biomathematics  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. .(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, ) . (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.

Figure 3 :
Figure 3: The densities of the  1 and  2 species are shown as contour plots in the  − -plane.The corresponding parameters are listed in Table1.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.

Figure 4 :
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.

Table 1 :
Parameters used in simulations.