Fish should not be in isolation: Calculating maximum sustainable yield using an ensemble model

Michael A. Spence Khatija Alliji Affiliation: Centre for Environment, Fisheries and Aquaculture Science, Pakefield Road, Lowestoft, Suffolk NR33 0HT, UK Hayley J. Bannister Affiliation: Centre for Environment, Fisheries and Aquaculture Science, Pakefield Road, Lowestoft, Suffolk NR33 0HT, UK Nicola D. Walker Affiliation: Centre for Environment, Fisheries and Aquaculture Science, Pakefield Road, Lowestoft, Suffolk NR33 0HT, UK Angela Muench Affiliation: Centre for Environment, Fisheries and Aquaculture Science, Pakefield Road, Lowestoft, Suffolk NR33 0HT, UK

Abstract

Many jurisdictions have a legal requirement to manage fish stocks to maximum sustainable yield (MSY). Generally, MSY is calculated on a single-species basis, however in reality, the yield of one species depends, not only on its own fishing level, but that of other species. We show that bold assumptions about the effect of interacting species on MSY are made when managing on a single-species basis, often leading to inconsistent and conflicting advice, demonstrating the requirement of a multispecies MSY (MMSY). Although there are several definitions of MMSY, there is no consensus. Furthermore, calculating a MMSY can be difficult as there are many models, of varying complexity, each with their own strengths and weaknesses, and the value if MMSY can be sensitive to the model used. Here, we use an ensemble model to combine different multispecies models, exploiting their individual strengths and quantifying their uncertainties and discrepancies, to calculate a more robust MMSY. We demonstrate this by calculating a MMSY for nine species in the North Sea. We found that it would be impossible to fish at single-species MSY and that MMSY led to higher yields and revenues than current levels.

Running title: Fish should not be in isolation

Keywords: Ensemble modelling; maximum sustainable yield, uncertainty analysis, multispecies modelling; Bayesian statistics; Nash equilibrium; emulators; Ecosystem based fisheries management

1 Introduction

The human population is growing, which has increased the demand for food production and security, which has led to an incompatibility of food production and conservation priorities. Both of which require a balance, to sustainably support an increasing population (Hilborn 2007). Marine fish are a valuable source of food and income for many countries, however the global yield has levelled off and begun to decline, since the 1990s (FAO 2009; Worm et al. 2009). There is now an urgent need to manage fish stocks sustainably so that the balance between food production, conservation and the socio-economics are considered. This will ideally lead to an increase in food production, whilst protecting fish stocks and jobs for future generations (Mesnil 2012).

Fisheries managers use maximum sustainable yield (MSY), which is intended to ensure the sustainability of fish stocks whist maximising food production, without compromising the reproductive potential of the stock (Hilborn 2007; Mesnil 2012). The legal requirement to manage fish stocks to MSY, was adopted in 1982 at the United Nations Convention on the Law of the Sea, the EU Regulation 1380/2013 and the 2002 UN world summit on sustainable development (A/CONF.199/20). The concept of MSY has been adopted by many fisheries management organisations throughout the 1900s, where mathematical and production models were used (Tsikliras & Froese 2018). Recently, MSY has become a widely used reference point in the assessment of fish stocks around the world (Hilborn & Walters 1992; Pauly & Froese 2014; Tsikliras & Froese 2018), here defined as:

Definition 1.

The fishing mortality that leads to the maximum sustainable yield of the $i$ th stock is

F_{MSY,i}(\bm{F}_{-i})={\text{arg}\sup}_{F_{i}}\left(f_{1,i}(F_{i},\bm{F_{-i}})\right),

where $F_{i}$ is the fishing mortality of the $i$ th species, $\bm{F}_{-i}$ is the fishing mortality of the other species and $f_{1,i}(F_{i},\bm{F_{-i}})$ is the $i$ th species’ long-term annual yield (see supplementary material).

MSY has historically been centred on single species MSY (Hart & Fay 2020, SS-MSY,), here defined as:

Definition 2.

The single species fishing mortality that leads to the maximum sustainable yield of the $i$ th stock, the single-species MSY (SS-MSY), is

F_{MSY,i}(\bm{F}_{-i})=F_{MSY,i},

$\forall{}\bm{F}_{-i}$ .

Although MSY is now a widely used concept it often applied to single stocks and when used in this way does not provide information on ecological interactions, resulting in significant ecosystem and fishing ramifications (Andersen et al. 2015; Säterberg et al. 2019). The use of MSY in fisheries management has been criticised for leading to significant changes in community structure, degradation of marine ecosystems and over-exploitation of fisheries resources (Larkin 1977; Hilborn 2007; Andersen et al. 2015). This has rendered MSY policy guidance as incomplete in terms of ecosystem sustainability (Gaichas 2008).

Proposition 1.

SS-MSY exists if and only if

\frac{\partial{}F_{MSY,i}(\bm{F}_{-i})}{\partial{}F_{j}}=0,

$\forall{}j\neq{}i$ .

Proof.

See supplementary material. ∎

Proposition 1 suggests that the fishing mortality on other species does not affect the value of $F_{MSY,i}$ , which does not seem plausible in reality and therefore Definition 2 is rarely met. To counteract this, there has been a push to move to ecosystem-based fisheries management (Pikitch et al. 2004; Link et al. 2011, EBFM,).

Although there is no universally agreed definition of multispecies MSY (Essington & Punt 2011; Norrström et al. 2017, MMSY,), there are a number of proposed alternatives such as the maximum sustainable yield of the community (Andersen 2019). When satisfied, the maximum sustainable yield of the community leads to over-exploitation of species with larger body sizes, leading to a decrease of predating pressure on fish stocks with smaller body sizes (Andersen 2019; Andersen et al. 2015; Szuwalski et al. 2016). Fishing in this way leads to the collapse of stocks of larger species, reducing the diversity of the ecosystem and decreasing the monetary value of the fishery (Andersen 2019; Andersen et al. 2015). Thorpe 2019 and Säterberg et al. 2019 extended this definition to include the risk of species collapse as a caveat to the maximum sustainable community yield. Alternatively the Nash equilibrium has been used to define MMSY (Norrström et al. 2017; Thorpe et al. 2017; Farcas & Rossberg 2016). The Nash equilibrium is a solution to multi-player games where no player can improve their payoff given fixed strategies played by the opponents (Nash 1951). This was applied, with the interest of maximising the yield of each fish stock, whilst taking account of the ecological impacts on other species (Norrström et al. 2017).

Multispecies reference points have been predicted by production models (Sissenwine & Shepherd 1987), however these models ignore the effects of ecological interactions between species (Norrström et al. 2017). Mechanistic multispecies models, henceforth known as simulators, are being increasingly used to support policy decisions, including fisheries and marine environmental polices (Hyder et al. 2015; Nielsen et al. 2018), and are able to capture these interactions. They describe how multiple species interact with their environment and one another, through mechanistic processes allowing them to better predict into the future (Hollowed et al. 2000). However, calculating reference points can be sensitive to the choice of simulator that was used to generate them (Collie et al. 2016; Essington & Plagányi 2013; Fulton et al. 2003; Hart & Fay 2020). This choice can be arbitrary as, although some simulators are better at describing some aspects of the system than others, in general no simulator is uniformly better than the others (Chandler 2013). Instead of choosing one simulator, it is possible to combine them using an ensemble model, allowing managers to maximise the predicting power of the simulators, whilst reducing the uncertainty and errors, improving their decision making. Spence et al. 2018 developed an ensemble model that treats the individual models as exchangeable and coming from a distribution. Their model exploits each simulators strengths, whilst discounting their weaknesses to give a combined solution.

In this paper, we demonstrate how multiple simulators can be combined to find a MMSY. We demonstrate it by finding a Nash equilibrium for nine species in the North Sea using the ensemble model developed in Spence et al. 2018. Although demonstrated with this definition of MMSY in the North Sea, the procedure can be used to optimise any objective function in any environment, including single-species reference points.

2 Methods

