Effects of Over-Harvesting and Drought on a Predator-Prey System with Optimal Control

Show more

1. Introduction

Predator-prey system is a dynamical system that explains the dynamics of ecological systems in which two species prey and predator do interact. The interactions can be through predation, competition, parasitism, mutualism and detritivory [1]. However, predation has been defined as the process in which species of one type consume species of another type, with condition that prey should be alive when it is captured [1][2]. Furthermore, [1]defined competition, parasitism, mutualism and detitrivory as follows: Competition is an interaction in which one organism consumes a resource that would have been available to and might have been consumed by another; parasitism is an interaction in which predators consume part of their prey, for example parasites and their hosts; mutualism is an interaction in which both species in the system benefits from each other; and detitrivory (commensalism) is an interaction in which a predator consumes prey which is dead already. In this study we consider prey and predator species that interact through predation.

The dynamic relationship between predators and their prey has been and continue to be one of the most important topic of discussion in ecology and mathematical ecology due to its universal existence and importance [3]. Therefore, the study of predator-prey dynamics helps scientists to understand interacting populations in the natural environment, which in turn guide them to develop solutions to problems affecting the growth rate of prey and predator species in the ecosystem [4]. The dynamics of predator-prey interactions is affected by factors such as over-harvesting, over predation, environmental pollution, mismanagement of the habitat, and hazards fire and drought. It has been recognized that many species in the ecosystem, have been pushed to extinction and many others are at the point of extinction due to interaction with the above mentioned factors [5]. From the bottom-up theory of community ecology, prey population density determines the predator’s population density since predators depend on prey as their only source of food. Thus, any change in prey population density affects the predator population density which means that a decrease in prey population density results in a decrease in predator population density and increase in prey population density results in increase in predator population density [5][6]. Therefore, in order to sustain the prey and predator species in the ecosystem, harvesting rate of prey population should not be either greater or equal or close to intrinsic growth rate of prey but should be at a rate at which the prey species will still survive in the system.

Harvesting in the predator-prey system is the killing of predator and prey species for the purpose of human consumption or trade. The problems of predator-prey systems with harvesting have been explored by many researchers, however most of them have put attention on the optimal harvesting guided entirely by economic profits from harvesting [7][8]. In fact, harvesting of species should be done by considering both economic and biological values of the population [9]. Recently, it has been observed that, the biological species of predator-prey system is harvested and sold for the purpose of making economic profit which slowly decreases the species and finally the ecosystem collapses [10].

Drought can be defined as the absence of water for a long time to result in depletion of soil water and damage to plants [11]. This condition results in the disturbance of the physical and biological structure of an ecological system which in turn alters the availability of resource and the physical environment of the habitat [12][13]. These ecological impacts caused by drought have affected the survival and reproduction of prey and predator species in the ecosystem [6]. For example, in the 1993 Serengeti ecosystem drought, about 30% of the wildebeest and about 40% to 50% of the total park population of large mammals died [5]. These external forces cause a rapid decrease of the population densities and probably the species can go to extinction if control measures are not taken into account. Thus, we need to control the system for the survival of the prey and predator species in the ecosystem. In other words, we need to act upon the problems over-harvesting and drought in order to guarantee that the system behaves as desired [14]. One way of doing so is by the use of optimal control theory which involves the formulation of objective function and dynamical constraints both with time dependent control mechanism.

Several studies have been conducted on optimal control strategies and management policies to keep and protect the ecological species such as those by [5][7][8][10][15][16][17][18][19]. In particular, [20]studied a diseased predator-prey system with stage structure. This study intends to apply optimal control theory to investigate optimal strategies for persistence of predator and prey species by considering the combined effects of harvesting and drought in the system. The optimal control theory is executed by first investigating the effects of over-harvesting and drought on predator-prey model, then extending the model to include time dependent control efforts and finally examining the impact of control efforts, suggested to alleviate the problem of species loss due to over-harvesting and drought. The suggested control efforts are creation of reserve areas with restriction of harvesting for controlling over-harvesting and construction of dams for controlling drought. However, the results obtained from numerical simulations of the model with control efforts showed that, the use of control efforts one at a time or two at a time encourages the increase of the prey and predator species to the optimal population size. Furthermore, the best result is attained when control efforts, creation of reserve areas with restriction of harvesting and construction of dams are simultaneously applied.

2. The Model with Harvesting and Drought

Consider two different populations consisting of prey and predator species. The population densities are respectively represented by x and y. The dynamics of the species interactions is modelled using Holling type II functional response with the following assumptions.

1) In the absence of predators, harvesting and drought, prey population is assumed to grow logistically to the carrying capacity.

2) The rate of increase of the predators depend on the amount of prey biomass predators convert as food.

3) Predators depend on prey population as their only source of food.

From Holling type II functional response, the model is represented by the system of equations as below:

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{{K}_{1}}\right)-\frac{{\alpha}_{1}xy}{x+A}-{\alpha}_{2}x-{d}_{1}x\\ \frac{\text{d}y}{\text{d}t}=-{\mu}_{1}y+\frac{{\beta}_{1}xy}{x+A}-{d}_{2}y,\\ x\left(t\right)>0,y\left(t\right)>0.\end{array}$ (2.1)

where $x\left(t\right)$ represents the prey population density at time t, $y\left(t\right)$ represents the predator population at time t, $r>0$ stands for intrinsic growth rate of prey population, ${K}_{1}$ is the carrying capacity of the environment to prey species, ${\alpha}_{1}>0$ represents the maximum per capita predation rate, ${\beta}_{1}>0$ is the conversion rate, ${\alpha}_{2}>0$ is the rate of harvesting, ${d}_{1}>0$ is the death rate of prey species due to drought, ${\mu}_{1}>0$ is the mortality rate of predator population, ${d}_{2}>0$ is the death rate of predator population due to drought and $A>0$ represents the half saturation of a predator.

2.1. Model Analysis

The solutions of the predator-prey model developed in (2.1) represent the populations of living individuals, and thus, they should be positive and bounded. The following lemma provides details:

Lemma: All the solutions of the system which start in ${\text{IR}}_{+}^{2}$ are uniformly bounded.

Proof: We define a function

$\sigma \left(t\right)=x\left(t\right)+\frac{{\alpha}_{1}y\left(t\right)}{{\beta}_{1}}$ (2.2)

where $\sigma \left(t\right)$ represents the total population of the prey-predator species. Differentiating (2.2) with respect to time gives

$\frac{\text{d}\sigma}{\text{d}t}=\frac{\text{d}x}{\text{d}t}+\frac{{\alpha}_{1}}{{\beta}_{1}}\frac{\text{d}y}{\text{d}t}$ (2.3)

Substitute (2.1) into (2.3) and simplify yields

