Nash equilibrium can resolve conflicting maximum sustainable yields in multi‐species fisheries management

&NA; The current fisheries management goals set by the European Commission states that fish stocks should be harvested to deliver maximum sustainable yields (MSY) and simultaneously, management should take ecosystem considerations into account. This creates unsolved trade‐offs for the management of the stocks. We suggest a definition of a multi‐species‐MSY (MS‐MSY) where no alternative fishing mortality (F) can increase yield (long term) for any ecologically interacting stock, given that the other stocks are fished at constant efforts (Fs). Such a MS‐MSY can be solved through the game theoretic concept of a Nash equilibrium and here we explore two solutions to this conflict in the Baltic Sea. We maximize the sustainable yield of each stock under two constraints: first, we harvest the other stocks at a fixed F (FNE); second, we keep the spawning stock biomasses of the other stocks fixed [biomass Nash equilibrium (BNE)]. As a case study, we have developed a multi‐species interaction stochastic operative model (MSI‐SOM), which contains a SOM for each of the three dominant species of the Baltic Sea, the predator cod (Gadus morhua), and its prey herring (Clupea harengus), and sprat (Sprattus sprattus). For our Baltic Sea case, MS‐MSYs exist under both the FNE and the BNE, but there is no guarantee that point solutions exists. We found that the prey species' spawning stock biomasses are additive in the cod growth function, which allowed for a point solution in BNE. In the FNE, the herring MSY was found to be relatively insensitive to the other species' fishing mortalities (F), which facilitated a point solution. The MSY targets of the BNE and the FNE differ slightly where the BNE gives higher predator yields and lower prey yields.


Introduction
The European Commission has decided to manage its fisheries sustainably with regards to the whole ecosystem (European Commission, 2002). This was followed by the member countries of the United Nations at the World summit on sustainable development at which they agreed to restore stocks to produce maximum sustainable yields (MSYs) by 2015 (United Nations, 2003). An ecosystem-based fishery management (EBFM) and a restoration to MSYs are now the current directive from the European Commission to the International council for the exploration of the sea (ICES, 2009). The problem is that multiple single-stock MSYs of ecologically linked species are conflicting. Although single stocks are harvested sustainably, it is likely that the multi-species complex is not (Pikitch et al., 2004;Mackinson et al., 2009). Although there has been a general acceptance of the rationale behind EBFM, there is no consensus regarding how it should be implemented (Essington and Punt, 2011). This has led to an increased interest in the so-called surplus models in order to provide a FMSY (a fishing mortality, F, leading to MSY) concept for an entire ecosystem. Fundamental to surplus models is the understanding of the drivers of stock productivity. These are sometimes described as a triad of exploitative, trophodynamic, and biophysical drivers (Link et al., 2010).
To develop biological reference points (BRP) for EBFM the state of a stock has to be assessed, sometimes incorporating trophodynamic and/or biophysical drivers (Link et al., 2010). Reference points have to be developed in relation to the drivers of productivity, usually in the terms of a FMSY analysis in which exploitation is varied and other drivers are kept constant (ICES, 2011b). Thus, BRP are specific to the drivers' state, which must be estimated directly or indirectly, e.g. as effects on individual growth or natural mortality (WD6, ICES, 2011a). A study of FMSY reference points using a multi-species model for the Baltic Sea reveals a dilemma for management (Gislason, 1999). To maximize total sustainable biomass yield when fishing on predator-prey species, the most rationale solution would be to eradicate the predator because of the higher productivity of lower trophic levels. If instead the net economic benefit had to be maximized, the fisheries on the predator should be sustained because of its higher market value compared to the forage fish (Clark, 1990). Maximization of monetary value is not in line with an MSY objective and eradication of predatory species challenges the concept of EBFM, not to mention the conservation of species, biodiversity, and ecosystem function in general. Hence, for the implementation of an EBFM, there is the urge to define an MSY for each species, which does not impair the productivity of interacting species. Such a framework of multispecies reference points (e.g. MS-MSY) would thus comprise a set of fishing mortalities or biomasses that meet the MS-MSY criteria. Currently, no definition of this type exists.
An outset definition of an MS-MSY fisheries can be that there is no alternative fishing mortality (F) that can increase yield (long term) for any of the ecologically interacting stocks, given that the other stocks are fished at a constant effort (F). A conflict framed in this way can be solved as a Nash equilibrium (Nash, 1951), which is one of the foundations in the scientific field of game theory. In this way, the interest of maintaining each stock at a maximum yield, while taking into account the conflict between predator-prey and competing species, is ensured.
The urge for MS-MSY reference points is an urge for a management tool. However, the game-theory applications in fisheries have so far been dealing with the incentive for joint management (Bailey et al., 2010). Seminal work analysed the two-player conflict over one resource (Mirman, 1979;Levhari and Mirman, 1980), which was extended to two players and two resources (Fisher and Mirman, 1996), and one resource with multiple players (Okuguchi, 1981) that could optionally form coalitions (Denisova and Garnaev, 2008). For the suggested MS-MSY, the MSY objective is already agreed on within the EU coalition, and a starting point for the analyses. In contrast to the conflict between harvesters, there are not necessarily any physical players with conflicting interests in MS-MSY conflicts. The conflict is inherent in the objective to simultaneously achieve MSY of multiple ecologically interlinked stocks. The MSY conflict can still be at hands of one managing agency or even one individual. To our knowledge, this game-theoretical approach has not been proposed previously to solve the multi-species MSY conflict or as a definition of MS-MSY reference points.
An MSY-analysis concept based on a stochastic operative model of a stock (SOM-MSY) has been developed earlier (Holmgren et al., 2012). The SOM-MSY returns targets and reference points for harvesting in relation to ecological and environmental drivers (Botsford et al., 1997). The SOM-MSY has been used for the ICES advice of the central Baltic Sea herring and Baltic Sea sprat in 2011 (ICES, 2011a). The method of using stochastic stock models is ranked as the second most preferred, after a management plan, for deriving reference points by ICES (2011b). In this article, we have further developed the SOM-MSY to a multi-species interaction-SOM (MSI-SOM), using the Baltic Sea as a case study. The Baltic Sea is characterized by few species compared with other marine systems. Cod is the main predator on the clupeids, herring, and sprat. This fact allows us to study the effects of both predator-prey interactions and indirect competition between prey species through their common predator in the MSY estimations. Accordingly, the aim of this article is to investigate the MS-MSY as a reference point concept with the three main Baltic Sea species, cod, herring, and sprat, being a model system. It is not, at this point, to produce recommendations of catch targets for the Baltic Sea.

