Mathematical modeling applied to renewable fishery management

Renewable resources are under extreme pressure worldwide in spite of taking efforts to design regulation to control the excessive and unsustainable exploitation. Fish is a renewable but finite resource. With the rapid growth of industrialization and population, the exploitation of fisheries has increased significantly. Although exploitation of resources is essential for the growth and development of a country, unplanned exploitation eventually leads to the extinction of the resources. Consequently, this will affect the growth and survival of species depending on the resource. From this point of view, it is a major challenge to manage renewable resources for the sustainable development of the country.

Mathematical modeling can play a significant role in the efficient and sustainable management of renewable resources. It is mainly used to describe the real phenomena leading to design better prediction, prevention, management and control techniques. Several well documented mathematical models regarding real life problems can be found in [1-7].

During the last few decades, mathematical models regarding renewable fishery resource management have been described. Mathematical modeling in harvesting of fisheries was studied first by Clark [5]. Biswas et al. presented a model for fishery resource with reserve area [1]. They studied the dynamics of a fishery resource in a two patch environment: a free fishing zone and a reserved zone where fishing is not allowed. Dubey et al. [6] presented a model for fishery resource with reserve area. They studied the dynamics of a fishery resource in a two patch environment: a free fishing zone and a reserved zone where fishing is not allowed. An optimal harvesting strategy is also derived using Pontryagin’s Maximum Principle.  Chaudhury [9] analyzed the dynamic optimization of combined harvesting of a two species fishery. Kar [10] presented a model for fishery resource with reserved area and facing prey-predator interaction. He considered that predation takes place in the unreserved zone. Local and global stability, optimal harvesting policy are also discussed. Roy et al. [11] investigated the effects of two predators on a prey population. They considered different types of functional responses to formulate the model.

From the literature discussed above and to the best of our views, a model for the cultivation of black tiger prawn was proposed and the work of Biswas et al. [1] was extended in this study. The dynamics of single species fishery in reserved and unreserved zone were discussed in the work but the effect of predation that the species might face in the unreserved zone were not studied. Also, it was considered that fish population can migrate from unreserved to reserved zone and vice versa. The effect of predation on the fish production was aimed to study and the migration from reserved to unreserved area was also restricted in this present work.  In this paper, a mathematical model of a prey predator fishery was proposed with the help of system of nonlinear differential equation. It was considered that the fish species in the unreserved area were related in prey predator relationship. Holling type II functional response was taken into account to study the interaction between prey and predator fish species. Harvesting was permissible only in the unreserved region. Predation and harvesting were restricted in the reserved regions. To analyze the model, the existence of equilibrium points, dynamical behavior of the points and also the stability and instability conditions were discussed. Finally, numerical simulations were carried out to verify the analytical result of our proposed model.

The model (1)-(3) had been analyzed in order to describe the dynamics of the fish species. For the analysis of the model the following studies were considered:

3.1 Boundedness of the model

To prove that the model system is biologically well posed the following Lemma was to be satisfied.

Lemma 1: The set

$\Omega =\left\{ \left( {{x}_{1}},{{x}_{2}},{{x}_{3}} \right):w\left( t \right)={{x}_{1}}\left( t \right)+{{x}_{2}}\left( t \right)+{{x}_{3}}\left( t \right),0<w\left( t \right)<\frac{\delta }{\eta } \right\}$

attracts all solutions initiating in the interior of the positive orthant, where η is a constant and

$\mu \text{=}\frac{{{k}_{1}}}{4{{r}_{1}}}{{({{r}_{1}}+\eta -\mu -qE)}^{2}}+\frac{{{k}_{2}}}{4{{r}_{2}}}{{({{r}_{2}}+\eta -\beta -\alpha )}^{2}}+\frac{{{k}_{3}}}{4{{r}_{3}}}{{({{r}_{3}}+\eta )}^{2}}$

Proof: Let, $w\left( t \right)={{x}_{1}}\left( t \right)+{{x}_{2}}\left( t \right)+{{x}_{3}}\left( t \right),\text{ }\eta >0$ be a constant. Then we can write,