We modelled nine species (see Table 1) in the North Sea using historical fishing mortality from 1984 until 2017 (ICES 2018a; ICES 2018c) and fixed fishing mortality, $\bm{F}=(F_{1},\ldots{},F_{9})^{\prime}$ , from 2017 to 2050, with $F_{i}\in[0,2]$ for $i=1,\ldots{},9$ . Our aim was to find $\bm{F}$ values that satisfy the Nash equilibrium (Nash 1951), with a probability that the spawning stock biomass (SSB) falls below $B_{lim}$ , the level of SSB at which recruitment becomes impaired, of 0.25 or less (see Section 4).

Table 1: A summary of the species in the model. The SS-MSY values were taken from ICES 2018a and ICES 2018c.

$i$	Species	Latin name	SS-FMSY	Price per tonne (£)
1	Sandeel	Ammodytes marinus	NA	1314.59
2	Norway pout	Trisopterus esmarkii	NA	151.96
3	Herring	Clupea harengus	0.33	528.34
4	Whiting	Merlangius merlangus	0.15	785.30
5	Sole	Solea solea	0.20	8387.12
6	Plaice	Pleuronectes platessa	0.21	1718.21
7	Haddock	Melanogrammus aeglefinus	0.19	1346.99
8	Cod	Gadus morhua	0.31	1745.22
9	Saithe	Pollachius virens	0.36	855.33

We defined our reference point as:

Definition 3.

$\bm{F}_{Nash}$ , is when

\forall{}i,F_{i}:f_{1,i}(F_{Nash,i},\bm{F}_{Nash,-i})\geq{}f_{1,i}(F_{i},\bm{F}_{Nash,-i})

and $Pr(B_{i}(\bm{F}_{Nash})<B_{lim,i}))<0.25$ , where $B_{i}(\bm{F})$ is the long-term SSB of the $i$ th species under future fishing mortality $\bm{F}$ .

We used the ensemble model yield and SSB in 2050 to be the long-term yield, $f_{1,i}(\bm{F})$ , and long-term SSB, $B_{i}(\bm{F})$ , for $i=1,\ldots{},9$ , respectively. As running simulators and the ensemble model is computationally expensive, we used a Gaussian process emulator (Kennedy & O’Hagan 2001) to describe $f_{1,i}(\bm{F})$ and the 25th percentile of the long-term SSB of the $i$ th species under future fishing mortality $\bm{F}$ , $f_{2,i}(\bm{F})$ (i.e. $Pr(B_{i}(\bm{F})<f_{2,i}(\bm{F}))=0.25$ ), which we iteratively updated after rounds of simulations, allowing us to efficiently search for $\bm{F}_{Nash}$ values. A round consisted of running each simulator and the ensemble model for $\bm{F}$ values to find the yield and SSB. We ran four rounds to find the $\bm{F}$ values that satisfied Definition 3. The first round of 196 $\bm{F}$ were chosen using Sobol’ sequences, a space filling algorithm (Sobol’ 1967). For the subsequent rounds we proposed 100 $\bm{F}$ values that we belied may be $\bm{F}_{Nash}$ values according to the Gaussian process emulator.

We found $\bm{F}_{Nash}$ values by the following steps:

1.

Generate $\bm{F}^{(l)}$ , for $l=1,\ldots{}196$ , using Sobol’ sequences.
2.

Evaluate the simulators and the ensemble model at each of the new scenarios to find the yield and the SSB.
3.

Emulate the predictions of the long-term yield, $f_{1,1:9}(\bm{F})$ , and the 25th percentile of the long-term SSB from the ensemble model, $f_{2,1:9}(\bm{F})$ .
4.

Find 100 potential Nash equilibria, $\bm{F}^{(l)}$ (for $l=197,\ldots{},296$ in the second round, $l=297,\ldots{},396$ in the third round and $l=397,\ldots{},496$ in the fourth round).
5.

Repeat step 2 to 4 twice.
6.

Evaluate the simulators and the ensemble model the scenarios $\bm{F}^{(l)}$ , for $l=397,\ldots{},496$ , to find the yield and the SSB.

After the fourth round (step 6), the final $\bm{F}_{Nash}$ values were all of $\bm{F}^{(l)}$ scenarios that satisfy Definition 3, for $l=397,\ldots{},496$ . Due to the uncertainty in the ensemble model, we had several final $\bm{F}_{Nash}$ values. To try and distinguish between these we calculated the expected revenue for each of them.

The rest of the methods are as follows: the simulators, the ensemble model and the Gaussian process emulator are described in Sections 2.1, 2.2 and 2.3 respectively; an algorithm of how we find potential $\bm{F}_{Nash}$ values is described in Section 2.4 and we conclude by describing how we calculated the revenue of the long-term yield in Section 2.5.

2.1 Simulators

Four multispecies simulators were used: EcoPath with EcoSim (EwE Mackinson et al. 2018), LeMans (Thorpe et al. 2015), mizer (Blanchard et al. 2014) and FishSUMs (Speirs et al. 2016). All of them the simulators were able to describe the dynamics of all nine species with the exception of FishSUMs, which did not model sole.

To keep the interpretation of fishing mortality the same across simulators, we used the single-species assessments fishing mortality at age to drive the dynamics of the simulators (ICES 2018a; ICES 2018c). For the size-based simulators, LeMans, mizer and FishSUMs, we calculated the length at age using their respective von Bertalanffy parameters, however for EwE we used the $\bar{F}$ values from the assessments. In the future (2018-2050), the age selectivity was the same as those in 2017 and species that appear in the models but not in the study were fished at their 2017 levels.

2.2 Ensemble model

The predicted yields and SSB from the four multispecies simulators were combined using the ensemble model of Spence et al. 2018. The ensemble model is described in Table 2 and the simulator specific values are described in Table 3. We fit two ensemble models, one for the yields ( $j=1$ ) and one for the SSB ( $j=2$ ). For the yields, the simulators and the observations are in natural log tonnes, with the observations, $\hat{\bm{y}}^{(t)}_{1}$ , coming from ICES 2017. In the SSB model ( $j=2$ ), the simulators predicted SSB and the observations were in natural log tonnes, for their respective species and years, with the observations, $\hat{\bm{y}}^{(t)}_{2}$ , coming from stock-assessments (ICES 2018a; ICES 2018c).

Due to the high dimensionality and correlation of the uncertain parameter space, we fitted the ensemble model using No U-turn Hamiltonian Monte Carlo (Hoffman & Gelman 2011) in the package Stan (Stan Development Team 2020). We ran the algorithm for 2000 iterations discarding the first 1000 as burn in.

Table 2: A summary of the variables in the ensemble model. The ensemble model is run for 1984–2050. For values of

n_{k}

M_{k}

and

T_{k}

see Table 3.

Variable	Dimensions	$t$	Description	Relationship
$\bm{y}^{(t)}_{j}$	$9$	1984–2050	The truth	$\bm{y}^{(t)}_{j}\sim{}N(\bm{y}^{(t-1)}_{j},{\Lambda}_{y,j})$
$\bm{\hat{y}}^{(t)}_{j}$	$9$	1984–2017	Noisy observation of $\bm{y}^{(t)}_{j}$	$\hat{\bm{y}}_{j}^{(t)}\sim{}N({\bm{y}}^{(t)}_{j},\Sigma_{y,j})$
$\bm{\delta}_{j}$	$9$	NA	Long-term shared discrepancy
$\bm{\eta}^{(t)}_{j}$	$9$	1984–2050	Short-term shared discrepancy	$\bm{\eta}_{j}^{(t)}\sim{}N(R_{\eta,j}\bm{\eta}_{j}^{(t-1)},\Lambda_{\eta,j}$ )
$\bm{\mu}^{(t)}_{j}$	$9$	1984–2050	Simulator consensus	$\bm{\mu}^{(t)}_{j}=\bm{y}^{(t)}_{j}+\bm{\delta}_{j}+\bm{\eta}^{(t)}_{j}$
$\bm{\gamma}_{k,j}$	$9$	NA	Simulator $k$ ’s long-term individual discrepancy	$\bm{\gamma}_{k,j}\sim{}N(\bm{0},C_{\gamma,j})$
$\bm{z}_{k,j}^{(t)}$	$9$	1984–2050	Simulator $k$ ’s short-term individual discrepancy	$\bm{z}_{k,j}^{(t)}\sim{}N(R_{k,j}\bm{z}_{k,j}^{(t-1)},\Lambda_{k,j})$
$\bm{x}_{k,j}^{(t)}$	$9$	1984–2050	Simulator $k$ ’s best guess	$\bm{x}_{k,j}^{(t)}=\bm{\mu}_{j}^{(t)}+\bm{\gamma}_{k,j}+\bm{z}_{k,j}^{(t)}$
$\bm{\hat{x}}_{k,j}^{(t)}$	$n_{k}$	$T_{k}$	The expectation of simulator $k$ ’s output $\bm{x}_{k,j}^{(t)}$	$\bm{\hat{x}}_{k,j}^{(t)}\sim{}N(M_{k}\bm{x}_{k,j}^{(t)},\Sigma_{k,j})$