Methods
The MSI-SOM The model consists of three SOMs for cod, sprat, and herring stocks, respectively. Each SOM has numbers-at-age (NAA) and weight-at-age (WAA) as dynamic variables. The changes in the dynamic variables are defined by a set of four stock functions: (i) a recruitment function, (ii) a weight-of-recruits function, (iii) a natural mortality function, and (iv) a body-growth function. This model structure has been used previously to describe stock dynamic processes (Beverton and Holt, 1957;Ault et al., 2005). Some of the functions included environmental forces as independent variables.
The interactions between the species consist of predator growth dependency on prey abundance, and prey mortality depending on predator abundance. These two interactions are the most important to consider in multi-species analyses of reference points (Collie and Gislason, 2001). In the Baltic Sea, there is evidence that herring and sprat compete for food and experience density-dependent individual growth (Casini et al., 2010;Lindegren et al., 2011). However, for the more general methodological scope of this article, it was not included in the MSI-SOM. Once the functions were estimated, the MSI-SOM was used for simulations of fisheries effects on the interacting stocks. The MSI-SOM can be run deterministically, or stochastically. The stochastic run includes random errors in the stock functions and in the environmental drivers.

Data
Time series of WAA and maturity ogives for herring and sprat, which are primarily used as input in the annual stock assessment (extended survivors analysis, XSA), were retrieved from ICES (2013a). Likewise, the output of the XSA provided the time-series of the NAA and spawning stock biomass (SSB) for herring and sprat (ICES, 2013a). For cod, WAA, maturity ogives, NAA and SSB were from the state-space assessment model SAM (ICES, 2013a). The total biomass of cod used in the cod growth function is the output of the XSA (ICES, 2012). The selectivity vectors for each species are the same as those used in the short-term forecast of WGBFAS (ICES, 2013a). The results from the multi-species SMS assessment were used to fit predation mortality on the clupeids as functions of cod SSB (ICES, 2013b).
Values of the reproductive volume (RV), the water volume with salinity, and dissolved oxygen required for cod egg survival (Vallin and Nissling, 2000) include the current and historic main cod spawning areas, i.e. the Bornholm deep, the Gdansk deep and the Gotland deep. The RV is referring to May conditions until 1991, after which it is adjusted to August to reflect the delayed spawning time of the cod. A gamma distribution was fitted to RV data from 1981 to 2010 (scale ¼ 48.7 and shape ¼ 3.08). This distribution was used in the stochastic runs of the model, whereas the mean was used in the deterministic runs (Supplementary Table S1). This accounted for the reduction in the peaks of RV compared with the period 1968-1980 ( Figure 1).
Salinity and temperatures at surface (10 m) and mid-water (50-60 m) in summer were from the Swedish Meteorological and Hydrological Institute and Helsinki Commission stations BY 5 (Bornholm basin, ICES SD 25) and BY 15 (Gotland deep, SD 28). A spectral Fourier analysis of the salinity and temperature was performed to separate correlated from uncorrelated inter-annual variation in the data. A Fourier function including the three dominant frequencies was fitted by least-squares to the data ( Figure  2). The last year estimate of the fitted functions was used as a constant in the deterministic runs, and as the mean of a normal distribution error in the stochastic runs (Supplementary Table  S1). The mean square errors of the residuals around the fitted Fourier functions were used as the variance of this error.

