# A Comparison of the Trojan Y Chromosome Strategy to Harvesting Models for
Eradication of Non-Native Species

###### Abstract.

The Trojan Y Chromosome Strategy (TYC) is a promising eradication method for
biological control of non-native species. The strategy works by manipulating
the sex ratio of a population through the introduction of supermales that
guarantee male offspring. In the current manuscript, we compare the TYC method
with a pure harvesting strategy. We also analyze a hybrid harvesting model
that mirrors the TYC strategy. The dynamic analysis leads to results on
stability of solutions and bifurcations of the model. Several conclusions
about the different strategies are established via optimal control methods. In
particular, the results affirm that either a pure harvesting or hybrid
strategy may work better than the TYC method at controlling a non-native
species population.

###### Key words and phrases:

mating system, stability and bifurcation, optimal control, biological
invasions, biological control

###### 1991 Mathematics Subject Classification:

Primary: 34C11, 34C23, 49J15; Secondary: 92D25, 92D40

Jingjing Lyu1, Pamela J. Schofield2, Kristen M. Reaver3,

Matthew Beauregard4 and Rana D. Parshad5

1) Department of Mathematics,

Chengdu University,

Chengdu, Sichuan 610106, China.

2) U.S. Geological Survey,

Wetland and Aquatic Research Center,

Gainesville, FL 32653, USA.

3) Cherokee Nation Technologies,

Wetland and Aquatic Research Center,

Gainesville, FL 32653, USA

4) Department of Mathematics,

Stephen F. Austin State University,

Nacogdoches, TX 75962 , USA

5) Department of Mathematics,

Iowa State University,

Ames, IA 50011, USA

## Recommendations for Resource Managers

  * •

Where harvesting is feasible, it is as effective if not more effective than
the classical TYC method. Therein managers may attempt harvesting female fish
while stocking males or harvesting both male and female fish.

  * •

Managers may attempt linear harvesting, saturating density dependent
harvesting and unbounded density dependent harvesting. Linear harvesting is
seen to be the most effective.

  * •

We caution against the outright use of harvesting due to various density
dependent effects that may arise. To this end hybrid models that involve a
combination of harvesting and TYC type methods might be a better strategy.

  * •

One may also use harvesting as a tool in mesocosm settings to predict the
efficacy of the TYC strategy in the wild.

##  1\. Introduction

###  1.1. Background

Biological invasions are the “uncontrolled spread and proliferation of species
to areas outside their native range” [1]. The rate of such invasions in the
United States continues to rise and, subsequently, the financial and
ecological damage caused by them [2, 3]. Non-native species can be difficult
to manage [4, 5]. For such reasons, the spread and control of non-native
species is an important and timely problem in spatial ecology and much work
has been devoted to this issue [6, 7, 8, 9, 10, 11, 12, 13, 14]. Current
eradication efforts for invasive aquatic species usually involve chemical
treatment, local harvesting, dewatering, ichthyocides, or a suitable
combination [15]. Unfortunately, all of these methods are known to negatively
impact native fauna [15].

An _alternative_ method that has been proposed for eradication of non-native
species is the Trojan Y Chromosome strategy (TYC) [17, 18]. Unlike other
technological approaches, the TYC strategy does not require within-chromosome
genetic modifications, rather, it involves a reassortment of (whole) pre-
existing sex chromosomes among individuals and is not considered a genetically
modified organism (GMO) [19]. These manipulations can cease at any time and
therefore the strategy is reversible. The strategy works by adding feminized
males and/or feminized super males (containing two Y chromosomes) to an
existing invasive population, to skew the sex ratio of subsequent generations
to contain an increasing number of males (i.e., fewer and fewer females in
each generation). The gradual reduction in females may lead to eventual
extinction of the population (see the right panel of Fig. 1). This strategy
has been of much interest lately [18, 20, 21, 22, 23, 24, 25, 26, 27, 29, 30,
31, 32].

Figure 1. Left panel: the female/male density change with time without TYC
strategy. Right Panel: the female/male density change with time by introducing
YY supermales (TYC strategy).

Fig. 1 is a simple demonstration of the theoretical power of the TYC strategy.
There is recent laboratory and field progress on TYC from two groups. (1) Dr.
Pamela Schofield’s laboratory, at the U.S. Geological Survey Wetland (USGS)
and Aquatic Research Center in Gainesville, Florida is working on developing
YY males for several fishes, including the fancy guppy fish (Poecilia
reticulata). (2) The Eagle Fish Genetics Lab of the Idaho Department of Fish
and Game has produced large stocks of YY males, for brook trout (_Salvelinus
fontinalis_) [34]. These have been released in the field with promising
preliminary results [35]. However, personal communication with chief
scientists at the above labs [36] suggests that the challenges in putting TYC
into practice are (1) the design and production of the $YY$ males and females
and (2) _in vitro_ testing of the strategy. The goal of this manuscript is to
compare harvesting strategies with TYC. We wanted to know whether TYC, a new
technology, could out compete the traditional strategy of harvesting in
reducing populations of non-native species. Harvesting has been used to reduce
populations of non-native species and subsequently moderate their negative
effects on environments [37, 38]. The use of harvesting is restricted to
certain field situations where the target organisms can be encountered,
detected and removed in a practical manner.

Our specific objectives were:

  1. (1)

Investigate alternate biological control strategies that _do not_ require $YY$
males or females.

  2. (2)

Compare and contrast such strategies to the TYC strategy, via an optimal
control approach.

  3. (3)

Develop such strategies so that they could be used in _conjunction_ with the
established TYC strategy as a hybrid strategy or become a _novel_ strategy in
itself.

  4. (4)

To use preliminary population data from mesocosm experiments to establish our
models. Note that these experiments contain wild type males and females only.
There are no YY males present in the mesocosm.

  5. (5)

To compare and contrast our harvesting strategies to the TYC strategy, using
realistic parameters that are _outputs of the above mentioned mesocosm
experiments_.

  6. (6)

Establish a mesocosm framework via harvesting - that could in turn be used to
test the efficacy of the TYC strategy in the wild.

For our mesocosm experiments we have chosen to work with guppy (_Poecilia
reticulata_) in the laboratory and get the data to model its population
dynamics herein. Our experiments contain only wild type male and female
guppies. Guppy is a tropical ornamental fish that is popular as an aquarium
pet. This species is introduced intermittently across the USA and has
established local reproducing populations in some areas [39]. Guppy is used
widely in laboratory studies and is especially well-suited to our experiments
due to its short generation time (ca. four weeks until sexual maturity),
sexual dimorphism, and peaceful nature.

###  1.2. Three-Variables TYC Model

The three-variables TYC model, in which only a YY supermale is introduced, is
described by a system of three ordinary differential equations for state
variables: a wild-type XX female ($f$), a wild-type XY male ($m$) and a YY
supermale ($s$):

(1) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta f,\\\
&\frac{dm}{dt}=\left(\frac{1}{2}fm+fs\right)\beta L-\delta m,\\\
&\frac{ds}{dt}=\mu-\delta s,\end{split}$

where $f,~{}m,$ and $s$ are the population densities; $\beta$ represents the
birth rate; $\delta$ is the death rate; $\mu$ denotes the non-negative
introduction rate of YY supermales. The logistic term is given by

(2) $L=1-\frac{f+m+s}{K},\quad\quad\quad\quad\quad\quad\quad\quad\quad$

where $K$ is the carrying capacity of the ecosystem. It is assumed that
$0<\beta$ and $\delta<1$ in this manuscrript.

###  1.3. Preliminary Data

Eleven months of preliminary population data were collected from a mesocosm
experiment conducted at the USGS facility in Gainesville, Florida (Fig. 2,
Fig. 3). Note these experiments consisted of wild type males and wild type
females only. There are no YY males present in the mesocosm. Our goal was to
infer life history parameters for _P. reticulata_ , such as birth and death
rates and carrying capacity, from this data - and then use that to simulate
TYC type models. The experiment was conducted in an indoor laboratory in one
rectangular fiberglass tank (114 x 56 x 61 cm) with aerated well water at a
depth of 30 cm (192L of total volume). Period water temperature fluctuated
with ambient indoor temperature and ranged from approximately 22-25 C. Living
aquatic vegetation (Hydrilla) was added to provide cover for the fish and to
seed live food sources such as microcrustaceans. The fish were fed three times
weekly with commercial flake food to supplement the live food sources. No
predators or inter-specific competitors were present in the mesocosms. A total
of 30 wild-type fancy guppies (15 XY males and 15 XX females) were initially
introduced to the mesocosm. They were allowed to reproduce in the mesocosm for
11 months. Population counts were made monthly (note that the guppy generation
time is $\approx$ 1 month) - see right panel in Fig. 3.

The population data is best fit to the basic population model where $s=0$
(without introduced YY super males) in (1). This enables the various
population parameters to be inferred, see right panel in Fig. 3. The best fit
parameters are $\beta=0.0057,~{}\delta=.0.0648,$ and $K=405.$

Figure 2. Male (top) and female (bottom) fancy guppy (Poecilia reticulata).
Photos were taken by Howard Jelk, USGS.

Figure 3. Left panel: This is a tank from the USGS facility where the mesocosm
experiment was run with fancy guppies (Poecilia reticulata). Right panel:
preliminary population data over 11-month period. The smooth curve is a result
of simulating the basic population model (1) (with $s=0$) with best-fit
parameters $\beta=0.0057,~{}\delta=0.0648,$ and $K=405$.

##  2\. Equilibrium and Stability Analysis

###  2.1. Equilibrium and Stability Analysis when $\boldsymbol{\mu>0}$

We relegate a significant amount of the local stability analysis to the
appendix see section 9. We begin by stating the following theorem,

######  Theorem 2.1.

Let $\mu>0$. The boundary equilibrium, $(0,0,\frac{\mu}{\delta})$, of system
(1) is locally stable if $\frac{1}{9}<\delta<1$ and unstable if
$0<\delta\leq\frac{1}{9}$.

###### Proof.

Now consider the boundary equilibrium of model (1), which satisfies $f^{*}=0$,
is $(0,0,\frac{\mu}{\delta})$. Now evaluate (21) at
$(0,0,\frac{\mu}{\delta})$. Then the characteristic equation about
$(0,0,\frac{\mu}{\delta})$ is

$\lambda^{3}+3\delta\lambda^{2}+3\delta^{2}\lambda+\delta^{2}=0,$

with corresponding eigenvalues

$\begin{split}&\lambda_{1}=(\delta^{3}-\delta^{2})^{\frac{1}{3}}-\delta,\\\
&\lambda_{2}=-\frac{1}{2}(\delta^{3}-\delta^{2})^{\frac{1}{3}}-\delta+\frac{\sqrt{3}}{2}(\delta^{3}-\delta^{2})^{\frac{1}{3}}i,\\\
&\lambda_{3}=-\frac{1}{2}(\delta^{3}-\delta^{2})^{\frac{1}{3}}-\delta-\frac{\sqrt{3}}{2}(\delta^{3}-\delta^{2})^{\frac{1}{3}}i.\end{split}$

Since $0<\delta<1$, then $\lambda_{1}<0$. The inequality
$-\frac{1}{2}(\delta^{3}-\delta^{2})^{\frac{1}{3}}-\delta<0$ implies that
$\delta>\frac{1}{9}$. According to the Routh Hurwitz stability criteria, the
system (1) under $f^{*}=0$ is locally stable if $\delta>\frac{1}{9}$. ∎

###  2.2. Equilibrium and Stability Analysis when $\boldsymbol{\mu=0}$

If $\mu=0$, then $s^{*}=0$. By direct computations, we can verify that the
equilibrium of model (1) is of the form

$(f^{*},m^{*},0),\quad\mathrm{with}\quad f^{*}=m^{*}.$

Let
$f^{*}_{\pm}=m^{*}_{\pm}=\frac{K}{4}\pm\frac{K}{4}\sqrt{1-\frac{16\delta}{\beta
K}}$. If $16\delta=\beta K$, then the model (1) has two equilibria for which
$f^{*}=\frac{K}{4}$ and $0$, respectively. If $16\delta<\beta K$, the model
has 3 different equilibria with $f^{*}=f^{*}_{+},~{}f^{*}_{-}$ and $0$. If
$16\delta>\beta K$, the model has only one trivial equilibrium, $f^{*}=0$.

In the case $16\delta=\beta K$, the equilibria are
$(\frac{K}{4},\frac{K}{4},0)$ and $(0,0,0)$, repectively. The corresponding
characteristic equation about $(\frac{K}{4},\frac{K}{4},0)$ is

$\lambda^{3}+(3\delta-\frac{K\beta}{16})\lambda^{2}+(3\delta^{2}-\frac{K\beta\delta}{8})\lambda-\frac{\delta^{2}(K\beta-16\delta)}{16}=0,$

with eigenvalues

$\begin{split}&\lambda_{1}=-\delta<0,\\\ &\lambda_{2}=-\delta<0,\\\
&\lambda_{3}=\frac{1}{16}K\beta-\delta=0.\end{split}$

By Routh Hurwitz criterion, $(\frac{K}{4},\frac{K}{4},0)$ is unstable.
Similarly, we find $(0,0,0)$ is locally stable.

In the case $16\delta<\beta K$, the equilibria are
$(f^{*}_{+},m^{*}_{+},0),(f^{*}_{-},m^{*}_{-},0)$ and $(0,0,0)$. The
characteristic equation about $(f^{*}_{+},m^{*}_{+},0)$ is

$(\lambda+\delta)^{2}\left(\lambda+\frac{1}{8}K\beta\sqrt{\frac{K\beta-16\delta}{k\beta}}+\frac{1}{8}K\beta-2\delta\right)=0$

The corresponding eigenvalues are $\lambda_{1}=\lambda_{2}=-\delta<0$, and
$\lambda_{3}=-\frac{1}{8}K\beta\sqrt{\frac{K\beta-16\delta}{k\beta}}-\frac{1}{8}K\beta+2\delta$.
Since $16\delta<\beta K$ then $\lambda_{3}<0$. This implies
$(f^{*}_{+},m^{*}_{+},0)$ is locally stable.