Table 3: A summary of the simulators, their outputs used in the case study, the simulator-specific values of

n_{k}

T_{k}

M_{k}

and

\Sigma_{k}


$k$	Simulator	Description	$n_{k}$	$T_{k}$	$M_{k}$	Reference for $\Sigma_{k}$
1	EcoPath with EcoSim (EwE)	An ecosystem model with 60 functional groups for the North Sea	$n_{1}=9$	$T_{1}=1991-2050$	$M_{1}=I_{9}$	Mackinson et al. 2018
2	LeMans	Abundance in length classes is modelled by species	$n_{2}=9$	$T_{2}=1986-2050$	$M_{2}=I_{9}.$	Thorpe et al. 2015
3	mizer	Total weight is modelled in weight classes by species	$n_{3}=9$	$T_{3}=1984-2050$	$M_{3}=I_{9}$	Spence et al. 2016
4	FishSUMs	Abundance in length classes is modelled by species	$n_{4}=8$	$T_{4}=1984-2050$ .	$M_{4}=\begin{pmatrix}1&0&0&0&0&0&0&0&0\\ 0&1&0&0&0&0&0&0&0\\ 0&0&1&0&0&0&0&0&0\\ 0&0&0&1&0&0&0&0&0\\ 0&0&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&1&0\\ 0&0&0&0&0&0&0&0&1\\ \end{pmatrix}$	Spence et al. 2018

2.3 Gaussian process emulator

The four simulators ran $m$ future fishing scenarios, $\bm{F}^{(l)}$ for $l=1,\ldots{},m$ , from 2018 until 2050. The ensemble model was evaluated at each of these future fishing scenarios to find the long-term yield, $f_{1,i}(\bm{F})$ , and long-term SSB. To find the Nash equilibrium we were required to evaluate $f_{1,i}(\bm{F})$ and the 25th percentile of the long-term SSB, $f_{2,i}(\bm{F})$ , of the $i$ th species at all $\bm{F}$ values. However, this was practically infeasible, as the simulators are relatively slow to run due to the computational complexity.

We used a Gaussian process emulator (Kennedy & O’Hagan 2001; Noè et al. 2019) to estimate $f_{1,i}(\bm{F})$ and $f_{2,i}(\bm{F})$ for all $\bm{F}$ values. If we let $\bm{f}_{j,i}=\left(f_{j,i}(\bm{F}^{(1)}),f_{j,i}(\bm{F}^{(2)}),\ldots,f_{j,i}(\bm{F}^{(m)})\right)^{\prime}$ , for $j=1$ and $2$ , then we say that

\bm{f}_{j,i}\sim{}GP(\bm{\eta}_{j,i},K_{j,i}),

where $\bm{\eta}_{j,i}$ was a generalised additive model (Wood 2017) and $K_{j,i}$ was the Matern covariance function (for more details see supplementary material), fitted using the DiceKriging package (Roustant et al. 2012) in R (R Core Team 2020).

2.4 Finding the Nash equilibrium

An algorithm to find the Nash equilibrium is to iteratively update $F_{i}$ by solving $F_{i}=F_{MSY,i}(\bm{F}_{-i})$ (Thorpe et al. 2017; Norrström et al. 2017). At each iteration, the long-yield term yield of the $i$ th species from the ensemble model was maximised by changing $F_{i}$ . In our case this meant maximising $f_{1,i}(F_{i},\bm{F}_{-i})$ , a stochastic function, such that $f_{2,i}(F_{i},\bm{F}_{-i})>B_{lim,i}$ . At each iteration, we sampled 500 potential $F_{i}$ values, using a Latin hypercube (McKay et al. 1979), and used the emulators to predict the long-term yield of the $i$ th species and the 25th percentile of the long-term SSB of all species. $F_{i}$ then became the potential new value with the largest estimated long-term yield for the $i$ th species such all of the species’ 25th percentile of their long-term SSB’s were above their respective $B_{lim}$ ’s. This is summarised in Algorithm 1.

Algorithm 1 A single iteration to find the Nash equilibrium.

LHC_{500}(0,2)

is 500 samples from a Latin hypercube of 1 dimension.

for

i

1:9

\bm{F}^{\prime}\sim{}LHC_{500}(0,2)

\bm{{\tilde{f}}}_{i,1}\sim{}GP(\bm{\eta}_{1,i}(\bm{F}^{\prime},\bm{F}_{-i}),K_{1,i})

\bm{{\tilde{f}}}_{1:9,2}\sim{}GP(\bm{\eta}_{2,i}(\bm{F}^{\prime},\bm{F}_{-i}),K_{2,i})

ll\leftarrow{}{\text{arg}\max}_{l}\left\{{{\tilde{f}}_{i,1,l}}:\tilde{f}_{i^{\prime},2,l}>B_{lim,i^{\prime}}\text{ for }i^{\prime}=1,\ldots{},9\right\}

F_{i}\leftarrow{}{F}_{ll}^{\prime}

end for

To initialise the algorithm, we sampled 10,000 $\bm{F}$ values, using a Latin hypercube design and used the emulators to estimate the long-term yield of each species and the 25th percentile of the long-term SSB. The initial $F_{i}$ value was the proposed fishing mortality for the $i$ th species that lead to the highest long-term yield such that the 0.25 percentile of the $i$ th species’ long-term SSB was above $B_{lim,i}$ .

The Nash equilibria was estimated with 100 samples from the posterior distribution of the ensemble model by repeating Algorithm 1 between 26 and 100 times, drawn at random. The resulting 100 samples were $\bm{F}_{Nash}$ values, which we ran the simulators and the ensemble model with.

2.5 Revenue of the long-term yield

For the $\bm{F}_{Nash}$ values, we calculated the expected revenue from the long-term yields for each Nash equilibrium found using Algorithm 1. To derive the revenue, we predicted the prices for each year until 2050 using a uni-variate Vector Auto-Regressive estimation model (VAR) and landings values per tonne per species (deflated) of the UK fleet in England from 1970-2018. The value of the landings per tonne are shown in Table 1.

3 Results

3.1 Simulator runs

Each simulator was run for 496 different fishing scenarios, $\bm{F}^{(l)}$ for $l=1,\ldots{}496$ . Figure 1 shows the historical yields and each of the simulators predicted yields for the period 1985 to 2017. Most of the simulators were able to qualitatively recreate the trends of the observed yields for most of the species, however no single simulator appears to be overall better than the others.

Refer to caption — Figure 1: Historical yield from observations (ICES 2017) and the simulators.

3.2 Ensemble outputs

We fitted the ensemble model and used it to describe, with uncertainty, what the yield and SSB would be under the future fishing scenarios. The median long-term yield for all of the scenarios is shown in Figure 2. Although the long-term yield and SSB were sensitive to the fishing mortality of that species, it was also sensitive to the fishing mortality of other species. Figure 3 shows the 5th and 25th percentile of the long-term SSB for cod and whiting for varying fishing mortality of cod respectively, with the solid lines being their $B_{lim}$ values. The long-term SSB’s of whiting and cod appear to be negatively correlated.

3.3 Gaussian process emulator

We fitted both the long-term yield and the 25th percentile of the long-term SSB for 100 iterations of the ensemble model. Table 4 shows the number of scenarios in each round and the number of species with acceptable risk, a probability that the long-runs SSB is above $B_{lim}$ of 0.75 or more. In later rounds the number of species with acceptable risk increases.

3.4 Nash equilibria

