Chaotic dynamics of a three species prey–predator competition model with bionomic harvesting due to delayed environmental noise as external driving force

Kalyan Das; M.N. Srinivas; M.A.S. Srinivas; N.H. Gazi

doi:10.1016/j.crvi.2012.06.001

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

Kalyan Das ¹ ; M.N. Srinivas ² ; M.A.S. Srinivas ³ ; N.H. Gazi ⁴

¹ National Institute of Food Technology Entrepreneurship and Management, Department of Mathematics, Plot No. 97, Sector 56, HSIIDC Industrial Estate, Kundli 131 028, Haryana, India
² School of Advanced Sciences, Department of Mathematics, VIT University, Vellore 632 014, Tamil Nadu, India
³ Department of Mathematics, Jawaharlal Nehru Technological University, Hyderabad, Andhra Pradesh, India
⁴ Department of Mathematics, Aliah University, DN-41, Salt Lake, Sector V, Kolkata 700 091, India

Comptes Rendus. Biologies, Volume 335 (2012) no. 8, pp. 503-513.

Résumé

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.

Métadonnées

Reçu le : 2012-01-23
Accepté le : 2012-06-04
Publié le : 2012-08-01

PMID

DOI : 10.1016/j.crvi.2012.06.001

Mots clés : Prey–predator, Local and global stability, Bionomic equilibrium, Stochastic delayed perturbation, Fourier transform methods, Chaos

Affiliations des auteurs :

Kalyan Das ¹ ; M.N. Srinivas ² ; M.A.S. Srinivas ³ ; N.H. Gazi ⁴

@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 (S₁) common to two predators (S₂) and (S₃) 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:

\frac{d x_{1}}{d t} = a_{1} x_{1} - α_{11} x_{1}^{2} + α_{12} x_{1} x_{2} - α_{13} x_{1} x_{3} - q_{1} E_{1} x_{1}

(1)

\frac{d x_{2}}{d t} = a_{2} x_{2} - α_{22} x_{2}^{2} - α_{21} x_{1} x_{2} - α_{23} x_{2} x_{3} - q_{2} E_{2} x_{2}

(2)

\frac{d x_{3}}{d t} = a_{3} x_{3} - α_{33} x_{3}^{2} + α_{31} x_{1} x_{3} - α_{32} x_{2} x_{3} - q_{3} E_{3} x_{3}

(3)

with the following notation:

• x_i(t): population density of the species S_i at time t; i = 1, 2, 3;
• a_i: the natural growth rates of S_i, i = 1, 2, 3;
• α_ii: the rate of decrease of S_i 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;
• K_i = a_i/α_ii: carrying capacities of species, i = 1, 2, 3.

Further, the variables x₁, x₂ and x₃ are non-negative and the model parameters a_i, K_i, α_ij, i = 1, 2, 3, j = 1, 2, 3.

• q₁: catch ability coefficient of prey species;
• q₂: catch ability coefficient of predator species;
• q₃: catch ability coefficient of predator species;
• E₁: effort applied to harvest the prey species (S₁);
• E₂: effort applied to harvest the predator species (S₂);
• E₃: effort applied to harvest the predator species (S₃).

Throughout our analysis, we take the initial conditions as

\begin{array}{l} a_{2} - q_{2} E_{2} > 0, \\ a_{3} - q_{3} E_{3} > 0, \\ a_{1} - q_{1} E_{1} > 0 \end{array}

(4)

3 Existence of equilibrium points

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

\begin{array}{l} 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 \end{array}

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

• the trivial state G₁ (0,0,0) (in the absence of all the species);
• the axial state $G_{2} ({\bar{x}}_{1}, {\bar{x}}_{2}, 0)$ (in the absence of second predator);
• the boundary state $G_{3} (x_{1}^{ϕ}, 0, x_{3}^{ϕ})$ (in the absence of first predator);
• the steady state of coexistence G₄ (x₁^*, x₂^*, x₃^*) (the interior equilibrium).

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

Case (ii): if ${\bar{x}}_{1}$ and ${\bar{x}}_{2}$ are the positive solutions of

\begin{array}{l} 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 \end{array}

Then

\begin{array}{l} {\bar{x}}_{1} = \frac{α_{22} (a_{1} - q_{1} E_{1}) - α_{12} (a_{2} - q_{2} E_{2})}{α_{21} α_{12} + α_{11} α_{22}} \\ {\bar{x}}_{2} = \frac{α_{21} (a_{1} - q_{1} E_{1}) + α_{11} (a_{2} - q_{2} E_{2})}{α_{21} α_{12} + α_{11} α_{22}} > 0 \end{array}

${\bar{x}}_{1}$ is positive provided the following inequality holds:

α_{22} (a_{1} - q_{1} E_{1}) > α_{12} (a_{2} - q_{2} E_{2})

Case (iii): if x₁^ϕ and x₃^ϕ are the positive solutions of

\begin{array}{l} 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_{3} x_{3} - α_{33} x_{3}^{2} + α_{31} x_{1} x_{3} - α_{32} x_{2} x_{3} - q_{3} E_{3} x_{3} = 0 \end{array}

Then

{x_{1}}^{ϕ} = \frac{α_{33} (a_{1} - q_{1} E_{1}) + α_{13} (a_{3} - q_{3} E_{3})}{α_{13} α_{31} + α_{11} α_{33}}

{x_{3}}^{ϕ} = \frac{α_{31} (a_{1} - q_{1} E_{1}) + α_{11} (a_{3} - q_{3} E_{3})}{α_{13} α_{31} + α_{11} α_{33}}

x₁^ϕ is positive provided the following inequality holds:

α_{33} (a_{1} - q_{1} E_{1}) > α_{13} (a_{3} - q_{3} E_{3})

Case (iv): if x₁^* and x₂^* are the positive solutions of

\begin{array}{l} 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 \end{array}

Then