$\frac{\text{d}\sigma}{\text{d}t}=\left(r-\left({\alpha}_{2}+{d}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}-\frac{{\alpha}_{1}}{{\beta}_{1}}\left({\mu}_{1}+{d}_{2}\right)y$ (2.4)

Choosing arbitrary constant say ${n}_{1}$ and applying it in (2.4) results into

$\begin{array}{l}\frac{\text{d}\sigma}{\text{d}t}=\left(r-\left({\alpha}_{2}+{d}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}-\frac{{\alpha}_{1}}{{\beta}_{1}}\left({\mu}_{1}+{d}_{2}\right)y+{n}_{1}\sigma \left(t\right)-{n}_{1}\sigma \left(t\right)\\ \frac{\text{d}\sigma}{\text{d}t}+{n}_{1}\sigma \left(t\right)=\left(r-\left({\alpha}_{2}+{d}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}-\frac{{\alpha}_{1}}{{\beta}_{1}}\left({\mu}_{1}+{d}_{2}\right)y+{n}_{1}\sigma \left(t\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\left(r-\left({\alpha}_{2}+{d}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}-\frac{{\alpha}_{1}}{{\beta}_{1}}\left({\mu}_{1}+{d}_{2}\right)y+{n}_{1}\left(x+\frac{{\alpha}_{1}}{{\beta}_{1}}y\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\left(r-\left({\alpha}_{2}+{d}_{1}+{n}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}-\frac{{\alpha}_{1}}{{\beta}_{1}}\left({\mu}_{1}+{d}_{2}-{n}_{1}\right)y\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\le \left(r-\left({\alpha}_{2}+{d}_{1}+{n}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}\end{array}$

Thus,

$\frac{\text{d}\sigma}{\text{d}t}+{n}_{1}\sigma \left(t\right)\le \frac{{K}_{1}}{4r}{\left(r-\left({\alpha}_{2}+{d}_{1}\right)+{n}_{1}\right)}^{2}$ (2.5)

where

$\begin{array}{l}\frac{{K}_{1}}{4r}{\left(r-\left({\alpha}_{2}+{d}_{1}\right)+{n}_{1}\right)}^{2}=max\left[\left(r-\left({\alpha}_{2}+{d}_{1}+{n}_{1}\right)\right)x-\frac{r{x}^{2}}{{K}_{1}}\right]\mathrm{.}\hfill \end{array}$

By letting

$\begin{array}{l}\frac{{K}_{1}}{4r}{\left(r-\left({\alpha}_{2}+{d}_{1}\right)+{n}_{1}\right)}^{2}={m}_{1}\hfill \end{array}$

we get

$\begin{array}{l}\frac{\text{d}\sigma}{\text{d}t}+{n}_{1}\sigma \left(t\right)\le {m}_{1}\hfill \end{array}$ (2.6)

Solving Equation (2.6) by integrating factor $I={\text{e}}^{{n}_{1}t}$ yields

$\sigma \left(t\right){\text{e}}^{{n}_{1}t}\le \frac{{m}_{1}}{{n}_{1}}+C{\text{e}}^{-{n}_{1}t}\mathrm{.}$ (2.7)

At $t=0$ , the Equation (2.7) becomes

$\sigma \left(t\right)=\frac{{m}_{1}}{{n}_{1}}+\left(\sigma \left(0\right)-\frac{{m}_{1}}{{n}_{1}}\right){\text{e}}^{-{n}_{1}t}\mathrm{.}$ (2.8)

Thus,

$0\le \sigma \left(t\right)\le \frac{{m}_{1}}{{n}_{1}}+\left(\sigma \left(0\right)-\frac{{m}_{1}}{{n}_{1}}\right){\text{e}}^{-{n}_{1}t}$ (2.9)

As $t\to \infty $ Equation (2.9) becomes

$\begin{array}{l}0\le \sigma \left(t\right)\le \frac{{m}_{1}}{{n}_{1}}\hfill \end{array}$ (2.10)

Therefore, $\sigma \left(t\right)$ is bounded and from positivity of x and y,

$\begin{array}{l}0\le x\le \frac{{m}_{1}}{{n}_{1}}\mathrm{,}\\ 0\le y\le \frac{{m}_{1}}{{n}_{1}}\mathrm{.}\end{array}$ (2.11)

2.1.1. Existence of Equilibrium Points

We start by setting the system of equations of model (2.1) equal to zero and by so doing we get the possible equilibrium points of the system as trivial equilibrium point $\left({A}_{0}\left(\mathrm{0,0}\right)\right)$ , predator free equilibrium point $\left({A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)\right)$ , prey free equilibrium point $\left({A}_{2}\left(\mathrm{0,}{y}^{\mathrm{*}}\right)\right)$ and ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ the co-existence or interior equilibrium point. Thereafter we establish the conditions for existence of each equilibrium point of the system. The condition of existence of the equilibrium point ${A}_{0}\left(\mathrm{0,0}\right)$ is trivial. The existence of the rest of the fixed equilibrium points is described as below:

1) Existence of ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ with ${x}^{*}>0$ .

Let $y=0$ , the system of equations of model (2.1) reduces to

$\begin{array}{l}{x}^{*}\left(r-\frac{r{x}^{*}}{{K}_{1}}-{\alpha}_{2}-{d}_{1}\right)=0,\hfill \end{array}$

which simplifies to

$\begin{array}{l}{x}^{*}=\frac{{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}{r}.\hfill \end{array}$

Thus,

${A}_{1}\left({x}^{*},0\right)={A}_{1}\left(\frac{{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}{r},0\right).$ (2.12)

From the expression of ${x}^{\mathrm{*}}$ we observe that harvesting and drought have negative impacts on the prey growth rate which eventually affects the prey population density. Moreover, for predator free equilibrium point ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ to exist $r-{\alpha}_{2}-{d}_{1}>0$ which implies that $r>{\alpha}_{2}+{d}_{1}$ . Therefore, in the absence of predators the intrinsic growth rate of prey population should be greater than the sum of harvesting and drought rates. Increasing harvesting and drought rates, the prey population decreases which affects the survival of predator species. This is from the fact that predators depend on prey as their only source of food.

2) Existence of ${A}_{2}\left(\mathrm{0,}{y}^{\mathrm{*}}\right)$ with ${y}^{*}>0$ .

Let ${x}^{*}=0$ , the system of equations of model (2.1) reduces to ${y}^{*}\left(-{\mu}_{1}-{d}_{2}\right)=0$ ,

from which we obtain ${y}^{*}=0$ . Thus,

${A}_{2}\left(0,{y}^{*}\right)={A}_{2}\left(0,0\right).$ (2.13)

The result above concur with the assumption that predators depend on prey as their only source of food. Thus, in the absence of prey, predator population becomes extinct.

3) Co-existence equilibrium point ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ .

We equate the system of equations of model (2.1) to zero, and solving for ${x}^{\mathrm{*}}$ and ${y}^{\mathrm{*}}$ yields

$\begin{array}{l}{x}^{*}=-\frac{A\left({d}_{2}+{\mu}_{1}\right)}{-{\beta}_{1}+{d}_{2}+{\mu}_{1}}=\frac{A\left({d}_{2}+{\mu}_{1}\right)}{{\beta}_{1}-{d}_{2}-{\mu}_{1}},\\ {y}^{*}=\frac{A{\beta}_{1}\left({K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)\left({\beta}_{1}-{d}_{2}-{\mu}_{1}\right)-Ar\left({d}_{2}+{\mu}_{1}\right)\right)}{{K}_{1}{\alpha}_{1}{\left({\beta}_{1}-{d}_{2}-{\mu}_{1}\right)}^{2}}.\end{array}$ (2.14)

which gives the co-existence point (interior point).

From the values of interior point in (2.14), we see that predator mortality rate and drought affect the predator’s birth rate ( ${\beta}_{1}$ ) negatively which in turn results into negative effects on predator population density. On the other hand, the co-existence equilibrium point exist if ${\beta}_{1}-{d}_{2}-{\mu}_{1}>0$ implying that ${\beta}_{1}>{d}_{2}+{\mu}_{1}$ . Therefore, in the presence of both populations birth rate of predators should be greater than the sum of rates of mortality and drought. Furthermore, the intrinsic growth rate of prey should be greater than the sum of harvesting rate and drought rate of prey for co-existence equilibrium point to exist. However, increasing death rate of predator population results in the increase of prey population density.

2.1.2. Local Stability Analysis of the Equilibrium Points

The stability properties of the equilibrium points are analyzed by computing the Jacobian matrix and determining the eigenvalues of the Jacobian matrix of each fixed point. The equilibrium points are asymptotically stable if the real parts of the eigenvalues of each Jacobian matrix is negative. From the system of equations of model (2.1) the general Jacobian matrix of the equations is given by:

$J\left({A}_{i}\right)=\left[\begin{array}{cc}\frac{\partial {g}_{1}}{\partial x}& \frac{\partial {g}_{1}}{\partial y}\\ \frac{\partial {g}_{2}}{\partial x}& \frac{\partial {g}_{2}}{\partial y}\end{array}\right]$

which provides

$J\left({A}_{i}\right)=\left[\begin{array}{cc}r\left(1-\frac{{x}^{\mathrm{*}}}{{K}_{1}}\right)-\frac{r{x}^{\mathrm{*}}}{{K}_{1}}-\frac{{\alpha}_{1}{y}^{\mathrm{*}}}{{x}^{\mathrm{*}}+A}+\frac{{\alpha}_{1}{x}^{\mathrm{*}}{y}^{\mathrm{*}}}{{\left({x}^{\mathrm{*}}+A\right)}^{2}}-{\alpha}_{2}-{d}_{1}& -\frac{{\alpha}_{1}{x}^{\mathrm{*}}}{{x}^{\mathrm{*}}+A}\\ \frac{{\beta}_{1}{y}^{\mathrm{*}}}{{x}^{\mathrm{*}}+A}-\frac{{\beta}_{1}{x}^{\mathrm{*}}{y}^{\mathrm{*}}}{{\left({x}^{\mathrm{*}}+A\right)}^{2}}& -{\mu}_{1}+\frac{{\beta}_{1}{x}^{\mathrm{*}}}{{x}^{\mathrm{*}}+A}-{d}_{2}\end{array}\right]\mathrm{.}$ (2.15)

1) For ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)={A}_{1}\left[\frac{{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}{r}\mathrm{,0}\right]$ ,

the Jacobian matrix evaluated at this equilibrium point is given by

$J\left({A}_{1}\right)=\left[\begin{array}{cc}-\left(r-{\alpha}_{2}-{d}_{1}\right)& -\frac{{K}_{1}{\alpha}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}{Ar+{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}\\ 0& {\beta}_{1}-\frac{Ar{\beta}_{1}}{Ar+{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}-{d}_{2}-{\mu}_{1}\end{array}\right].$

The eigenvalues of the Jacobian matrix $J\left({A}_{1}\right)$ are $-\left(r-{\alpha}_{2}-{d}_{1}\right)$ and

${\beta}_{1}-\frac{Ar{\beta}_{1}}{Ar+{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}-{d}_{2}-{\mu}_{1}$ .

Therefore, the equilibrium point ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ is locally asymptotically stable if

it satisfies the following condition: ${\beta}_{1}-\frac{Ar{\beta}_{1}}{Ar+{K}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}-{d}_{2}-{\mu}_{1}<0$ .

which is negative if ${\beta}_{1}<{d}_{2}+{\mu}_{1}$ . Thus, the predator free equilibrium point is locally asymptotically stable if predator population goes to extinction.

2) The corresponding Jacobian matrix of the equations evaluated at the equilibrium point ${A}_{2}\left(0,{y}^{*}\right)={A}_{2}\left(0,0\right)$ is

$J\left({A}_{2}\right)=\left[\begin{array}{cc}r-{\alpha}_{2}-{d}_{1}& 0\\ 0& -{\mu}_{1}-{d}_{2}\end{array}\right]$

The eigenvalues for the Jacobian matrix $J\left({A}_{2}\right)$ are $r-{\alpha}_{2}-{d}_{1}$ and $-\left({\mu}_{1}+{d}_{2}\right)$ . If $r>{\alpha}_{2}+{d}_{1}$ , then the equilibrium point ${A}_{2}\left(\mathrm{0,}{y}^{\mathrm{*}}\right)$ is unstable saddle point. That is, to say the prey population is growing. However, if $r<{\alpha}_{2}+{d}_{1}$ then the prey population goes to extinction which results to extinction of predator population under the assumption that predators depends on prey as their only source of food. Thus, the prey free equilibrium point is locally asymptotically stable if prey population goes to extinction.

3) For co-existence fixed point ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$

Consider the following Jacobian matrix

$J\left({A}_{3}\right)=\left(\begin{array}{cc}{a}_{11}& {a}_{12}\\ {a}_{21}& {a}_{22}\end{array}\right)$

where

$\begin{array}{l}{a}_{11}=\frac{\left({d}_{2}+{\mu}_{1}\right)\left(\left(-Ar+{K}_{1}\left(r-{\alpha}_{2}\right)\right){\beta}_{1}+\left({\alpha}_{2}{K}_{1}-\left(A+{K}_{1}\right)r\right)\left({d}_{2}+{\mu}_{1}\right)-M\right)}{{K}_{1}{\beta}_{1}\left({\beta}_{1}-{d}_{2}-{\mu}_{1}\right)}\\ {a}_{12}=-\frac{{\alpha}_{1}\left({d}_{2}+{\mu}_{1}\right)}{{\beta}_{1}}\\ {a}_{21}=\frac{{K}_{1}\left(r-{\alpha}_{2}\right){\beta}_{1}+\left({\alpha}_{2}{K}_{1}-\left(A+{K}_{1}\right)r\right)\left({d}_{2}+{\mu}_{1}\right)-{K}_{1}{d}_{1}\left({\beta}_{1}-{d}_{2}-{\mu}_{1}\right)}{{\alpha}_{1}{K}_{1}}\\ {a}_{22}=0\end{array}$

and $M={K}_{1}{d}_{1}\left({\beta}_{1}-{d}_{2}-{\mu}_{1}\right)$

The stability of this equilibrium point is stated using the trace/determinant technique as follows: Suppose the Jacobian matrix evaluated at the co-existence equilibrium point has the characteristic polynomial equation

${\lambda}^{2}-P\lambda +Q=0,$ (2.16)

Such that $P=tr\left(J\left({A}_{3}\right)\right)={a}_{11}+{a}_{22}={a}_{11}$ and $Q=Det\left(J\left({A}_{3}\right)\right)={a}_{11}{a}_{22}-{a}_{21}{a}_{12}=-{a}_{21}{a}_{12}.$

Then, the co-existence equilibrium point is locally asymptotically stable or stable spiral if $P<0$ and $Q>0$ . However, the interior equilibrium point is a center (neutrally stable) if $P=0$ and $Q>0$ . Moreover, if the co-existence equilibrium point is locally asymptotically stable or stable spiral or neutrally stable, it implies that there is a stable dynamic relationship between predator and prey.

2.1.3. Global Stability Analysis

The global stability of equilibrium points ${A}_{1}$ and ${A}_{3}$ is shown by linearizing the system of equations of model (2.1) and defining appropriate Lyapunov function to separately describe each equilibrium point. The linearization process is done using Jacobian technique such that

$\frac{\text{d}{X}_{i}}{\text{d}t}=J\left({A}_{i}\right){X}_{i},$ (2.17)

where $J\left({A}_{i}\right)$ is the Jacobian matrix and ${X}_{i}$ is a small perturbation on ${x}_{i}$ . Therefore the system of equations of model (2.1) reduces to the following linear system of equations;

$\begin{array}{l}\frac{\text{d}X}{\text{d}t}=\left[r\left(1-\frac{{x}^{*}}{{K}_{1}}\right)-\frac{r{x}^{*}}{{K}_{1}}-\frac{{\alpha}_{1}{y}^{*}}{{x}^{*}+A}+\frac{{\alpha}_{1}{x}^{*}{y}^{*}}{{\left({x}^{*}+A\right)}^{2}}-{\alpha}_{2}-{d}_{1}\right]X-P\\ \frac{\text{d}Y}{\text{d}t}=\left[\frac{{\beta}_{1}{y}^{*}}{{x}^{*}+A}-\frac{{\beta}_{1}{x}^{*}{y}^{*}}{{\left({x}^{*}+A\right)}^{2}}\right]X+\left[-{\mu}_{1}+\frac{{\beta}_{1}{x}^{*}}{{x}^{*}+A}-{d}_{2}\right]Y\end{array}$ (2.18)

where $P=\frac{{\alpha}_{1}{x}^{*}}{{x}^{*}+A}Y$ .

The Lyapunov function is chosen as

$V\left(X,Y\right)=\frac{{X}^{2}}{2}+\frac{{Y}^{2}}{2}.$ (2.19)

The function $V\left(X,Y\right)$ is a positive definite function since $V\left(X,Y\right)>0$ for any values of $\left(X\mathrm{,}Y\right)$ and it is minimum at the origin. That is, $V\left(0,0\right)=0$ . The time derivative of $V\left(X,Y\right)$ is given by

$\frac{\text{d}V\left(X,Y\right)}{\text{d}t}=\frac{\partial V}{\partial X}\frac{\text{d}X}{\text{d}t}+\frac{\partial V}{\partial Y}\frac{\text{d}Y}{\text{d}t}.$ (2.20)

By substitution of (2.18) into (2.20) and differentiate V partially, we obtain the following

$\begin{array}{c}\frac{\text{d}V\left(X,Y\right)}{\text{d}t}=X\left[\left(r\left(1-\frac{{x}^{*}}{{K}_{1}}\right)-\frac{r{x}^{*}}{{K}_{1}}-\frac{{\alpha}_{1}{y}^{*}}{{x}^{*}+A}+\frac{{\alpha}_{1}{x}^{*}{y}^{*}}{{\left({x}^{*}+A\right)}^{2}}-{\alpha}_{2}-{d}_{1}\right)X-S\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+Y\left[\left(\frac{{\beta}_{1}{y}^{*}}{{x}^{*}+A}-\frac{{\beta}_{1}{x}^{*}{y}^{*}}{{\left({x}^{*}+A\right)}^{2}}\right)X+\left(-{\mu}_{1}+\frac{{\beta}_{1}{x}^{*}}{{x}^{*}+A}-{d}_{2}\right)Y\right]\end{array}$ (2.21)

where $S=\frac{{\alpha}_{1}{x}^{*}}{{x}^{*}+A}Y$ .

1) For fixed point ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$

Here, Equation (2.12) is used in Equation (2.21) to get

$\begin{array}{c}\frac{\text{d}V\left(X,Y\right)}{\text{d}t}=-\left(r-{\alpha}_{2}-{d}_{1}\right){X}^{2}-\frac{{K}_{1}{\alpha}_{1}\left(r-{\alpha}_{2}-{d}_{1}\right)}{\left(A+K\right)r-{K}_{1}{\alpha}_{2}-K{d}_{1}}XY\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({\beta}_{1}-\frac{Ar{\beta}_{1}}{\left(A+{K}_{1}\right)r-{K}_{1}{\alpha}_{2}-{K}_{1}{d}_{1}}-\left({d}_{2}+{\mu}_{1}\right)\right){Y}^{2}\end{array}$ (2.22)

From (2.22), the equilibrium point ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ is globally asymptotically stable if the following condition is satisfied: That is,

$\left({\beta}_{1}-\frac{Ar{\beta}_{1}}{\left(A+K\right)r-{K}_{1}{\alpha}_{2}-{K}_{1}{d}_{1}}-\left({d}_{2}+{\mu}_{1}\right)\right){Y}^{2}<0,$

which is negative if ${\beta}_{1}<{d}_{2}+{\mu}_{1}$ , implying that, the equilibrium point ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ is globally stable if predators go to extinction.

2) For equilibrium point ${A}_{2}\left(\mathrm{0,}{y}^{\mathrm{*}}\right)$ .

In this subsection, we use Equation (2.13) in Equation (2.21) to obtain

$\frac{\text{d}V\left(X,Y\right)}{\text{d}t}=\left(r-{\alpha}_{2}-{d}_{1}\right){X}^{2}-\left({d}_{2}+{\mu}_{1}\right){Y}^{2}$ (2.23)

From Equation (2.23), the fixed point ${A}_{2}\left(\mathrm{0,}{y}^{\mathrm{*}}\right)$ is globally asymptotically stable if $\left(r-{\alpha}_{2}-{d}_{1}\right){X}^{2}<0$ , which implies $r<{\alpha}_{2}+{d}_{1}$ . Thus, the prey free equilibrium point is globally asymptotically stable if prey population goes to extinction.

3) For steady state point ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ .

Here we substitute the values of interior point (2.14) into Equation (2.21) to obtain

$\frac{\text{d}V\left(X\mathrm{,}Y\right)}{\text{d}t}={a}_{11}{X}^{2}+\left({a}_{12}+{a}_{21}\right)XY\mathrm{,}$ (2.24)

with usual notations for ${a}_{11}\mathrm{,}{a}_{12}$ and ${a}_{21}$ . Therefore the point is globally stable if the condition below holds,

$\frac{\text{d}V\left(X,Y\right)}{\text{d}t}=\left({a}_{11}{X}^{2}+\left({a}_{12}+{a}_{21}\right)XY\right)<0.$ (2.25)

3. The Model with Control

We introduce into system of equations of model (2.1), time dependent control efforts on harvesting and drought to alleviate the loss of species in the predator-prey system. Let ${u}_{1}\left(t\right)$ represents over-harvesting control strategy (Creation of reserve areas with restriction of harvesting) and ${u}_{2}\left(t\right)$ represents the drought control strategy (construction of dams). Thus, the system of equations of model (2.1) becomes:

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{{K}_{1}}\right)-\frac{{\alpha}_{1}xy}{x+A}-\left(1-{u}_{1}\left(t\right)\right){\alpha}_{2}x-\left(1-{u}_{2}\left(t\right)\right){d}_{1}x\\ \frac{\text{d}y}{\text{d}t}=-{\mu}_{1}y+\frac{{\beta}_{1}xy}{x+A}-\left(1-{u}_{2}\left(t\right)\right){d}_{2}y\end{array}$ (3.1)

3.1. Analysis of the Optimal Control

Here we construct an objective function that provides the optimal population size of the prey-predator species at minimum costs for over-harvesting and drought control strategies. The objective function takes the following form, that is

$J\left({u}_{i}\right)=\mathrm{max}\left[G\left({x}_{i}\left({T}_{1}\right),{T}_{1}\right)-{\displaystyle {\int}_{{T}_{0}}^{{T}_{1}}}\left(F\left({x}_{i}\left(t\right),{u}_{i}\left(t\right),t\right)\right)\text{d}t\right]$ (3.2)

subject to

$\frac{\text{d}{x}_{i}}{\text{d}t}={g}_{i}\left(t\mathrm{,}{u}_{i}\left(t\right)\mathrm{,}{x}_{i}\left(t\right)\right)\mathrm{,}$

where ${x}_{i}\left({T}_{0}\right)={x}_{i0}$ and $0\le {u}_{i}\le 1$ for $t\in \left[{T}_{0}\mathrm{,}{T}_{1}\right]$ . The terms $G\left({x}_{i}\left({T}_{1}\right)\right)$ and $F\left({x}_{i}\left(t\right)\mathrm{,}{u}_{i}\left(t\right)\mathrm{,}t\right)$ represent the predator-prey populations to be optimized at terminal time of control and total cost of control respectively.

Therefore, from (3.2) the objective function becomes:

$J\left(U\right)=\underset{{u}_{1},{u}_{2}}{\mathrm{max}}\left[{D}_{1}x\left({T}_{1}\right)+{D}_{2}y\left({T}_{1}\right)-{\displaystyle {\int}_{0}^{{T}_{1}}}\left(\frac{{C}_{1}{u}_{1}^{2}}{2}+\frac{{C}_{2}{u}_{2}^{2}}{2}\right)\text{d}t\right]$ (3.3)

subject to

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{{K}_{1}}\right)-\frac{{\alpha}_{1}xy}{x+A}-\left(1-{u}_{1}\left(t\right)\right){\alpha}_{2}x-\left(1-{u}_{2}\left(t\right)\right){d}_{1}x\\ \frac{\text{d}y}{\text{d}t}=-{\mu}_{1}y+\frac{{\beta}_{1}xy}{x+A}-\left(1-{u}_{2}\left(t\right)\right){d}_{2}y\end{array}$

$x\left({T}_{0}\right)={x}_{0}$ , $y\left({T}_{0}\right)={y}_{0}$ and $0\le {u}_{1}\le \mathrm{1,0}\le {u}_{2}\le 1$ for $t\in \left[\mathrm{0,}{T}_{1}\right]$ ; ${u}_{1}\mathrm{,}{u}_{2}\in U$ .

However, the terms ${D}_{1}x\left({T}_{1}\right)$ and ${D}_{2}y\left({T}_{1}\right)$ represent the prey and predator

populations to be optimized at terminal time of control and $\frac{{C}_{1}{u}_{1}^{2}}{2}$ and $\frac{{C}_{2}{u}_{2}^{2}}{2}$

are total control costs for over-harvesting and drought respectively. The costs weights ${C}_{1}\mathrm{,}{C}_{2}$ and state weights ${D}_{1}\mathrm{,}{D}_{2}$ are all positive constants. In here, the aim is to pick up ${u}_{1}^{\mathrm{*}}$ and ${u}_{2}^{\mathrm{*}}$ such that,

$J\left({u}_{1}^{\mathrm{*}}\mathrm{,}{u}_{2}^{\mathrm{*}}\right)=\underset{U}{max}\left(J\left({u}_{1}\mathrm{,}{u}_{2}\right)\right)\mathrm{,}$ (3.4)

with $0\le {u}_{1}\le \mathrm{1,0}\le {u}_{2}\le 1$ for $t\in \left[\mathrm{0,}{T}_{1}\right]$ .

3.1.1. Existence of Optimal Control

The aim is to show that, the optimal control problem formulated in (3.3) has at least one solution before trying to solve for the optimal control values. The following theorem explains. Theorem: Given an optimal control problem in (3.3) with ${u}_{1}$ and ${u}_{2}$ as control variables, then there exist ${u}_{1}^{\mathrm{*}}\mathrm{,}{u}_{2}^{\mathrm{*}}\in U$ (optimal control set) such that,

$J\left({u}_{1}^{\mathrm{*}}\mathrm{,}{u}_{2}^{\mathrm{*}}\right)=\underset{U}{max}\left(J\left({u}_{1}\mathrm{,}{u}_{2}\right)\right)\mathrm{,}$ (3.5)

Proof: The proof for the existence of optimal control provided by [5][21][22]is also appropriate in here.

1) Model Equation (3.1) with control are linear in control variables.

2) The control $U$ is convex, closed and bounded set.

3) The integrand $-\left(\frac{{C}_{1}{u}_{1}^{2}}{2}+\frac{{C}_{2}{u}_{2}^{2}}{2}\right)$ of the objective function (3.3) is concave