Out of the 100 potential $\bm{F}_{Nash}$ values found in the fourth round, 39 of them satisfied Definition 3. For these 39 $\bm{F}_{Nash}$ values, we calculated the revenue of the long-term yield. Figure 4 shows the marginal distributions of the accepted $\bm{F}_{Nash}$ values, with the solid line showing the $\bm{F}_{Nash}$ value that led to the highest revenue. The revenues of all the $\bm{F}_{Nash}$ values were between £1.7 billion and £2.2 billion, larger than the revenue in 2017, £1.3 billion. See Table S1 in the supplementary material for the 39 $\bm{F}_{Nash}values$ .

Table 4: The number of scenarios that have acceptable risk to the species long-term SSB. Acceptable risk to a species is that the long-term SSB is above

B_{lim}

with a probability of 0.75 or more.

# species	Rnd 1	Rnd 2	Rnd 3	Rnd 4
0	$0$	$0$	$0$	$0$
1	$14$	$0$	$0$	$0$
2	$45$	$0$	$0$	$0$
3	$42$	$0$	$0$	$0$
4	$37$	$1$	$0$	$0$
5	$31$	$12$	$1$	$0$
6	$19$	$33$	$2$	$0$
7	$8$	$43$	$28$	$13$
8	$0$	$11$	$46$	$48$
9	$0$	$0$	$23$	$39$

4 Discussion

In this paper, we showed that using a SS-MSY is only possible under very strict assumptions. We demonstrated how to calculate reference points using a specific definition of MMSY, the Nash equilibrium (Farcas & Rossberg 2016; Norrström et al. 2017; Thorpe et al. 2017), with a caveat for the risk of species collapse for nine species in the North Sea. We did this by combining multiple simulators using an ensemble model, removing the arbitrary choices of which simulator to use to calculate the reference points. We found that the Nash equilibrium led to higher fishing mortality rates than SS-MSY, leading to an increase in the long-term yields and revenue. To our knowledge, ensemble modelling has never been used to calculate multispecies reference points before.

We found that the $\bm{F}_{Nash}$ values were generally higher than SS-MSY. Fishing predators at higher levels can relieve stress on prey, leading to an increase in prey, which can be exploited by the fishery (Andersen et al. 2015). These interactions are not accounted for when calculating SS-MSY, therefore adopting a MMSY means it is possible to have a higher yield (Beddington & Cooke 1982; Norrström et al. 2017). Furthermore the revenue generated from MMSY was greater than current levels, allowing an economic gain for fishers and their families, something that is not always the case for SS-MSY (Giron-Nava et al. 2019).

A common criteria when defining SS-MSY is that the SSB is larger than $B_{lim}$ with a probability greater than 0.95 (ICES 2018b). We demonstrated that the SSB of a single species not only depends on its own fishing mortality, but also the fishing mortality of the other species. For example, only a small range of cod fishing mortality would lead to both whiting and cod’s long-term SSB being above $B_{lim}$ with a probability greater than 0.75. However, no combination of $\bm{F}$ values would result in both whiting and cod’s long-term SSB being above $B_{lim}$ with a probability greater than 0.95 (Figure 3). This effect between whiting and cod was also found by EwE (Mackinson et al. 2009) and the stochastic multispecies model (Lewy & Vinther 2004; Kempf et al. 2010). As it is impossible to find $\bm{F}$ values that satisfy the 0.95 probability criteria for all nine species in this study, we reduced our criteria to 0.75. In general uncertainty is subjective, specific to the decision maker, the study, the information and the simulators (Gelman et al. 2013). In this study our certainty is limited to the simulators used, and could be reduced if they were improved, however, we were able quantify this uncertainty in a robust and interpretable manner (Harwood & Stokes 2003). Currently when calculating SS-MSY, large amounts of uncertainty are ignored, e.g. species interactions, and thus estimations of probability are not robust, which makes the 0.95 caveat rather arbitrary.

Generally, fisheries managers select a single simulator for a species, from a set of competing simulators to calculate reference points (ICES 2018b, e.g), however, not including species interactions can lead to inconsistencies in the reference points. Should a manager chose a simulator for defining an MMSY, they would have to decide which simulator based on an arbitrary choice, and the values of the reference points are sensitive to the simulator (see supplementary material Figures S1-S4) (Gaichas 2008). Furthermore, choosing a simulator, without accounting for its model discrepancy (Kennedy & O’Hagan 2001), can lead to biased advice. For example, if we chose LeMans, sandeel yields would be consistently under-estimated, however correcting for the discrepancy would lead to estimations that were closer to the true yields (Figure 1). In general, no simulator is uniformly better than the others (Chandler 2013). In our example, mizer captures the dynamics of the saithe yields, however it does not capture the cod yields as well (Figure 1). We combined four different simulators, accounting for their discrepancies and uncertainties, to define the reference points, suggesting their values are no longer sensitive to the simulator selection (Spence et al. 2018).

Due to the robust quantification of uncertainty in the ensemble model, we found 39 different $\bm{F}_{Nash}$ values. In practice, a manager would have to decide which of the $\bm{F}_{Nash}$ values is the ‘best’, which is dependent on their needs and priorities. For example, they may want to maximise the total revenue or to minimise the risk to a specific species. In general, the manager’s utility can be computed for the different MMSYs and then they can decide which of them is the ‘best’. We calculated the revenue for each $\bm{F}_{Nash}$ value, and, if we were to give advice, we would select the $\bm{F}_{Nash}$ value that would lead to the highest revenue, as shown in Figure 4.

The results should be interpreted in the light of the limitations of the four simulators used in this study, which were the only ones available. Using as many simulators as possible would improve the robustness of the results, however it would be more beneficial to use better or improved simulators. By robustly quantifying the uncertainty, the ensemble model uses all of the information from the simulators. If there was no, or very little, information in all the simulators then the ensemble model would give very uncertain predictions. Currently the ensemble model of Spence et al. 2018 assumes that the discrepancies of the simulators are the same in the future as they are in the past, for example a simulator that was uncertain at predicting the past would also be uncertain when predicting the future. More work is required to find the predictive power of these simulators, so we can include this information in the ensemble model.

When calculating the Nash equilibrium, we would like to use a sequential algorithm, such as in Norrström et al. 2017, which would require running the simulators and the ensemble model many times. Currently this is not feasible due to computational and time constraints caused by the simulators. To limit the number of simulator runs required, we used a Gaussian process emulator to predict, with uncertainty, what the ensemble model would say for all future scenarios. Gaussian process emulators have been used in other fields when simulators are expensive to run (Vernon et al. 2014; Kennedy et al. 2006, e.g). This allowed us to limit simulator runs to fishing scenarios that may be close to the Nash equilibrium and result in acceptable risk (Table 4), or where the emulator was unsure of the outcome. Although replacing the ensemble model with a Gaussian process leads to uncertainty in the final $\bm{F}_{Nash}$ values, this would not matter in practice, as the uncertainty Gaussian process will be small. Using the ensemble model and the Gaussian process emulator allows for the calculation of MMSY in a robust and timely manner.

The Nash equilibrium was calculated for nine species in the North Sea, with caveats for the risk of stock collapse, although the methods described would be applicable for any definition of MMSY, or even SS-MSY, at any location. Alternative objectives could be ecosystem based yield (Steele et al. 2011) or to aim to either maximise profits in the fisheries (e.g. using an MEY approach (Dichmont et al. 2010; Pascoe et al. 2018; Guillen et al. 2013)) or focus on the efficiency of the fishing practice. The latter, for example, could aim to define the reference points based on marginal value of yields, apply pareto-efficiency criteria or include joint-technology in production (i.e. mixed fisheries considerations). In this paper, the Nash equilibrium was chosen as it is a way of combining the SS-MSY with the MMSY as aligned concepts of MSY and EBFM (Norrström et al. 2017).

5 Conclusion

In this study we calculated MMSY in the North Sea using an ensemble model, demonstrating that it can lead to sustainable yields whilst ensuring ecosystem health is not diminished. This approach can be adopted by fisheries scientists and mangers worldwide, taking account of structural uncertainties and removing arbitrary modelling decisions, leading to more robust, and therefore better science and management. The reference points can be applied to problems in other fields such as climate science, epidemiology or systems biology. Using the methods described in this paper, we were able to provide a practical tool to optimise any objective function for use by scientists and managers alike.