{x_{1}}^{*} = \frac{N_{1}}{D}, {x_{2}}^{*} = \frac{N_{2}}{D}, {x_{3}}^{*} = \frac{N_{3}}{D}

where

\begin{array}{l} N_{1} = (a_{1} - q_{1} E_{1}) [α_{22} α_{33} - α_{23} α_{32}] \\ + (a_{2} - q_{2} E_{2}) [α_{13} α_{32} - α_{12} α_{33}] \\ + (a_{3} - q_{3} E_{3}) [α_{12} α_{23} - α_{13} α_{22}] \end{array}

\begin{array}{l} N_{2} = (a_{1} - q_{1} E_{1}) [α_{31} α_{23} - α_{21} α_{33}] \\ + (a_{2} - q_{2} E_{2}) [α_{11} α_{33} - α_{13} α_{31}] \\ + (a_{3} - q_{3} E_{3}) [α_{11} α_{23} - α_{13} α_{21}] \end{array}

\begin{array}{l} N_{3} = (a_{1} - q_{1} E_{1}) [α_{31} α_{22} - α_{21} α_{32}] \\ + (a_{2} - q_{2} E_{2}) [α_{11} α_{32} - α_{12} α_{31}] \\ + (a_{3} - q_{3} E_{3}) [α_{11} α_{22} - α_{12} α_{21}] \end{array}

\begin{array}{l} D = α_{11} [α_{22} α_{33} - α_{23} α_{32}] \\ + α_{22} [α_{21} α_{33} - α_{31} α_{23}] \\ + α_{13} [α_{31} α_{22} - α_{21} α_{32}] \end{array}

For x₁^*, x₂^*, x₃^* to be positive, the following inequalities hold:

\begin{array}{l} α_{22} α_{33} > α_{23} α_{32}; α_{13} α_{32} > α_{12} α_{33}; \\ α_{31} α_{23} > α_{21} α_{33}; α_{31} α_{22} > α_{21} α_{32}; \\ α_{21} α_{33} > α_{31} α_{32}; α_{31} α_{22} > α_{21} α_{32} \end{array}

(a_{1} - q_{1} E_{1}) [α_{31} α_{23} - α_{21} α_{33}] + (a_{2} - q_{2} E_{2}) [α_{11} α_{33} - α_{13} α_{31}] > (a_{3} - q_{3} E_{3}) [α_{11} α_{23} - α_{13} α_{21}]

