Outline
Comptes Rendus

Biological modelling/Biomodélisation
Chaotic dynamics of a three species prey–predator competition model with bionomic harvesting due to delayed environmental noise as external driving force
Comptes Rendus. Biologies, Volume 335 (2012) no. 8, pp. 503-513.

Abstract

We consider a biological economic model based on prey–predator interactions to study the dynamical behaviour of a fishery resource system consisting of one prey and two predators surviving on the same prey. The mathematical model is a set of first order non-linear differential equations in three variables with the population densities of one prey and the two predators. All the possible equilibrium points of the model are identified, where the local and global stabilities are investigated. Biological and bionomical equilibriums of the system are also derived. We have analysed the population intensities of fluctuations i.e., variances around the positive equilibrium due to noise with incorporation of a constant delay leading to chaos, and lastly have investigated the stability and chaotic phenomena with a computer simulation.

Metadata
Received:
Accepted:
Published online:
DOI: 10.1016/j.crvi.2012.06.001
Keywords: Prey–predator, Local and global stability, Bionomic equilibrium, Stochastic delayed perturbation, Fourier transform methods, Chaos

Kalyan Das 1; M.N. Srinivas 2; M.A.S. Srinivas 3; N.H. Gazi 4

1 National Institute of Food Technology Entrepreneurship and Management, Department of Mathematics, Plot No. 97, Sector 56, HSIIDC Industrial Estate, Kundli 131 028, Haryana, India
2 School of Advanced Sciences, Department of Mathematics, VIT University, Vellore 632 014, Tamil Nadu, India
3 Department of Mathematics, Jawaharlal Nehru Technological University, Hyderabad, Andhra Pradesh, India
4 Department of Mathematics, Aliah University, DN-41, Salt Lake, Sector V, Kolkata 700 091, India
@article{CRBIOL_2012__335_8_503_0,
     author = {Kalyan Das and M.N. Srinivas and M.A.S. Srinivas and N.H. Gazi},
     title = {Chaotic dynamics of a three species prey{\textendash}predator competition model with bionomic harvesting due to delayed environmental noise as external driving force},
     journal = {Comptes Rendus. Biologies},
     pages = {503--513},
     publisher = {Elsevier},
     volume = {335},
     number = {8},
     year = {2012},
     doi = {10.1016/j.crvi.2012.06.001},
     language = {en},
}
TY  - JOUR
AU  - Kalyan Das
AU  - M.N. Srinivas
AU  - M.A.S. Srinivas
AU  - N.H. Gazi
TI  - Chaotic dynamics of a three species prey–predator competition model with bionomic harvesting due to delayed environmental noise as external driving force
JO  - Comptes Rendus. Biologies
PY  - 2012
SP  - 503
EP  - 513
VL  - 335
IS  - 8
PB  - Elsevier
DO  - 10.1016/j.crvi.2012.06.001
LA  - en
ID  - CRBIOL_2012__335_8_503_0
ER  - 
%0 Journal Article
%A Kalyan Das
%A M.N. Srinivas
%A M.A.S. Srinivas
%A N.H. Gazi
%T Chaotic dynamics of a three species prey–predator competition model with bionomic harvesting due to delayed environmental noise as external driving force
%J Comptes Rendus. Biologies
%D 2012
%P 503-513
%V 335
%N 8
%I Elsevier
%R 10.1016/j.crvi.2012.06.001
%G en
%F CRBIOL_2012__335_8_503_0
Kalyan Das; M.N. Srinivas; M.A.S. Srinivas; N.H. Gazi. Chaotic dynamics of a three species prey–predator competition model with bionomic harvesting due to delayed environmental noise as external driving force. Comptes Rendus. Biologies, Volume 335 (2012) no. 8, pp. 503-513. doi : 10.1016/j.crvi.2012.06.001. https://comptes-rendus.academie-sciences.fr/biologies/articles/10.1016/j.crvi.2012.06.001/

Version originale du texte intégral

1 Introduction

The prey–predator system is an important feature in population dynamics and has been studied by many authors [1–4]. It is already known in the classical predator–prey system where each individual predator possesses the ability to attack prey. Even since research in the discipline of theoretical ecology was initiated by Lotka [5] and Volterra [6], several mathematicians and ecologists contributed to the growth of this area of knowledge and this has been extensively reported in the treatises of Meyer [7], Cushing [8], Conlinvaux [9], Freedman [10], Kapur [11,12]. Recently, the optimal management of renewable resources which has a direct impact on sustainable development, has been studied extensively by Chaudhuri [13], Kar and Swarnakamal [14], and Clark [15]. Now, people are facing the problems due to shortage of resource management. Extensive and unregulated harvesting of marine fish has lead to the extinction of several fish species. These problems are seen in marine reserved zones, where fishing and other activities are prohibited. This marine reserve not only protects species inside the reserve area, but also increases fish abundance in adjacent areas. The model of ecological system reflecting this has been given by Kar and Swarnakamal [14] and Rui Zhang et al. [16]. The ecological interactions are broadly classified as prey-predation, competitions, neutralism, mutualism and so on.

Harvesting has a strong impact on the dynamic evolution of a population. We know, depending on the nature of the applied harvesting strategy, the long-run stationary density of the population may be significantly smaller than the long-run stationary density of a population in the absence of harvesting [15].