$\begin{align}  & \frac{dw}{dt}+\eta w=\left( {{r}_{1}}+\eta -\mu -qE \right){{x}_{1}}-\frac{{{r}_{1}}}{{{k}_{1}}}{{x}_{1}}^{2}+\left( {{r}_{2}}+\eta -\alpha -\beta  \right){{x}_{2}}-\frac{{{r}_{2}}}{{{k}_{2}}}{{x}_{2}}^{2} \\  & -\left( m-n \right)\frac{{{x}_{1}}{{x}_{3}}}{a+{{x}_{1}}}+\left( {{r}_{3}}+\eta -{{\beta }_{1}} \right){{x}_{3}}-\frac{{{r}_{3}}}{{{k}_{3}}}{{x}_{3}}^{2} \\ \end{align}$

Since m is the depletion rate coefficient of prey due to its intake by the predator and n is the growth rate coefficient of predator due to its interaction with their prey, so it is assumed that m≥n. Now η is chosen such that 0<η<β₁.

$\begin{align}  & \frac{dw}{dt}+\eta w\le \frac{{{k}_{1}}}{4{{r}_{1}}}{{\left( {{r}_{1}}+\eta -\mu -qE \right)}^{2}}+\frac{{{k}_{2}}}{4{{r}_{2}}}{{\left( {{r}_{2}}+\eta -\alpha -\beta  \right)}^{2}} \\ & +\frac{{{k}_{3}}}{4{{r}_{3}}}{{\left( {{r}_{3}}+\eta  \right)}^{2}} \\\end{align}$

By using differential inequality, we get,

$0<w({{x}_{1}}(t),{{x}_{2}}(t),{{x}_{3}}(t))\le \frac{\delta }{\eta }(1-{{e}^{-\eta t}})+({{x}_{1}}(0),{{x}_{2}}(0),{{x}_{3}}(0)){{e}^{-\eta t}}$

Taking limit as t→∞, we get, 0<w(t)≤δ/η.

3.2 Positivity of the solution of the model

Lemma 2:  For${{x}_{1}}\left( 0 \right)\ge 0,{{x}_{2}}\left( 0 \right)\ge 0,{{x}_{3}}\left( 0 \right)\ge 0$, the solutions ${{x}_{1}}\left( t \right),{{x}_{2}}\left( t \right),{{x}_{3}}\left( t \right)$ are all non-negative for all t≥0

Proof: For positivity, equation (2.1) can be written as, $\frac{d{{x}_{1}}}{dt}\ge -\left( \mu +\sigma +qE \right){{x}_{1}}$ $\Rightarrow \frac{d{{x}_{1}}}{{{x}_{1}}}\ge -kdt$ where,$k=\left( \mu +\sigma +qE \right)$ $\Rightarrow {{x}_{1}}\ge {{c}_{1}}{{e}^{-kt}}$, where c₁ is an integrating constant. Applying the initial condition, at $t=0,{{x}_{1}}\left( 0 \right)\ge 0$, we get, ${{x}_{1}}\left( 0 \right)={{c}_{1}}$. Putting the value of c₁ in the equation, we get, ${{x}_{1}}\left( t \right)\ge {{x}_{1}}\left( 0 \right){{e}^{-kt}}$. When $t\to \infty ,{{x}_{1}}\left( t \right)\ge 0$. Therefore x₁(t) is positive for all t≥0. Again equation (2.2) can be written as $\frac{d{{x}_{2}}}{dt}\ge -\left( \alpha +\beta  \right){{x}_{2}}$ $\Rightarrow \frac{d{{x}_{2}}}{{{x}_{2}}}\ge -Ldt$, where $L=\left( \alpha +\beta  \right)$. Integrating we get, ${{x}_{2}}\left( t \right)\ge {{c}_{2}}{{e}^{-Lt}}$, where c₂ is an integrating constant. Applying the initial condition, at $t=0,{{x}_{2}}\left( 0 \right)\ge 0$, we get, ${{x}_{2}}\left( 0 \right)={{c}_{2}}$. Putting the value of c₂, we get, ${{x}_{2}}\left( t \right)\ge {{x}_{2}}\left( 0 \right){{e}^{-Lt}}$. When $t\to \infty ,\text{ }{{x}_{2}}\left( t \right)\ge 0$. Hence x₂(t) is positive for all t≥0. To prove that x₃(t) is positive, equation (2.3) is written as, $\frac{d{{x}_{3}}}{dt}\ge -{{\beta }_{1}}{{x}_{3}}$ $\Rightarrow \frac{d{{x}_{3}}}{{{x}_{3}}}\ge -{{\beta }_{1}}dt$. Integrating ${{x}_{3}}(t)\ge {{c}_{3}}{{e}^{-{{\beta }_{1}}t}}$, where c₃ is an integrating constant. Now, at $t=0,{{x}_{3}}\left( 0 \right)\ge 0.$