in $U$ .

3.1.2. Characterization of the Optimal Control

The calculation of the optimal control strategy is done by the application of Pontryagin’s maximum principle (PMP) which provides the Hamiltonian H and necessary conditions with which the optimal control and the co-state variables must satisfy. From the objective function (3.3), the Hamiltonian H is given by

$\begin{array}{c}H=-\left(\frac{{C}_{1}{u}_{1}^{2}}{2}+\frac{{C}_{2}{u}_{2}^{2}}{2}\right)+{R}_{1}\left[rx\left(1-\frac{x}{{K}_{1}}\right)-\frac{{\alpha}_{1}xy}{x+A}-\left(1-{u}_{1}\left(t\right)\right){\alpha}_{2}x-\left(1-{u}_{2}\left(t\right)\right){d}_{1}x\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{R}_{2}\left[-{\mu}_{1}y+\frac{{\beta}_{1}xy}{x+A}-\left(1-{u}_{2}\left(t\right)\right){d}_{2}y\right],\end{array}$ (3.6)

where ${R}_{1}$ and ${R}_{2}$ represents the co-state variables or adjoint variables. Using Pontryagin’s maximum principle (PMP) [23]and the existence result for the optimal control variables, we state a proposition as follows:

Proposition: Let ${u}_{1}^{\mathrm{*}}$ and ${u}_{2}^{\mathrm{*}}$ be optimal control variables that optimize the objective function in (3.3), then there exist adjoint variables ${R}_{1}$ and ${R}_{2}$

satisfying the adjoint condition $\frac{\text{d}{R}_{i}}{\text{d}t}=\frac{\partial H}{\partial {x}_{i}}$ with ${R}_{i}\left(T\right)={D}_{i}$ ; $i=1,2$ .

Proof: Using the adjoint condition, the set of adjoint equations becomes:

$\begin{array}{l}\frac{\text{d}{R}_{1}}{\text{d}t}=-\frac{\partial H}{\partial x}=-{R}_{1}[r\left(1-\frac{x}{{K}_{1}}\right)-\frac{rx}{{K}_{1}}-\frac{{\alpha}_{1}y}{x+A}+\frac{{\alpha}_{1}xy}{{\left(x+A\right)}^{2}}-\left(1-{u}_{1}\left(t\right)\right){\alpha}_{2}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\begin{array}{c}\\ \end{array}-\left(1-{u}_{2}\left(t\right)\right){d}_{1}]-{R}_{2}\left[\frac{{\beta}_{1}y}{x+A}-\frac{{\beta}_{1}xy}{{\left(x+A\right)}^{2}}\right]\\ \frac{\text{d}{R}_{2}}{\text{d}t}=-\frac{\partial H}{\partial y}=-{R}_{1}\left[-\frac{{\alpha}_{1}x}{x+A}\right]-{R}_{2}\left[-{\mu}_{1}+\frac{{\beta}_{1}x}{x+A}-\left(1-{u}_{2}\left(t\right)\right){d}_{2}\right]\end{array}$ (3.7)

with transversality conditions:

$\begin{array}{l}{R}_{1}\left({T}_{1}\right)=\frac{\text{d}}{\text{d}x}\left({D}_{1}x\left({T}_{1}\right)\right)={D}_{1},\\ {R}_{2}\left({T}_{1}\right)=\frac{\text{d}}{\text{d}y}\left({D}_{2}y\left({T}_{1}\right)\right)={D}_{2}.\end{array}$ (3.8)

Using optimality condition,

$\frac{\partial H}{\partial {u}_{i}}=0$ at ${u}_{i}^{\mathrm{*}}$ , we have $\frac{\partial H}{\partial {u}_{1}}$ at ${u}_{1}^{\mathrm{*}}$ and $\frac{\partial H}{\partial {u}_{2}}$ at ${u}_{2}^{\mathrm{*}}$ . Thus, $\frac{\partial H}{\partial {u}_{1}}=-{C}_{1}{u}_{1}^{*}+{R}_{1}{\alpha}_{2}x=0$ which gives

${u}_{1}^{*}=\frac{{R}_{1}{\alpha}_{1}x}{{C}_{1}}$ (3.9)

Also,

$\frac{\partial H}{\partial {u}_{2}}=-{C}_{2}{u}_{2}^{*}+{R}_{1}{d}_{1}x+{R}_{2}{d}_{2}y=0$

Thus,

${u}_{2}^{*}=\frac{{R}_{1}{d}_{1}x+{R}_{2}{d}_{2}y}{{C}_{2}}$ (3.10)

Solving for x, y, ${R}_{1}$ and ${R}_{2}$ , we substitute ${u}_{1}^{\mathrm{*}}$ and ${u}_{2}^{\mathrm{*}}$ in the following system of equations below:

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{{K}_{1}}\right)-\frac{{\alpha}_{1}xy}{x+A}-\left(1-{u}_{1}\left(t\right)\right){\alpha}_{2}x-\left(1-{u}_{2}\left(t\right)\right){d}_{1}x,x\left(0\right)={x}_{0}\\ \frac{\text{d}y}{\text{d}t}=-{\mu}_{1}y+\frac{{\beta}_{1}xy}{x+A}-\left(1-{u}_{2}\left(t\right)\right){d}_{2}y,y\left(0\right)={y}_{0}\\ \frac{\text{d}{R}_{1}}{\text{d}t}=-{R}_{1}\left[r\left(1-\frac{x}{{K}_{1}}\right)-\frac{rx}{{K}_{1}}-\frac{{\alpha}_{1}y}{x+A}+\frac{{\alpha}_{1}xy}{{\left(x+A\right)}^{2}}-\left(1-{u}_{1}\left(t\right)\right){\alpha}_{2}-{W}_{1}\right]\end{array}$