Generally, a bionomic model consists of a biological or biophysical type which describes the behaviour of a living system, and an economic model which relates the biological system to market prices and resources with institutional constraints. Bio-economic models contain a single mathematical equation to represent a biological process. The logistic equation is the most commonly used function to capture the essential features of population densities in fishery and forestry management. However, there exists an increasing trend towards simulation models which are developed by biologists and agricultural scientists. These types of models also model the approximate dynamical behaviour in real situations and their complexity may preclude their use directly as part of optimal control models.

As harvesting has a strong impact on the dynamic evolution of a population, depending on the nature of applied harvesting strategy, the long-run stationary density of the population may be significantly smaller than the long-run stationary density of a population in the absence of harvesting. In the absence of harvesting, a population can be free of extinction risk; however, harvesting can lead to the incorporation of a positive extinction probability and to potential extinction in finite time. If a population is subject to a positive extinction rate then harvesting can drive the population density to a dangerous low level at which extinction becomes sure, no matter how the harvester affects the population afterwards. The exploitation of biological resources and harvest of population species are commonly practiced in fishery, forestry and wildlife management. Fishery, an ancient human tradition, has satisfied the food needs of mankind for thousands of years and has become economically, socially and culturally fundamental. Nowadays, fishes are in real trouble as their populations are being depleted to dangerous low levels, and it is necessary to discuss further in order to understand short- and long-term exploitation patterns in fishery management [17].

The problems of predator–prey systems in the presence of harvesting have been discussed by many authors; most of them have focused attention on optimal exploitation guided entirely by profits from harvesting. Brauer and Soudack [18,19] studied a class of predator–prey models under a constant rate of harvesting and under a constant quota of harvesting of both species simultaneously. The prey–predator model with harvesting was also studied by Dai and Tang [20], Myerscough et al. [21], Kar and Chaudhuri [22,23].

Recently, stability analysis of three species with environmental fluctuations was investigated by Gazi et al. [24]. Local stability analysis for a two-species ecological mutualism model has been presented by Reddy et al. [25] and stability analysis of neutral species was carried out by Reddy et al. [26]. Also, the stability analysis of prey, predator and super-predator was carried out by Reddy et al. [2,27]. In 2006, Carletti [28] considered a delay differential equations model with bacteriophage infection and discussed the robustness of the positive equilibrium with respect to stochastic perturbations of the environment using two different approaches. He investigated the analytical estimates of the population intensities with fluctuations by Fourier transform methods. Extensive numerical simulation suggested that a noisy environment for the bacteria population has much more destabilizing behaviour on the concentrations at the equilibrium point than a noisy environment for the phage.

The present investigation is an analytical study of three species system comprising a Prey (S1) common to two predators (S2) and (S3) which are in competition with each other. Some of the equilibrium points are identified based on the model equations and these are given in four distinct classes:

  • • in the absence of all species;
  • • in the absence of second predator;
  • • in the absence of first predator;
  • • in the presence of all the species (the co-existent state).

This article is organized in the following way. In the next section, we study the existence of local stability of the equilibria and their dependence on the harvesting efforts. We have concentrated much more on the interior equilibrium of the system as we are interested in the existence of the species. Next, we derived the criteria for the global stabilities of the system. Taking simple economic considerations into account, we have discussed the possibilities of the existence of a bionomic equilibrium when the system is exploited. We also computed the population intensities of fluctuations, i.e., variances around the positive equilibrium due to incorporation of noise which leads to chaos in reality [29]. Some numerical results are presented. The problem ends with a brief description of the principal results obtained here. In particular, our present model could be applied to fish species like goldband fish as a prey, sharks as a first predator and baleen whales or dolphins as second predator in the real ecological arena.

2 Mathematical model

We consider three species multisystem model by the following set of three non-linear ordinary differential equations:

dx1dt=a1x1α11x12+α12x1x2α13x1x3q1E1x1(1)
dx2dt=a2x2α22x22α21x1x2α23x2x3q2E2x2(2)
dx3dt=a3x3α33x32+α31x1x3α32x2x3q3E3x3(3)
with the following notation:
  • xi(t): population density of the species Si at time t; i = 1, 2, 3;
  • ai: the natural growth rates of Si, i = 1, 2, 3;
  • αii: the rate of decrease of Si due to its own insufficient resources, i = 1, 2, 3;
  • • α1n: the rate of decrease of the prey due to inhibition by the predators, n = 2, 3;
  • • αm1: the rate of increase of the predator due to its successful attacks on the prey, m = 2, 3;
  • Ki = aiii: carrying capacities of species, i = 1, 2, 3.

Further, the variables x1, x2 and x3 are non-negative and the model parameters ai, Ki, αij, i = 1, 2, 3, j = 1, 2, 3.

  • q1: catch ability coefficient of prey species;
  • q2: catch ability coefficient of predator species;
  • q3: catch ability coefficient of predator species;
  • E1: effort applied to harvest the prey species (S1);
  • E2: effort applied to harvest the predator species (S2);
  • E3: effort applied to harvest the predator species (S3).

Throughout our analysis, we take the initial conditions as

a2q2E2>0,a3q3E3>0,a1q1E1>0(4)

3 Existence of equilibrium points

The steady state equations of (1)–(3) are

a 1 x 1 α 11 x 1 2 α 12 x 1 x 2 α 13 x 1 x 3 q 1 E 1 x 1 = 0 a 2 x 2 α 22 x 2 2 + α 21 x 1 x 2 α 23 x 2 x 3 q 2 E 2 x 2 = 0 a 3 x 3 α 33 x 3 2 + α 31 x 1 x 3 α 32 x 2 x 3 q 3 E 3 x 3 = 0