(a_{1} - q_{1} E_{1}) [α_{31} α_{22} - α_{21} α_{32}] + (a_{3} - q_{3} E_{3}) [α_{11} α_{22} - α_{12} α_{21}] > (a_{2} - q_{2} E_{2}) [α_{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 G₄ (x₁^*, x₂^*, x₃^*) is:

A = [\begin{array}{c} - α_{11} {α_{1}}^{*} & - α_{12} {α_{1}}^{*} & - α_{13} {α_{1}}^{*} \\ - α_{21} {α_{2}}^{*} & - α_{22} {α_{2}}^{*} & - α_{23} {α_{2}}^{*} \\ - α_{31} {α_{3}}^{*} & - α_{32} {α_{3}}^{*} & - α_{33} {α_{3}}^{*} \end{array}]

We will point out here that, although (0,0,0) is defined for the system, it cannot be linearised there. So local stability of G₁ (0,0,0) cannot be studied. At $G_{2} ({\bar{x}}_{1}, {\bar{x}}_{2}, 0)$ , the characteristic equation of $A ({\bar{x}}_{1}, {\bar{x}}_{2}, 0)$ is $λ^{2} + λ (α_{11} + α_{22}) + (α_{11} α_{22} + α_{12} α_{21}) = 0$ since $λ_{1} + λ_{2} = - (α_{11} + α_{12}) < 0$ and λ₁λ₂ = α₁₁α₂₂ + α₁₂α₂₁ > 0. Thus $G_{2} ({\bar{x}}_{1}, {\bar{x}}_{2}, 0)$ is locally asymptotically stable.

At G₃ (x₁^ϕ, 0, x₃^ϕ), the characteristic equation of A (x₁^ϕ, x₂^ϕ, 0) is $λ^{2} + λ (α_{11} + α_{33}) + (α_{11} α_{33} + α_{31} α_{13}) = 0$ , since $λ_{1} + λ_{2} = - (α_{11} + α_{33}) < 0$ and λ₁λ₂ = α₁₁α₃₃ + α₃₁α₁₃ > 0.

Thus G₃ (x₁^ϕ, 0, x₃^ϕ) is locally asymptotically stable.

The characteristic equation A (x₁^*, x₂^*, x₃^*) is

λ^{3} + a_{1} λ^{2} + a_{2} λ + a_{3} = 0

(5)

where

a_{1} = α_{33} {x_{3}}^{*} + α_{22} {x_{2}}^{*} + α_{11} {x_{1}}^{*} > 0

\begin{array}{l} a_{2} = (α_{22} α_{33} - α_{23} α_{32}) {x_{2}}^{*} {x_{3}}^{*} \\ + (α_{11} α_{33} + α_{13} α_{31}) {x_{1}}^{*} {x_{3}}^{*} \\ + (α_{12} α_{21} + α_{11} α_{22}) {x_{1}}^{*} {x_{2}}^{*} > 0 \end{array}

a_{3} = (\begin{array}{l} α_{11} α_{22} α_{33} \\ + α_{13} α_{31} α_{22} \\ - α_{13} α_{21} α_{32} \\ + α_{12} α_{21} α_{33} \\ + α_{12} α_{23} α_{31} \\ - α_{11} α_{23} α_{32} \end{array}) {x_{1}}^{*} {x_{2}}^{*} {x_{3}}^{*}

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

a_{1} > 0; a_{3} > 0; a_{1} a_{2} - a_{3} > 0

Clearly a₁ > 0 and a₃ > 0, and after some manipulation, it is very easy to check that a₁a₂ − a₃ > 0. Hence G₄(x₁*, x₂*, x₃*) 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 point $G_{2} ({\bar{x}}_{1}, {\bar{x}}_{2}, 0)$ is globally asymptotically stable.

Proof: Let us consider the following:

\begin{array}{l} V (x_{1}, x_{2}) \\ = (x_{1} - {x_{1}}^{*} - {x_{1}}^{*} \ln (\frac{x_{1}}{{x_{1}}^{*}})) \\ + l_{1} (x_{2} - {x_{2}}^{*} - {x_{2}}^{*} \ln (\frac{x_{2}}{{x_{2}}^{*}})) \end{array}

where l₁ is a suitable constant.

Differentiating V with respect to time t along the solutions of model (1)–(3) and by choosing $l_{1} = \frac{α_{12}}{α_{21}}$ , a little algebraic manipulation yields

\frac{d V}{d t} < - α_{11} {(x_{1} - {x_{1}}^{*})}^{2} - \frac{α_{12} α_{22}}{α_{21}} {(x_{2} - {x_{2}}^{*})}^{2} < 0

This shows that $\frac{d V}{d t}$ is negative definite, and hence by the Lyapunov theorem on stability, it follows that the positive equilibrium $G_{2} ({\bar{x}}_{1}, {\bar{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 G₃(x₁^ϕ, 0, x₃^ϕ) is globally asymptotically stable.

Proof: Let us consider

\begin{array}{l} V (x_{1}, x_{3}) \\ = (x_{1} - {x_{1}}^{*} - {x_{1}}^{*} \ln (\frac{x_{1}}{{x_{1}}^{*}})) + l_{2} (x_{3} - {x_{3}}^{*} - {x_{3}}^{*} \ln (\frac{x_{3}}{{x_{3}}^{*}})) \end{array}

where l₂ is a suitable constant. Differentiating V with respect to time t along the solutions of model (1)–(3) and by choosing

l_{2} = \frac{α_{13}}{α_{31}}

, a little algebraic manipulation yields

\begin{array}{l} \frac{d V}{d t} < - α_{11} {(x_{1} - {x_{1}}^{ϕ})}^{2} \\ - \frac{α_{13} α_{33}}{α_{31}} {(x_{3} - {x_{3}}^{ϕ})}^{2} < 0 \end{array}

This shows that $\frac{d V}{d t}$ is negative definite, and hence by Lyapunov theorem on stability, it follows that the positive equilibrium G₃(x₁^ϕ, 0, x₃^ϕ) is globally asymptotically stable with respect to all solutions initiating in the interior of the positive quadrant.

Theorem 3. The equilibrium point G₄(x₁*, x₂*, x₃*) is globally asymptotically stable.

Proof: Let us consider

\begin{array}{l} V (x_{1}, x_{2}, x_{3}) \\ = (x_{1} - {x_{1}}^{*} - {x_{1}}^{*} \ln (\frac{x_{1}}{{x_{1}}^{*}})) \\ + l_{1} (x_{2} - {x_{2}}^{*} - {x_{2}}^{*} \ln (\frac{x_{2}}{{x_{2}}^{*}})) \\ + l_{2} (x_{3} - {x_{3}}^{*} - {x_{3}}^{*} \ln (\frac{x_{3}}{{x_{3}}^{*}})) \end{array}

where l₁ and l₂ 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 G₄(x₁*, x₂*, x₃*). Differentiating V with respect to time t along the solutions of model (1)–(3) and by choosing

l_{1} = \frac{α_{12}}{α_{21}}

, and

l_{2} = \frac{α_{13}}{α_{31}}

with a little algebraic manipulation, we get

\frac{d V}{d t} < 0

. Since

\frac{d V}{d t} < 0

in some neighbourhood (x₁*, x₂*, x₃*), therefore the interior equilibrium point (x₁*, x₂*, x₃*) 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 $\frac{d x_{i}}{d t} = 0$ , i = 1, 2, 3.

Let, c₁ be the harvesting cost per unit effort for prey species, c₂ be the harvesting cost per unit effort for first predator species (x₂), c₃ be the harvesting cost per unit effort for the second predator species (x₃), p₁ be the price per unit biomass of the prey, p₂ be the price per unit biomass of the first predator species (x₂), p₃ be the price per unit biomass of the second predator species (x₃).

Therefore, net revenue or economic rent at any time given by R = R₁ + R₂ + R₃ where

\begin{array}{l} R_{1} = (p_{1} q_{1} x_{1} - c_{1}) E_{1}; \\ R_{2} = (p_{2} q_{2} x_{2} - c_{2}) E_{2}; \\ R_{3} = (p_{3} q_{3} x_{3} - c_{3}) E_{3} \end{array}

here R₁ represents Net Revenue for the prey; R₂ represents Net Revenue for first predator species (x₂); R₃ represents Net Revenue for second predator species (x₃).

The bionomic equilibrium $({(x_{1})}_{\infty}, {(x_{2})}_{\infty}, {(x_{3})}_{\infty}, {(E_{1})}_{\infty}, {(E_{2})}_{\infty}, {(E_{3})}_{\infty})$ is given by the following equations:

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

(6)

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

(7)

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

(8)

\begin{array}{l} R = (p_{1} q_{1} x_{1} - c_{1}) E_{1} + (p_{2} q_{2} x_{2} - c_{2}) E_{2} \\ + (p_{3} q_{3} x_{3} - c_{3}) E_{3} = 0 \end{array}

(9)

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

Case (i): if c₁ > p₁q₁x₁, c₂ > p₂q₂x₂, c₃ > p₃q₃x₃, then the cost is greater than revenue for three species and the whole system will be closed.

Case (ii): if c₁ < p₁q₁x₁, c₂ < p₂q₂x₂, c₃ < p₃q₃x₃, then the revenues for all the species being positive, then the whole system will be in operation.

{(x_{1})}_{\infty} = \frac{c_{1}}{p_{1} q_{1}}, {(x_{2})}_{\infty} = \frac{c_{2}}{p_{2} q_{2}}, {(x_{3})}_{\infty} = \frac{c_{3}}{p_{3} q_{3}},

Now substitute ${(x_{1})}_{\infty}, {(x_{2})}_{\infty}, {(x_{3})}_{\infty}$ in Eqs. (6)–(9), we get

{(E_{1})}_{\infty} = \frac{1}{q_{1}} [a_{1} - \frac{α_{11} c_{1}}{p_{1} q_{1}} - \frac{α_{12} c_{2}}{p_{2} q_{2}} - \frac{α_{13} c_{3}}{p_{3} q_{3}}]

{(E_{2})}_{\infty} = \frac{1}{q_{2}} [a_{2} - \frac{α_{22} c_{2}}{p_{2} q_{2}} + \frac{α_{21} c_{1}}{p_{1} q_{1}} - \frac{α_{23} c_{3}}{p_{3} q_{3}}]

{(E_{3})}_{\infty} = \frac{1}{q_{3}} [a_{3} - \frac{α_{33} c_{3}}{p_{3} q_{3}} + \frac{α_{31} c_{1}}{p_{3} q_{3}} - \frac{α_{32} c_{2}}{p_{2} q_{2}}]

Now (E₁)_∞ > 0, if

a_{1} > \frac{α_{11} c_{1}}{p_{1} q_{1}} + \frac{α_{12} c_{2}}{p_{2} q_{2}} + \frac{α_{13} c_{3}}{p_{3} q_{3}}

(10)

(E₂)_∞ > 0, if

a_{2} + \frac{α_{21} c_{1}}{p_{1} q_{1}} > \frac{α_{22} c_{2}}{p_{2} q_{2}} + \frac{α_{23} c_{3}}{p_{3} q_{3}}

(11)

and (E₃)_∞ > 0, if

a_{3} + \frac{α_{31} c_{1}}{p_{3} q_{3}} > \frac{α_{33} c_{3}}{p_{3} q_{3}} + \frac{α_{32} c_{2}}{p_{2} q_{2}}

(12)

Thus the Non-trivial bionomic equilibrium point $({(x_{1})}_{\infty}, {(x_{2})}_{\infty}, {(x_{3})}_{\infty}, {(E_{1})}_{\infty}, {(E_{2})}_{\infty}, {(E_{3})}_{\infty})$ 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 G₄ 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:

\begin{array}{l} \frac{d x_{1}}{d t} = a_{1} x_{1} - α_{11} {x_{1}}^{2} - α_{12} x_{1} x_{2} - α_{13} x_{1} x_{3} \\ - q_{1} E_{1} x_{1} + η_{1} ξ_{1} (t) \end{array}

(13)

\begin{array}{l} \frac{d x_{1}}{d t} = a_{2} x_{2} - α_{22} {x_{2}}^{2} - α_{21} x_{1} (t - τ) x_{2} (t - τ) \\ - α_{23} x_{2} x_{3} - q_{2} E_{2} x_{2} + η_{2} ξ_{2} (t) \end{array}

(14)

\begin{array}{l} \frac{d x_{3}}{d t} = a_{3} x_{3} - α_{33} {x_{3}}^{2} + α_{31} x_{1} (t - τ) x_{3} (t - τ) \\ - α_{32} x_{2} x_{3} - q_{3} E_{3} x_{3} + η_{3} ξ_{3} (t) \end{array}

(15)

where x₁(t) represents prey species, x₂(t) represents predator-one species, x₃(t) represents predator-two species, α₁, α₂, α₃ are real constants and ξ(t) = [ξ₁(t), ξ₂(t), ξ₃(t)] is a three-dimensional Gaussian white noise process satisfying

E [ξ_{i} (t)] = 0; i = 1,2,3

(16)

E [ξ_{i} (t) ξ_{j} (t^{'})] = δ_{i j} δ (t - t^{'}); i = j = 1,2,3

(17)

where δ_ij is the Kronecker symbol; δ is the δ-dirac function.

Let

\begin{array}{l} x_{1} (t) = u_{1} (t) + S^{*}; \\ x_{2} (t) = u_{2} (t) + P^{*}; \\ x_{3} (t) = u_{3} (t) + T^{*}; \end{array}

(18)

\frac{d x_{1}}{d t} = \frac{d u_{1} (t)}{d t}; \frac{d x_{2}}{d t} = \frac{d u_{2} (t)}{d t}; \frac{d x_{3}}{d t} = \frac{d u_{3} (t)}{d t}

(19)

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

\begin{array}{l} \frac{d u_{1} (t)}{d t} = a_{1} u_{1} (t) + a_{1} S^{*} - α_{11} {u_{1}}^{2} (t) - α_{11} {(S^{*})}^{2} \\ - 2 α_{11} u_{1} (t) S^{*} - α_{12} u_{1} (t) u_{2} (t) + α_{12} u_{1} (t) P^{*} \\ - α_{12} u_{2} (t) S^{*} - α_{12} S^{*} P^{*} - α_{13} u_{1} (t) u_{3} (t) \\ - α_{13} u_{1} (t) T^{*} - α_{13} u_{3} (t) S^{*} + α_{13} S^{*} T^{*} \\ - q_{1} E_{1} u_{1} (t) - q_{1} E_{1} S^{*} + η_{1} ξ_{1} (t) \end{array}

(20)

The linear part of (20) is

\begin{array}{l} \frac{d u_{1} (t)}{d t} = - α_{11} u_{1} (t) S^{*} - α_{12} u_{2} (t) S^{*} \\ - α_{13} u_{3} (t) S^{*} + η_{1} ξ_{1} (t) \end{array}

(21)

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

\begin{array}{l} \frac{d u_{2} (t)}{d t} = a_{2} u_{2} (t) + a_{2} P^{*} - α_{22} {u_{2}}^{2} (t) - α_{22} {(P^{*})}^{2} \\ - 2 α_{22} u_{2} (t) P^{*} + α_{21} u_{1} (t - τ) u_{2} (t - τ) P^{*})^{2} \\ + α_{21} u_{1} (t - τ) P^{*} + α_{21} u_{2} (t - τ) S^{*} + α_{21} S^{*} P^{*} \\ - α_{23} u_{2} (t) u_{3} (t) - α_{23} u_{2} (t) T^{*} - α_{23} u_{3} (t) P^{*} \\ + α_{23} P^{*} T^{*} - q_{2} E_{2} u_{2} (t) - q_{2} E_{2} P^{*} + η_{2} ξ_{2} (t) \end{array}

(22)

The linear part of Eq. (22) is

\begin{array}{l} \frac{d u_{2} (t)}{d t} = - α_{22} u_{2} (t) P^{*} + α_{21} u_{1} (t - τ) P^{*} \\ - α_{23} u_{3} (t) P^{*} + η_{2} ξ_{2} (t) \end{array}

(23)

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

\begin{array}{l} \frac{d u_{3} (t)}{d t} = a_{3} u_{3} (t) + a_{3} T^{*} - α_{33} {u_{3}}^{2} (t) - α_{33} {(T^{*})}^{2} \\ - 2 α_{33} u_{3} (t) T^{*} + α_{31} u_{1} (t - τ) u_{3} (t - τ) \\ + α_{31} u_{1} (t - τ) T^{*} + α_{31} u_{3} (t - τ) S^{*} + α_{31} S^{*} T^{*} \\ - α_{32} u_{2} (t) u_{3} (t) - α_{32} u_{2} (t) T^{*} - α_{32} u_{3} (t) P^{*} \\ + α_{32} P^{*} T^{*} - q_{3} E_{3} u_{3} (t) - q_{3} E_{3} T^{*} + η_{3} ξ_{3} (t) \end{array}

(24)

The linear part of (24) is

\begin{array}{l} \frac{d u_{3} (t)}{d t} = - α_{33} u_{3} (t) T^{*} + α_{31} u_{1} (t - τ) T^{*} \\ - α_{32} u_{2} (t) T^{*} + η_{3} ξ_{3} (t) \end{array}

(25)

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

\begin{array}{l} i ω {\tilde{u}}_{1} (ω) = - α_{11} S^{*} {\tilde{u}}_{1} (ω) - α_{12} S^{*} {\tilde{u}}_{2} (ω) \\ - α_{13} S^{*} {\tilde{u}}_{3} (ω) + η_{1} {\tilde{ξ}}_{1} (ω) \end{array}

(or)

\begin{array}{l} η_{1} {\tilde{ξ}}_{1} (ω) = (i ω + α_{11} S^{*}) {\tilde{u}}_{1} (ω) + α_{12} S^{*} {\tilde{u}}_{2} (ω) \\ + α_{13} S^{*} {\tilde{u}}_{3} (ω) \\ i ω {\tilde{u}}_{2} (ω) = - α_{22} P^{*} {\tilde{u}}_{2} (ω) + α_{21} P^{*} e^{- i ω τ} {\tilde{u}}_{1} (ω) \\ - α_{23} P^{*} {\tilde{u}}_{3} (ω) + η_{2} {\tilde{ξ}}_{2} (ω) \end{array}

(26)

(or)

\begin{array}{l} η_{2} {\tilde{ξ}}_{2} (ω) = - α_{21} P^{*} e^{- i ω τ} {\tilde{u}}_{1} (ω) + (i ω + α_{22} P^{*}) {\tilde{u}}_{2} (ω) \\ + α_{23} P^{*} {\tilde{u}}_{3} (ω) \\ i ω {\tilde{u}}_{3} (ω) = - α_{33} T^{*} {\tilde{u}}_{3} (ω) + α_{31} T^{*} e^{- i ω τ} {\tilde{u}}_{1} (ω) \\ - α_{32} T^{*} {\tilde{u}}_{2} (ω) + η_{3} {\tilde{ξ}}_{3} (ω) \end{array}

(27)

(or)

\begin{array}{l} η_{3} {\tilde{ξ}}_{3} (ω) = - α_{31} T^{*} e^{- i ω τ} {\tilde{u}}_{1} (ω) + α_{32} T^{*} {\tilde{u}}_{2} (ω) \\ + (i ω + α_{33} T^{*}) {\tilde{u}}_{3} (ω) \end{array}

(28)

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

M (ω) \tilde{u} (ω) = \tilde{ξ} (ω)

(29)

where

M (ω) = (\begin{array}{c} A (ω) & B (ω) & C (ω) \\ D (ω) & E (ω) & F (ω) \\ G (ω) & H (ω) & I (ω) \end{array});

\tilde{u} (ω) = [\begin{array}{c} {\tilde{u}}_{1} (ω) \\ {\tilde{u}}_{2} (ω) \\ {\tilde{u}}_{3} (ω) \end{array}]; \tilde{ξ} (ω) = [\begin{array}{c} {\tilde{ξ}}_{1} (ω) \\ {\tilde{ξ}}_{2} (ω) \\ {\tilde{ξ}}_{3} (ω) \end{array}]

\begin{array}{l} A (ω) = i ω + α_{11} S^{*}; B (ω) = α_{12} S^{*}; \\ C (ω) = α_{13} S^{*}; D (ω) = - α_{21} P^{*} e^{- i ω τ}; \\ E (ω) = i ω + α_{22} P^{*}; G (ω) = - α_{31} T^{*} e^{- i ω t}; \\ H (ω) = α_{32} T^{*}; I (ω) = i ω + α_{33} T^{*} \end{array}

(30)

Eq. (29) can also be written as

\tilde{u} (ω) = {[M (ω)]}^{- 1} \tilde{ξ} (ω)

Let [M(ω)]⁻¹ = K(ω), therefore

\tilde{u} (ω) = K (ω) \tilde{ξ} (ω)

(31)

where

K (ω) = \frac{A d j M (ω)}{|M (ω)|}

(32)

if the function Y(t) has a zero mean value, then the fluctuation intensity (variance) of its components in the frequency interval [ω, ω + dω] is S_Y(ω)dω where S_Y(ω) is spectral density of Y and is defined as

S_{Y} (ω) = \lim_{\tilde{T} \to \infty} \frac{{|\tilde{Y} (ω)|}^{2}}{\tilde{T}}

(33)

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

C_{Y} (τ) = \frac{1}{2 π} \int_{- \infty}^{\infty} S_{Y} (ω) e^{i ω τ} d ω

(34)

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

{σ_{Y}}^{2} = C_{Y} (0) = \frac{1}{2 π} \int_{- \infty}^{\infty} S_{Y} (ω) d ω

(35)

and the autocorrelation function is the normalized autocovariance

P_{Y} (τ) = \frac{C_{Y} (τ)}{C_{Y} (0)}

(36)

For a Gaussian white noise process, it is

\begin{array}{l} S_{ξ_{i} ξ_{j}} (ω) = \lim_{\hat{T} \to + \infty} \frac{E [{\tilde{ξ}}_{i} (ω) {\tilde{ξ}}_{j} (ω)]}{\hat{T}} \\ = \lim_{\hat{T} \to + \infty} \frac{1}{\hat{T}} \int_{- \frac{\hat{T}}{2}}^{\frac{\hat{T}}{2}} \int_{- \frac{\hat{T}}{2}}^{\frac{\hat{T}}{2}} E [{\tilde{ξ}}_{i} (t) {\tilde{ξ}}_{j} (t^{'})] e^{- i ω (t - t^{'})} d t d t^{'} \\ = δ_{i j} \end{array}

(37)

From (31), we have

{\tilde{u}}_{i} (ω) = \sum_{j = 1}^{3} K_{i j} (ω) {\tilde{ξ}}_{j} (ω); i = 1,2,3

(38)

From (33) we have

S_{u_{i}} (ω) = \sum_{j = 1}^{3} η_{j} {|K_{i j} (ω)|}^{2}; i = 1,2,3

(39)

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

{σ_{u_{i}}}^{2} = \frac{1}{2 π} \sum_{j = 1}^{3} \int_{- \infty}^{\infty} η_{j} {|K_{i j} (ω)|}^{2} d ω; i = 1,2,3

(40)

and by (31) and (32), we obtain three variances of u_i (i = 1, 2, 3) of the model system (13)–(15) as follows:

\begin{array}{l} {σ_{u_{1}}}^{2} = \\ \frac{1}{2 π} \int_{- \infty}^{\infty} [η_{1} {|\frac{A d j (1)}{|M (ω)|}|}^{2} + η_{2} {|\frac{A d j (2)}{|M (ω)|}|}^{2} + η_{3} {|\frac{A d j (3)}{|M (ω)|}|}^{2}] d ω \\ {σ_{u_{2}}}^{2} = \\ \frac{1}{2 π} \int_{- \infty}^{\infty} [η_{1} {|\frac{A d j (4)}{|M (ω)|}|}^{2} + η_{2} {|\frac{A d j (5)}{|M (ω)|}|}^{2} + η_{3} {|\frac{A d j (6)}{|M (ω)|}|}^{2}] d ω \\ {σ_{u_{3}}}^{2} = \\ \frac{1}{2 π} \int_{- \infty}^{\infty} [η_{1} {|\frac{A d j (7)}{|M (ω)|}|}^{2} + η_{2} {|\frac{A d j (8)}{|M (ω)|}|}^{2} + η_{3} {|\frac{A d j (9)}{|M (ω)|}|}^{2}] d ω \end{array}

(41)

where |M(ω)| = |R(ω)| + i|I(ω)|.

The real part:

\begin{array}{l} |M (ω)| = R^{2} (ω) = - α_{33} ω^{2} T^{*} - ω^{2} α_{22} P^{*} - ω^{2} α_{11} S^{*} \\ + α_{11} α_{22} α_{33} S^{*} P^{*} - α_{11} α_{23} α_{32} S^{*} T^{*} P^{*} \\ + α_{12} α_{21} α_{33} S^{*} P^{*} T^{*} \cos ω τ + ω α_{12} α_{21} S^{*} P^{*} \sin ω τ \\ - α_{12} α_{31} α_{23} S^{*} T^{*} P^{*} \cos ω τ - α_{13} α_{21} α_{32} S^{*} T^{*} P^{*} \cos ω τ \\ + ω α_{13} α_{31} S^{*} T^{*} \sin ω τ + α_{13} α_{31} α_{22} S^{*} P^{*} T^{*} \end{array}

(42)

The imaginary part:

\begin{array}{l} |M (ω)| = I^{2} (ω) = - ω^{3} + ω α_{22} α_{33} P^{*} T^{*} + ω α_{11} α_{33} S^{*} T^{*} \\ + ω α_{11} α_{22} S^{*} P^{*} - ω α_{23} α_{32} T^{*} P^{*} \\ + ω α_{12} α_{21} S^{*} P^{*} \cos ω τ - α_{12} α_{21} α_{33} S^{*} P^{*} T^{*} \sin ω τ \\ + α_{12} α_{31} α_{23} S^{*} T^{*} P^{*} \sin ω τ + α_{13} α_{21} α_{32} S^{*} T^{*} P^{*} \sin ω τ \\ + ω α_{13} α_{31} S^{*} T^{*} \cos ω τ - α_{22} α_{13} α_{31} S^{*} P^{*} T^{*} \sin ω τ \end{array}

(43)

and

\begin{array}{l} {|A d j (1)|}^{2} = {X_{1}}^{2} + {Y_{1}}^{2}; {|A d j (2)|}^{2} = {X_{2}}^{2} + {Y_{2}}^{2}; \\ {|A d j (3)|}^{2} = {X_{3}}^{2} + {Y_{3}}^{2}; {|A d j (4)|}^{2} = {X_{4}}^{2} + {Y_{4}}^{2}; \\ {|A d j (5)|}^{2} = {X_{5}}^{2} + {Y_{5}}^{2}; {|A d j (6)|}^{2} = {X_{6}}^{2} + {Y_{6}}^{2}; \\ {|A d j (7)|}^{2} = {X_{7}}^{2} + {Y_{7}}^{2}; {|A d j (8)|}^{2} = {X_{8}}^{2} + {Y_{8}}^{2}; \\ {|A d j (9)|}^{2} = {X_{9}}^{2} + {Y_{9}}^{2} \end{array}

where

X_{1} = - ω^{2} + α_{22} α_{23} P^{*} T^{*} - α_{23} α_{32} T^{*} P^{*}

Y_{1} = ω α_{33} T^{*} + ω α_{22} P^{*}

\begin{array}{l} X_{2} = α_{21} α_{33} P^{*} T^{*} \cos ω τ - α_{23} α_{31} T^{*} P^{*} \cos ω τ \\ + ω α_{21} P^{*} \sin ω τ \end{array}

\begin{array}{l} Y_{2} = ω α_{21} P^{*} \cos ω τ - ω α_{21} α_{33} P^{*} T^{*} \sin ω τ \\ + α_{23} α_{31} T^{*} P^{*} \end{array}

\begin{array}{l} X_{3} = - α_{21} α_{32} P^{*} T^{*} \cos ω τ + α_{31} α_{22} T^{*} P^{*} \cos ω τ \\ + ω α_{31} T^{*} \sin ω τ \end{array}

\begin{array}{l} Y_{3} = ω α_{31} T^{*} \cos ω τ + α_{21} α_{32} P^{*} T^{*} \sin ω τ \\ - α_{31} α_{22} T^{*} P^{*} \sin ω τ \end{array}

X_{4} = α_{12} α_{13} S^{*} T^{*} - α_{32} α_{13} S^{*} T^{*}; Y_{4} = ω α_{12} S^{*}

X_{5} = - ω^{2} + α_{11} α_{33} S^{*} T^{*} + α_{31} α_{13} S^{*} T^{*} \cos ω τ

Y_{5} = ω α_{33} T^{*} + ω α_{11} S^{*} - α_{31} α_{13} S^{*} T^{*} \sin ω τ

X_{6} = α_{11} α_{32} S^{*} T^{*} + α_{31} α_{12} S^{*} T^{*} \cos ω τ

Y_{6} = ω α_{32} T^{*} - α_{31} α_{12} S^{*} T^{*} \sin ω τ

X_{7} = α_{12} α_{23} S^{*} P^{*} - α_{22} α_{13} S^{*} P^{*}; Y_{7} = - ω α_{13} S^{*}

X_{8} = α_{12} α_{23} S^{*} P^{*} + α_{21} α_{13} S^{*} P^{*} \cos ω τ

Y_{8} = ω α_{23} P^{*} - α_{21} α_{13} S^{*} P^{*}

X_{9} = ω^{2} + α_{11} α_{22} S^{*} P^{*} + α_{12} α_{21} S^{*} P^{*} \cos ω τ

Y_{9} = ω α_{22} P^{*} + ω α_{11} S^{*} - α_{12} α_{21} S^{*} P^{*} \sin ω τ

Thus Eq. (41) becomes

{σ_{u_{1}}}^{2} = \frac{1}{2 π} \{\int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} [\begin{array}{l} η_{1} ({X_{1}}^{2} + {Y_{1}}^{2}) \\ + η_{2} ({X_{4}}^{2} + {Y_{4}}^{2}) \\ + η_{3} ({X_{7}}^{2} + {Y_{7}}^{2}) \end{array}] d ω\}

{σ_{u_{1}}}^{2} = \frac{1}{2 π} \{\int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} [\begin{array}{l} η_{1} ({X_{2}}^{2} + {Y_{2}}^{2}) + \\ η_{2} ({X_{5}}^{2} + {Y_{5}}^{2}) + \\ η_{3} ({X_{8}}^{2} + {Y_{8}}^{2}) \end{array}] d ω\}

{σ_{u_{1}}}^{2} = \frac{1}{2 π} \{\int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} [\begin{array}{l} η_{1} ({X_{3}}^{2} + {Y_{3}}^{2}) \\ + η_{2} ({X_{6}}^{2} + {Y_{6}}^{2}) \\ + η_{3} ({X_{9}}^{2} + {Y_{9}}^{2}) \end{array}] d ω\}

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

(i) if η₁ = 0, η₂ = 0, then

{σ_{u_{1}}}^{2} = \frac{η_{3} ({X_{7}}^{2} + {Y_{7}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

{σ_{u_{2}}}^{2} = \frac{η_{3} ({X_{8}}^{2} + {Y_{8}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

{σ_{u_{3}}}^{2} = \frac{η_{3} ({X_{9}}^{2} + {Y_{9}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

(ii) if η₂ = 0, η₃ = 0, then

{σ_{u_{1}}}^{2} = \frac{η_{1} ({X_{1}}^{2} + {Y_{1}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

{σ_{u_{2}}}^{2} = \frac{η_{1} ({X_{2}}^{2} + {Y_{2}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

{σ_{u_{3}}}^{2} = \frac{η_{1} ({X_{3}}^{2} + {Y_{3}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

(iii) if η₃ = 0, η₁ = 0, then

{σ_{u_{1}}}^{2} = \frac{η_{2} ({X_{4}}^{2} + {Y_{4}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

{σ_{u_{2}}}^{2} = \frac{η_{2} ({X_{5}}^{2} + {Y_{5}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} d ω

{σ_{u_{3}}}^{2} = \frac{η_{2} ({X_{6}}^{2} + {Y_{6}}^{2})}{2 π} \int_{- \infty}^{\infty} \frac{1}{R^{2} (ω) + I^{2} (ω)} 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:

\begin{array}{l} x_{1} = 10; x_{2} = 10; x_{3} = 20; \\ a_{1} = 8; α_{11} = 0.05; α_{12} = 0.6; α_{13} = 0.7; q_{1} = 0.15; \\ E_{1} = 10; a_{2} = 2; α_{22} = 0.7; \\ α_{21} = 0.17; α_{23} = 0.13; q_{2} = 0.01; E_{2} = 15; \\ a_{3} = 1; α_{33} = 0.4; α_{31} = 0.15; α_{32} = 0.18; \\ q_{3} = 0.03; E_{3} = 5; \end{array}

Fig. 1 shows that the variation of population against time, initially with x₁ = 10; x₂ = 10; x₃ = 20; in three different graphs.

7.2 Example 2

We take the following parameters:

\begin{array}{l} x_{1} = 10; x_{2} = 10; x_{3} = 20; a_{1} = 8; α_{11} = 0.05; \\ α_{12} = 0.6; α_{13} = 0.7; q_{1} = 0.15; E_{1} = 10; \\ a_{2} = 2; α_{22} = 0.7; α_{21} = 0.17; α_{23} = 0.13; \\ q_{2} = 0.01; E_{2} = 15; a_{3} = 1; α_{33} = 0.4; \\ α_{31} = 0.15; α_{32} = 0.18; q_{3} = 0.03; E_{3} = 5; \end{array}

Fig. 2 shows that the variation of population against time, initially with x₁ = 10; x₂ = 10; x₃ = 20; in a single graph.

7.3 Example 3

We take the following parameters:

\begin{array}{l} x_{1} = 10; x_{2} = 15; x_{3} = 20; a_{1} = 8; α_{11} = 0.05; \\ α_{12} = 0.6; α_{13} = 0.7; q_{1} = 0.15; E_{1} = 10; \\ a_{2} = 2; α_{22} = 0.7; α_{21} = 0.17; α_{23} = 0.13; \\ q_{2} = 0.01; E_{2} = 15; a_{3} = 1; α_{33} = 0.4; \\ α_{31} = 0.15; α_{32} = 0.18; \\ q_{3} = 0.03; E_{3} = 5; \end{array}

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

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 x₁ = 10; x₂ = 15; x₃ = 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 x₁ = 10; x₂ = 15; x₃ = 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

\begin{array}{l} V (x_{1}, x_{2}, x_{3}) \\ = (x_{1} - {x_{1}}^{*} - {x_{1}}^{*} \ln (\frac{x_{1}}{{x_{1}}^{*}})) \\ + l_{1} (x_{2} - {x_{2}}^{*} - {x_{2}}^{*} \ln (\frac{x_{2}}{{x_{2}}^{*}})) \\ + l_{2} (x_{3} - {x_{3}}^{*} - {x_{3}}^{*} \ln (\frac{x_{3}}{{x_{3}}^{*}})) \end{array}

where l₁ and l₂ 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 G₄(x₁*, x₂*, x₃*). The time derivative of V along the trajectories of (1)–(3) is

\begin{array}{l} \frac{d V}{d t} = (\frac{x_{1} - {x_{1}}^{*}}{x_{1}}) \frac{d x_{1}}{d t} + l_{1} (\frac{x_{2} - {x_{2}}^{*}}{x_{2}}) \frac{d x_{2}}{d t} \\ + l_{2} (\frac{x_{3} - {x_{3}}^{*}}{x_{3}}) \frac{d x_{3}}{d t} \end{array}

Choosing $l_{1} = \frac{α_{12}}{α_{21}}$ , $l_{2} = \frac{α_{13}}{α_{31}}$ after a little algebraic manipulation, we get,

\begin{array}{l} \frac{d V}{d t} = - α_{11} {(x_{1} - {x_{1}}^{*})}^{2} \\ - \frac{α_{12} α_{22}}{α_{21}} {(x_{2} - {x_{2}}^{*})}^{2} - \frac{α_{33} α_{13}}{α_{31}} {(x_{3} - {x_{3}}^{*})}^{2} \\ - \frac{α_{23} α_{12}}{α_{21}} (x_{2} - {x_{2}}^{*}) (x_{3} - {x_{3}}^{*}) - \frac{α_{13} α_{32}}{α_{31}} (x_{2} - {x_{2}}^{*}) (x_{3} - {x_{3}}^{*}) \end{array}

\begin{array}{l} \frac{d V}{d t} < - α_{11} {(x_{1} - {x_{1}}^{*})}^{2} - \frac{α_{12} α_{22}}{α_{21}} {(x_{2} - {x_{2}}^{*})}^{2} \\ - \frac{α_{33} α_{13}}{α_{31}} {(x_{3} - {x_{3}}^{*})}^{2} \\ - (\frac{α_{23} α_{12}}{α_{21}} + \frac{α_{13} α_{32}}{α_{31}}) [\frac{{(x_{2} - {x_{2}}^{*})}^{2}}{2} + \frac{{(x_{3} - {x_{3}}^{*})}^{2}}{2}] \end{array}

Now it is very clear that $\frac{d V}{d t} < 0$ .

Since $\frac{d V}{d t} < 0$ in some neighbourhood (x₁*, x₂*, x₃*) therefore the interior equilibrium point (x₁*, x₂*, x₃*) is globally asymptotically stable.

Bibliographie

[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

Commentaires - Politique

Ces articles pourraient vous intéresser

Combined harvesting of a stage structured prey–predator model incorporating cannibalism in competitive environment

Kunal Chakraborty; Kunal Das; Tapan Kumar Kar

C. R. Biol (2013)

A focus on long-run sustainability of a harvested prey predator system in the presence of alternative prey

T.K. Kar; S.K. Chattopadhyay

C. R. Biol (2010)

Stability and bifurcation of a prey–predator model with time delay

T.K. Kar; Ashim Batabyal

C. R. Biol (2009)