$\begin{array}{l}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-{R}_{2}\left[\frac{{\beta}_{1}y}{x+A}-\frac{{\beta}_{1}xy}{{\left(x+A\right)}^{2}}\right],{R}_{1}\left({T}_{1}\right)={D}_{1}\\ \frac{\text{d}{R}_{2}}{\text{d}t}=-{R}_{1}\left[-\frac{{\alpha}_{1}x}{x+A}\right]-{R}_{2}\left[-{\mu}_{1}+\frac{{\beta}_{1}x}{x+A}-\left(1-{u}_{2}\left(t\right)\right){d}_{2}\right],{R}_{2}\left({T}_{2}\right)={D}_{2},\end{array}$ (3.11)

where, ${W}_{1}=\left(1-{u}_{2}\left(t\right)\right){d}_{1}$ .

Taking back the values obtained in (3.11) into the expressions of ${u}_{1}^{\mathrm{*}}$ and ${u}_{2}^{\mathrm{*}}$ in (3.9) and (3.10) we get the characterization of the optimal control as:

$\begin{array}{l}{u}_{1}^{*}=\mathrm{min}\left\{1,\mathrm{max}\left(0,\frac{{R}_{1}{\alpha}_{2}x}{{C}_{1}}\right)\right\}\\ {u}_{2}^{*}=\mathrm{min}\left\{1,\mathrm{max}\left(0,\frac{{R}_{1}{d}_{1}x+{R}_{2}{d}_{2}y}{{C}_{2}}\right)\right\},\end{array}$ (3.12)