System (1)–(3) possesses the following steady states:

  • • the trivial state G1 (0,0,0) (in the absence of all the species);
  • • the axial state G2(x¯1,x¯2,0) (in the absence of second predator);
  • • the boundary state G3x1ϕ,0,x3ϕ (in the absence of first predator);
  • • the steady state of coexistence G4 (x1*, x2*, x3*) (the interior equilibrium).

Case (i): G1 (0,0,0) i.e. the population is extinct and this state always exists.

Case (ii): if x¯1 and x¯2 are the positive solutions of

a1x1α11x12α12x1x2α13x1x3q1E1x1=0a2x2α22x22+α21x1x2α23x2x3q2E2x2=0

Then

x¯1=α22(a1q1E1)α12(a2q2E2)α21α12+α11α22x¯2=α21(a1q1E1)+α11(a2q2E2)α21α12+α11α22>0

x¯1 is positive provided the following inequality holds:

α22(a1q1E1)>α12(a2q2E2)

Case (iii): if x1ϕ and x3ϕ are the positive solutions of

a1x1α11x12α12x1x2α13x1x3q1E1x1=0a3x3α33x32+α31x1x3α32x2x3q3E3x3=0

Then

x1ϕ=α33a1q1E1+α13a3q3E3α13α31+α11α33
x3ϕ=α31a1q1E1+α11a3q3E3α13α31+α11α33

x1ϕ is positive provided the following inequality holds:

α33(a1q1E1)>α13(a3q3E3)

Case (iv): if x1* and x2* are the positive solutions of

a1x1α11x12α12x1x2α13x1x3q1E1x1=0a2x2α22x22+α21x1x2α23x2x3q2E2x2=0a3x3α33x32+α31x1x3α32x2x3q3E3x3=0

Then

x1*=N1D,x2*=N2D,x3*=N3D
where
N1=a1q1E1α22α33α23α32+a2q2E2α13α32α12α33+a3q3E3α12α23α13α22
N2=a1q1E1α31α23α21α33+a2q2E2α11α33α13α31+a3q3E3α11α23α13α21
N3=a1q1E1α31α22α21α32+a2q2E2α11α32α12α31+a3q3E3α11α22α12α21
D=α11α22α33α23α32+α22α21α33α31α23+α13α31α22α21α32

For x1*, x2*, x3* to be positive, the following inequalities hold:

α22α33>α23α32;α13α32>α12α33;α31α23>α21α33;α31α22>α21α32;α21α33>α31α32;α31α22>α21α32
a1q1E1α31α23α21α33+a2q2E2α11α33α13α31>a3q3E3α11α23α13α21
a1q1E1α31α22α21α32+a3q3E3α11α22α12α21>a2q2E2α11α32α12α31

4 Stability

4.1 Local stability analysis

To compute the local stability behaviour, we have studied the variational matrix corresponding to the interior equilibrium. The variational matrix of the system (1)–(3) at G4 (x1*, x2*, x3*) is:

A=α11α1*α12α1*α13α1*α21α2*α22α2*α23α2*α31α3*α32α3*α33α3*

We will point out here that, although (0,0,0) is defined for the system, it cannot be linearised there. So local stability of G1 (0,0,0) cannot be studied. At G2(x¯1,x¯2,0), the characteristic equation of A(x¯1,x¯2,0) is λ2+λα11+α22+α11α22+α12α21=0 since λ1+λ2=α11+α12<0 and λ1λ2 = α11α22 + α12α21 > 0. Thus G2(x¯1,x¯2,0) is locally asymptotically stable.

At G3 (x1ϕ, 0, x3ϕ), the characteristic equation of A (x1ϕ, x2ϕ, 0) is λ2+λα11+α33+α11α33+α31α13=0, since λ1+λ2=α11+α33<0 and λ1λ2 = α11α33 + α31α13 > 0.

Thus G3 (x1ϕ, 0, x3ϕ) is locally asymptotically stable.

The characteristic equation A (x1*, x2*, x3*) is

λ3+a1λ2+a2λ+a3=0(5)
where
a1=α33x3*+α22x2*+α11x1*>0
a2=(α22α33α23α32)x2*x3*+(α11α33+α13α31)x1*x3*+(α12α21+α11α22)x1*x2*>0
a3=α11α22α33+α13α31α22α13α21α32+α12α21α33+α12α23α31α11α23α32x1*x2*x3*

By the Routh-Hurwitz criterion, it follows that all eigen values of (5) will have negative real parts, if and only if,

a1>0;a3>0;a1a2a3>0

Clearly a1 > 0 and a3 > 0, and after some manipulation, it is very easy to check that a1a2 − a3 > 0. Hence G4(x1*, x2*, x3*) is locally asymptotically stable.

4.2 Global stability analysis

In this section, we have considered the global stability of the system (1)–(3) by constructing a suitable Lyapunov function.

Theorem 1. The equilibrium pointG2(x¯1,x¯2,0)is globally asymptotically stable.

Proof: Let us consider the following:

V(x1,x2)=x1x1*x1*lnx1x1*   +l1x2x2*x2*lnx2x2*
where l1 is a suitable constant.

Differentiating V with respect to time t along the solutions of model (1)–(3) and by choosing l1=α12α21, a little algebraic manipulation yields