Similarly, we find the characteristic equation about $(f^{*}_{-},m^{*}_{-},0)$
is

$(\lambda+\delta)^{2}\left[\lambda-\left(\frac{1}{8}K\beta\sqrt{\frac{K\beta-16\delta}{K\beta}}-\frac{1}{8}K\beta+2\delta\right)\right]=0.$

It is easy to verify the third eigenvalue
$\left(\frac{1}{8}K\beta\sqrt{\frac{K\beta-16\delta}{K\beta}}-\frac{1}{8}K\beta+2\delta\right)>0$
under $16\delta<\beta K$. Therefore $(f^{*}_{-},m^{*}_{-},0)$ is unstable.

In the case of $16\delta>\beta K$, the only equilibrium is $(0,0,0)$, which is
locally stable.

######  Remark 1.

The three and four species TYC models are now known to blow-up in finite time
[28]. Thus one must exhibit caution when dealing with such models by
restricting initial conditions, as well as the rate of feminized
male/feminized super male introduction $\mu$. In the current manuscript we
restrict initial data and the size of $\mu$, s.t. we always have positive
bounded solutions, and all of the ensuing optimal control theory can be
applied.

###  2.3. Optimal Control Analysis

The goal of this section is to investigate the mechanisms in our TYC system of
equations, that, if controlled, could lead to optimal levels of both wild type
female and male densities. We assume that the introduction rate $\mu$ is not
known _a priori_ and enter the system as a time-dependent control. The
response for the range is $0\leq\mu(t)<\infty$.

Consider the following objective function

$J_{0}(\mu)=\int^{T}_{0}-(f+m)-\frac{1}{2}\mu^{2}dt,$

subject to the governing equations and initial conditions. Optimal strategies
are derived for the objective function, where we want to minimize both female
and male populations and also minimizing the YY males introduction rate $\mu$.
Optimal controls are searched for within the set $U_{0}$, namely,

$U_{0}=\\{\mu~{}|~{}\mu\ \mbox{measurable},\ 0\leq\mu<\infty,\ t\in[0,T],\
\forall T\\}.$

The goal is to seek an optimal $\mu^{*}(t)$ such that,

$\begin{split}J_{0}(\mu^{*})&=\underset{\mu}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}\mu^{2}dt\\\
&=\underset{\mu}{\min}\int^{T}_{0}f+m+\frac{1}{2}\mu^{2}dt.\end{split}$

Consider the following existence theorem,

######  Theorem 2.2.

Consider the optimal control problem (1) with $\mu=\mu(t)$. There exists
$\mu^{*}(t)\in U_{0}$ such that

(3)
$J_{0}(\mu^{*})=\underset{\mu}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}\mu^{2}dt\\\
$

###### Proof.

The compactness (closed and bounded in the ODE case) of the functional $J$
follows from the global boundedness of the state variables and the control
$\mu$. Also the functional $J$ is concave in the argument $\mu$. This is
easily verified via standard application [40]. These facts in conjunction give
the existence of an optimal control. ∎

We use Pontryagin’s maximum principle to derive the necessary conditions on
the optimal control. The Hamiltonian for $J_{0}$ is given by

$H_{0}=-(f+m)-\frac{1}{2}\mu^{2}+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}+\lambda_{3}s^{\prime}.$

We use the Hamiltonian to find a differential equation of the adjoint
$\lambda_{i},i=1,2,3.$ Namely,

$\begin{split}&\lambda_{1}^{\prime}(t)=1-\lambda_{1}[\frac{m\beta}{2}(1-\frac{f+m+s}{K})-\frac{fm\beta}{2K}-\delta]\\\
&\quad\quad\quad-\lambda_{2}[(\frac{m}{2}+s)\beta(1-\frac{f+m+s}{K})-\frac{\beta}{K}(\frac{fm}{2}+fs)],\\\
&\lambda_{2}^{\prime}(t)=1-\lambda_{1}[\frac{f\beta}{2}(1-\frac{f+m+s}{K})-\frac{fm\beta}{2K}]\\\
&\quad\quad\quad-\lambda_{2}[\frac{f\beta}{2}(1-\frac{f+m+s}{K})-\frac{fm\beta}{2K}-\delta],\\\
&\lambda_{3}^{\prime}(t)=\frac{\lambda_{1}fm\beta}{2K}+\lambda_{2}[\frac{\beta}{K}(\frac{fm}{2}+fs)-\beta
f(1-\frac{f+m+s}{K})]+\lambda_{3}\delta,\end{split}$

with the transversality condition given as

$\lambda_{1}(T)=\lambda_{2}(T)=\lambda_{3}(T)=0.$

Now considering the optimality conditions, the Hamiltonian function is
differentiated with respect to control variable $\mu$ resulting in

$\frac{\partial H_{0}}{\partial\mu}=\lambda_{3}-\mu.$

Then a compact way of writing the optimal control $\mu$ is

(4) $\mu^{*}(t)=\max(0,\lambda_{3}).$

The following theorem encapsulates the above.

######  Theorem 2.3.

An optimal control $\mu^{*}\in U_{0}$ for the system (1) that maximizes the
objective functional $J$ is characterized by (4).

##  3\. Mirroring TYC startegy via Harvesting

A key issue in the implementation of the TYC strategy is the “production” of
Trojan YY males.

######  Remark 2.

Here we present a model that can be reared effectively in the laboratory,
while mirroring the TYC strategy - but that does not require YY males/females.
Essentially the TYC strategy works by reducing the number of females in each
generation, whilst increasing the number of males. We implement this in the
current model where we remove a specified fraction of females whilst adding in
males. This strategy is called female harvesting male stocking (FHMS) and
attempts to mirror the manipulation of the sex-ratio via the TYC strategy. We
will introduce several forms of harvesting/stocking.

Consider the following model (We also note this model as Model 1)

(5) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta f-\eta_{1}f,\\\
&\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta m+\eta_{2}m,\end{split}$

with $L=1-\frac{f+m}{K}.$ Here, $\eta_{1}$ and $\eta_{2}$ are non-negative
parameters the specify the removal rate of females and addition rate of males,
respectively. It is assumed that $\delta>\eta_{2}$ if $\eta_{2}\not=0$.

###  3.1. Equilibria and Stability Analysis

An equilibrium of $f^{*}=0$ of this system requires $m^{*}=0$. It is easy to
verify that the trivial equilibrium $(0,0)$ is locally stable since both
eigenvalues of the linearized systems are negative.

Let $f^{*}>0$ and consider the three cases:

  1. (1)

$\eta_{1}>0$ and $\eta_{2}=0$,

  2. (2)

$\eta_{1}=0$ and $\eta_{2}>0$,

  3. (3)

$\eta_{1}>0$ and $\eta_{2}>0$.

####  3.1.1. Case 1: $\eta_{1}>0,\eta_{2}=0$.

In this situation the equilibrium of (5) is

$m^{*}=\left(1+\frac{\eta_{1}}{\delta}\right)f^{*}$

where $f^{*}$ satisfies

$\frac{1}{2}\left(2+\frac{\eta_{1}}{\delta}\right)\beta f^{2}-\frac{1}{2}\beta
Kf+K\delta=0.$

Define

$\Delta=(\frac{1}{2}\beta K)^{2}-2(2+\frac{\eta_{1}}{\delta})\beta K\delta.$

If $\Delta=0$, that is, $\beta K=8(2\delta+\eta_{1})$, then
$f^{*}=\frac{K}{2(2+\frac{\eta_{1}}{\delta})}=\frac{4\delta}{\beta}$ and
$m^{*}=\frac{4(\delta+\eta_{1})}{\beta}$. In such a case, the corresponding
characteristic equation is

(6) $\lambda^{2}+\frac{2\delta(\delta+\eta_{1})\lambda}{2\delta+\eta_{1}}=0.$

Since $0$ is a root of (6) then by Routh Hurwitz Criterion,
$\left(\frac{4\delta}{\beta},\frac{4(\delta+\eta_{1})}{\beta}\right)$ is not
stable.

If $\Delta>0$, that is, $\beta K>8(2\delta+\eta_{1})$, then we have the
equilibrium solutions $(f^{*}_{+},m^{*}_{+})$ and $(f^{*}_{-},m^{*}_{-})$
where

$\displaystyle f^{*}_{\pm}$ $\displaystyle=$ $\displaystyle\frac{(\beta
K\pm\sqrt{\beta K(\beta
K-16\delta-8\eta_{1}})\delta}{2(2\delta+\eta_{1})\beta},$ $\displaystyle
m^{*}_{\pm}$ $\displaystyle=$ $\displaystyle\frac{(\delta+\eta_{1})(\beta
K\pm\sqrt{\beta K(\beta K-16\delta-8\eta_{1})})}{2(2\delta+\eta_{1})\beta}.$

Notice that both $f^{*}_{+}$ and $f^{*}_{-}$ are positive. The characteristic
equation about $(f^{*}_{+},m^{*}_{+})$ is given by

(7) $\lambda^{2}+k_{41}\lambda+k_{40}=0,$

where

$\displaystyle k_{41}$ $\displaystyle=$
$\displaystyle\delta\frac{(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}}{2(2\delta+\eta_{1})^{2}}+\delta\frac{K\beta\delta+K\beta\eta_{1}-8\delta^{2}-12\delta\eta_{1}-4\eta_{1}^{2}}{2(2\delta+\eta_{1})^{2}},$
$\displaystyle k_{40}$ $\displaystyle=$
$\displaystyle\delta\frac{(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}}{4(2\delta+\eta_{1})}+\delta\frac{K\beta\delta+K\beta\eta_{1}-16\delta^{2}-24\delta\eta_{1}-8\eta_{1}^{2}}{4(2\delta+\eta_{1})}.$

By Routh Hurwitz criterion, the system with characteristic equation (7) is
locally stable if and only if both coefficients satisfy $k_{4i}>0$ for
$i=0,1$. Since $\beta K>8(2\delta+\eta_{1})$, then
$\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}>0$. Likewise,

$\displaystyle
K\beta\delta+K\beta\eta_{1}-8\delta^{2}-12\delta\eta_{1}-4\eta_{1}^{2}$
$\displaystyle\quad>8\delta(2\delta+\eta_{1})+8\eta_{1}(2\delta+\eta_{1})-8\delta^{2}-12\delta\eta_{1}-4\eta_{1}^{2}$
$\displaystyle\quad=16\delta^{2}+8\delta\eta_{1}+16\eta_{1}\delta+8\eta_{1}^{2}-8\delta^{2}-12\delta\eta_{1}-4\eta_{1}^{2}$
$\displaystyle\quad=8\delta^{2}+12\delta\eta_{1}+4\eta_{1}^{2}>0.$

Similarly,

$\displaystyle
K\beta\delta+K\beta\eta_{1}-16\delta^{2}-24\delta\eta_{1}-8\eta_{1}^{2}$
$\displaystyle\quad>8\delta(2\delta+\eta_{1})+8\eta_{1}(2\delta+\eta_{1})-16\delta^{2}-24\delta\eta_{1}-8\eta_{1}^{2}$
$\displaystyle\quad=16\delta^{2}+8\delta\eta_{1}+16\eta_{1}\delta+8\eta_{1}^{2}-16\delta^{2}-24\delta\eta_{1}-8\eta_{1}^{2}=0.$

This implies that $k_{41}>0$ and $k_{40}>0$. Subsequently, the equilibrium
$(f^{*}_{+},m^{*}_{+})$ is locally stable. Now consider the characteristic
equation about $(f^{*}_{-},m^{*}_{-})$, namely,

$\lambda^{2}+k_{51}\lambda+k_{50}=0.$

with

$\displaystyle k_{51}$ $\displaystyle=$
$\displaystyle\delta\frac{K\beta\delta+K\beta\eta_{1}-8\delta^{2}-12\delta\eta_{1}-4\eta_{1}^{2}-(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}}{2(2\delta+\eta_{1})^{2}},$
$\displaystyle k_{50}$ $\displaystyle=$
$\displaystyle\delta\frac{K\beta\delta+K\beta\eta_{1}-16\delta^{2}-24\delta\eta_{1}-8\eta_{1}^{2}-(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}}{4(2\delta+\eta_{1})}.$

The signs of $k_{51}$ and $k_{50}$ need to be known to quantify the stability
of the equilibrium solution. Firstly, consider the numerator of $k_{50}$,

$\displaystyle
K\beta\delta+K\beta\eta_{1}-16\delta^{2}-24\delta\eta_{1}-8\eta_{1}^{2}-(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}$
$\displaystyle\quad=(\delta+\eta_{1})(K\beta-(16\delta+8\eta_{1}))-(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}.$

To compare $(\delta+\eta_{1})(K\beta-(16\delta+8\eta_{1}))$ and
$(\delta+\eta_{1})\sqrt{K\beta(K\beta-16\delta-8\eta_{1})}$, it is enough to
compare $(K\beta-(16\delta+8\eta_{1}))^{2}$ and
$(\sqrt{K\beta(K\beta-16\delta-8\eta_{1})})^{2}$. Namely,

$\begin{split}&(K\beta-(16\delta+8\eta_{1}))^{2}-(\sqrt{K\beta(K\beta-16\delta-8\eta_{1})})^{2}\\\
&\quad=(K\beta)^{2}-2K\beta(16\delta+8\eta_{1})+(16\delta+8\eta_{1})^{2}-(K\beta)^{2}+K\beta(16\delta+8\eta_{1})\\\
&\quad=(16\delta+8\eta_{1})^{2}-K\beta(16\delta+8\eta_{1})\\\
&\quad=(16\delta+8\eta_{1})(16\delta+8\eta_{1}-K\beta)<0.\end{split}$

This is less than zero because $\beta K>8(2\delta+\eta_{1}).$ Therefore,
$k_{50}<0$ and $(f^{*}_{-},m^{*}_{-})$ is unstable.

