# Bioeconomic analysis of harvesting within a predator-prey system: A case
study in the Chesapeake Bay fisheries

N. Kolling, J. Herrmann, M.B. Hallare and I.N. Panayotova

(July 2022)

###### Abstract

Sustainable use of biological resources is very important as over-exploitation
on a long run may lead to the stock depletion which in turn may threaten
biodiversity. The Chesapeake Bay is an extremely complex ecosystem, and
sustainable harvesting of its fisheries is very important both for the
ecosystem biodiversity and for the economic prosperity of the area. Here we
use mathematical modeling to study the population dynamics with harvesting of
two very important fish in the Chesapeake Bay, menhaden as a prey and striped
bass as a predator. We start by fitting the generalized Lotka-Volterra model
to real population data for the two species obtained from the fisheries in the
Bay. We derive conditions for the existence of the bio-economic equilibrium
and investigate the stability and the resilience of the biological system. We
study the maximum sustainable yield, maximum economic yield, and resilience
maximizing yield policies and their effects on the fisheries long time
sustainability and particularly on the menhaden-bass population dynamics. This
study may be used by policy-makers to balance the economic and ecological
harvesting goals while managing the populations of Atlantic menhaden and
striped bass in the Chesapeake Bay fisheries.

Keywords Mathematical Modeling, Differential Equations, Prey-Predator
Dynamics, harvesting policies, maximum sustainable yield, maximum economic
yield, resilience maximizing yield.

##  1 Introduction

The Atlantic menhaden (Brevoortia tyrannus) is widely considered to be one of
the most important species of fish in the Chesapeake Bay. Atlantic menhaden
grows up to 18 inches, can live up to around 10 to 12 years, and can weigh up
to a pound [1, 7]. They generally move around in schools and their spawning
_can_ occur year round, but will generally take place during cooler months [1,
2]. The Atlantic menhaden is a filter-feeder, and its diet consists mostly of
various types of plankton and a bit of plant materials as well [7, 16, 22].
The Atlantic menhaden is an extremely important fish in the Chesapeake Bay,
providing both biological and economic benefits. It serves as middle grounds
in many various food chains in the Bay and is a key prey species for many
other fish, such as the striped bass [42], the Bluefish [15], the Weakfish
[20], and many other species [7, 22]. Along with this, the oil produced from
the Atlantic menhaden has both economic and health benefits. It is a good
source of Omega-3 fatty acids, could be used to lower health care costs and
help with cardiovascular diseases [22].

However, the overall population of the Atlantic menhaden has been on a
significant decline since roughly 1978, and continues to remain at a low mean
[4]. Species within the Chesapeake Bay and menhaden in particular, are
seriously affected by problems such as the air pollution, the number invasive
species in the area, the wetland destruction from development, deforestation,
climate change and the rising sea levels, runoff pollution, and chemical
contaminants [6, 8, 23, 34]. But one of the biggest threat to the Atlantic
menhaden that directly affects their constantly decreasing population is
overharvesting [4]. The Omega Protein Corporation has been in the spotlight
for many years for ignoring harvest caps placed on the Atlantic menhaden [4,
26, 35, 37]. Even with all these factors that contribute to increasing the
risk of further depletion of menhaden as a species, in the ASMFC 2019 Stock
Assessment Report it was stated that ”the Atlantic menhaden species is not
overfished and overfishing is not occurring” [29].

Another very important fish within the Chesapeake Bay is the striped bass
(Morone saxatilis). The adults can range typically anywhere from 2 to 3 feet
in length and 10 to 30 lbs in weight [10, 39]. The striped bass feeds on
various species, such as the Atlantic menhaden, crustaceans, and anchovies
[10, 18, 42], and are fed on by more predatory species in the Chesapeake Bay,
such as sharks [12], ospreys [14], and other larger species of fish [10].
Striped bass spawning will generally take place from April through June, where
the females will lay their eggs typically in fresh waters. From there, the
larvae hatch after a few days, and the new young will spend a few years around
the areas they hatched in [10, 27, 39]. The most important roles of the
striped bass in the Bay is its key role in the fisheries food chain, prey to
larger predators and a key predator to smaller species, as well as its role in
supporting one of the Bay’s most popular and biggest commercial and
recreational fisheries [10]. A study conducted by Southwick Associates for the
McGraw Center for Conservation Leadership analyzed the economic contributions
of both the recreational and commercial fisheries for the striped bass [33].
With the exception of Maryland, fisheries in eleven different states along the
Atlantic coast have the majority of landings, GDP, number of jobs supported,
and salaries and wage amounts that mostly come from the recreational
fisheries. Furthermore, this study have also found the largely positive impact
the recreational fishery has on the GDP of the East Coast of the United States
[33].

Despite the importance of the striped bass, similar to the Atlantic menhaden,
the fish faces issues that put the species at great risk. The striped bass
populations have been at risk for extinction many times in the past several
decades. The bass population faced a significant decline in its numbers
throughout the 1970s and 1980s. The largest factor that was considered was
overharvesting, which put the bass population at even more risk to other
dangers, such as pollution, acidity in bass spawning areas, and low dissolved
oxygen in the bass’s habitats [9]. In the past couple of years, the bass have
yet again been facing a decline in their numbers. Overfishing is a major
factor in this cause, yet other issues have just as much of a negative impact
on the bass as overfishing does, such as lower levels of the Atlantic
menhaden, which plays a big role in the Bass’s diet, low oxygen levels,
climate change, and even poor recreational handling of the fish has caused for
released Bass to die, roughly up to 9% of released bass [5]. According to the
66th Northeast Regional Stock Assessment Workshop Assessment Report by
Northeast Fisheries Science Center, the striped bass stock was overfished and
experienced overfishing in 2017 [24].

Renewable resources and ”common property” fisheries are frequently exploited
at levels that do not account for its effect on the ecosystem and long term
sustainability, which may lead to extinction of some or all overexploited
species. Reducing overexploitation and ensuring fisheries conservation for
future generations is a main challenge for policy-makers. Mathematical models
could help in understanding the complexities involved in this endeavor. Here
we use the Lotka-Volterra predator-prey model to study the population dynamics
of Atlantic menhaden and striped bass and the possibility of balancing
economic gains with ecosystem conservation. We use real time series population
data of the two fish to fit a mathematical model which is then used to study
the effects of harvesting on the economic rent function and long time
sustainability of the fish population. We find the bioeconomic equilibria of
the model as well as the stability of the equilibria. Along with this, we
analyze the effects of applying the maximum sustainable yield (MSY) policy and
the maximum economic yield (MEY) policy on our system and see how these
strategies might be useful to both maximize profit from harvesting and keep a
sustainable fish population. We also study the effects of resilience
maximizing yield (RMY) policy when applied to the Atlantic menhaden-striped
bass system. Such policy is a management approach that prioritizes resilience,
a measure of how fast the system stabilizes back to biological equilibrium
after disruption, and allows fisheries to preserve the longevity and stability
of the fish population over long time. Finally, we perform sensitivity
analysis to measure how small perturbations of our parameter values affect the
population dynamics of the two species in time. This analysis allows to
identify the most sensitive parameters of the system which could be then used
for management and conservation decisions by policy-makers.

##  2 Mathematical Model

###  2.1 Two-species prey-predator model

Let us consider the generalized Lotka-Volterra model of the population
dynamics of a prey $x=x(t)$ and a predator $y=y(t)$ given by the following
system of differential equations