where ${R}_{1}$ and ${R}_{2}$ are the solutions of the system of adjoint Equation (3.7).

4. Numerical Results and Discussion

We verify the equilibrium points of system equations of model (2.1), using phase diagrams representation. The aim is to show the dynamical behavior of the equilibrium points discussed in the theoretical part. Using parameter values as in Table 2, the phase diagrams for equilibrium points ${A}_{2}\left(0,{y}^{*}\right)={A}_{2}\left(0,0\right)$ , ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ and ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ are drawn as follows:

Figure 1 indicates that, in the presence of over-harvesting and drought the dynamical behavior of equilibrium points ${A}_{2}\left(0,{y}^{*}\right)$ and ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ is unstable saddle and the dynamical behavior of co-existence equilibrium point ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ is a spiral unstable surrounded by a stable limit cycle. The equilibrium points ${A}_{2}\left(0,{y}^{*}\right)={A}_{2}\left(0,0\right)$ and ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ are clearly shown in a top sub-figure of Figure 1, while the equilibrium point ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ is well shown by a bottom sub-figure of Figure 1.

4.1. The Effects of Varying Harvesting Rate on Prey and Predator Population Densities

Figure 2(a) shows the effects of varying harvesting rate to prey population density. From the result we observe that, as the harvesting rate of prey increases, the prey population decreases. The decrease in prey population proves that harvesting has strong negative impact on predator-prey population. Figure 2(b) shows the effects of varying harvesting rate of prey on predator population density. The result indicates that, as harvesting rate of prey increases, predator population decreases. However, the decrease in predator population while they are not affected by over-harvesting directly, is due to the fact that predators depend on prey as their only source of food. Thus, any effect on prey species results into indirect effect on predator species.