So, ${{x}_{3}}\left( 0 \right)={{c}_{3}}.$

Putting the value of c₃, we obtain, ${{x}_{3}}\left( t \right)\ge {{x}_{3}}\left( 0 \right){{e}^{-{{\beta }_{1}}t}}$. At $t\to \infty ,\text{ }{{x}_{3}}\left( t \right)\ge 0$. Therefore, x₃(t) is positive for all t≥0. Hence this completes the proof.

3.3 Existence of equilibria

Let $\left( {{{\bar{x}}}_{1}},{{{\bar{x}}}_{2}},0 \right)$ be the positive solution of the equations (1) and (2)

${{\varphi }_{1}}{{x}_{1}}-\frac{{{r}_{1}}{{x}_{1}}^{2}}{{{k}_{1}}}=0$   (5)

$\frac{d{{x}_{2}}}{dt}={{\varphi }_{2}}{{x}_{2}}-\frac{{{r}_{2}}{{x}_{2}}^{2}}{{{k}_{2}}}+\sigma {{x}_{1}}\text{=0         }$  (6)   

where, ${{\varphi }_{1}}={{r}_{1}}-\mu -\sigma -qE$and ${{\varphi }_{2}}={{r}_{2}}-\alpha -\beta $. From (5), we get, ${{x}_{1}}=0,{{\bar{x}}_{1}}=\frac{{{\varphi }_{1}}{{k}_{1}}}{{{r}_{1}}}$. Putting the value of ${{\bar{x}}_{1}}$in (6), we get,

${{\bar{x}}_{2}}=\frac{{{\varphi }_{2}}\pm \sqrt{\varphi _{2}^{2}+\frac{4\sigma {{r}_{2}}{{k}_{1}}{{\varphi }_{1}}}{{{r}_{1}}{{k}_{2}}}}}{\frac{2{{r}_{2}}}{{{k}_{2}}}}$

Again, let ($x_{1}^{*},x_{2}^{*},x_{3}^{*}$ ) be the positive solution of the equations

${{\varphi }_{1}}{{x}_{1}}-\frac{{{r}_{1}}{{x}_{1}}^{2}}{{{k}_{1}}}-\frac{m{{x}_{1}}{{x}_{3}}}{a+{{x}_{1}}}=0$  (7)

${{\varphi }_{2}}{{x}_{2}}-\frac{{{r}_{2}}{{x}_{2}}^{2}}{{{k}_{2}}}+\sigma {{x}_{1}}\text{=0         }$ (8)   

${{\varphi }_{3}}{{x}_{3}}+\frac{n{{x}_{1}}{{x}_{3}}}{a+{{x}_{1}}}-\frac{{{r}_{3}}{{x}_{3}}^{2}}{{{k}_{3}}}=0$  (9) 

where, ${{\varphi }_{3}}={{r}_{3}}-{{\beta }_{1}}$. Equation (9) yields

