- Email: [email protected]
Hydra effects in stable food chain models
BioSystems 185 (2019) 104018
Contents lists available at ScienceDirect
BioSystems journal homepage: www.elsevier.com/locate/biosystems
Hydra effects in stable food chain models a,b
Debprasad Pal a b c
c,⁎
, Bapan Ghosh , Tapan Kumar Kar
T
a
Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Botanic Garden, Howrah 711103, West Bengal, India Department of Mathematics, Bethune College, 181, Bidhan Sarani, Kolkata 700006, West Bengal, India Department of Mathematics, National Institute of Technology Meghalaya, Bijni Complex, Shillong 793003, Meghalaya, India
ARTICLE INFO
ABSTRACT
Keywords: Population dynamics Food chain Harvesting Hydra effect Dynamical system Stability
In this paper, we explore the occurrence of the hydra effect in food chains, a popular research theme in the current decade. The hydra effect, one of the paradoxical results in theoretical and applied ecology refers to the fact where increasing mortality rate on a population enhances its own stock. The main focus is to propose a dynamical system model of food chain showing a stable steady state and estimate the variation of stock of targeted species with increasing mortality. In our model, the per capita growth rate of any predator trophic level does not depend upon its density. The prey–predator model incorporating such a feature for predator growth is referred to as ‘pure predator system’ (see Sieber and Hilker (2012), J. Math. Biol. (2012) 64: 341–360, Journal of Mathematical Biology). Keeping the above feature in mind, we study a Rosenweig-MacArthur food chain model with logistic prey growth and Holling type II functional responses. It is shown that hydra effect at stable state appears on (a) prey in a four-trophic system, (b) first predator in a five-trophic system, and (c) prey and second predator in a six-trophic system. Xiao and Cao (2009) (Mathematical and Computer Modelling 50 (2009) 360–379) established that limit cycle may be observed due to harvesting in a system with the ratio-dependent prey-predator system (example of a non “non-pure predator system”). Therefore, if harvesting causes instability on some range of mortality rate, the hydra effect cannot occur at a stable state. Some results show that the unique stable steady state in our model remains stable under harvesting of either trophic level. As a whole, our investigations have some contribution in understanding population interactions, fishery management and biological pest control tactic.
1. Introduction
hydra effect in more details and claimed that the hydra effect may also occur at a stable fixed point in unstructured population models. Interested readers should go through his contribution to acquire adequate knowledge about hydra effect and the conditions under which the hydra effect is more likely to happen. With the progress in understanding the hydra effect in population dynamics, there is some evidence regarding the existence of the hydra effect where species coexist in the form of oscillation. In the context of pioneering work of Sieber and Hilker (2012) by the investigation of Gause-type predator–prey systems with Holling type II and III functional responses, the time-average density of predator increases with the mortality rate when the prey–predator communities are evolved in a cyclic mode. Later, Ghosh et al. (2014) found that both the mean density (temporal average of biomass over a stable limit cycle) and the unstable equilibrium abundance of the predator become larger when fishing strength on the predator increases. Several prey–predator systems with stage-structured for prey were investigated by Abrams and Quince (2005) and showed that predator
The ecological interactions among populations are very complex and can lead to many paradoxical results such as: hydra effect (Matsuda and Abrams, 2004; Sieber and Hilker, 2012; Abrams, 2015), species enrichment (Abrams and Roth, 1994; Weide et al., 2019), species extinction in marine protected areas (Takashina et al., 2012), etc. Ricker (1954) was the first to discover that increase of mortality may increase the mean density of a population through the study of a discrete singlespecies model. Unfortunately, no much attention is paid to further research on this topic for a long time. The idea became more popular due to Abrams and Matsuda (2005), Matsuda and Abrams (2004), when they termed the event as hydra effect. Schröder et al. (2014) demonstrated that the hydra effect is a common feature in a stage-structured community in the presence of cyclic dynamics of the species. They believed that mere population increased due to higher mortality in an unstructured system may not be considered as hydra effect. In response to Schröder et al. (2014), Abrams (2015) reviewed the definition of
⁎
Corresponding author. E-mail address: [email protected] (B. Ghosh).
https://doi.org/10.1016/j.biosystems.2019.104018 Received 13 December 2018; Received in revised form 7 August 2019; Accepted 13 August 2019 Available online 19 August 2019 0303-2647/ © 2019 Elsevier B.V. All rights reserved.
BioSystems 185 (2019) 104018
D. Pal, et al.
biomass increases at stable state when its mortality rate was increased. Recently, Cortez and Abrams (2016) investigated different systems: (i) a prey–predator system, (ii) One predator and two-prey system and (iii) A simple food chain model with three species. In all the cases, they proved the existence of hydra effect at a stable steady state. In their models, the positive density-dependent effect on the hydra effect species is incorporated. A prey–predator system with generalist predator and stagestructure for both the populations was developed by Costa et al. (2017), and proved that augmentation (population introduction or recruitment into a system) of adult predator causes an increase in stock for the adult prey. Very recently, Costa and dos Anjos (2018) studied a predator–prey model with Allee effect and mutual interference among predators, and discover the existence of double hydra effect at stable states. We have already mentioned that the work of Ricker (1954) was the first to show the positive impact of mortality in a discrete model. Liz and Ruiz-Herrera (2012) proposed a single-species discrete-time model with harvesting and showed that both the stable and unstable equilibria achieve higher value while harvesting rate is increased. A study on a stage-structure single-species model in discrete time with harvesting showed that harvesting the adult may enhance its steady state density, which establishes the existence of hydra effect (Neverova et al., 2018). Based on a simultaneous work Weide et al. (2019) concluded that a hydra effect occurs in a discrete type prey–predator system as the mean predator population size increases for increasing its mortality. In the above literature survey, we found that the hydra effect appears at a stable steady state when the continuous system has stagestructure (Costa et al., 2017), the system has one-prey two-predator (Cortez and Abrams, 2016) or the system involves an Allee effect in predator (Costa and dos Anjos, 2018). A food chain with a very different per-capita growth rate can also possess a hydra effect (Cortez and Abrams, 2016). However, well known prey–predator systems such as Lotka-Volterra or Rosenzweig-MacArthur (RM) type models do not exhibit any hydra effect at a stable state. In fact, no hydra effect is observed in a food chain with three or four trophic levels when predation involves only linear functional responses (Ghosh and Kar, 2013). Even though Cortez and Abrams (2016) showed the existence of the hydra effect in a tri-trophic food chain system with linear interaction, but they included the nonlinear mortality terms for the middle and toppredator. Thus from the above discussion it seems to us that the existence of the hydra effect in the paper by Cortez and Abrams (2016) is due to the presence of nonlinear mortality terms. Moreover, the existence of the hydra effect was not established in a linear food chain model with arbitrary trophic levels (Legović et al., 2010) and RM type tri-trophic food chain model (Ghosh et al., 2018). Ghosh et al. (2018) did not prove analytically whether no hydra effect is observed at a stable state in the tri-trophic food chain. Thus the above literature survey motivates us to examine the existence of hydra effect in RM type higher trophic levels food chain with Holling type functional response. Due to the presence of Holling type functional response, such models with more than three trophic levels also experience both oscillation and chaos (Abrams and Roth, 1994; Gross et al., 2005). Also, we intend to examine the impact of harvesting on stability exchanges in the long food chain system numerically. The rest part of the paper is organized as follows: In Section 2 we first propose a food chain model with six trophic levels. The mathematical theory, methods and simulation techniques to study the model are also briefly explained therein. The positive impacts on biomass due to mortality are analyzed in Section 3. The discussion of obtained results and the possible applications of the same are presented in the last section.
and Holling type II functional response and found several distinct dynamics. Sieber and Hilker (2012) studied the RM model and observed the hydra effect accompanying unstable equilibrium. Analytical results through graphical representation on the RM model ensured that no hydra effect occurs with stable equilibrium (Ghosh et al., 2014). On the other hand, Ghosh et al. (2018) could not identify any hydra effect on RM type tri-trophic food chain at a stable state through numerical simulation with particular parameter sets, which does not produce any generic outcome. Therefore, a stronger result is necessary if the hydra effect occurs in such a food chain. Also, to the best of the authors’ knowledge, no attempt is found on four or higher trophic food chain, to examine the existence of hydra effect at stable state. Thus a novel question arises whether the hydra effect, at a stable state, occurs in a long food chain. For seeking an answer to this question, in the present article, we will investigate food chains up to six trophic levels. Such high numbers of trophic levels have been reported for many real ecosystems. Some authors have studied longer food chain for different purposes. A general food chain models with an arbitrary trophic level is investigated to estimate the species biomass and yield by Legović et al. (2010), Wollrab et al. (2012). A statistical model is studied by Pauly et al. (2002) to estimate stock with four trophic level fish species (peruvian anchoveta, pelagic fish, groundfish and invertebrates). Wilen and Wilen (2012) also proposed a typical food chain with four trophic levels and concluded that fishing can down the food chain, but it is not the only cause. Pauly and Watson (2003) investigated a six-trophic food chain. However, none of these articles studied the occurrence of hydra effect in a larger food chain. We first consider an RM type resourceconsumer model with six trophic levels. This food chain contains only one prey and five specialist predator trophic levels. The model reads as follows:
x dx = rx(1 ) dt K dy1 1 xy1 = 1 dt h1 + x dy2 2y y = 2 1 2 dt h2 + y1 dy3 3y y = 3 2 3 dt h3 + y2 dy4 4 y3 y4 = 4 dt h 4 + y4 dy5 5y y = 5 4 5 dt h5 + y4
1 xy1 q0 e0 x h1 + x 2 y1 y2 m1 y1 q1 e1 y1 h2 + y1 3 y2 y3 m2 y2 q2 e2 y2 h3 + y3 4 y3 y4 m3 y3 q3 e3 y3 h4 + y4 5 y4 y5 m4 y4 q4 e4 y4 h5 + y4
m5 y5
q5 e5 y5 ,
where x and yi (i = 1, 2,…, 5) are the biomass of the prey and ith predator, respectively at any time t. Here r and K are respectively the intrinsic growth rate and the environmental carrying capacity of the prey species. αi and βi are, respectively, the attack rate and conversion coefficient due to predation by the ith predator to its adjacent lower trophic level. The parameter hi is the half-saturation constants linked with the ith predator. mi is the density-independent specific mortality rate of the ith predator, qi (resp . ei) is the catchability coefficient (resp. fishing effort) for harvesting the ith predator while q0 and e0 represent the same for the prey populations. Both the natural mortality and fishing mortality proportionally depend on the population biomass. Thus, the model could be simplified by adding the mortality and harvesting terms together for any trophic level. However, the specific mortality rate is a demographic parameter whereas harvesting effort is a control parameter (independent parameter) which can be adjusted by a regulatory agency. Thus, we prefer to keep both the terms explicitly as this framework would allow us to analyze the system dynamics without harvesting. When harvesting is introduced, all the measurements such as species biomass and yield are functions of effort. Hence, the variation of species stock is studied with respect to the harvesting effort. We are interested to know whether there is any positive effect on
2. Model and methodology A large number of articles are devoted to studying RosenzweigMacArthur (RM) type prey–predator model with logistic prey growth 2
BioSystems 185 (2019) 104018
D. Pal, et al.
stable biomass of the target species (for fishing) with increasing effort on any trophic level. In order to determine the coexisting equilibrium points (x *, y1* , y2* , y3* , y4* , y5* ) (represent the population biomass) for the model (1), we solve the following nonlinear system:
x* 1 y1* q0 e0 ) h1 + x * K 2 y2* 1 1x* (m1 + q1 e1) h2 + y1* h1 + x *
r (1
2 2 y1*
3 y3*
h2 + y1* 3 3 y2*
4 y4* 5 y5*
h4 + y3*
(m5 + q5 e5)
h5 + y4*
(i) The food chain possesses a unique equilibrium under certain condition. (ii) The coexisting equilibrium is asymptotically stable without harvesting.
(m4 + q4 e4 ) = 0,
h5 + y4*
5 5 y4*
= 0,
(m3 + q3 e3) = 0,
h4 + y3*
4 4 y3*
= 0,
(m2 + q2 e2) = 0,
h3 + y2*
h3 + y2*
To achieve the stability condition we use Routh-Hurwitz criteria by computing the determinant of Hurwitz matrices. Asymptotic stability is ensured when the determinants of the all the Hurwitz matrices are positive. The details theory for the same is available in Gantmacher (2000). It can be seen later that the results for food chains with different trophic levels are distinct. Hence, we present the conditions for different food chains separately in subsequent analysis. We will show that these stability conditions are very much useful to calculate the rate of change of equilibrium stock due to harvesting intensity. Two important assumptions, in order to analyze the system analytically, are:
= 0.
(2)
Thus we compute the equilibrium value of a species increases/decreases with the effort to verify the existence of the hydra effect. If the system is asymptotically stable, then from the first principal minor of the Jacobian matrix a11 < 0. The other entries satisfy aij < 0 if i < j and aij > 0 otherwise. Without looking at the explicit form of aij, we will be able to know the sign of these quantities quickly. In particular, Ghosh et al. (2018) mathematically shown that a11 is negative once the positive equilibrium exists. In addition, we make use of the effective numerical scheme to calculate the equilibrium value of the stock. The stability behavior of the equilibrium is to be identified from the Jacobian matrix. Instead of verifying several conditions to due Routh-Hurwitz criteria, we directly compute the eigenvalue of the Jacobian matrix evaluated at the equilibrium. When effort is introduced into the model, we compute the equilibrium and corresponding eigenvalue for each value of effort. Instability can occur at equilibrium if one of the eigenvalues achieves positive real part. Finally, we will plot the stock and eigenvalues with respect to the effort. In all our simulation we fix the catchability coefficient as 0.05 for harvesting either species. An explanatory Matlab code for generating the figure for the four-trophic system is provided as a Supplementary material for interested readers.
Because of the long food chain, it is difficult to determine the form of the equilibrium in term of parameters explicitly. The number of coexisting equilibrium and existence conditions for a tri-trophic food chain are presented by Ghosh et al. (2018). Shortly, we would see that our analysis can be performed without finding the explicit form of the equilibrium in term of demographic parameters and effort. In the next step we establish the stability condition of the model at the equilibrium. The Jacobian matrix at (x *, y1* , y2* , y3* , y4* , y5* ) corresponding to the model is
a11 a12 0 0 0 0 a21 a22 a23 0 0 0 0 a32 a33 a34 0 0 = , 0 0 a43 a44 a45 0 0 0 0 a54 a55 a56 0 0 0 0 a65 0
A(x*, y1*, y2*, y3*, y4*, y5*)
where
a11 = a21 =
1 x *y1*
1 1 h1 y1*
2 y1*
a43 =
h2 + y1*
a65 =
,
a32 =
,
(h3 + y2* )2
3 3 h3 y3* , (h3 + y2*) 2 4 y3*
a45 =
a55 =
a22 =
3 y2* y3*
h4 + y3*
a34 =
a44 =
,
a54 =
,
a56 =
5 y4* y5*
(h5 + y4* ) 2
1x* , h1 + x *
a12 =
,
(h1 + x *)2
a23 =
a33 =
rx * , K
(h1 + x *) 2
2 y1* y2*
,
(h2 + y1* ) 2 2 2 h2 y2*
(h2 + y1* )2 3 y2*
h3 + y2*
In order to find the possibility of hydra effect, we systematically examine the food chains with three, four, five and six trophic levels both analytically and numerically.
,
, 3.1. Food chain with three trophic levels
4 y3* y4*
(h4 + y3* )2 4 4 h 4 y4*
(h4 + y3* )2 5 y4*
h5 + y4*
3. Model analysis
,
Ghosh et al. (2018) have observed that stock of prey and top-predator decreases at a stable state due to harvesting the respective species. However, this outcome is based on simulation for a particular set of model parameters. Therefore, we first study the tri-trophic food chain analytically to verify if the hydra effect appears at a stable steady state. The unique (assumption) equilibrium state E (x *, y1* , y2* ) satisfies the following nonlinear system:
,
,
5 5 h5 y5* . (h5 + y4* ) 2
1
The eigenvalues of the Jacobian matrix satisfy the equation 6
+ Q1
5
+ Q2
4
x* ) K * 1 x *y1
rx*(1
+ Q3
3
+ Q4
2
h1 + x * 2 2 y1* y2*
+ Q5 + Q6 = 0.
at the equilibrium, where the quantities Qi are easily be determined for a specific system from the principal minors of corresponding the Jacobian matrix at the equilibrium.
h2 + y1*
1 x *y1*
h1 + x * y 2 1* y2* h2 + y1*
q0 e0 x *
= 0,
(m1 + q1 e1) y1* = 0,
(m2 + q2 e2) y2*
= 0.
The Jacobian matrix at the equilibrium E (x *, y1* , y2* ) is 3
(3)
BioSystems 185 (2019) 104018
D. Pal, et al.
3.1.2. Harvesting first predator When predator is harvested then its equilibrium biomass is not a function of e1. Hence, no hydra effect is possible.
a11 a12 0 J = a21 a22 a23 , 0 a32 0 where
a11 =
a22 =
1 x *y1*
(h1 +
rx * , K
x *) 2
2 y1* y2*
(h2 + y1* )2
,
1x*
a12 =
a23 =
h1 + x *
2 y1*
h2 + y1*
,
,
a32 =
a21 =
3.1.3. Harvesting top predator Again differentiating each of the equation with respect to e2 we obtain the matrix equation
1 1 h1 y1* , (h1 + x *)2
2 2 h2 y2*
(h2 + y1* )2
. a11 a12 0 a21 a22 a23 0 a32 0
The eigenvalues of the Jacobian matrix satisfy the equation 3
+ Q1
2
(4)
+ Q2 + Q3 = 0,
where
Q1 =
de2 dy *2
=
0 0 . q2 y2* (7)
de2 (a11 + a22),
Q2 =
a12 a21 + a11 a22
a23 a32 ,
Q3 = a11 a23 a32.
One can notice that the column vector of the right hand side of the matrix equation (7) has only one non-zero element and located in a row equal to the specific harvested trophic level. We use this concept through the paper without writing the matrix equation explicitly. Cramer's rule yields
We calculate the variation of stock under harvesting until the system is asymptotically stable. According to Routh-Hurwitz criteria, the following three conditions are to be maintained for stability: (i) Q1 > 0, (ii) Q3 > 0 and (iii) Q1Q2 − Q3 > 0. Now
Q1 Q2
dx* de2 dy1*
Q3 = =
a12 q2 y2* dx* = < 0, de2 a11 a32 dy1* q y* = 2 2 > 0, de2 a32 dy *2 (a12 a21 a11 a22) q2 y2* = = de2 a11 a23 a32
(a11 + a22)( a12 a21 + a11 a22 a23 a32) a11 a23 a32 (a11 + a22)(a11 a22 a12 a21) + a22 a23 a32.
As a11 < 0 and a22a23a32 < 0 and Q1 > 0, Q1Q2 − Q3 > 0, we found a11a22 − a12a21 > 0. Now we employ effort on individual trophic level successively to estimate the variation of stock with effort.
(a11 a22
a12 a21) q2 y2* Q3
.
dy*
Since a11a22 − a12a21 > 0, we obtain de2 < 0 . Thus hydra effect 2 does not occur for harvesting the top predator. A major conclusion is that the hydra effect does not appear, at a stable steady state, in the tri-trophic food chain.
3.1.1. Harvesting prey When prey species is harvested all catchability coefficient is zero except q0. Differentiating each equation of the system (3) with respect to e0 we can obtain the following equations:
dy * dx* + a12 1 q0 x * = 0, de0 de0 dy* dy * dx* a21 + a22 1 + a23 2 = 0, de 0 de 0 de 0 dy1* a32 = 0. de 0
3.2. Food chain with four trophic levels
a11
We now analyze the food chain with four trophic levels. Let E (x *, y1* , y2* , y3*) is the only positive equilibrium. Then the eigenvalues of the Jacobian matrix satisfy the equation (5)
4
By making use the Jacobian matrix, the above system of equations (5) can be written in matrix form as
a11 a12 0 a21 a22 a23 0 a32 0
dx* de0 dy1* de0 dy *2
=
q0 x * 0 . 0
dy*1
de 0
+ Q2
2
+ Q3 + Q4 = 0,
(8)
= (a11 + a22 + a33), = a12 a21 + a11 a22 a23 a32 + a11 a33 + a22 a33 a34 a43, = ( a11 a23 a32 a12 a21 a33 + a11 a22 a33 a11 a34 a43 a22 a34 a43), = (a12 a21 a11 a22) a34 a43.
According to Routh-Hurwitz criteria, the following three conditions are to be satisfied for obtaining negative real part of the eigenvalues: (i) Q1 > 0, (ii) Q4 > 0, (iii) Q1Q2 − Q3 > 0 and (iv) Q1 Q2 Q3 Q32 > Q12 Q4 . As a consequence of the last two inequalities, we obtain Q3 > 0. Moreover, Q4 > 0 leads to the fact that
(6)
We now solve the above system of non-homogeneous equations (6) by Cramer's rule. The last equation yields
3
where
Q1 Q2 Q3 Q4
de0
+ Q1
= 0 . This is expected as the equilibrium
a12 a21
dy
stock of the middle trophic level, which is found from dt2 = 0 , does not depend on the prey harvesting. Ghosh et al. (2018) proved that a11 is necessarily negative when the coexisting equilibrium exists. In fact, for asymptotic stability the first principal minor needs to be negative, and otherwise Q3 would be neq x* dx* gative. Hence the first equation leads de = a0 < 0 . It is now easy to 0 11 check from the second equation that the equilibrium stock of the top predator decreases with effort. Thus we conclude that the hydra effect does not occur when prey species is harvested in the tri-trophic food chain.
a11 a22 < 0.
(9)
Likewise, tri-trophic food chain we can easily obtain a non-homogeneous system to compute the rate of change of species biomass with effort. We are not explicitly providing the matrix equation for the individual case. 3.2.1. Harvesting prey When prey is harvested, we can obtain the following inequalities for stock variation of individual species: 4
BioSystems 185 (2019) 104018
D. Pal, et al.
matrix evaluated at the equilibrium is negative. Hence, the steady-state of the unharvested system is asymptotically stable. We would like to determine the variation of stable steady stock of all the populations when prey is harvested. It is calculated that top-predator persists when effort ranges from e0 = 0 to e0 = 1.97 and the equilibrium within this effort range is hyperbolic. The variation of stock of all the species and maximum real part of all the eigenvalues are shown in Fig. 1. The corresponding coexisting stock within this effort is stable as the real part of all the eigenvalues is negative.1 It suggests that the prey harvesting does not destabilize the system which was asymptotically stable before harvesting. The change in the stock of prey species with respect to the effort is shown in the top left panel of Fig. 1. We can observe that prey biomass, at a stable state, increases with effort. Therefore, a hydra effect is exhibited when prey species is harvested. It is also examined that harvesting any other species does not have any positive impact towards the increase of stock. Moreover, there is no change in stability behavior when the harvesting is considered on the other trophic level. Thus, harvesting has a more stabilizing effect under such a parameter configuration.
a22 q0 x * dx* = > 0, de0 a11 a22 a12 a21 dy1* a21 q0 x * = < 0, de 0 a11 a22 a12 a21 dy *2 = 0, de 0 dy3* a21 a32 q0 x * = < 0. de 0 (a11 a22 a12 a21) a34 It claims that the hydra effect always occurs at stable state when prey species is harvested. 3.2.2. Harvesting first predator In this case we can obtain the following expressions:
a12 q1 y1* dx* = > 0, de1 a12 a21 a11 a22 dy1* a11 q1 y1* = < 0, de1 a12 a21 + a11 a22 dy *2 = 0, de1 dy3* a11 a32 q1 y1* = < 0. de1 (a12 a21 a11 a22) a34
3.3. Food chain with five trophic levels In the four-trophic system, we have detected the hydra effect only when prey species was harvested. In this subsection, we examine whether the hydra effect can appear only on prey species or it may appear on predator trophic levels in a community with five species. The unique interior equilibrium E (x *, y1* , y2* , y3* , y4* ) satisfies the characteristic equation
Since the biomass of the first predator decreases with effort, no hydra effect appears when the first predator is harvested. 3.2.3. Harvesting second predator Biomass of the target species is independent of effort, no hydra effect is observed.
5
de3
=
(
(a12 a21
+ Q3
2
+ Q4 + Q5 = 0,
(a11 + a22 + a33 + a44),
Q2 =
a12 a21
a23 a32 + a22 a33
a34 a43 + (a22 + a33) a44 a45 a54 ,
( a11 a23 a32 + a11 a22 a33
a11 a34 a43
a22 a34 a43 + a11 a22 a44
a23 a32 a44 + a11 a33 a44 + a22 a33 a44 Q4 Q5
)
+ a33 q3 y3*
a34 a43 (a11 a23 a32 + a33 (a12 a21
3
Q1 =
Q3 =
Finally,
=
+ Q2
+ a11 (a22 + a33 + a44 )
a12 a23 q3 y3* dx* = > 0, de3 (a12 a21 a11 a22) a43 dy1* a11 a23 q3 y3* = < 0, de3 (a12 a21 a11 a22) a43 dy *2 q y* = 3 3 > 0. de3 a43
dy3*
4
where
3.2.4. Harvesting top predator For harvesting the top predator, we can obtain
a11 a23 a32 a12 a21 a11 a22
+ Q1
a11 a22)) q3 y3*
a11 a22) a34 a43
a12 a21 (a33 + a44 )
(a11 + a22 + a33) a45 a54), = (a23 a32 a22 a33) a45 a54 + a12 a21 (a34 a43 a33 a44 + a45 a54 ) a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 + (a22 + a33) a45 a54), = (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 a54 .
The Routh-Hurwitz criteria for stability are: (i) Q1 > 0, (ii) Q5 > 0, (iii) Q1Q2 − Q3 > 0, (iv) Q1 Q2 Q3 Q32 Q12 Q4 + Q1 Q5 > 0 and (v) Q1 Q2 Q3 Q4 Q32 Q4 Q12 Q42 Q1 Q22 Q5 + Q2 Q3 Q5 + 2Q1 Q4 Q5 Now Q5 > 0 reduces to
.
In this case, we are not able to prove analytically the sign of the numerator for the above expression. However, many simulations reveal that the biomass of the top species decreases with effort. In conclusion, we observe that the hydra effect certainly occurs at stable state when prey is harvested. Harvesting any other trophic level decreases its density. One of the assumptions in our analysis was that harvesting can be continued until the equilibrium is stable. We do not have much information when the system loses its stability or it remains stable under harvesting. We provide some additional information for the above queries numerically. Moreover, we wish to visualize the change in the stock of all the trophic levels due to harvesting by a simulation technique. Considering parameter set r = 0.4, K = 20, α1 = 0.65, α2 = 0.8, α3 = 0.85, β1 = 0.8, β2 = 0.25, β3 = 0.11, m1 = 0.16, m2 = 0.1, m3 = 0.04, h1 = 8, h2 = 10 and h3 = 18 leads unique positive equilibrium P(8.5712, 5.8274, 2.1557, 0.7976) in the unharvested system. It is observed that the real part of all the eigenvalues of the Jacobian
Q52 > 0 .
a11 a23 a32 + a12 a21 a33 a11 a22 a33 > 0, i. e., a11 a23 a32 a11 a22 a33 > 0. Since, a11 < 0, we can write
a23 a32
a22 a33 < 0.
(10)
The analytical results for harvesting prey species are provided in Appendix A.1. It shows that increasing effort on prey harvesting reduces its equilibrium stock. Therefore, the hydra effect does not appear. We now employ effort on the first predator to examine the existence of the hydra effect. 1 Unlike the case of a smooth curve representing the prey stock, Max ( Re( i )) is a non-smooth curve. It is not surprising as Max ( Re( i )) collects the real part from eigenvalues with varying effort. More clearly, Max ( Re( i )) is not computed from only one eigenvalue.
5
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 1. The change in biomass of all the trophic levels is explored when the prey species is exploited in the food chain with three predator trophic levels. Biomass of the prey, first predator, second predator and top-predator increases, decreases, remains constant and decreases, respectively with effort. It is clear that the equilibrium value for the second predator does not depend on effort as its pre capita growth does not depends on its density. Hydra effect occurs on prey species. The change in the maximum value of the real part of all four eigenvalues shows that the system is stable under harvesting.
Harvesting first predator The rate of change in stock for different trophic levels under the harvesting of first predator is calculated as follows:
could not detect any hydra effect under the harvesting of the top-predator. For explanatory, we select the parameter set as r = 0.2, K = 20, α1 = 1.2, α2 = 1.5, α3 = 1.7, α4 = 1.9, β1 = 0.92, β2 = 0.8, β3 = 0.7, β4 = 0.4, m1 = 0.1, m2 = 0.07, m3 = 0.04, m4 = 0.02, h1 = 50, h2 = 55, h3 = 60 and h4 = 65 to show the existence of hydra effect when first predator is harvested. A unique coexisting equilibrium (13.3180, 5.8682, 5.3649, 1.7568, 2.0263), which is asymptotically stable, is obtained in the unharvested food chain with five trophic levels. All the populations persist when the first predator is harvested with fishing effort lesser than 1.47. Fig. 2 shows the equilibrium biomass y1* of the first predator increases with effort on first predator. We have computed the real part of the eigenvalues and observe that maximum value of the real parts is negative for any effort within e1 = 0 and e1 = 1.47. Thus, the equilibrium under harvesting is asymptotically stable and results in a hydra effect at a stable steady state. As before we have verified that harvesting either trophic level does not destabilize the coexisting steady state.
a12 a33 q1 y1* dx* = < 0, de1 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy1* a11 a33 q1 y1* = > 0, de1 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy *2 a11 a32 q1 y1* = < 0, de1 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy3* = 0, de1 dy *4 a11 a32 a43 q1 y1* = < 0. de1 a11 a23 a32 a45 + a12 a21 a33 a45 a11 a22 a33 a45 The above expression shows that the equilibrium stock of the first predator increases when it is harvested with increasing effort. Thus, the hydra effect, at stable equilibrium, appears within the whole range of effort where all populations coexist. The mathematical calculations in estimating the stock for harvesting the second predator are postponed to Appendix A.2. The main key to understand the ecological result depends upon sign of (a12a21 − a11a22). Several numerical simulations could not detect any hydra effect. On the other hand, one can easily prove that the stable stock of the third predator is independent of effort when it is harvested. We have further represented the analytical expression in Appendix A.3 when the top-predator is harvested. However, the complex form of the expressions prevents to establish any concrete conclusion on the variation of stock with effort. All parameter sets, we have selected so far,
3.4. Food chain with six trophic levels So far we have been able to detect the hydra effect only on a single trophic level in the food chain with four or five species. Now we analyze a food chain with six trophic levels. The eigenvalues of the Jacobian matrix corresponding to the food chain satisfy the equation 6
+ Q1
where
6
5
+ Q2
4
+ Q3
3
+ Q4
2
+ Q5 + Q6 = 0,
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 2. Impacts on biomass at stable steady state for all the trophic levels are shown when the first predator harvested in a food chain with five trophic levels. Hydra effect is exhibited on the first predator. The bottom panel claims that all the eigenvalues has a negative real part. It also shows that harvesting does not destabilize the dynamics of the equilibrium.
Q1 = Q2 = Q3 =
(a11 + a22 + a33 + a44 + a55), a12 a21 a23 a32 a34 a43 + a33 a44 a45 a54 + (a33 + a44) a55 + a22 (a33 + a44 + a55) + a11 (a22 + a33 + a44 + a55) a56 a65 , ( a22 a34 a43 + a22 a33 a44
a22 a45 a54
Q6 = (a12 a21 ( a34 a43 + a33 a44 ) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a56 a65. According to Routh-Hurwitz criteria, the following conditions are to be satisfied for obtaining negative real part of the eigenvalues: (i) Q1 > 0, (ii) Q6 > 0, (iii) Q1Q2 − Q3 > 0, (iv) Q1 Q2 Q3 Q32 Q12 Q4 + Q1 Q5 > 0 , Q1 Q2 Q3 Q4 Q32 Q4 Q12 Q42 Q1 Q22 Q5 + Q2 Q3 Q5 + 2Q1 Q4 Q5 (v) Q52 + Q12 Q2 Q6 Q1 Q3 Q6 > 0 and (vi) Q1 Q2 Q3 Q4 Q5 Q32 Q4 Q5 Q12 Q42 Q5 Q1 Q22 Q52 + Q2 Q3 Q52 + 2Q1 Q4 Q52 Q53 Q1 Q2 Q32 Q6 + Q33 Q6 + Q12 Q3 Q4 Q6 + 2Q12 Q2 Q5 Q6
a33 a45 a54 + a22 a33 a55
a34 a43 a55 + a22 a44 a55 + a33 a44 a55
a23 a32 (a44 + a55)
a12 a21 (a33 + a44 + a55)
(a22 + a33 + a44) a56 a65 + a11 ( a23 a32
a34 a43 + a33 a44
+ a22 (a33 + a44 + a55) Q4 = a23 a32 a45 a54 + a22 a33 a44 a55
a45 a54 + (a33 + a44 ) a55
a56 a65)),
a22 a33 a45 a54
a22 a34 a43 a55
a11 (a33 a45 a54 + a34 a43 a55
3Q1 Q3 Q5 Q6 Q13 Q62 > 0. Now Q6 > 0 implies
a23 a32 a44 a55
a12 a21 (a33 a44
a33 a44 a55
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44 ) < 0. (11)
+ a23 a32 (a44 + a55) + a22 (a34 a43
a33 a44 + a45 a54
+ a44 ) a56 a65 + (a23 a32 + a34 a43
(a33 + a44) a55))
The left hand side quantity will appear in the denominator for all cases in the expression of rate of change of the stock due to harvesting effort. From Q6 > 0 we can also obtain
a11 (a22 + a33
a22 a33
(a22 + a33) a44) a56 a65 + a12 a21 (a34 a43
a33 a44 + a45 a54
(a33 + a44) a55 + a56 a65),
Q5 =
((a22 a34 a43 + a23 a32 a44
a12 a21 (a33 a44 i. e., (a12 a21
a22 a33 a44 ) a56 a65
Thus,
+ a12 a21 (a33 a45 a54 + a34 a43 a55 a33 a44 a55 + (a33 + a44 ) a56 a65) + a11 ((a34 a43
a34 a43) + a11 (a22 a34 a43 a22 a33 a44 ) < 0, a11 a22 )(a33 a44 a34 a43 ) < 0.
a12 a21
a33 a44) a56 a65
a11 a22 < 0.
(12)
Now we examine the impacts on stock levels when individual trophic level.
+ a23 a32 (a45 a54 a44 a55 + a56 a65) a22 (a33 a45 a54 + a34 a43 a55 a33 a44 a55 + (a33 + a44 ) a56 a65)))
3.4.1. Harvesting prey When prey is harvested, the rate of change in the stock of individual
and 7
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 3. The variation of the biomass of all the trophic levels is shown when prey species is harvested in the six-species model. Clearly, the hydra effect occurs on the prey species. The bottom panel represents the largest of real part among all eigenvalues.
species with effort e0 are presented as:
a12 a23 a44 q2 y2* dx* = > 0, de2 a12 a21 (a33 a 44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44) dy1* a11 a23 a44 q2 y2* = < 0, de2 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy 2* (a12 a21 a11 a22) a44 q2 y2* = de2 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )
(a22 a34 a 43 + a23 a32 a 44 a22 a33 a44) q0 x * dx* = > 0, de0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a 43 + a23 a32 a 44 a22 a33 a44) dy1* a21 (a34 a43 a33 a 44 ) q0 x * = < 0, de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44 ) * dy 2 a21 a32 a44 q0 x * = > 0, de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44 ) dy*3 a21 a32 a43 q0 x * = < 0, de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a 43 + a23 a32 a 44 a22 a33 a44) dy*4 = 0, de 0 dy 5* a21 a32 a43 a54 q0 x * = < 0. de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44 )) a56
> 0, (using Eq. (12)) dy*3 de2 dy*4 de2 dy*5 de2
=
(a12 a21 a12 a21 (a33 a 44
a11 a22) a 43 q2 y2*
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)
< 0,
= 0, =
(a12 a21 (a12 a21 (a33 a 44
a11 a22) a 43 a54 q2 y2*
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)) a56
< 0.
Our above analysis reveals that the hydra effect is ensured when the second predator is harvested. The mathematical calculations are provided in Appendix B.2 when third predator is subjected to harvested. In this situation, we face difficulty to understand the variation of the biomass of the third predator. Several numerical simulations do not show any increase in biomass on the third predator. Thus, it is unlike to experience the hydra effect on the third predator. However, we must admit that an affirmative answer is yet to investigate. It is again cleared from Appendix B.3 that stock of the forth predator is independent of e4. Hence, we cannot expect any hydra effect. When the top predator is removed, it is found that variation of biomass for both the third and top predator will be the same (See Appendix B.4). However, analytically we could not establish any affirmative statement. Simulations predict that no hydra effect is observed when the top species is harvested. For further explanation, in our six species food chain, parameters
The sign of the numerator in the above quantities can easily be determined by identifying the sign of the individual entries of the Jacobian matrix. It is cleared that prey biomass increases with effort and results in a hydra effect. In the case of harvesting the first predator, the hydra effect, at stable state, is not appeared as the biomass of the first predator follows a decreasing trend for all effort range. The mathematical derivations are detailed in Appendix B.1. 3.4.2. Harvesting second predator Here the variation of stock with effort e2 can be seen as:
8
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 4. In the six-species model, the biomass of all the species is presented when the second predator is harvested. The bottom panel claims that all the eigenvalues have a negative real part.
are chosen as r = 0.2, K = 30, α1 = 1.2, α2 = 1.5, α3 = 1.7, α4 = 1.9, α5 = 10, β1 = 0.9, β2 = 0.8, β3 = 0.7, β4 = 0.4, β5 = 0.1, m1 = 0.1, m2 = 0.07, m3 = 0.04, m4 = 0.02, m5 = 0.01, h1 = 40, h2 = 45, h3 = 50, h4 = 55 and h5 = 60. This set produces a stable steady state (7.9195, 5.8783, 2.6622, 2.1264, 0.6061, 0.05024). For increasing effort for prey species within the range (0, 0.475) increases the prey density (Fig. 3). The top-predator population decreases with effort at equilibrium (Fig. 3). As the top-predator biomass is positive within the said effort range, biomass of all the species are positive at equilibrium. The steady state is asymptotically stable as the maximum real part of all the eigenvalues are negative (bottom panel). Hence, the hydra effect appears at stable state when prey species is harvested. Computation reveals that the effort range for harvesting the second predator should lie within (0, 0.41) in order to maintain the coexistence of all the species. Fig. 4 shows that the biomass of the second predator is increased with effort. The biomass of the top-predator tends to zero as e2 → 0.41. Then all the eigenvalues have a negative real part in the aforesaid effort interval as can be seen from the figure. Hence, the hydra effect is observed at a stable state when the second predator is harvested. Unlike the case of the previous two models, we here able to identify the existence of the hydra effect in multiple trophic levels. We have also verified that harvesting either trophic level under this parameter restriction preserves the asymptotic local stability behavior of the coexisting equilibrium.
by Sieber and Hilker (2012) to prove the existence of a hydra effect. However, that effect at stable state was not detected either in a RM models (Sieber and Hilker, 2012) or in tri-trophic food chains (Ghosh et al., 2018). The purpose of this article was to verify the existence of the hydra effect in the food chains. We have proposed a linear food chain model with logistic prey growth and Holling type II functional response. We posed the question of whether the hydra effect, at a stable state, can be a special feature in “pure predator systems” with higher trophic levels. An affirmative answer is found through the study of food chains up to six trophic levels successively. Although prey-predator models with two species are investigated in establishing many ecological phenomenons, multiple trophic levels are common in a real ecosystem. For instance, four- (Pauly et al., 2002) and five-trophic (Pauly and Watson, 2003) system for fish dynamics are available in literature. Therefore, proposing and analyzing ecological systems with higher trophic levels system are justified. We observed that hydra effect, at stable state, occurs on (a) prey in a four-trophic system (Fig. 1), (b) first predator in a five-trophic system (Fig. 2) and (c) prey (Fig. 3) and second predator (Fig. 4) in a six-trophic system. To the best of authors’ knowledge, we are the first to establish the hydra effect in relatively simple food chain models involving unique stable fixed point. Therefore, hydra effect may be observed even in the absence of nonlinear mortality term among different predator trophic levels (Cortez and Abrams, 2016), stage structure of community (Abrams and Quince, 2005; Costa et al., 2017) or Allee effect among species (Costa and dos Anjos, 2018). The existence of the hydra effect as well as the variation of stable stock under harvesting effort is shown in Table 1, when an individual
4. Results and conclusion Gause type prey-predator systems under cyclic dynamics are studied 9
BioSystems 185 (2019) 104018
D. Pal, et al.
number of hydra effect species is ([n/2] − 1). Also, two consecutive trophic levels in our models cannot experience the hydra effect. These statements can be considered as conjectures. This study is important in fishery management including biological conservation for estimating stock and yield correctly. Most importantly, knowledge of the hydra effect might help to implement proper action in the pest control program. Because, an attempt for reducing pest or removing/culling arbitrary species, may lead to a higher density of pest populations (occurrence of hydra effect) in complex pest-predator (natural enemy) systems. In this case, farmers and gardeners should look for alternative action plans beforehand. In some situations, we could not derive results analytically about the hydra effect. However, our predictions are reported above. All the stability conditions are not utilized for estimating the variation of stock. The complete scenario will be cleared when the remaining RouthHurwitz criteria can be applied effectively. One can take the work forward to prove these issues. Thus, a general framework could be established analytically to identify the hydra affect species in a food chain with an arbitrary number of trophic levels. This is an interesting and challenging problem as a future work from mathematical and ecological viewpoints. Studying the unstructured food chain models, in the present contribution, is the first step and hence we intend to work on the same theme for stage-structured prey-predator models in a continuous and discrete-time framework. Further, the existence of the hydra effect in a spatially structured community could equally be interesting as a future perspective.
Table 1 The change of biomass of all the trophic levels is shown when a single trophic is targeted for harvesting. We also reported if the targeted species is a hydra affect species. Here the symbol ↑ (resp. ↓) indicates that the biomass of a particular trophic level is increased (resp. decreased) with effort, whereas no variation of biomass is represented by the symbol −. The observation based on numerical simulation is indicated by *. Model
Target species
x*
y* 1
y* 2
Tri-trophic
x(t) y1(t) y2(t)
↓ – ↓
– – ↑
↓ ↓ ↓
Four-trophic
x(t) y1(t) y2(t) y3(t)
↑ ↑ – ↑
↓ ↓ – ↓
– – – ↑
↓ ↓ ↓ *↓
Five-trophic
x(t) y1(t) y2(t) y3(t) y4(t)
↓ ↓ ↓ – ↓
↑ ↑ ↑ – ↑
↓ ↓ *↓ – *↓
– – – – ↑
↓ ↓ *↓ ↓ *↓
Six-trophic
x(t) y1(t) y2(t) y3(t) y4(t) y5(t)
↑ ↑ ↑ ↑ – ↑
↓ ↓ ↓ ↓ – ↓
↑ ↑ ↑ ↑ – ↑
↓ ↓ ↓ *↓ – *↓
– – – – ↑
y* 3
y* 4
y* 5
Hydra effect NO NO NO YES NO NO *NO NO YES *NO NO *NO
↓ ↓ ↓ *↓ ↓ *↓
YES NO YES *NO NO *NO
trophic level is harvested in different food chains. Interestingly, one can observe that the rule of biomass variation maintains a reverse order from one trophic level to another one when the hydra effect appears. This order does not match the conventional rule. Legović et al. (2010) generalized a rule for the change of biomass due to harvesting in the ecosystem by investigating a simple prey-predator model. Our study suggests that Legović et al. (2010)’ results may not hold for some specific systems in the presence of a hydra effect. From our experience, we predict that a food chain with even (resp. odd) number (n) of trophic levels, hydra effect occurs at odd (resp. even) trophic levels except possibly the top three species. Moreover, the
Conflicts of interest The authors declare no conflicts of interest. Acknowledgment We appreciate the substantial comments and suggestions received from both the eminent reviewers and the Editor, which certainly helped to improve the content of the manuscript.
Appendix A. Analysis of food chain with five trophic levels Here we present some of the mathematical results and their ecological relevances when an individual trophic level is harvested. A.1 Harvesting prey Prey harvesting generates
(a23 a32 a22 a33) q0 x * dx* = < 0, de0 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy1* a21 a33 q0 x * = > 0, de 0 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy *2 a21 a32 q0 x * < 0, = de 0 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy3* = 0, de 0 dy *4 a21 a32 a43 q0 x * = < 0. de 0 a11 a23 a32 a45 + a12 a21 a33 a45 a11 a22 a33 a45 Here prey harvesting reduces its own biomass, and hence no hydra effect is possible. A.2 Harvesting second predator The variation of biomass at the equilibrium of the species are determined from the followings:
10
BioSystems 185 (2019) 104018
D. Pal, et al.
a12 a23 q2 y2* dx* = < 0, de2 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy1* a11 a23 q2 y2* = > 0, de2 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy *2 (a12 a21 a11 a22) q2 y2* = , de2 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy3* = 0, de2 dy *4 ( a12 a21 + a11 a22) a43 q2 y2* = . de2 (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 We are not able to determine the sign of (a12a21 − a11a22). But numerically, we have observed that the sign is negative. Hence it is more likely that no hydra effect appears when the second predator is harvested. A.3 Harvesting top predator When top predator is harvested, we observe that
a12 a23 a34 q4 y4* dx* = < 0, de4 a11 a23 a32 a54 + a12 a21 a33 a54 a11 a22 a33 a54 dy1* a11 a23 a34 q4 y4* = > 0, de4 a11 a23 a32 a54 + a12 a21 a33 a54 a11 a22 a33 a54 dy *2 ( a12 a21 + a11 a22) a34 q4 y4* = , de4 (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a54 dy3* q y* = 4 4 > 0. de4 a54 Now
dy*4 de4
=
a22 a33 a44 )) q4 y4*
(a12 a21 ( a34 a43 + a33 a44 ) + a11 (a22 a34 a43 + a23 a32 a44 (a11 a23 a32 + a12 a21 a33
a11 a22 a33) a45 a54
.
Again the analytical result is not available to know if the biomass of the top predator increases for some effort range when it is harvested. Thus, we cannot arrive at any specific conclusion. However, many simulations show a decreasing trend in biomass of the target species. Appendix B. Analysis of food chain with six trophic levels We derive mathematical results in support of the statements provided in Section 3.4. B.1 Harvesting first predator In this case, we have
dx* = de1 dy1* = de1 dy *2 = de1 dy3* = de1 dy *4 = de1 dy5* = de1
a12 (a34 a43 a12 a21 (a33 a44
a33 a44 ) q1 y1*
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a11 (a34 a43 a33 a44 ) q1 y1*
a22 a33 a44 )
a12 a21 ( a34 a43 + a33 a44) + a11 (a22 a34 a43 + a23 a32 a44 a11 a32 a44 q1 y1*
a22 a33 a44)
a12 a21 (a33 a44 a12 a21 (a33 a44
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a11 a32 a43 q1 y1* a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)
> 0, < 0,
> 0,
a22 a33 a44 )
< 0,
0, a11 a32 a43 a54 q1 y1* (a12 a21 (a33 a44
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)) a56
< 0.
The above inequalities hold true as (a34a43 − a33a44) < 0 along with a11 < 0 and a12 < 0. Hence, the hydra effect is not appeared. B.2 Harvesting third predator When e3 is the control parameter, we can obtain
11
BioSystems 185 (2019) 104018
D. Pal, et al.
a12 a23 a34 q3 y3* dx* = > 0, de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy1* a11 a23 a34 q3 y3* = < 0, de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy *2 (a12 a21 a11 a22 ) a34 q3 y3* = > 0, (using Eq.(12)) de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy3* (a11 a23 a32 + a12 a21 a33 a11 a22 a33) q3 y3* , = de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44) dy *4 = 0, de3 dy5* (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a54 q3 y3* = . de3 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a56 We find difficult to understand the variation of biomass for the third predator. Numerical results do not show any increase in biomass for third predator and top predator biomass decreases with effort. It is unlike to happen hydra effect on the third predator. B.3 Harvesting fourth predator Since the biomass of this trophic level is independent of e4, we cannot expect any hydra effect. Besides, we would like to know the variation of the biomass of other trophic levels under harvesting. To estimate the variation we represent the following system of equations:
dy * dx* + a12 1 de4 de4 dy1* dx* a21 + a22 de4 de4 dy1* dy2* a32 + a33 de4 de4 dy* dy*2 a43 + a44 3 de4 de4 dy*3 dy *4 a54 + a55 de4 de4 dy *4 a65 de4 a11
= 0, + a23 + a34 + a45 + a56
dy *2
= 0,
de4 dy*3
= 0,
de4 dy *4 de4 dy5* de4
= 0, q4 y4* = 0, = 0.
dy*
As de 4 = 0 , the first four equations form a homogeneous system with coefficient matrix non-singular. Thus the sub-system with the first four 4 equations has a unique solution. As a result, there is no change in biomass for all the trophic levels below y4(t). Finally, the fifth equation reveals that top predator's biomass decreases with effort. Although this analysis is not provided explicitly for food chain with smaller trophic levels. However, the same analysis can be applied easily. B.4 Harvesting top predator When the top predator is removed, we have
a12 a23 a34 a45 q5 y5* dx* = > 0, de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a65 dy1* a11 a23 a34 a45 q5 y5* = < 0, de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44)) a65 dy *2 ( a12 a21 + a11 a22 ) a34 a45 q5 y5* = > 0, de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a65 dy3* (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 q5 y5* = , de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a65 dy *4 q5 y5* = > 0, de5 a65 dy5* de5
=
(
(a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 a54 a12 a21 (a33 a 44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44)
a56 a65
+ a55 ) q5 y5*
.
It is cleared that the variation of biomass for both the third and top predator will be the same. However, analytically we could not establish any affirmative statement. Simulations confirm that no hydra effect is observed when top species is harvested. Appendix C. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.biosystems.2019.104018.
12
BioSystems 185 (2019) 104018
D. Pal, et al.
References
Matsuda, H., Abrams, P.A., 2004. Effects of predator prey interactions and adaptive change on sustainable yield. Can. J. Fish. Aquat. Sci. 61 (2), 175–184. Neverova, G., Abakumov, A., Yarovenko, I., Frisman, E.Y., 2018. Mode change in the dynamics of exploited limited population with age structure. Nonlinear Dyn. 94 (2), 827–844. Pauly, D., Christensen, V., Guénette, S., Pitcher, T.J., Sumaila, U.R., Walters, C.J., Watson, R., Zeller, D., 2002. Towards sustainability in world fisheries. Nature 418 (6898), 689. Pauly, D., Watson, R., 2003. Counting the last fish. Sci. Am. 289 (1), 42–47. Ricker, W.E., 1954. Stock and recruitment. J. Fish. Board Can. 11 (5), 559–623. Schröder, A., van Leeuwen, A., Cameron, T.C., 2014. When less is more: positive population-level effects of mortality. Trends Ecol. Evol. 29 (11), 614–624. Sieber, M., Hilker, F.M., 2012. The hydra effect in predator-prey models. J. Math. Biol. 64 (1-2), 341–360. Takashina, N., Mougi, A., Iwasa, Y., 2012. Paradox of marine protected areas: suppression of fishing may cause species loss. Popul. Ecol. 54 (3), 475–485. Weide, V., Varriale, M.C., Hilker, F.M., 2019. Hydra effect and paradox of enrichment in discrete-time predator-prey models. Math. Biosci. 310, 120–127. Wilen, C.D., Wilen, J.E., 2012. Fishing down the food chain revisited: Modeling exploited trophic systems. Ecol. Econ. 79, 80–88. Wollrab, S., Diehl, S., De Roos, A.M., 2012. Simple rules describe bottom-up and topdown control in food webs with alternative energy pathways. Ecol. Lett. 15 (9), 935–946. Xiao, M., Cao, J., 2009. Hopf bifurcation and non-hyperbolic equilibrium in a ratio-dependent predator-prey model with linear harvesting rate: analysis and computation. Math. Comput. Model. 50 (3-4), 360–379.
Abrams, P., 2015. The hydra effect is no myth. New Sci. 226 (3023), 28–29. Abrams, P.A., Matsuda, H., 2005. The effect of adaptive change in the prey on the dynamics of an exploited predator population. Can. J. Fish. Aquat. Sci. 62 (4), 758–766. Abrams, P.A., Quince, C., 2005. The impact of mortality on predator population size and stability in systems with stage-structured prey. Theor. Popul. Biol. 68 (4), 253–266. Abrams, P.A., Roth, J.D., 1994. The effects of enrichment of three-species food chains with nonlinear functional responses. Ecology 75 (4), 1118–1130. Cortez, M.H., Abrams, P.A., 2016. Hydra effects in stable communities and their implications for system dynamics. Ecology 97 (5), 1135–1145. Costa, M.I.d.S., dos Anjos, L., 2018. Multiple hydra effect in a predator–prey model with allee effect and mutual interference in the predator. Ecol. Model. 373, 22–24. Costa, M.I.D.S., Esteves, P.V., Faria, L.D.B., dos Anjos, L., 2017. Prey dynamics under generalist predator culling in stage structured models. Math. Biosci. 285, 68–74. Gantmacher, F.R., 2000. The Theory of Matrices, vol. 131 American Mathematical Soc. Ghosh, B., Kar, T., 2013. Possible ecosystem impacts of applying maximum sustainable yield policy in food chain models. J. Theor. Biol. 329, 6–14. Ghosh, B., Kar, T., Legović, T., 2014. Relationship between exploitation, oscillation, msy and extinction. Math. Biosci. 256, 1–9. Ghosh, B., Pal, D., Legović, T., Kar, T., 2018. Harvesting induced stability and instability in a tri-trophic food chain. Math. Biosci. 304, 89–99. Gross, T., Ebenhöh, W., Feudel, U., 2005. Long food chains are in general chaotic. Oikos 109 (1), 135–144. Legović, T., Klanjšček, J., Geček, S., 2010. Maximum sustainable yield and species extinction in ecosystems. Ecol. Model. 221 (12), 1569–1574. Liz, E., Ruiz-Herrera, A., 2012. The hydra effect, bubbles, and chaos in a simple discrete population model with constant effort harvesting. J. Math. Biol. 65 (5), 997–1016.
13
Contents lists available at ScienceDirect
BioSystems journal homepage: www.elsevier.com/locate/biosystems
Hydra effects in stable food chain models a,b
Debprasad Pal a b c
c,⁎
, Bapan Ghosh , Tapan Kumar Kar
T
a
Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Botanic Garden, Howrah 711103, West Bengal, India Department of Mathematics, Bethune College, 181, Bidhan Sarani, Kolkata 700006, West Bengal, India Department of Mathematics, National Institute of Technology Meghalaya, Bijni Complex, Shillong 793003, Meghalaya, India
ARTICLE INFO
ABSTRACT
Keywords: Population dynamics Food chain Harvesting Hydra effect Dynamical system Stability
In this paper, we explore the occurrence of the hydra effect in food chains, a popular research theme in the current decade. The hydra effect, one of the paradoxical results in theoretical and applied ecology refers to the fact where increasing mortality rate on a population enhances its own stock. The main focus is to propose a dynamical system model of food chain showing a stable steady state and estimate the variation of stock of targeted species with increasing mortality. In our model, the per capita growth rate of any predator trophic level does not depend upon its density. The prey–predator model incorporating such a feature for predator growth is referred to as ‘pure predator system’ (see Sieber and Hilker (2012), J. Math. Biol. (2012) 64: 341–360, Journal of Mathematical Biology). Keeping the above feature in mind, we study a Rosenweig-MacArthur food chain model with logistic prey growth and Holling type II functional responses. It is shown that hydra effect at stable state appears on (a) prey in a four-trophic system, (b) first predator in a five-trophic system, and (c) prey and second predator in a six-trophic system. Xiao and Cao (2009) (Mathematical and Computer Modelling 50 (2009) 360–379) established that limit cycle may be observed due to harvesting in a system with the ratio-dependent prey-predator system (example of a non “non-pure predator system”). Therefore, if harvesting causes instability on some range of mortality rate, the hydra effect cannot occur at a stable state. Some results show that the unique stable steady state in our model remains stable under harvesting of either trophic level. As a whole, our investigations have some contribution in understanding population interactions, fishery management and biological pest control tactic.
1. Introduction
hydra effect in more details and claimed that the hydra effect may also occur at a stable fixed point in unstructured population models. Interested readers should go through his contribution to acquire adequate knowledge about hydra effect and the conditions under which the hydra effect is more likely to happen. With the progress in understanding the hydra effect in population dynamics, there is some evidence regarding the existence of the hydra effect where species coexist in the form of oscillation. In the context of pioneering work of Sieber and Hilker (2012) by the investigation of Gause-type predator–prey systems with Holling type II and III functional responses, the time-average density of predator increases with the mortality rate when the prey–predator communities are evolved in a cyclic mode. Later, Ghosh et al. (2014) found that both the mean density (temporal average of biomass over a stable limit cycle) and the unstable equilibrium abundance of the predator become larger when fishing strength on the predator increases. Several prey–predator systems with stage-structured for prey were investigated by Abrams and Quince (2005) and showed that predator
The ecological interactions among populations are very complex and can lead to many paradoxical results such as: hydra effect (Matsuda and Abrams, 2004; Sieber and Hilker, 2012; Abrams, 2015), species enrichment (Abrams and Roth, 1994; Weide et al., 2019), species extinction in marine protected areas (Takashina et al., 2012), etc. Ricker (1954) was the first to discover that increase of mortality may increase the mean density of a population through the study of a discrete singlespecies model. Unfortunately, no much attention is paid to further research on this topic for a long time. The idea became more popular due to Abrams and Matsuda (2005), Matsuda and Abrams (2004), when they termed the event as hydra effect. Schröder et al. (2014) demonstrated that the hydra effect is a common feature in a stage-structured community in the presence of cyclic dynamics of the species. They believed that mere population increased due to higher mortality in an unstructured system may not be considered as hydra effect. In response to Schröder et al. (2014), Abrams (2015) reviewed the definition of
⁎
Corresponding author. E-mail address: [email protected] (B. Ghosh).
https://doi.org/10.1016/j.biosystems.2019.104018 Received 13 December 2018; Received in revised form 7 August 2019; Accepted 13 August 2019 Available online 19 August 2019 0303-2647/ © 2019 Elsevier B.V. All rights reserved.
BioSystems 185 (2019) 104018
D. Pal, et al.
biomass increases at stable state when its mortality rate was increased. Recently, Cortez and Abrams (2016) investigated different systems: (i) a prey–predator system, (ii) One predator and two-prey system and (iii) A simple food chain model with three species. In all the cases, they proved the existence of hydra effect at a stable steady state. In their models, the positive density-dependent effect on the hydra effect species is incorporated. A prey–predator system with generalist predator and stagestructure for both the populations was developed by Costa et al. (2017), and proved that augmentation (population introduction or recruitment into a system) of adult predator causes an increase in stock for the adult prey. Very recently, Costa and dos Anjos (2018) studied a predator–prey model with Allee effect and mutual interference among predators, and discover the existence of double hydra effect at stable states. We have already mentioned that the work of Ricker (1954) was the first to show the positive impact of mortality in a discrete model. Liz and Ruiz-Herrera (2012) proposed a single-species discrete-time model with harvesting and showed that both the stable and unstable equilibria achieve higher value while harvesting rate is increased. A study on a stage-structure single-species model in discrete time with harvesting showed that harvesting the adult may enhance its steady state density, which establishes the existence of hydra effect (Neverova et al., 2018). Based on a simultaneous work Weide et al. (2019) concluded that a hydra effect occurs in a discrete type prey–predator system as the mean predator population size increases for increasing its mortality. In the above literature survey, we found that the hydra effect appears at a stable steady state when the continuous system has stagestructure (Costa et al., 2017), the system has one-prey two-predator (Cortez and Abrams, 2016) or the system involves an Allee effect in predator (Costa and dos Anjos, 2018). A food chain with a very different per-capita growth rate can also possess a hydra effect (Cortez and Abrams, 2016). However, well known prey–predator systems such as Lotka-Volterra or Rosenzweig-MacArthur (RM) type models do not exhibit any hydra effect at a stable state. In fact, no hydra effect is observed in a food chain with three or four trophic levels when predation involves only linear functional responses (Ghosh and Kar, 2013). Even though Cortez and Abrams (2016) showed the existence of the hydra effect in a tri-trophic food chain system with linear interaction, but they included the nonlinear mortality terms for the middle and toppredator. Thus from the above discussion it seems to us that the existence of the hydra effect in the paper by Cortez and Abrams (2016) is due to the presence of nonlinear mortality terms. Moreover, the existence of the hydra effect was not established in a linear food chain model with arbitrary trophic levels (Legović et al., 2010) and RM type tri-trophic food chain model (Ghosh et al., 2018). Ghosh et al. (2018) did not prove analytically whether no hydra effect is observed at a stable state in the tri-trophic food chain. Thus the above literature survey motivates us to examine the existence of hydra effect in RM type higher trophic levels food chain with Holling type functional response. Due to the presence of Holling type functional response, such models with more than three trophic levels also experience both oscillation and chaos (Abrams and Roth, 1994; Gross et al., 2005). Also, we intend to examine the impact of harvesting on stability exchanges in the long food chain system numerically. The rest part of the paper is organized as follows: In Section 2 we first propose a food chain model with six trophic levels. The mathematical theory, methods and simulation techniques to study the model are also briefly explained therein. The positive impacts on biomass due to mortality are analyzed in Section 3. The discussion of obtained results and the possible applications of the same are presented in the last section.
and Holling type II functional response and found several distinct dynamics. Sieber and Hilker (2012) studied the RM model and observed the hydra effect accompanying unstable equilibrium. Analytical results through graphical representation on the RM model ensured that no hydra effect occurs with stable equilibrium (Ghosh et al., 2014). On the other hand, Ghosh et al. (2018) could not identify any hydra effect on RM type tri-trophic food chain at a stable state through numerical simulation with particular parameter sets, which does not produce any generic outcome. Therefore, a stronger result is necessary if the hydra effect occurs in such a food chain. Also, to the best of the authors’ knowledge, no attempt is found on four or higher trophic food chain, to examine the existence of hydra effect at stable state. Thus a novel question arises whether the hydra effect, at a stable state, occurs in a long food chain. For seeking an answer to this question, in the present article, we will investigate food chains up to six trophic levels. Such high numbers of trophic levels have been reported for many real ecosystems. Some authors have studied longer food chain for different purposes. A general food chain models with an arbitrary trophic level is investigated to estimate the species biomass and yield by Legović et al. (2010), Wollrab et al. (2012). A statistical model is studied by Pauly et al. (2002) to estimate stock with four trophic level fish species (peruvian anchoveta, pelagic fish, groundfish and invertebrates). Wilen and Wilen (2012) also proposed a typical food chain with four trophic levels and concluded that fishing can down the food chain, but it is not the only cause. Pauly and Watson (2003) investigated a six-trophic food chain. However, none of these articles studied the occurrence of hydra effect in a larger food chain. We first consider an RM type resourceconsumer model with six trophic levels. This food chain contains only one prey and five specialist predator trophic levels. The model reads as follows:
x dx = rx(1 ) dt K dy1 1 xy1 = 1 dt h1 + x dy2 2y y = 2 1 2 dt h2 + y1 dy3 3y y = 3 2 3 dt h3 + y2 dy4 4 y3 y4 = 4 dt h 4 + y4 dy5 5y y = 5 4 5 dt h5 + y4
1 xy1 q0 e0 x h1 + x 2 y1 y2 m1 y1 q1 e1 y1 h2 + y1 3 y2 y3 m2 y2 q2 e2 y2 h3 + y3 4 y3 y4 m3 y3 q3 e3 y3 h4 + y4 5 y4 y5 m4 y4 q4 e4 y4 h5 + y4
m5 y5
q5 e5 y5 ,
where x and yi (i = 1, 2,…, 5) are the biomass of the prey and ith predator, respectively at any time t. Here r and K are respectively the intrinsic growth rate and the environmental carrying capacity of the prey species. αi and βi are, respectively, the attack rate and conversion coefficient due to predation by the ith predator to its adjacent lower trophic level. The parameter hi is the half-saturation constants linked with the ith predator. mi is the density-independent specific mortality rate of the ith predator, qi (resp . ei) is the catchability coefficient (resp. fishing effort) for harvesting the ith predator while q0 and e0 represent the same for the prey populations. Both the natural mortality and fishing mortality proportionally depend on the population biomass. Thus, the model could be simplified by adding the mortality and harvesting terms together for any trophic level. However, the specific mortality rate is a demographic parameter whereas harvesting effort is a control parameter (independent parameter) which can be adjusted by a regulatory agency. Thus, we prefer to keep both the terms explicitly as this framework would allow us to analyze the system dynamics without harvesting. When harvesting is introduced, all the measurements such as species biomass and yield are functions of effort. Hence, the variation of species stock is studied with respect to the harvesting effort. We are interested to know whether there is any positive effect on
2. Model and methodology A large number of articles are devoted to studying RosenzweigMacArthur (RM) type prey–predator model with logistic prey growth 2
BioSystems 185 (2019) 104018
D. Pal, et al.
stable biomass of the target species (for fishing) with increasing effort on any trophic level. In order to determine the coexisting equilibrium points (x *, y1* , y2* , y3* , y4* , y5* ) (represent the population biomass) for the model (1), we solve the following nonlinear system:
x* 1 y1* q0 e0 ) h1 + x * K 2 y2* 1 1x* (m1 + q1 e1) h2 + y1* h1 + x *
r (1
2 2 y1*
3 y3*
h2 + y1* 3 3 y2*
4 y4* 5 y5*
h4 + y3*
(m5 + q5 e5)
h5 + y4*
(i) The food chain possesses a unique equilibrium under certain condition. (ii) The coexisting equilibrium is asymptotically stable without harvesting.
(m4 + q4 e4 ) = 0,
h5 + y4*
5 5 y4*
= 0,
(m3 + q3 e3) = 0,
h4 + y3*
4 4 y3*
= 0,
(m2 + q2 e2) = 0,
h3 + y2*
h3 + y2*
To achieve the stability condition we use Routh-Hurwitz criteria by computing the determinant of Hurwitz matrices. Asymptotic stability is ensured when the determinants of the all the Hurwitz matrices are positive. The details theory for the same is available in Gantmacher (2000). It can be seen later that the results for food chains with different trophic levels are distinct. Hence, we present the conditions for different food chains separately in subsequent analysis. We will show that these stability conditions are very much useful to calculate the rate of change of equilibrium stock due to harvesting intensity. Two important assumptions, in order to analyze the system analytically, are:
= 0.
(2)
Thus we compute the equilibrium value of a species increases/decreases with the effort to verify the existence of the hydra effect. If the system is asymptotically stable, then from the first principal minor of the Jacobian matrix a11 < 0. The other entries satisfy aij < 0 if i < j and aij > 0 otherwise. Without looking at the explicit form of aij, we will be able to know the sign of these quantities quickly. In particular, Ghosh et al. (2018) mathematically shown that a11 is negative once the positive equilibrium exists. In addition, we make use of the effective numerical scheme to calculate the equilibrium value of the stock. The stability behavior of the equilibrium is to be identified from the Jacobian matrix. Instead of verifying several conditions to due Routh-Hurwitz criteria, we directly compute the eigenvalue of the Jacobian matrix evaluated at the equilibrium. When effort is introduced into the model, we compute the equilibrium and corresponding eigenvalue for each value of effort. Instability can occur at equilibrium if one of the eigenvalues achieves positive real part. Finally, we will plot the stock and eigenvalues with respect to the effort. In all our simulation we fix the catchability coefficient as 0.05 for harvesting either species. An explanatory Matlab code for generating the figure for the four-trophic system is provided as a Supplementary material for interested readers.
Because of the long food chain, it is difficult to determine the form of the equilibrium in term of parameters explicitly. The number of coexisting equilibrium and existence conditions for a tri-trophic food chain are presented by Ghosh et al. (2018). Shortly, we would see that our analysis can be performed without finding the explicit form of the equilibrium in term of demographic parameters and effort. In the next step we establish the stability condition of the model at the equilibrium. The Jacobian matrix at (x *, y1* , y2* , y3* , y4* , y5* ) corresponding to the model is
a11 a12 0 0 0 0 a21 a22 a23 0 0 0 0 a32 a33 a34 0 0 = , 0 0 a43 a44 a45 0 0 0 0 a54 a55 a56 0 0 0 0 a65 0
A(x*, y1*, y2*, y3*, y4*, y5*)
where
a11 = a21 =
1 x *y1*
1 1 h1 y1*
2 y1*
a43 =
h2 + y1*
a65 =
,
a32 =
,
(h3 + y2* )2
3 3 h3 y3* , (h3 + y2*) 2 4 y3*
a45 =
a55 =
a22 =
3 y2* y3*
h4 + y3*
a34 =
a44 =
,
a54 =
,
a56 =
5 y4* y5*
(h5 + y4* ) 2
1x* , h1 + x *
a12 =
,
(h1 + x *)2
a23 =
a33 =
rx * , K
(h1 + x *) 2
2 y1* y2*
,
(h2 + y1* ) 2 2 2 h2 y2*
(h2 + y1* )2 3 y2*
h3 + y2*
In order to find the possibility of hydra effect, we systematically examine the food chains with three, four, five and six trophic levels both analytically and numerically.
,
, 3.1. Food chain with three trophic levels
4 y3* y4*
(h4 + y3* )2 4 4 h 4 y4*
(h4 + y3* )2 5 y4*
h5 + y4*
3. Model analysis
,
Ghosh et al. (2018) have observed that stock of prey and top-predator decreases at a stable state due to harvesting the respective species. However, this outcome is based on simulation for a particular set of model parameters. Therefore, we first study the tri-trophic food chain analytically to verify if the hydra effect appears at a stable steady state. The unique (assumption) equilibrium state E (x *, y1* , y2* ) satisfies the following nonlinear system:
,
,
5 5 h5 y5* . (h5 + y4* ) 2
1
The eigenvalues of the Jacobian matrix satisfy the equation 6
+ Q1
5
+ Q2
4
x* ) K * 1 x *y1
rx*(1
+ Q3
3
+ Q4
2
h1 + x * 2 2 y1* y2*
+ Q5 + Q6 = 0.
at the equilibrium, where the quantities Qi are easily be determined for a specific system from the principal minors of corresponding the Jacobian matrix at the equilibrium.
h2 + y1*
1 x *y1*
h1 + x * y 2 1* y2* h2 + y1*
q0 e0 x *
= 0,
(m1 + q1 e1) y1* = 0,
(m2 + q2 e2) y2*
= 0.
The Jacobian matrix at the equilibrium E (x *, y1* , y2* ) is 3
(3)
BioSystems 185 (2019) 104018
D. Pal, et al.
3.1.2. Harvesting first predator When predator is harvested then its equilibrium biomass is not a function of e1. Hence, no hydra effect is possible.
a11 a12 0 J = a21 a22 a23 , 0 a32 0 where
a11 =
a22 =
1 x *y1*
(h1 +
rx * , K
x *) 2
2 y1* y2*
(h2 + y1* )2
,
1x*
a12 =
a23 =
h1 + x *
2 y1*
h2 + y1*
,
,
a32 =
a21 =
3.1.3. Harvesting top predator Again differentiating each of the equation with respect to e2 we obtain the matrix equation
1 1 h1 y1* , (h1 + x *)2
2 2 h2 y2*
(h2 + y1* )2
. a11 a12 0 a21 a22 a23 0 a32 0
The eigenvalues of the Jacobian matrix satisfy the equation 3
+ Q1
2
(4)
+ Q2 + Q3 = 0,
where
Q1 =
de2 dy *2
=
0 0 . q2 y2* (7)
de2 (a11 + a22),
Q2 =
a12 a21 + a11 a22
a23 a32 ,
Q3 = a11 a23 a32.
One can notice that the column vector of the right hand side of the matrix equation (7) has only one non-zero element and located in a row equal to the specific harvested trophic level. We use this concept through the paper without writing the matrix equation explicitly. Cramer's rule yields
We calculate the variation of stock under harvesting until the system is asymptotically stable. According to Routh-Hurwitz criteria, the following three conditions are to be maintained for stability: (i) Q1 > 0, (ii) Q3 > 0 and (iii) Q1Q2 − Q3 > 0. Now
Q1 Q2
dx* de2 dy1*
Q3 = =
a12 q2 y2* dx* = < 0, de2 a11 a32 dy1* q y* = 2 2 > 0, de2 a32 dy *2 (a12 a21 a11 a22) q2 y2* = = de2 a11 a23 a32
(a11 + a22)( a12 a21 + a11 a22 a23 a32) a11 a23 a32 (a11 + a22)(a11 a22 a12 a21) + a22 a23 a32.
As a11 < 0 and a22a23a32 < 0 and Q1 > 0, Q1Q2 − Q3 > 0, we found a11a22 − a12a21 > 0. Now we employ effort on individual trophic level successively to estimate the variation of stock with effort.
(a11 a22
a12 a21) q2 y2* Q3
.
dy*
Since a11a22 − a12a21 > 0, we obtain de2 < 0 . Thus hydra effect 2 does not occur for harvesting the top predator. A major conclusion is that the hydra effect does not appear, at a stable steady state, in the tri-trophic food chain.
3.1.1. Harvesting prey When prey species is harvested all catchability coefficient is zero except q0. Differentiating each equation of the system (3) with respect to e0 we can obtain the following equations:
dy * dx* + a12 1 q0 x * = 0, de0 de0 dy* dy * dx* a21 + a22 1 + a23 2 = 0, de 0 de 0 de 0 dy1* a32 = 0. de 0
3.2. Food chain with four trophic levels
a11
We now analyze the food chain with four trophic levels. Let E (x *, y1* , y2* , y3*) is the only positive equilibrium. Then the eigenvalues of the Jacobian matrix satisfy the equation (5)
4
By making use the Jacobian matrix, the above system of equations (5) can be written in matrix form as
a11 a12 0 a21 a22 a23 0 a32 0
dx* de0 dy1* de0 dy *2
=
q0 x * 0 . 0
dy*1
de 0
+ Q2
2
+ Q3 + Q4 = 0,
(8)
= (a11 + a22 + a33), = a12 a21 + a11 a22 a23 a32 + a11 a33 + a22 a33 a34 a43, = ( a11 a23 a32 a12 a21 a33 + a11 a22 a33 a11 a34 a43 a22 a34 a43), = (a12 a21 a11 a22) a34 a43.
According to Routh-Hurwitz criteria, the following three conditions are to be satisfied for obtaining negative real part of the eigenvalues: (i) Q1 > 0, (ii) Q4 > 0, (iii) Q1Q2 − Q3 > 0 and (iv) Q1 Q2 Q3 Q32 > Q12 Q4 . As a consequence of the last two inequalities, we obtain Q3 > 0. Moreover, Q4 > 0 leads to the fact that
(6)
We now solve the above system of non-homogeneous equations (6) by Cramer's rule. The last equation yields
3
where
Q1 Q2 Q3 Q4
de0
+ Q1
= 0 . This is expected as the equilibrium
a12 a21
dy
stock of the middle trophic level, which is found from dt2 = 0 , does not depend on the prey harvesting. Ghosh et al. (2018) proved that a11 is necessarily negative when the coexisting equilibrium exists. In fact, for asymptotic stability the first principal minor needs to be negative, and otherwise Q3 would be neq x* dx* gative. Hence the first equation leads de = a0 < 0 . It is now easy to 0 11 check from the second equation that the equilibrium stock of the top predator decreases with effort. Thus we conclude that the hydra effect does not occur when prey species is harvested in the tri-trophic food chain.
a11 a22 < 0.
(9)
Likewise, tri-trophic food chain we can easily obtain a non-homogeneous system to compute the rate of change of species biomass with effort. We are not explicitly providing the matrix equation for the individual case. 3.2.1. Harvesting prey When prey is harvested, we can obtain the following inequalities for stock variation of individual species: 4
BioSystems 185 (2019) 104018
D. Pal, et al.
matrix evaluated at the equilibrium is negative. Hence, the steady-state of the unharvested system is asymptotically stable. We would like to determine the variation of stable steady stock of all the populations when prey is harvested. It is calculated that top-predator persists when effort ranges from e0 = 0 to e0 = 1.97 and the equilibrium within this effort range is hyperbolic. The variation of stock of all the species and maximum real part of all the eigenvalues are shown in Fig. 1. The corresponding coexisting stock within this effort is stable as the real part of all the eigenvalues is negative.1 It suggests that the prey harvesting does not destabilize the system which was asymptotically stable before harvesting. The change in the stock of prey species with respect to the effort is shown in the top left panel of Fig. 1. We can observe that prey biomass, at a stable state, increases with effort. Therefore, a hydra effect is exhibited when prey species is harvested. It is also examined that harvesting any other species does not have any positive impact towards the increase of stock. Moreover, there is no change in stability behavior when the harvesting is considered on the other trophic level. Thus, harvesting has a more stabilizing effect under such a parameter configuration.
a22 q0 x * dx* = > 0, de0 a11 a22 a12 a21 dy1* a21 q0 x * = < 0, de 0 a11 a22 a12 a21 dy *2 = 0, de 0 dy3* a21 a32 q0 x * = < 0. de 0 (a11 a22 a12 a21) a34 It claims that the hydra effect always occurs at stable state when prey species is harvested. 3.2.2. Harvesting first predator In this case we can obtain the following expressions:
a12 q1 y1* dx* = > 0, de1 a12 a21 a11 a22 dy1* a11 q1 y1* = < 0, de1 a12 a21 + a11 a22 dy *2 = 0, de1 dy3* a11 a32 q1 y1* = < 0. de1 (a12 a21 a11 a22) a34
3.3. Food chain with five trophic levels In the four-trophic system, we have detected the hydra effect only when prey species was harvested. In this subsection, we examine whether the hydra effect can appear only on prey species or it may appear on predator trophic levels in a community with five species. The unique interior equilibrium E (x *, y1* , y2* , y3* , y4* ) satisfies the characteristic equation
Since the biomass of the first predator decreases with effort, no hydra effect appears when the first predator is harvested. 3.2.3. Harvesting second predator Biomass of the target species is independent of effort, no hydra effect is observed.
5
de3
=
(
(a12 a21
+ Q3
2
+ Q4 + Q5 = 0,
(a11 + a22 + a33 + a44),
Q2 =
a12 a21
a23 a32 + a22 a33
a34 a43 + (a22 + a33) a44 a45 a54 ,
( a11 a23 a32 + a11 a22 a33
a11 a34 a43
a22 a34 a43 + a11 a22 a44
a23 a32 a44 + a11 a33 a44 + a22 a33 a44 Q4 Q5
)
+ a33 q3 y3*
a34 a43 (a11 a23 a32 + a33 (a12 a21
3
Q1 =
Q3 =
Finally,
=
+ Q2
+ a11 (a22 + a33 + a44 )
a12 a23 q3 y3* dx* = > 0, de3 (a12 a21 a11 a22) a43 dy1* a11 a23 q3 y3* = < 0, de3 (a12 a21 a11 a22) a43 dy *2 q y* = 3 3 > 0. de3 a43
dy3*
4
where
3.2.4. Harvesting top predator For harvesting the top predator, we can obtain
a11 a23 a32 a12 a21 a11 a22
+ Q1
a11 a22)) q3 y3*
a11 a22) a34 a43
a12 a21 (a33 + a44 )
(a11 + a22 + a33) a45 a54), = (a23 a32 a22 a33) a45 a54 + a12 a21 (a34 a43 a33 a44 + a45 a54 ) a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 + (a22 + a33) a45 a54), = (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 a54 .
The Routh-Hurwitz criteria for stability are: (i) Q1 > 0, (ii) Q5 > 0, (iii) Q1Q2 − Q3 > 0, (iv) Q1 Q2 Q3 Q32 Q12 Q4 + Q1 Q5 > 0 and (v) Q1 Q2 Q3 Q4 Q32 Q4 Q12 Q42 Q1 Q22 Q5 + Q2 Q3 Q5 + 2Q1 Q4 Q5 Now Q5 > 0 reduces to
.
In this case, we are not able to prove analytically the sign of the numerator for the above expression. However, many simulations reveal that the biomass of the top species decreases with effort. In conclusion, we observe that the hydra effect certainly occurs at stable state when prey is harvested. Harvesting any other trophic level decreases its density. One of the assumptions in our analysis was that harvesting can be continued until the equilibrium is stable. We do not have much information when the system loses its stability or it remains stable under harvesting. We provide some additional information for the above queries numerically. Moreover, we wish to visualize the change in the stock of all the trophic levels due to harvesting by a simulation technique. Considering parameter set r = 0.4, K = 20, α1 = 0.65, α2 = 0.8, α3 = 0.85, β1 = 0.8, β2 = 0.25, β3 = 0.11, m1 = 0.16, m2 = 0.1, m3 = 0.04, h1 = 8, h2 = 10 and h3 = 18 leads unique positive equilibrium P(8.5712, 5.8274, 2.1557, 0.7976) in the unharvested system. It is observed that the real part of all the eigenvalues of the Jacobian
Q52 > 0 .
a11 a23 a32 + a12 a21 a33 a11 a22 a33 > 0, i. e., a11 a23 a32 a11 a22 a33 > 0. Since, a11 < 0, we can write
a23 a32
a22 a33 < 0.
(10)
The analytical results for harvesting prey species are provided in Appendix A.1. It shows that increasing effort on prey harvesting reduces its equilibrium stock. Therefore, the hydra effect does not appear. We now employ effort on the first predator to examine the existence of the hydra effect. 1 Unlike the case of a smooth curve representing the prey stock, Max ( Re( i )) is a non-smooth curve. It is not surprising as Max ( Re( i )) collects the real part from eigenvalues with varying effort. More clearly, Max ( Re( i )) is not computed from only one eigenvalue.
5
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 1. The change in biomass of all the trophic levels is explored when the prey species is exploited in the food chain with three predator trophic levels. Biomass of the prey, first predator, second predator and top-predator increases, decreases, remains constant and decreases, respectively with effort. It is clear that the equilibrium value for the second predator does not depend on effort as its pre capita growth does not depends on its density. Hydra effect occurs on prey species. The change in the maximum value of the real part of all four eigenvalues shows that the system is stable under harvesting.
Harvesting first predator The rate of change in stock for different trophic levels under the harvesting of first predator is calculated as follows:
could not detect any hydra effect under the harvesting of the top-predator. For explanatory, we select the parameter set as r = 0.2, K = 20, α1 = 1.2, α2 = 1.5, α3 = 1.7, α4 = 1.9, β1 = 0.92, β2 = 0.8, β3 = 0.7, β4 = 0.4, m1 = 0.1, m2 = 0.07, m3 = 0.04, m4 = 0.02, h1 = 50, h2 = 55, h3 = 60 and h4 = 65 to show the existence of hydra effect when first predator is harvested. A unique coexisting equilibrium (13.3180, 5.8682, 5.3649, 1.7568, 2.0263), which is asymptotically stable, is obtained in the unharvested food chain with five trophic levels. All the populations persist when the first predator is harvested with fishing effort lesser than 1.47. Fig. 2 shows the equilibrium biomass y1* of the first predator increases with effort on first predator. We have computed the real part of the eigenvalues and observe that maximum value of the real parts is negative for any effort within e1 = 0 and e1 = 1.47. Thus, the equilibrium under harvesting is asymptotically stable and results in a hydra effect at a stable steady state. As before we have verified that harvesting either trophic level does not destabilize the coexisting steady state.
a12 a33 q1 y1* dx* = < 0, de1 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy1* a11 a33 q1 y1* = > 0, de1 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy *2 a11 a32 q1 y1* = < 0, de1 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy3* = 0, de1 dy *4 a11 a32 a43 q1 y1* = < 0. de1 a11 a23 a32 a45 + a12 a21 a33 a45 a11 a22 a33 a45 The above expression shows that the equilibrium stock of the first predator increases when it is harvested with increasing effort. Thus, the hydra effect, at stable equilibrium, appears within the whole range of effort where all populations coexist. The mathematical calculations in estimating the stock for harvesting the second predator are postponed to Appendix A.2. The main key to understand the ecological result depends upon sign of (a12a21 − a11a22). Several numerical simulations could not detect any hydra effect. On the other hand, one can easily prove that the stable stock of the third predator is independent of effort when it is harvested. We have further represented the analytical expression in Appendix A.3 when the top-predator is harvested. However, the complex form of the expressions prevents to establish any concrete conclusion on the variation of stock with effort. All parameter sets, we have selected so far,
3.4. Food chain with six trophic levels So far we have been able to detect the hydra effect only on a single trophic level in the food chain with four or five species. Now we analyze a food chain with six trophic levels. The eigenvalues of the Jacobian matrix corresponding to the food chain satisfy the equation 6
+ Q1
where
6
5
+ Q2
4
+ Q3
3
+ Q4
2
+ Q5 + Q6 = 0,
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 2. Impacts on biomass at stable steady state for all the trophic levels are shown when the first predator harvested in a food chain with five trophic levels. Hydra effect is exhibited on the first predator. The bottom panel claims that all the eigenvalues has a negative real part. It also shows that harvesting does not destabilize the dynamics of the equilibrium.
Q1 = Q2 = Q3 =
(a11 + a22 + a33 + a44 + a55), a12 a21 a23 a32 a34 a43 + a33 a44 a45 a54 + (a33 + a44) a55 + a22 (a33 + a44 + a55) + a11 (a22 + a33 + a44 + a55) a56 a65 , ( a22 a34 a43 + a22 a33 a44
a22 a45 a54
Q6 = (a12 a21 ( a34 a43 + a33 a44 ) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a56 a65. According to Routh-Hurwitz criteria, the following conditions are to be satisfied for obtaining negative real part of the eigenvalues: (i) Q1 > 0, (ii) Q6 > 0, (iii) Q1Q2 − Q3 > 0, (iv) Q1 Q2 Q3 Q32 Q12 Q4 + Q1 Q5 > 0 , Q1 Q2 Q3 Q4 Q32 Q4 Q12 Q42 Q1 Q22 Q5 + Q2 Q3 Q5 + 2Q1 Q4 Q5 (v) Q52 + Q12 Q2 Q6 Q1 Q3 Q6 > 0 and (vi) Q1 Q2 Q3 Q4 Q5 Q32 Q4 Q5 Q12 Q42 Q5 Q1 Q22 Q52 + Q2 Q3 Q52 + 2Q1 Q4 Q52 Q53 Q1 Q2 Q32 Q6 + Q33 Q6 + Q12 Q3 Q4 Q6 + 2Q12 Q2 Q5 Q6
a33 a45 a54 + a22 a33 a55
a34 a43 a55 + a22 a44 a55 + a33 a44 a55
a23 a32 (a44 + a55)
a12 a21 (a33 + a44 + a55)
(a22 + a33 + a44) a56 a65 + a11 ( a23 a32
a34 a43 + a33 a44
+ a22 (a33 + a44 + a55) Q4 = a23 a32 a45 a54 + a22 a33 a44 a55
a45 a54 + (a33 + a44 ) a55
a56 a65)),
a22 a33 a45 a54
a22 a34 a43 a55
a11 (a33 a45 a54 + a34 a43 a55
3Q1 Q3 Q5 Q6 Q13 Q62 > 0. Now Q6 > 0 implies
a23 a32 a44 a55
a12 a21 (a33 a44
a33 a44 a55
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44 ) < 0. (11)
+ a23 a32 (a44 + a55) + a22 (a34 a43
a33 a44 + a45 a54
+ a44 ) a56 a65 + (a23 a32 + a34 a43
(a33 + a44) a55))
The left hand side quantity will appear in the denominator for all cases in the expression of rate of change of the stock due to harvesting effort. From Q6 > 0 we can also obtain
a11 (a22 + a33
a22 a33
(a22 + a33) a44) a56 a65 + a12 a21 (a34 a43
a33 a44 + a45 a54
(a33 + a44) a55 + a56 a65),
Q5 =
((a22 a34 a43 + a23 a32 a44
a12 a21 (a33 a44 i. e., (a12 a21
a22 a33 a44 ) a56 a65
Thus,
+ a12 a21 (a33 a45 a54 + a34 a43 a55 a33 a44 a55 + (a33 + a44 ) a56 a65) + a11 ((a34 a43
a34 a43) + a11 (a22 a34 a43 a22 a33 a44 ) < 0, a11 a22 )(a33 a44 a34 a43 ) < 0.
a12 a21
a33 a44) a56 a65
a11 a22 < 0.
(12)
Now we examine the impacts on stock levels when individual trophic level.
+ a23 a32 (a45 a54 a44 a55 + a56 a65) a22 (a33 a45 a54 + a34 a43 a55 a33 a44 a55 + (a33 + a44 ) a56 a65)))
3.4.1. Harvesting prey When prey is harvested, the rate of change in the stock of individual
and 7
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 3. The variation of the biomass of all the trophic levels is shown when prey species is harvested in the six-species model. Clearly, the hydra effect occurs on the prey species. The bottom panel represents the largest of real part among all eigenvalues.
species with effort e0 are presented as:
a12 a23 a44 q2 y2* dx* = > 0, de2 a12 a21 (a33 a 44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44) dy1* a11 a23 a44 q2 y2* = < 0, de2 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy 2* (a12 a21 a11 a22) a44 q2 y2* = de2 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )
(a22 a34 a 43 + a23 a32 a 44 a22 a33 a44) q0 x * dx* = > 0, de0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a 43 + a23 a32 a 44 a22 a33 a44) dy1* a21 (a34 a43 a33 a 44 ) q0 x * = < 0, de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44 ) * dy 2 a21 a32 a44 q0 x * = > 0, de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44 ) dy*3 a21 a32 a43 q0 x * = < 0, de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a 43 + a23 a32 a 44 a22 a33 a44) dy*4 = 0, de 0 dy 5* a21 a32 a43 a54 q0 x * = < 0. de 0 a12 a21 (a33 a44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44 )) a56
> 0, (using Eq. (12)) dy*3 de2 dy*4 de2 dy*5 de2
=
(a12 a21 a12 a21 (a33 a 44
a11 a22) a 43 q2 y2*
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)
< 0,
= 0, =
(a12 a21 (a12 a21 (a33 a 44
a11 a22) a 43 a54 q2 y2*
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)) a56
< 0.
Our above analysis reveals that the hydra effect is ensured when the second predator is harvested. The mathematical calculations are provided in Appendix B.2 when third predator is subjected to harvested. In this situation, we face difficulty to understand the variation of the biomass of the third predator. Several numerical simulations do not show any increase in biomass on the third predator. Thus, it is unlike to experience the hydra effect on the third predator. However, we must admit that an affirmative answer is yet to investigate. It is again cleared from Appendix B.3 that stock of the forth predator is independent of e4. Hence, we cannot expect any hydra effect. When the top predator is removed, it is found that variation of biomass for both the third and top predator will be the same (See Appendix B.4). However, analytically we could not establish any affirmative statement. Simulations predict that no hydra effect is observed when the top species is harvested. For further explanation, in our six species food chain, parameters
The sign of the numerator in the above quantities can easily be determined by identifying the sign of the individual entries of the Jacobian matrix. It is cleared that prey biomass increases with effort and results in a hydra effect. In the case of harvesting the first predator, the hydra effect, at stable state, is not appeared as the biomass of the first predator follows a decreasing trend for all effort range. The mathematical derivations are detailed in Appendix B.1. 3.4.2. Harvesting second predator Here the variation of stock with effort e2 can be seen as:
8
BioSystems 185 (2019) 104018
D. Pal, et al.
Fig. 4. In the six-species model, the biomass of all the species is presented when the second predator is harvested. The bottom panel claims that all the eigenvalues have a negative real part.
are chosen as r = 0.2, K = 30, α1 = 1.2, α2 = 1.5, α3 = 1.7, α4 = 1.9, α5 = 10, β1 = 0.9, β2 = 0.8, β3 = 0.7, β4 = 0.4, β5 = 0.1, m1 = 0.1, m2 = 0.07, m3 = 0.04, m4 = 0.02, m5 = 0.01, h1 = 40, h2 = 45, h3 = 50, h4 = 55 and h5 = 60. This set produces a stable steady state (7.9195, 5.8783, 2.6622, 2.1264, 0.6061, 0.05024). For increasing effort for prey species within the range (0, 0.475) increases the prey density (Fig. 3). The top-predator population decreases with effort at equilibrium (Fig. 3). As the top-predator biomass is positive within the said effort range, biomass of all the species are positive at equilibrium. The steady state is asymptotically stable as the maximum real part of all the eigenvalues are negative (bottom panel). Hence, the hydra effect appears at stable state when prey species is harvested. Computation reveals that the effort range for harvesting the second predator should lie within (0, 0.41) in order to maintain the coexistence of all the species. Fig. 4 shows that the biomass of the second predator is increased with effort. The biomass of the top-predator tends to zero as e2 → 0.41. Then all the eigenvalues have a negative real part in the aforesaid effort interval as can be seen from the figure. Hence, the hydra effect is observed at a stable state when the second predator is harvested. Unlike the case of the previous two models, we here able to identify the existence of the hydra effect in multiple trophic levels. We have also verified that harvesting either trophic level under this parameter restriction preserves the asymptotic local stability behavior of the coexisting equilibrium.
by Sieber and Hilker (2012) to prove the existence of a hydra effect. However, that effect at stable state was not detected either in a RM models (Sieber and Hilker, 2012) or in tri-trophic food chains (Ghosh et al., 2018). The purpose of this article was to verify the existence of the hydra effect in the food chains. We have proposed a linear food chain model with logistic prey growth and Holling type II functional response. We posed the question of whether the hydra effect, at a stable state, can be a special feature in “pure predator systems” with higher trophic levels. An affirmative answer is found through the study of food chains up to six trophic levels successively. Although prey-predator models with two species are investigated in establishing many ecological phenomenons, multiple trophic levels are common in a real ecosystem. For instance, four- (Pauly et al., 2002) and five-trophic (Pauly and Watson, 2003) system for fish dynamics are available in literature. Therefore, proposing and analyzing ecological systems with higher trophic levels system are justified. We observed that hydra effect, at stable state, occurs on (a) prey in a four-trophic system (Fig. 1), (b) first predator in a five-trophic system (Fig. 2) and (c) prey (Fig. 3) and second predator (Fig. 4) in a six-trophic system. To the best of authors’ knowledge, we are the first to establish the hydra effect in relatively simple food chain models involving unique stable fixed point. Therefore, hydra effect may be observed even in the absence of nonlinear mortality term among different predator trophic levels (Cortez and Abrams, 2016), stage structure of community (Abrams and Quince, 2005; Costa et al., 2017) or Allee effect among species (Costa and dos Anjos, 2018). The existence of the hydra effect as well as the variation of stable stock under harvesting effort is shown in Table 1, when an individual
4. Results and conclusion Gause type prey-predator systems under cyclic dynamics are studied 9
BioSystems 185 (2019) 104018
D. Pal, et al.
number of hydra effect species is ([n/2] − 1). Also, two consecutive trophic levels in our models cannot experience the hydra effect. These statements can be considered as conjectures. This study is important in fishery management including biological conservation for estimating stock and yield correctly. Most importantly, knowledge of the hydra effect might help to implement proper action in the pest control program. Because, an attempt for reducing pest or removing/culling arbitrary species, may lead to a higher density of pest populations (occurrence of hydra effect) in complex pest-predator (natural enemy) systems. In this case, farmers and gardeners should look for alternative action plans beforehand. In some situations, we could not derive results analytically about the hydra effect. However, our predictions are reported above. All the stability conditions are not utilized for estimating the variation of stock. The complete scenario will be cleared when the remaining RouthHurwitz criteria can be applied effectively. One can take the work forward to prove these issues. Thus, a general framework could be established analytically to identify the hydra affect species in a food chain with an arbitrary number of trophic levels. This is an interesting and challenging problem as a future work from mathematical and ecological viewpoints. Studying the unstructured food chain models, in the present contribution, is the first step and hence we intend to work on the same theme for stage-structured prey-predator models in a continuous and discrete-time framework. Further, the existence of the hydra effect in a spatially structured community could equally be interesting as a future perspective.
Table 1 The change of biomass of all the trophic levels is shown when a single trophic is targeted for harvesting. We also reported if the targeted species is a hydra affect species. Here the symbol ↑ (resp. ↓) indicates that the biomass of a particular trophic level is increased (resp. decreased) with effort, whereas no variation of biomass is represented by the symbol −. The observation based on numerical simulation is indicated by *. Model
Target species
x*
y* 1
y* 2
Tri-trophic
x(t) y1(t) y2(t)
↓ – ↓
– – ↑
↓ ↓ ↓
Four-trophic
x(t) y1(t) y2(t) y3(t)
↑ ↑ – ↑
↓ ↓ – ↓
– – – ↑
↓ ↓ ↓ *↓
Five-trophic
x(t) y1(t) y2(t) y3(t) y4(t)
↓ ↓ ↓ – ↓
↑ ↑ ↑ – ↑
↓ ↓ *↓ – *↓
– – – – ↑
↓ ↓ *↓ ↓ *↓
Six-trophic
x(t) y1(t) y2(t) y3(t) y4(t) y5(t)
↑ ↑ ↑ ↑ – ↑
↓ ↓ ↓ ↓ – ↓
↑ ↑ ↑ ↑ – ↑
↓ ↓ ↓ *↓ – *↓
– – – – ↑
y* 3
y* 4
y* 5
Hydra effect NO NO NO YES NO NO *NO NO YES *NO NO *NO
↓ ↓ ↓ *↓ ↓ *↓
YES NO YES *NO NO *NO
trophic level is harvested in different food chains. Interestingly, one can observe that the rule of biomass variation maintains a reverse order from one trophic level to another one when the hydra effect appears. This order does not match the conventional rule. Legović et al. (2010) generalized a rule for the change of biomass due to harvesting in the ecosystem by investigating a simple prey-predator model. Our study suggests that Legović et al. (2010)’ results may not hold for some specific systems in the presence of a hydra effect. From our experience, we predict that a food chain with even (resp. odd) number (n) of trophic levels, hydra effect occurs at odd (resp. even) trophic levels except possibly the top three species. Moreover, the
Conflicts of interest The authors declare no conflicts of interest. Acknowledgment We appreciate the substantial comments and suggestions received from both the eminent reviewers and the Editor, which certainly helped to improve the content of the manuscript.
Appendix A. Analysis of food chain with five trophic levels Here we present some of the mathematical results and their ecological relevances when an individual trophic level is harvested. A.1 Harvesting prey Prey harvesting generates
(a23 a32 a22 a33) q0 x * dx* = < 0, de0 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy1* a21 a33 q0 x * = > 0, de 0 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy *2 a21 a32 q0 x * < 0, = de 0 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy3* = 0, de 0 dy *4 a21 a32 a43 q0 x * = < 0. de 0 a11 a23 a32 a45 + a12 a21 a33 a45 a11 a22 a33 a45 Here prey harvesting reduces its own biomass, and hence no hydra effect is possible. A.2 Harvesting second predator The variation of biomass at the equilibrium of the species are determined from the followings:
10
BioSystems 185 (2019) 104018
D. Pal, et al.
a12 a23 q2 y2* dx* = < 0, de2 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy1* a11 a23 q2 y2* = > 0, de2 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy *2 (a12 a21 a11 a22) q2 y2* = , de2 a11 a23 a32 + a12 a21 a33 a11 a22 a33 dy3* = 0, de2 dy *4 ( a12 a21 + a11 a22) a43 q2 y2* = . de2 (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 We are not able to determine the sign of (a12a21 − a11a22). But numerically, we have observed that the sign is negative. Hence it is more likely that no hydra effect appears when the second predator is harvested. A.3 Harvesting top predator When top predator is harvested, we observe that
a12 a23 a34 q4 y4* dx* = < 0, de4 a11 a23 a32 a54 + a12 a21 a33 a54 a11 a22 a33 a54 dy1* a11 a23 a34 q4 y4* = > 0, de4 a11 a23 a32 a54 + a12 a21 a33 a54 a11 a22 a33 a54 dy *2 ( a12 a21 + a11 a22) a34 q4 y4* = , de4 (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a54 dy3* q y* = 4 4 > 0. de4 a54 Now
dy*4 de4
=
a22 a33 a44 )) q4 y4*
(a12 a21 ( a34 a43 + a33 a44 ) + a11 (a22 a34 a43 + a23 a32 a44 (a11 a23 a32 + a12 a21 a33
a11 a22 a33) a45 a54
.
Again the analytical result is not available to know if the biomass of the top predator increases for some effort range when it is harvested. Thus, we cannot arrive at any specific conclusion. However, many simulations show a decreasing trend in biomass of the target species. Appendix B. Analysis of food chain with six trophic levels We derive mathematical results in support of the statements provided in Section 3.4. B.1 Harvesting first predator In this case, we have
dx* = de1 dy1* = de1 dy *2 = de1 dy3* = de1 dy *4 = de1 dy5* = de1
a12 (a34 a43 a12 a21 (a33 a44
a33 a44 ) q1 y1*
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a11 (a34 a43 a33 a44 ) q1 y1*
a22 a33 a44 )
a12 a21 ( a34 a43 + a33 a44) + a11 (a22 a34 a43 + a23 a32 a44 a11 a32 a44 q1 y1*
a22 a33 a44)
a12 a21 (a33 a44 a12 a21 (a33 a44
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a11 a32 a43 q1 y1* a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)
> 0, < 0,
> 0,
a22 a33 a44 )
< 0,
0, a11 a32 a43 a54 q1 y1* (a12 a21 (a33 a44
a34 a43) + a11 (a22 a34 a43 + a23 a32 a44
a22 a33 a44)) a56
< 0.
The above inequalities hold true as (a34a43 − a33a44) < 0 along with a11 < 0 and a12 < 0. Hence, the hydra effect is not appeared. B.2 Harvesting third predator When e3 is the control parameter, we can obtain
11
BioSystems 185 (2019) 104018
D. Pal, et al.
a12 a23 a34 q3 y3* dx* = > 0, de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy1* a11 a23 a34 q3 y3* = < 0, de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy *2 (a12 a21 a11 a22 ) a34 q3 y3* = > 0, (using Eq.(12)) de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 ) dy3* (a11 a23 a32 + a12 a21 a33 a11 a22 a33) q3 y3* , = de3 a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44) dy *4 = 0, de3 dy5* (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a54 q3 y3* = . de3 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a56 We find difficult to understand the variation of biomass for the third predator. Numerical results do not show any increase in biomass for third predator and top predator biomass decreases with effort. It is unlike to happen hydra effect on the third predator. B.3 Harvesting fourth predator Since the biomass of this trophic level is independent of e4, we cannot expect any hydra effect. Besides, we would like to know the variation of the biomass of other trophic levels under harvesting. To estimate the variation we represent the following system of equations:
dy * dx* + a12 1 de4 de4 dy1* dx* a21 + a22 de4 de4 dy1* dy2* a32 + a33 de4 de4 dy* dy*2 a43 + a44 3 de4 de4 dy*3 dy *4 a54 + a55 de4 de4 dy *4 a65 de4 a11
= 0, + a23 + a34 + a45 + a56
dy *2
= 0,
de4 dy*3
= 0,
de4 dy *4 de4 dy5* de4
= 0, q4 y4* = 0, = 0.
dy*
As de 4 = 0 , the first four equations form a homogeneous system with coefficient matrix non-singular. Thus the sub-system with the first four 4 equations has a unique solution. As a result, there is no change in biomass for all the trophic levels below y4(t). Finally, the fifth equation reveals that top predator's biomass decreases with effort. Although this analysis is not provided explicitly for food chain with smaller trophic levels. However, the same analysis can be applied easily. B.4 Harvesting top predator When the top predator is removed, we have
a12 a23 a34 a45 q5 y5* dx* = > 0, de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a65 dy1* a11 a23 a34 a45 q5 y5* = < 0, de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44)) a65 dy *2 ( a12 a21 + a11 a22 ) a34 a45 q5 y5* = > 0, de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a65 dy3* (a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 q5 y5* = , de5 (a12 a21 (a33 a44 a34 a43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a44 )) a65 dy *4 q5 y5* = > 0, de5 a65 dy5* de5
=
(
(a11 a23 a32 + a12 a21 a33 a11 a22 a33) a45 a54 a12 a21 (a33 a 44 a34 a 43) + a11 (a22 a34 a43 + a23 a32 a44 a22 a33 a 44)
a56 a65
+ a55 ) q5 y5*
.
It is cleared that the variation of biomass for both the third and top predator will be the same. However, analytically we could not establish any affirmative statement. Simulations confirm that no hydra effect is observed when top species is harvested. Appendix C. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.biosystems.2019.104018.
12
BioSystems 185 (2019) 104018
D. Pal, et al.
References
Matsuda, H., Abrams, P.A., 2004. Effects of predator prey interactions and adaptive change on sustainable yield. Can. J. Fish. Aquat. Sci. 61 (2), 175–184. Neverova, G., Abakumov, A., Yarovenko, I., Frisman, E.Y., 2018. Mode change in the dynamics of exploited limited population with age structure. Nonlinear Dyn. 94 (2), 827–844. Pauly, D., Christensen, V., Guénette, S., Pitcher, T.J., Sumaila, U.R., Walters, C.J., Watson, R., Zeller, D., 2002. Towards sustainability in world fisheries. Nature 418 (6898), 689. Pauly, D., Watson, R., 2003. Counting the last fish. Sci. Am. 289 (1), 42–47. Ricker, W.E., 1954. Stock and recruitment. J. Fish. Board Can. 11 (5), 559–623. Schröder, A., van Leeuwen, A., Cameron, T.C., 2014. When less is more: positive population-level effects of mortality. Trends Ecol. Evol. 29 (11), 614–624. Sieber, M., Hilker, F.M., 2012. The hydra effect in predator-prey models. J. Math. Biol. 64 (1-2), 341–360. Takashina, N., Mougi, A., Iwasa, Y., 2012. Paradox of marine protected areas: suppression of fishing may cause species loss. Popul. Ecol. 54 (3), 475–485. Weide, V., Varriale, M.C., Hilker, F.M., 2019. Hydra effect and paradox of enrichment in discrete-time predator-prey models. Math. Biosci. 310, 120–127. Wilen, C.D., Wilen, J.E., 2012. Fishing down the food chain revisited: Modeling exploited trophic systems. Ecol. Econ. 79, 80–88. Wollrab, S., Diehl, S., De Roos, A.M., 2012. Simple rules describe bottom-up and topdown control in food webs with alternative energy pathways. Ecol. Lett. 15 (9), 935–946. Xiao, M., Cao, J., 2009. Hopf bifurcation and non-hyperbolic equilibrium in a ratio-dependent predator-prey model with linear harvesting rate: analysis and computation. Math. Comput. Model. 50 (3-4), 360–379.
Abrams, P., 2015. The hydra effect is no myth. New Sci. 226 (3023), 28–29. Abrams, P.A., Matsuda, H., 2005. The effect of adaptive change in the prey on the dynamics of an exploited predator population. Can. J. Fish. Aquat. Sci. 62 (4), 758–766. Abrams, P.A., Quince, C., 2005. The impact of mortality on predator population size and stability in systems with stage-structured prey. Theor. Popul. Biol. 68 (4), 253–266. Abrams, P.A., Roth, J.D., 1994. The effects of enrichment of three-species food chains with nonlinear functional responses. Ecology 75 (4), 1118–1130. Cortez, M.H., Abrams, P.A., 2016. Hydra effects in stable communities and their implications for system dynamics. Ecology 97 (5), 1135–1145. Costa, M.I.d.S., dos Anjos, L., 2018. Multiple hydra effect in a predator–prey model with allee effect and mutual interference in the predator. Ecol. Model. 373, 22–24. Costa, M.I.D.S., Esteves, P.V., Faria, L.D.B., dos Anjos, L., 2017. Prey dynamics under generalist predator culling in stage structured models. Math. Biosci. 285, 68–74. Gantmacher, F.R., 2000. The Theory of Matrices, vol. 131 American Mathematical Soc. Ghosh, B., Kar, T., 2013. Possible ecosystem impacts of applying maximum sustainable yield policy in food chain models. J. Theor. Biol. 329, 6–14. Ghosh, B., Kar, T., Legović, T., 2014. Relationship between exploitation, oscillation, msy and extinction. Math. Biosci. 256, 1–9. Ghosh, B., Pal, D., Legović, T., Kar, T., 2018. Harvesting induced stability and instability in a tri-trophic food chain. Math. Biosci. 304, 89–99. Gross, T., Ebenhöh, W., Feudel, U., 2005. Long food chains are in general chaotic. Oikos 109 (1), 135–144. Legović, T., Klanjšček, J., Geček, S., 2010. Maximum sustainable yield and species extinction in ecosystems. Ecol. Model. 221 (12), 1569–1574. Liz, E., Ruiz-Herrera, A., 2012. The hydra effect, bubbles, and chaos in a simple discrete population model with constant effort harvesting. J. Math. Biol. 65 (5), 997–1016.
13