dVdt<α11x1x1*2α12α22α21x2x2*2<0

This shows that dVdt is negative definite, and hence by the Lyapunov theorem on stability, it follows that the positive equilibrium G2(x¯1,x¯2,0) is globally asymptotically stable with respect to all solutions initiating in the interior of the positive quadrant.

Theorem 2. The equilibrium point G3(x1ϕ, 0, x3ϕ) is globally asymptotically stable.

Proof: Let us consider

V(x1,x3)=x1x1*x1*lnx1x1*+l2x3x3*x3*lnx3x3*
where l2 is a suitable constant. Differentiating V with respect to time t along the solutions of model (1)–(3) and by choosing l2=α13α31, a little algebraic manipulation yields
dVdt<α11x1x1ϕ2α13α33α31x3x3ϕ2<0

This shows that dVdt is negative definite, and hence by Lyapunov theorem on stability, it follows that the positive equilibrium G3(x1ϕ, 0, x3ϕ) is globally asymptotically stable with respect to all solutions initiating in the interior of the positive quadrant.

Theorem 3. The equilibrium point G4(x1*, x2*, x3*) is globally asymptotically stable.

Proof: Let us consider

V(x1,x2,x3)=x1x1*x1*lnx1x1*+l1x2x2*x2*lnx2x2*+l2x3x3*x3*lnx3x3*
where l1 and l2 are suitable constants to be determined in the subsequent steps. It can be easily verified that the function V is zero at the equilibrium point G4(x1*, x2*, x3*). Differentiating V with respect to time t along the solutions of model (1)–(3) and by choosing l1=α12α21, and l2=α13α31 with a little algebraic manipulation, we get dVdt<0. Since dVdt<0 in some neighbourhood (x1*, x2*, x3*), therefore the interior equilibrium point (x1*, x2*, x3*) is globally asymptotically stable.

5 Quantitative bionomic aspect of the model

It is mainly associated with study of the dynamics of living resources using economic models. The bionomic equilibrium is said to be achieved when the total revenue obtained by selling the harvested biomass equals the total cost utilized in harvesting. As we have already discussed, a biological equilibrium is given by dxidt=0, i = 1, 2, 3.

Let, c1 be the harvesting cost per unit effort for prey species, c2 be the harvesting cost per unit effort for first predator species (x2), c3 be the harvesting cost per unit effort for the second predator species (x3), p1 be the price per unit biomass of the prey, p2 be the price per unit biomass of the first predator species (x2), p3 be the price per unit biomass of the second predator species (x3).

Therefore, net revenue or economic rent at any time given by R = R1 + R2 + R3 where

R1=(p1q1x1c1)E1;R2=(p2q2x2c2)E2;R3=(p3q3x3c3)E3
here R1 represents Net Revenue for the prey; R2 represents Net Revenue for first predator species (x2); R3 represents Net Revenue for second predator species (x3).

The bionomic equilibrium (x1),(x2),(x3),E1,E2,E3 is given by the following equations:

a1x1α11x12α12x1x2α13x1x3q1E1x1=0(6)
a2x2α22x22α21x1x2α23x2x3q2E2x2=0(7)
a3x3α33x32α31x1x3α32x2x3q3E3x3=0(8)
R=p1q1x1c1E1+p2q2x2c2E2+p3q3x3c3E3=0(9)

In order to determine the bionomic equilibrium, we come across the following cases.

Case (i): if c1 > p1q1x1, c2 > p2q2x2, c3 > p3q3x3, then the cost is greater than revenue for three species and the whole system will be closed.

Case (ii): if c1 < p1q1x1, c2 < p2q2x2, c3 < p3q3x3, then the revenues for all the species being positive, then the whole system will be in operation.

x1=c1p1q1,x2=c2p2q2,x3=c3p3q3,

Now substitute x1,x2,x3 in Eqs. (6)–(9), we get

E1=1q1a1α11c1p1q1α12c2p2q2α13c3p3q3
E2=1q2a2α22c2p2q2+α21c1p1q1α23c3p3q3
E3=1q3a3α33c3p3q3+α31c1p3q3α32c2p2q2

Now (E1) > 0, if

a1>α11c1p1q1+α12c2p2q2+α13c3p3q3(10)

(E2) > 0, if

a2+α21c1p1q1>α22c2p2q2+α23c3p3q3(11)
and (E3) > 0, if
a3+α31c1p3q3>α33c3p3q3+α32c2p2q2(12)

Thus the Non-trivial bionomic equilibrium point (x1),(x2),(x3),E1,E2,E3 exists, if conditions (10)–(12) hold.

6 The Stochastic delayed model

Now, we have extended the above model (1)–(3) incorporation the effect of random fluctuations. The sensitive parameters of the system fluctuate about their average values due to these random fluctuations. We incorporated such randomness to the model by adding white noise. The main assumption that leads us to extend the deterministic model (1)–(3) to a stochastic counterpart is that it is reasonable to conceive the open sea as a noisy environment with chaos [30]. There are a number of ways in which environmental noise may be incorporated in system (1)–(3). Generally, the environmental noise is distinguished from demographic or internal noise, for the variation over time is due to different causes. External noise may arise either from random fluctuations of the parameters around some known mean values or from stochastic fluctuations of the population densities around some constant values.