Lastly, if $\Delta<0$, that is, $\beta K<8(2\delta+\eta_{1})$, then the only
equilibrium is trivial which is locally stable.

####  3.1.2. Case 2: $\eta_{1}=0,\eta_{2}>0$.

In the case of $\eta_{1}=0,\eta_{2}>0$, the model simplifies to

(8) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta f,\\\
&\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta m+\eta_{2}m,\end{split}$

which is symmetric to the case $\eta_{1}>0,\eta_{2}=0$ if we replace
$\eta_{2}$ by $-\eta_{1}$. Therefore, the following corollary can be
developed.

######  Corollary 1.

The boundary equilibrium $(0,0)$ is locally stable. There are 3 cases for the
interior equilibria,

  1. (i)

If $\beta K=8(2\delta-\eta_{2})$, the equilibrium of (8) is
$\left(\frac{4(\delta-\eta_{2})}{\beta},\frac{4\delta}{\beta}\right)$ and it
is locally stable.

  2. (ii)

If $\beta K>8(2\delta-\eta_{2})$, the equilibria of (8) are

$\left(\frac{(\delta-\eta_{2})(\beta K\pm\sqrt{\beta K(\beta
K-16\delta+8\eta_{2})})}{2(2\delta+\eta_{1})\beta},\frac{(\beta
K\pm\sqrt{\beta K(\beta
K-16\delta+8\eta_{2}})\delta}{2(2\delta-\eta_{2})\beta}\right)$ which the
postive and negative branches are locally stable and unstable, respectively.

  3. (iii)

If $\beta K<8(2\delta-\eta_{2})$, the only equilibrium of (8) is $(0,0)$.

####  3.1.3. Case 3: $\eta_{1}>0,\eta_{2}>0$.

In the case of $\eta_{1}>0,\eta_{2}>0$, the equilibrium of model (5) is given
by

$m^{*}=\frac{\delta+\eta_{1}}{\delta-\eta_{2}}f^{*},$

where $f^{*}$ satisfies

$(-2\beta\delta-\beta\eta_{1}+\beta\eta_{2})f^{2}+(K\beta\delta-K\beta\eta_{2})f-2K\delta^{2}+4K\delta\eta_{2}-2K\eta_{2}^{2}=0.$

If $\beta K=8(2\delta+\eta_{1}-\eta_{2})$, then
$f^{*}=\frac{K\beta(\delta-\eta_{2})}{2\beta(2\delta+\eta_{1}-\eta_{2})}$ and
$m^{*}=\frac{K\beta(\delta+\eta_{1})}{2\beta(2\delta+\eta_{1}-\eta_{2})}$. The
corresponding characteristic equation is

$\lambda\left(\lambda+\frac{2(\delta-\eta_{2})(\delta+\eta_{1})}{2\delta+\eta_{1}-\eta_{2}}\right)=0.$

It is easy to see
$(\frac{K\beta(\delta-\eta_{2})}{2\beta(2\delta+\eta_{1}-\eta_{2})},\frac{K\beta(\delta+\eta_{1})}{2\beta(2\delta+\eta_{1}-\eta_{2})})$
is unstable because one of the eigenvalues is $0$.

If $\beta K>8(2\delta+\eta_{1}-\eta_{2})$, then we have two positive
equilibria, $(f^{*}_{+},m^{*}_{+}$ and $(f^{*}_{-},m^{*}_{-})$, where

$\begin{split}&f^{*}_{\pm}=\frac{(K\beta\pm\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})})(\delta-\eta_{2})}{2\beta(2\delta+\eta_{1}-\eta_{2})},\\\
&m^{*}_{\pm}=\frac{(K\beta\pm\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})})(\delta+\eta_{1})}{2\beta(2\delta+\eta_{1}-\eta_{2})}.\end{split}$

The characteristic equation about $(f^{*}_{+},m^{*}_{+})$ is

$\lambda^{2}+k_{61}\lambda+k_{60}=0$

with

$\begin{split}k_{61}=&\frac{-8\delta^{3}-12\delta^{2}\eta_{1}+12\delta^{2}\eta_{2}-4\delta\eta_{1}^{2}+16\delta\eta_{1}\eta_{2}-4\delta\eta_{2}^{2}+4\eta_{1}^{2}\eta_{2}-4\eta_{1}\eta_{2}^{2}}{2(2\delta+\eta_{1}-\eta_{2})^{2}}\\\
&+\frac{\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})}{2(2\delta+\eta_{1}-\eta_{2})^{2}}\\\
&+\frac{K\beta(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})}{2(2\delta+\eta_{1}-\eta_{2})^{2}}\end{split}$

and

$\begin{split}k_{60}=&\frac{-16\delta^{3}-24\delta^{2}\eta_{1}+24\delta^{2}\eta_{2}-8\delta\eta_{1}^{2}+32\delta\eta_{1}\eta_{2}-8\delta\eta_{2}^{2}+8\eta_{1}^{2}\eta_{2}-8\eta_{1}\eta_{2}^{2}}{4(2\delta+\eta_{1}-\eta_{2})}\\\
&+\frac{\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})}{4(2\delta+\eta_{1}-\eta_{2})}\\\
&+\frac{K\beta(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})}{4(2\delta+\eta_{1}-\eta_{2})}.\end{split}$

Since $\beta K>8(2\delta+\eta_{1}-\eta_{2})$ and $\delta>\eta_{2}$, then

$\begin{split}&\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})\\\
&\quad=\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}(\delta+\eta_{1})(\delta-\eta_{2})>0.\end{split}$

Likewise, the remaining two terms in the numerator of $k_{61}$ can be combined
and shown to be greater than zero, namely,

$\begin{split}&-8(\delta+\eta_{1})(\delta-\eta_{2})(\delta+\frac{1}{2}\eta_{1}-\frac{1}{2}\eta_{2})+K\beta(\delta+\eta_{1})(\delta-\eta_{2})\\\
&\quad=(\delta+\eta_{1})(\delta-\eta_{2})(K\beta-8(\delta+\frac{1}{2}\eta_{1}-\frac{1}{2}\eta_{2}))\\\
&\quad>(\delta+\eta_{1})(\delta-\eta_{2})(K\beta-16(\delta+\frac{1}{2}\eta_{1}-\frac{1}{2}\eta_{2}))\\\
&\quad=(\delta+\eta_{1})(\delta-\eta_{2})(K\beta-8(2\delta+\eta_{1}-\eta_{2}))>0.\end{split}$

Therefore $k_{61}>0$. In addition, since

$\begin{split}&-16\delta^{3}-24\delta^{2}\eta_{1}+24\delta^{2}\eta_{2}-8\delta\eta_{1}^{2}+32\delta\eta_{1}\eta_{2}-8\delta\eta_{2}^{2}+8\eta_{1}^{2}\eta_{2}-8\eta_{1}\eta_{2}^{2}\\\
&+K\beta(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})\\\
&=(\delta+\eta_{1})(\delta-\eta_{2})(\beta
K-8(2\delta+\eta_{1}-\eta_{2}))>0,\end{split}$

then $k_{60}>0.$ Therefore, $(f^{*}_{+},m^{*}_{+})$ is locally stable.

As for $(f^{*}_{-},m^{*}_{-})$, the corresponding characteristic equation is

$\lambda^{2}+k_{71}\lambda+k_{70}=0$

where

$\begin{split}k_{71}=&\frac{-8\delta^{3}-12\delta^{2}\eta_{1}+12\delta^{2}\eta_{2}-4\delta\eta_{1}^{2}+16\delta\eta_{1}\eta_{2}-4\delta\eta_{2}^{2}+4\eta_{1}^{2}\eta_{2}-4\eta_{1}\eta_{2}^{2}}{2(2\delta+\eta_{1}-\eta_{2})^{2}}\\\
&+\frac{\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}(-\delta^{2}-\delta\eta_{1}+\delta\eta_{2}+\eta_{1}\eta_{2})}{2(2\delta+\eta_{1}-\eta_{2})^{2}}\\\
&+\frac{K\beta(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})}{2(2\delta+\eta_{1}-\eta_{2})^{2}}\end{split}$

and

$\begin{split}k_{70}=&\frac{-16\delta^{3}-24\delta^{2}\eta_{1}+24\delta^{2}\eta_{2}-8\delta\eta_{1}^{2}+32\delta\eta_{1}\eta_{2}-8\delta\eta_{2}^{2}+8\eta_{1}^{2}\eta_{2}-8\eta_{1}\eta_{2}^{2}}{4(2\delta+\eta_{1}-\eta_{2})}\\\
&+\frac{\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}(-\delta^{2}-\delta\eta_{1}+\delta\eta_{2}+\eta_{1}\eta_{2})}{4(2\delta+\eta_{1}-\eta_{2})}\\\
&+\frac{K\beta(\delta^{2}+\delta\eta_{1}-\delta\eta_{2}-\eta_{1}\eta_{2})}{4(2\delta+\eta_{1}-\eta_{2})}.\end{split}$

The denominator of $k_{70}$ and $k_{71}$ are positive. The numerator of
$k_{70}$ can be simplified to

(9)
$\begin{split}(\delta+\eta_{1})(\delta-\eta_{2})(K\beta-16\delta-8\eta_{1}+8\eta_{2}-\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}).\end{split}$

The sign of (9) depends on
$(K\beta-16\delta-8\eta_{1}+8\eta_{2})-\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})}$.
It is enough to compare $(K\beta-16\delta-8\eta_{1}+8\eta_{2})^{2}$ and
$(\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})})^{2}$ because both of
them are positive. Therefore,

$\begin{split}&(K\beta-16\delta-8\eta_{1}+8\eta_{2})^{2}-(\sqrt{K\beta(K\beta-16\delta-8\eta_{1}+8\eta_{2})})^{2}\\\
&\quad=(K\beta-16\delta-8\eta_{1}+8\eta_{2})(K\beta-16\delta-8\eta_{1}+8\eta_{2}-K\beta)\\\
&\quad=(K\beta-16\delta-8\eta_{1}+8\eta_{2})(-16\delta-8\eta_{1}+8\eta_{2})\\\
&\quad=-8((\delta+\eta_{1})+(\delta-\eta_{2}))(K\beta-16\delta-8\eta_{1}+8\eta_{2})<0.\end{split}$

This shows that $k_{70}<0$. Subsequently, $(f^{*}_{-},m^{*}_{-})$ is unstable
by Routh Hurwitz criterion.

The local stability for the equilibria in each case is now understood. The
global stability of the trivial equilibrium $(0,0)$ is stated in the following
theorem.

######  Theorem 3.1.

Consider the model (5), the equilibrium $(0,0)$ is globally asymptotically
stable if $\beta K<2\delta+\eta_{1}-\eta_{2}.$

###### Proof.

Consider the Lyapunov function $V=f+m$. It is left to show that
$\frac{dV}{dt}<0$ for all $(f,m)\not=(0,0)$. Consider

$\begin{split}\frac{dV}{dt}&=\frac{df}{dt}+\frac{dm}{dt}\\\ &=fm\beta
L-(\delta+\eta_{1})f-(\delta-\eta_{2})m\\\ &\leq
fm\beta-(\delta+\eta_{1})f-(\delta-\eta_{2})m,\quad\mathrm{since}\quad L\leq
1\\\
&\leq[\beta-\frac{1}{m}(\delta+\eta_{1})-\frac{1}{f}(\delta-\eta_{2})]fm\\\
&\leq[\beta-\frac{1}{K}(\delta+\eta_{1})-\frac{1}{K}(\delta-\eta_{2})]fm.\end{split}$

It is enough to show
$\beta-\frac{1}{K}(\delta+\eta_{1})-\frac{1}{K}(\delta-\eta_{2})<0$. By direct
calculation, we have $\beta K<2\delta+\eta_{1}-\eta_{2}$, which proves theorem
(3.1). ∎

Various other forms of harvesting for FHMS are also tried, but we maintain the
same general structure that is females are harvested via a term
$-\eta_{1}G_{1}(f)$, whilst males are added to the system via a term
$\eta_{2}G_{2}(m)$. The results spanning various forms of $G_{1},G_{2}$ are
relegated to the appendix, section 10.

###  3.2. Optimal Control Analysis

The goal of this section is to further investigate controls the female-male
system. In particular, assume that the removal rate $\eta_{1}$ or the addition
rate $\eta_{2}$ are not known _a priori_ and enter the system as time-
dependent controls. The responses are over the range
$0\leq\eta_{1},\eta_{2}\leq 1$. For clarity, the model is restated with the
temporal dependence in $\eta_{i}$,

(10) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}(t)f,\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m+\eta_{2}(t)m.\end{split}$

Consider the objective function

$J_{1}(\eta_{1},\eta_{2})=\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.$

subject to the governing equations in (10) and specified initial conditions.
Optimal strategies are derived for the following objective function, where we
minimize both female and male populations while minimizing the harvesting and
addition rates. Optimal controls are searched for within the set $U_{1}$,
namely,

$U_{1}=\\{(\eta_{1},\eta_{2})|\eta_{i}\ \mbox{measurable},\ 0\leq\eta_{1}\leq
1,0\leq\eta_{2}<\delta,\ t\in[0,T],\ \forall T\\}.$

The goal is to seek an optimal $(\eta_{1}^{*},\eta_{2}^{*})$ such that

$J_{1}(\eta_{1}^{*},\eta_{2}^{*})=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}(-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt$

The following existence theorem is stated.

######  Theorem 3.2.

Consider the optimal control problem (10). There exists
$(\eta_{1}^{*},\eta_{2}^{*})\in U_{1}$ such that

$\begin{split}J_{1}(\eta_{1}^{*},\eta_{2}^{*})&=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt\\\
&=\underset{(\eta_{1},\eta_{2})}{\min}\int^{T}_{0}f+m+\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.\\\
\end{split}$

###### Proof.

Since $f+m\leq K,~{}0\leq\eta_{1}\leq 1$ and $0\leq\eta_{2}\leq\delta\leq 1$,
then the compactness of the functional $J$ follows from the global boundedness
of the state variables and the controls $\eta_{1}$ and $\eta_{2}$. Also the
functional $J$ is concave in the both of its arguments $\eta_{1}$ and
$\eta_{2}$. This is easily verified via standard application [40].
Subsequently, these facts in conjunction provide the existence of an optimal
control. ∎

We use the Pontryagin’s maximum principle to derive the necessary conditions
on the optimal controls. Consider the Hamiltonian for $J_{1}$, namely,

$H_{1}=-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}.$