Fitting the stock functions Recruitment
We evaluate three different types of recruitment functions: Beverton-Holt, Ricker, and the Quadratic functions, where the dependent variable is recruits per SSB (Cadima, 2003). Each recruitment function type was tested with environmental influence as additive factors or as interacting with SSB or as a combination of both. Since the Quadratic recruitment function intersects with the x-axis, potentially negative recruitment is truncated to zero in our simulations.
For cod recruitment, we explored the influence by the RV (see Supplementary Table S1 for the units of the parameters). We also included average parental weight (PW) to account for the effect of different egg buoyancy due to female weight (Vallin and Nissling, 2000). The interactive effect of PW with SSB was included as a multiplicative effect both in the numerator and in the denominator.
For clupeid recruitment, we investigated if it is influenced by the declining salinity in the Baltic Sea observed since the 1980s (H€ anninen et al., 2000). We also explore the impact of changing temperature (Margonski et al., 2010). Summer temperature and salinity from 10-m depth were used as a covariate in the clupeid recruitment functions. The interaction of salinity and temperature with SSB is put in the denominator in order to get increased competition with decreased salinity and temperature.

Weight of recruits
The weight of recruits was estimated as a function of PW, i.e. linking the weight to succeeding cohorts (Geffen, 2009). The weights of recruits were estimated as linear functions of average PW. The average PW was calculated by taking the sum of WAA, weighed by each age's share of the total SSB.

Natural mortality
Natural mortality for the clupeids was estimated with speciesspecific general linear models with age as a category variable and cod SSB as age-specific covariate. A similar approach was adopted by Pope and Macer (1991) who used a linear relationship between the average number of predators and predation mortality. The natural mortality for cod is constant ¼ 0.2.

Body growth
When we tested different functions of body growth of the species, we always included body weight as an independent variable to get an asymptotic growth (Lester et al., 2004). The cod body growth functions evaluated included sprat SSB and herring SBB (prey abundance), number of 0-year-old cod and 1-year-old cod (cannibalism), and total biomass of cod (density dependence by food depletion within a year, interference competition, or possibly cannibalism). The growth functions of the clupeids were fitted including effects of summer salinity and temperature measured at 50-60 m. It has been proposed that the availability of the clupeids' preferred food is linked to salinity and temperature (Möllmann et al., 2004(Möllmann et al., , 2005. The clupeids are naturally marine species so a lower salinity is also assumed to be directly detrimental to growth (Casini et al., 2011). An increasing temperature, on the other hand, is assumed to improve primary production and ultimately clupeid growth. Food competition and density-dependent growth of sprat and herring (Casini et al., 2010;Lindegren et al., 2011) was not included in the MSI-SOM due to the more general methodological scope of this article.
The growth functions used the difference in WAA of a cohort, i.e. growth (G) for age a at year y ¼ W a þ 1,y þ 1 À W a,y , as the dependent variable for the clupeids. For cod, trials with growth as the dependent variable exhibited cone-shaped residuals (not shown here). Residuals were uncorrelated with predicted values when growth per unit weight was chosen as the dependent variable, i.e. G a,y /W a,y .