In this section, we computed the population intensities of fluctuations, i.e., variances around the positive equilibrium G4 due to noise, according to the method introduced by Nisbet and Gurney [31]. Now, we assumed the presence of randomly fluctuating driving forces on the deterministic growth of the prey, predator-one and predator-two populations at time t, so that the system (1)–(3) is the following:

dx1dt=a1x1α11x12α12x1x2α13x1x3q1E1x1+η1ξ1t(13)
dx1dt=a2x2α22x22α21x1tτx2tτα23x2x3q2E2x2+η2ξ2t(14)
dx3dt=a3x3α33x32+α31x1tτx3tτα32x2x3q3E3x3+η3ξ3t(15)
where x1(t) represents prey species, x2(t) represents predator-one species, x3(t) represents predator-two species, α1, α2, α3 are real constants and ξ(t) = [ξ1(t), ξ2(t), ξ3(t)] is a three-dimensional Gaussian white noise process satisfying
Eξit=0;i=1,2,3(16)
Eξitξjt=δijδtt;i=j=1,2,3(17)
where δij is the Kronecker symbol; δ is the δ-dirac function.

Let

x1(t)=u1(t)+S*;x2(t)=u2(t)+P*;x3(t)=u3(t)+T*;(18)
dx1dt=du1(t)dt;dx2dt=du2(t)dt;dx3dt=du3(t)dt(19)

Using (18) and (19), Eq. (13) becomes

du1tdt=a1u1t+a1S*α11u12tα11S*22α11u1tS*α12u1tu2t+α12u1tP*α12u2tS*α12S*P*α13u1tu3tα13u1tT*α13u3tS*+α13S*T*q1E1u1tq1E1S*+η1ξ1t(20)

The linear part of (20) is

du1tdt=α11u1tS*α12u2tS*α13u3tS*+η1ξ1t(21)

Again using (18) and (19), Eq. (14) becomes

du2tdt=a2u2t+a2P*α22u22tα22P*22α22u2tP*+α21u1tτu2tτP*)2+α21u1tτP*+α21u2tτS*+α21S*P*α23u2tu3tα23u2tT*α23u3tP*+α23P*T*q2E2u2tq2E2P*+η2ξ2t(22)

The linear part of Eq. (22) is

du2(t)dt=α22u2(t)P*+α21u1(tτ)P*α23u3(t)P*+η2ξ2(t)(23)

Using Eqs. (18) and (19), Eq. (15) becomes

du3tdt=a3u3t+a3T*α33u32tα33T*22α33u3tT*+α31u1tτu3tτ+α31u1tτT*+α31u3tτS*+α31S*T*α32u2tu3tα32u2tT*α32u3tP*+α32P*T*q3E3u3tq3E3T*+η3ξ3t(24)

The linear part of (24) is

du3(t)dt=α33u3(t)T*+α31u1(tτ)T*α32u2(t)T*+η3ξ3(t)(25)

Taking the Fourier transform on both sides of (21), (23), (25) we get:

iωu˜1(ω)=α11S*u˜1(ω)α12S*u˜2(ω)α13S*u˜3(ω)+η1ξ˜1(ω)
(or)
η1ξ˜1(ω)=iω+α11S*u˜1(ω)+α12S*u˜2(ω)+α13S*u˜3(ω)iωu˜2(ω)=α22P*u˜2(ω)+α21P*eiωτu˜1(ω)α23P*u˜3(ω)+η2ξ˜2(ω)(26)
(or)
η2ξ˜2(ω)=α21P*eiωτu˜1(ω)+iω+α22P*u˜2(ω)+α23P*u˜3(ω)iωu˜3(ω)=α33T*u˜3(ω)+α31T*eiωτu˜1(ω)α32T*u˜2(ω)+η3ξ˜3(ω)(27)
(or)
η3ξ˜3(ω)=α31T*eiωτu˜1(ω)+α32T*u˜2(ω)+iω+α33T*u˜3(ω)(28)

The matrix form of Eqs. (26)–(28) is

Mωu˜ω=ξ˜ω(29)
where
Mω=A(ω)B(ω)C(ω)D(ω)E(ω)F(ω)G(ω)H(ω)I(ω);
u˜ω=u˜1(ω)u˜2(ω)u˜3(ω);ξ˜ω=ξ˜1ωξ˜2ωξ˜3ω
A(ω)=iω+α11S*;B(ω)=α12S*;C(ω)=α13S*;D(ω)=α21P*eiωτ;E(ω)=iω+α22P*;G(ω)=α31T*eiωt;H(ω)=α32T*;I(ω)=iω+α33T*(30)

Eq. (29) can also be written as

u˜ω=Mω1ξ˜ω

Let [M(ω)]−1 = K(ω), therefore

u˜ω=K(ω)ξ˜ω(31)
where
K(ω)=AdjMωMω(32)
if the function Y(t) has a zero mean value, then the fluctuation intensity (variance) of its components in the frequency interval [ω, ω + ] is SY(ω) where SY(ω) is spectral density of Y and is defined as
SY(ω)=limT˜Y˜ω2T˜(33)

If Y has a zero mean value, the inverse transform of SY(ω) is the autocovariance function

CYτ=12πSYωeiωτdω(34)

The corresponding variance of fluctuations in Y(t) is given by

σY2=CY0=12πSYωdω(35)
and the autocorrelation function is the normalized autocovariance
PYτ=CYτCY0(36)

For a Gaussian white noise process, it is