The Hamiltonian is used to establish a differential equation of the adjoint
$\lambda_{i},i=1,2.$ That is,

$\begin{split}&\lambda_{1}^{\prime}(t)=1-\lambda_{1}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\eta_{1}]-\lambda_{2}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}],\\\
&\lambda_{2}^{\prime}(t)=1-\lambda_{1}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}]-\lambda_{2}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta+\eta_{2}].\end{split}$

with the transversality condition of $\lambda_{1}(T)=\lambda_{2}(T)=0.$ In
consideration of the optimality conditions, the Hamiltonian is differentiated
with respect to the control variables $\eta_{1}$ and $\eta_{2}$ resulting in

$\begin{split}&\frac{\partial
H_{1}}{\partial\eta_{1}}=-f\lambda_{1}-\eta_{1},\\\ &\frac{\partial
H_{1}}{\partial\eta_{2}}=m\lambda_{2}-\eta_{2}.\end{split}$

We find a characterization of $\eta_{1}$ by considering three cases,

  1. (1)

If $\frac{\partial H_{1}}{\partial\eta_{1}}<0$ then $\eta_{1}^{*}=0$. This
implies that $-f\lambda_{1}<0$;

  2. (2)

If $\frac{\partial H_{1}}{\partial\eta_{1}}=0$ then
$\eta_{1}^{*}=-f\lambda_{1}$. This implies that $0\leq-f\lambda_{1}\leq 1$;

  3. (3)

If $\frac{\partial H_{1}}{\partial\eta_{1}}>0$ then $\eta_{1}^{*}=1$. This
implies that $1<-f\lambda_{1}$.

Notice the recurrence of the expression $-f\lambda_{1}$. This expression is
strictly less than 0 when the control is at the lower bound. In constrast,
when the control is at the upper bound, then $-f\lambda_{1}>1$. Thus a compact
way of writing the optimal control $\eta_{1}$ is

(11) $\eta_{1}(t)^{*}=\min(1,\max(0,-f\lambda_{1})).$

Similarly, for $\eta_{2}$,

(12) $x\eta_{2}(t)^{*}=\min(1,\max(0,m\lambda_{2})).$

######  Theorem 3.3.

An optimal control $(\eta_{1}^{*},\eta_{2}^{*})\in U_{1}$ for the system (5)
that maximizes the objective functional $J$ is characterized by (11)-(12).

######  Remark 3.

The optimal control anlysis of the nonlinear forms of FHMS model is relegated
to the appendix, section 10.

##  4\. Optimal Strategy Comparisons between Classic TYC Model and FHMS Model

In this section numerical simulations are used to compare the optimal controls
for the classic TYC model to the FHMS model. The wild-type female and male
populations are compared for each control. The Guppy fish population data from
USGS is utilized to run the simulations with our best fit parameters, that is,
$\beta=0.005774,~{}\delta=.0.0648,$ and $K=405.0705.$

Table 1. Parameters used for numerical simulations Parameter | Description | Value
---|---|---
$\beta$ | Birth rate |  0.0057
$\delta$ | Death rate |  0.0648
$K$ | Carrying capacity | 405
$T$ | Terminal time | 200

(a)

(b)

Figure 4. (A) Female (red), male (green) and supermale (blue) densities and
optimal control $\mu(t)$ in (1) change with time $t$. (B) Female (red) and
male (green) densities and optimal controls $\eta_{1}(t)$ and $\eta_{2}(t)$ in
(10) change with time $t$.

(a)

(b)

Figure 5. (A) In this simulation we look at the objective function towards the
optimal implement for each model, we clearly see the cost function of TYC
model (Model 0) is much larger than FHMS model (Model 1). (B) Since the
constraints and scales for $\mu$ and $\eta_{1},\eta_{2}$ are different, it’s
hard to compare the effectiveness for each strategy, so here we also look at
the cost deducting the controls, that is, $\int^{T}_{0}f+m$ $dt$ under optimal
control.

In the classical TYC case, supermales are introduced at high rate for longer
time as shown in Fig. 4. Once the females are brought below a particular
threshold, the introduction of supermales decreases and finally turns off as
extinction nears. Similarly to the TYC strategy, the female removal rate
$\eta_{1}$ keeps high for a certain time to reduce the female population and
gradually declines as the entire population approaches extinction.

######  Remark 4.

From Fig. 4, we see that introducing males to the entire population is not
helpful for the eradication. Note, the optimal control $\eta_{2}$ turns out to
be essentially 0.

From Fig. 5, it is clear that the cost function of FHMS strategy at each time
is much less than the classical TYC strategy. We also consider that the range
for $\mu(0<\mu<\infty)$ and $\eta_{1},\eta_{2}(0\leq\eta_{1},\eta_{2}\leq 1)$
are different, we calculate the the cost excluding the controls under the
optimal control in each model as shown in the Fig. 5. And we get the following
result:

Table 2. Comparison of TYC and FHMS strategy results Results |  TYC |  FHMS
---|---|---
Objective function | 270820 | 84168
Cost deducting controls | 106782 | 84156
Female population in final time | 0.0184 | 0
Male population in final time | 0.005 | 0.0028
Female Approximate Eradication Time (month ) | 140 | 60
Male Approximate Eradication Time (month) | 156 | 110

######  Remark 5.

Note the approximate eradication time is defined as the time that the
population is less than some $\epsilon$ where $0<\epsilon<1.$ Through this
paper, we set $\epsilon=0.5$.

##  5\. Alternative Harvesting Approach (FHMH)

In this section we try an alternative harvesting approach.

######  Remark 6.

When we compare the linear harvesting model to the classical TYC model, it is
observed via Fig. 4, that the optimal levels $\eta_{2}$ of adding males in the
harvesting model, is virtually zero. This is also seen in our various
nonlinear models of harvesting, where we harvest females and introduce males,
see section 10. Motivated by this observation, we analyze the optimal
harvesting levels when both females and males are removed/harvested from the
population. We call these _female harvesting and male harvesting_ models
(FHMH).

The general model is given by

(13) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}G_{1}(f),\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m-\eta_{2}G_{2}(m),\end{split}$

where $G_{1}(f)$ and $G_{2}(m)$ are nonnegative functions. Again, assume that
$\eta_{1}$ and $\eta_{2}$ are nonnegative.

###  5.1. Model 4: $\boldsymbol{G_{1}(f)=f,G_{2}(m)=m}$

For clarity we restate the model:

(14) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta f-\eta_{1}f,\\\
&\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta m-\eta_{2}m.\end{split}$

####  5.1.1. Global stability

The global stability of the trivial equilibrium $(0,0)$ is stated in the
following theorem.

######  Theorem 5.1.

The trivial equilibrium is globally asymptotically stable in (14) if $\beta
K<2\delta+\eta_{1}+\eta_{2}$.

###### Proof.

Again, consider the Lyapunov function $V=f+m$. It is left to to show
$\frac{dV}{dt}<0$ for all $(f,m)\not=(0,0)$. Taking the derivative yields,

$\begin{split}\frac{dV}{dt}&=\frac{df}{dt}+\frac{dm}{dt}\\\ &=fm\beta
L-(\delta+\eta_{1})f-(\delta+\eta_{2})m\\\ &\leq
fm\beta-(\delta+\eta_{1})f-(\delta+\eta_{2})m,\quad\mathrm{since}\quad L\leq
1\\\
&\leq[\beta-\frac{1}{m}(\delta+\eta_{1})-\frac{1}{f}(\delta+\eta_{2})]fm\\\
&\leq[\beta-\frac{1}{K}(\delta+\eta_{1})-\frac{1}{K}(\delta+\eta_{2})]fm.\end{split}$

It is enough to show
$\beta-\frac{1}{K}(\delta+\eta_{1})-\frac{1}{K}(\delta+\eta_{2})<0$. By direct
calculation, we have $\beta K<2\delta+\eta_{1}+\eta_{2}$, which completes the
proof. ∎

######  Remark 7.

If one considers the condition for global asymptotic stability via theorem 5.1
that is $\beta K<2\delta+\eta_{1}+\eta_{2}$, to the condition for global
asymptotic stability via theorem 3.1 that is $\beta
K<2\delta+\eta_{1}-\eta_{2}$, we see this is stronger, and the weaker
condition via theorem 5.1 enables global asymptotic stability - or extinction.
Thus there is more merit in harvesting both males and females. Although
counter intuitive to TYC methodology, the reason is that when we introduce
normal males, there is always a chance they will mate with females - producing
more females.

####  5.1.2. Optimal control analysis

Assume that the removal rates are not known _a priori_ and enter the system as
time-dependent controls. Optimal controls are sought within the range
$0\leq\eta_{1},\eta_{2}\leq 1$. Consider the objective function,

$J_{4}(\eta_{1},\eta_{2})=\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt$

subject to (14) and with initial conditions $f(t_{0})=f_{0},m(t_{0})=m_{0}$.
Optimal controls are sought that minimize the populations and remain in the
set $U_{4}$, namely,

$U_{4}=\\{(\eta_{1},\eta_{2})|\eta_{i}\ \mbox{measurable},\ 0\leq\eta_{1}\leq
1,0\leq\eta_{2}\leq 1,\ t\in[0,T],\ \forall T\\}.$

Optimal functions $(\eta_{1}^{*},\eta_{2}^{*})$ are sought such that,

(15)
$\begin{split}J_{4}(\eta_{1}^{*},\eta_{2}^{*})&=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}(-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt\\\
&=\underset{(\eta_{1},\eta_{2})}{\min}\int^{T}_{0}f+m+\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.\\\
\end{split}$

Consider the following existence theorem.

######  Theorem 5.2.

Consider the optimal control problem (14). There exists
$(\eta_{1}^{*},\eta_{2}^{*})\in U_{4}$ such that

(16)
$J_{4}(\eta_{1}^{*},\eta_{2}^{*})=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt$

###### Proof.

Since $f+m\leq K,~{}0\leq\eta_{1}\leq 1$ and $0\leq\eta_{2}\leq 1$, the
compactness of the functional $J$ follows from the global boundedness of the
state variables and the controls $\eta_{1}$ and $\eta_{2}$. Also the
functional $J_{4}$ is concave in both arguments, $\eta_{1}$ and $\eta_{2}$.
These facts in conjunction give the existence of an optimal control. ∎

Pontryagin’s maximum principle is used to derive the necessary conditions on
the optimal controls. Consider the Hamiltonian for $J_{4}$, namely,

$H_{4}=-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}.$

The Hamiltonian is used to find a differential equation of the adjoint
$\lambda_{i},i=1,2.$ That is,

$\begin{split}&\lambda_{1}^{\prime}(t)=1-\lambda_{1}\left[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\eta_{1}\right]-\lambda_{2}\left[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}\right],\\\
&\lambda_{2}^{\prime}(t)=1-\lambda_{1}\left[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}\right]-\lambda_{2}\left[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\eta_{2}\right]\end{split}$

with the transversality condition $\lambda_{1}(T)=\lambda_{2}(T)=0$.

The Hamiltonian function is differentiated with respect to the control
variables $\eta_{1}$ and $\eta_{2}$ resulting in

$\begin{split}&\frac{\partial
H_{4}}{\partial\eta_{1}}=-f\lambda_{1}-\eta_{1},\\\ &\frac{\partial
H_{4}}{\partial\eta_{2}}=-m\lambda_{2}-\eta_{2}\end{split}$

As shown previously, a compact way of writing the optimal control for
$\eta_{i}$ is

(17) $\displaystyle\eta_{1}(t)^{*}$ $\displaystyle=$
$\displaystyle\min(1,\max(0,-f\lambda_{1})),$ (18)
$\displaystyle\eta_{2}(t)^{*}$ $\displaystyle=$
$\displaystyle\min(1,\max(0,-m\lambda_{2})).$

######  Theorem 5.3.

An optimal control $(\eta_{1}^{*},\eta_{2}^{*})\in U_{4}$ for the system (14)
that maximizes the objective functional $J$ is characterized by (17)-(18).

Various other forms of harvesting for FHMH are also tried, but we maintain the
same general structure that is females are harvested via a term
$-\eta_{1}G_{1}(f)$ and males are harvested to the system via a term
$-\eta_{2}G_{2}(m)$. The results spanning various forms of $G_{1},G_{2}$ are
relegated to the appendix, section 11.

###  5.2. Numerical Simulations and Comparisons

In this section, we will numerically simulate the optimal strategy and its
corresponding optimal states for the models 4-6 - the FHMH class of models.
For models 5,6 the reader is referred to section 11. The parameters for
simulating are given by:

Table 3. Parameters used for numerical simulations Parameter | Description | Value
---|---|---
$\beta$ | Birth rate |  0.0057
$\delta$ | Death rate |  0.0648
$K$ | Carrying capacity | 405
$T$ | Terminal time | 200
$d_{1}$ | Parameter In $G_{1}$ | 1
$d_{2}$ | Parameter In $G_{2}$ | 1

(a)

(b)

Figure 6. (A) Female (red) and male (green) densities change with time $t$
with optimal control $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) in Model
4. As we can see the controls and densities are indistinguishable. (B)
Objective function $J_{4}$ decrease as increasing time with optimal control.

(a)

(b)

Figure 7. (A) Female (red) and male (blue) densities change with time $t$ with
optimal control $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) in Model 5. (B)
Objective function $J_{5}$ varies as increasing time with optimal control.

(a)

(b)