4.2. The Effects of Varying Harvesting and Drought Rates on Prey and Predator Population Densities

Figure 3(a) is the effects of increasing harvesting and drought rates of prey at a

Figure 1. Phase diagrams showing dynamical behavior of ${A}_{2}\left(0,{y}^{*}\right)$ , ${A}_{3}\left({x}^{\mathrm{*}}\mathrm{,}{y}^{\mathrm{*}}\right)$ and ${A}_{1}\left({x}^{\mathrm{*}}\mathrm{,0}\right)$ .

time, on prey population density. The result shows that, there is a rapid decrease in prey population. The rapid decrease in prey population is due to the combinations of the problems (harvesting and drought) on prey species both with strong negative impact on the prey species.

Figure 3(b) shows the effects of increasing harvesting and drought rates of prey at a time, on predator population density. From the result we observe that, there is a rapid decrease in predator population density. The rapid decrease in

Figure 2. Simulations of a predator-prey model showing the impacts of varying harvesting rate to prey and predator population densities.

predator population density while the effects are on prey individuals is due to dependence of predators on prey as their only source of food. Thus, decreasing prey population rapidly results into rapid decrease in predator population density.

Figure 3. Simulations of a predator-prey model showing the impacts of varying harvesting and drought rates on prey and predator densities.

