The Dynamics of Biological Models with Optimal Harvesting

This paper aims to introduce a concept of an equilibrium point of a dynamical system which will call it almost global asymptotically stable. A biological prey-predator model is also analyzed with a modification function growth in prey species. The conditions of the local stable and existence of all its equilibria are given. After that the model is extended to an optimal control problem to obtain an optimal harvesting strategy. The discrete time version of Pontryagin's maximum principle is applied to solve the optimality problem. The characterization of the optimal harvesting variable and the adjoint variables are derived. Finally numerical simulations of various set of values of parameters are provided to confirm the theoretical findings.


Introduction
The theory of mathematical models plays an important role for studying populations behavior. These models can be described in continuous time case or in discrete time case by a system of ordinary differential equations or a system of deference equations respectively, this description depends on the study problem . Discrete time systems are suitable for populations that reproduce at specific times each month or year or each circle, this can be seen in many insects populations, marine fish, and plants. There are variety of studies in the literature that analyzed and investigated the dynamical behavior of this kind of models, we refer to [1,2,3,4],and the references therein.
The essential concept in mathematical modeling is the stability of an equilibrium point so that an equilibrium point is called globally asymptotically stable if the solution approaches to this equilibrium regardless of the initial condition, while it is called locally asymptotically stable if there exists a neighborhood of this equilibrium such that from every initial condition within this neighborhood the solution approaches to it. Some authors in the literatures( [5,6,7,8]) proved that the unique positive equilibrium point in their models is global asymptotically stable even one or more than one boundary equilibrium point always exits. They ignored or excluded the boundary points from the domain of the function in their biological models.In these cases they violated the definition of the global asymptotically stable. To treat this situation we introduce a definition of an equilibrium point which we call it almost global asymptotically stable, simply means an equilibrium point is global asymptotically stable in the interior of the domain. A system of difference equations may have and show a rich and more complicated dynamical behaviors even for a simple one dimensional system. For example the logistic equation which is very known equation its equilibria can vary from stability behavior to chaotic behavior [9,10,11]. Many researchers investigated and analyzed different kind of two or more than two dimensional models in ecology [12,13,14], they derived conditions for local and global stability of solutions as well as the existence of periodic solutions [15,16]. Some authors have discussed and assumed that the life of populations have two stages immature and adults. However stage-structured or age structured models are are considered in the literatures [17,18,19]. In this study we will introduce a concept of an equilibrium point of a dynamical system as well as biological model, prey predator models with modification growth function in prey species is investigated in details. The general form is given by The continuous time version of model (1.1) can be written in the following form Where x t , (x(t)), y t (y(t)) and h(x t ) are the prey population density, the predator population density and the predator functional response at time t respectively. The parameters r, c, a,and d are model parameters supposing only positive values.These parameters are the intrinsic growth rate of prey species, the mortality rate of predators species, the maximum per capita killing rate,and the conversion rate predator respectively.While b, m, k and n are positive constants. This paper is organized as following: In section 2 we will discuss the model (1.1) with absent the predator species when n = m = b = 1. In section 3 the preypredator system is analyzed and all behavior of its equilibria are investigated.
In section 4 the model is extended to an optimal control problem.The discrete time version of Pontryagin's maximum principle is applied to solve the optimality problem. In section 5 numerical simulations is provided to confirm the theoretical results. Finally conclusion is given.