Acknowledgements

The work was funded by the Department for Environment, Food and Rural Affairs (Defra). We would like to thank Robert Thorpe, Michaela Schratzberger and Paul Dolder for comments on earlier versions of the manuscript.

Authors contribution

MAS, HJB and KA conceived the ideas and designed the methodology; NDW, AM and MAS extracted data for the study; MAS, KA and HJB ran simulators; MAS and KA led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Data availability statement

Data sharing is not applicable to this article as no new data were created; rather, data were acquired from existing published sources (all sources are cited in the text), or are described, figured and tabulated within the manuscript or supplementary information of this article.

References

Andersen (2019) Andersen, K.H. (2019) Fish Ecology, Evolution, and Exploitation A New Theoretical Synthesis. Princeton University Press.
Andersen et al. (2015) Andersen, K.H., Brander, K. & Ravn-Jonsen, L. (2015) Trade-offs between objectives for ecosystem management of fisheries. Ecological Applications, 25, 1390–1396.
Beddington & Cooke (1982) Beddington, J. & Cooke, J. (1982) Harvesting from a prey-predator complex. Ecological Modelling, 14, 155 – 177. Ecology, Renewable Resources and Optimal Control.
Blanchard et al. (2014) Blanchard, J.L., Andersen, K.H., Scott, F., Hintzen, N.T., Piet, G. & Jennings, S. (2014) Evaluating targets and trade-offs among fisheries and conservation objectives using a multispecies size spectrum model. Journal of Applied Ecology, 51, 612–622.
Chandler (2013) Chandler, R.E. (2013) Exploiting strength, discounting weakness: combining information from multiple climate simulators. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371, 20120388.
Collie et al. (2016) Collie, J.S., Botsford, L.W., Hastings, A., Kaplan, I.C., Largier, J.L., Livingston, P.A., Plagányi, E., Rose, K.A., Wells, B.K. & Werner, F.E. (2016) Ecosystem models for fisheries management: finding the sweet spot. Fish and Fisheries, 17, 101–125.
Dichmont et al. (2010) Dichmont, C.M., Pascoe, S., Kompas, T., Punt, A.E. & Deng, R. (2010) On implementing maximum economic yield in commercial fisheries. Proceedings of the National Academy of Sciences, 107, 16–21.
Essington & Punt (2011) Essington, T. & Punt, A. (2011) Implementing ecosystem-based fisheries management: Advances, challenges and emerging tools. Fish and Fisheries, 12.
Essington & Plagányi (2013) Essington, T.E. & Plagányi, E.E. (2013) Pitfalls and guidelines for “recycling” models for ecosystem-based fisheries management: evaluating model suitability for forage fish fisheries. ICES Journal of Marine Science, 71, 118–127.
FAO (2009) FAO (2009) How to Feed the World in 2050 - Food and Agriculture organization. Http://www.fao.org/docrep/pdf/012/ak542e/ak542e00.pdf.
Farcas & Rossberg (2016) Farcas, A. & Rossberg, A.G. (2016) Maximum sustainable yield from interacting fish stocks in an uncertain world: two policy choices and underlying trade-offs. ICES Journal of Marine Science, 73, 2499–2508.
Fulton et al. (2003) Fulton, E., Smith, A. & Johnson, C. (2003) Effect of complexity of marine ecosystem models. Marine Ecology Progress Series, 253, 1–16.
Gaichas (2008) Gaichas, S.K. (2008) A context for ecosystem-based fishery management: Developing concepts of ecosystems and sustainability. Marine Policy, 32, 393 – 401.
Gelman et al. (2013) Gelman, A., Carlin, J., Stern, H., Dunson, D., Vehtari, A. & Rubin, D. (2013) Bayesian Data Analysis. Chapman and Hall/CRC, third edition edition.
Giron-Nava et al. (2019) Giron-Nava, A., Johnson, A.F., Cisneros-Montemayor, A.M. & Aburto-Oropeza, O. (2019) Managing at maximum sustainable yield does not ensure economic well-being for artisanal fishers. Fish and Fisheries, 20, 214–223.
Guillen et al. (2013) Guillen, J., Macher, C., Merzéréaud, M., Bertignac, M., Fifas, S. & Guyader, O. (2013) Estimating MSY and MEY in multi-species and multi-fleet fisheries, consequences and limits: an application to the bay of biscay mixed fishery. Marine Policy, 40, 64 – 74.
Hart & Fay (2020) Hart, A.R. & Fay, G. (2020) Applying tree analysis to assess combinations of ecosystem-based fisheries management actions in management strategy evaluation. Fisheries Research, 225, 105466.
Harwood & Stokes (2003) Harwood, J. & Stokes, K. (2003) Coping with uncertainty in ecological advice: lessons from fisheries. Trends in Ecology & Evolution, 18, 617 – 622.
Hilborn (2007) Hilborn, R. (2007) Defining success in fisheries and conflicts in objectives. Marine Policy, 31, 153 – 158.
Hilborn & Walters (1992) Hilborn, R. & Walters, C.J. (1992) Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty. Springer Science.
Hoffman & Gelman (2011) Hoffman, M. & Gelman, A. (2011) The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo. Journal of Machine Learning Research, 15.
Hollowed et al. (2000) Hollowed, A.B., Bax, N., Beamish, R., Collie, J., Fogarty, M., Livingston, P., Pope, J. & Rice, J.C. (2000) Are multispecies models an improvement on single-species models for measuring fishing impacts on marine ecosystems? ICES Journal of Marine Science, 57, 707–719.
Hyder et al. (2015) Hyder, K., Rossberg, A.G., Allen, J.I., Austen, M.C., Barciela, R.M., Bannister, H.J., Blackwell, P.G., Blanchard, J.L., Burrows, M.T., Defriez, E., Dorrington, T., Edwards, K.P., Garcia-Carreras, B., Heath, M.R., Hembury, D.J., Heymans, J.J., Holt, J., Houle, J.E., Jennings, S., Mackinson, S., Malcolm, S.J., McPike, R., Mee, L., Mills, D.K., Montgomery, C., Pearson, D., Pinnegar, J.K., Pollicino, M., Popova, E.E., Rae, L., Rogers, S.I., Speirs, D., Spence, M.A., Thorpe, R., Turner, R.K., van der Molen, J., Yool, A. & Paterson, D.M. (2015) Making modelling count - increasing the contribution of shelf-seas community and ecosystem models to policy development and management. Marine Policy, 61, 291–302.
ICES (2017) ICES (2017) Official Nominal Catches. http://ices.dk/marine-data/dataset-collections/Pages/Fish-catch-and-stock-assessment.aspx.
ICES (2018a) ICES (2018a) Herring Assessment Working Group for the Area South of 62 N (HAWG). Technical report, ICES Scientific Reports. ACOM:07. 960 pp, ICES, Copenhagen.
ICES (2018b) ICES (2018b) ICES Advice basis. Technical report, International Council for Exploration of the Seas.
ICES (2018c) ICES (2018c) Report of the Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak. Technical report, ICES Scientific Reports. ACOM:22. pp, ICES, Copenhagen.
Kempf et al. (2010) Kempf, A., Dingsør, G.E., Huse, G., Vinther, M., Floeter, J. & Temming, A. (2010) The importance of predator-prey overlap: predicting North Sea cod recovery with a multispecies assessment model. ICES Journal of Marine Science, 67, 1989–1997.
Kennedy et al. (2006) Kennedy, M.C., Anderson, C.W., Conti, S. & O’Hagan, A. (2006) Case studies in Gaussian process modelling of computer codes. Reliability Engineering & System Safety, 91, 1301 – 1309. The Fourth International Conference on Sensitivity Analysis of Model Output (SAMO 2004).
Kennedy & O’Hagan (2001) Kennedy, M.C. & O’Hagan, A. (2001) Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63, 425–464.
Larkin (1977) Larkin, P. (1977) An epitaph for the concept of maximum sustained yield. Transactions of the American Fisheries Society, 106, 1–11.
Lewy & Vinther (2004) Lewy, P. & Vinther, M. (2004) A stochastic age-length-structured multispecies model applied to north sea stocks. Technical report, ICES.
Link et al. (2011) Link, J.S., Bundy, A., Overholtz, W.J., Shackell, N., Manderson, J., Duplisea, D., Hare, J., Koen-Alonso, M. & Friedland, K.D. (2011) Ecosystem-based fisheries management in the Northwest Atlantic. Fish and Fisheries, 12, 152–170.
Mackinson et al. (2009) Mackinson, S., Deas, B., Beveridge, D. & Casey, J. (2009) Mixed-fishery or ecosystem conundrum? multispecies considerations inform thinking on long-term management of north sea demersal stocks. Canadian Journal of Fisheries and Aquatic Sciences, 66, 1107–1129.
Mackinson et al. (2018) Mackinson, S., Platts, M., Garcia, C. & Lynam, C. (2018) Evaluating the fishery and ecological consequences of the proposed North Sea multi-annual plan. PLOS ONE, 13, 1–23.
McKay et al. (1979) McKay, M.D., Beckman, R.J. & Conover, W.J. (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21, 239–245.
Mesnil (2012) Mesnil, B. (2012) The hesitant emergence of maximum sustainable yield (MSY) in fisheries policies in europe. Marine Policy, 36, 473 – 480.
Nash (1951) Nash, J. (1951) Non-cooperative games. Annals of Mathematics, 54, 286–295.
Nielsen et al. (2018) Nielsen, J.R., Thunberg, E., Holland, D.S., Schmidt, J.O., Fulton, E.A., Bastardie, F., Punt, A.E., Allen, I., Bartelings, H., Bertignac, M., Bethke, E., Bossier, S., Buckworth, R., Carpenter, G., Christensen, A., Christensen, V., Da-Rocha, J.M., Deng, R., Dichmont, C., Doering, R., Esteban, A., Fernandes, J.A., Frost, H., Garcia, D., Gasche, L., Gascuel, D., Gourguet, S., Groeneveld, R.A., Guillén, J., Guyader, O., Hamon, K.G., Hoff, A., Horbowy, J., Hutton, T., Lehuta, S., Little, L.R., Lleonart, J., Macher, C., Mackinson, S., Mahevas, S., Marchal, P., Mato-Amboage, R., Mapstone, B., Maynou, F., Merzéréaud, M., Palacz, A., Pascoe, S., Paulrud, A., Plaganyi, E., Prellezo, R., van Putten, E.I., Quaas, M., Ravn-Jonsen, L., Sanchez, S., Simons, S., Thébaud, O., Tomczak, M.T., Ulrich, C., van Dijk, D., Vermard, Y., Voss, R. & Waldo, S. (2018) Integrated ecological-economic fisheries models-evaluation, review and challenges for implementation. Fish and Fisheries, 19, 1–29.
Noè et al. (2019) Noè, U., Lazarus, A., Gao, H., Davies, V., Macdonald, B., Mangion, K., Berry, C., Luo, X. & Husmeier, D. (2019) Gaussian process emulation to accelerate parameter estimation in a mechanical model of the left ventricle: a critical step towards clinical end-user relevance. Journal of The Royal Society Interface, 16, 20190114.
Norrström et al. (2017) Norrström, N., Casini, M. & Holmgren, N. (2017) Nash equilibrium can resolve conflicting maximum sustainable yields in multi-species fisheries management. ICES Journal of Marine Science, 74, 78–90.
Ok (2017) Ok, E.A. (2017) Real Analysis with Economic Applications. Princeton University Press.
Pascoe et al. (2018) Pascoe, S., Hutton, T. & Hoshino, E. (2018) Offsetting externalities in estimating MEY in multispecies fisheries. Ecological Economics, 146, 304 – 311.
Pauly & Froese (2014) Pauly, D. & Froese, R. (2014) Fisheries Management. American Cancer Society.
Pikitch et al. (2004) Pikitch, E.K., Santora, C., Babcock, E.A., Bakun, A., Bonfil, R., Conover, D.O., Dayton, P., Doukakis, P., Fluharty, D., Heneman, B., Houde, E.D., Link, J., Livingston, P.A., Mangel, M., McAllister, M.K., Pope, J. & Sainsbury, K.J. (2004) Ecosystem-Based Fishery Management. Science, 305, 346–347.
R Core Team (2020) R Core Team (2020) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Roustant et al. (2012) Roustant, O., Ginsbourger, D. & Deville, Y. (2012) DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization. Journal of Statistical Software, 51, 1–55.
Säterberg et al. (2019) Säterberg, T., Casini, M. & Gardmark, A. (2019) Ecologically Sustainable Exploitation Rates-A multispecies approach for fisheries management. Fish and Fisheries, 20, 952–961.
Sissenwine & Shepherd (1987) Sissenwine, M.P. & Shepherd, J.G. (1987) An alternative perspective on recruitment overfishing and biological reference points. Canadian Journal of Fisheries and Aquatic Sciences, 44, 913–918.
Sobol’ (1967) Sobol’, I. (1967) On the distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics, 7, 86 – 112.
Speirs et al. (2016) Speirs, D., Greenstreet, S. & Heath, M. (2016) Modelling the effects of fishing on the North Sea fish community size composition. Ecological Modelling, 321, 35–45.
Spence et al. (2016) Spence, M.A., Blackwell, P.G. & Blanchard, J.L. (2016) Parameter uncertainty of a dynamic multispecies size spectrum model. Canadian Journal of Fisheries and Aquatic Sciences, 73, 589–597.
Spence et al. (2018) Spence, M.A., Blanchard, J.L., Rossberg, A.G., Heath, M.R., Heymans, J.J., Mackinson, S., Serpetti, N., Speirs, D.C., Thorpe, R.B. & Blackwell, P.G. (2018) A general framework for combining ecosystem models. Fish and Fisheries, 19, 1031–1042.
Stan Development Team (2020) Stan Development Team (2020) RStan: the R interface to Stan. R package version 2.19.3.
Steele et al. (2011) Steele, J.H., Gifford, D.J. & Collie, J.S. (2011) Comparing species and ecosystem-based estimates of fisheries yields. Fisheries Research, 111, 139 – 144.
Szuwalski et al. (2016) Szuwalski, C., Burgess, M., Costello, C. & Gaines, S. (2016) High fishery catches through trophic cascades in China. Proceedings of the National Academy of Sciences, 114, 201612722.
Thorpe (2019) Thorpe, R.B. (2019) What is multispecies msy? a worked example from the north sea. Journal of Fish Biology, 94, 1011–1018.
Thorpe et al. (2017) Thorpe, R.B., Jennings, S. & Dolder, P.J. (2017) Risks and benefits of catching pretty good yield in multispecies mixed fisheries. ICES Journal of Marine Science, 74, 2097–2106.
Thorpe et al. (2015) Thorpe, R.B., Le Quesne, W.J.F., Luxford, F., Collie, J.S. & Jennings, S. (2015) Evaluation and management implications of uncertainty in a multispecies size-structured model of population and community responses to fishing. Methods in Ecology and Evolution, 6, 49–58.
Tsikliras & Froese (2018) Tsikliras, A. & Froese, R. (2018) Maximum Sustainable Yield, pp. 108–115. Elsevier.
Vernon et al. (2014) Vernon, I., Goldstein, M. & Bower, R. (2014) Galaxy formation : Bayesian history matching for the observable universe. Statistical science, 29, 81–90.
Wood (2017) Wood, S.N. (2017) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC, second edition edition.
Worm et al. (2009) Worm, B., Hilborn, R., Baum, J.K., Branch, T.A., Collie, J.S., Costello, C., Fogarty, M.J., Fulton, E.A., Hutchings, J.A., Jennings, S., Jensen, O.P., Lotze, H.K., Mace, P.M., McClanahan, T.R., Minto, C., Palumbi, S.R., Parma, A.M., Ricard, D., Rosenberg, A.A., Watson, R. & Zeller, D. (2009) Rebuilding global fisheries. Science, 325, 578–585.

S1 MSY

S1.1 Function definition

Let $y_{i}(F_{1},F_{2},\ldots{}F_{n})$ be a continuous function such that

f_{1,i}:D\xrightarrow{}\mathbb{R}_{\geq 0}

with

D=\left\{(x_{1}\times{}x_{2}\times\ldots\times{}x_{n})\in\mathbb{R}^{n}_{\geq 0}\right\}.

S1.2 Proof of Proposition 1

Proof.

As $f_{1,i}(F_{i},\bm{F}_{-i})$ is a continuous function, then $F_{MSY,i}(\bm{F}_{-i})={\text{arg}\sup}_{F_{i}}\left(f_{1,i}(F_{i},\bm{F}_{-i})\right)$ is also a continuous function due to the maximum theorem (Ok 2017). Suppose

\frac{\partial{}F_{MSY,i}(\bm{F}_{-i})}{\partial{}F_{j}}=0

then

$\displaystyle F_{MSY,i}$	$\displaystyle=$	$\displaystyle F_{MSY,i}(\bm{F}_{-i,j},F_{j})$
	$\displaystyle=$	$\displaystyle\lim_{\delta\to 0}F_{MSY,i}(\bm{F}_{-i,j},F_{j}+\delta)$
	$\displaystyle=$	$\displaystyle F_{MSY,i}.$

Now suppose

\frac{\partial{}F_{MSY,i}(\bm{F}_{-i})}{\partial{}F_{j}}>0

then

$\displaystyle F_{MSY,i}$	$\displaystyle=$	$\displaystyle F_{MSY,i}(\bm{F}_{-i,j},F_{j})$
	$\displaystyle<$	$\displaystyle\lim_{\delta\to 0}F_{MSY,i}(\bm{F}_{-i,j},F_{j}+\delta)$
	$\displaystyle=$	$\displaystyle F_{MSY,i}^{\prime}$

hence $F_{MSY,i}\neq{}F_{MSY,i}^{\prime}$ . Alternatively suppose

\frac{\partial{}F_{MSY,i}(\bm{F}_{-i})}{\partial{}F_{j}}<0

then

$\displaystyle F_{MSY,i}$	$\displaystyle=$	$\displaystyle F_{MSY,i}(\bm{F}_{-i,j},F_{j})$
	$\displaystyle>$	$\displaystyle\lim_{\delta\to 0}F_{MSY,i}(\bm{F}_{-i,j},F_{j}+\delta)$
	$\displaystyle=$	$\displaystyle F_{MSY,i}^{\prime}$

hence $F_{MSY,i}\neq{}F_{MSY,i}^{\prime}$ . Hence

\frac{\partial{}F_{MSY,i}(\bm{F}_{-i})}{\partial{}F_{j}}=0

$\forall{}j\neq{i}$ if Definition 2 is to exist. ∎

S2 Gaussian process emulator

A stochastic process $f_{j,i}(\bm{F})$ is said to be a Gaussian process if the random vector, $\bm{f}_{j,i}=\left(f_{j,i}(\bm{F}^{(1)}),f_{j,i}(\bm{F}^{(2)}),\ldots,f_{j,i}(\bm{F}^{(n)})\right)^{\prime}$ , for $j=1$ and $2$ and $i=1,\ldots{}9$ , has the distribution

\bm{f}_{j,i}\sim{}N(\bm{\eta}_{j,i},K_{j,i}).

Similarly to a multivariate Gaussian, completely specified by a mean vector and a covariance matrix, the Gaussian Process is parametried by a mean and a covariance function with

\eta_{j,i}(\bm{F})=E(f_{j,i}(\bm{F}))

and

k_{j,i}(\bm{F}^{(l)},\bm{F}^{(l^{\prime})})=Cov\left(f_{j,i}(\bm{F}^{(l)}),f_{j,i}(\bm{F}^{(l^{\prime})})\right)

respectively, returning the mean of a random variable and the covariance between two random variables, as function of the inputs only (Noè et al. 2019). In this work we consider $\eta_{j,i}$ to be a generalised additive model (Wood 2017), see Section S2.1. We used the covariance $k_{j,i}(\bm{F}^{(l)},\bm{F}^{(l^{\prime})})=C_{j,i,1}(F_{1}^{(l)},F^{(l^{\prime})}_{1})\otimes{}\ldots\otimes{}C_{j,i,9}(F_{9}^{(l)},F_{9}^{(l^{\prime})})$ with a Matèrn covariance function,

	$\displaystyle C_{j,i,d}(F_{d}^{(l)},F_{d}^{(l^{\prime})})=$		$\displaystyle\sigma^{2}\left(1+\sqrt{5}\frac{\|F_{d}^{(l^{\prime})}-F_{d}^{(l)}\|}{\rho_{j,i,d}}+\frac{5}{3}\left(\frac{\|F_{d}^{(l^{\prime})}-F_{d}^{(l)}\|}{\rho_{j,i,d}}\right)^{2}\right)$
			$\displaystyle\times\exp\left(\frac{-\sqrt{5}\|F_{d}^{(l^{\prime})}-F_{d}^{(l)}\|}{\rho_{j,i,d}}\right),$

for $d=1\ldots{}9$ .

Denote the observed data $\mathcal{D}=\left\{(\bm{F}^{(1)},y_{j,i}^{(1)}),\ldots{},(\bm{F}^{(m)},y_{j,i}^{(m)})\right\}$ to be training data, with inputs $\bm{F}^{(l)}$ and outputs $y_{j,i}^{(l)}$ for $l=1,\ldots{},m$ . The outputs are denoted $\bm{y}_{j,i}=(y_{j,i}^{(1)},\ldots,y_{j,i}^{(m)})^{\prime}$ . Conditioning the Gaussian process on the observed data

f_{j,i}(\bm{F})\sim{}GP(\tilde{f}_{j,i}(\bm{F}),s(\bm{F},\bm{F}^{\prime}))

with

\tilde{f}_{j,i}(\bm{F})=\eta_{j,i}(\bm{F})+\bm{k}(\bm{F})^{\prime}(K+\sigma^{2}I)^{-1}(\bm{y}_{j,i}-\bm{\eta}_{j,i})

and

s(\bm{F},\bm{F}^{\prime})=k(\bm{F},\bm{F}^{\prime})-\bm{k}(\bm{F})^{\prime}(K+\sigma^{2}I)^{-1}\bm{k}(\bm{F}^{\prime}),

where $\bm{k}(\bm{F})=(k(\bm{F},\bm{F}^{(1)}),\ldots{}k(\bm{F},\bm{F}^{(m)}))^{\prime}$ , $K=\left[k(\bm{F}^{(l)},\bm{F}^{(l^{\prime})})\right]^{m}_{l,l^{\prime}=1}$ is the training covariance, $\bm{\eta}_{j,i}=(\eta_{j,i}(\bm{F}^{(1)}),\ldots,\eta_{j,i}(\bm{F}^{(m)}))^{\prime}$ and $I$ is the identity matrix of dimensions $m$ (Noè et al. 2019).

S2.1 Generalised additive models

The mean function from the Gaussian process emulator was a cubic spline such that

s(x)=\sum_{h=1}^{H}1_{x\geq\lambda_{h}}\beta_{h}(x-\lambda_{k})^{3},

where $H$ is the number of ‘knots’ and $\lambda_{k}$ is the location of the $k$ th ‘knot’.

Sandeel

The yield for sandeel was

\eta_{1,1}(\bm{F})=\beta_{1,1}+s(F_{1})+s(F_{3})+s(F_{5}),

and the SSB was

\eta_{2,1}(\bm{F})=\beta_{2,1}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{5})+s(F_{8})+s(F_{9}).