Sξiξjω=limTˆ+Eξ˜iωξ˜jωTˆ=limTˆ+1TˆTˆ2Tˆ2Tˆ2Tˆ2Eξ˜itξ˜jt'eiω(tt)dtdt=δij(37)

From (31), we have

u˜iω=j=13Kijωξ˜jω;i=1,2,3(38)

From (33) we have

Suiω=j=13ηjKijω2;i=1,2,3(39)

Hence by (35) and (39), the intensities of fluctuations in the variable ui; i = 1, 2, 3 are given by

σui2=12πj=13ηjKijω2dω;i=1,2,3(40)
and by (31) and (32), we obtain three variances of ui (i = 1, 2, 3) of the model system (13)–(15) as follows:
σu12=12πη1Adj(1)M(ω)2+η2Adj(2)M(ω)2+η3Adj(3)M(ω)2dωσu22=12πη1Adj(4)M(ω)2+η2Adj(5)M(ω)2+η3Adj(6)M(ω)2dωσu32=12πη1Adj(7)M(ω)2+η2Adj(8)M(ω)2+η3Adj(9)M(ω)2dω(41)
where |M(ω)| = |R(ω)| + i|I(ω)|.

The real part:

M(ω)=R2ω=α33ω2T*ω2α22P*ω2α11S*+α11α22α33S*P*α11α23α32S*T*P*+α12α21α33S*P*T*cosωτ+ωα12α21S*P*sinωτα12α31α23S*T*P*cosωτα13α21α32S*T*P*cosωτ+ωα13α31S*T*sinωτ+α13α31α22S*P*T*(42)

The imaginary part:

M(ω)=I2ω=ω3+ωα22α33P*T*+ωα11α33S*T*+ωα11α22S*P*ωα23α32T*P*+ωα12α21S*P*cosωτα12α21α33S*P*T*sinωτ+α12α31α23S*T*P*sinωτ+α13α21α32S*T*P*sinωτ+ωα13α31S*T*cosωτα22α13α31S*P*T*sinωτ(43)
and
Adj(1)2=X12+Y12;Adj(2)2=X22+Y22;Adj(3)2=X32+Y32;Adj(4)2=X42+Y42;Adj(5)2=X52+Y52;Adj(6)2=X62+Y62;Adj(7)2=X72+Y72;Adj(8)2=X82+Y82;Adj(9)2=X92+Y92
where
X1=ω2+α22α23P*T*α23α32T*P*
Y1=ωα33T*+ωα22P*
X2=α21α33P*T*cosωτα23α31T*P*cosωτ+ωα21P*sinωτ
Y2=ωα21P*cosωτωα21α33P*T*sinωτ+α23α31T*P*
X3=α21α32P*T*cosωτ+α31α22T*P*cosωτ+ωα31T*sinωτ
Y3=ωα31T*cosωτ+α21α32P*T*sinωτα31α22T*P*sinωτ
X4=α12α13S*T*α32α13S*T*;Y4=ωα12S*
X5=ω2+α11α33S*T*+α31α13S*T*cosωτ
Y5=ωα33T*+ωα11S*α31α13S*T*sinωτ
X6=α11α32S*T*+α31α12S*T*cosωτ
Y6=ωα32T*α31α12S*T*sinωτ
X7=α12α23S*P*α22α13S*P*;Y7=ωα13S*
X8=α12α23S*P*+α21α13S*P*cosωτ
Y8=ωα23P*α21α13S*P*
X9=ω2+α11α22S*P*+α12α21S*P*cosωτ
Y9=ωα22P*+ωα11S*α12α21S*P*sinωτ

Thus Eq. (41) becomes

σu12=12π1R2(ω)+I2(ω)η1X12+Y12+η2X42+Y42+η3X72+Y72dω
σu12=12π1R2(ω)+I2(ω)η1X22+Y22+η2X52+Y52+η3X82+Y82dω
σu12=12π1R2(ω)+I2(ω)η1X32+Y32+η2X62+Y62+η3X92+Y92dω

If we are interested in the dynamics of system (13)–(15) with either η1 = 0 or η2 = 0 or η3 = 0, then the population variances are the following:

(i) if η1 = 0, η2 = 0, then

σu12=η3X72+Y722π1R2(ω)+I2(ω)dω
σu22=η3X82+Y822π1R2(ω)+I2(ω)dω
σu32=η3X92+Y922π1R2(ω)+I2(ω)dω

(ii) if η2 = 0, η3 = 0, then

σu12=η1X12+Y122π1R2(ω)+I2(ω)dω
σu22=η1X22+Y222π1R2(ω)+I2(ω)dω
σu32=η1X32+Y322π1R2(ω)+I2(ω)dω

(iii) if η3 = 0, η1 = 0, then

σu12=η2X42+Y422π1R2(ω)+I2(ω)dω
σu22=η2X52+Y522π1R2(ω)+I2(ω)dω
σu32=η2X62+Y622π1R2(ω)+I2(ω)dω

The expressions in (41) gives three variances of the populations. The integrations over the real line can be evaluated which gives the variances of the populations.

7 Numerical simulation results

In this section, we substantiate as well as augment our analytical findings through numerical simulations considering the following examples.

7.1 Example 1

We take the following parameters:

x1=10;x2=10;x3=20;a1=8;α11=0.05;α12=0.6;α13=0.7;q1=0.15;E1=10;a2=2;α22=0.7;α21=0.17;α23=0.13;q2=0.01;E2=15;a3=1;α33=0.4;α31=0.15;α32=0.18;q3=0.03;E3=5;