Single species model
Definition 1. Consider the following nonlinear discrete dynamical system x t+1 = f (x t ), where f : D ⊂ R n → R n . An equilibrium point x e , that means f (x e ) = x e is said to be almost global asymptotically stable in D if it is global asymptotically stable in D − ∂D. For the continuous time case the definition will be as follows: Consider the following nonlinear continuous dynamical system x (t) = f (x(t)), where f : D ⊂ R n → R n . An equilibrium point x e , that means f (x e ) = 0 is said to be almost global asymptotically stable in D if it is global asymptotically stable in D − ∂D.
Remark 1. 1-It is clear that every global asymptotically stable is almost global asymptotically stable, and the converse is not true. For this one can see in ( [7] authors proved that the positive equilibrium point S 3 = (x * 1 , x * 2 ) there, is globally asymptotically stable only in the interior of the domain which means it is almost global asymptotically stable. However it is not global asymptotically stable because the trivial equilibrium point always exists in their model. 2-Clearly that if x e is almost global asymptotically stable then it is local asymptotically stable,but the converse is not true.For example : Consider the the following system . If s = 3.1 then x * = 0.558 and y * = 0.7646. The point x * = 0.558 is locally stable which is not almost global asymptotically stable point. Now we will investigate the dynamics of single species model of (1.1) in the absent of the predator species and n = m = b = 1. Thus the model will be as the following The model(2.1) has two equilibria , namely, the trivial equilibrium point x 1 = 0, and the unique positive point x 2 = ln( r−1 k ). The trivial equilibrium point always exists, while the positive equilibrium exists when r > k + 1. The following lemma gives the behavior of its equilibria.

The equilibrium point
Proof. It is clear that f (x 1 ) = r k+1 , so that the results in 1 can be easily obtained.
r−1 < k < r − 1 and the results can be got. Now we will consider a situations that population is exposition to harvest by a constant rate harvesting which is proportional to the respective population size therefore the model (2.1) including the harvesting will be as the following : Where h is a positive constant representing the intensity of removing due to hunting or removal. It is obvious that one cannot remove more than the population density therefore h ≤ h max < 1, h max is the maximum removing amount. The model(2.2) has also two equilibria , the trivial equilibrium point x o ,which always exists,and the unique positive equilibrium x h = ln( r−(1+h) k(1+h) ) exists only when r−(1+h) k(1+h) > 1. Next lemma describes the behavior of the equilibria of model(2.2).
, therefore all results can be obtained.

Two Species Model,Prey-Predator Model
In this section we will study in details the dynamics of the two species model discrete time case of model (1.1) with n = m = b = 1. Thus the system can be written as The all parameter a, r, c, d and k are defined the same as before. By solving the following algebraic equation one can get all equilibrium points of the model(3.1): Therefore we have the following lemma. 3. The unique positive equilibrium point e 2 = (x * , y * ) = ( 1+c d , r−(1+ke x * ) a(1+ke x * ) ),which exists if r > 1 + ke x * .
In order to investigate the dynamic behavior of the model(3.1) one has to compute the general Jacobian matrix of the model (3.1) at point (x,y). This is given by: − ay, J 12 = −ax, J 21 = dy, and J 22 = dx − c.
Next theorems give the local stability of e 0 , and e 1 respectively. where I 1 = ((r − 1)e − 2r r−1 , (r − 1)), Proof. The Jacobian matrix at the point e 1 is In order to discuss the dynamic behavior of the unique positive equilibrium, the next lemma is needed. Proof. see [19]. Where Proof. The Jacobian matrix at the unique positive equilibrium point is given by So that the characteristic polynomial of J e 3 is with d > N . According to lemma (4) the proof is finished.

An optimal harvesting approach
In this section we will extend the model(1.1) to an optimal control problem and will discuss the optimal harvesting management of renewable resources.We assume that the population is harvested or removed with the harvesting rate h t , which represents our control variable. For the single species the model (2.1) including the harvesting effect becomes : The x t , r, and k are defined as before. In this problem the free terminal value problem is discussed and the terminal time T is specified. The aim is to maximize the following objective functional where c 1 h t x t represents the amount of money that one has to obtain, and c 2 h 2 t is the cost of catching and supporting the animal. c 1 and c 2 are positive constants. The control variable is subject to the constraint Now according to the discrete version of Pontryagin's maximum principle [20], the Hamiltonian functional for this problem is given by − h t , is the adjoint variable or shadow price [21]. Then the characterization of the optimal control solution is In order to extend the two species model to an optimal control problem, the model (3.1)with control harvesting variable will become as the following: Our aim in this problem is to get an optimal harvesting amount for that we will maximize the following cost functional All terms and parameters are as before.So that the Hamiltonian functional will be as the following: Where λ 1 and λ 2 are the adjoint functions that satisfy : Furthermore the characterization of the optimal harvesting solution h * satisfies: An iterative method in [20] is used to get the optimal control with corresponding optimal state solutions of the above optimal control problems at time t by maximizing the Hamiltonian functional at that t numerically.

Numerical Simulations
In this section we will illustrate the theoretical findings numerically for various set of parameters. For the optimal control problem an iterative numerical method in [20] is used to determine the optimal solutions with corresponding state solutions. For the control problem of single species we choose this set of values of parameters r = 1.999; k = .8; c1 = 0.1; c2 = 0.01, x 0 = 0.1 and T = 80 so that the total optimal objective functional Jopt is found equal to 0.0412. In Figure 4 the prey population density with control , without control and with constant harvesting is plotted. Figure 5 shows the control variable as a function of time. For the optimal control problem of two species model we choose this set of values of parameters a = 0.1; k = 2.1; c = 0.5; d = 2.9; r = 5.2; c1 = 0.025; c2 = 0.08, (x 0 , y 0 ) = (0.5, 0.8) and T = 80. So that the total optimal objective functional Jopt is found equal to 0.04121. Figures 6-7 shows the prey population and the predator population with control ,without control and with constant control respectively, and Figure 8 shows the control variable as a function of time. Finally Table 1 contains the total optimal objective functional and other different total harvesting amount strategies of both control problems by using the same values of the parameters in each problem.

Conclusion
In this paper the definition of almost global asymptotically stable of an equilibrium point is introduced with examples. We have also investigate biological models in discrete time case for one species with a modification growth function and for two species, prey-predator, model.We have studied the dynamics behavior of the equilibrium points for each model. after that we have extended these model to an optimal control problems. We have used the Pontaygin's principle maximum to get the optimal solutions numerically. Finally one can investigate the dynamics behavior of the continuous time case as well as the other values of n, m, and b in future work.