Norway pout

The yield for Norway pout was

\eta_{1,2}(\bm{F})=\beta_{1,2}+s(F_{2})+s(F_{3})+s(F_{5}),

and the SSB was

\eta_{2,2}(\bm{F})=\beta_{2,2}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{8})+s(F_{9}).

Herring

The yield for herring was

\eta_{1,3}(\bm{F})=\beta_{1,3}+s(F_{3})+s(F_{5})+s(F_{8})+s(F_{9}),

and the SSB was

\eta_{2,3}(\bm{F})=\beta_{1,3}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{5})+s(F_{6})+s(F_{8})+s(F_{9}).

Whiting

The yield for whiting was

\eta_{1,4}(\bm{F})=\beta_{1,4}+s(F_{3})+s(F_{4}),

and the SSB was

\eta_{2,4}(\bm{F})=\beta_{2,4}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{5})+s(F_{7})+s(F_{8})+s(F_{9}).

Sole

The yield for sole was

\eta_{1,5}(\bm{F})=\beta_{1,5}+s(F_{3})+s(F_{4})+s(F_{5}),

and the SSB was

\eta_{2,5}(\bm{F})=\beta_{2,5}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{5})+s(F_{6})+s(F_{7})+s(F_{8})+s(F_{9}).

Plaice

