Complex dynamic of plankton–fish interaction with quadratic harvesting and time delay

Modeling Earth Systems and Environment, Nov 2016

The present paper deals with a food chain model consisting of three species phytoplankton, zooplankton and fish. We have divided the present paper into two parts. In the first part, we have assumed that the fish population is harvested using a non-linear harvesting function. Considering this rate of harvesting ‘E’ as a control parameter, we have estimated different ranges of harvesting parameter for maintaining the sustainability in the plankton ecosystem. Moreover, the bifurcation analysis of the system is carried out using normal form theorem by taking ‘E’ as bifurcation parameter. In the second part, a digestion delay corresponding to zooPlankton–fish interaction is introduced for more realistic consideration of the real world problem. Taking harvesting parameter in the stability range, the effect of time delay on the given system is investigated. This research demonstrate that for a certain range of delay, system enters into the excited state with the existence of stability switches which seems new findings for the Plankton–fish system. Explicit results are derived for stability and direction of the bifurcating periodic solution by using normal form theory and center manifold arguments. To validate our analytical findings numerical simulations are also executed.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs40808-016-0248-x.pdf

Complex dynamic of plankton–fish interaction with quadratic harvesting and time delay

Model. Earth Syst. Environ. Complex dynamic of plankton-fish interaction with quadratic harvesting and time delay Amit Sharma 0 1 2 Anuj Kumar Sharma 0 1 2 Kulbhushan Agnihotri 0 1 2 0 LRDAV College , Jagraon, Punjab , India 1 DAV Institute of Engineering and Technology , Jalandhar, Punjab , India 2 SBSSTC , Ferozpur, Punjab , India The present paper deals with a food chain model consisting of three species phytoplankton, zooplankton and fish. We have divided the present paper into two parts. In the first part, we have assumed that the fish population is harvested using a non-linear harvesting function. Considering this rate of harvesting 'E' as a control parameter, we have estimated different ranges of harvesting parameter for maintaining the sustainability in the plankton ecosystem. Moreover, the bifurcation analysis of the system is carried out using normal form theorem by taking 'E' as bifurcation parameter. In the second part, a digestion delay corresponding to zooPlankton-fish interaction is introduced for more realistic consideration of the real world problem. Taking harvesting parameter in the stability range, the effect of time delay on the given system is investigated. This research demonstrate that for a certain range of delay, system enters into the excited state with the existence of stability switches which seems new findings for the Plankton-fish system. Explicit results are derived for stability and direction of the bifurcating periodic solution by using normal form theory and center manifold arguments. To validate our analytical findings numerical simulations are also executed. Plankton; Fish; Quadratic harvesting; Time delay; Stability switches; Normal form; Center manifold theorem Introduction A major concern in population ecology is to understand how a population of a given species influences the dynamics of populations of other species. The dynamic relationship between predator and their prey has long been and continue to be one of the dominant themes in mathematical ecology due to its universal existence and importance (Berryman and Millstein 1989). The pioneering work of May (1976) has established some mathematical models based on certain ecological principles to explore the complexity of the ecological system. The research of the last two decades have demonstrated that very complex dynamics can arise in three or more species food chain models (Hastings and Powell 1991; Klebanoff and Hastings 1993; Rai and Upadhyay 2004; Gakkhar and Naji 2005; Upadhyay and Chattopadhyay 2005) including quasi-periodic or chaos. A great number of theoretical studies indicates that, indeed, plankton systems are, in principle, capable to generate their own chaos. Preypredator models (e.g., phytoplankton-zooplankton interaction) in seasonally varying environments have been proved to be chaotic (Inoue and Kamifukumoto 1984; Klebanoff and Hastings 1994; Kuznetsov and Rinaldi 1996). Chaotic behavior in tritrophic food chain interactions have been shown in Gilpin (1979), Hogeweg and Hesper (1978), Scheffer (1991). It is well established that toxin has great impact on phytoplankton-zooplankton interaction and can be used as a mechanism of controlling complexity (bloom) of the plankton ecosystem (Pal et al. 2007; Chakarborty et al. 2008; Sarkar and Chattopadhyay 2003; Mukhopadhyay and Bhattacharyya 2006; Hansen 1995; Ives 1987; Buskey and Hyatt 1995). Recently, Upadhyay and Chattopadhyay (2005), Upadhyay et al. (2007), Upadhyay and Rao (2009) have studied a series of mathematical models representing the prey(TPP)-specialist predator(Zooplankton)-generalist predator(Fish) interaction in the context of deterministic chaos in the plankton system. They have studied the role of toxin as a control parameter in stabilizing the plankton system by using different response functions for fish grazing. Further, the introduction of strictly planktivorous fish to lakes can alter plankton communities via cascading interactions in food webs, where strong fish predation on zooplankton leads to reduction of their grazing pressure on algae. Many limnological studies have focused on this dramatic impact of fish on plankton communities (Carpenter et al. 1987; Leavitt and Findlay 1994; Pace et al. 1999). Furthermore, enhancement of piscivorous fish by stocking and catch restrictions has been attempted in several eutrophic or hypertrophic waters, and has resulted in the elimination of excess algae (Carpenter et al. 1987; Rinaldi and Solidoro 1998; Benndorf et al. 1984; McQueen et al. 1989; Lazzaro 1987). Recently, research carried by authors in Strock et al. (2013) reveals that the introduction of white perch in eutrophic lake reduces the algal pigments concentration and results high label of stability in trophic interaction due to the cascading effects of white perch. Attayde et al. (2010) in their model, also supports the idea that omnivory decreases the amplitude of limit cycles and increases the persistence in plankton dynamics. Keeping in view the above findings, here, we have considered a tri-trophic food chain model of Plankton–fish interaction. We have assumed that the generalist predator (fish) is harvested using a nonlinear harvesting function (Feng 2014). Our major concern in this research to determine the impact of non-linear harvesting on the Plankton– fish interaction by predicting different ranges of harvesting parameter. So far many researcher (Upadhyay et al. 2007; Upadhyay and Rao 2009, references there in) have studied that toxicated phytoplankton may be used as controller and can stabilize the Plankton–fish dynamic. In this paper, we have introduced quadratic harvesting of the fish population and proposed different ranges of harvesting for regulating agencies so that they can utilize it for future harvesting along with the stabilization of the ecosystem. The organization of our paper is as follows: In Sect. 2, we have developed a mathematical model followed by its properties in Sect. 2.1. The stability and bifurcation analysis of the given model system with rate of harvesting as a bifurcation parameter is presented in Sect. 3. In Sect. 4, we have assumed that there is a gestation delay in fish population and considering this time delay as a bifurcation parameter, we have discussed its possible impact on the Plankton–fish dynamic. The stability and bifurcation properties are provided in Sect. 5. The justification of our analytical findings by numerical simulation and concluding remarks are given in Sects. 6 and 6.2. The mathematical model Let p(t) be population density of the toxin producing phytoplankton (prey) which is predated by individuals of specialist predator zooplankton of population density z(t) and this zooplankton population, in turn, serves as a favorite food for the generalist predator (fish of mollusca) population of size f(t). This tri-trophic interaction is represented by the following system of ordinary differential equation as, >>>>>8 ddpt ¼ rp d1p2 a1 a þp p z; >< dz p z p >> dt ¼ c1a1 a þ p z d2z a2 b þ z f h a þ p z; ð1Þ :>>>> ddft ¼ d3f þ c2a2 f Ef 2; where r be the growth rate of the TPP species, a1 be the maximum ingestion rate by zooplankton population, c1 is the fraction of biomass converted to zooplankton, d1, d2, d3 be the mortality rates of phytoplankton, zooplankton and fish population respectively. Let a2 be the rate at which fish graze on specialist predator (zooplankton) and c2 is the corresponding biomass conversion by fish predator. a, b and d are the half saturation constants. The effect of harmful phytoplankton species on zooplankton is modeled by Holling-II type response function ‘ðaþp pÞ’. It is assumed that the fish population is harvested quadratically at constant rate of harvesting E. Stability properties of equilibria The given model system (1) has the following equilibrium points in the closed first octant R3þ ¼ ðp; z; f Þ : p 0; z 0; f 0 R0 ¼ ð0; 0; 0Þ, a trivial equilibrium always exist and unstable, R1 ¼ ðc1a1adh2 d2 ; a11 ðr d1p Þða þ p Þ; 0Þ, the boundary equilibrium exists when the conditions ðc1a1 hÞ [ d2 and p \ dr1 are satisfied and The interior equilibrium point R ¼ ðp ; z ; f Þ of (1) exists if and only if conditions p \ dr1 and d þzz \ cd2a32 are satisfied and is given by, d2 ¼ baþ2fz , z ¼ ðr d1pa1Þðaþp Þ and In this paper, our main interest is to investigate the stability of the interior equilibrium R for coexistence of the population. Using the transformation, p ¼ x1 þ p , z ¼ x2 þ z and f ¼ x3 þ f , the linearization of given system (1) in the neighborhood of R gives the following system, 3 2 2 þb030x2 þ b210x1x2 þ b021x2x3 þ Theorem 2.1 Let B2 [ 0. Then (3) will have a pair of purely imaginary roots provided H B1B2 B3 ¼ 0 and system (1) will undergo a hopf-bifurcation if ½ddE ðHÞ E¼E0 6 ¼ 0 where E0 is the value of E correspond to H ¼ 0. Proof is obvious and is left for the reader. a210 ¼ ða2þapa1Þ3, b100 ¼ ðca1þa1pazÞ2 c020 ¼ ð2dcþ2az2Þd3f , c002 ¼ 2E, c011 ¼ ðddþc2za2Þ2, c030 ¼ 6ðcd2þaz2dÞf4 , 2dc2a2. c021 ¼ ðdþz Þ3 Rewriting the given system in the form, X_ ¼ MX þ R, we have and R is the matrix containing the nonlinear terms of (2). The characteristic equation corresponding to the matrix M is, k3 þ B1k2 þ B2k þ B3 ¼ 0: Analysis of bifurcating solutions In the last section, we have chosen E as the bifurcation parameter and deduced that the system undergoes a Hopfbifurcation as E passes through some threshold value E0. In this section we perform a detailed analysis about the bifurcating solutions (Hassard et al. 1981). Let E1 ¼ E0 E and denote the matrix M as MðEÞ ¼ FX ðX ; EÞ, where the given system is of the form X_ ¼ FðX; EÞ. Let q(E) and q^ E ð Þ be the eigenvectors of M(E) and MtðEÞ respectively corresponding to simple eigenvalues k and k. So we have MðEÞq ¼ kðEÞq, For an imaginary value of k ¼ ir0 and q ¼ ðq1; q2; q3Þ, we have q1 ¼ a010b001 ¼ f1 þ if2, q2 ¼ a100b001 þ ix0b001 ¼ g1 þ ig2, q3 ¼ a100b001 a010b100 x20 ix0ða100 þ b001Þ ¼ h1 þ ih2, We normalize q^ relative to q, so that \q^; q [ ¼ 1: q^1ðf1 þ if2Þ þ q^2ðg1 þ ig2Þ þ q^3ðh1 þ ih2Þ ¼ 1: Solving (6) and (7), we have where k is a constant. So, q^1, q^2 and q^3 can be determined by taking complex conjugates of (8). Using the notation in Hassard et al. (1981), we substitute At bifurcation parameter E1 ¼ 0, we have, u_ ¼ ix0uðtÞ þ q^ð0Þ:fð2ReðuqÞ þ Wð0ÞÞ ¼ ix0 þ q^ð0Þ:fðu;uÞ; f01 ¼ u2ða200q12 þ a110q1q2Þ þ uuð2a200q1q1 þ a110ðq1q2 þ q1q2ÞÞ þ u3 a3600 q13 þ a2610 q12q2 þ ; f02 ¼ u2ðb200q12 þ b020q22 þ b110q1q2 þ b011q2q3Þ þ uuð2b200q1q1 þ 2b020q2q2 þ b110ðq1q2 þ q1q2Þ þ b011q2q3 þ b011q2q3Þ þ u2ðb200q12 þ b020q22 þ b110q1q2Þ þ u3ðb300q13 þ b030q23 þ b210q12q2 þ b021q22q3Þ þ ; q^1ðf1 if2Þ þ q^2ðg1 ig2Þ þ q^3ðh1 ih2Þ ¼ 0: ð7Þ f03 ¼ u2ðc020q22 þ c002q32 þ c011q2q3Þ þ uuð2c020q2q2 þ 2c002q3q3 þ c011ðq2q3 þ q2q3ÞÞ þ u2ðc020q22 þ c002q32 þ c011q2q3Þ þ u2u c030 q2q2 þ c021 q2q3 q2q2q3 2 6 þ 3 þ uu2 c030 q2q22 þ c021 q22q3 q2q2q3 6 6 þ 3 Now taking the dot product of f0 and q^ð0Þ and expanding we have, u_ ¼ ix0u þ q^1f01 þ q^2f02 þ q^3f03; ¼ ix0u þ 2g20u2 þ 2g02u2 þ g11uu þ 6g30u3 þ 61g03u3 1 1 1 From (9), (10), (11) and (12), we obtain g20 ¼ 2½q^1ða200q1 þ a110q1q2Þ þ q^2ðb200q12 þ b020q22 2 þ b110q1q2 þ b011q2q3Þ 2 2 þ q^3ðc020q2 þ c002q3 þ c011q2q3Þ ; g02 ¼ 2½q^1ða200q12 þ a110q1q2Þ þ q^2ðb200q12 þ b020q22 þ b110q1q2 þ b011q2q3Þ þ q^3ðc020q22 þ c002q32 þ c011q2q3Þ ; g11 ¼ q^1ð2a200q1q1 þ a110ðq1q2 þ q1q2ÞÞ þ q^2ð2b200q1q1 þ 2b020q2q2 þ b110ðq1q2 þ q1q2Þ þ b011q2q3 þ b011q2q3Þ þ q^3ð2c020q2q2 þ 2c002q3q3 þ c011ðq2q3 þ q2q3ÞÞ; h g21 ¼ 2 q^1 a3200 q21q1 þ a2610 q21q2 þ a2310 q1q1q2 þq^2 b3200 q21q1 þ b0230 q22q2 þ q^3 c030 q2q2 þ c021 q2q3 q2q2q3 : ð13Þ 2 6 þ 3 Using the notation of Hassard et al. (1981) we write ing20g11 2jg11j2 jg032j2o where gij are given by (13). Theorem 3.1 From the above obtained quantities, it can be found that bifurcation is supercritical (stable) if l2 [ 0 and subcritical (unstable) if l2\0. Coexistence in the presence of digestion delay In most of ecosystems, both natural and manmade, population of one species does not respond instantaneously to the interactions of other species. Sharma et al. (2014) mentioned that animals must take some times to digest their foods before further activities and responses. To incorporate this idea in the modeling approach, we have considered a gestation delay in the zooPlankton–fish interaction. This delay represents the time lag required for gestation of predator which is based on the assumption that the rate of change of predator depends on the number of prey and of predators present at some previous time (Sarwardi et al. 2012). Further time delays in ecological system has destabilizing effects (Sharma et al. 2014; Cushing 1977; Gopalsamy 1992; Kuang 1993; Macdonald 1989; Beretta and Kuang 2002; Das and Ray 2008; Sharma et al. 2015). So, the main aim of this section is to determined whether the dynamic of given system under investigation in stability region for some parameter values changes in the presence of digestion delay or not. The model is governed by the following system of the nonlinear delay differential equation, >>>>>><8> dddzpt ¼ rp d1pp2 a1 a þp p z; z p >> dt ¼ c1a1 a þ p z d2z a2 b þ z f h d þ p z; ð14Þ >>>>>: ddft ¼ d3f þ c2a2 sÞ Ef 2; b þzðtzðt sÞsÞ f ðt The initial condition of the system (14) take the form pðhÞ ¼ /1ðhÞ, zðhÞ ¼ /2ðhÞ and f ðhÞ ¼ /3ðhÞ, h 2 ½ s; 0 , /ið0Þ 0, /iðhÞ 2 Cð½ s; 0 ; R3þÞ. Theorem 4.1 The system (14) remained uniformly bounded in P, where P ¼ ½ðpðtÞ; zðtÞ; f ðt sÞÞ : 0 pðtÞ 1g ðdr21 2ð1cþ2ac22Ea2ÞÞ and g ¼ minð1; r; d2 Proof From (14) for all t 2 ½0; 1Þ, dp d1 dt rpðtÞð1 r pðtÞÞ; or it can be obtain that, lim supt!1 pðtÞ dr1. Let us consider a time dependent function, Using (14) in above expression, we can obtain 2ð1 þ c2a2Þ ; where g ¼ minð1; r; d2 c2a2Þ: c2a2E Now, applying the theorem of differential inequalities (Birkhoff and Rota 1882), we obtain 0\WðtÞ As t ! 1, we have 0 Now, using Taylor’s expension, (14) can be reduced to the following system of differential equations, þc002x32 þ c0011x2ðt sÞx3ðt sÞ þ The corresponding characteristic equation is, DðR ; irÞ ¼ e irsð B4ir2 þ B5ir þ B6Þ; Separating the real and imaginary parts, we obtain B4r2ÞcosðrsÞ þ B5rsinðrsÞ; r3 ¼ B5rcosðrsÞ þ ðB6r B4ÞsinðrsÞ; Eliminating s from above equations, it can be obtained that Setting z ¼ r2, (21) can be written as, Hðr; sÞ ¼ z3 þ pz2 þ qz þ r ¼ 0: 2 2B1B3 þ 2B4B6 þ B5; Now from the sign of B1, B2 and B3, it can be obtained that Eq. (22) possesses two positive roots say, z0 which implies (21) has two roots r ¼ pffizffi0ffiffi. Substituting this value in (19) and (20), the corresponding critical value of delay s ¼ sk at which R stability switches occurs are given by, B6r2ÞðB2r r3Þ þ B5r ðB1r2 B1r2ÞðB6r2 B4Þ þ B5r ðB2r Remark The local asymptotical stability of the non-delayed system around R using Routh–hurwitz criterion is already discussed in above section and it follows that the necessary and sufficient conditions for all roots of (3) having negative real parts ðs ¼ 0Þ are given by B2 [ 0 and B1B2 B3 [ 0, which means that (14) is remained asymptotically stable around R in the absence of gestation delay. To obtain the effect of digestion delay on the stability of equilibrium point R , let us consider the case when s 6¼ 0 and the characteristic equation can be written as DðR ; sÞ ¼ k3 þ B1k2 þ B2k þ B3 þ e ksðB4k2 þ B5k þ B6Þ: Now for k ¼ 0, DðR ; sÞ ¼ B3 B6 6¼ 0, it means that k ¼ 0 is not a root of the characteristic equation (16) if the interior equilibrium point R exists. And the conditions for stability of R when s varies are given by the following theorem. Theorem 4.2 If all roots of the Eq. (16) has negative real parts and, When Hðr; sÞ ¼ 0 has no positive root, the interior equilibrium point R is asymptotically stable for an arbitrary delay s. When Hðr; sÞ ¼ 0 has at least one positive root say r0, there exists a critical delay s0 [ 0 such that the interior equilibrium point R is asymptotically stable for seð0; s0Þ. When Hðr; sÞ ¼ 0 has at least two positive roots say r then there exist positive integer n, such that the equilibrium R switches n times from stability to instability to stability and so on such that R is locally asymptotically stable whenever se½0; s0þ [ [ ðsn 1; snþÞ and is unstable when The model system (14) undergoes a Hopf-bifurcation around R for every s ¼ sk . Where Hðr; sÞ ¼ f12ðr; sÞ þ f22ðr; sÞ f1ðr; sÞ ¼ ðB3 f2ðr; sÞ ¼ ðB4 Proof Let ir0ðr0 [ 0Þ is a pair of purely imaginary roots of (16), then we have G2ðn; r; sÞ ¼ 1 ðB4 B6r2ÞðB2r r3Þ þ B5r ðB1r2 B3Þ, we have r0 ðB3 B1r2ÞðB6r2 B4Þ þ B5r ðB2r r3Þ two sequence of delays skþ ¼ hþrþþ2kp and sk ¼ h þr12kp for which there are two purely imaginary roots of the (16). We will now study how the real parts of the roots of (16) vary as s varies in a small neighborhood of sþ and s . Let k ¼ n þ ir be a root of (16), then substituting k ¼ n þ ir in (18) and separating real and imaginary parts we have, 3nr2 þ B1ðn2 r2Þ þ B2n þ B3 e nsðB4 þ B5n þ B6ðn2 3n2r þ 2B1nr þ B2r þ e nsðB4 þ B5n þ B6ðn2 e nsðB5r þ 2B6nrÞcosrs; ð0; r ; s Þ such that nþðskþÞ ¼ 0, rþðskþÞ ¼ rþ, j ddnsþ js¼skþ This completes the proof of theorem. Stability of bifurcating periodic solution In this section we are employing the explicit formulae suggested by Hassard et al. Hassard et al. (1981) to determine the direction, stability and the period of the periodic solutions bifurcating from interior equilibrium R . Following the procedure of the computation in Hassard et al. (1981), we compute (see Appendix A for details of the computation) where the first four coefficients that we need for determining properties of the Hopf bifurcation are of the form 8 2 >>> g20 ¼ d1 s0fD11 þ q1 D21 þ q2 D31g; > > >>> 1 >><> g11 ¼ d1 s0fD12 þ q1 D22 þ q2 D32g; 2 >>>> g02 ¼ d1 s0fD13 þ q1 D23 þ q2 D33g; > > >>> 2 >:> g21 ¼ d1 s0fD14 þ q1 D24 þ q2 D34g; b2 ¼ 2Refc1ð0Þg; Imfc1ð0Þg þ l2ImfdkdðsskÞg; Theorem 5.1 l2 determines the direction of the Hopf bifurcation: if l2 [ ðl2\0Þ, then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for s2 [ s20 ðs2\s20 Þ; b2 determines the stability of bifurcating periodic solutions: the bifurcating periodic solutions are orbitally asymptotically stable (unstable) if b2\0ðb2 [ 0Þ; and T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2 [ 0ðT2\0Þ: Numerical simulation Impact of nonlinear harvesting The primary interest here is to find the impact of quadratic harvesting of the fish population on the stability of the given Plankton–fish system. We obtain numerically that when fish population is harvested at low rate (E¼ 0:0001), given system (1) exhibits chaotic dynamic with initial data pð0Þ ¼ 10, z 0 ð Þ ¼ 5 and f ð0Þ ¼ 3 as shown in Fig. 1. When, we further increase the value of harvesting parameter E from 0.0001 to 0.09, Figs. 2 and 3 and 4 depicts that stability exchange takes place by changing the state of the system from chaos to limit cycle and finally to stable focus. In ecological sense, it means that when the fish population is maintained between E 2 ð0:0001; 0:1Þ, it has stabilizing effect on the phytoplankton (TPP)-zooPlankton–fish dynamic which is well agreed with the experiment findings of Carpenter et al. (1987), Leavitt and Findlay (1994), Pace et al. (1999), Strock et al. (2013), Attayde et al. (2010). Figure 5 depicts the existence of limit cycles around R at E ¼ 0:15 and E ¼ 1:0 and ecologically it implies the prevalence of bloom like situations in plankton ecosystem. These results indicate that the enhancement in harvesting of fish population induces instability in the plankton system. Moreover, in case of excessive harvesting, the entire fish population will be elimination and phytoplankton-zooplankton population coexist with the existence of limit cycles as depicted in Fig. 6. The numerical results predicts some estimation of the harvesting parameter that will be helpful for regulating agencies to formulate some harvesting policies of fisheries for coexistence of population in ecology. The equilibrium level of all population at different harvesting rates along with nature of stability around R is given in Table 1. litno20 a u op15 P Fig. 2 Change in stability from chaos to limit cycle at E ¼ 0:01 It is already well established that toxin has stabilizing effect on the tri-trophic interaction system (Upadhyay and Chattopadhyay 2005; Upadhyay et al. 2007; Upadhyay and Rao 2009). But, in this research, we have shown that toxin is not the only mechanism of inducing stability in the Plankton–fish system rather harvesting the generalist Fig. 1 System exhibits chaotic oscillation around R at E ¼ 0:0001 Fig. 3 Change in stability from limit cycle to a stable focus at E ¼ 0:09 Fig. 4 Time series graph shows the converges of solution trajectories to R at E ¼ 0:1 predator (fish) at appropriate level of harvesting can induce stability and trigger new mechanism of controlling complexity such as chaos in plankton system. Thus, our findings show that fish may have stabilizing effects if it can be harvested in the proposed harvesting range. Delay induced instability In the second part of this paper, we have considered the digestion delay in zooPlankton–fish interaction. We consider the given system in stable state by taking E ¼ 0:09 and keeping other parameters unchanged. From the sign of p ¼ 3:4729, q ¼ 3:2571 and r ¼ 0:0808, Eq. (22) has two purely imaginary roots ir0 with rþ ¼ 0:1612 and r ¼ 0:0356. Now from (23) the critical values of parameter delay at which stability exchange take place are given by s0þ ¼ 1:013 and s0 ¼ 49:981 such that R remain stable for s 2 ½0; 1:0Þ (see Fig. 7) litn20 o a u op15 P 200 400 600 800 1000 1200 1400 1600 1800 2000 Time and is unstable when s 2 ð1; 49:981Þ. Thus numerical simulation shows that the interior equilibrium converges in the range 0\s\s0þ and it loses stability when s passes through its critical value s0 and a Hopf bifurcation occur which is depicted in Fig. 8. When the parameter s varies further from critical value s0, it is observed that periodic oscillation appeared around E till s lies in the range s0þ\s\s0 (see Fig. 9) and regains its stability for s [ s0 (see Figs. 10, 11). The equilibrium point R remains locally asymptotically stable whenever the delay parameter lies in the range s0 \s\s1þ implies R again switches from stability to instability as s passes through s ¼ s1þ. Thus it is seen that a finite number of stability switches for the given system occurs when delay varies continuously. Thus it can be seen that there exist intervals of stability and instability as delay parameter varies through its critical value. Ecologically, it signifies the appearance and disappearance of blooms in plankton ecosystem. Finally the stability determining quantities for Hopfbifurcating periodic solutions at s ¼ s0 are given by c1ð0Þ ¼ 1:6307e þ 009 6:7954e þ 009i, l2 ¼ 1:5409e þ011, b2 ¼ 3:2614e þ 009, T2 ¼ 3:7051e þ 009. Hence, we conclude that the Hopf-bifurcation is subcritical in nature as well as the bifurcating periodic orbits are unstable and decreases as delay increases further through its critical value ‘s ¼ s0’. Discussion and conclusion The consequences of toxin liberation by phytoplankton are of great interest due to its possible detrimental effect on fisheries. The role of top predator on the stability of tritrophic food chain has been remained an interesting area of research. In this paper the Plankton–fish interaction were Fig. 5 Limit cycles oscillations of solution curves around R at E ¼ 0:15 and E ¼ 1:0 iton20 a l u p o P15 non linear (quadratic) harvesting. We have analyzed the given system with various ranges of harvesting parameter keeping other parametric values fixed. It is found numerically that the chaotic oscillation has occurred in the range 0:0001 E\0:009 and when values of E is taken in 0:009 E 0:01, these oscillations disappeared with the existence of limit cycle. Again, for 0:09 E 0:1, the system converges to a stable focus from limit cycle which shows that the harvesting of fish population induces the stability in the plankton food chain and hence stimulate a possible mechanism for controlling blooms in plankton ecosystem. Our findings has shown that introduction of nonlinear harvesting of fish population in Plankton–fish interaction can be utilized as a controlling agent. In the second part of the model system we have assumed that the interaction of species in real life are not instantaneous rather it is delayed by time known as time delay. In this paper we have assumed that the zooPlankton–fish interaction are delayed by some time lag known as digestion delay. The main interest of this part is to observe the effect of this delay on the stability of the given system brought by stocking of the fish population. It is found that the given system remained asymptotically stable for delay n20 o ilt a u p oP15 ino 20 lt a u p o P15 Fig. 6 Time series graph shows that phytoplankton-zooplankton coexist but fish population almost dies out at E ¼ 30 studied in the context of obtaining the impact of quadratic harvesting and time delay. We have assumed that the top predator in the tri-trophic food chain is harvested through Table 1 Impact of Non-linear Harvesting on stability of given system around interior equilibrium R Rate of harvesting, E Nature of system System oscillate chaotically around R System oscillate around R ¼(29,3,4) System converge to R ¼(15.2,12.4,3.1) System converge to R ¼(14.7,2.6,2.8) system oscillate around R ¼(14.7,12.6,2.8) system oscillate around R ¼(14,16,1.6) system oscillate around R ¼(4,7,0.00) in pz-plane Fig. 7 Stability of interior equilibrium R at s ¼ 0:1 1000 1500 2000 2500 3000 3500 4000 4500 5000 Fig. 8 Hopf bifurcation in the form of limit cycle oscillation at s ¼ 1:0 in the range 0 s so (see Fig. 3). When delay ‘s’ passes through its critical value ‘s0’ it is obtained that a phenomenon of switching of stability arises i.e. there exist a sequence of two time delays s0þ ðk ¼ 0; 1; 2; . . .Þ and s0 ðk ¼ 0; 1; 2; . . .Þ under certain parametric restrictions for which the interior equilibrium point E is stable whenever s 2 ½0; s0þÞ [ ðs0 ; s1þÞ [ [ ðsm 1; smþÞ, and unstable when s 2 ½s0þ; s0 Þ [ ðs1þ; s1 Þ [ [ ðsmþ 1; sm 1Þ[ ðsmþ; 1Þ. The above switching behavior of the system is shown in numerical simulation performed in the above section and it seems to be new findings in delayed Plankton–fish interactions. Thus it can be concluded that digestion delay has destabilizing effect on the TPPzooPlankton–fish interaction. Finally it can be predicted h s i 3 F that induction of stability by harvesting of top predator in plankton food chain can be destabilized by digestion delay. Appendix A Let s ¼ sk þ l , uiðtÞ ¼ uiðstÞ and utðhÞ ¼ uðt þ hÞ for h 2 ½ 1; 0 and dropping the bars for simplification of notations, system (15) becomes a functional differential equation in C ¼ Cð½ 1; 0 ; R3Þ as u_ðtÞ ¼ LlðutÞ þ f ðl; utÞ; where uðtÞ ¼ ðu1ðtÞ; u2ðtÞÞT 2 R2 and Ll : C ! R3; R C ! R2 are given, respectively, by 2 h s i F 0 20 2 h s i F 0 20 Fig. 9 Limit cycle at s ¼ 40. Fig. 10 Switch of stability from limit cycle to stable focus around s ¼ 50 f ðl; /Þ ¼ ðsk þ lÞ 2 a200/12ð0Þ þ a110/ð0Þ/2ð0Þ By Riesz representation theorem, there exists a function /ðh; lÞ of bounded variation for h 2 ½ 1; 0 such that 2 h s i F 0 20 ion15 lt a u p o P 10 dgðh; lÞ/ðhÞ In fact, we can choose /ðh; lÞ ¼ ðsk þ lÞ½A11dðhÞ A22dðh þ 1Þ ; where d denote the Dirac delta function. For / 2 Cð½ 1; 0 ; R3Þ, define 8 d/ðhÞ < AðlÞ/ðhÞ ¼ dh : R 01 d1ðs; lÞ/ðsÞ h ¼ 0 RðlÞ/ðhÞ ¼ Then system (26) is equivalent to Solving the above equations, we get 0 b001a010 From the definition of A , we obtain Z 0 A ðlÞwðsÞ ¼ 1 dgðt;0Þwð tÞ ¼ A1T1wð0Þ þ A2T2wð 1Þ: Since qðhÞ is the eigenvector of A corresponding to -ir0sk, then we have A q ð0Þ ¼ ir0skq ðhÞ; where I is identity matrix of order 2, that is, u_ðtÞ ¼ AðlÞut þ RðlÞut; For w 2 C1ð½0;1 ;ðRÞ Þ, define and a bilinear inner product c001 þ c0001e ir0sk þ ir0 Z 0 Z h \wðsÞ;/ðhÞ[ ¼ wð0Þ/ð0Þ 1 n¼0 wðn hÞdgðhÞ/ðnÞdn: Solving the above equations, we get 0 Að0ÞqðhÞ ¼ ixskqðhÞ; It follow from the definition of A(0), Llð/Þ and gðh;lÞ that ½ðA11 þ A22e ir0s0Þ ir0I qð0Þ ¼ 0; where I is identity matrix of order 2, that is, 0 0a100 ir0 a010 @B b100 b010 ir0 b001 0 c010 þ c0100e ir0sk c001 þ c0001e ir0sk ir0 0 1 1 001 B@ a CA ¼ @B0CA; b 0 ¼ d1½1þaa0 þab0 ð1;aa0;b0Þheihr0skÞdgðhÞð1;a;bÞT 1 ¼ d1½1þaa0 þab0 þskb0ðac0010 þbc0001Þe ir0sk : In order to have ðq ;qÞ ¼ 1, we can choose d1 ¼ 1 C 1 þ aa0 þ bb0 þ b0skeix0skðac0010 þ bc0001Þ, A such that ðw;A/Þ ¼ ðA w;/Þ. or ir0ðq ;qÞ ¼ ðq ;AqÞ ¼ ðA q ;qÞ ¼ ð ir0q ;qÞ ¼ ir0 ðq ;qÞ i.e. ðq ;qÞ ¼ 0. Next we will compute the coordinates to describe the center manifold C0 at l ¼ 0. Let ut be the solution of (31) when l ¼ 0. Define zðtÞ ¼ \q ; ut [ ; Wðt; hÞ ¼ utðhÞ 2RefzðtÞqðhÞg: ð33Þ On the center manifold C0, we have, Wðt; hÞ ¼ WðzðtÞ; zðtÞ; hÞ, where D13 ¼ a200 þ a110a; D11 ¼ a200 þ a110a; D12 ¼ 2a200 þ 2ReðaÞa110; D14 ¼ ð2W210ð0Þ þ 4W111ð0ÞÞa200 þ a110ðaW210ð0Þ Comparing its coefficients with (35), we find ð35Þ 8 2 >>> g20 ¼ d1 s0fD11 þ q1 D21 þ q2 D31g; > > >>> 1 >>>< g11 ¼ d1 s0fD12 þ q1 D22 þ q2 D32g; Since, W20 and W11 are in g21, we still to compute them. Now from (31) and (33), we have z and z are local coordinates for center manifold C0 in the direction of q and q respectively. Here we consider only real solutions as W is real if ut is real. For the solution ut 2 C0 of (31), since l ¼ 0, we have This equation can be rewritten as z_ðtÞ ¼ ir0skz þ gðz; zÞ; From (33) and (34), we have þ ð1;q1;q2ÞTe ix0skhz þ :::: Now from (28) and (35), it follows that Substituting (40) into (39) and comparing the coefficients, we get ðAð0Þ 2ix0skIÞW20ðhÞ ¼ H20ðhÞ; Að0ÞW11ðhÞ ¼ H11ðhÞ; From (39) and for h 2 ½ 1;0Þ H20ðhÞ ¼ g20qðhÞ g02qðhÞ; H11ðhÞ ¼ g11qðhÞ g11qðhÞ; From the definition of A(0), (41) and (43), we obtain W_20ðhÞ ¼ 2ir0skW20ðhÞ þ g20qðhÞ þ g02qðhÞ; Solving it and for qðhÞ ¼ ð1;q1;q2ÞTeir0skh, we have Similarly, from (41) and (44) it follows that, ig11 qð0Þeir0skh þ r0sk ig11 qð0Þe ir0skh þ E2; W11ðhÞ ¼ r0sk ð42Þ From (33), we have ð43Þ Thus we can obtain, 2D11 3z2 2D12 3 ð44Þ f0 ¼ 2sk64D21 572 þ sk64D22 57zz þ 2 D11 3 H20ð0Þ ¼ g20qð0Þ g02qð0Þ þ 2sk64 D21 57 H11ð0Þ ¼ g11qð0Þ g11qð0Þ þ sk64 D22 57 Since ir0sk is the eigenvalue of A(0) corresponding to eigenvector q(0), then z2 z2 H20ðhÞ 2 þ H11ðhÞzz þ H02ðhÞ 2 þ 2 2 ¼ qð0Þ g20 z2 þ g11zz þ g02 z2 þ ... 2 2 qð0Þ g20 z2 þ g11zz þ g02 z2 þ :::: þ f0ðz;zÞ; Z 0 1 d1ðhÞW20ðhÞ ¼ 2ir0skW20ð0Þ H20ð0Þ; where 1ðhÞ ¼ 1ð0;hÞ. From (39), we know when h ¼ 0, d1ðhÞW11ðhÞ ¼ H11ð0Þ; Substituting (45) and (51) into (47), we find f2ir0skI Z 0 e2ir0skhd1ðhÞgE1 ¼ 2sk64D21 57 1 a010 0 2ir0 b010 b001 c010 c0010e 2ir0s0 c001 c0001e 2ir0s0 E1 ¼ 264 D21 57 E2 ¼ 46 D22 5746 b100 b2 ¼ 2Refc1ð0Þg; Finally, we can compute the following quantities: IfcðÞg þ lIfdkdðsskÞg Attayde JL , Van Nes EH , Araujo AIL , Corso G , Scheffer M ( 2010 ) Omnivory by planktivores stabilizes plankton dynamics, but May either promote or reduce algal biomass . Ecosystem 13 ( 3 ): 410 - 420 Benndorf J , Kneschke H , Kossatz K , Penz E ( 1984 ) Manipulation of the pelagical food web by stocking with predaceous fishes . Int Revue Ges Hydrobiol 69 : 407 - 428 Beretta E , Kuang Y ( 2002 ) Geometric stability switch criteria in delay differential systems with delay dependant parameters . SIAM J Math Anal 33 : 1144 - 65 Berryman AA , Millstein JA ( 1989 ) Are ecological systems chaotic and if not , why not? Trends Ecol Evol 4 : 26 - 28 Birkhoff G , Rota GS ( 1882 ) Ordinary differential equations . Ginn , Boston Buskey E , Hyatt C ( 1995 ) Effects of Texas (USA) brown tide alga on planktonic grazers . Mar Ecol Prog Ser 126 : 285 - 292 Carpenter SR , Kitchell JR , Hodgson J , Cochran P , Elser J , Elser M , Lodge DM , Kretchmer D , He X ( 1987 ) Regulation of lake Rai V , Upadhyay RK ( 2004 ) Chaotic population dynamics and biology of the top pradator . Chaos Solitons Fractals 21 : 1195 - 1204 Rinaldi S , Solidoro C ( 1998 ) Chaos and peak-to-peak dynamics in a plankton-fish model . Theor Popul Biol 54 : 62 - 77 Sarkar RR , Chattopadhyay J ( 2003 ) Occurence of planktonic blooms under environmental fluctuations and its possible control mechanism-mathematical models and experimental observations . J Theor Biol 224 : 501 - 516 Sarwardi S , Haque M , Mandal PK ( 2012 ) Ratio-dependent predatorprey model of interacting population with delay effect . Nonlinear Dyn 69 : 817 - 836 Scheffer M ( 1991 ) Should we expect strange attractors behind plankton dynamics and if so, should we bother ? J Plankton Res 13 : 1291 - 1305 Sharma A , Sharma AK , Agnihotri K ( 2014 ) The dynamic of plankton-nutrient interaction with discrete delay . Appl. Math. Comput . 231 : 503 - 515 Sharma A , Sharma AK , Agnihotri K ( 2015 ) Analysis of a toxin producing phytoplankton-zooplankton interaction with holling IV type scheme and time delay . Nonlinear Dyn Int J Nonlinear Dyn Chaos Eng 81 ( 1 ): 13 - 25 Strock KE , Saros JE , Simon KS , McGowan S , Kinnison MT ( 2013 ) Cascading effects of generalist fish introduction in oligotrophic lakes . Hydrobiologia 711 : 99 - 113 Upadhyay RK , Chattopadhyay J ( 2005 ) Chaos to order: role of toxin producing phytoplankton in aquatic systems . Nonlin Anal Model Control 10 : 383 - 96 Upadhyay RK , Naji RK , Kumari N ( 2007 ) Dynamical complexity in some ecological models: effect of toxin production by phytoplankton . Nonlinear Anal Model Control 12 ( 1 ): 123 - 138 Upadhyay RK , Rao V ( 2009 ) Short-term recurrent chaos and role of toxin producing phytoplankton (TPP) on chaotic dynamics in aquatic systems . Chaos Solitons Fractals 39 : 1550 - 1564


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs40808-016-0248-x.pdf

Amit Sharma, Anuj Kumar Sharma, Kulbhushan Agnihotri. Complex dynamic of plankton–fish interaction with quadratic harvesting and time delay, Modeling Earth Systems and Environment, 2016, 204, DOI: 10.1007/s40808-016-0248-x