\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2012 (2012), No. 76, pp. 1--11.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2012 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2012/76\hfil Bioeconomical Ricker's model] {Bioeconomical Ricker's model of marine protected areas} \author[C. N. Christou, L. V. Idels \hfil EJDE-2012/76\hfilneg] {Cameron N. Christou, Lev V. Idels} % in alphabetical order \address{Cameron N. Christou \newline Department of Mathematics, University of British Columbia, Vancouver, B.C., V6T 1Z2, Canada} \email{cnchrist@math.ubc.ca} \address{Lev V. Idels \newline Department of Mathematics, Vancouver Island University\\ 900 Fifth St. Nanaimo BC, V9S5S5, Canada} \email{lev.idels@viu.ca} \thanks{Submitted November 29, 2011. Published May 14, 2012.} \subjclass[2000]{34C60, 37N25, 49K15, 92B05} \keywords{Marine protected areas; source-sink models; fishery models, \hfill\break\indent bioeconomical analysis, Pontryagin's maximum principle, nonlinear differential equations} \begin{abstract} Marine protected areas (MPA) become part of modern fishery management to safeguard marine life and sustain ecosystem processes. Based on a classical Ricker's model, the deterministic nonlinear sink-source model of MPA is presented. Sufficient conditions for the existence of positive bounded steady-states are obtained. The present value of net revenue is maximized subject to population dynamics in the fishing zone and in the protected area. The analysis has shown that there is an optimal equilibrium solution, and the best harvesting policy to attain this equilibrium position is a bang-bang control effort. It was demonstrated numerically by comparing the optimal harvesting rate against a constant harvesting rate, and the fast convergence to the optimal equilibrium guarantees greater profits under the optimal harvesting strategy. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}[theorem]{Remark} \newtheorem{example}[theorem]{Example} \allowdisplaybreaks \section{Introduction} Traditional methods of controlling fishing include: \begin{itemize} \item limiting the effort (restricting the fishing season); \item limiting the catch (restricting the number or biomass of fish captured); \item limiting the fish size (restricting the size of fish captured); \item limiting the location (restricting the places available for fishing). \end{itemize} The last option includes the establishment of marine reserves, areas permanently closed to all fishing. These can supplement or potentially replace other methods of regulation. For this reason, the reserve concept has recently attracted great interest in the community of fisheries scientists and managers. Marine protected areas have been promoted as conservation and fishery management tools, and at present, there are over 1300 marine reserves in the world \cite{Ami,Gell,Graf,Ko,Pi,So, Sum}. Since 1999, the number and extent of, and the amount of research on, MPA has grown rapidly. Areas protected from fishing provide multiple benefits, including: reduced likelihood of a stock collapsing, enhanced spawning biomass, improved recruitment survival of fish to a more mature state, allowing increased harvest and recovery of the habitat. Much of the empirical case study research on the ecological impact of MPA indicate that within the no-take area boundaries, fish both grow and attain a broader age and size distribution \cite{Ko,Pi,So}. There is robust demonstration of conservation benefits of the MPA, but fishery benefits and design of reserves remain controversial. Most of the researchers suggest that in the long term no-fishing zones sustain both fisheries and fish populations. It is also known that highly migratory species are among the most overexploited finfish species. Many species, especially adults such as cod, show directed migration on an annual basis; thus adults could be very vulnerable at certain periods of the year when most of them would be likely to move outside of a protected area \cite{Do2,Gell,Graf}. The biological interconnections between patches may be of many varieties and we will use a sink-source model based on the Ricker's growth curve. In recent publications \cite{Arm,Do,Dub,Kara}, MPAs are modeled by using logistic differential equations. More complex models were studied in \cite{B1,B2,B3,Je,I}. Suppose that two regions $R$ (reserve area) and $U$ (unprotected or fishing zone) have habitat areas: $A_1$ is the area of a protected zone and $A_2$ is the area of the fishing zone. Let $x(t)$ denotes the biomass of fish in the reserve $R$, and $y(t)$ is the fish biomass in the unprotected (fishing) area. Assume that fishing takes place only in the region $U$, with the region $R$ established as a marine reserve or protected area. The two equations in model \eqref{rickermodel} are linked by the movement terms, which depend on a difference of biomass densities per unit of the habitat areas. The depletion of fish in region $U$ will generally decrease the fish density there and cause fish to move away from a higher density in the protected region $R$. We introduced the migration rate \[ d_i = \frac{\sigma}{A_i},\quad \sigma > 0, \] where $\sigma$ is a migration coefficient. Then the coupled system of nonlinear differential equations \begin{equation} \label{e1} \begin{gathered} \frac{dx}{dt} = -m_1 x(t) + d_2 y(t) - d_{1} x(t) + F_1 (x(t)), \\ \frac{dy}{dt} = -m_2 y(t) + d_{1} x(t)- d_2 y(t) + F_2 (y(t)) - E(t)qy(t) \end{gathered} \end{equation} could represent the population dynamics in each region, where $m_i$ is a natural mortality rate, \[ F_i (x) = p_i xe^{-x},\quad p_i > 0. \] is the birth function. Finally, the system under study has the form \begin{equation}\label{rickermodel} \begin{gathered} \frac{dx}{dt} = -m_1 x + d_2 y - d_{1} x + p_{1}x e^{-x}, \\ \frac{dy}{dt} = -m_2 y + d_{1} x- d_2 y + p_2ye^{-y} - E(t)qy. \end{gathered} \end{equation} In system \eqref{rickermodel} a positive constant $q$ is called a catchability coefficient, whereas the effort function $E(t)$ is a dynamic variable. MPA design is a joint economic and biological problem, and MPA creation alters economic incentives. In fisheries the optimal control has often been applied to compute bieconomic optima via locally-optimal linear feedback controls, and mostly the optimal solutions obtained are bang-bang controls \cite{Arm,Bea,conrad,Do,Dub,Je,Mch,Ne}. \section{Equilibrium analysis} To find all equilibria of system \eqref{rickermodel} we set \begin{equation}\label{eqsystem} \begin{gathered} m_1 x + d_{1} x - d_2 y - p_1 xe^{- x} = 0, \\ m_2 y - d_{1} x + d_2 y + Eqy - p_2 ye^{- y} = 0. \end{gathered} \end{equation} Apart of a trivial solution, system \eqref{eqsystem} has a positive equilibrium $(x^{*}, y^{*}).$ From system \eqref{eqsystem} we have two curves: \begin{equation}\label{equilib2} \begin{gathered} L_{1}: \quad y =\phi_{1}(x)= \frac{m_{1}+d_1-p_{1}e^{-x} }{d_2} x,\\ L_2: \quad x=\phi_2(y)= \frac{m_2+d_2+Eq-p_2e^{-y}}{d_1} y. \end{gathered} \end{equation} \begin{theorem}\label {thm3.1} System \eqref{rickermodel} has a unique internal positive equilibrium if the following two conditions hold: (i) $m_1+d_{1}\leq p_1$ and (ii) $m_2+d_2+Eq\leq p_2$. \end{theorem} \begin{proof} To prove the existence of a nontrivial equilibrium, firstly note that curves $L_1$ and $L_2$ have asymptotes $$ y=\frac{m_1 +d_1}{d_2}x\quad\text{and}\quad y=\frac{d_1}{m_2+d_2+Eq}x $$ correspondingly. Clearly, $$ \frac{m_{1}+d_1}{d_2}>\frac{d_1}{m_2+d_2+Eq} \hspace{1mm}, $$ thus, for sufficiently large $x$, points on the curve $L_1$ lie above the corresponding points of the curve $L_2$. On the other hand, in the neighborhood of the origin, each of the conditions (i)--(ii) guarantees that points on the curve $L_2$ lie above the points of the curve $L_1$. Therefore a positive internal equilibrium of system \eqref{rickermodel} exists. To prove that this equilibrium is a unique point, firstly, we note that from system \eqref{equilib2} \begin{equation}\label{divequilib2} \begin{gathered} L_1: \quad \frac{dy}{dx} =\frac{y}{x}+ \frac{p_1}{d_2} x \exp(- x )>\frac{y}{x},\\ L_2:\quad \frac{dx}{dy} =\frac{x}{y}+\frac{p_2}{d_1} y \exp(- y )>\frac{x}{y}, \end{gathered} \end{equation} since $x > 0$ on $L_1$ and $ y >0$ on $L_2$. Let $\theta$ be a polar angle of the point on the curve $L_1$ (with $Ox$ as a polar axis), $$ \theta=\arctan\frac{y}{x}=\arctan\frac{\phi_1}{x}. $$ If $x$ moves from $0$ to $\infty$, then $$ \frac{d\theta}{dx}=\frac{1}{1+(\frac{y}{x})^2}\frac{d}{dx} \frac{\phi_1}{x}=\frac{x}{x^2+y^2} \big(\frac{d\phi_1}{dx}-\frac{\phi_1}{x}\big). $$ The latter equality and inequalities \eqref{divequilib2} guarantee $\frac{d\theta}{dx}>0$; thus $\theta(x)$ is a monotone increasing function with $$ \theta_0 < \theta(x) < \arctan\frac{m_1+d_1}{d_2}. $$ Similarly, let $\vartheta$ be a polar angle of the point on the curve $L_2$ with the same polar axis. If $y$ moves from $0$ to $\infty$, then $\vartheta (y)$ is a monotone decreasing function, with $$ \arctan\frac{d_1}{m_2+d_2+Eq}< \vartheta(y) < \vartheta_0. $$ At the equilibrium point $\theta=\vartheta$, increase of the function $\theta$ and decrease of the function $\vartheta$ guarantee the uniqueness of a nontrivial equilibrium. Let none of the conditions (i)-(ii) hold; then for any $x>0$ and $y>0$, \begin{equation}\label{impossib} \begin{gathered} \arctan\frac{m_1+d_1}{d_2}>\theta \geq \theta_0=\arctan\frac{m_1+d_1-p_1}{d_2}>\\ \vartheta_0=\arctan\frac{d_1}{m_2+d_2+Eq-p_2}\geq \vartheta > \arctan\frac{d_1}{m_2+d_2+Eq}\nonumber. \end{gathered} \end{equation} Therefore, the equality $\theta=\vartheta$ is impossible, and the positive equilibrium does not exist. Theorem \ref{thm3.1} is proved. \end{proof} \begin{remark} \label{rmk2.1} \rm Note that nonlinear system \eqref{eqsystem} can not be solved analytically, the latter fact makes bioeconomical analysis of the model more difficult than the corresponding analysis of the system with logistic growth \begin{gather*} m_1 x + d_{1} x - d_2 y - p_1 x^{2} = 0, \\ m_2 y - d_{1} x + d_2 y + Eqy - p_2 y^{2} = 0. \end{gather*} in \cite{Do,Dub,Graf,Kara}. \end{remark} \begin{theorem} \label{thm2.2} Any solution of system \eqref{rickermodel} with nonnegative initial functions and positive initial conditions ultimately enters the square region $$ K=\{0< x(t) \leq a \text{ and } 0 0$ and $y(t)>0$ for $t>0$, provided that $x(0)> 0$ and $y(0)> 0$. Let $ z=x+y $, then from system \eqref{rickermodel} we have \[ \frac{dz}{dt} =-m_1x-m_2y-Eqy + F_1(x)+F_2(y). \] For $0< \epsilon < \min (m_1, m_2 )$ we have \begin{align*} \frac{dz}{dt}+\epsilon z &= -(m_1-\epsilon)x-(m_2-\epsilon)y-Eqy+ F_1(x))+F_2(y) \\ &< F_1(x)+ F_2(y) \leq \max\{F_1,F_2\}\\ &\leq \max \{p_1, p_2\}=W. \end{align*} Therefore, $\frac{dz}{dt}+\epsilon z < W $, or \begin{equation}\label{eval} z(t)< \frac{W}{\epsilon}+ \big(z(0)+\frac{W}{\epsilon}\big) e^{-\epsilon t}. \end{equation} Thus all solutions of \eqref{rickermodel} are bounded. \end{proof} \section{Optimal harvesting} In the presence of an MPA, the objective is to optimally exploit the available resources. The discounted revenue of the fishing industry is given in \cite{Kot} \[ J = \int_0^{\infty}e^{-\delta t}\left [pqy(t) - c\right ]E(t) \mathrm{d}t, \] where $\delta$ is the instantaneous rate of annual discount, $p$ is the unit price of the biomass being harvested, and $c$ is the cost per unit of harvested biomass. The cost of fishing is assumed to decrease with the size of the stock, i.e. $c'=\frac{dc}{dy} \leq 0$. Thus, the objective is to maximize the revenue functional $J$ with respect to the harvesting rate $E(t)$. In order to prevent extinction of the species, for the given size of the MPA, determined by the parameters $A_1$ and $A_2$, a critical harvesting rate $E_c$ is defined as \[ E_c = \sup\{ E : m_1 + d_{1} \leq p_1 \text{ or} m_2 + d_2 + Eq \leq p_2 \} \] or \[ E_c = \sup\{ E:(m_1 + d_{1} - p_1)(m_2 + d_2 + Eq - p_2) < d_{1}d_2\} \] Therefore, conditions (i)--(ii) in Theorem \ref{thm3.1} are satisfied. Note that if condition (i) of Theorem \ref{thm3.1} is satisfied, then $E_c = \infty$ since any fishing effort will lead to non-trivial solutions. Also note that $E_c > 0$, otherwise the fish population will become extinct by natural causes. Therefore, a maximal harvesting rate is chosen so that $0 \leq E \leq E_{\rm max} < E_c$. In this manner, the fish population should be protected from over-exploitation. Hence the set of allowable harvesting rate functions, $E(t)$, belongs to the set \[ E = \{ E(t) : 0 \leq E(t) \leq E_{\rm max} \}. \] The objective is to find $\sup_{E} J$ subject to the dynamical system \eqref{rickermodel}. An application of the Pontryagin's Maximum Principle yields the Hamiltonian \begin{equation}\label{ham1} \begin{aligned} H &= e^{-\delta t}(pqy - c)E(t) + \lambda_1\left(-\left(m_1+d_{1}\right)x + d_2y + p_1xe^{-x}\right) \\ &\quad + \lambda_2\left(-\left(m_2 + d_2+Eq\right)y + d_{1}x + p_2ye^{-y} \right). \end{aligned} \end{equation} The PMP also produces the equations of motion for the shadow prices \begin{equation}\label{costate1} \begin{gathered} -\dot{\lambda}_1 = \frac{\partial H}{\partial x} = \lambda_1\left[-\left(m_1+d_{1}\right)+p_1(1-x)e^{-x}\right]+ \lambda_2d_{1} \\ -\dot{\lambda}_2=\frac{\partial H}{\partial y} = e^{-\delta t}\left(pq - c'\right)E(t) + \lambda_1d_2 + \lambda_2\left[-\left(m_2 + d_2+Eq\right)+p_2(1-y)e^{-y}\right].\nonumber \end{gathered} \end{equation} Consider the change of variables: $\lambda_1 = e^{-\delta t}\Gamma_1$, $\lambda_2=e^{-\delta t}\Gamma_2$. Substitution of $\Gamma_1$ and $\Gamma_2$ in \eqref{ham1} yields \begin{equation}\label{ham2} \begin{aligned} H &= e^{-\delta t}\left[(pqy - c)E(t) + \Gamma_1\left(-\left(m_1+d_{1}\right)x + d_2y + p_1xe^{-x}\right)\right] \\ &\quad + e^{-\delta t}\left[ \Gamma_2\left(-\left(m_2 + d_2+Eq\right)y + d_{1}x + p_2ye^{-y} \right)\right]. \end{aligned} \end{equation} Equivalently, substitution of $\Gamma_1$ and $\Gamma_2$ in \eqref{costate1} produces \begin{equation}\label{CO1} \begin{gathered} -e^{-\delta t}\dot{\Gamma}_1 + \delta e^{-\delta t}\Gamma_1 = e^{-\delta t}\Gamma_1\left[-\left(m_1+d_{1}\right)+p_1(1-x)e^{-x}\right] + e^{-\delta t}\Gamma_2d_{1}\\ \Rightarrow \; \dot{\Gamma}_1 = \Gamma_1\left[\delta +m_1+d_{1}-p_1(1-x)e^{-x}\right] - \Gamma_2d_{1}, \end{gathered} \end{equation} and \begin{equation} \label{CO2} \begin{gathered} \begin{aligned} -e^{-\delta t}\dot{\Gamma}_2 + \delta e^{-\delta t}\Gamma_2 &= e^{-\delta t}\left(pq - c'\right)E(t) + e^{-\delta t}\Gamma_1d_2 \\ + e^{-\delta t}\Gamma_2 \left[-\left(m_2 + d_2\right)+p_2(1-y)e^{-y}-Eq\right] \end{aligned} \\ \Rightarrow\; \dot{\Gamma}_2 = (c'-pq)E - \Gamma_1d_2 + \Gamma_2[\delta +m_2 + d_2- p_2(1-y)e^{-y}+Eq]. \end{gathered} \end{equation} Equations \eqref{CO1} and \eqref{CO2} represent the weighted shadow prices \cite{Kot}. The objective now is to maximize the Hamiltonian \eqref{ham2} over the set $E$ subject to equations \eqref{CO1} and \eqref{CO2}. Note that the Hamiltonian is linearly dependent on $E$; therefore, the optimal harvesting rate must be a bang-bang control in the set $E$ . Explicitly, the optimal harvesting rate is \[ E_{\rm opt}(t) = \begin{cases} E_{\rm max} & T> 0\\ 0 & T < 0 \\ \tilde{E}& T=0 , \end{cases} \] where $\tilde{E}$ is a singular harvesting rate and $T=pqy - c - qy\Gamma_2$. This singular harvesting rate occurs because the value of $E(t)$ is not determined when $T = 0$. The possible end-states of the system \eqref{rickermodel} under the optimal harvesting rate $E_{\rm opt}(t)$ are the same as those mentioned in \cite{Kara}. These states are \begin{itemize} \item No harvesting \item Maximum harvesting \item Singular rate harvesting \item A ``bang-bang'' cycle of harvesting rates. \end{itemize} If ``no harvesting'' is the optimal end-state of the system, then it must be the case that it is never profitable to fish for any biomass up to the natural carrying capacity of the unprotected fishing zone, i.e. $pqy - c \leq 0$ for all positive values of $y$ up to the carrying capacity of the zone. If this was not the case, then there would exist a bang-bang cycle that would produce a greater profit. If ``maximum harvesting'' is the optimal end-state of the system; then it must be the case that fishing at a constant rate of $E(t) = E_{\rm max}$ drives the system \eqref{rickermodel} to a non-trivial equilibrium point $(x^*, y^*)$, and that $pqy^* - c(y^*) \geq 0$. If this was not the case, then the profit would be negative, and greater profit would be obtained by not harvesting at all. For the singular harvesting rate i.e. $T=0$, function $\Gamma_2$ satisfies \begin{equation}\label{SING} \Gamma_2 = p - \frac{c}{qy}. \end{equation} If the end-state is a singular rate, then equation \eqref{SING} holds for an extended period of time, and \begin{equation}\label{singdot} \dot{\Gamma}_2 = \Big(\frac{c'}{qy} - \frac{c}{qy^2}\Big)\dot{y}. \end{equation} Equating \eqref{CO2} and \eqref{singdot} produce an expression for function $\Gamma_1$ which can be solved explicitly. \begin{equation}\label{singdot2} \Gamma_1 = \frac{1}{f_{y}}\Big(\delta\big(p-\frac{c}{qy}\big) - (pq-c')E - \big(p-\frac{c}{qy}\big){g_{y}} + \frac{c'g}{qy} - \frac{cg}{qy^2}\Big), \end{equation} where \[ f = -m_1 x + d_2 y - d_{1} x + p_1xe^{-x},\quad g = -m_2 y - d_2 y +d_{1} x + p_2ye^{-y} - Eqy. \] Since the system stays in this singular state, function $\Gamma_1$ also satisfies \begin{equation} \label{gam1} \begin{split} \dot{\Gamma}_1 & = -\frac{f_{yy} \dot{y}}{f^{2}_{y}} \Big(\delta\big(p-\frac{c}{qy}\big) - (pq-c')E - \big(p-\frac{c}{qy}\big)g_{y} + \frac{c'g}{qy} - \frac{cg}{qy^2}\Big)\\ &\quad +\frac{1}{f_{y}}\Big(\frac{-\delta c'\dot{y}}{qy} + \frac{\delta c \dot{y}}{qy^2} + c''E\dot{y}-g_{yy}\dot{y}\big(p-\frac{c}{qy}\big)\\ &\quad +2\frac{c'\dot{y}g_{y}}{qy} - 2\frac{c\dot{y}g_{y}}{qy^2} +\frac{c''\dot{y}g}{qy} - 2\frac{c'g\dot{y}}{qy^2} + 2\frac{cg\dot{y}}{qy^3} \Big) \end{split} \end{equation} Equating \eqref{CO1} and \ref{gam1} implicitly define a singular harvesting rate for every point in the first quadrant of the $(x,y)$-plane. However, due to the imposed restriction $0 \leq E(t) \leq E_{\rm max}$, it may not be possible to employ a singular harvesting rate to certain points in the plane. If system \eqref{rickermodel} reaches an equilibrium, while fishing at a singular harvesting rate, $(x^*,y^*)$; then the harvesting rate is defined by \begin{equation}\label{equiharvest} \tilde{E} = \frac{1}{q}\Big( -m_2 - d_2 + d_{1}\frac{x^*}{y^*} + p_2e^{-y^*}\Big). \end{equation} \section{Optimal harvesting: equilibrium solutions}\label{optiharvest} The objective of this section is to optimize the profit earned while keeping the total biomass of system \eqref{rickermodel} at equilibrium. If system \eqref{rickermodel} is in stable equilibrium, then the harvesting rate must be constant and equal to $\tilde{E}$ defined by \eqref{equiharvest}. \textbf{Part 1.} For a given harvesting rate to be optimal, in the sense that profit is maximized, the following condition must be satisfied: \[ \frac{\partial H}{\partial E} = 0, \] where $H$ is given by equation \eqref{ham2}. Therefore, equation \[ \frac{\partial H}{\partial E} = pqy - c - \Gamma_2 qy = 0 \] must be satisfied and \begin{equation}\label{gamma2} \Gamma_2 = p - \frac{c}{qy}. \end{equation} Under an equilibrium solution, it is also the case that $\dot{\Gamma}_1 = \dot{\Gamma}_2 = 0$; i.e., \begin{gather}\label{gamma1} \Gamma_2d_{1}-\Gamma_1 [\delta +m_1+d_{1}-p_1(1-x)e^{-x} ]=0,\\ \label{gamma3} (pq - c')E(t) + \Gamma_1d_2 - \Gamma_2 [\delta +m_2 + d_2-p_2(1-y)e^{-y}+Eq]=0. \end{gather} Substitution of equations \eqref{singdot2}, \eqref{equiharvest} and \eqref{gamma2} in \eqref{gamma1} and \eqref{gamma3} define a system of equations that can be solved for the optimal equilibrium solution, $(x^*, y^*)$. Alternatively, the optimal equilibrium solution can be found by locating the critical points of the profit function \[ (pqy - c)\tilde{E}, \] where $\tilde{E}$ is given by equations \eqref{equiharvest}, and \[ y = \frac{m_1x + d_{1}x - p_1xe^{-x}}{d_2}. \] The latter is a necessary condition for system \eqref{rickermodel} to be at the steady state. Clearly, either of these two methods will produce the same result. \textbf{Part 2.} The objective now is to reach the optimal equilibrium state in an optimal way. As shown in the Section 3, this objective can be achieved by using a bang-bang control policy for the harvesting rate $E(t)$. Until the system reaches the optimal equilibrium point, we define \[ E(t) = \begin{cases} E_{\rm max}, &\text{if } pqy(t) - cy(t) - \Gamma_2qy(t) > 0\\ 0, & \text{otherwise } \end{cases} \] and then $E(t)$ is a constant at the equilibrium harvesting rate. Note that \begin{equation} \Gamma_2 = p - \frac{c(y^*)}{qy^*},\nonumber \end{equation} where $(x^*, y^*)$ is the optimal equilibrium point. Alternatively, since an internal equilibrium point of system \eqref{rickermodel} is globally asymptotically stable under constant harvesting rates, the optimal equilibrium point can be reached by employing the constant optimal harvesting rate for all time. \section{Numerical simulations} The following simulations were produced using the system parameter values: $m_1 = 0.5$, $m_2 = 0.2$, $\sigma = 0.4$, $A_1 = 0.2$, $A_2=0.8$, $p_1 = p_2 = 1$, $q = 0.1$, $p=7$, $E_{\rm max} = 7$, and $c = 0.4 + \frac{1}{1+y}$. \begin{example} \label{examp5.1} \rm Using the methods described in Section \ref{optiharvest}, the optimal equilibrium solution is found to be $(x^*, y^*) = (0.4044, 1.4822)$. Keeping for all time the system \eqref{rickermodel} at this equilibrium position produces a maximum profit of 0.17085. This information is presented in Figure \ref{optequibfig}. \begin{figure}[ht] \begin{center} \includegraphics[trim= 2.5cm 7.5cm 3cm 7.5cm, clip, scale=0.5]{fig1} %OptEquib \end{center} \caption{Finding the optimal equilibrium point for maximizing profit}\label{optequibfig} \end{figure} Starting from the initial conditions $(x_0, y_0) = (0,2)$, the trajectory of the system is shown in Figure \ref{trajectoryfig}. \begin{figure}[ht] \begin{center} \includegraphics[trim= 2.5cm 7.5cm 3cm 7.5cm, clip, scale=0.5]{fig2} % trajectory.eps \end{center} \caption{Trajectory of the system given the initial condition $(0, 2)$}\label{trajectoryfig} \end{figure} The full phase portrait of this system is shown in Figure \eqref{phasefig}. \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig3} %optim \end{center} \caption{Phase portrait of the system}\label{phasefig} \end{figure} The switching times of the function and the value of the harvesting rate become apparent in view of the phase portrait. Also, the time where the singular solution is the optimal harvesting rate is also visible. \end{example} \begin{example} \label{examp5.2} \rm As a second example of the optimal harvesting strategy, consider the trajectory followed by the system \eqref{rickermodel} from the initial position $(x_0, y_0) = (2, 1.2)$. The trajectory is shown in Figure \eqref{trajec2fig}. The benefits of the optimal harvesting strategy are clearly demonstrated in this example. \begin{figure}[ht] \begin{center} \includegraphics[trim= 2.5cm 7.5cm 3cm 7.5cm, clip, scale=0.5]{fig4} %trajectory2.eps \end{center} \caption{Trajectory of the system given the initial condition $(2,1.2)$}\label{trajec2fig} \end{figure} \end{example} \section{Conclusions} It is difficult to solve the marine protected areas problems empirically, because they require population data of large spatial and temporal extent. One of the benefits of a theoretical modeling of MPAs is insight into the type of informative fisheries data that should be collected in order that the best design can be established. Based on a classical Ricker's model, the deterministic nonlinear sink-source model of MPAs is presented. Sufficient conditions for the existence of the positive bounded steady-states are obtained. The present value of net revenue is maximized subject to population dynamics in the fishing zones and in the protected areas. The impact of MPA can be felt in the critical value of the harvesting rate. For most parameter choices, it is unreasonable to create an MPA big enough so that condition (i) of Theorem \ref{thm3.1} is satisfied, and hence a critical value of the harvesting rate $E$ must not be exceeded in order to ensure that the fish species does not go extinct. The analysis has shown that there is an optimal equilibrium solution, and the best harvesting policy to attain this equilibrium position is a bang-bang control effort. This was demonstrated numerically by comparing the optimal harvesting rate against constant harvesting rates. The convergence to the optimal equilibrium solution was always faster, which implies greater profits under the optimal harvesting strategy. Also included is a phase portrait, which shows a more complete picture of the system being studied. We have attempted to illustrate the effect of the optimal harvesting strategy more clearly by picking parameters and initial conditions that better demonstrate the benefit of the optimal harvesting policy. Unfortunately, field data required to adequately populate mathematical models is not classified and scattered in various scientific papers and reports. Access to the specific fishery data will help to develop models that can predict the features required for optimizing the economic benefits of MPA. Contrariwise, one of the benefits of a theoretical modeling of MPA is insight into the type of data that should be collected in order that the best harvest strategies can be established in the fishing zones. There are many assumptions and simplifications in our model that are open to objections. \begin{thebibliography}{99} \bibitem{Ami} D. Ami, P. Cartigny, A. Rapaport; \emph{Can marine protected areas enhance both economical and biological situations?} C.R. Biologies \textbf{328} (2005) 357-366. \bibitem{Arm} C. Armstrong; \emph{A note on the ecological-economic modelling of marine reserves in fisheries}, Ecol. Economics \textbf{62} (2007) 242-250. \bibitem{Bea} A. Beattie, U. Sumaila, V. Christensen, D. Pauly; \emph{A model for the bionomic evaluation of MPA size and placement in the Northern Sea}, Nat. Resource Mod. \textbf{15} (2002) 413-37. \bibitem {Bell} R. Bellman, K. L. Cooke; \emph{Differential-Difference Equations}. Mathematics in Science and Engineering, Volume 28, 1963. \bibitem{B1} L. Berezansky, L. Idels, M. Kipnis; \emph{Mathematical Model of Marine Protected Area}, IMA Journal of Applied Math, 76 (2) 2011 312-325. \bibitem{B2} L. Berezansky, L. Idels, L. Troib; \emph{Global Dynamics of One Class of Nonlinear Nonautonomous Systems with Time-Varying Delays}, Nonlinear Analysis: Theory, Methods \& Applications, Volume 74, Issue 18, December 2011, Pages 7499-7512. \bibitem{B3} L. Berezansky, L. Idels, L. Troib; \emph{Global Dynamics of Nicholson-Type Delay Systems with Applications}, Nonlinear Analysis Series B: Real World Applications, 12 (2011) 436-445. \bibitem{conrad} J. Conrad; \emph{The bioeconomics of marine sanctuaries}, J. of Bioeconomics \textbf{1} (1999), 205-217. \bibitem{Do} L. Doyen, C. Bene; \emph{Sustainability of fisheries through marine reserves: a robust modeling analysis}, J. of Env. Management \textbf{69} (2003) 1-13. \bibitem{Do2} L. Doyen, M. Lara, J. Ferraris, D. Pelletier, Sustainability of exploited marine ecosystems through protected areas: A viability model and a coral reef case study Ecol. Mod. \textbf{208} (2007) 353-366. \bibitem{Dub}B. Dubey, P. Chandra, P. Sinha; \emph{A model for fishery resource with reserve area}, Nonlinear Analysis: Real World Appl. \textbf{4} (2003) 625-637. \bibitem{Gell} F. Gell, C. Roberts; \emph{Benefits beyond boundaries: the fishery effects of marine reserves}, Trends in Ecology and Evolution \textbf{18} (2003) 448-455. \bibitem{Graf} R. Grafton, T. Kompas, D. Lindenmayer; \emph{Marine reserves with ecological uncertainty}, Bull. of Math. Bio. \textbf{67} (2005) 957-971. \bibitem{Je} M. Jerry, N. Raissi; \emph{Optimal strategy for structured model of fishing population}, C.R. Biologies \textbf{328} (2005) 351-356. \bibitem{I} L. Idels, M. Kipnis; \emph{Stability Criteria for a Nonautonomous Nonlinear System with Delay}, Applied Mathematical Modelling, 33 (5) 2293-2297. \bibitem{Ko} J. Lubchenco, S. Palumbi, S. Gaines, S. Andelman; \emph{Plugging a hole in the ocean: the emerging science of marine reserves}, Ecol. Appl. \textbf{13} (2003) 3-7. \bibitem{Kara} T. Kar, H. Masuda; \emph{A Bioeconomical model of a single-species fishery with a marine reserve}, J. of Environ. Management \textbf{86 } (2008) 171-180. \bibitem{Kot} M. Kot; \emph{Elements of Mathematical Ecology}, Cambridge Univ Press, Cambridge, 2001. \bibitem{Mch} R. Mchich, P. Auger, N. Raissi; \emph{The stabilizability of a controlled system describing the dynamics of a fishery}, C. R. Biologies \textbf{328} (2005) 337-350. \bibitem{Ne} G. Neubert; \emph{Marine reserves and optimal harvesting}, Ecol. Lett. \textbf{6} (2003) 843-849. \bibitem{Pi} J. Pitchford, E. Codling, D. Psarra; \emph{Uncertainty and sustainability in fisheries and the benefit of marine protected areas}, Ecol. Mod. \textbf{207} (2007) 286-292. \bibitem{So} J. Sobel, G. Dahlgren; \emph{Marine Reserves: A guide to science, design, and use}, Island Press, Washington, 2004. \bibitem{Sum} U. Sumaila, S. Guenette, J. Alder, R. Chuenpagdee; \emph{Addressing ecosystem effects of fishing using marine protected areas}, ICES J. of Marine Sci. \textbf{57} (2000) 752-760. \end{thebibliography} \end{document}