Figure 8. (A) Female (red) and male (green) densities over time $t$ with
optimal control $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) in Model 6. We
can see that the densities and controls are indistinguishable. (B) Objective
function $J_{6}$ decrease as increasing time with optimal control.

Table. 4 compares all harvesting models including nonlinear forms of
harvesting for FHMS (Model 2-3) and FHMH (Model 5-6), which are refered to
section 10 and 11, respectively.

Table 4. Comparisons among Model 1-Model 6

Results | Model 1 | Model 2 | Model 3 | Model 4 | Model 5 | Model 6
---|---|---|---|---|---|---
Total cost | 84168 | 185660 | 31050 | 14250 | 93600 | 15790
Females population at final time | 0 | 171 | 0 | 0 | 3.2291 | 0
Males population at final time | 0.0028 | 180 | 0.0002 | 0 | 3.2966 | 0
Female Approximate Eradication Time | 60 | / | 2 | 7 | / | 9
Male Approximate Eradication Time | 110 | / | 73 | 7 | / | 8

##  6\. Conclusions

In this paper we derive and analyze optimal controls for TYC strategy. The
population can be driven to extinction with an optimal $\mu(t)$ that requires
a large initial introduction of YY super males. Our simulations, see Fig. 4,
show that one requires to introduce supermales at a continuous rate of
12-14$\%$ of the carrying capacity for at least a month, before this
introduction can catapult the population densities toward controllable levels,
for which extinction can occur. The difficulty in creating the supermale
population motivates alternative approaches, such as harvesting, that attempt
to mirror the TYC strategy. However, the effectiveness of driving a population
to extinction seems to depend delicately on the form of harvesting. With this
in mind we introduce two classes of models,

  * •

Harvesting females while stocking males - FHMS models. Herein we introduce
three subclasses, linear harvesting, saturating density dependent harvesting
and unbounded density dependent harvesting.

  * •

Harvesting females and harvesting males - FHMH models. Herein again we
introduce three subclasses, linear harvesting, saturating density dependent
harvesting and unbounded density dependent harvesting.

Both classes of models can yield extinction - and are generally more effective
than the TYC strategy. To this end we compare them via a rigorous optimal
control approach, where our metrics of judging how good or bad a model is, by
looking at various criteria such as (1) costs of putting in the controls (2)
time it takes to eradicate females (3) time it takes to eradicate males. Also
note, we use population parameters that are actual outputs from mesocosm
experiments at the USGS laboratory in Gainesville, Florida. To the best of our
knowledge this is the first study in the literature which attempts to
determine the efficay of the TYC strategy, where best fit parameters from
population experiments have been used to tune the models. Based on these
criteria we summarize some of our key results.

Fig. 6 indicates that extinction occurs with our optimal controls of Model 4.
Notice that a large amount of harvesting must first occur to ensure extinction
occurs. This is reasonable to expect since in order for the population to be
attracted to the extinction state then the populations must be driven through
substantial harvesting to low enough densities. Model 2 and Model 5, see Fig.
7, are not able to achieve eradication (for the parameters considered),
although Model 5 female and male fish numbers are lower at the terminal time.
Hence, this indicates that if harvesting is modeled by rational functions,
that is at some saturating rate, then using this form of harvesting is
disadvantageous. Furthermore, the total costs of Model 2 and 5 are _185660_
and _93600_ , respectively, indicating the hybrid harvesting strategy FHMS,
which mirrors the TYC strategy is less effective as an eradication and control
tool than its counterpart FHMH models, for such saturating forms of
harvesting.

A similar situation is observed in comparing Model 6 to Model 3. Fig. 8
indicates that extinction can occur when optimal harvesting controls are used.
Again, the total costs of Models 3 and 6 are _31050_ and _15790_ ,
respectively which indicates that FHMH model is quantifiably better than FHMS
model, even when we have unbounded harvesting functions. Fig. 9 implies that
extinction cannot occur if harvesting obeys the saturating form found in Model
2. However, females and males can be driven to extinction with our optimal
controls for Model 3, depicted in Fig. 9, where we use unbounded harvesting
functions. This is further evidence that harvesting via a saturating
harvesting term such as in Model 2 indicates unfavorable results. However, if
harvesting is modeled by a power function then favorable results are found. In
fact, the total cost of _31050_ is lower for Model 3 than that of Model 1 -
and so we can conclude that within the class of FHMS models, Model 3 is the
best; the total cost of _14250_ is lower for Model 4 than that of Model 6 -
and so we can conclude that within the class of FHMH models, Model 4 is the
best. Also we can probably say the class of FHMH models is better than the
class of FHMS models.

There is one _very interesting_ exception to this. If we compare Model 6 to
Model 3, we see that although Model 6 has a lower cost, females approximately
eradicate at the $2^{nd}$ month in the Model 3, whereas it takes about 9
months to make females extinct in the Model 6. Thus if we are concerned with
female eradication speed, for the realistic parameters that we have gotten out
of our population experiments, Model 3 actually does better than Model 6.

Note, in various situations the efficacy of harvesting can be indeed
questioned. For example, there is a large body of literature that points to
various changes in population parameters as well as body size in fish, under
pressure of harvesting. See [41], where it is observed that Lake trout
populations exhibit age reduction at first maturity and clear increases in
fecundity when harvested. Also, lake trout body size increased in populations
where exploitation caused density reduction [41]. Thus in reality, there is
strong evidence that fishes will exhibit a density-dependent response to
harvest, as seen in these examples - questioning the efficacy of harvest,
unless models show otherwise. Clearly their life history parameters change, as
their density is affected by harvest, resulting in say increased fecundity.
One then must consider density dependent life history parameters, such as
$\delta,\beta$, when we consider our harvesting class of models. To this end
we have begun mesocosm experiments to try and derive the right forms of
density dependence under harvesting. Modeling harvesting processes via such
density dependent population parameters is our next immediate future goal. The
use of harvesting is restricted to certain field situations where the target
organisms can be encountered, detected and removed in a practical manner.
Thus, harvesting may not be possible in some situations, such as in vast
landscapes that are difficult or impossible to access or when working with
organisms that are cryptic or difficult to locate. Thus, we acknowledge that
while it may not be possible in all habitats or with all species, in cases
where harvesting is feasible it may be an efficient strategy to reduce or
eliminate non-native species populations. Furthermore there is literature on
selective harvesting, such as in harvesting females versus males, and its
population consequences [33]. Another question is what our results predict
about the use of the TYC strategy as an eradication tool. This is very
subjective and depends on the resources available to the manager. If one looks
at the average rate of harvesting, in our optimally-controlled scenarios it
seems the average rate of harvest can vary between $0.03$ for Model 1 (only
harvest females) to $0.125$ for Model 4 (harvest both females and males).
Note, the natural death rate of the fish ($\delta$) is $0.065$. Thus it really
depends if one _can_ harvest at these high rates - and what the associated
costs are to sustain such a high harvesting rate - again calling into question
the efficacy of the harvesting strategy. For example, Model 4 clearly seems
the best model in terms of cost and eradication time, but what we don’t
account for here is the cost of harvesting at approx. $12\%$, when the natural
death rate is about $6\%$ \- and how feasible this is in the first place. A
viable alternative, depending again on what is feasible from a management
point of view, could be to try various _hybrid_ models, that invoke TYC type
dynamics coupled with a certain amount of harvesting. Comparing such hybrid
models to pure TYC or pure harvesting strategies, are among some of our future
directions.

In essence these harvesting models are a means of mimicking the effect of the
TYC strategy, without the use of the YY super males. The TYC strategy remains
a powerful method for invasive species control [17, 18, 19, 34, 35]. Our
future goal is to use these results to guide further experiments into
harvesting strategies. In a closed mesocosm environment this is easily doable
as we can always count total populations of the fish periodically. These
counts can then be used to harvest/restock exactly those many fish based on
the form of harvesting to be used. If favorable results are obtained in the
mesocosm, this becomes a confidence booster for natural resource managers to
try TYC type strategies coupled with some harvesting, or hybrid strategies for
the control of non-native species in the wild. Our results also establish a
framework where via harvesting in mesocosm settings - we could actually
predict the efficacy of the powerful TYC strategy in the wild.

##  7\. Figure Legends

Figure 1:
Left panel: the female/male density change with time without TYC strategy.
Right Panel: the female/male density change with time by introducing YY
supermales (TYC strategy).

Figure 2:
Male (top) and female (bottom) fancy guppy (Poecilia reticulata). Photos were
taken by Howard Jelk, USGS.

Figure 3:
Left panel: This is a tank from the USGS facility where the mesocosm
experiment was run with fancy guppies (Poecilia reticulata). Right panel:
preliminary population data over 11-month period. The smooth curve is a result
of simulating the basic population model (1) (with $s=0$) with best-fit
parameters $\beta=0.0057,~{}\delta=0.0648,$ and $K=405$.

Figure 4:
(A) Female (red), male (green) and supermale (blue) densities and optimal
control $\mu(t)$ in (1) change with time $t$. (B) Female (red) and male
(green) densities and optimal controls $\eta_{1}(t)$ and $\eta_{2}(t)$ in (10)
change with time $t$.

Figure 5:
(A) In this simulation we look at the objective function towards the optimal
implement for each model, we clearly see the cost function of TYC model (Model
0) is much larger than FHMS model (Model 1). (B) Since the constraints and
scales for $\mu$ and $\eta_{1},\eta_{2}$ are different, it’s hard to compare
the effectiveness for each strategy, so here we also look at the cost
deducting the controls, that is, $\int^{T}_{0}f+m$ $dt$ under optimal control.

Figure 6:
(A) Female (red) and male (green) densities change with time $t$ with optimal
control $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) in Model 4. As we can
see the controls and densities are indistinguishable. (B) Objective function
$J_{4}$ decrease as increasing time with optimal control.

Figure 7:
(A) Female (red) and male (blue) densities change with time $t$ with optimal
control $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) in Model 5. (B)
Objective function $J_{5}$ varies as increasing time with optimal control.

Figure 8:
(A) Female (red) and male (green) densities over time $t$ with optimal control
$\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) in Model 6. We can see that the
densities and controls are indistinguishable. (B) Objective function $J_{6}$
decrease as increasing time with optimal control.

Figure 9:
Female (red) and male (green) densities change with time $t$ and the optimal
controls $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) for (A) Model 2 and
(B) Model 3. This figure can be found in appendix.

Figure 10:
Objective function $J_{2}$ and $J_{3}$ change with time $t$ with optimal
controls $\eta_{1}^{*}$ and $\eta_{2}^{*}$ for (A) Model 2 and (B) Model 3.
This figure is in appendix.

##  8\. Acknowledgements

JL and RP would like to acknowledge valuable support from the NSF via
DMS-1715377 and DMS-1839993. MB would like to acknowledge valuable support
from the NSF via DMS-1715044. U.S. Geological Survey researchers were
supported by the Greater Everglades Priority Ecosystem Sciences program. Any
use of trade, firm or product names is for descriptive purposes only and does
not imply endorsement by the U.S. Government.