${{x}_{3}}=0,\text{ }a{{\varphi }_{3}}+{{\varphi }_{3}}{{x}_{1}}-\frac{a{{r}_{3}}{{x}_{3}}}{{{k}_{3}}}-\frac{{{r}_{3}}{{x}_{1}}{{x}_{3}}}{{{k}_{3}}}+n{{x}_{1}}=0$ $\Rightarrow {{x}_{1}}^{*}=\frac{a{{r}_{3}}{{x}_{3}}^{*}-a{{\varphi }_{3}}{{k}_{3}}}{{{\varphi }_{3}}{{k}_{3}}+n{{k}_{3}}-{{r}_{3}}{{x}_{3}}^{*}}$.

Using the value of ${{x}_{1}}^{*}$ in equation (8), we get

$\frac{{{r}_{2}}x_{2}^{{{*}^{2}}}}{{{k}_{2}}}-{{\varphi }_{2}}x_{2}^{*}-\sigma (\frac{a{{r}_{3}}x_{3}^{*}-a{{\varphi }_{3}}{{k}_{3}}}{{{\varphi }_{3}}{{k}_{3}}+n{{k}_{3}}-{{r}_{3}}x_{3}^{*}})=0$

This equation will have a positive solution if ${{\varphi }_{2}}\ge 0\Rightarrow {{r}_{2}}-\alpha -\beta \ge 0$. Finally from equation (3.3), we obtain

${{x}_{3}}^{*}=\frac{a+{{x}_{1}}^{*}}{m}\left( {{\varphi }_{1}}-\frac{{{r}_{1}}{{x}_{1}}^{*}}{{{k}_{1}}} \right)$

Hence the equilibrium points of the system are ${{P}_{1}}(0,0,0),{{P}_{2}}({{\bar{x}}_{1}},{{\bar{x}}_{2}},0)\text{ and  }{{P}_{3}}\left( {{x}_{1}}^{*},{{x}_{2}}^{*},{{x}_{3}}^{*} \right)$.

3.3 Stability analysis at steady state

The Jacobean of the system is                                                                                                                 

$J\left( {{x}_{1}},{{x}_{2}},{{x}_{3}} \right)=\left( \begin{align}  & {{\text{r}}_{1}}\text{-}\frac{2{{r}_{1}}{{x}_{1}}}{{{k}_{1}}}\text{-}\left( \mu +\sigma +qE \right)\text{-}\frac{am{{x}_{3}}}{{{\left( a+{{x}_{1}} \right)}^{2}}}\text{            0                   -}\frac{m{{x}_{1}}}{a+{{x}_{1}}} \\ & \text{ }\sigma \text{                              }{{\text{r}}_{2}}\text{-}\frac{2{{r}_{2}}{{x}_{2}}}{{{k}_{2}}}\text{-}\beta -\alpha \text{                             0} \\ & \text{ }\frac{an{{x}_{3}}}{{{\left( a+{{x}_{1}} \right)}^{2}}}\text{                      0                                     }{{\text{r}}_{3}}-\frac{2{{r}_{3}}{{x}_{3}}}{{{k}_{3}}}\text{-}{{\beta }_{1}} \\ \end{align} \right)$ $\therefore \text{ }J\left( 0,0,0 \right)=\left( \begin{align}  & {{\text{r}}_{1}}-\left( \mu +\sigma +qE \right)\text{           0                       0} \\ & \text{ }\sigma \text{                              }{{\text{r}}_{2}}\text{-}\beta -\alpha \text{               0} \\ & \text{   0                                0                     }{{\text{r}}_{3}}-{{\beta }_{1}} \\\end{align} \right)$

Then the characteristics equation of the matrix with eigenvalue λ is

$\left| J-\lambda I \right|=\left| \begin{align}  & {{\text{r}}_{1}}-\left( \mu +\sigma +qE \right)\text{-}\lambda \text{           0                         0} \\ & \text{ }\sigma \text{                              }{{\text{r}}_{2}}\text{-}\beta -\alpha -\lambda \text{               0} \\ & \text{   0                                0                     }{{\text{r}}_{3}}-{{\beta }_{1}}-\lambda  \\\end{align} \right|=0$