The yield for plaice was

\eta_{1,6}(\bm{F})=\beta_{1,6}+s(F_{1})+s(F_{3})+s(F_{6})+s(F_{7}),

and the SSB was

\eta_{2,6}(\bm{F})=\beta_{2,6}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4}).

Haddock

The yield for haddock was

\eta_{1,7}(\bm{F})=\beta_{1,7}+s(F_{3})+s(F_{4})+s(F_{7})+s(F_{8}),

and the SSB was

\eta_{2,7}(\bm{F})=\beta_{2,7}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{5})+s(F_{6})+s(F_{7})+s(F_{8}).

Cod

The yield for cod was

\eta_{1,8}(\bm{F})=\beta_{1,8}+s(F_{3})+s(F_{5})+s(F_{8}),

and the SSB was

\eta_{2,8}(\bm{F})=\beta_{2,8}+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{5})+s(F_{8}).

Saithe

The yield for saithe was

\eta_{1,9}(\bm{F})=\beta_{1,9}+s(F_{3})+s(F_{5})+s(F_{8})+s(F_{9}),

and the SSB was

\eta_{2,9}(\bm{F})=\beta_{2,9}+s(F_{1})+s(F_{2})+s(F_{3})+s(F_{4})+s(F_{8})+s(F_{9}).