Fig. 1 shows that the variation of population against time, initially with x1 = 10; x2 = 10; x3 = 20; in three different graphs.

Fig. 1

Results of Example 1.

7.2 Example 2

We take the following parameters:

x1=10;x2=10;x3=20;a1=8;α11=0.05;α12=0.6;α13=0.7;q1=0.15;E1=10;a2=2;α22=0.7;α21=0.17;α23=0.13;q2=0.01;E2=15;a3=1;α33=0.4;α31=0.15;α32=0.18;q3=0.03;E3=5;

Fig. 2 shows that the variation of population against time, initially with x1 = 10; x2 = 10; x3 = 20; in a single graph.

Fig. 2

Results of Example 2.

7.3 Example 3

We take the following parameters:

x1=10;x2=15;x3=20;a1=8;α11=0.05;α12=0.6;α13=0.7;q1=0.15;E1=10;a2=2;α22=0.7;α21=0.17;α23=0.13;q2=0.01;E2=15;a3=1;α33=0.4;α31=0.15;α32=0.18;q3=0.03;E3=5;

Fig. 3 shows that the variation of population amongst three species prey, predator-one, and predator-two initially with x1 = 10; x2 = 15; x3 = 20; in three dimensions.

Fig. 3

Results of Example 3.

Fig. 4 shows that the population oscillation gives chaos against time under random environmental noise of three species prey, predator-one, and predator-two with initial population sizes x1 = 10; x2 = 15; x3 = 20.

Fig. 4

Example 3, population oscillation gives chaos against time.

Fig. 5 shows that the variation of population against time under random environmental noise of species prey, predator-one and predator-two and the fourth is the phase portrait of the same initially with x1 = 10; x2 = 15; x3 = 20.

Fig. 5

Example 3, variation of population against time under random environmental noise.

8 Concluding remarks

In this article, a mathematical model has been proposed and analysed to study the dynamics of fishery resources. It is assumed that fish populations are growing logistically in reality. Using the stability theory of ordinary differential equations [32], it is known that the interior equilibrium should exist under certain conditions and it is globally asymptotically stable.

From Figs. 1–3, we may conclude that the given biological system is stable and in the case of fish (in particular fish species like goldband fish as a prey, shark as a first predator and baleen whale or dolphin as second predator), Eqs. (10)–(12) are satisfied, then the fishermen will get profit and further, they can run fisheries. Usually fishermen may harvest the predator population to increase profit even though harvesting of predator species is very costly and sometimes difficult [33–35]. However, in the case of prey species fishermen will not get that much profit when compared to predators. Even though, some fishermen will concentrate on the harvesting of prey species, since prey fishes (like goldband small fishes) play a major vital role while curing of the diseases related to eye. Also from the Figs. 4 and 5, we have seen that the system (13)–(15) is chaotic [29,30], in nature as it is a-periodic along with high sensitivity to the initial conditions for solution trajectories in bounded region. Recently, ratio-dependent system interactions have attracted many researchers, as these models produce richer dynamics in real media. The dynamics of ratio-dependent fishery models are also an important area of research which is left for future investigations for good innovative researcher.

Disclosure of interest

The authors declare that they have no conflicts of interest concerning this article.

Acknowledgements

The authors are grateful to the reviewers for their critical evaluations and suggestions.

Appendix A Proof of Theorem 3

Let us consider the Lyapunov function

V(x1,x2,x3)=x1x1*x1*lnx1x1*+l1x2x2*x2*lnx2x2*+l2x3x3*x3*lnx3x3*
where l1 and l2 are suitable constants to be determined in the subsequent steps. It can be easily verified that the function V is zero at the equilibrium point G4(x1*, x2*, x3*). The time derivative of V along the trajectories of (1)–(3) is
dVdt=x1x1*x1dx1dt+l1x2x2*x2dx2dt+l2x3x3*x3dx3dt

Choosing l1=α12α21, l2=α13α31 after a little algebraic manipulation, we get,

dVdt=α11x1x1*2α12α22α21x2x2*2α33α13α31x3x3*2α23α12α21x2x2*x3x3*α13α32α31x2x2*x3x3*
dVdt<α11x1x1*2α12α22α21x2x2*2α33α13α31x3x3*2α23α12α21+α13α32α31x2x2*22+x3x3*22

Now it is very clear that dVdt<0.

Since dVdt<0 in some neighbourhood (x1*, x2*, x3*) therefore the interior equilibrium point (x1*, x2*, x3*) is globally asymptotically stable.


References

[1] T.K. Kar; K. Chakraborty Bioeconomic modelling of a prey predator system using differential algebraic equations, Int. J. Eng. Sci. Technol., Volume 2 (2010) no. 1, pp. 13-34

[2] A. Hastings Global stability of two species systems, J. Math. Boil., Volume 5 (1978), pp. 399-403

[3] X.Z. He Stability and delays in a predator-prey system, J. Math. Anal. Appl., Volume 198 (1996), pp. 355-370

[4] B.S. Goh Global stability in two species interactions, J. Math. Biol., Volume 3 (1976), pp. 313-318

[5] A.J. Lokta Elements of Physical biology, Williams and Wilkins, Baltimore, 1925