## References

  * [1] M. A. Lewis, S. V. Petrovskii and J. R. Potts.  _The mathematics behind biological invasions_.  Springer, vol 44, 2016.
  * [2] D. Pimentel, R. Zuniga and D. Morrison.  _Update on the environmental and economic costs associated with alien-invasive species in the united states_.  Ecological Economics, vol 52, no.3, pp. 273-288, 2005.
  * [3] D. Simberloff.  _Introduced species: the threat to biodiversity and what can be done_.  http://www.actionbioscience.org/biodiversity/simberloff.html .
  * [4] J. E. Hill and C. E. Cichra.  _Eradication of a reproducing population of convict cichlids cichlasoma nigrofasciatum (cichlidae) in north-central florida_.  Florida Scientist, vol 68, pp. 65-74, 2005.
  * [5] H. F. Weinberger.  _Diffusion-induced blow up in a system with equal diffusions_.  Journal of Differential Equations, vol 154, pp. 225-237, 1999.
  * [6] M. Arim, S. Abades, P. Neill, M. Lima and P. Marquet.  _Spread Dynamics of invasive species_.  Proceedings of the National Academy of Sciences, vol 103, no.2, pp 374-378, 2006.
  * [7] I. Averill and Y. Lou.  _On several conjectures from evolution of dispersal_.  , Journal of Biological Dynamics, vol.6, no.2, pp 117-130, 2012.
  * [8] C. J. Bampfylde and M. A. Lewis.  _Biological control through intraguild predation:case studies in pest control, invasive species and range expansion_.  Bulletin of Mathematical Biology, vol 69, pp 1031-1066, 2007.
  * [9] J. S. Clark, M. Lewis and L. Horvath.  _Invasion by extremes: Population spread with variation in dispersal and reproduction_.  The American Naturalist, vol 157, no.5, 2001.
  * [10] Y. Lou and D. Munther.  _Dynamics of a three species competition model_.  Discrete $\&$ Continuous Dynamical Systems-A, vol.32, pp 3099-3131, 2012.
  * [11] J. H. Myers, D. Simberloff, A. M. Kuris and J. R. Carey.  _Eradication revisited: dealing with exotic species_.  Trends in Ecology $\&$ Evolution, vol.15, pp. 316-320, 2000.
  * [12] A. Okubo, P. K. Maini, M. H. Williamson and J. D. Murray.  _The spread of the grey squirrel in Britain_.  Proceedings of the Royal Society of London, Series B, vol.238, pp. 113-125.
  * [13] N. Shigesada and K. Kawasaki.  _“Biological invasions:Theory and practice”._ Oxford University Press, Oxford, 1997.
  * [14] R. Van Driesche and T. Bellows.  “Biological Control”,  Kluwer Academic Publishers, Massachusetts, 1996.
  * [15] P. Schofield and W. Loftus.  _Nonnative fishes in Florida freshwaters: a literature review and synthesis_.  Reviews in Fish Biology and Fisheries, vol. 25, no.1, pp. 117-145.
  * [16] J. Hood.  _Asian carp: State’s fish kill in Chicago Sanitary and Ship Canal yield only 1 Asian carp_ ,  www.articles.chicagotribune.com.
  * [17] J. B. Gutierrez and J. Teem.  _A model describing the effect of sex-reversed YY fish in an established wild population: the use of a Trojan Y-Chromosome to cause extinction of an introduced exotic species_.  Journal of Theoretical Biology, vol 241, no.22, pp 333-341, 2006.
  * [18] J. L. Teem, J. B. Gutierrez and R. D. Parshad.  _A comparison of the Trojan Y Chromosome model and daughterless carp eradication strategies_.  Biological Invasions, DOI 10.1007/s10530-013-0475-2, Published online May 2013.
  * [19] D. J. Schill, K. A. Meyer and M. J. Hansen.  _Simulated effects of YY-male stocking and manual suppression for eradicating nonnative brook trout populations_.  North American Journal of Fisheries Management, 37:1054–1066, 2017.
  * [20] S. Cotton and C. Wedekind.  _Control of introduced species using Trojan sex chromosomes_.  Trends in Ecology $\&$ Evolution, vol 22, pp 441-443, 2007.
  * [21] S. Cotton and C. Wedekind.  _Introduction of Trojan sex chromosomes to boost population growth_.  Journal of Theoretical Biology, vol 249, pp 153-161, 2007.
  * [22] S. Cotton and C. Wedekind.  _Population consequences of environmental sex reversal_.  Conservation Biology, vol 23, pp 196-206, 2009.
  * [23] J. B. Gutierrez, M. K. Hurdal, R. D. Parshad and J. L. Teem.  _Analysis of the Trojan Y Chromosome model for eradication of invasive species in a dendritic riverine system_.  Journal of Mathematical Biology, vol 64, no.1-2, pp 319-340, 2012.
  * [24] R. D. Parshad.  _On the long time behavior of a PDE model for invasive species control_.  International Journal of Mathematical Analysis, vol 5, no.40, pp 1991-2015, 2011.
  * [25] R. D. Parshad and J. B. Gutierrez.  _On the well posedness of the Trojan Y Chromosome model_.  Boundary Value Problems, Volume 2010, Article ID 405816, pp 1-29, 2010.
  * [26] R. D. Parshad and J. B. Gutierrez.  _On the global attractor of the Trojan Y-Chromosome model_.  Communications on Pure and Applied Analysis, vol 10, pp 339-359, 2011.
  * [27] R. D. Parshad, S. Kouachi and J. B. Gutierrez.  _Global existence and asymptotic behavior of a model for biological control of invasive species via supermale introduction_.  Communications in Mathematical Sciences, vol 11, no.4, pp 951-972, 2013.
  * [28] R.D. Parshad, M. Beauregard, E. Takyi, T. Griffin and L. Bobo,  _Large and Small Data Blow-Up Solutions in the Trojan Y Chromosome Model_ ,  Journal of Mathematical Biology, In Reveiw 2019, arXiv:1907.06079.
  * [29] N. Perrin.  _Sex reversal: a fountain of youth for sex chromosomes_?  Evolution, vol 63, pp. 3043-3049, 2009.
  * [30] X. Wang, R. D. Parshad and J. Walton.  _The Stochastic Trojan Y Chromosome model for eradication of an invasive species_.  Journal of Biological Dynamics, vol 10, issue 1, pp 179-199, 2016.
  * [31] X. Wang, J. R. Walton, R. D. Parshad, K. Storey and M. Boggess.  _Analysis of Trojan Y-Chromosome eradication strategy_.  Journal of Mathematical Biology, vol 68, issue 7, pp. 1731-1756, 2014.
  * [32] X. Zhao, B. Liu and N. Duan.  _Existence of global attractor for the Trojan Y Chromosome model_.  Electronic Journal of Qualitative Theory of Differential Equations, No.36, 16, 2012.
  * [33] J. M. Milner, E.B. Nilsen and H.P. Andreassen,  _Demographic side effects of selective hunting in ungulates and carnivores_.  Conservation biology, 21(1), 36-47, 2007.
  * [34] D. J. Schill, J. A. Heindel, M. R. Campbell, K. A. Meyer and E. R. Mamer.  _Production of a YY Male Brook Trout Broodstock for Potential Eradication of Undesired Brook Trout Populations_.  North American Journal of Aquaculture, 78(1), 72-83, 2016.
  * [35] P. A. Kennedy, K. A. Meyer, D. J. Schill, M. R. Campbell and N.V. Vue.  _Survival and reproduction success of hatchery YY male brook trout stocked in Idaho streams_.  Transactions of the American Fisheries Society, vol. 147, pp. 419-430, 2018.
  * [36] D. Schill, personal communication, November $4^{th}$, 2016.
  * [37] S. J. Green, N. K. Dulvy, A. M. Brooks, J. L. Akins, A. B. Cooper, S. Miller and I. M. Côté.  _Linking removal targets to the ecological effects of invaders: a predictive model and fieldtest._ Ecological Applications, 24(6), pp.1311-1322, 2014.
  * [38] J. P. Parkes, G. Nugent and B. Warburton  _Commercial exploitation as a pest control tool for introduced mammals in New Zealand._ Wildlife Biology 2: 171-177, 1996.
  * [39] L. Nico, M. Neilson and B. Loftus.  _Poecilia reticulata Peters, 1859: U.S. Geological Survey, Nonindigenous Aquatic Species Database, Gainesville, FL_.  https://nas.er.usgs.gov/queries/FactSheet.aspx?SpeciesID=863, Access Date: 7/30/2019
  * [40] W. H. Fleming and R. W. Rishel.  _Deterministic and Stochastic Optimal Control_.  Springer Verlag, New York, 1975.
  * [41] Syslo, J. M., Guy, C. S., Bigelow, P. E., Doepke, P. D., Ertel, B. D. and Koel, T. M.  _Response of non-native lake trout (Salvelinus namaycush) to 15 years of harvest in Yellowstone Lake, Yellowstone National Park_.  Canadian Journal of Fisheries and Aquatic Sciences, 68(12), 2132-2145, 2011.

##  9\. Appendix A: Stability Analysis of classical TYC model

###  9.1. Equilibrium and Stability Analysis when $\boldsymbol{\mu>0}$

Let $(f^{*},m^{*},s^{*})$ represent an equilibrium of model (1). We are
interested in the region $\Omega=\\{(f,m,s)~{}|~{}0\leq f\leq K,~{}0\leq m\leq
K,~{}0\leq s\leq K,~{}0\leq f+m+s\leq K\\}$. The dynamics are investigated for
the model with a positive introduction rate, $\mu>0$.

####  9.1.1. Equilibrium analysis

An equilibrium with $f^{*}=0$ denotes a steady state of model (1) at which the
wild-type XX females are eradicated. There is only one equilibrium that
satisfies $f^{*}=0$, that is, $\left(0,0,\dfrac{\mu}{\delta}\right)$.

In contrast, assume $f^{*}>0$ then by direct calculation the equilibrium
satisfies

$\begin{split}&f^{*}=K-m^{*}-s^{*}-\frac{2K\delta}{m^{*}\beta},\\\
&s^{*}=\frac{\mu}{\delta},\end{split}$

where $m^{*}$ satisfies

(19) $am^{3}+bm^{2}+cm+d=0.$

The constants $a,~{}b,~{}c$ and $d$ are given by

$\begin{split}&a=a(s^{*})=2\beta,\\\ &b=b(s^{*})=3\beta s^{*}-K\beta,\\\
&c=c(s^{*})=2\beta s^{*2}+2K\delta-2K\beta s^{*},\\\ &d=d(s^{*})=4K\delta
s^{*}.\end{split}$

Clearly $a,d>0$ since $\mu,~{}\beta$, and $K$ are positive. Let

$\displaystyle s_{c}^{+}$ $\displaystyle=$
$\displaystyle\frac{K}{2}+\sqrt{\frac{K^{2}}{4}-\frac{K\delta}{\beta}},$
$\displaystyle s_{c}^{-}$ $\displaystyle=$
$\displaystyle\frac{K}{2}-\sqrt{\frac{K^{2}}{4}-\frac{K\delta}{\beta}},$
$\displaystyle s_{b}$ $\displaystyle=$ $\displaystyle\frac{K}{3}.$

It is easy to check that $s_{c}^{\pm}$ and $s_{b}$ are the nonzero roots for
$c(s^{*})=0$ and $b(s^{*})=0$, respectively. Define
$\alpha=\frac{K\beta}{\delta}$. For illustrative purposes, the population data
provides the following estimates: $\alpha\approx 30.375$, $s_{c}^{+}\approx
404.925$, $s_{c}^{-}\approx.075014$, and $s_{b}\approx 135$.

######  Lemma 9.1.

  1. (a)

If $\alpha>\frac{9}{2}$ then $s_{c}^{-}<s_{b}<s_{c}^{+}$;

  2. (b)

If $4<\alpha<\frac{9}{2}$ then $s_{b}<s_{c}^{-}<s_{c}^{+}$;

  3. (c)

If $\alpha=4$ then $s_{b}<s_{c}^{-}=s_{c}^{+}$;

  4. (d)

If $\alpha=\frac{9}{2}$ then $s_{b}=s_{c}^{-}<s_{c}^{+}$.

###### Proof.

It is easy to be verified by direct calculations. ∎

Next, consider the discriminant of (19), that is,

$\Delta=b^{2}c^{2}-4ac^{3}-4b^{3}d-27a^{2}d^{2}+18abcd.$

######  Proposition 1.

Assume that $\Delta>0$.

  1. (a)

If $\alpha>\frac{9}{2}$, the (19) has 2 positive roots when
$\frac{\mu}{\delta}<s_{c}^{+}$ and no positive root when
$\frac{\mu}{\delta}\geq s_{c}^{+}.$

  2. (b)

If $s^{*}\in[s_{b},s_{c}^{-}]\cup[s_{c}^{+},\infty)$, then (19) has 2 positive
roots when $\frac{\mu}{\delta}\in(0,s_{b})\cup(s_{c}^{-},s_{c}^{+})$ and no
positive root when
$\frac{\mu}{\delta}\in[s_{b},s_{c}^{-}]\cup[s_{c}^{+},\infty)$.

###### Proof.

Define

(20) $h(m)=am^{3}+bm^{2}+cm+d.$

Since $a>0,d>0$ and $\Delta>0$, the cubic function (20) has 3 distinct nonzero
real roots.

  1. (a)

If $\alpha>\frac{9}{2}$, by lemma (9.1), we have $s_{c}^{-}<s_{b}<s_{c}^{+}$.

    1. (i)

Assume $s^{*}=\frac{\mu}{\delta}\in(0,s_{c}^{-})$. Then $b<0$ and $c>0$. The
signs of coefficients of $h(m)$ and $h(-m)$ are

$\displaystyle\mbox{sign}(a,b,c,d)$ $\displaystyle=$ $\displaystyle(+,-,+,+),$
$\displaystyle\mbox{sign}(-a,b,-c,d)$ $\displaystyle=$
$\displaystyle(-,-,-,+).$

Thus, two sign changes in $h(m)$ and one sign changes in $h(-m)$. By
Descartes’ rule of signs, $h(m)=0$ has either 2 or 0 positive roots and
exactly 1 negative root. Since $h(m)=0$ has three distinct nonzero real roots,
the roots of (20) has to be 2 positive roots and 1 negative root.

    2. (ii)

Assume $s^{*}=\frac{\mu}{\delta}\in[s_{c}^{-},s_{b})$. If $s^{*}=s_{c}^{-}$,
then $c=0,b<0$; if $s^{*}\in(s_{c}^{-},s_{b})$, then $c<0,b<0$. By Descartes’
rule of signs, it implies $h(m)$ has 2 positive roots and 1 negative one.

    3. (iii)

Assume $s^{*}=\frac{\mu}{\delta}\in[s_{b},s_{c}^{+})$. If $s^{*}=s_{b}$, then
$b=0,c<0$; if $s^{*}\in(s_{b},s_{c}^{+})$, then $b>0,c<0$. By a similar
argument, $h(m)$ has 2 positive roots and 1 negative root.

    4. (iv)

Assume $s^{*}=\frac{\mu}{\delta}\in[s_{c}^{+},\infty)$. Consequently, $b,c>0$.
We have no positive roots for $h(m)$ because there is no sign change in the
coefficient of $h(m)$.

  2. (b)

If $\alpha\leq\frac{9}{2}$, we always have $s_{b}\leq s_{c}^{-}$. Similarly,

    1. (i)

If $s^{*}\in(0,s_{b})\cup(s_{c}^{-},s_{c}^{+})$, then either $b<0,c>0$ or
$b>0,c<0$. There are two sign changes for the coefficient of $h(m)$, therefore
$h(m)$ has 2 positive roots.

    2. (ii)

If $s^{*}\in[s_{b},s_{c}^{-}]\cup[s_{c}^{+},\infty)$, $b\geq 0,c\geq 0$, there
is no sign change for the coefficient of $h(m)$ and therefore $h(m)$ has no
positive root.

∎

######  Proposition 2.

Assume that $\Delta=0$.

  1. (a)

If $\alpha>\frac{9}{2}$, then (19) has 2 positive roots when
$\frac{\mu}{\delta}<s_{c}^{+}$ and no positive root when
$\frac{\mu}{\delta}\geq s_{c}^{+}.$

  2. (b)

If $s^{*}\in[s_{b},s_{c}^{-}]\cup[s_{c}^{+},\infty)$, then (19) has 2 positive
roots when $\frac{\mu}{\delta}\in(0,s_{b})\cup(s_{c}^{-},s_{c}^{+})$ and no
positive root when
$\frac{\mu}{\delta}\in[s_{b},s_{c}^{-}]\cup[s_{c}^{+},\infty)$.

###### Proof.

If $\Delta=0,h(m)=0$ has a repeated root and all its roots are real. The
results follow by Descartes’ rule of signs. ∎

######  Proposition 3.

Assume that $\Delta<0$. The equation (19) has no positive root.