$\left\\{\begin{array}[]{lll}\dfrac{dx}{dt}&=&x(r-ax-by)\\\
\dfrac{dy}{dt}&=&y(-e-cy+dx).\end{array}\right.$ (1)

Here $r$ is the intrinsic growth rate of the prey, $e$ is the death rate of
the predator in absence of prey, $a$ and $c$ are the self-limitation
parameters for the prey and the predator respectively, $b$ is the predator
effect of $y$ on $x$, and $d$ is the preying effect of $x$ on $y$.

###  2.2 Data and Methodology

We used the Northeast Area Monitoring and Assessment Program (NEAMAP) survey
data of the Atlantic menhaden and striped bass conducted by the Multi-species
Research Group at the Virginia Institute of Marine Science [40, 41]. The data
used was the numerical index, a metric representative of population, collected
every Fall from 2008 to 2020. We used the gauseR package [21] to fit the model
(1) to the data. The gauseR package is a collection of various functions that
are made specifically for modeling and fitting curves to predator-prey
dynamics by creating time-lagged abundances, estimating per-capita growth
rates, and parameter tuning via nonlinear least squares method.

###  2.3 Results and Prediction analysis

After fitting the Lotka-Volterra model to the numerical index data for
menhaden and bass, the fitted model is graphed in Figure 1 by the red and blue
curves correspondingly. We can see the phase-shifting behaviour of the
menhaden and bass which portray the natural fluctuations commonly seen within
predator-prey dynamics. The system’s parameters obtained after the fitting are
given in Table 1, with $r^{2}$ of 0.72.

Parameter | Description | Value
---|---|---
r | Intrinsic growth rate of Atlantic menhaden | 0.416
a | Self-limitation of Atlantic menhaden | 0.534
b | Effect of Striped bass on Atlantic menhaden | 0.515
e | Death rate of striped bass | 0.397
c | Self-limitation of striped bass | 0.262
d | Effect of Atlantic menhaden on striped bass | 1.042
Table 1: Table of the model’s parameter values obtained from the fitting.
Figure 1: Fitted curves to the given data for the population of menhaden and
bass from 2008 to 2020.

We also used the model and the obtained parameters to predict the future
population dynamics of menhaden and bass up to year 2050, which is shown in
Figure 2. The system converges to equilibrium past 2040, and using the
obtained parameters from Table 1, we calculate the numerical index equilibrium
$E=(x^{*},y^{*})=(0.4187,0.2958)$ for menhaden and bass respectively.

Figure 2: Prediction of menhaden-bass population dynamics up to 2050.

###  2.4 Resilience of the menhaden-bass system

Resilience ($R$) measures the stability of a system and is defined as the
inverse of the time required for the system to recover from a perturbation and
come back to its original steady state. According to Pimm and Lawton [25], the
time required for the system to return to equilibrum is given by

$\tau=\displaystyle\frac{-1}{Re(\lambda_{m})}\quad\mbox{and}\quad
R=\dfrac{1}{\tau},$

where $Re(\lambda_{m})$ is the real part of the leading eigenvalue of the
system, and $\tau$ is time. Note that we can talk about resilience only if the
system is stable, as real part of the leading eigenvalue will be negative,
hence the resilience $\tau$ will be a positive value. As higher the resilience
of the system is, as smaller is the time for its return to equilibrium if the
system is perturbed from its steady state. In our case study of menhaden and
striped bass as a prey and predator, we found the eigenvalues of the fitted
model to the given observational data corresponding to the co-existence
equilibrium to be

$\lambda_{1}\approx-0.1506+0.2473i,\qquad\lambda_{2}\approx-0.1506-0.2473i$
(2)

As the real part has a negative sign, this means that the co-existence
equilibrium is a stable one and hence we can calculate the resilience of this
system and the return time to equilibrium as

$R\approx 0.1506\quad\mbox{and}\quad\tau\approx\frac{1}{0.1506}\approx 6.64.$

This means that if the system is perturbed from its steady state it will
return to its equilibrium after roughly 6.64 years.

##  3 Ecological and Economic Considerations

Understanding the interactions between species in the Chesapeake Bay is
important, especially when it comes to predator-prey dynamics. However, even
more important is to understand these dynamics when we introduce a human
factor in the form of harvesting into the system. Not only are the Atlantic
menhaden and the striped bass two very important species in the Chesapeake Bay
in terms of their roles in the ecological food chain, but they have proven to
be key economic factors as well. Harvesting these species for economic gain
makes sense, but when does harvesting becomes excessive and hence detrimental
to these species? How does such harvesting affect the predator-prey dynamics?
How can we balance harvesting for economic gains and long-term sustainability
of these species? What are some policies that could be implemented to achieve
such a balance?

In this section we employ mathematical modeling and techniques to study the
complexity of the problem and shed light on the possibility to balance
economically effective policies with the fisheries conservation for future
generations.

###  3.1 Bionomic Model

Let us consider the system where both the prey and the predator are harvested
based on the catch per unit effort hypothesis. Incorporating this into model
(1) yields to a system of two ordinary differential equations. Adding an
economic rent function, representing the profit made from harvesting and
selling, we get the following bioeconomic model, or simply bionomic model:

$\left\\{\begin{array}[]{lll}\displaystyle\frac{dx}{dt}&=&x(r-ax-
by)-q_{1}E_{1}x\\\ \displaystyle\frac{dy}{dt}&=&y(-e-cy+dx)-q_{2}E_{2}y\\\
\pi(E_{1},E_{2})&=&(p_{1}q_{1}x-c_{1})E_{1}+(p_{2}q_{2}y-c_{2})E_{2}\end{array}\right.$
(3)

with initial data values $x(0)\geq 0,y(0)\geq 0$. Here, $E_{1}$ and $E_{2}$
are the harvesting efforts of the menhaden and bass respectively, $q_{1}$ and
$q_{2}$ represent the catchability coefficients of each species. In the right
hand side of $\pi$ are the fishing cost per unit effort for both species,
$c_{1}$ and $c_{2}$, and the price per unit biomass of the species, $p_{1}$
and $p_{2}$.

###  3.2 Positivity and Boundedness of the Solutions

In this section, we determine a region $\Omega$ where the state variables are
positive and bounded. We show that the set defined as follows

$\Omega=\left\\{(x,y)\|\quad 0\leq x\leq\frac{r}{a},0\leq
x+\frac{y}{\gamma}\leq L\right\\},$ (4)

for some $\gamma$ and $L$, which will be defined later, is an invariant set,
meaning that every solution with initial conditions in this set, remains in it
for all $t\geq 0$. Figure 3 illustrates the region, using the parameter values
from Table 1.

Figure 3: The region $\Omega$ obtained using the parameter values from Table
1. Only the first quadrant is biologically relevant.

####  3.2.1 Positivity of the solutions

Rewriting equation one of model (3) in the form

$\frac{dx}{x}=\phi_{1}dt,$

where $\phi_{1}(x,y)=r-ax-by$, and integrating over $[0,t]$ we obtain

$x(t)=x(0)\exp{\left[\int_{0}^{y}\phi_{1}(x,y)ds\right]}.$

Since $x(0)\geq 0$, we can conclude that $x(t)\geq 0$ for all $t\geq 0$.
Similarly, from the second equation of (3), we have

$\frac{dy}{y}=\phi_{2}dt,$

where $\phi_{2}=dx-e-cy$. Integrating both sides over $[0,t]$, we obtain

$y(t)=y(0)\exp{\left[\int_{0}^{y}\phi_{2}(x,y)ds\right]},$

and as the initial condition is positive, $y(0)\geq 0$, then solution
$y(t)\geq 0$ for all $t\geq 0$ which proves the positivity of the solutions in
the defined region $\Omega$.

####  3.2.2 Boundedness of the solutions

Next, we prove that all solutions remain bounded in $\Omega$, which means that
as long as the initial conditions are in $\Omega$, the solutions to (3) remain
in $\Omega$ for all $t$.

As all coefficients are non-negative, we have $\dfrac{dx}{dt}=x(r-ax-
by)-q_{1}E_{1}x\leq x(r-ax)$. Now solving $\dfrac{dx}{dt}=x(r-ax)$ by
separation of variables and using the partial fractions decomposition method,
we get:

$\left(\frac{1}{rx}+\frac{a}{r}\frac{1}{r-ax}\right){dx}={dt}.$

After integrating both sides and using some manipulations, we can express $x$:

$x=\frac{Cre^{rt}}{1+aCe^{rt}}=\frac{Cr}{aC+e^{-rt}}.$

Now letting $t\rightarrow\infty,\text{implies that }x\rightarrow\frac{r}{a}.$
Therefore $x(t)\leq\frac{r}{a}$ is bounded in $\Omega$.

Now we want to show that $y(t)$ is also bounded as $t\rightarrow\infty$. To
prove that, let’s introduce a new variable, $w(t)$, defined as

$w(t)=x(t)+\frac{y(t)}{\gamma},$ (5)

where $\gamma$ is given by $\gamma=\frac{d}{b}$ or equivalently $d=\gamma b$.
Using the first two equations in (3) we obtain the following differential
equation for the new variable $w$

$\frac{dw}{dt}=\frac{dx}{dt}+\frac{1}{\gamma}\frac{dy}{dt}=x(r-ax)-q_{1}E_{1}x-\frac{1}{\gamma}(ey+cy^{2}+q_{2}E_{2}y).$

Using the fact that $-cy^{2}<0$, we find the differential inequality:

$\frac{dw}{dt}\leq x(r-ax)-q_{1}E_{1}x-\frac{y}{\gamma}(e+q_{2}E_{2}).$

Observe that if we let $f(x)=x(r-ax)$, then $f(x)$ obtains its maximum at
$x=\dfrac{r}{2a}$ and that maximum is $f(\dfrac{r}{2a})=\dfrac{r^{2}}{4a}$.
Substituting this back into the differential inequality, along with letting
$\delta=\min{(q_{1}E_{1},q_{2}E_{2},e)}$, we get:

$\displaystyle\frac{dw}{dt}$ $\displaystyle\leq
x(r-ax)-q_{1}E_{1}x-\frac{y}{\gamma}(e+q_{2}E_{2})$
$\displaystyle\leq\frac{r^{2}}{4a}-\delta(x+\frac{y}{\gamma})=\frac{r^{2}}{4a}-\delta
w.$

Thus, the new differential inequality for $w$ is as follows

$\frac{dw}{dt}+\delta w\leq\frac{r^{2}}{4a}.$

Using the integrating factor, multiplying both sides by $e^{\delta t}$ and
integrating, after some manipulations we get

$w(t)\leq\frac{w(0)}{e^{\delta t}}-\frac{r^{2}}{4a\delta e^{\delta
t}}+\frac{r^{2}}{4a\delta}=\frac{1}{e^{\delta
t}}(w(0)-\frac{r^{2}}{4a\delta})+\frac{r^{2}}{4a\delta}.$

From here, as $t\rightarrow\infty$, $e^{\delta t}$ $\rightarrow\infty$, which
in turn means that the term $\frac{1}{e^{\delta
t}}(w(0)-\frac{r^{2}}{4a\delta})\rightarrow{0}$ as $t\rightarrow\infty$,
reducing the above inequality to

$\lim_{t\to\infty}w(t)\leq\frac{r^{2}}{4a\delta}.$

From here $w=x+\frac{y}{\gamma}$ is bounded and hence are solutions in the
region $\Omega$, where $L=\frac{r^{2}}{4a\delta}$. Therefore, the region
$\Omega$ is an invariant set, and any solution that has initial conditions in
$\Omega$ remains in $\Omega$. With this, all solutions of the model (3) are
nonnegative and bounded in $\Omega$.

###  3.3 Bionomic Equilibrium Existence

The bionomic equilibrium is obtained when $dx/dt=dy/dt=0$ and when the total
revenue earned by selling the harvested biomass is equal to the total cost for
the effort in harvesting the biomass. Here, we want to find conditions under
which bionomic equilibria exist. Suppose the bionomic equilibria is
$P(x_{B},y_{B},E_{1B},E_{2B})$. Then $P$ will satisfy the system

$\left\\{\begin{array}[]{rrr}r-ax-by-q_{1}E_{1}&=&0\\\ -e-cy+dx-
q_{2}E_{2}&=&0\\\
(p_{1}q_{1}x-c_{1})E_{1}+(p_{2}q_{2}y-c_{2})E_{2}&=&0.\end{array}\right.$ (6)

Note that we did not consider the biological equilibrium $E_{0}(x,y)=(0,0)$
because in this case, the overall net revenue $\pi$ becomes negative. To find
a solution to (6) we will consider four different cases based on the
harvesting efforts.

#### Case 1. Harvesting the predator only ($E_{1}=0$)

If the fishing cost per unit effort for the prey (menhaden) is greater than
the revenue obtained from this fishing effort, $c_{1}>p_{1}q_{1}x$, the
fishermen will naturally withdraw from this fishery and as a result the
fishing effort for menhaden will be zero $E_{1}=0$, but the effort for striped
bass will be not zero $E_{2}\neq 0$. Then for the third equation in (6) to be
true, we can take $y_{B}=c_{2}/p_{2}q_{2}$. From the first equation of (6) we
can express $x_{B}$:

$x_{B}=\displaystyle\frac{r}{a}-\frac{bc_{2}}{ap_{2}q_{2}}.$

While from the second equation we can express $E_{2B}$ such as

$E_{2B}=-\displaystyle\frac{e}{q_{2}}-\frac{cc_{2}}{p_{2}q_{2}^{2}}+\frac{d}{q_{2}}\left(\frac{r}{a}-\frac{bc_{2}}{ap_{2}q_{2}}\right).$
(7)

######  Theorem 3.1.

If there are no harvesting efforts on the prey, the bionomic equilibrium

$P_{1}\left(x_{B},y_{B},E_{1B},E_{2B}\right)=P_{1}\left(\displaystyle\frac{r}{a}-\frac{bc_{2}}{ap_{2}q_{2}},\frac{c_{2}}{p_{2}q_{2}},0,E_{2B}\right),$

where $E_{2B}$ is given by (7), exists provided that the following two
conditions hold:

$\displaystyle\frac{r}{a}>\frac{bc_{2}}{ap_{2}q_{2}}$
$\frac{d}{q_{2}}\left(\frac{r}{a}-\frac{bc_{2}}{ap_{2}q_{2}}\right)>\displaystyle\frac{e}{q_{2}}+\frac{cc_{2}}{p_{2}q_{2}^{2}}.$

#### Case 2. Harvesting the prey only ($E_{2}=0$)

Likewise, if the fishing cost per unit effort for the predator (striped bass)
is greater than the revenue obtained from this fishing effort,
$c_{2}>p_{2}q_{2}x$, the fishermen will naturally withdraw from this fishery
and as a result the fishing effort for the striped bass will be zero
$E_{2}=0$, but the effort for menhaden will be not zero $E_{1}\neq 0$. Then
for the third equation in (6) to be true, we can take
$x_{B}=c_{1}/p_{1}q_{1}$. From the second equation of (6) we can express
$y_{B}$:

$y_{B}=\displaystyle\frac{dc_{1}}{cp_{1}q_{1}}-\frac{e}{c}.$

While from the first equation we can express $E_{1B}$ such as

$E_{1B}=\displaystyle\frac{r}{q_{1}}-\frac{ac_{1}}{p_{1}q_{1}^{2}}-\frac{b}{q_{1}}\left(\frac{dc_{1}}{cp_{1}q_{1}}-\frac{e}{c}\right).$
(8)

######  Theorem 3.2.

If there are no harvesting efforts on the predator, the bionomic equilibrium

$P_{2}\left(x_{B},y_{B},E_{1B},E_{2B}\right)=P_{2}\left(\frac{c_{1}}{p_{1}q_{1}},\displaystyle\frac{dc_{1}}{cp_{1}q_{1}}-\frac{e}{c},E_{1B},0\right),$

where $E_{1B}$ is given by (8), exists provided that the following two
conditions hold:

$\displaystyle\frac{dc_{1}}{cp_{1}q_{1}}>\frac{e}{c}$
$\displaystyle\frac{r}{q_{1}}-\frac{ac_{1}}{p_{1}q_{1}^{2}}>\frac{b}{q_{1}}\left(\frac{dc_{1}}{cp_{1}q_{1}}-\frac{e}{c}\right).$

#### Case 3. No harvesting of both the prey and the predator ($E_{1}=E_{2}=0$)

If $c_{1}>p_{1}q_{1}x$ and $c_{2}>p_{2}q_{2}y$, then the fishing cost is
greater than the revenue for menhaden and striped bass and the fisheries will
be closed. We can not have a bionomic equilibrium in this case, however the
biological equilibrium $(x^{*},y^{*})$ will still exist if $rd>ae$:

$x^{*}=\frac{rc+be}{ac+bd},\qquad y^{*}=\frac{rd-ae}{ac+bd}.$

#### Case 4. Simultaneous harvesting of both the prey and the predator
($E_{1}\neq 0,E_{2}\neq 0$)

If $c_{1}\leq p_{1}q_{1}x$ and $c_{2}\leq p_{2}q_{2}y$, then the cost of
fishing is less or equal to the revenue for both the menhaden and the striped
bass, and the fisheries have an incentive to harvest both species, and thus
$E_{1}\neq 0$ and $E_{2}\neq 0$. Therefore, to find the bionomic equilibrium
of (6) and in order for the third equation to be true, we can set $x_{B}$ and
$y_{B}$ as

$x_{B}=\frac{c_{1}}{p_{1}q_{1}},\quad y_{B}=\frac{c_{2}}{p_{2}q_{2}}$

respectively. From these, we can express $E_{1}$ and $E_{2}$ from the first
and second equations in (6) respectively:

$\begin{array}[]{ccc}E_{1B}&=&\displaystyle\frac{r}{q_{1}}-\frac{ac_{1}}{p_{1}q_{1}^{2}}-\frac{bc_{2}}{p_{2}q_{1}q_{2}}\\\
E_{2B}&=&\displaystyle\frac{-e}{q_{2}}-\frac{cc_{2}}{p_{2}q_{2}^{2}}+\frac{dc_{1}}{p_{1}q_{1}q_{2}}.\end{array}$
(9)

######  Theorem 3.3.

If both the prey and the predator are harvested, the bionomic equilibrium

$P_{3}\left(x_{B},y_{B},E_{3B},E_{3B}\right)=P_{1}\left(\displaystyle\frac{c_{1}}{p_{1}q_{1}},\frac{c_{2}}{p_{2}q_{2}},E_{1B},E_{2B}\right),$

where $E_{1B}$ and $E_{2B}$ are given by (9), exists provided that the
following two conditions hold:

$\displaystyle\frac{r}{q_{1}}>\frac{ac_{1}}{p_{1}q_{1}^{2}}+\frac{bc_{2}}{p_{2}q_{1}q_{2}}$
$\displaystyle\frac{dc_{1}}{p_{1}q_{1}q_{2}}>\frac{e}{q_{2}}+\frac{cc_{2}}{p_{2}q_{2}^{2}}.$

###  3.4 Effect of harvesting on stability of the prey-predator system

To study the stability of the biological system when either the prey, or the
predator, or both prey and predator are harvested we will use the Jacobian of
the system (6):

$J(x,y)=\left(\begin{array}[]{ccc}r-2ax-by-q_{1}E_{1}&-bx\\\ dy&-e-2cy+dx-
q_{2}E_{2}\end{array}\right).$ (10)

#### Case 1. Harvesting the predator only

If only the predator is harvested, the co-existence equilibrium
$(\hat{x},\hat{y})$ is the solution to the system:

$\left\\{\begin{array}[]{ccc}ax+by&=&r\\\ dx-cy&=&e+q_{2}E_{2}.\\\
\end{array}\right.$ (11)

Using Cramer’s rule, the solutions are

$\hat{x}=\displaystyle\frac{b(e+q_{2}E_{2})+rc}{ac+bd},$ (12)

$\hat{y}=\displaystyle\frac{rd-a(e+q_{2}E_{2})}{ac+bd}.$ (13)

Note that $\hat{x}>0$ is always positive, hence for the co-existence
equilibrium to exist we need to require $\hat{y}>0$, which provides a
constraint on the harvesting effort for the predator:

$E_{2}<\frac{rd-ae}{aq_{2}}.$ (14)

Substituting $(\hat{x},\hat{y})$ into the Jacobian matrix (10) and using the
properties of the equilibrium solution, we get the simplified Jacobian matrix:

$J(\hat{x},\hat{y})=\left(\begin{array}[]{ccc}-a\hat{x}&-b\hat{x}\\\
d\hat{y}&-c\hat{y}\end{array}\right)=A.$ (15)

Its characteristic equation is the quadratic equation

$\lambda^{2}-trace(A)\lambda+det(A)=0,$

where $trace(A)=-(a\hat{x}+c\hat{y})$ and $det(A)=(ac+bd)\hat{x}\hat{y}$. Note
that the trace is always negative, while the determinant is always positive
which implies that the eigenvalues will have negative real parts. Hence, the
system is asymptotically stable for any harvesting effort $E_{2}$ that
satisfies condition (14).

Using the obtained model parameters from Table 1 and assuming $q_{2}=1$, this
means that the co-existence equilibrium for menhaden and striped bass is
asymptotically stable as long as the effort is $E_{2}\lessapprox 0.4147$. Note
that it is not required for the effort to be positive, negative values of
$E_{2}$ mean that harvesting effort is decreased from the current bass
harvesting. To analyse the resilience of the menhaden-bass system we calculate
the return time to equilibrium shown in Figure 4. One interesting findings is
the fact that the return time of the menhaden-bass system to equilibrium
decreases monotonically as a function of $E_{2}$ linear-like up until
$E_{2}\approx 0.3102$. Hence, even with additional effort, when harvesting the
predator only, the menhaden-bass system is more resilient to small
perturbations in their populations. For $E_{2}>0.3102$, the return time to
equilibrium increases significantly and eventually grows asymptotically
towards infinity as $E_{2}$ approaches the upper bound of $\sim 0.4147$.

Figure 4: Return time to equilibrium $\tau$ versus the effort $E_{2}$ when
harvesting the predator only.

#### Case 2. Prey harvesting only

If only the predator is harvested, the co-existence equilibrium
$(\bar{x},\bar{y})$ is the solution to the system:

$\left\\{\begin{array}[]{ccc}r-ax-by&=&q_{1}E_{1}\\\ -e-cy+dx&=&0\\\
\end{array}\right.$ (16)

Using Cramer’s rule, the solutions are

$\bar{x}=\displaystyle\frac{c(r-q_{1}E_{1})+be}{ac+bd},$ (17)
$\bar{y}=\displaystyle\frac{d(r-q_{1}E_{1})-ae}{ac+bd}.$ (18)

Note that $\bar{x}>0$ requires that $r>q_{1}E_{1}$, or equivalently,
$E_{1}<\displaystyle\frac{r}{q_{1}}$ which is a natural condition for
existence of the menhaden, as otherwise, if $r<q_{1}E_{1}$, this means that
the menhaden are being farmed to extinction without predation from the bass.
Now, for the co-existence equilibrium to be biologically relevant, we need to
require $\bar{y}>0$, which provides a constraint on the prey’s harvesting
effort:

$E_{1}<\frac{r}{q_{1}}-\frac{ae}{dq_{1}}.$ (19)

Observe also that from the condition (19) for the effort $E_{1}$ follows that
$E_{1}<\displaystyle\frac{r}{q_{1}}$ and hence $\bar{x}>0$. Substituting
$(\bar{x},\bar{y})$ into the Jacobian matrix (10) and using the properties of
the equilibrium solution, we get the simplified Jacobian matrix:

$J(\bar{x},\bar{y})=\left(\begin{array}[]{ccc}-a\bar{x}&-b\bar{x}\\\
d\bar{y}&-c\bar{y}\end{array}\right)=\bar{A}.$ (20)

Similar considerations about the trace and determinant of $\bar{A}$ as before
lead to the conclusion that the coexistence equilibrium $(\bar{x},\bar{y})$ is
asymptotically stable provided condition (19) for the prey effort $E_{1}$ is
satisfied. Using the fitted model’s parameters from Table 1 and assuming
$q_{1}=1$, this gives that the equilibrium is asymptotically stable as long as
the effort on menhaden is bounded by $E_{1}\lessapprox 0.2125$. Note that
negative effort means that the menhaden harvesting is decreased from the
current level. To analyse the resilience of the menhaden-bass system we
calculate the return time to equilibrium shown in Figure 5. The return time to
equilibrium increases monotonically as a function of $E_{1}$ linear-like up
until $E_{1}\approx 0.1824$. After that it increases significantly and
eventually grows asymptotically towards infinity as $E_{1}$ approaches the
upper bound of $\sim 0.2125$.

Figure 5: Return time to equilibrium $\tau$ versus the effort $E_{1}$ when
harvesting the prey only.

#### Case 3. Combined harvesting of both the prey and the predator

Let us now apply a combined harvesting effort $E_{1}=E_{2}=E$ of both the prey
and predator. The coexistence equilibrium $(x^{*},y^{*})$ is a solution to the
system:

$\left\\{\begin{array}[]{ccc}ax+by&=&r-q_{1}E\\\ dx-cy&=&e+q_{2}E.\\\
\end{array}\right.$ (21)

Here, the solutions are

$x^{*}=\displaystyle\frac{(cr+be)-(cq_{1}-bq_{2})E}{ac+bd},$ (22)
$y^{*}=\displaystyle\frac{(rd-ae)-(dq_{1}+aq_{2})E}{ac+bd}.$ (23)

Therefore, in order for the coexistence equilibrium to exist we need to
require $x^{*}>0$ and $y^{*}>0$. If $cq_{1}<bq_{2}$, then the $x^{*}$ will be
always positive for a positive effort. If $cq_{1}>bq_{2}$, we get the
following constraint on the harvesting effort $E$:

$0<E<\frac{cr+be}{cq_{1}-bq_{2}}.$

Requiring $y^{*}>0$ provides another constraint on the effort:

$0<E<\displaystyle\frac{rd-ae}{dq_{1}+aq_{2}}.$

Hence, for the coexistence equilibrium to exist, the effort must satisfy the
following condition:

$E<\min\left\\{\frac{cr+be}{cq_{1}-bq_{2}},\frac{rd-
ae}{dq_{1}+aq_{2}}\right\\}$ (24)

Substituting $(\hat{x},\hat{y})$ into the Jacobian matrix, we get the
following:

$J(x^{*},y^{*})=\left(\begin{array}[]{ccc}r-2ax^{*}-by^{*}-q_{1}E&-bx^{*}\\\
dy^{*}&-e-2cy^{*}+dx^{*}-q_{2}E\end{array}\right)=\left(\begin{array}[]{ccc}-ax^{*}&-bx^{*}\\\
dy^{*}&-cy^{*}\end{array}\right)=A^{*}.$

Therefore, the coexistence equilibrium $(x^{*},y^{*})$ is asymptotically
stable as long as condition (24) for the combined harvesting effort is
satisfied. As $q_{2}\geq\frac{c}{b}\approx 0.5087$, $q_{2}$ must be less than
$\sim 0.5087$ in order for $E\geq 0$. Using our obtained parameters and
assuming $q_{1}=1,q_{2}=0.5$, we get that $\dfrac{cr+be}{cq_{1}-bq_{2}}\approx
69.6549$ and $\dfrac{rd-ae}{dq_{1}+aq_{2}}\approx 0.1692$. Thus, the
equilibrium is asymptotically stable as long as $E\lessapprox 0.1692$. Note
that for stability it is not required for the effort to be positive; a
negative effort means that the harvesting is decreased from the current level.
To analyse the resilience of the menhaden-bass system we calculate the return
time to equilibrium shown in Figure 6. The return time to equilibrium
increases monotonically as a function of $E$ linear-like up until $E\approx
0.1692$. After that it increases exponentially and eventually grows
asymptotically towards infinity as $E$ approaches the upper bound of $\sim
0.1692$.

Figure 6: Return time to equilibrium $\tau$ versus the effort $E$ when
harvesting both the prey and the predator with a same effort $E$.

##  4 Analysis of Harvesting Policies

###  4.1 Maximum Sustainable Yield

We now want to analyze how existing fisheries management policies might affect
the biological sustainability and economic gain within our bionomic system
(3). We will first consider the implementation of an maximum sustainable yield
(MSY) policy on our system. MSY is a theoretical management approach, as it is
the supposed highest catch of a species that maximizes the economic rent while
simultaneously ensuring that the species population is still at a sustainable
level [13]. The concept of MSY originated in Schaefer’s 1954 paper [28], and
has since been widely considered and researched [3, 17, 19, 30, 31, 32, 36,
38]. Applying MSY to our bionomic system allows for a conceptually easy
comparison to other policies we consider.

As we wish to harvest the fisheries such that we can get maximum yield
sustainably, our goal is to maximize the yield function
$Y(E)=q_{1}Ex^{*}+q_{2}Ey^{*}$ with respect to $E$ when both species are
harvested with the same effort $E$. Here $(x^{*},y^{*})$ is the equilibrium
point of the prey-predator system.

To illustrate how MSY policy could be applied let us use Figure 7. The yield
function is in green and the prey and predator equilibria curves, which are
linearly dependent on effort, are given in red and blue respectively. The
effort level corresponding to MSY is $E_{MSY}$. If the predator equilibrium
curve is similar to the bottom blue line, then the effort level corresponding
to the maximum of the yield curve, denoted as $E_{Y_{m}}$, will result in a
negative numerical index, or population, of the predator. This means that the
predator will go extinct if this effort level is used, and therefore $E_{MSY}$
will be at a point before the x-intercept of the equilibrium curve, resulting
in an MSY value lower than the theoretical maximum of $Y(E)$. If the predator
curve follows the middle blue line, then harvesting at $E_{Y_{m}}$ will result
in a positive, but near 0, equilibrium population of the predator. This means
that while sustainable, the predator population is at high existential risk
due to the practical difficulties of low populations in reality. Lastly, if
the predator’s equilibrium curve follows the top blue line, then the predator
population becomes much less at risk, and MSY becomes the maximum of the yield
function and can be obtained with implementing the effort level of
$E_{Y_{m}}$.

Figure 7: A general case for two species MSY analysis with multiple predator
equilibria lines.

#### Case 1. Harvesting the predator only

Now we will consider each of the three scenarios: when only the predator is
harvested, only the prey is harvested and when both species are harvested. Let
us assume that only the predator is harvested with an effort $E$. As we
already calculated the equilibrium $(\hat{x},\hat{y})$ in (12) and (13), the
yield at the equilibrium as a function of the effort defined as follows

$Y(E)=E\hat{y}(E)=E\left(\displaystyle\frac{rd-
ae}{ac+bd}-\frac{aq}{ac+bd}E\right),$ (25)

where $E_{2}=E$ and $q_{2}=q$. As the yield function is a quadratic function
with respect to $E$, $E_{Y_{m}}=\displaystyle\frac{rd-ae}{2aq}$. Next, we can
calculate both $\hat{x}$ and $\hat{y}$ with respect to $E$ as follows

$\hat{x}(E)=\displaystyle\frac{be+rc}{ac+bd}+\frac{bq}{ac+bd}E,$
$\hat{y}(E)=\displaystyle\frac{rd-ae}{ac+bd}-\frac{aq}{ac+bd}E,$

where both $\hat{x}$ and $\hat{y}$ are linearly dependent on $E$, with
$\hat{x}$ increasing and $\hat{y}$ decreasing as functions of $E$.

Figure 8: Yield function and the species equilibria curves for menhaden and
bass when only the predator is harvested and parameters from Table 1 are used,
with $q=1$.

Note that predator harvesting reduces the population of the predator however
increasing the prey’s population at equilibrium. This is illustrated in Figure
(8) using the values from Table 1. As $\hat{y}(E_{Y_{m}})$ and
$\hat{y}(E_{Y_{m}})$ are both positive, the maximum yield can be attained by
harvesting at $E_{Y_{m}}$, with corresponding equilibrium values of
$\hat{x}(E_{MSY})=\hat{x}(E_{Y_{m}})=\displaystyle\frac{be+rc}{ac+bd}+\frac{bq}{ac+bd}E_{Y_{m}}$
and $\hat{Y}(E_{MSY})=\hat{y}(E_{Y_{m}})=\displaystyle\frac{rd-ae}{2(ac+bd)}$.
The predator’s population in this case is half of its population in the
absence of harvesting. In order for MSY to exist, we need $rd>ae$, therefore
we may state the following theorem about MSY in this case.

######  Theorem 4.1.

In the system (3), if $rd>ae$, MSY will exist and be obtained as a result of
harvesting the predator with effort level $E_{MSY}$.

######  Corollary 4.1.1.

For the menhaden-bass system under consideration, harvesting the menhaden at
the MSY level can occur as Theorem 4.1 is satisfied as, using the parameters
found in Table 1, $rd\approx 0.4335$ and $ae\approx 0.2120$, and thus $rd>ae$.

#### Case 2. Harvesting the prey only

Next, let us consider the case when only the prey is harvested with an effort
$E$. For our equilibrium $(\hat{x},\hat{y})$ from (17) and (18), defining the
yield at equilibrium as a function of effort follows

$Y(E)=E\hat{x}(E)=E\left(\displaystyle\frac{cr+be}{ac+bd}-\frac{cq}{ac+bd}E\right),$
(26)

where $E_{1}=E$ and $q_{1}=q$. As the yield function is a quadratic function
with respect to $E$, $E_{Y_{m}}=\displaystyle\frac{cr+be}{2cq}$. Next, we can
calculate both $\hat{x}$ and $\hat{y}$ with respect to $E$ as follows

$\hat{x}(E)=\displaystyle\frac{cr+be}{ac+bd}-\frac{cq}{ac+bd}E,$
$\hat{y}(E)=\displaystyle\frac{rd-ae}{ac+bd}-\frac{dq}{ac+bd}E,$

where both $\hat{x}$ and $\hat{y}$ are linearly dependent on $E$, with both
$\hat{x}$ and $\hat{y}$ decreasing as functions of $E$.

Figure 9: Yield function and the species equilibria curves for menhaden and
bass when only the prey is harvested and parameters from Table 1 are used,
with $q=1$.

Figure (9) shows the yeld function and the equilibria curves for menhaden and
bass using the parameters from Table 1. It could be seen that prey harvesting
reduces the population of both the prey and predator at equilibrium. As
$\hat{y}(E_{Y_{m}})$ is negative, harvesting at $E_{Y_{m}}$ will result in the
extinction of the predator and thus we have
$E_{MSY}<E_{Y_{m}}<\frac{r}{q_{1}}-\frac{ae}{dq_{1}}$ using condition 19.
Observe that there is a trade off when the effort $E_{MSY}$ increases close to
the effort that provides maximum yield $E_{Y_{m}}$, the larger MSY gets the
smaller predator’s equilibrium population gets and eventually decreases to 0.
The corresponding equilibrium values are
$\hat{x}(E_{MSY})=\displaystyle\frac{cr+be}{ac+bd}-\frac{cq}{ac+bd}E_{MSY}$
and $\hat{y}(E_{MSY})=\displaystyle\frac{rd-
ae}{ac+bd}-\frac{dq}{ac+bd}E_{MSY}$. Hence, in order for MSY to exist, we need
$rd>ae$ and therefore we may state the following theorem about MSY in this
case.

######  Theorem 4.2.

In the system (3), if $rd>ae$, MSY will exist and could be obtained as a
result of harvesting the prey with some effort level
$E_{MSY}<E_{Y_{m}}<\frac{rd-ac}{dq_{1}}$.

######  Corollary 4.2.1.

For the menhaden-bass system under consideration, harvesting the striped bass
at the MSY level can occur as Theorem 4.2 is satisfied as, using the
parameters found in Table 1, $rd\approx 0.4335$ and $ae\approx 0.2120$, and
thus $rd>ae$.

#### Case 3. Combined harvesting of both the prey and the predator

Let us assume combined harvesting of both the prey and the predator with
effort $E$. For our equilibrium $(\hat{x},\hat{y})$ from (22) and (23),
defining the yield at equilibrium as a function of effort follows

$Y(E)=E(q_{1}\hat{x}+q_{2}\hat{y})=E\left(\displaystyle\frac{cr+be+rd-
ae}{ac+bd}-\frac{cq_{1}-bq_{2}+dq_{1}+aq_{2}}{ac+bd}E\right)$ (27)

where $E$ is the effort level of harvesting effort for both prey and predator.
As $Y(E)$ is quadratic, $E_{Y_{m}}=\displaystyle\frac{cr+be+rd-
ae}{2(cq_{1}-bq_{2}+dq_{1}+aq_{2})}$. Next, we can calculate both $\hat{x}$
and $\hat{y}$ with respect to E as follows

$\hat{x}(E)=\displaystyle\frac{cr+be}{ac+bd}-\frac{cq_{1}-bq_{2}}{ac+bd}E,$
$\hat{y}(E)=\displaystyle\frac{rd-ae}{ac+bd}-\frac{dq_{1}+aq_{2}}{ac+bd}E,$

where both $\hat{x}$ and $\hat{y}$ are linearly dependent on $E$, with
$\hat{y}$ is decreasing as the effort increases. For $\hat{x}$, there are two
cases, depending on the values of $cq_{1}$ and $bq_{2}$. In the first case,
when $cq_{1}<bq_{2}$, $\hat{x}(E)$ is increasing. In the second case, where
$cq_{1}>bq_{2}$, $\hat{x}(E)$ is decreasing.

Figure 10: Yield function and the species equilibria curves for menhaden and
bass when both prey and predator are harvested and parameters from Table 1 are
used, with $q_{1}=1$ and $q_{2}=0.5$.

Figure (10) shows the yield curve and the equilibria curves for menhaden and
bass using the parameters from Table 1. Simultaneous prey and predator
harvesting reduces the population of both the prey and predator at
equilibrium. As $\hat{y}(E_{Y_{m}})$ is negative, harvesting at $E_{Y_{m}}$
will result in the extinction of the predator. Thus
$E_{MSY}<E_{Y_{m}}<\min\left\\{\frac{cr+be}{cq_{1}-bq_{2}},\frac{rd-
ae}{dq_{1}+aq_{2}}\right\\}$, using condition 24. Although there is an upper
bound on $E_{MSY}$, there is uncertainty as to what its exact value should be.
Here, the same trade off seen from the harvesting the prey only case occurs,
in that as $E_{MSY}$ increases to $E_{Y_{m}}$, the larger MSY gets, but at the
cost of the predator’s equilibrium population which could get extinct. The
corresponding equilibrium values are
$\hat{x}(E_{MSY})=\displaystyle\frac{cr+be}{ac+bd}-\frac{cq_{1}-bq_{2}}{ac+bd}E_{MSY}$
and $\hat{y}(E_{MSY})=\displaystyle\frac{rd-
ae}{ac+bd}-\frac{dq_{1}+aq_{2}}{ac+bd}E_{MSY}$. Hence, in order for MSY to
exist, we need $rd>ae$ and therefore we may state the following theorem about
MSY in this case.

######  Theorem 4.3.

In the system (3), if $rd>ae$, MSY will exist and be obtained as a result of
harvesting both the prey and predator with some effort level such that
$E_{MSY}<E_{Y_{m}}<\min\left\\{\frac{cr+be}{cq_{1}-bq_{2}},\frac{rd-
ae}{dq_{1}+aq_{2}}\right\\}$.

######  Corollary 4.3.1.

For the menhaden-bass system under consideration, harvesting both the menhaden
and striped bass at the MSY level can occur as Theorem 4.3 is satisfied as,
using the parameters found in Table 1, $rd\approx 0.4335$ and $ae\approx
0.2120$, and thus $rd>ae$.

Note that harvesting both menhaden and bass with the same effort at MSY level
affects the predator population quite severely as they are at higher risk of
extinction while the prey is practically not affected.

###  4.2 Maximum Economic Yield

Although the implementation of an MSY policy may theoretically be useful in
sustaining fish populations while maximizing economic rent, it’s not without
its downsides. A potential issue with harvesting at MSY is that due to
inaccurate data and models, harvesting at the theoretical effort corresponding
to MSY may result in over harvesting - leading to at risk or extinct species.
This is where a different policy, maximum economic yield (MEY), can help
minimize risk in this regard. Most define MEY as the largest difference
between total revenue earned from fishing and total costs of fishing. In other
words, MEY is the largest net profit, or rent, that is earned from fishing. In
order to understand why MEY may be the better option over MSY, it’s important
to analyze how they theoretically compare to one another.

#### Case 1. Harvesting the predator only

MEY is defined as the largest difference between the total revenue and total
costs. Here, we need three equations:

$R(E)=Y(E)\hat{p}=E(\frac{rd-ae}{ac+bd}-\frac{aq}{ac+bd}E)\hat{p},$ (28)

$C(E)=\hat{c}E,$ (29)

$Rent(E)=R(E)-C(E),$ (30)

where $R(E)$ is our revenue function, composed of the yield equation for
harvesting the predator only defined in (25) times $\hat{p}$, which is the
price per unit of fish, $C(E)$ is our cost function, representing the total
costs of fishing, defined by the cost per unit of effort times harvesting
effort, and $Rent(E)$ is the difference between total revenue and total costs,
which is also the objective in MEY policy. Observe that all are functions of
the harvesting effort $E$. To find the largest difference between revenue and
costs, we need to find the critical points. If we plug in the equations (28)
and (29) into our rent function and simplify, we get

$Rent(E)=\frac{\hat{p}(rd-
ae)}{ac+bd}E-\frac{\hat{p}(aq)}{ac+bd}E^{2}-\hat{c}E,$

and finding the critical points by solving for $E$ at

$\frac{d}{dE}(Rent(E))=\frac{\hat{p}(rd-
ae)}{ac+bd}-\frac{2\hat{p}aq}{ac+bd}E-\hat{c}=0,$

we find that the harvesting effort at MEY is

$E_{MEY}=\frac{\hat{p}(rd-ae)-\hat{c}(ac+bd)}{2\hat{p}aq}.$ (31)

In order to compare both $E_{MSY}$ and $E_{MEY}$ for the harvesting of the
predator only, we need for them to have a common denominator. This can be done
by modifying $E_{MSY}$ to include price:

$E_{MSY}\frac{\hat{p}}{\hat{p}}=\frac{rd-
ae}{2aq}\frac{\hat{p}}{\hat{p}}=\frac{\hat{p}(rd-ae)}{2\hat{p}aq}.$ (32)

Now comparing both $E_{MSY}$ and $E_{MEY}$, we now see that the only
difference between them is the subtracted term $\hat{c}(ac+bd)$ within
$E_{MEY}$. Because $\hat{c}(ac+bd)>0$, this shows that $E_{MEY}<E_{MSY}$,
meaning that, when harvesting the predator only, it takes less effort to
achieve maximum economic yield than it does to achieve maximum sustainable
yield. Additionally, it follows that, although MEY is less than MSY, MEY is
still sustainable.

#### Case 2. Harvesting the prey only

In the case of harvesting the prey only, our rent function becomes

$Rent(E)=\frac{(cr+be)\hat{p}}{ac+bd}E-\frac{\hat{p}cq}{ac+bd}E^{2}-\hat{c}E.$

This time our yield function comes from (26). Solving for the critical points
by solving for E at

$\frac{d}{dE}(Rent(E))=\frac{(cr+be)\hat{p}}{ac+bd}-\frac{2\hat{p}cq}{ac+bd}E-\hat{c}=0$

gives us the MEY harvesting effort

$E_{MEY}=\frac{(cr+be)\hat{p}-\hat{c}(ac+bd)}{2\hat{p}cq}.$ (33)

Similarly to the case of harvesting the predator only, we want to compare
$E_{MSY}$ and $E_{MEY}$. By modifying $E_{MSY}$:

$E_{MSY}\frac{\hat{p}}{\hat{p}}=\frac{cr+be}{2cq}\frac{\hat{p}}{\hat{p}}=\frac{\hat{p}(cr+be)}{2\hat{p}cq},$
(34)

we again see that the only difference between $E_{MSY}$ and $E_{MEY}$ is the
subtracted term $\hat{c}(ac+bd)$. From this, we can see that $E_{MEY}<E_{MSY}$
for harvesting the prey only.

#### Case 3. Combined harvesting of both the prey and predator

In the case of combined harvesting of both prey and predator, our rent
function becomes

$Rent(E)=\frac{\hat{p}(cr+be+rd-
ae)}{ac+bd}E-\frac{\hat{p}(cq_{1}-bq_{2}+dq_{1}+aq_{2})}{ac+bd}E^{2}-\hat{c}E,$

with the yield function coming from (27). Solving for critical points by
solving for E at

$\frac{d}{dE}(Rent(E))=\frac{\hat{p}(cr+be+rd-
ae)}{ac+bd}-\frac{2\hat{p}(cq_{1}-bq_{2}+dq_{1}+aq_{2})}{ac+bd}E-\hat{c}=0$

gives us the MEY harvesting effort

$E_{MEY}=\frac{\hat{p}(cr+be+rd-
ae)-\hat{c}(ac+bd)}{2\hat{p}(cq_{1}-bq_{2}+dq_{1}+aq_{2})}.$ (35)

After the modification of $E_{MSY}$

$E_{MSY}\frac{\hat{p}}{\hat{p}}=\frac{cr+be+rd-
ae}{2(cq_{1}-bq_{2}+dq_{1}+aq_{2})}\frac{\hat{p}}{\hat{p}}=\frac{\hat{p}(cr+be+rd-
ae)}{2\hat{p}(cq_{1}-bq_{2}+dq_{1}+aq_{2})},$ (36)

we see that, when comparing $E_{MSY}$ and $E_{MEY}$, that $E_{MEY}<E_{MSY}$
for the combined harvesting of both prey and predator. Hence, even though
sustainable, MSY policy does not result in a higher net yield compared to MEY
policy and consequently, is not economically justified as it does not increase
profit while reducing the species populations.

###  4.3 Resilience Maximizing Yield

Resilience maximizing yield (RMY) is a management approach that prioritizes
resilience, a measure of how fast the system stabilizes back to biological
equilibrium after disruption. Doing so allows fisheries to preserve the
longevity and stability of the fish populations under consideration. Here, we
continue to analyze additional effort onto the current implementation in
practice. The plots feature varying levels of the relevant $q_{i}$, the
catchability coefficient representative of how “hard” the fish is to catch, to
illustrate that the general behavior of the resilience curves will remain the
same. However, it also shows that varying the $q_{i}$ can have significant
economical implications when considering RMY. Vertical dashed lines indicate
when the resilience of the system is equal to 0, meaning that the system can
no longer recover as one of the species has gone extinct. Stars indicate the
maximum of the resilience curves at effort levels $E_{RMY}$, with their
counterparts plotted on the corresponding yield curves, with coordinates
($E_{RMY}$, $Y(E_{RMY})$). Note that the resilience curves in all cases have
initial value of approximately the original resilience of the system without
harvesting from 2.

#### Case 1. Harvesting the predator only

Figure 11: Yield curves (left) and resilience (right) for harvesting the
predator only given by 25, with varying values of $q_{2}$

When harvesting the striped bass only, the resilience curves increase linearly
until they reach the star at an effort level dependent on $q_{2}$, and then
rapidly drop to 0 as illustrated in Figure 11.The yield curves are all
downward facing non-negative parabolas. As the resilience curves increase
initially, resilience can be maximized by harvesting at the corresponding
effort levels, at $E_{RMY}$. This means that any effort up to $E_{RMY}$ will
increase the system’s biological stability, and any effort past $E_{RMY}$ will
begin to harm it. Note that all of the $E_{RMY}$’s (effort which maximizes the
resilience) are greater than their $E_{R_{m}}$ counterparts (effort that
maximizes the yield).

#### Case 2. Harvesting the prey only

Figure 12: Yield curves (left) and resilience (right) for harvesting the prey
only given by 26, with varying values of $q_{1}$

When harvesting the menhaden only, the resilience curves are strictly
decreasing as seen from Figure 12, with corresponding downward facing non-
negative parabola yield curves. As the resilience curves are strictly
decreasing, resilience is maximized at $E_{1}=0$, showing that any additional
effort on harvesting the prey harms the system’s resilience. However,
additional yield and profit can be obtained at the sacrifice of the system’s
biological stability.

#### Case 3. Combined harvesting of both the prey and the predator

Figure 13: Yield curves (left) and resilience (right) for harvesting both
species given by 27, with varying values of $q_{2}$ when $q_{1}=1,0.5$ Figure
14: Yield curves (left) and resilience (right) for harvesting both species
given by 27, with varying values of $q_{1}$ when $q_{2}=1,0.5$

When harvesting both menhaden and bass at the same effort level, we see that
the behavior of the resilience and yield curves are no longer uniform with
$q_{i}$. In Figure 13, all of the resilience curves are strictly decreasing.
Hence, to maximize the resilience the effort must be zero $E_{RMY}$. When
$q_{1}=1$, the yield parabolas are strictly negative when considering
additional effort for the given $q_{2}$ values, signifying that additional
yield cannot be gained, and that no further effort towards harvesting both
species should be applied as it will decrease the resilience of the prey-
predator system. When $q_{1}=0.5$, some additional yield could be gained at
the expense of the resilience of the system for certain $q_{2}$ values. In
Figure 14, $E_{RMY}$ is only non zero when $q_{2}$ is high and $q_{1}$ is low,
meaning these are the conditions on $q_{i}$ in which additional harvesting
should be done to increase resilience. Here, $E_{RMY}<E_{R_{m}}$, meaning that
the effort needed for gaining maximum resilience is less than the effort
needed to gain maximum yield. In other words, additional effort towards
simultaneous harvesting results in an increase in yield at the cost of
resilience.

##  5 Sensitivity analysis

One important question to pose is if the solutions to the system of
differential equations are stable with respect to small perturbations of the
parameters of the system. To answer this question we perform sensitivity
analysis. Such analysis gives an idea of how small changes in the parameters
of the system affect the behavior of the solutions in time. Here we use the
direct differential method which is a local method for sensitivity analysis.

Assume that the mathematical model, consisting of the first two differential
equations in (3) has the form

$\frac{du}{dt}=f(u,p),$ (37)

where $u=u(t,p)=[x(t,p)\quad y(t,p)]$ is the two component vector solution of
the system, $p=(p_{1},\dots,p_{n})$ is the vector of all parameters, and the
$f$ is the two component vector of the right-hand sides. We are looking to
find the rate of change of $u$ with respect to $p_{j}$, i.e. $\dfrac{\partial
u}{\partial p_{j}}$. We can find the time evolution of this quantity by using
the given differential equation as follows

$\frac{d}{dt}\left(\frac{\partial u_{i}}{\partial p_{j}}\right)=\frac{\partial
f_{i}}{\partial u_{i}}\frac{\partial u_{i}}{\partial p_{j}}+\frac{\partial
f_{i}}{\partial p_{j}}=J\cdot S_{j}+F_{j},$ (38)

where matrix $J$, the Jacobian, is given by

$J=\left(\begin{array}[]{cccccc}\frac{\partial f_{1}}{\partial
u_{1}}&\frac{\partial f_{1}}{\partial u_{2}}&...&\frac{\partial
f_{1}}{\partial u_{k}}\\\ \frac{\partial f_{2}}{\partial u_{1}}&\frac{\partial
f_{2}}{\partial u_{2}}&...&\frac{\partial f_{2}}{\partial u_{k}}\\\
...&...&...&...\\\ \frac{\partial f_{k}}{\partial u_{1}}&\frac{\partial
f_{k}}{\partial u_{2}}&...&\frac{\partial f_{k}}{\partial u_{k}}\\\
\end{array}\right),$

vector $S_{j}$, the sensitivity vector for parameter $p_{j}$, and vector
$F_{j}$, the derivative of the right-hand side of the system with respect to
the parameter $p_{j}$, are given by the vectors:

$S_{j}=\left(\frac{\partial y_{i}}{\partial
p_{j}}\right)=\left(\begin{array}[]{c}\frac{\partial y_{1}}{\partial p_{j}}\\\
\frac{\partial y_{2}}{\partial p_{j}}\\\ ...\\\ \frac{\partial y_{k}}{\partial
p_{j}}\end{array}\right),F_{j}=\left(\begin{array}[]{c}\frac{\partial
f_{1}}{\partial p_{j}}\\\ \frac{\partial f_{2}}{\partial p_{j}}\\\ ...\\\
\frac{\partial f_{k}}{\partial p_{j}}\end{array}\right).$

The Jacobian for our system from (3) is given by the matrix

$J(x,y)=\left(\begin{array}[]{ccc}r-2ax-by-q_{1}E_{1}&-bx\\\ dy&-e-2cy+dx-
q_{2}E_{2}\end{array}\right).$

As the system (3) has 8 different parameters, let us arrange all parameters as
a vector
$p=\left(r\>\>\>a\>\>\>b\>\>\>E_{1}\>\>\>e\>\>\>c\>\>\>d\>\>\>E_{2}\right).$
Then the full sensitivity matrix $S$ is given by

$S=\left(\begin{array}[]{cccccccc}\frac{\partial x}{\partial r}&\frac{\partial
x}{\partial a}&\frac{\partial x}{\partial b}&\frac{\partial x}{\partial
E_{1}}&\frac{\partial x}{\partial e}&\frac{\partial x}{\partial
c}&\frac{\partial x}{\partial d}&\frac{\partial x}{\partial E_{2}}\\\
\frac{\partial y}{\partial r}&\frac{\partial y}{\partial a}&\frac{\partial
y}{\partial b}&\frac{\partial y}{\partial E_{1}}&\frac{\partial y}{\partial
e}&\frac{\partial y}{\partial c}&\frac{\partial y}{\partial d}&\frac{\partial
y}{\partial E_{2}}\\\ \end{array}\right).$ (39)

And the full derivative of the right-hand side with respect to the parameters
is given by the matrix

$F=\left(\begin{array}[]{cccccccccccccc}x&-x^{2}&-xy&-x&0&0&0&0\\\
0&0&0&0&-y&-y^{2}&xy&-y\\\ \end{array}\right).$ (40)

Hence, for example, the system of differential equations that investigates the
rate of change of the dependent variables $x$ and $y$ with respect to the
parameter $r$ is given by four differential equations:

$\left\\{\begin{array}[]{lll}\dfrac{dx}{dt}&=&x(r-ax-by)-q_{1}E_{1}x\\\
\dfrac{dy}{dt}&=&y(-e-cy+dx)-q_{2}E_{2}y\\\ \dfrac{ds_{11}}{dt}&=&(r-2ax-by-
E_{1})s_{11}-bxs_{21}+x\\\ \dfrac{ds_{21}}{dt}&=&dys_{11}-(e+2cy-dx-
q_{2}E_{2})s_{21}\\\ \end{array}\right.$

Observe that the first two differential equations are exactly the differential
equations in the system (3) while the third and fourth equations come from
(38) by using the Jacobian and the first column of (40). Since different
parameters can have different units and thus can have different orders of
magnitude, we calculate the relative sensitivity matrix
$\bar{S}_{ij}=\dfrac{\partial u_{i}}{\partial p_{j}}\dfrac{p_{j}}{u_{i}}$.
Then the sensitivity of each parameter $p_{j}$ is given by the sensitivity
index defined as the magnitude of the corresponding norm of the $i,j$-th
element in $\bar{S}$:

$\left\|\bar{S}_{ij}\right\|_{2}=\sqrt{\sum_{k=1}^{N}\left(\dfrac{\partial
u_{i}}{\partial p_{j}}(t_{k})\dfrac{p_{j}}{u_{i}(t_{k})}\right)^{2}}.$

For the parameter set given in Table 1, we computed the relative sensitivity
vector for each of the steady-states of the variables $x(t)$ and $y(t)$ (for
$t>50$) which are given in descending order below.

Menhaden’s Norms | Bass’s Norms
---|---
| Parameter p | $\left\|\frac{dx}{dp}\right\|_{2}$
---|---
d | 0.7932
e | 0.6523
r | 0.3477
a | 0.2068
b, c | 0.1409
$E_{2}$ | $1.6430E_{2}$
$E_{1}$ | $0.8359E_{1}$
| Parameter p | $\left\|\frac{dy}{dp}\right\|_{2}$
---|---
r | 1.9572
a, d | 1.1640
e | 0.9572
b | 0.7932
c | 0.2068
$E_{1}$ | $4.7049E_{1}$
$E_{2}$ | $2.4111E_{2}$

Figure 15 illustrates the sensitivity of striped bass (the top graph) and
menhaden (bottom graph) in respect to each of the eight parameters by showing
the rate of change of each variable in respect to the parameter over time.
From the sensitivity tables and Figure 15 we see that the striped bass is most
sensitive to perturbations in the intrinsic growth rate of the menhaden $r$.
As illustrated in Figure 16, even small increase in the growth rate of
menhaden increases the steady-state of bass significantly from about 0.3 to
0.46. From the sensitivity table for menhaden we see that it is the most
sensitive to changes in the parameter $d$, the preying effect of menhaden on
bass. A marginal increase of $d$ decreases significantly the steady-state of
menhaden but increases the steady-state of bass’ population. The next most
sensitive parameter for menhaden is the growth rate of bass $e$. Even a
marginal increase of $e$, as shown in Figure 16, results in higher steady-
state for menhaden but lower steady-state for the bass.

Other important parameters as identified from the sensitivity tables are the
harvesting coefficients $E_{1}$ and $E_{2}$. The predator (striped bass) is
more sensitive to the prey’s harvesting effort $E_{1}$ while the prey
(menhaden) is more sensitive to the predator’s harvesting effort $E_{2}$. As
illustrated in Figure 16, a marginal increase of ${E_{2}}$ results in an
increase of the prey’s steady-state, and a decrease in the predator’s steady-
state; further, it affects the predator’s population more than the prey’s
population, meaning that the predator’s population is more sensitive to
changes in the harvesting of the predator. Also in Figure 16 we can see that a
marginal increase of ${E_{1}}$ results in a small decrease in the steady-state
of the prey’s population but results in a much serious decrease in the steady-
state of predator’s population. This means that a reduction in the harvesting
of the prey (menhaden) would be also beneficial to the predator’s (striped
bass) long term sustainability.

Figure 15: The rate of change of each species’ population with respect to the
system parameters over time.

##  6 Conclusions

The Atlantic menhaden and striped bass are two of the most important fish in
the Chesapeake Bay. Yet, both of these fish populations are at risk throughout
the last few decades due to several issues, one of the biggest being
overharvesting. There is still uncertainty among policymakers about the
management approach which would preserve these species. Here we applied
mathematical modeling to study the current predator-prey dynamics between the
species, and particularly the long-term effects of three harvesting policies:
maximum sustainable yield (MSY), maximum economic yield (MEY), and resilience
maximizing yield (RMY). Specifically, we examine how each harvesting policy
affects the system’s dynamics, and explore the trade-off between maximizing
economic yield and maintaining long-term sustainability within each policy.

We used the Lotka-Volterra predator-prey model to analyze the population
dynamics of menhaden and bass and the possibility of balancing economic gains
with ecosystem conservation. We began by fitting the generalized Lotka-
Volterra model to the time series data of the menhaden and bass’s numerical
indices to find the model’s parameter values and calculate the resilience of
the system. Next, we included a harvesting effort component into the model and
found the bioeconomic equilibria, as well as derived conditions for their
stability. We have proven the positivity and boundedness of the solutions
within an invariant region. This shows that the model is well posed, and if
the initial conditions of the system are within this region, then all
solutions will remain there. We studied the effects of maximum sustainable
yield (MSY), maximum economic yield (MEY), and resilience maximizing yield
(RMY) policies on the sustainability of the fish populations.

Our analysis has shown that for all three harvesting options, harvesting the
predator only, the prey only, or both the predator and prey, additional
harvesting at MSY and MEY results in an ecologically sustainable system.
Notably, only the case of harvesting the predator only will allow for being
able to harvest at $E_{Y_{m}}$, resulting in achieving the maximum theoretical
yield, as the other two harvesting cases, prey only and both prey and
predator, will have the predator reach extinction from implementing
$E_{Y_{m}}$. Even though relatively sustainable, MSY policy does not result in
a higher yield compared to MEY policy; hence, it is not economically justified
as it does not increase the profit, but cuts the bass population at
equilibrium by at least half in all harvesting cases. When harvesting the prey
only or both prey and predator, the equilibrium population of the striped bass
decreases at much higher rate compared to menhaden with respect to $E$.
Harvesting both menhaden and bass at MEY levels is efficient and balances the
economic profit and sustainability of the two species, but can have varying
effects on the system’s resilience to perturbations in their populations.

As currently ongoing harvesting effort is already included in the $r$
parameter value (the intrinsic growth rate) of the fitted model, any
additional harvesting $E$ will change the resilience, depending on the
harvesting case, and further by the values of $q_{i}$. When harvesting the
predator only, resilience is initially increasing with $E$, showing that
additional effort on harvesting the predator can increase the system’s
resilience, while also providing yield. When harvesting the prey only,
resilience is maximized when the effort is strictly zero, showing that
additional effort on harvesting the prey harms the system’s resilience. In
this case, additional yield and profit can be obtained by sacrificing of the
system’s biological stability. When harvesting both the predator and prey
simultaneously, the resilience curves and corresponding analysis become
dependent on $q_{i}$, and similar scenarios to the first two cases are
observed with different $q_{i}$ values.

Figure 16: State variables, menhaden in red and bass in blue, time evolution
for small perturbations in the model parameters values from Table 1.

Finally, we performed sensitivity analysis of the mathematical model with
respect to small perturbations in its parameters. Sensitivity analysis
revealed that the bass is most sensitive to the intrinsic growth rate $r$ of
menhaden and even a marginal increase of the growth rate of the prey would
benefit the predator and increase significantly its steady-state on a long
run. The sensitivity analysis also has shown that both menhaden and bass are
very sensitive to the preying effect coefficient $d$. Even a marginal increase
in $d$ increases the steady-state of the bass but decreases the steady-state
of menhaden. Our sensitivity analysis also shows that harvesting the prey
affects much more negatively the predator than the prey while harvesting the
predator decreases the predator’s population but simultaneously increases
significantly the prey’s population. These sensitivity results shed light on
the effects of the system’s parameters on the long-term population dynamics of
menhaden and bass and in combination with the results on the different
harvesting policies may be very beneficial for policy-makers while deciding on
future fishing quotes for these species.

## Acknowledgments

We wish to acknowledge the support of the NSF-CURM grant DMS-1722563.

## References

  * [1] Ahrenholz, D.W., 1991. Population Biology and Life History of the North American Menhadens, _Brevoortia_ spp. Marine Fisheries Review. 53(4), 3-19. http://hdl.handle.net/1834/26530
  * [2] Ahrenholz, D.W., et al., 1987. Population and Fishery Characteristics of Atlantic Menhaden, _Brevoortia tyrannus_. Fishery Bulletin. 85(3), 569-600.
  * [3] Casas-Valdez, M., et al., 2005. Estimation of maximum sustainable yield of _Gelidium robustum_ seaweed fishery in Mexico. Journal of the Marine Biological Association of the United Kingdom. 85(4), 775-778. https://doi.org/10.1017/S0025315405011689
  * [4] Chesapeake Bay Foundation. Atlantic Menhaden. https://www.cbf.org/about-the-bay/more-than-just-the-bay/chesapeake-wildlife/menhaden/ (accessed 13 November 2021)
  * [5] Chesapeake Bay Foundation. Striped Bass (Rockfish). https://www.cbf.org/about-the-bay/more-than-just-the-bay/chesapeake-wildlife/rockfish/ (accessed 20 November 2021)
  * [6] Chesapeake Bay Foundation. The Issues. https://www.cbf.org/issues/ (accessed 24 October 2021)
  * [7] Chesapeake Bay Program. Atlantic Menhaden. https://www.chesapeakebay.net/discover/field-guide/entry/atlantic_menhaden (accessed 13 November 2021)
  * [8] Chesapeake Bay Program. Learn the Issues. https://www.chesapeakebay.net/issues (accessed 24 October 2021)
  * [9] Chesapeake Bay Program. Learn the Issues: Striped Bass. https://www.chesapeakebay.net/issues/striped_bass (accessed 24 October 2021)
  * [10] Chesapeake Bay Program. Striped Bass. https://www.chesapeakebay.net/S=0/fieldguide/critter/striped_bass (accessed 24 October 2021)
  * [11] Durell, E.Q., Weedon, C. 2020. Striped Bass Seine Survey Juvenile Index Web Page. Maryland Department of Natural Resources. https://dnr.maryland.gov/fisheries/pages/striped-bass/juvenile-index.aspx
  * [12] Ellis, Julia K., 2003. Diet of the Sandbar Shark, Carcharhinus plumbeus, in Chesapeake Bay and Adjacent Waters. Dissertations, Theses, and Masters Projects. William and Mary. Paper 1539617797. https://dx.doi.org/doi:10.25773/v5-h091-j251
  * [13] Fath, Brian. 2018. Encyclopedia of Ecology, second ed. Elsevier, pp. 108-115.
  * [14] Glass, A.K., Watts, B.D., 2009. Osprey diet composition and quality in high- and low-salinity areas of lower Chesapeake Bay. The Journal of Raptor Research. 43(1), 27-36. https://doi.org/10.3356/JRR-08-22.1
  * [15] Harding, J.M., Mann, R., 2001. Diet and habitat use by bluefish, _Pomatomus saltatrix_ , in a Chesapeake Bay estuary. Environmental Biology of Fisheries. 60, 401-409. https://doi.org/10.1023/A:1011071127930
  * [16] June, F.C., Carlson, F.T., 1970. Food of young Atlantic Menhaden, _Brevoortia tyrannus_ , in relation to metamorphosis. Fishery Bulletin. 68(3).
  * [17] Karim, E., et al., 2020. Maximum Sustainable Yield Estimates of Marine Captured Shrimp Fishery of the Bay of Bengal, Bangladesh by Using Surplus Production Model. Thalassas: An International Journal of Marine Sciences. 36, 471-480. https://doi.org/10.1007/s41208-020-00198-9
  * [18] Manooch III, C.S., 1973. Food Habits of Yearling and Adult Striped Bass, Morone saxatilis (Walbaum), from Albemarle Sound, North Carolina. Chesapeake Science. 14(2), 73-86. https://doi.org/10.2307/1350872
  * [19] Maunder, M.N., 2002. The relationship between fishing methods, fisheries management and the estimation of maximum sustainable yield. Fish and Fisheries. 3(4), 251-260. https://doi.org/10.1046/j.1467-2979.2002.00089.x
  * [20] Merriner, J.V., 1975. Food Habits of the weakfish, Cynoscion regalis, in North Carolina waters. Chesapeake Science. 16(1,) 74-76. https://doi.org/10.2307/1351090
  * [21] Mühlbauer, L.K., et al., 2020. gauseR: Simple methods for fitting Lotka-Volterra models describing Gause’s ”Struggle for Existence”. Ecology and Evolution. 10(23), 13275-13283. https://doi.org/10.1002/ece3.6926
  * [22] National Oceanic and Atmospheric Administration. Atlantic Menhaden. https://www.fisheries.noaa.gov/species/atlantic-menhaden (accessed 14 January 2022)
  * [23] National Park Service. Environmental Factors. https://www.nps.gov/chba/learn/environmental-factors.htm (accessed 14 January 2022)
  * [24] Northeast Fisheries Science Center (NEFSC). 2019. 66th Northeast Regional Stock Assessment Workshop (66th SAW) Assessment Report. US Dept Commer, Northeast Fish Sci Cent Ref Doc. 19-08; 1170 p. Available from: http://www.nefsc.noaa.gov/publications/
  * [25] Pimm, S.L., J.H. Lawton, 1978. On feeding on more than one trophic level. Nature. 275, 542-544. https://doi.org/10.1038/275542a0
  * [26] Ress, D. 2020. Chesapeake Bay’s menhaden catch cut drastically, along with Omega Protein quota. Daily Press. https://www.dailypress.com/news/dp-nw-menhaden-20200428-gxtukdev4vd5rmyff7xlr7grhm-story.html
  * [27] Robichaud-LeBlanc, K.A., et al., 1996. Spawning and early life history of a northern population of striped bass (_Morone saxatilis_) in the Miramichi River estuary, Gulf of St. Lawrence. Canadian Journal of Zoology. 74(9), 1645-1655. https://doi.org/10.1139/z96-182
  * [28] Schaefer, Milner B., 1954. Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Inter-American Tropical Tuna Commission. 1(2), 23-56. http://hdl.handle.net/1834/21257
  * [29] SEDAR. 2020. SEDAR 69 – Atlantic Menhaden Benchmark Stock Assessment Report. SEDAR, North Charleston SC. 691 pp. available online at: http://sedarweb.org/sedar-69
  * [30] Shelton, P.A., Morgan, M.J., 2014. Impact of maximum sustainable yield-based fisheries management frameworks on rebuilding North Atlantic cod stocks. Journal of Northwest Atlantic Fishery Science. 46, 15-25. https://doi:10.2960/J.v46.m697
  * [31] Shibata, Y., Nagao, J., et al., 2020. Estimating the maximum sustainable yield of snow crab (_Chionoecetes opilio_) off Tohoku, Japan via a state-space stock assessement model with time-varying natural mortality. Population Ecology. 63(1), 41-60. https://doi.org/10.1002/1438-390X.12068
  * [32] Siyal, F.K., et al., 2013. Maximum Sustainable Yield Estimates of Silver Pomfret, _Pampus argenteus_ (Family: Strometidae) Fishery in Pakistan. Pakistan Journal of Zoology. 45(2), 447-452.
  * [33] Southwick Associates, The Economic Contributions of Recreational and Commercial Striped Bass Fishing, 2019-04, https://mcgrawconservation.org/wp-content/uploads/McGraw-Striped-Bass-Report-FINAL.pdf
  * [34] The National Wildlife Federation. Protecting the Chesapeake Bay. https://www.nwf.org/Our-Work/Waters/Great-Waters-Restoration/Chesapeake-Bay. (accessed January 15 2022)
  * [35] Theodore Roosevelt Conservation Partnership, 2019. Omega Protein Breaks the Rules in Chesapeake Bay. https://www.trcp.org/2019/09/19/omega-protein-breaks-rules-chesapeake-bay/ (accessed 11 March 2022)
  * [36] Thorpe, R.B., 2019. What is a multispecies MSY? A worked example from the North Sea. Journal of Fish Biology. 94(6), 1011-1018. https://doi.org/10.1111/jfb.13967
  * [37] Tolliver, L. 2019. Menhaden fishery to be shut down in Virginia if company continues overfishing, secretary of commerce says. The Virginia Pilot. https://www.pilotonline.com/life/wildlife-nature/vp-nw-menhaden-shutdown-20191219-x6i6af3karbnnn5abksmlrzcku-story.html (accessed 19 March 2022)
  * [38] Ulrich, C., et al., 2017. Achieving maximum sustainable yield in mixed fisheries: a management approach for the North Sea demersal fisheries. ICES Journal of Marine Science. 74(2), 566-575. https://doi.org/10.1093/icesjms/fsw126
  * [39] US Fish and Wildlife Service. Striped Bass. https://www.fws.gov/species/maryland-bass-hybrid-morone-saxatilis (accessed 24 October 2021)
  * [40] [dataset] Virginia Institute of Marine Science. Species of Interest: Atlantic Menhaden. NEAMAP Abundance Index. https://www.vims.edu/research/departments/fisheries/programs/multispecies_fisheries_research/speciesofinterest/atlantic_menhaden.php
  * [41] [dataset] Virginia Institute of Marine Science. Species of Interest: Striped Bass. NEAMAP Abundance Index. https://www.vims.edu/research/departments/fisheries/programs/multispecies_fisheries_research/speciesofinterest/striped-bass.php
  * [42] Walter III, J.F., Austin, H.M., 2003. Diet composition of large striped bass (_Morone saxatilis_) in Chesapeake Bay. Fishery Bulletin. 101(2), 414-423. https://scholarworks.wm.edu/vimsarticles/572


Proficiency Test: Someone attaches the pen drive / camera to the cd player.
Main Test: Someone attaches / detaches the pen drive to the cd player.
Proficiency Test: Someone reveals the battery / truth.
Main Test: Someone reveals / covers the battery.
Figure 16: Sample instances from the Change of State (Action) tests.
---
Proficiency Test: The skier / jet skies.
Main Test: Initially, the skier is in a higher / lower position.
Proficiency Test: Someone digs the seeds / body out.
Main Test: Initially, the seeds are inside / outside.
Proficiency Test: Someone empties the cup / trash.
Main Test: Initially, the cup is full / empty.
Proficiency Test: Someone pulls the spray / phone out.
Main Test: Initially, the spray is inside / outside.
Proficiency Test: Someone pulls a baby baby toy / gun.
Main Test: Initially, a baby toy is further away / closer.
Figure 17: Sample instances from the Change of State (Pre-State) test.
---
Proficiency Test: Someone connects the chair and the base / dots.
Main Test: At the end, the chair and the base is attached / detached.
Proficiency Test: Someone folds the paper / laundry.
Main Test: At the end, the paper is folded / unfolded.
Proficiency Test: Someone uncovers the bulb / truth.
Main Test: At the end, the bulb is visible / hidden.
Proficiency Test: Someone plugs the usb / plug.
Main Test: At the end, the usb is inside / outside.
Proficiency Test: Someone soaks the wakame / car.
Main Test: At the end, the wakame is wet / dry.
Figure 18: Sample instances from the Change of State (Post-State) test.
---
Proficiency Test: Someone wraps the pizza / presents.
Main Test: Initially, the pizza is unwrapped / wrapped. Then, someone wraps /
unwraps the pizza. At the end, the pizza is wrapped / unwrapped.
Proficiency Test: Someone fills the glass / gap.
Main Test: Initially, the glass is empty / full. Then, someone fills / empties
the glass. At the end, the glass is full / empty.
Proficiency Test: Someone chops an onion / apple.
Main Test: Initially, an onion is whole / in pieces. Then, someone chops /
connects an onion. At the end, an onion is in pieces / whole.
Proficiency Test: Someone wipes the screen / floor.
Main Test: Initially, the screen is dirty / clean. Then, someone wipes /
stains the screen. At the end, the screen is clean / dirty.
Proficiency Test: Someone fixes the new string / problem.
Main Test: Initially, the new string is separated / attached. Then, someone
fixes / detaches the new string. At the end, the new string is attached /
separated.
Figure 19: Sample instances from the Change of State (Reverse) test.

###  C.5 Rare Actions

In the Rare Actions test, we investigate the ability of VidLMs to identify
novel compositions and recognise unusual events, such as a computer keyboard
is being cut using a chainsaw by someone described. These events are described
by a verb-noun pair, e.g. “cutting a keyboard". We choose foils from more
likely events taking place in videos.

#### Data sources.

We leverage the RareAct dataset (Miech et al., 2020), which consists of videos
accompanied by action-object pairs describing events within the videos. These
action-object pairs are extracted by analysing co-occurrence statistics from
the widely used HowTo100M (Miech et al., 2019) dataset for VidLM pretraining.
To enrich this dataset, we generate simple captions based on the action-object
pairs. For instance, given the action-object pair cut-keyboard, we create the
descriptive caption cutting a keyboard. We avoid subject phrases in captions
since some videos do not comprise a human being as an actor.

#### Foiling method.

This test offers two subtests: the Action Replacement and the Object
Replacement subtest. In the Action Replacement subtest, we substitute the
original action with a more plausible alternative that can be applied to the
given object, e.g. type on for the previous keyboard example. To generate
foils in this subtest, we employ T5 (Raffel et al., 2020), as it can produce
foil candidates with compound verbs, e.g., talk on, place at, etc. We discard
foil candidates with some general actions (e.g. use, have etc.) or actions
that imply some form of touching (e.g. hold, reach etc.). We then perform an
NLI and GRUEN score filtering. In this stage, we perform a manual intervention
and abandon the low quality candidates. As for the Object replacement subtest,
we focus on replacing the object in the action-object pair. For instance,
revisiting the previous example, we replace the object keyboard with bread.
Here, we prefer to use a set of token-based MLMs (Devlin et al., 2019; Lan et
al., 2020; Liu et al., 2019). To further enhance the quality of the foils, we
opt for an ensembling approach in the object replacement test. Particularly,
we use three MLMs which are BERT, RoBERTa and ALBERT (Devlin et al., 2019; Liu
et al., 2019; Lan et al., 2020). We also run an object detector called End-to-
End Object Detection model (DETR) (Carion et al., 2020) and discard the foil
candidates that contain the detected objects. We use `facebook/detr-
resnet-101` DETR in Huggingface’s transformers package (Wolf et al., 2019). We
uniformly sample $K=8$ frames and run the object detector on these sampled
frames. We classify an object detected if its confidence threshold exceeds
$0.80$.

#### Proficiency Test.

The proficiency test of the Rare Actions focuses on the existence of objects
in the given input videos, similar to the VALSE existence instrument. This
time we do not use negated foils: We replace the correct object with another.
To create captions, we create a statement about the existence of the ground
truth object, e.g. “there is at least one keyboard.". To create foils, we
randomly sample from the objects appear in ground-truth captions, e.g. “there
are some flowers". Similar to the main test, we implement the same object
detection filtering process.

#### In-Depth Results.

Table 9 presents the model performance achieved on the Rare Actions test.

Table 9: Rare Actions subtest results using pairwise accuracy ($acc_{r}$) metric. P, T and P+T stand for the scores achieved on the proficiency tests, the main tests only and the combined tests, respectively. Model | Action Replacement | Object Replacement | All
---|---|---|---
P | T | P+T | P | T | P+T | P | T | P+T
Random | 50.00 | 50.00 | 25.00 | 50.00 | 50.00 | 25.00 | 50.00 | 50.00 | 25.00
GPT-2 | 61.25 | 17.98 | 9.85 | 55.20 | 34.39 | 26.16 | 58.35 | 25.85 | 17.67
OPT | 60.99 | 20.51 | 10.79 | 56.79 | 27.60 | 19.36 | 58.97 | 23.91 | 14.90
CLIP | 92.81 | 93.74 | 86.95 | 92.63 | 94.08 | 88.73 | 92.72 | 93.90 | 87.80
BLIP2 | 93.21 | 64.98 | 60.19 | 94.51 | 84.83 | 81.65 | 93.83 | 74.50 | 70.48
ClipBERT | 38.75 | 27.96 | 7.99 | 40.75 | 52.60 | 20.81 | 39.71 | 39.78 | 14.14
UniVL | 77.23 | 78.83 | 60.59 | 77.75 | 77.17 | 59.25 | 77.48 | 78.04 | 59.89
VideoCLIP | 84.15 | 75.10 | 62.98 | 83.82 | 80.64 | 72.40 | 83.99 | 77.75 | 67.5
FiT | 89.21 | 86.55 | 76.30 | 90.17 | 92.49 | 85.55 | 89.67 | 89.40 | 80.73
CLIP4Clip | 82.95 | 93.07 | 78.16 | 82.94 | 95.23 | 79.19 | 82.95 | 94.10 | 78.65
VIOLET | 86.68 | 85.49 | 73.50 | 87.57 | 87.86 | 75.87 | 87.10 | 86.63 | 74.64
X-CLIP | 82.69 | 86.28 | 70.57 | 85.12 | 84.97 | 74.13 | 83.85 | 85.65 | 72.28
MCQ | 91.48 | 86.82 | 79.09 | 91.18 | 90.75 | 85.69 | 91.34 | 88.70 | 82.26
Singularity | 92.14 | 85.22 | 78.16 | 93.21 | 92.77 | 88.44 | 92.65 | 88.84 | 83.09
UniPerceiver | 57.12 | 55.12 | 31.02 | 59.39 | 62.71 | 38.72 | 58.21 | 58.76 | 34.71
Merlot Reserve | 83.75 | 96.27 | 80.42 | 83.81 | 84.39 | 74.56 | 83.78 | 80.57 | 77.61
VindLU | 94.01 | 92.14 | 86.42 | 94.36 | 94.08 | 89.60 | 94.18 | 93.07 | 87.94

##### Unimodal Results.

Unimodal baselines GPT-2 and OPT (Radford et al., 2019; Zhang et al., 2022)
perform very poorly on the main tests as expected since the captions describe
less likely events. We observe that the models consistently perform better in
object replacement tests in comparison to action replacement. This is inline
with previous work where the models are more biased towards nouns and they
fail to process verbs sufficiently (Momeni et al., 2023; Park et al., 2022;
Lei et al., 2022).

##### Image-Language Model Results.

CLIP consistently outperforms random in all individual and combined tests,
showcasing its effectiveness in understanding and executing the given tasks.
In contrast, BLIP2, while excelling in proficiency tests, demonstrates lower
performance in the main test and combined test compared to CLIP, implying
potential limitations or challenges in action replacement subtest for BLIP2.
Additionally, CLIP outperforms all VidLMs except VindLU, which demonstrates
its ability to handle novel compositions is superior to these models.

##### Video-Language Model Results.

All VidLMs, consistently outperform random in all individual and combined
tests, showcasing their effectiveness in understanding and performing the
given tasks. The high accuracy scores across proficiency, main, and combined
tests suggest the capability of these models in achieving accurate results in
Rare Actions subtests. However, UniPerceiver struggles with accurately
comprehending and replacing actions or objects in the given scenarios compared
to other VidLMs, potentially indicating a need for improvements in its ability
to understand nuanced visual and contextual information.

#### Test Examples.

In Figure 20 and Figure 21, we show some sample validated examples from the
Rare Actions – Action Replacement and Object Replacement subtests.

---
Proficiency Test: there is at least one fridge / chocolate
Main Test: washing some fridges / eating from a fridge
Proficiency Test: there is at least one banana / blender
Main Test: weighing / eating a banana
Proficiency Test: there is at least one shirt / are some peppers
Main Test: rolling / putting a shirt
Proficiency Test: there are some spoons / is at least one car
Main Test: throwing / serving with some spoons
Proficiency Test: there is at least one laptop / car
Main Test: measuring / accessing a laptop
Figure 20: Sample instances from the Rare Actions (Action Replacement) test.
---
Proficiency Test: there are some corns / is at least one blender
Main Test: peeling some corns / a lemon
Proficiency Test: there is at least one phone / egg
Main Test: hammering a phone / some nails
Proficiency Test: there is at least one door / towel
Main Test: spraying a door / an area
Proficiency Test: there is at least one phone / bicycle
Main Test: microwaving a phone / some food
Proficiency Test: there is at least one chocolate / are some corns
Main Test: drinking a chocolate / some tea
Figure 21: Sample instances from the Rare Actions (Object Replacement) test.

###  C.6 Spatial Relations

In the Spatial Relations test, we investigate the capabilities of VidLMs to
understand spatial as well as spatio-temporal relations in a video (e.g.
moving something towards or from something).

#### Data sources

We create the foils starting from the Something-Something V2 dataset Goyal et
al. (2017a), which contains a collection of $220,847$ labelled clips of 174
pre-defined basic actions performed by humans with common objects, such as
putting something into something or turning something upside down. Since this
dataset is often used to pretrain VidLMs models, we leverage only the video-
caption pairs present in the validation set.

#### Foiling method.

In order to create a foil, we generate several candidates for each caption by
replacing a preposition in the caption with others drawn from the set of
prepositions in the validation set itself. This ensures that foils and
captions express relations with prepositions from the same distribution. Each
candidate is scored for perplexity using
T5-base141414https://huggingface.co/t5-base. We test also with T5-large, but
we obtain similar results, therefore we use the base version for better
efficiency.

In order to facilitate the model in the scoring stage, we prepend a subject to
the candidate before feeding it into the model. Therefore a sentence like,
rolling a can onto a flat surface becomes someone rolling a can onto a flat
surface. We observe that with this format, the model produces better
perplexity scores, probably because the data T5 was trained on tends to
contain full sentences with explicit subjects, rather than verb phrases only.
However, this is carried out only for scoring purposes. Once a candidate is
selected, the version included in the main test does not have any prepended
subject.

We select the top $10$ best-scoring candidates for each sample and we run them
through a state-of-the-art NLI model based on
RoBERTa151515https://huggingface.co/ynie/roberta-large-
snli_mnli_fever_anli_R1_R2_R3-nli, using the caption as premise and the
candidate foil as a hypothesis. Similarly to Parcalabescu et al. (2022), we
keep only the candidates predicted as neutral or contradiction. The candidate
foils are checked for grammaticality, by computing the GRUEN score and
filtering out candidates with scores lower than $0.6$. Finally, we adjust the
distribution of prepositions in the foils to match the caption distribution.
This is achieved by removing foil candidates containing prepositions with
frequencies that exceed the frequency of the caption distribution. This
mitigates distribution biases arising in the foil generation process, which
could be exploited by the tested models.

#### Proficiency Test.

We focus on the object identification task for the proficiency test of Spatial
Relations, similar to Change of State proficiency tests. We use RoBERTa (Liu
et al., 2019) to generate the foil by either masking the subject or the object
of the caption, again, depending on the transitivity of the verb as we did in
Change of State.

#### In-Depth Results.

We report in Table 10 the results for Spatial Relations.

Table 10: Spatial Relations results using pairwise accuracy ($acc_{r}$) metric. P, T and P+T stand for the scores achieved on the proficiency tests, the main tests only and the combined tests, respectively. Model | Spatial Relations
---|---
P | T | P+T
Random | 50.00 | 50.00 | 25.00
GPT-2 | 49.1 | 72.8 | 43.0
OPT | 59.0 | 84.7 | 55.7
CLIP | 78.6 | 58.3 | 44.8
BLIP2 | 91.1 | 86.0 | 79.4
ClipBERT | 44.0 | 65.1 | 30.0
UniVL | 62.5 | 51.7 | 33.2
VideoCLIP | 67.9 | 54.7 | 39.7
FiT | 70.5 | 51.9 | 38.7
CLIP4Clip | 79.8 | 56.7 | 44.2
VIOLET | 73.3 | 50.4 | 38.7
X-CLIP | 74.8 | 56.2 | 43.5
MCQ | 79.4 | 48.9 | 39.4
Singularity | 80.7 | 46.8 | 38.9
UniPerceiver | 45.5 | 48.0 | 20.1
Merlot Reserve | 63.1 | 41.9 | 29.2
VindLU | 83.2 | 45.6 | 39.4

##### Unimodal Results.

Unimodal models, namely GPT-2 and OPT show higher performance on the main test
(T) than on the proficiency test (P) where they are close to random chance.
The quite decent result on the main test (T) suggests that the models can
still exploit some spurious correlations in the text. This may be due to the
design choice to include only in-distribution relations in the foil
construction which although, filtered for grammaticality may still be detected
as less likely than the actual caption. This can certainly be seen as a
potential limitation of the Spatial Relations test. However, this is less
visible in the combined test (P+T) where both perform better than the random
baseline.

##### Image-Language Model Results.

Both CLIP and BLIP2 perform well on the proficiency test (P). This is expected
as this test heavily relies on object identification and it is known that ILMs
are usually trained with object-centric texts. On the other hand, we observe a
drop in performance on the main test (T) for both models and, to a greater
extent on CLIP. This is also expected since the task is designed to exploit
temporal information, which cannot be exploited by these models.

Despite that, BLIP2 keeps performing decently in both the main test (T) and
the combined one (P+T), setting the best performance on the spatial relations
among the evaluated models. This raises questions regarding the capability of
VidLMs to ground visio-temporal and textual information.

##### Video-Language Model Results.

Apart from UniPerceiver and ClipBERT (which perform below chance), all the
VidLMs perform decently on the proficiency test (P). However, their
performance consistently drops for both the main (T) and the combined test
(P+T). CLIP4Clip is the best-performing VidLM however, its performance is far
from both unimodal and ILMs. This suggests a consistent lack of grounding in
VidLMs concerning spatial relations and the need for VidLMs to design better
strategies to properly leverage the temporal information in video-language
tasks.

#### Test Examples.

In Figure 22, we show some sample validated examples from the Spatial
Relations test.

---
Proficiency Test: rolling a can on a flat surface / screen
Main Test: rolling a can on / onto a flat surface
Proficiency Test: pouring water onto a chair / table
Main Test: pouring water onto / out of a chair
Proficiency Test: pulling keys out of a lock / drawer
Main Test: pulling keys out of / on a lock
Proficiency Test: pushing small sunscreen lotion from right to left / bottom
Main Test: pushing small sunscreen lotion from / on right to left
Proficiency Test: holding a pen over scissors / paper
Main Test: holding a pen over / with scissors
Figure 22: Sample instances from the Spatial Relations test. The descriptions
are quite similar except for a small lexical variation (in blue) for the
caption and (in orange) for the foil.