[6] V. Volterra Leçons sur la théorie mathématique de la lutte pour la vie, Gauthier-Villars, Paris, 1931

[7] W.J. Meyer, Concepts of Mathematical Modelling, Mc Graw-Hill, 1985.

[8] J.M. Cushing Integro-differential equations and Delay models in Population dynamics, Lecture Notes in Biomathematics, vol. 20, Springer-Verlag, Heidelberg, 1997

[9] P. Colinvaux Ecology, John Wiley and Sons Inc., New York, 1986

[10] H.I. Freedman Deterministic Mathematical Models in Population Ecology, Marcel Dekker, New York, 1980

[11] J.N. Kapur, Mathematical Modelling, Weley-Eastern, 1998.

[12] J.N. Kapur, Mathematical Models in Biology and Medicine Affiliated, East-West, 1985.

[13] K.S. Chaudhuri A bio economic mode of harvesting of a multi species fishery, Ecol. Model, Volume 32 (1986), pp. 267-279

[14] T.K. Kar; M. Swarnakamal Influence of prey reserve in a prey-predator fishery, Non-Linear Anal., Volume 65 (2006), pp. 1725-1735

[15] C.W. Clark Mathematical bioeconomics: the optimal management of renewable resources, Wiley, New York, 1976

[16] R. Zhang; S. Junfang; Y. Haixia Analysis of a prey-predator fishery model with prey reserve, Appl. Math. Sci., Volume 1 (2007) no. 50, pp. 2481-2492

[17] S.M. Garcia, C. Newton, Current situation, trends and prospects in world capture fisheries, in: E.K. Pikitch, D.D. Huppert, M.P. Sissenwine (Eds.), Global trends in fisheries management, American Fisheries Society Symposium, 20 Bethesda MD, 1997, pp. 3–27.

[18] F. Brauer; A.C. Soudack Stability regions and transition phenomena for harvested predator–prey systems, J. Math. Biol., Volume 7 (1979), pp. 319-337

[19] F. Brauer; A.C. Soudack Stability regions in predator–prey systems with constant rate prey harvesting, J. Math. Biol., Volume 8 (1979), pp. 55-71

[20] G. Dai; M. Tang Coexistence region and global dynamics of a harvested predator–prey systems, SIAM J. Appl. Math., Volume 58 (1998), pp. 193-210

[21] M.R. Myerscough; B.F. Gray; W.L. Hogarth; J. Norbury An analysis of an ordinary differential equations model for a two species predator–prey system with harvesting and stocking, J. Math. Biol., Volume 30 (1992), pp. 389-411

[22] T.K. Kar; K.S. Chaudhuri On non-selective harvesting of a multispecies fishery, Int. J. Math. Educ. Technol., Volume 33 (2002) no. 4, pp. 543-556

[23] T.K. Kar; K.S. Chaudhuri Regulation of a prey–predator fishery by taxation: a dynamic reaction models, J. Biol. Syst., Volume 11 (2003) no. 2, pp. 173-187

[24] N.H. Gazi; K. Das; Z. Mukandavire; C. Chiyaka; P. Das A study of Cholera model with environmental fluctuations, Int. J. Math. Models Methods Appl. Sci., Volume 3 (2010) no. 4, pp. 150-155

[25] B. Ravindra Reddy; K. Lakshmi Narayan; N.C. Pattabhiramacharyulu A model of two mutually interacting species with limited resources for both the species, Int. J. Eng. Res. Indu. Appls, Volume 2 (2009) no. II, pp. 281-291

[26] K. Shiva Reddy; K. Lakshmi Narayan; N.C. Pattabhiramacharyulu A Three Species Eco system consisting of A Prey And Two Predators, Int. J. Math. Sci. Eng. Appls, Volume 14 (2010) no. 4, pp. 129-145

[27] K. Shiva Reddy; K. Lakshmi Narayan; N.C. Pattabhiramacharyulu A three species ecosystem consisting of a prey, predator and super predator, Math. Appl. Sci. Technol., Volume 2 (2010) no. 1, pp. 95-107

[28] M. Carletti Numerical simulation of a Campbell-like stochastic delay model for bacteriophage infection, Math. Med. Biol., Volume 23 (2006), pp. 297-310

[29] R.L. Deveney, An Introduction to Chaotic Dynamical Systems, Wiseley Publishing Company, Inc., 1986.

[30] M.J. Stutzer Chaotic dynamics and bifurcation in a macro model, J. Econ. Dyn. Control, Volume 2 (1980), pp. 353-376

[31] R.M. Nisbet; W.S.C. Gurney The systematic formulation of population models for insects with dynamically varying instar duration, Theoret. Popul. Biol., Volume 23 (1982) no. 1, pp. 114-135

[32] G.F. Simmons Differential Equations with Applications and Historical Notes, Tata McGraw-Hill, New Delhi, 1974

[33] T.K. Kar; K.S. Chaudhuri Harvesting in a two-prey one-predator fishery: a bioeconomic model, ANZIAM J., Volume 45 (2004), pp. 443-456

[34] D. Purohit; K.S. Chaudhuri A bioeconomic model of non selective harvesting of two competing fish species, ANZIAM J., Volume 46 (2004), pp. 299-308

[35] K. Shiva Reddy; N.C. Pattabhiramacharyulu A three species ecosystem consisting of two predators competing for a prey, J. Adv. Appl. Sci. Res., Volume 2 (2011) no. 3, pp. 208-218


Comments - Policy