$\Rightarrow {{\text{r}}_{1}}-\left( \mu +\sigma +qE \right)\text{-}\lambda =0,\text{ }{{\text{r}}_{2}}\text{-}\beta -\alpha -\lambda =0\text{ or }{{\text{r}}_{3}}-{{\beta }_{1}}-\lambda =0$

$\Rightarrow {{\lambda }_{1}}={{\text{r}}_{1}}-\left( \mu +\sigma +qE \right),\text{ }{{\lambda }_{2}}={{\text{r}}_{2}}\text{-}\beta -\alpha \text{  and }{{\lambda }_{3}}={{\text{r}}_{3}}-{{\beta }_{1}}$

The eigenvalue λ of the matrix determines the stability of the states. Depending on λ, the stability conditions are: 1. if the eigenvalue λ>0, then the steady state is unstable, 2. if the eigenvalue λ<0, then the system is stable.

In the following lemma, we show that ${{P}_{3}}(x_{1}^{*},x_{2}^{*},x_{3}^{*})$ is locally asymptotically stable.

Lemma 3. The equilibrium point ${{P}_{3}}(x_{1}^{*},x_{2}^{*},x_{3}^{*})$ is always locally asymptotically stable.

Proof: Let, the three functions be $u=u\left( {{x}_{1}},{{x}_{2}},{{x}_{3}} \right),v=v\left( {{x}_{1}},{{x}_{2}},{{x}_{3}} \right)\text{ and }w=w\left( {{x}_{1}},{{x}_{2}},{{x}_{3}} \right)$ of the system (2.1)-(2.3). Then the Jacobean at ${{P}_{3}}(x_{1}^{*},x_{2}^{*},x_{3}^{*})$ is

${{J}_{3}}(x_{1}^{*},x_{2}^{*},x_{3}^{*})=\left( \begin{align}  & {{\text{r}}_{1}}\text{-}\frac{2{{r}_{1}}{{x}_{1}}^{*}}{{{k}_{1}}}\text{-}\left( \mu +\sigma +qE \right)\text{-}\frac{am{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}}\text{      0                       }\frac{m{{x}_{1}}^{*}}{a+{{x}_{1}}^{*}} \\ & \text{     }\sigma \text{                                        }{{\text{r}}_{2}}\text{-}\frac{2{{r}_{2}}{{x}_{2}}^{*}}{{{k}_{2}}}\text{-}\beta -\alpha \text{               0} \\ & \text{ }\frac{an{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}}\text{                                        0                }{{\text{r}}_{3}}-\frac{2{{r}_{3}}{{x}_{3}}^{*}}{{{k}_{3}}}\text{-}{{\beta }_{1}} \\ \end{align} \right)$

Then the characteristics equation of the matrix with eigenvalue is λ.

|J-λI|=0

$\left| \begin{align}  & {{\text{r}}_{1}}\text{-}\frac{2{{r}_{1}}{{x}_{1}}^{*}}{{{k}_{1}}}\text{-}\left( \mu +\sigma +qE \right)\text{-}\frac{am{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}}\text{-}\lambda \text{     0                         -}\frac{m{{x}_{1}}^{*}}{a+{{x}_{1}}^{*}} \\ & \text{ }\sigma \text{                                              }{{\text{r}}_{2}}\text{-}\frac{2{{r}_{2}}{{x}_{2}}^{*}}{{{k}_{2}}}\text{-}\beta -\alpha -\lambda\text{            0} \\ & \text{ }\frac{an{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}}\text{                                         0                    }{{\text{r}}_{3}}-\frac{2{{r}_{3}}{{x}_{3}}^{*}}{{{k}_{3}}}\text{-}{{\beta }_{1}}-\lambda  \\\end{align} \right|=0$