4.3. The Effects of Varying Drought Rates on Prey and Predator Population Densities

From Figure 4(a) we observe that, the increase of drought, results into decrease of prey population density. The decrease of prey population, is due to the fact

Figure 4. Simulations of a prey-predator model showing the impacts of varying drought rates on prey and predator densities.

that drought affect the population negatively. From Figure 4(b) we observe that, the increase of drought to predator population decreases predator population density. The decrease in predator population shows that, drought has strong negative impact on predator-prey system.

4.4. Numerical Results of the Model with Control

The optimal control strategy can numerically be solved by several numerical techniques using parameter values. A forward-backward sweep method (FBSM) is one of the numerical techniques that can be used to solve an optimal control problem. The technique uses Runge-Kutta 4th order numerical method (RK4) in solving the optimal control problem. The following scholars [5][18][23], suggested the method to be executed as follows:

1) Make an initial guess for control vector and use initial conditions for state vector to solve for the state variables forward in time using Runge-Kutta 4th order numerical method (RK4).

2) Using the new sate values, transversality conditions and the guessed values for control vector, solve for adjoint vector backward in time using RK4.

3) The obtained values for state and adjoint variables are entered on the characterization of the optimal control (3.12) to update the control vector which becomes a new value for the control vector.

4) Stop the process when the values of the control variables in the current and previous iterations are sufficiently close.

The investigation of the impacts of the control efforts on the predator-prey system is studied numerically through different combinations of the control strategies. The following Table 1 and Table 2 are respectively indicating different combinations of the control strategies and parameter values of the predator-prey model.

Table 1. Control strategies for prey-predator model (3.1).

Table 2. Parameter values for prey-predator model (3.1).

According to [5], the constants for state and control variables are chosen depending on their relative importance and relative implications of cost used for controlling the problems. Thus, the initial state variables and other constants for state and control variables are chosen as follows: ${C}_{1}=110,{C}_{2}=46,{D}_{1}=45$ and ${D}_{2}=36$ and the initial state variables are $x\left(0\right)=48$ and $y\left(0\right)=24$ . The simulations for different strategies are carried out as below:

4.4.1. Control Strategy I: Creation of Reserve Areas with Restriction of Harvesting for Controlling Over-Harvesting

From Figure 5(a), the impact of control strategy for controlling over-harvesting

Figure 5. Simulations of a predator-prey model affected by over-harvesting and drought showing the impacts of creating reserve areas with restriction of harvesting to prey and predator population densities.

u_{1} is shown by using it to optimize the objective function
$J\left(U\right)$ while setting the control strategy for drought u_{2} equal to zero. The result shows that, at the final time of control the prey population increases from 24 to about 128 individuals. These results indicate that, over-harvesting has strong negative impact on the population dynamics of prey-predator population densities.

In Figure 5(b), optimal control strategy for over-harvesting, u_{1} is used to optimize the objective function
$J\left(U\right)$ while optimal control strategy for drought, u_{2} is set to zero. The result shows that, Despite of predators having enough food due to increased prey species as shown in Figure 5(a), there is a very small increase in a predator population density. The very small increase in predator population is caused by the effects of drought which is in this case uncontrolled.

4.4.2. Control Strategy II: Construction of Dams for Mitigating Drought Effects

From Figure 6(a) above the impact of construction of dams u_{2} for mitigating drought effects is shown by using the control strategy to optimize the objective function
$J\left(U\right)$ while setting the control strategy for over-harvesting u_{1} equal to zero. The result shows that after control there is increase in small number in prey population from about 24 (final state before control) to about 30 species only, at the final time of control. The increase in small number of prey population density while drought has been controlled is due to increased predators (Figure 6(b)) after controlling drought, and over-harvesting which is in this case uncontrolled.

In Figure 6(b) we show the impact of construction of dams, u_{2} to predator-prey population densities. We observe that, the predator population density has increased from 24 individuals to about 76 individuals. The strategy is much successful in optimizing predator population. This is due to the assumption that, drought is the only external factor that affects predator population density.

4.4.3. Control Strategy III: Combination of Creation of Reserve Areas with Restriction of Harvesting and Construction of Dams