Selecting functions and running the MSI-SOM
Alternative recruitment functions and growth functions were evaluated using the Akaike information criteria (AIC) (Burnham et al., 2011). We selected the single best function for inclusion in the SOM model. If the best function had a possibility to produce a population with unlimited growth or negative recruitment at extreme densities, we individually evaluated the succeeding functions for inclusion in the SOM (see Supplementary Tables S3, S4, and S7 for the criteria of non-inclusions). The functions for natural mortality of the clupeids and the functions relating weight of recruits to PW were fitted to data without being challenged by alternative models. The mean squares of the model errors (MSE) of all the fitted functions were incorporated into the SOM by adding them as normal distributed error terms to the functions with 0 mean and the variance equal to the MSE.
The program R (http://www.r-project.org/) was used for all statistical analyses. The quadratic recruitment function, the herring and sprat mortality functions, and the recruit-weight function were fitted using the linear model module (lm). The nonlinear least squares model (nls) was used for the Berverton-Holt and Ricker recruitment functions. With the fitted functions in place, the MSI-SOM was run in R as a deterministic or a stochastic version. For the MSY analyses, the model was run with fixed fishing mortalities. The first 100 generations of data were discarded, whereas the consecutive 50 generations were used as the basis for the mean values of yield and SSB. The MSY was defined as the peak of the yield curve, stepping through a range of runs with 0.1 increments of Fs.
Nash equilibrium MSY Nash (1951) formulated a solution to multi-player games, who states that an equilibrium of strategies exists if no player can improve its payoff given fixed strategies played by the opponents. When we apply this approach to a multi-species MSY-oriented fishery, we secure the interest of maximizing the yield (in mass) of each fishery taking into account the ecological impacts on other stocks. Technically, we can analyse the Nash equilibrium by keeping fishing mortalities or stock biomasses of interacting stocks constant. Prior to our analyses, we want to keep as many options as possible open and analyse the properties of both methods. When keeping biomasses constant this condition can be defined as when all species (i) obey where Ài denotes the other species. This means that there is an F that returns a higher yield Y than any alternative F, given constant biomasses B of the other species. Technically, the Nash equilibrium is conducted by solving SSB at MSY (B MSY ) for one species while stepping through a range of constant Bs of the other species. Keeping the Bs constant should be regarded as a technique to find MS-MSYs, when applied in reality the Bs will vary. We denote this biomass Nash equilibrium (BNE), and hence the associated yield, Y BNE , biomass B BNE and fishing mortality F BNE . Even under BNE when the opposing stock biomasses are constant, the target is still FMSY. An alternative condition can be defined when the other species are harvested at a fixed fishing mortality: We denote this condition fishing mortality Nash equilibrium, FNE. Because the FNE is not controlled with constant biomasses B, the biomasses of all species are allowed to change dynamically during calculations. We have calculated FMSYs over the historical ranges of Fs and present a graphical analysis of the FNE. This, however, is computationally intensive so to increase the applicability of the method the FNE can also be solved readily through iteration: Initial values are chosen for F C (cod), F H (herring), and F S (sprat), and the F vector that maximizes yield are calculated (vector of F 0 ). The F 0 vector is then used as input to the next iteration step. When the F vector does not change between iterations the FNE is found. The BNE and FNE analyses are performed without stochastic noise in the MSI-SOM model or in the environmental variables. Stochasticity in both the MSI-SOM and the environmental variables is included in the simulations when the reference points are evaluated. The resulting FNE reference points are compared with the reference points of systems where one of the clupeid species is fished at a high fishing mortality (F ¼ 1.0). This approximates a predator-prey (pp) pair-FNE and the differences in the results compared with the triad-FNE reveal how the number of prey species affects the reference points. More specifically, for the case of the Baltic Sea, we want to investigate how the different predation mortality of the clupeids and their effect on predator growth affects the reference points.

Recruitment functions Cod
Based on AIC model selection, the best model to define cod recruitment is the Beverton-Holt with RV as an additive factor and PW (PW C ) interacting with SSB C (Supplementary Table S2): The residuals of the function have no strong trend and exhibit a similar distribution over the range of the fitted values. The RV and the PW C add substantial explanatory value to the function (Supplementary Table S2 and Figure S1).

Herring
The best applicable model to describe herring recruitment is Beverton-Holt with salinity (S B10 ) interacting with temperature (T B10 ) and herring SSB (SSB H ). Three other models had better fit, but they had positive density dependence, which we regarded as biologically implausible (Supplementary Table S3): The residuals of the function are well behaved and the influence of salinity and temperature on the function is weak (Supplementary Table S3 and Figure S2).

Sprat
The model chosen to define sprat recruitment is the quadratic with temperature (T G10 ) as an additive factor (Supplementary Table S4):

R SSB S
¼ 22:6 Á T G10 À 9:03 ð Þ À 9:58 Á 10 À05 Á SSB S This means that the summer surface temperature must be above 9 C for successful recruitment. The influence of temperature on the function is strong (Supplementary Table S4 and Figure S3).

Mortality functions Herring
Herring mortality is affected by cod SSB when the herring is young but the influence of cod decreases markedly with increasing age (Figure 3). For the three last age classes, the linear relationship with cod is negative so for these age classes we use the average mortality with no influence from cod SBB. The effect of age and the interaction between age and cod SSB are not significant (Table 1).

Sprat
Sprat mortality is linked to the cod SSB (Figure 4), and all age classes are similarly affected. The effect of age and the interaction between age and cod SSB are not significant ( Table 2). The parameter values of the age-specific mortality functions are in Table 3.

Body growth functions Cod
The best fitting body growth function for cod is growth per unit weight as a function of sprat and herring SSB, the total cod biomass (TBM), and weight (W C ; Supplementary Table S5).

Herring
The clupeid growth functions are functions of salinity, temperature, and weight where weight can be an additive effect and interact with S and T (Supplementary Tables S6 and S7). The chosen model for herring growth includes salinity and temperature as additive terms and salinity as a factor interacting with herring weight (W H ): G ¼ À1:24 Á 10 À2 þ 3:55 Á 10 À03 Á S B50 À 1:18 Á 10 À03 Á T B50 The model produces well-behaved residuals and predicts growth well (Supplementary Figure S4), and all effects in the model are significant (Table 5).

Sprat
The chosen model for sprat growth includes salinity and weight with salinity as an additive factor: G ¼ 3:30 Á 10 À03 þ 1:03 Á 10 À03 Á S G60 À 0:295 Á W S The model produces well-behaved residuals and predicts growth well (Supplementary Figure S5), and both the salinity and the weight effects in the model are significant (Table 6).

Weight of recruit functions
The average weight of recruits (RW) were positively correlated with average PW for all three species, but only significantly so for the clupeids (Table 7). The functions are as follows: RW H ¼ 0:59 Á 10 À03 þ 0:347 Á PW H (11) RW S ¼ À9:02 Á 10 À05 þ 0:571 Á PW S Nash equilibrium reference points      ). This enables the use of a weighed SSB sum of the clupeids, here measured in herring SSB equivalents. The intersection between BMSY-isoclines plotted in the phase portrait of SSB for clupeids and cod defines the condition in (Equation 1). The BMSY-isocline for herring equivalents is plotted as a function of SSB of cod ( Figure 5). The BMSY-isocline for cod is plotted as a function of SSB of herring equivalents, but for this isocline the independent herring equivalents are on the y-axis that makes it possible to view the two isoclines in the same graph. The isoclines intersect at a BMSY of 295 t tons for cod and at a BMSY of 717 t tons in herring equivalents. The intersection defines the BNE since it is the best response (MSY) given a constant biomass of the other species ( Figure 5). The clupeid biomass can be divided into herring and sprat biomasses through FMSY analyses run with the BMSY of cod (295 t tons). For herring the BMSY is 485 t tons and for sprat the BMSY is 784 t tons (

FNE
We investigated the MSYs and associated F and SSB within the historical F-space of the species triad ( Figure 6). As expected, there is a conflict between MSYs of the three species. From the cod fishery perspective, the fishing mortality on sprat and herring should be as small as possible. From the herring fishery perspective, the fishing mortality should be as large as possible on cod and sprat. In fact, a lower cod biomass, obtained by either directly increasing the fishing on cod or indirectly by increasing the fishing on sprat (and thus affecting cod growth), would decrease predation mortality on herring. The same applies to the sprat fishery ( Figure 6, lower row). The Nash equilibrium can be found by iteratively reading the FMSY range of a species given the possible range of the other two. The iteration preferably starts with herring since it has the  Figure 5. BMSY-isoclines for cod and clupeids in terms of herring SSB equivalents (sprat SSB weighed by its impact on cod in relation to herring) with corresponding FMSYs. The F-values for sprat and herring given the cod SSB, and the F-values for cod given the herring equivalent SSB can be read on the right y-axis. The BNE is defined by the intersection of the BMSY-isoclines. The dotted vertical line indicates the corresponding FMSYs of the BNE. The pp analyses set a high fishing mortality for one of the prey species which decimates the population and approximates a two species (pp) system. The symbol -means an F of 1.0 for herring or sprat in the FMSY column and low BMSY and MSY values in the respective columns. SSBs and yields are shown in thousand tons.
narrowest FMSY range, and continues with cod that has the narrowest FMSY range given the range of herring. With two defined ranges we can read the FMSY range for sprat, and then reiterate in the order herring-cod-sprat until one FMSY-value (with a desired precision) is achieved for each species (Figure 6, upper row). The herring FMSY is very insensitive to the applied F of the other two species, the sensitivity increases with increased fishing mortality on sprat. The cod FMSY is only sensitive to herring F when F is relatively low for both herring and sprat. When sprat F is relatively high, the cod FMSY becomes more sensitive also to sprat F. Sprat FMSY shows a number of concave curves with rather uniform sensitivity over the investigated ranges.
The cod BMSY is not monotonically changing with increasing sprat F. For a given herring F it first declines with increasing sprat F, but around F ¼ 0.75 the cod BMSY increases. This can be understood with the sharp decrease in cod FMSY with increasing sprat F > 0.75. Hence, in a herring-dominated community, FMSY for cod should decrease with a resulting increase in BMSY (compared to F ¼ 0.75 for sprat; Figure 6).
We also investigated the pp-FNEs with one prey species available and the other under high fishing pressure. With two prey species available, the predator fisheries MSY is higher than if only one prey species is available (Table 8). For the prey species fisheries, the effect of the presence of a second prey species is the opposite. Their MSY increases when they are alone with the predator. The difference is small for herring, but for sprat there is a predatory release due to a smaller cod stock. The cod FMSY of the triad-FNE is higher than in the pp-FNE, with a higher MSY and higher BMSY. In contrast to the predator, the prey species' FMSYs are higher in the pp-FNE than in the triad-FNE. The clupeid stocks respond differently to this change. The herring BMSY increases when going from the pp-FNE to the triad-FNE, whereas

Cod
Herring Sprat the sprat BMSY decreases. This is probably a result of the cod's different predation profile on herring and sprat. All ages of sprat are subject to predation, whereas only the young herring is subject to cod predation. The MSY increases for both clupeids in the pp-FNE compared to the triad-FNE. It is valid to ask how these NE reference points compare to single species (SS) reference points. For the predator (here cod), the reference point will depend on the assumed or actual food situation accounted for in the SS-FMSY analysis, usually given by the WAA of cod in the catches and surveys. Based on the results from our model, cod SS-FMSY can be anything between 0.38 and 0.60. Similarly, the prey species can achieve a range of SS-FMSYs depending on the natural mortality assumed or estimated in the SS-FMSY analysis. In this case, 0.29-0.31 for herring and 0.35-0.75 for sprat. The ranges are given in Figure 6 that depicts the FMSYs for the ranges of F observed historically for the three stocks. The drawback with an SS analysis is that the dynamic effects of the interaction between the species are not taken into account. Hence, there exists no fixed target FMSY, it will change with changes in the input assumptions.

Comparison of the MS-FMSYs
The BNE gives higher MSY for cod but lower MSY for herring and sprat compared to the FNE (Table 8). The FMSYs of the BNE is lower for all three species compared to the FMSYs of the FNE. We ascribe this effect to the lack of stock biomass resilience in the BNE analysis. The higher FMSYs of the FNE reduce the BMSYs for cod and herring but not for sprat compared with BNE.
When evaluating the FNE and BNE reference points with stochastic simulations, it is clear that all populations exhibit large variation in SSBs when harvested at the equilibrium points (Figure 7). At both the FNE and the BNE cod is continually above B lim and B pa reference levels and sprat is also above B lim and B pa most of the time (ICES, 2013b). The herring is instead continually below current B pa and below B lim most of the time when harvested at the BNE. At the FNE state, herring SSBs are above B pa a few times but also sometimes below B lim and the variation in SSB appear to be larger at the FNE than at the BNE. Sprat exhibits the largest variation in SSB (400-2,000 thousand tons at the FNE) with slightly lower SSBs at the BNE than at the FNE. The variation in cod SSB appears to be similar in FNE and BNE.

Discussion
The arguments for multi-species approach vs. SS when setting fisheries targets Harvesting targets based on SS analyses can vary widely depending on the assumptions of predation and body growth, something we have demonstrated here on the simple Baltic Sea cod-herringsprat system. In a short-term perspective, current natural mortalities and predator WAA can be used, but this would require a reestimate of the reference points when the conditions change (Collie and Gislason, 2001). To our knowledge, the dynamics of such a management strategy with a moving target has not been analysed, with the risk that stock productivity is suboptimal or stocks even jeopardized. SS analyses simply does not incorporate the mutual dynamic effects a predator and prey have on each other which, as shown here, can be used to define an equilibrium point. Christensen (1996) reviewed the management of pp systems and concluded that "SS models will not suffice." The same conclusion was reached in a sensitivity analysis of BRPs (Collie and Gislason, 2001). They analysed the effects of historical variation in predation mortality on sprat and food availability for cod. The MSY reference points for the sprat were very sensitive to predation mortality, whereas the MSY reference points for cod were less sensitive to food abundance. The multi-species approach must be an integrated part of ecosystem-based management. An unbalanced fishing effort on the species in the ecosystem can have undesired effects. Intense fishing on predator species can increase the competition between prey species, of which the subordinate species will decline in abundance or disappear (Botsford et al., 1997). In our multispecies approach, we see that the higher F of cod in the FNE (compared with BNE) is accompanied by higher Fs also on the prey species, hence reducing the risk of elevated competition between the prey species.

Alternative definitions of multi-species targets
The quest for suitable targets in fisheries where stocks are ecologically linked is not new. Several authors have suggested the objective of maximizing total economic yield for multi-species fisheries. In a simple pp system where the predator's growth rate is entirely determined by its intake rate of prey, Clark (1990) analysed the total sustainable yield or revenue from harvesting. He shows that to maximize the total sustainable yield or revenue either the predator should be eliminated and only the prey exploited, or, if the predator is sufficiently valuable, only the predator should be exploited. The objective to maximize the profit changes the ecosystem to an "agricultural configuration" in which only the profitable species are supported and their predators and competitors are eliminated (Christensen and Walters, 2004). In some cases, supposedly when market values are not too different, total economic yield can balance the harvest rates between the species compared to a strategy of maximizing total yield biomass (Gislason, 1999).
The objective to maximize the total yield biomass has been suggested as a multi-species target. It usually has the consequence of fishing down the food chain (Gislason, 1999;Christensen, 1996) but can be combined with precautionary restrictions so that stocks are not jeopardized (Guillen et al., 2013). Both maximizing economic yield and total yield biomass have no inherent mechanisms to prevent unbalanced exploitation among stocks in the way the NE approach does. May et al. (1978) proposed a precautionary approach to multispecies fisheries at which the predator would be fished at MSY and the prey species should be fished so that they are not depleted (Christensen, 1996). By keeping the biomass of the prey species intact when analysing the MSY of the predator, it is basically the first step of a BNE analysis. The second step of finding the BNE is to analyse the MSY of the prey given different BMSYs of the predator. If the second step is not taken, the predator MSY would be somewhat arbitrary.
A similar idea to deal with conflicts in multi-species MSYs is suggested by Beddington and Cooke (1982) where the difference between the amount of prey eaten by predators in the natural state and in a state where the predator is depleted gives a surplus of prey which can be harvested. Their analysis is based on finding the theoretical maximum MSY for a species in a multi-species system given fixed harvest regimes for the other species. This procedure is similar to the initial procedure of the FNE analysis. In FNE, however, a similar type of procedure is done for each species in relation to the other species which is then further analysed to find the state at which the targets of all species exist simultaneously. Beddington and Cooke (1982) show that given that the predator is harvested at constant fishing effort, a higher stable yield of prey can be taken than in the case of no harvest of the predator. However, not all prey surplus can be fully exploited due to the system loosing stability. Since the prey fishery in their system is not aiming for MSY, it is not in line with the MSY directives. By analyzing the yield planes of both the BNE and the FNE solutions, it is clear that our approach gives the same tendency as in the analysis of Beddington and Cooke (1982): a decreasing cod population size (directly in BNE or indirectly through a high F in FNE) results in increasing SSBs and yields at MSY for the clupeids.
Multi-species reference points have also been determined from the analysis of production models, i.e. models that combines the production from body mass growth of the fish and population growth in numbers (Sissenwine and Shepherd, 1987). However, the natural mortality of a prey caused by ecological interactions with a predator, and the effect of prey size and abundance on body growth rates of the predator are not addressed. These models, therefore, miss to incorporate the dynamic interactions of the pp system investigated here.
The MS-MSY is a game of continuous strategies of different Fs, and the 'sustainability' of the MSY analyses is one of the premises, which prevents strategic temporal resource allocation at the detriment of others (Sumaila, 1999). The previous gametheoretical studies on fisheries have exclusively investigated the conditions for cooperative management framed as a tragedy of the commons and a prisoners' dilemma game and resolving conflicts between fleets (Sumaila, 1999;Hannesson, 2011;Bischi et al., 2013). However, in a study of international conflicts in the Baltic Sea fisheries, a secondary aim was to maximize the payoff in a 50-year perspective, which has the potential to be sustainable (Nieminen et al., 2015). The economic payoff was maximized by a coalition between the participating countries (grand coalition), in which the fishing mortalities were allowed to vary over time and be freely distributed across countries and stocks. Average fishing mortalities turn out balanced (FMSY % 0.25, 0.3 and 0.23 for cod, herring and sprat, respectively), which depend on the higher payoffs of the predator. Other settings of the economic parameters may lead to intentional fishing down the food-chain or eradication of less profitable species. The aim to maximize the pay-off in a coalition setting does not guarantee the sustainability of all stocks and thereby does not comply with a MS-MSY.

Recommendations
Our result points out the state of the three dominating species in the Baltic Sea at which they simultaneously produce MSYs. The presented FMSYs, associated yields, SSBs and trigger points only apply to the system being in any of the two MSY states described. The analysis does not provide a strategy for how to advance the stocks to the state where they are capable of producing MSY. To apply the proposed FMSYs on the current stocks could prevent the stocks from reaching BMSY, as long as this has not been analysed.
The proposed joint MS-MSYs are conditional on the estimated resilience of the stocks. Although these are based on ICES data and up-to-date knowledge of stock biology and environmental drivers, resilience is sensitive to density-dependent processes and our knowledge about these can improve significantly with new research. Even if we have picked the most likely functions based on AIC, some have alternatives that are almost as likely. We therefore suggest that MS-MSYs are combined with an adaptive management approach with a regular update of stock biology processes. In addition, current model parameters are estimated on historical data during which the Baltic Sea ecosystem has undergone dramatic changes. It is not certain that the functions are valid for the entire data period.
Quite naturally, we found that MSY and BMSY of the predator increase with the number of available prey species. Similar results have been shown in other studies (Smith et al., 2011). For the Baltic Sea herring, it has previously been shown that a higher cod abundance leads to a lower MSY and FMSY (Holmgren et al., 2012). The same causal relationship is found here for both herring and sprat.
Both Nash equilibrium methods deliver MSYs for all species under the constraint that all other species must simultaneously be at MSY. They would both protect the interest of maximizing the utility of each fishery and have a preserving effect on the stocks. Looking at the stochastic variation in SSBs at the suggested FMSYs, it is evident that FNE maintains herring and sprat at higher SSBs at the same time as cod is kept at slightly lower biomasses compared to BNE.
From a fisheries perspective, it is certainly easier to control fishing mortality than to control biomass levels and since the BNE assumes constant SSBs in one of the interacting trophic levels the basic assumptions of the equilibrium may not hold as much applicability as the FNE. From the Nash equilibrium theory, one can argue that FNE is a true equilibrium of strategies because stocks are actually managed by fishing constraints (e.g. regulation of the fishing mortality), and not biomass per se. These circumstances, therefore, suggest the use of FNE in favour of BNE, but see Holmgren et al. (2014) for a discussion on biomass vs. effort target.
The FNE is not always necessarily a point solution where all species simultaneously can be harvested at MSY. Alternative settings can result in a system where the FMSY planes have several intersections. Ecosystems with more interacting species are also more likely to have multiple intersections of MSY planes. In such cases, the FNE will provide ranges of FMSYs instead (not shown). The BNE method will provide a point solution at the intersection of the BMSY-isoclines as long as the biomasses of the prey influence predator growth additively. Other interactions between or within the two trophic levels, or an increased number of species, will most likely provide MSY ranges. The FNE (and BNE) can in principle be applied on larger ecosystems, but we should expect outputs of MSY reference ranges, rather than reference points to be more common as the number of species increases. MSY ranges may even be preferred by management bodies, as it has become in the Baltic Sea and the North Sea (European Commission, 2015).

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.

Funding
This work was financed by the Swedish Research Council FORMAS project "Towards sustainable fisheries by ecosystembased management."