$\begin{align}  & \Rightarrow \left( {{\text{r}}_{1}}\text{-}\frac{2{{r}_{1}}{{x}_{1}}^{*}}{{{k}_{1}}}\text{-}\left( \mu +\sigma +qE \right)\text{-}\frac{am{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}}\text{-}\lambda  \right)\left( {{\text{r}}_{2}}\text{-}\frac{2{{r}_{2}}{{x}_{2}}^{*}}{{{k}_{2}}}\text{-}\beta -\alpha -\lambda  \right) \\ & \left( {{\text{r}}_{3}}-\frac{2{{r}_{3}}{{x}_{3}}^{*}}{{{k}_{3}}}\text{-}{{\beta }_{1}}-\lambda  \right)+\left( \frac{m{{x}_{1}}^{*}}{a+{{x}_{1}}^{*}} \right)\left( \frac{an{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}} \right)\left( {{\text{r}}_{2}}\text{-}\frac{2{{r}_{2}}{{x}_{2}}^{*}}{{{k}_{2}}}\text{-}\beta -\alpha -\lambda \right) \\\end{align}=0$$\Rightarrow \left( {{A}_{1}}\text{-}\lambda  \right)\left( {{A}_{2}}-\lambda  \right)\left( {{A}_{3}}-\lambda  \right)+{{A}_{4}}\left( {{A}_{2}}-\lambda  \right)=0$

$\Rightarrow {{\lambda }^{3}}+{{a}_{1}}{{\lambda }^{2}}+{{a}_{2}}\lambda +{{a}_{3}}=0$ (10)

where,

${{a}_{1}}=-\left( {{A}_{1}}+{{A}_{2}}-{{A}_{3}} \right),{{a}_{2}}=\left( {{A}_{1}}{{A}_{2}}+{{A}_{2}}{{A}_{3}}+{{A}_{3}}{{A}_{1}}+{{A}_{4}} \right),{{a}_{3}}=-\left( {{A}_{1}}{{A}_{2}}{{A}_{3}}+{{A}_{2}}{{A}_{4}} \right)$,

$\begin{align}  & {{A}_{1}}={{\text{r}}_{1}}\text{-}\frac{2{{r}_{1}}{{x}_{1}}^{*}}{{{k}_{1}}}\text{-}\left( \mu +\sigma +qE \right)\text{-}\frac{am{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}},\text{ }{{A}_{2}}={{\text{r}}_{2}}\text{-}\frac{2{{r}_{2}}{{x}_{2}}^{*}}{{{k}_{2}}}\text{-}\beta -\alpha ,\text{ } \\ & {{A}_{3}}={{\text{r}}_{3}}-\frac{2{{r}_{3}}{{x}_{3}}^{*}}{{{k}_{3}}}\text{-}{{\beta }_{1}},\text{ }{{A}_{4}}=\left( \frac{m{{x}_{1}}^{*}}{a+{{x}_{1}}^{*}} \right)\left( \frac{an{{x}_{3}}^{*}}{{{\left( a+{{x}_{1}}^{*} \right)}^{2}}} \right) \\\end{align}$

By Routh-Hurwitch criterion, all the eigenvalues of (3.7) have negative real roots if and only if a₁>0, a₃>0 and a₁a₂>a₃.

Then the equilibrium point ${{P}_{3}}(x_{1}^{*},x_{2}^{*},x_{3}^{*})$ is locally asymptotically stable.

Considering two ecosystems, a mathematical model of fishery has been formulated in this paper. We have analyzed the behavior of the model, focusing on the parameters that are mainly responsible for the production and reduction of black tiger prawn. The numerical results reveal that the high mortality rate, fishing rate and predation rate of black tiger prawn in the unreserved area have a decreasing effect on the species in the reserved area. It also shows that the number of black tiger prawn increases in the reserved area due to the high migration from the unreserved area. When predation increases, the number of black tiger prawn in both reserved and unreserved area reduces by a significant amount. It is also evident from the simulations that the production of black tiger prawn in the reserved area decreases due to death from different diseases and for being stolen in absence of security. The main conclusion based on the result is that it is possible to maximize the production of black tiger prawn by proper management and thus it can play a significant role in the economy of a country.