###### Proof.

If $\Delta=0,$ then $h(m)=0$ has one real root and two non-real complex
conjugate roots. The results follow by Descartes’ rule of signs. ∎

####  9.1.2. Stability analysis of equilibria

The associated Jacobian matrix of the model (1) is given by

(21) $\mathrm{J}=\left[\begin{array}[]{ccc}J_{11}&J_{12}&J_{13}\\\
J_{21}&J_{22}&J_{23}\\\ 0&0&J_{33}\end{array}\right]$

where

$\displaystyle J_{11}$ $\displaystyle=$ $\displaystyle\frac{1}{2}m^{*}\beta
L^{*}-\frac{1}{2}\frac{f^{*}m^{*}\beta}{K}-\delta,$ $\displaystyle J_{12}$
$\displaystyle=$ $\displaystyle\frac{1}{2}f^{*}\beta
L^{*}-\frac{1}{2}\frac{f^{*}m^{*}\beta}{K},$ $\displaystyle J_{13}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\frac{f^{*}m^{*}\beta}{K},$
$\displaystyle J_{21}$ $\displaystyle=$
$\displaystyle(\frac{1}{2}m^{*}+s^{*})\beta
L^{*}-\frac{(\frac{1}{2}f^{*}m^{*}+f^{*}s^{*})\beta}{K},$ $\displaystyle
J_{22}$ $\displaystyle=$ $\displaystyle\frac{1}{2}f^{*}\beta
L^{*}-\frac{(\frac{1}{2}f^{*}m^{*}+f^{*}s^{*})\beta}{K}-\delta,$
$\displaystyle J_{23}$ $\displaystyle=$ $\displaystyle f^{*}\beta
L^{*}-\frac{(\frac{1}{2}f^{*}m^{*}+f^{*}s^{*})\beta}{K},$ $\displaystyle
J_{33}$ $\displaystyle=$ $\displaystyle-\delta,$

and

$L^{*}=1-\frac{f^{*}+m^{*}+s^{*}}{K}.$

The corresponding characteristic function is given by

$\lambda^{3}+k_{12}\lambda^{2}+k_{11}\lambda+k_{10}=0$

with

$\begin{split}&k_{12}=-(J_{11}+J_{22}+J_{33}),\\\
&k_{11}=J_{11}J_{22}+J_{11}J_{33}-J_{12}J_{21}+J_{22}J_{33},\\\
&k_{10}=J_{12}J_{21}J_{33}-J_{11}J_{22}J_{33}.\end{split}$

######  Theorem 9.2.

Let $\mu>0$. The interior equilibrium $(f^{*},m^{*},\frac{\mu}{\delta})$ where
$f^{*},m^{*}>0$ is locally stable if

(22) $k_{10}>0,~{}k_{11}>0,~{}k_{12}>0,~{}k_{11}k_{12}-k_{10}>0.$

###### Proof.

It follows from Routh Hurwitz stability criteria. ∎

##  10\. Appendix B: Nonlinear Forms of Harvesting for FHMS Models

In the previous two sections linear harvesting functions considered only.
Here, we compare the previous results to nonlinear harvesting functions. That
is, we examine the control on the system

(23) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}G_{1}(f),\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m+\eta_{2}G_{2}(m),\end{split}$

where $G_{1}(f)$ and $G_{2}(m)$ are non-negative, possibly nonlinear,
functions such that $\\{f,m\\}\in\\{0\leq f,m\leq K$ and $0\leq f+m\leq K\\}$.
In the previous section $G_{1}(f)=f$ and $G_{2}(m)=m$; that model will be
called Model 1 in the forthcoming discussion.

###  10.1. Model 2: Nonlinear Harvesting when
$\boldsymbol{G_{1}(f)=\frac{f}{f+d_{1}},~{}G_{2}(m)=\frac{m}{m+d_{2}}}$

Consider the model

(24) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}\frac{f}{f+d_{1}},\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m+\eta_{2}\frac{m}{m+d_{2}},\end{split}$

where $d_{1}>0,d_{2}>0$. The nonlinear harvesting function are rational
functions and attempt to model a saturation of harvesting at large
populations.

######  Theorem 10.1.

The trivial equilibrium of (24) is globally asymptotically stable if $\beta
K<2\delta+\frac{\eta_{1}}{K+d_{1}}-\frac{\eta_{2}}{d_{2}}$ and
$\delta>\frac{\eta_{2}}{d_{2}}$.

###### Proof.

Consider the Lyapunov function $V=f+m$. It is left to show that
$\frac{dV}{dt}<0$ for all $(f,m)\not=(0,0)$. Consider

$\begin{split}\frac{dV}{dt}&=\frac{df}{dt}+\frac{dm}{dt}\\\ &=fm\beta
L-(\delta+\frac{\eta_{1}}{f+d_{1}})f-(\delta-\frac{\eta_{2}}{m+d_{2}})m\\\
&\leq
fm\beta-(\delta+\frac{\eta_{1}}{f+d_{1}})f-(\delta-\frac{\eta_{2}}{m+d_{2}})m,\quad\mathrm{since}\quad
L\leq 1\\\ &\leq
fm\beta-(\delta+\frac{\eta_{1}}{K+d_{1}})f-(\delta-\frac{\eta_{2}}{d_{2}})m\\\
&\leq[\beta-\frac{1}{m}(\delta+\frac{\eta_{1}}{K+d_{1}})-\frac{1}{f}(\delta-\frac{\eta_{2}}{d_{2}})]fm\\\
&\leq[\beta-\frac{1}{K}(\delta+\frac{\eta_{1}}{K+d_{1}})-\frac{1}{K}(\delta-\frac{\eta_{2}}{d_{2}})]fm.\\\
\end{split}$

It is enough to show
$\beta-\frac{1}{K}(\delta+\frac{\eta_{1}}{K+d_{1}})-\frac{1}{K}(\delta-\frac{\eta_{2}}{d_{2}})<0$.
Now $\delta-\frac{\eta_{2}}{d_{2}}>0$ and, by direct calculation, $\beta
K<2\delta+\frac{\eta_{1}}{K+d_{1}}-\frac{\eta_{2}}{d_{2}}$. Subsequently,
$\frac{dV}{dt}<0.$ ∎

As in the previous section, assume $\eta_{1},\eta_{2}$ are time dependent and
set the objective function as

$J_{2}(\eta_{1},\eta_{2})=\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt,$

subject to (24) and with initial conditions $f(t_{0})=f_{0},m(t_{0})=m_{0}$.
Optimal controls are sought within the set $U_{2}$ where

$U_{2}=\\{(\eta_{1},\eta_{2})~{}|~{}\eta_{i}\ \mbox{measurable},\
0\leq\eta_{1}\leq 1,0\leq\eta_{2}\leq 1,\ t\in[0,T],\ \forall T\\}.$

The goal is to seek an optimal $(\eta_{1}^{*},\eta_{2}^{*})$ such that

$\begin{split}J_{2}(\eta_{1}^{*},\eta_{2}^{*})&=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt\\\
&=\underset{(\eta_{1},\eta_{2})}{\min}\int^{T}_{0}f+m+\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.\end{split}$

Consider the following existence theorem.

######  Theorem 10.2.

Consider the optimal control problem (24). There exists
$(\eta_{1}^{*},\eta_{2}^{*})\in U_{2}$ such that

$J_{2}(\eta_{1}^{*},\eta_{2}^{*})=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.$

###### Proof.

The proof is similar to Theorem 3.2 and is omitted for brevity. ∎

Consider the Hamiltonian for $J_{2}$,

$H_{2}=-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}.$

The necessary conditions on the optimal controls are derived by applying
Pontryagin’s maximum principle. Namely,

$\begin{split}\lambda_{1}^{\prime}(t)=&1-\lambda_{1}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\frac{\eta_{1}}{f+d_{1}}+\frac{\eta_{1}f}{(f+d_{1})^{2}}]\\\
&-\lambda_{2}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}],\\\
\lambda_{2}^{\prime}(t)=&1-\lambda_{1}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}]\\\
&-\lambda_{2}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta+\frac{\eta_{2}}{m+d_{2}}-\frac{\eta_{2}m}{(m+d_{2})^{2}}],\end{split}$

with the transversality condition given as $\lambda_{1}(T)=\lambda_{2}(T)=0$.
Differentiating with respect to our controls yields,

$\begin{split}&\frac{\partial
H_{2}}{\partial\eta_{1}}=-\frac{f\lambda_{1}}{f+d_{1}}-\eta_{1},\\\
&\frac{\partial
H_{2}}{\partial\eta_{2}}=\frac{m\lambda_{2}}{m+d_{2}}-\eta_{2}.\end{split}$

Hence, the optimal $\eta_{1}^{*},$ is characterized by the three cases:

  1. (1)

If $\frac{\partial H_{2}}{\partial\eta_{1}}<0$ then $\eta_{1}^{*}=0$. This
implies that $-\frac{f\lambda_{1}}{f+d_{1}}<0;$

  2. (2)

If $\frac{\partial H_{2}}{\partial\eta_{1}}=0$ then
$\eta_{1}^{*}=-\frac{f\lambda_{1}}{f+d_{1}}$. This implies that
$0\leq-\frac{f\lambda_{1}}{f+d_{1}}\leq 1$;

  3. (3)

If $\frac{\partial H_{2}}{\partial\eta_{1}}>0$ then $\eta_{1}^{*}=1$. This
implies that $1<-\frac{f\lambda_{1}}{f+d_{1}}.$

Similar results are established for the optimal $\eta_{2}^{*}$. A compact way
of writing the optimal control $\eta_{1},\eta_{2}$ is

(25) $\displaystyle\eta_{1}^{*}$ $\displaystyle=$
$\displaystyle\min\left(1,\max\left(0,-\frac{f\lambda_{1}}{f+d_{1}}\right)\right)$
(26) $\displaystyle\eta_{2}^{*}$ $\displaystyle=$
$\displaystyle\min\left(1,\max\left(0,\frac{m\lambda_{2}}{m+d_{2}}\right)\right).$

######  Theorem 10.3.

An optimal control $(\eta_{1}^{*},\eta_{2}^{*})\in U_{2}$ for the system (24)
that maximizes the objective functional $J_{2}$ is characterized by (25)-(26).

###  10.2. Model 3: Nonlinear Harvesting when
$\boldsymbol{G_{1}(f)=f^{\frac{3}{2}},~{}G_{2}(m)=m^{\frac{3}{2}}}$

Consider the the model

(27) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}f^{\frac{3}{2}},\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m+\eta_{2}m^{\frac{3}{2}}.\end{split}$

The nonlinear harvesting functions attempt to model a sharp increase in the
ability to harvest fish at low populations. Of course, in contrast to the
previous Model 2, for large populations the harvesting terms has no upper
bound. Such terms could express chemical use that are density dependent [15].
In order to compare Models 2 and 3 see Fig. 10 \- Fig. 10.

######  Theorem 10.4.

The trivial equilibrium is globally asymptotically stable in (27) if $\beta
K<2\delta-\eta_{2}\sqrt{K}$ and $\delta>\eta_{2}\sqrt{K}$.

###### Proof.

Consider the Lyapunov function $V=f+m$. It remains to show that
$\frac{dV}{dt}<0$ for all $(f,m)\not=(0,0)$. Consider

$\begin{split}\frac{dV}{dt}&=\frac{df}{dt}+\frac{dm}{dt}\\\ &=fm\beta
L-(\delta+\eta_{1}\sqrt{f})f-(\delta-\eta_{2}\sqrt{m})m\\\ &\leq
fm\beta-(\delta+\eta_{1}\sqrt{f})f-(\delta-\eta_{2}\sqrt{m})m,\quad\mathrm{since}\quad
L\leq 1\\\ &\leq fm\beta-\delta f-(\delta-\eta_{2}\sqrt{K})m\\\
&\leq(\beta-\frac{\delta}{m}-\frac{(\delta-\eta_{2}\sqrt{K})}{f})fm\\\
&\leq(\beta-\frac{\delta}{K}-\frac{(\delta-\eta_{2}\sqrt{K})}{K})fm.\\\
\end{split}$

It is enough to show
$\beta-\frac{\delta}{K}-\frac{(\delta-\eta_{2}\sqrt{K})}{K}<0$ and
$\delta-\eta_{2}\sqrt{K}>0$. By direct calculation, we have $\beta
K<2\delta-\eta_{2}\sqrt{K}$ and $\delta>\eta_{2}\sqrt{K}$, which complete the
proof. ∎

As in the previous section, we assume $\eta_{1}$ and $\eta_{2}$ are time
dependent and set the objective function

$J_{3}(\eta_{1},\eta_{2})=\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2}).dt$

subject to (27) and with initial conditions $f(t_{0})=f_{0},m(t_{0})=m_{0}$.
The optimal controls are searched for within the set $U_{3}$ where

$U_{3}=\\{(\eta_{1},\eta_{2})~{}|~{}\eta_{i}\ \mbox{measurable},\
0\leq\eta_{1}\leq 1,0\leq\eta_{2}\leq 1,\ t\in[0,T],\ \forall T\\}.$

The goal is to seek an optimal $(\eta_{1}^{*},\eta_{2}^{*})$ such that

$\begin{split}J_{3}(\eta_{1}^{*},\eta_{2}^{*})&=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt\\\
&=\underset{(\eta_{1},\eta_{2})}{\min}\int^{T}_{0}f+m+\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.\\\
\end{split}$

Again, we prove the following existence theorem.

######  Theorem 10.5.

Consider the optimal control problem (27). There exists
$(\eta_{1}^{*},\eta_{2}^{*})\in U_{3}$ s.t.

$J_{3}(\eta_{1}^{*},\eta_{2}^{*})=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.$

###### Proof.

The proof is similar to Theorem 3.2 and is omitted for brevity. ∎

The Hamiltonian for $J_{3}$ is

$H_{3}=-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}.$

The necessary conditions on the optimal controls are determined by applying
Pontryagin’s maximum principle, that is,