S3 Results

S3.1 Simulator runs

Figures S1-S4 show the long-term yields from the simulators.

S3.2 Spawning stock biomass

Figure S5 shows the 25th percentile of the long-term SSB. The solid line is the $B_{lim}$ for each species.

S3.3 Reference points

Table S1 shows the 39 Nash equilibria and their expected long-term revenue that satisfy have acceptable risk to the species long-term SSB.

Table S1: The 39 values of

\bm{F}_{Nash}

and their expected long-term revenue that we found in this study.

Sandeel	N.pout	Herring	Whiting	Sole	Plaice	Haddock	Cod	Saithe	Revenue (£billions)
$1.05$	$1.49$	$0.46$	$0.82$	$0.31$	$0.48$	$0.94$	$0.62$	$1.10$	$2.16$
$1.10$	$1.47$	$0.44$	$0.87$	$0.31$	$0.48$	$0.76$	$0.63$	$1.16$	$2.15$
$1.11$	$1.44$	$0.39$	$0.85$	$0.37$	$0.51$	$0.82$	$0.63$	$1.13$	$2.12$
$1.04$	$1.42$	$0.40$	$0.78$	$0.31$	$0.49$	$0.97$	$0.64$	$1.09$	$2.11$
$1.39$	$1.41$	$0.38$	$0.86$	$0.31$	$0.44$	$0.86$	$0.66$	$0.97$	$2.10$
$1.06$	$1.38$	$0.41$	$0.77$	$0.27$	$0.50$	$0.80$	$0.64$	$0.83$	$2.09$
$0.93$	$1.40$	$0.46$	$0.82$	$0.35$	$0.51$	$0.77$	$0.62$	$1.16$	$2.09$
$1.31$	$1.44$	$0.44$	$0.82$	$0.37$	$0.39$	$0.87$	$0.65$	$0.98$	$2.09$
$1.01$	$1.38$	$0.48$	$0.78$	$0.30$	$0.46$	$0.87$	$0.60$	$0.98$	$2.08$
$1.12$	$1.39$	$0.41$	$0.81$	$0.33$	$0.53$	$0.72$	$0.65$	$0.98$	$2.08$
$1.10$	$1.53$	$0.47$	$0.74$	$0.37$	$0.50$	$0.94$	$0.61$	$1.11$	$2.08$
$1.08$	$1.44$	$0.44$	$0.77$	$0.35$	$0.48$	$0.79$	$0.62$	$0.96$	$2.07$
$1.04$	$1.16$	$0.36$	$0.79$	$0.36$	$0.55$	$0.90$	$0.65$	$1.24$	$2.05$
$0.91$	$1.39$	$0.42$	$0.74$	$0.27$	$0.40$	$0.95$	$0.58$	$0.93$	$2.05$
$1.10$	$1.39$	$0.47$	$0.76$	$0.31$	$0.51$	$0.69$	$0.60$	$0.93$	$2.05$
$1.20$	$1.42$	$0.43$	$0.76$	$0.36$	$0.51$	$0.73$	$0.63$	$0.94$	$2.05$
$1.14$	$1.45$	$0.46$	$0.79$	$0.43$	$0.42$	$0.84$	$0.63$	$1.00$	$2.04$
$1.02$	$1.36$	$0.41$	$0.76$	$0.32$	$0.18$	$1.00$	$0.64$	$1.12$	$2.03$
$0.98$	$1.36$	$0.38$	$0.79$	$0.32$	$0.23$	$0.83$	$0.66$	$1.07$	$2.03$
$1.01$	$1.43$	$0.45$	$0.74$	$0.39$	$0.48$	$0.74$	$0.59$	$0.91$	$2.03$
$1.33$	$1.38$	$0.40$	$0.78$	$0.31$	$0.48$	$0.69$	$0.59$	$1.02$	$2.02$
$1.16$	$1.44$	$0.42$	$0.82$	$0.14$	$0.49$	$0.66$	$0.63$	$1.15$	$2.02$
$1.01$	$1.41$	$0.42$	$0.70$	$0.36$	$0.56$	$0.76$	$0.56$	$1.03$	$2.01$
$1.04$	$1.39$	$0.42$	$0.75$	$0.35$	$0.44$	$0.63$	$0.62$	$0.99$	$2.01$
$1.05$	$1.24$	$0.40$	$0.74$	$0.33$	$0.43$	$0.70$	$0.65$	$1.24$	$2.01$
$0.90$	$1.35$	$0.38$	$0.71$	$0.34$	$0.44$	$0.79$	$0.65$	$0.83$	$2.01$
$0.90$	$1.30$	$0.46$	$0.72$	$0.36$	$0.49$	$0.72$	$0.57$	$0.96$	$1.98$
$0.94$	$1.46$	$0.40$	$0.65$	$0.27$	$0.48$	$0.83$	$0.63$	$1.14$	$1.98$
$1.20$	$1.37$	$0.38$	$0.69$	$0.41$	$0.59$	$0.79$	$0.56$	$0.94$	$1.98$
$1.08$	$1.20$	$0.43$	$0.70$	$0.41$	$0.45$	$0.80$	$0.63$	$0.83$	$1.97$
$1.07$	$1.43$	$0.39$	$0.65$	$0.33$	$0.43$	$0.63$	$0.55$	$1.20$	$1.95$
$0.78$	$1.35$	$0.39$	$0.71$	$0.43$	$0.42$	$1.02$	$0.66$	$1.19$	$1.92$
$1.44$	$1.47$	$0.35$	$0.73$	$0.32$	$0.46$	$0.92$	$0.64$	$1.26$	$1.91$
$1.49$	$1.35$	$0.38$	$0.78$	$0.33$	$0.55$	$0.69$	$0.63$	$1.23$	$1.90$
$1.66$	$1.43$	$0.41$	$0.80$	$0.35$	$0.69$	$0.77$	$0.61$	$1.00$	$1.88$
$1.68$	$1.42$	$0.38$	$0.82$	$0.33$	$0.52$	$0.90$	$0.65$	$1.18$	$1.83$
$0.71$	$1.40$	$0.49$	$0.58$	$0.38$	$0.47$	$0.64$	$0.55$	$0.87$	$1.81$
$1.38$	$0.69$	$0.30$	$0.72$	$0.42$	$0.57$	$0.96$	$0.67$	$1.44$	$1.75$
$0.66$	$1.28$	$0.23$	$0.69$	$0.36$	$0.51$	$0.67$	$0.60$	$1.25$	$1.75$

S3.4 Value of the yield

Figure S6 shows the value of the yield for the 40 final Nash equilibria.