From Figure 7(a), the effects of both control Strategies u_{1} and u_{2} are shown by using them simultaneously to optimize the objective function
$J\left(U\right)$ . The result of control showed that, prey population increases from 24 to about 150 species. Thus, applying both control strategies at a time to predator-prey densities, encourages the growth of prey in the predator-prey system.

From Figure 7(b), both control strategies are simultaneously applied to optimize the objective function $J\left(U\right)$ . At the final time of control, We observe that predator population increases from 42 to about 84 individuals. This information from the Figure implies that, the use of both control strategies at a time on prey-predator system, promote the growth of predator population in the predator-prey system.

Figure 6. Simulations of a predator-prey model affected by over-harvesting and drought showing the impacts of construction of dams to prey and predator population densities.

5. Conclusions

The aim of this study was to investigate the effects of over-harvesting and drought, and to assess the impact of control measures to these external forces. We developed a predator-prey model which incorporates harvesting and drought without control. Thereafter we carried out some analysis (theoretically and numerically) of the developed model to describe how harvesting and drought affect the prey and predator species. Furthermore, the model was modified to include the time dependent control efforts on over-harvesting and on drought. The numerical analysis to investigate the impact of the time dependent control

Figure 7. Simulations of a predator-prey model affected by over-harvesting and drought showing the impacts of combining control strategies to prey and predator population densities.

efforts on the predator-prey model was also carried out.

However, from the numerical simulations we observe that, using control strategies one at a time or two at a time encourages the increase of prey and predator population densities to the optimal population size. Furthermore, the best results is attained when the control strategies are applied at a time. Thus, in order to keep and protect the prey and predator species from over-harvesting and drought, the stakeholders should choose to apply for control strategies to both over-harvesting and drought simultaneously.

References

[1] Begon, M., Townsend, C.R. and Harper, J.L. (2006) ECOLOGY: From Individuals to Ecosystem. CPI Bath Press, UK.

[2] Taylor, R.J. (1984) Predation. Champman and Hall, New York.

https://doi.org/10.1007/978-94-009-5554-7

[3] Berryman, A.A. (1992) The Origins and Evolution of Predator Prey Theory. Ecology, 73, 1530-1535.

https://doi.org/10.2307/1940005

[4] Dubey, B. (2007) A Prey-Predator Model with a Reserved Area. Nonlinear Analysis: Modeling and Control, 12, 479-494.

[5] Sagamiko, T.D., Shaban, N., Nahonyo, C.L. and Makinde, O.D. (2015) Optimal Control of a Threatened Wildebeest-Lion Prey-Predator System in the Serengeti Ecosystem. Open Journal of Ecology, 5, 110-119.

https://doi.org/10.4236/oje.2015.54010

[6] Sperry, J.H. and Weatherhead, P.J. (2008) Prey-Mediated Effects of Drought on Condition and Survival of a Terrestrial Snake. Ecological Society of America, 89, 2770-2776.

[7] Chakraborty, K., Das, S. and Kar, T.K. (2011) Optimal Control of Effort of a Stage Structured Prey Predator Fishery Model with Harvesting. Nonlinear Analysis: Real World Applications, 12, 3452-3467.

https://doi.org/10.1016/j.nonrwa.2011.06.007

[8] Chakraborty, K., Chakraborty, M. and Kar, T.K. (2011) Optimal Control of Harvesting and Bifurcation of a Prey-Predator Model with Stage Structure. Applied Mathematics and Computation, 217, 8778-8792.

https://doi.org/10.1016/j.amc.2011.03.139

[9] Kaiyuan, L. and Lansun, C. (2007) Harvesting Control for a Stage-Structured Predator-Prey Model with Ivlev’s Functional Response and Impulsive Stocking on Prey. Discrete Dynamics in nature and Society, 2007, Article ID: 86482.

https://doi.org/10.1155/2007/86482

[10] Kar, T.K. and Ghosh, B. (2010) Bifurcations and Feedback Control of a Stage Structure Exploited Prey Predator System. International Journal of Engineering, Science and Technology, 2, 131-141.

[11] Hanson, P.J. and Weltzin, J.F. (2000) Drought Disturbance from Climate Change: Response of United States Forests. Science of the Total Environment, 262, 205-220.

https://doi.org/10.1016/S0048-9697(00)00523-4

[12] Pickett, S.T.A. and White, P.S. (1985) The Ecology of Natural Disturbance and Patch Dynamics. Academic Press, London.

[13] Pickett, S.T.A., Wu, J.G. and Cadenasso, M.L. (1999) Patch Dynamics and the Ecology of Disturbed Ground. Elsevier, Amsterdam, 707-722.

[14] Fernandez-Cara, E. and ZuaZua, E. (2003) Control Theory: History, Mathematical Achievements and Perspectives. SeMA Journal, 26, 79-140.

[15] Bera, S.P., Maiti, A. and Samanta, G.P. (2014) A Prey-Predator Model with Infection in Both Prey and Predator. Filomat, 29,1753-1767.

[16] Kar, T. and Ghosh, B. (2012) Sustainability and Optimal Control of an Exploited Prey-Predator System through Provision of Alternative Food to Predator. Biosystems, 109, 220-232.

https://doi.org/10.1016/j.biosystems.2012.02.003

[17] Hugo, A., Massawe, E.S. and Makinde, O.D. (2012) An Eco-Epidemiological Mathematical Model with Treatment and Disease Infection in both Prey and Predator Population. Journal of Ecology and the Natural Environment, 4, 266-279.

[18] Bodine, E., Gross, L. and Lenhart, S. (2008) Optimal Control Applied to a Model for Species Augmentation. Mathematical Bioscience and Engineering, 5, 669-680.

https://doi.org/10.3934/mbe.2008.5.669

[19] Goh, B., Leitman, G. and Vicent, T.L. (1974) Optimal Control of a Prey Predator System. Mathematical Bioscience, 19, 263-286.

https://doi.org/10.1016/0025-5564(74)90043-1

[20] Naji, R.K. and Dina, S.A.J. (2011) The Dynamics of Stage Structured Prey-Predator Model Involving Parasitic Infectious Disease. Applied Mathematics: An International Journal, 6, 529-551.

[21] Stoddart, A.W.J. (1967) Existence of Optimal Controls. Pacific Journal of Mathematics, 20, 167-177.

https://doi.org/10.2140/pjm.1967.20.167

[22] Fleming, W.H. and Rishel, R.W. (1975) Deterministic and Stochastic Optimal Control. Vol. 268, Springer-Verlag, New York.

https://doi.org/10.1007/978-1-4612-6380-7

[23] Lenhart, S. and Workman, J. (2007) Optimal Control Applied to Biological Models. Champman and Hall/CRC, London.

[24] Schaller, G.B. (1972) The Serengeti Lion: A Study of Predator-Prey Relations. University of Chicagopress, Chicago.

[25] Kar, T.K. (2010) A Dynamic Reaction Model of a Prey-Predator System with Stage-Structure for Predator. Modern Applied Science, 4, 183-195.

https://doi.org/10.5539/mas.v4n5p183