$\begin{split}\lambda_{1}^{\prime}(t)=&1-\lambda_{1}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\frac{3}{2}\eta_{1}f^{\frac{1}{2}}]\\\
&-\lambda_{2}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}],\\\
\lambda_{2}^{\prime}(t)=&1-\lambda_{1}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}]\\\
&-\lambda_{2}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta+\frac{3}{2}\eta_{2}m^{\frac{1}{2}}],\end{split}$

with the transversality condition of $\lambda_{1}(T)=\lambda_{2}(T)=0$.
Differentiating the Hamiltonian with respect to the control variables yields,

$\begin{split}&\frac{\partial
H_{3}}{\partial\eta_{1}}=-\lambda_{1}f^{\frac{3}{2}}-\eta_{1},\\\
&\frac{\partial
H_{3}}{\partial\eta_{2}}=\lambda_{2}m^{\frac{3}{2}}-\eta_{2}.\end{split}$

Hence, the optimal control for $\eta_{1}^{*}$ and $\eta_{2}^{*}$ is written
compactly as

(28) $\displaystyle\eta_{1}^{*}$ $\displaystyle=$
$\displaystyle\min(1,\max(0,-\lambda_{1}f^{\frac{3}{2}})),$ (29)
$\displaystyle\eta_{2}^{*}$ $\displaystyle=$
$\displaystyle\min(1,\max(0,~{}\lambda_{2}m^{\frac{3}{2}})).$

######  Theorem 10.6.

An optimal control $(\eta_{1}^{*},\eta_{2}^{*})\in U_{3}$ for the system (27)
that maximizes the objective functional $J_{3}$ is characterized by (28)-(29).

###  10.3. Numerical Simulations and Comparisons

In this section, we will numerically simulate the optimal strategy and its
corresponding optimal states for each model. Table (5) details the parameters
used in the simulations.

Table 5. Parameters used for numerical simulations Parameter | Description | Value
---|---|---
$\beta$ | Birth rate | 0.005774
$\delta$ | Death rate | 0.0648
$K$ | Carrying capacity | 405.0705
$T$ | Terminal time | 200
$d_{1}$ | parameter in $G_{1}$ | 1
$d_{2}$ | parameter in $G_{2}$ | 1

(a) Model 2

(b) Model 3

Figure 9. Female (red) and male (green) densities change with time $t$ and the
optimal controls $\eta_{1}^{*}$ (blue) and $\eta_{2}^{*}$ (red) for (A) Model
2 and (B) Model 3.

(a) Model 2

(b) Model 3

Figure 10. Objective function $J_{2}$ and $J_{3}$ change with time $t$ with
optimal controls $\eta_{1}^{*}$ and $\eta_{2}^{*}$ for (A) Model 2 and (B)
Model 3.

##  11\. Appendix C: Nonlinear forms of Harvesting for FHMH Models

###  11.1. Model 5:
$\boldsymbol{G_{1}(f)=\frac{f}{f+d_{1}},G_{2}(m)=\frac{m}{m+d_{2}}}$

The model is given by

(30) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}\frac{f}{f+d_{1}},\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m-\eta_{2}\frac{m}{m+d_{2}},\end{split}$

where $d_{1}>0,d_{2}>0$. Again, the harvesting functions attempt to model a
saturation of harvesting for large populations.

####  11.1.1. Global stability

The following theorem provides the criteria for global stability of the
trivial equilibrium.

######  Theorem 11.1.

The trivial equilibrium is globally asymptotically stable in (30) if $\beta
K<2\delta+\frac{\eta_{1}}{K+d_{1}}+\frac{\eta_{2}}{K+d_{2}}$.

###### Proof.

Consider the Lyapunov function $V=f+m$. It is left to show that
$\frac{dV}{dt}<0$ for all $(f,m)\not=(0,0)$. Consider

$\begin{split}\frac{dV}{dt}&=\frac{df}{dt}+\frac{dm}{dt}\\\ &=fm\beta
L-(\delta+\frac{\eta_{1}}{f+d_{1}})f-(\delta+\frac{\eta_{2}}{m+d_{2}})m\\\
&\leq
fm\beta-(\delta+\frac{\eta_{1}}{f+d_{1}})f-(\delta+\frac{\eta_{2}}{m+d_{2}})m,\quad\mathrm{since}\quad
L\leq 1\\\
&=[\beta-\frac{1}{m}(\delta+\frac{\eta_{1}}{f+d_{1}})-\frac{1}{f}(\delta+\frac{\eta_{2}}{m+d_{2}})]fm\\\
&\leq[\beta-\frac{1}{m}(\delta+\frac{\eta_{1}}{K+d_{1}})-\frac{1}{f}(\delta+\frac{\eta_{2}}{K+d_{2}})]fm\\\
&\leq[\beta-\frac{1}{K}(\delta+\frac{\eta_{1}}{K+d_{1}})-\frac{1}{K}(\delta+\frac{\eta_{2}}{K+d_{2}})]fm.\\\
\end{split}$

It is enough to show
$\beta-\frac{1}{K}(\delta+\frac{\eta_{1}}{K+d_{1}})-\frac{1}{K}(\delta+\frac{\eta_{2}}{K+d_{2}})<0$.
By direct calculation, we have $\beta
K<2\delta+\frac{\eta_{1}}{K+d_{1}}+\frac{\eta_{2}}{K+d_{2}}$, which proves
theorem (11.1). ∎

####  11.1.2. Optimal control

We again assume $\eta_{1},\eta_{2}$ are time dependent and set the objective
function as

$J_{5}(\eta_{1},\eta_{2})=\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.$

subject to (30) and with initial conditions $f(t_{0})=f_{0},m(t_{0})=m_{0}$.
Optimal controls are sought within $U_{5}$, that is,

$U_{5}=\\{(\eta_{1},~{}\eta_{2})~{}|~{}\eta_{i}\ \mbox{measurable},\
0\leq\eta_{1}\leq 1,0\leq\eta_{2}\leq 1,\ t\in[0,T],\ \forall T\\}.$

Hence, optimal $(\eta_{1}^{*},\eta_{2}^{*})$ are sought such that,

$\begin{split}J_{5}(\eta_{1}^{*},\eta_{2}^{*})&=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt\\\
&=\underset{(\eta_{1},\eta_{2})}{\min}\int^{T}_{0}f+m+\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.\\\
\end{split}$

The following existence theorem is given.

######  Theorem 11.2.

Consider the optimal control problem (24). There exists
$(\eta_{1}^{*},\eta_{2}^{*})\in U_{5}$ such that

$J_{5}(\eta_{1}^{*},\eta_{2}^{*})=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.$

###### Proof.

The proof is similar to Theorem 5.2 and is omitted for brevity. ∎

Consider the Hamiltonian for $J_{5}$,

(31)
$H_{5}=-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}.$

Pontryagin’s maximum principle is applied to determine necessary conditions on
the optimal controls. The differential equations for the adjoint $\lambda_{i}$
is

$\begin{split}\lambda_{1}^{\prime}(t)=&1-\lambda_{1}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\frac{\eta_{1}}{f+d_{1}}+\frac{\eta_{1}f}{(f+d_{1})^{2}}]\\\
&-\lambda_{2}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}].\\\
\lambda_{2}^{\prime}(t)=&1-\lambda_{1}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}]\\\
&-\lambda_{2}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\frac{\eta_{2}}{m+d_{2}}+\frac{\eta_{2}m}{(m+d_{2})^{2}}].\end{split}$

with the transversality condition $\lambda_{1}(T)=\lambda_{2}(T)=0.$
Differentiating the Hamiltonian,

$\begin{split}&\frac{\partial
H_{5}}{\partial\eta_{1}}=-\frac{f\lambda_{1}}{f+d_{1}}-\eta_{1},\\\
&\frac{\partial
H_{5}}{\partial\eta_{2}}=\frac{m\lambda_{2}}{m+d_{2}}-\eta_{2}.\end{split}$

Hence, the optimal control $\eta_{1},\eta_{2}$ is written compactly as,

(32) $\displaystyle\eta_{1}^{*}$ $\displaystyle=$
$\displaystyle\min\left(1,\max\left(0,-\frac{f\lambda_{1}}{f+d_{1}}\right)\right),$
(33) $\displaystyle\eta_{2}^{*}$ $\displaystyle=$
$\displaystyle\min\left(1,\max\left(0,-\frac{m\lambda_{2}}{m+d_{2}}\right)\right).$

######  Theorem 11.3.

An optimal control $(\eta_{1}^{*},\eta_{2}^{*})\in U_{5}$ for the system (30)
that maximizes the objective functional $J_{2}$ is characterized by (32)-(33).

###  11.2. Model 6:
$\boldsymbol{G_{1}(f)=f^{\frac{3}{2}},G_{2}(m)=m^{\frac{3}{2}}}$

Consider an alternative model of the harvesting given by,

(34) $\begin{split}&\frac{df}{dt}=\frac{1}{2}fm\beta L-\delta
f-\eta_{1}f^{\frac{3}{2}},\\\ &\frac{dm}{dt}=\frac{1}{2}fm\beta L-\delta
m-\eta_{2}m^{\frac{3}{2}}.\end{split}$

####  11.2.1. Global stability

The global stability of the trivial equilibrium for model (34) is proven in
the following theorem.

######  Theorem 11.4.

The trivial equilibrium $(0,0)$ for (34) is globally asymptotically stable if
$\beta K<2\delta$.

###### Proof.

Again, consider the Lyapunov function $V=f+m$. It is left to show that
$\frac{dV}{dt}<0$ for all $(f,m)\not=(0,0)$. Consider

$\begin{split}\frac{dV}{dt}&=\frac{df}{dt}+\frac{dm}{dt}\\\ &=fm\beta
L-(\delta+\eta_{1}\sqrt{f})f-(\delta+\eta_{2}\sqrt{m})m\\\ &\leq
fm\beta-(\delta+\eta_{1}\sqrt{f})f-(\delta+\eta_{2}\sqrt{m})m,\quad\mathrm{since}\quad
L\leq 1\\\ &\leq fm\beta-\delta f-\delta m\\\
&\leq(\beta-\frac{\delta}{m}-\frac{\delta}{f})fm\\\
&\leq(\beta-\frac{\delta}{K}-\frac{\delta}{K})fm.\\\ \end{split}$

It is enough to show $\beta-\frac{\delta}{K}-\frac{\delta}{K}<0$. By direct
calculation, we have $\beta K<2\delta$, which proves theorem (11.4). ∎

####  11.2.2. Optimal control

The parameters $\eta_{1}$ and $\eta_{2}$ are assumed to be time dependent and
define the objective function

$J_{6}(\eta_{1},\eta_{2})=\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt$

subject to (34) and with initial conditions $f(t_{0})=f_{0},m(t_{0})=m_{0}$.
The optimal controls are sought within the set $U_{6}$ where

$U_{6}=\\{(\eta_{1},\eta_{2})~{}|~{}\eta_{i}\ \mbox{measurable},\
0\leq\eta_{1}\leq 1,0\leq\eta_{2}\leq 1,\ t\in[0,T],\ \forall T\\}.$

The goal is obtain optimal $(\eta_{1}^{*},\eta_{2}^{*})$ such that

$\begin{split}J_{6}(\eta_{1}^{*},\eta_{2}^{*})&=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt\\\
&=\underset{(\eta_{1},\eta_{2})}{\min}\int^{T}_{0}f+m+\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.\\\
\end{split}$

The following existence theorem is given.

######  Theorem 11.5.

Consider the optimal control problem (27). There exists
$(\eta_{1}^{*},\eta_{2}^{*})\in U_{6}$ such that

(35)
$J_{6}(\eta_{1}^{*},\eta_{2}^{*})=\underset{(\eta_{1},\eta_{2})}{\max}\int^{T}_{0}-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})dt.$

###### Proof.

The proof is similar to Theorem 5.2 and is omitted for brevity. ∎

As in the previous sections, the Hamiltonian for $J_{6}$ is considered, that
is,

$H_{6}=-(f+m)-\frac{1}{2}(\eta_{1}^{2}+\eta_{2}^{2})+\lambda_{1}f^{\prime}+\lambda_{2}m^{\prime}.$

Again, we establish differential equations for the adjoint,

$\begin{split}\lambda_{1}^{\prime}(t)=&1-\lambda_{1}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\frac{3}{2}\eta_{1}f^{\frac{1}{2}}]\\\
&-\lambda_{2}[\frac{m\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}],\\\
\lambda_{2}^{\prime}(t)=&1-\lambda_{1}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}]\\\
&-\lambda_{2}[\frac{f\beta}{2}(1-\frac{f+m}{K})-\frac{fm\beta}{2K}-\delta-\frac{3}{2}\eta_{2}m^{\frac{1}{2}}],\end{split}$

with the transversality condition $\lambda_{1}(T)=\lambda_{2}(T)=0$. The
derivatives of the Hamiltonian are,

$\begin{split}&\frac{\partial
H_{6}}{\partial\eta_{1}}=-\lambda_{1}f^{\frac{3}{2}}-\eta_{1},\\\
&\frac{\partial
H_{6}}{\partial\eta_{2}}=-\lambda_{2}m^{\frac{3}{2}}-\eta_{2}.\end{split}$

It is simple to verify that the optimal controls $\eta_{1}^{*}$ and
$\eta_{2}^{*}$ are compactly given by

(36) $\displaystyle\eta_{1}^{*}$ $\displaystyle=$
$\displaystyle\min(1,\max(0,-\lambda_{1}f^{\frac{3}{2}})),$ (37)
$\displaystyle\eta_{2}^{*}$ $\displaystyle=$
$\displaystyle\min(1,\max(0,-\lambda_{2}m^{\frac{3}{2}})).$

######  Theorem 11.6.

An optimal control $(\eta_{1}^{*},\eta_{2}^{*})\in U_{6}$ for the system (34)
that maximizes the objective functional $J_{6}$ is characterized by (36)-(37).