Maximum sustainable yield from interacting fish stocks in an uncertain world: two policy choices and underlying trade-offs

The case of fisheries management illustrates how the inherent structural instability of ecosystems can have deep-running policy implications. We contrast ten types of management plans to achieve maximum sustainable yields (MSY) from multiple stocks and compare their effectiveness based on a management strategy evaluation (MSE) that uses complex food webs in its operating model. Plans that target specific stock sizes ($B_{\text{MSY}}$) consistently led to higher yields than plans targeting specific fishing pressures ($F_{\text{MSY}}$). A new self-optimising control rule, introduced here for its robustness to structural instability, led to intermediate yields. Most plans outperformed single-species management plans with pressure targets set without considering multispecies interactions. However, more refined plans to"maximise the yield from each stock separately", in the sense of a Nash equilibrium, produced total yields comparable to plans aiming to maximise total harvested biomass, and were more robust to structural instability. Our analyses highlight trade-offs between yields, amenability to negotiations, pressures on biodiversity, and continuity with current approaches in the European context. Based on these results, we recommend directions for developments of EU fisheries policy.


Four questions about multispecies MSY
Fisheries management aiming to attain maximum sustainable yield (MSY) from multiple interacting stocks is considerably more complicated than single-stock management (Pope 1976;Bulgakova and Kizner 1986;Collie et al. 2003;Walters et al. 2005;Matsuda and Abrams 2006;Geček and Legović 2012;Houle et al. 2013;Voss et al. 2014;Thorpe et al. 2016). A priori it is not even clear what the best translation of the MSY objective for a single, isolated stock is to cases with multispecies interactions, i.e. with feeding and competitive interactions among species (for simplicity, we use "stock" interchangeable with "fish species" in this paper). Even when a management objective considering multispecies interactions is defined, attaining it can be difficult because these interactions are generally not well known. It is therefore not surprising if legislation aiming at MSY acknowledges the role of multispecies interactions in principle, but tends to play it down. Examples are § §301.a.3, 303.b.12 of the Magnuson-Stevens Act (U.S. Department of Commerce 2007) or Article 9.3.b of the Common Fisheries Policy (EU 2013), but see, §13(2) of New Zealand's Fisheries Act (Parliamentary Council Office 2014). In the current practice of fisheries management, multispecies interactions among managed stocks play only a minor role. Noteworthy exceptions are cases where stock assessments take account of the dependence of predation mortality on the abundances of other species (Gjøsaeter et al. 2015;ICES 2014b).
To contribute to the development of management practices mindful of multispecies interactions we ask here: Q1: How can the singe-species MSY objective be translated into the multispecies case? Q2: Which strategies are suited to achieve such objectives? Q3: How do these options compare with regards to the yields they achieve, the degree of collaboration among players required to reach the objectives, pressures on biodiversity, and their political acceptability?
Q4: How much can be gained from multispecies management in comparison with management disregarding ecological interactions?
Three conceivable answers to each Q1 and Q2 are given in the following two sections. Hence there are two high-level policy choices to be made. The first choice is of the management objective, the second of the type of strategy employed to achieve it. Each of the resulting nine variants we call a management plan.
To answer Q3, we performed formal management strategy evaluations (MSE, Walters and Hilborn 1976;Hilborn and Walters 1992) of these plans. That is, the outcomes of the plans were evaluated by applying them to fisheries in hypothetical ecosystems, described by an operating model, that are not fully known to the manager. To answer Q4, evaluation scores were compared with those for management that aims at MSY for the same set of stocks while disregarding multispecies interactions, modelled after current EU practice.

First policy choice: the management objective
There is no unique way of translating the single-species MSY objective to the multispecies case. Maximisation of yield from one stock will generally require different strategies than maximisation of yield from another. In a simple predator-prey system, for example, the maximisation of yield from the prey requires culling the predator, while by not exploiting the prey yield from the predator can be maximised (May et al. 1979;Clark 1990).
The three types of high-level objective we consider here are: -Nash Pressure: To fish all exploited stocks at such rates that changes in the exploitation rate of any single stock cannot increase the long-term yield from that stock.
-Nash State: To keep all exploited stocks at such sizes (e.g. in terms of spawning stock biomass) that changes in the size of any single stock cannot increase the long-term yield from that stock.
-Total Yield : Maximisation of the summed long-term yield from all exploited stocks.
These objectives are generally equivalent only in absence of multispecies interactions. The objectives Nash Pressure and Nash State are two ways of implementing the idea of "maximisation of yield from each stock separately". They correspond to the Nash equilibrium outcome in game theory (Osborne and Rubinstein 1994), where no player of a game could improve its gains by changing its moves. Nash equilibria are traditionally understood as arising naturally when players are not collaborating. For the Nash Pressure objective the players are hypothetical fleets, each targeting one specific stock, and the permitted moves changes in their exploitation rates. For the Nash State objective the players are hypothetical managers of individual stocks and their moves are the stock sizes they aim at. The Total Yield objective corresponds a situation where players chose their moves such that the total gain by all players is maximised. Attaining this objective generally requires collaboration or enforcement through governing institutions. Figure 1 depicts these three hypothetical games and objectives.

Second policy choice: the strategy
The management strategies that we consider differ by the structure of the corresponding harvest control rules (Deroba and Bence 2008). All rules ultimately prescribe the total catch from each stock for a given time period, e.g. a year. This corresponds to total allowable catches (TAC) if one assumes that allowances are fully used. Combined with estimates of stock sizes, rules for catches can be formulated in terms of rules for fishing mortality rates. Depending on how these rates are determined from managers' knowledge of stock sizes and interactions, we consider: -Pressure Target Control, where fishing mortality rates are kept fixed at target values deemed to be consistent with the selected objective.

TOTAL YIELD
NASH STATE NASH PRESSU RE Figure 1: Illustration of the first high-level policy choice as an option between three hypothetical games. The games correspond to the three high-level management objectives we consider: Nash Pressure, Nash State and Total Yield. The Nash Pressure players are managers of hypothetical fleets targeting a single stock, each aiming to maximise yield from that stock by adjusting its exploitation rate. The Nash State players are hypothetical managers of individual stocks, aiming to maximise yield from that stock by adjusting its equilibrium abundance. The last game, maximisation of total yield, requires collaboration between all players. Infographic design by Georgia Bayliss-Brown.
-State Target Control, where fishing mortalities are continuously modified such as to adjust stock sizes to target values deemed to be consistent with the selected objective.
-Self-Optimising Control, where neither targets for fishing mortalities nor stock sizes are fixed, and instead target fishing mortalities are varied according to homogeneous linear functions of stock sizes. These functions are chosen such that the expected resulting equilibrium state is consistent with the selected objective.
Formal mathematical definitions of these objectives are given in Supplementary Material, SM S1.

Simplifications, complications, negotiations
To simplify analysis and discussion, and to isolate the particular implications of multispecies interactions, a number of complications relevant in the practice of management are not covered in our analysis. First, plans do not generally include constraints resulting from conservation objectives. There are several conceivable implementations of such constraints, among others those associated with F PA , B PA , or B trigger reference points (ICES 2014a; Task Force on Multiannual Plans 2014), and we would not want to confound the side effects of specific implementations with the implications of the policy choices studied here. Instead, we determine biodiversity impacts of the unconstrained management plants. Second, neither technical interactions among fleets, that is, complications due to limitations to stock-selectivity of fishing, nor the resulting issues surrounding by-catch and discards are considered. Third, we disregard differences in market value of species and other factors differentiating between biomass and economic yields. Forth, environmental stochasticity, including recruitment variability, is not considered, despite its importance for stock assessments and decision making. Real-world management must take all these and many other complications into account and adapt objectives and strategies accordingly. The plans we consider could be modified to this end, at moderate cost in complexity and presumably with little effect on the trends we find. However, managing each of these complications involves value judgements, and presumably requires negotiations among stakeholders. Rather than prejudicating such negotiations, one can ask how easy the negotiations would be when based on the different objectives we investigate. By definition, long-term yield is maximal for each player under the Nash Nash Pressure objective, and so will not decline much when the player fishes at slightly different rates-and similarly for Nash State (impacts on yields from other stocks, however, can be larger!). Stakeholders will therefore be open to negotiations of their targets. This intuition is captured by proposals to define MSY as a "range" of targets rather than a particular "point" (Task Force on Multiannual Plans 2014), and to use this flexibility to take technical constraints into account. Achievement of the Total Yield objective generally involves trading the yield from one stock off against that from another, generating winners and losers among fleets (Voss et al. 2014). Negotiations aiming for Total Yield can therefore be more difficult. The question arises what the loss in yields would be when opting for a Nash equilibrium to avoid these difficulties. This is addressed below.
As a technical simplification, required for computational reasons, our models do not structure populations by size or age. For food-web interactions, corresponding population dynamics can be derived from size-structured models using projection methods (Rossberg and Farnsworth 2011;Rossberg 2012b). However, this means we disregard size-selectivity effects. We consider this a legitimate simplification, because, in a first approximation, optimisation of size selectivity and of fishing mortality for given size selectivity are separate problems (Law and Grey 1988; Getz 2012; Scott and Sampson 2011).

Structural instability of natural and modelled ecological communities
A particular concern of this study are the management implications of structural instability of ecological communities (Rossberg 2013, Chapter 18). By structural instability we mean that, in a community of interacting populations, small changes in environmental conditions or external pressures can lead to large changes in population sizes (Yodzis 1988). Structural instability manifests itself, among others, through the difficulties fisheries modellers experience in parameterizing models of interacting species such that these reproduce observed community structure and dynamics. Small deviations of model structure or parameters from reality can lead to very different model states. While the phenomenon as such is well known (Andersen and Ursin 1977;Flora et al. 2011;Hill et al. 2007;Novak et al. 2011;Yodzis 1988), the likely underlying ecological mechanism has only recently been revealed (Rossberg 2013, Chapters 14-18). Because structural instability increases with species richness (Novak et al. 2011), naturally assembled communities eventually saturate at richness levels where structural instability is so large that the perturbation from one species' invasion leads, on average, to extirpation of one other species. Structural instability poses a dilemma to every fisheries modeller interested in the implications of multispecies interactions. If models have less structural instability than in nature, pressure-state relations might not be correctly represented, leading to incorrect projections for MSY. If, on the other hand, the model's structural instability is as large as in nature, the model will be difficult to parameterize. Besides, comparing a model's structural instability to reality is difficult. In our MSE, these problems are avoided by abandoning the idea of using an operating model representative of a specific natural community-for our general analysis, this might not even be desirable. Instead, model communities are constructed through a random process that mimics natural assembly and turnover of aquatic food webs. This leads to the kind of community saturation thought to be responsible for structural instability in nature (Borrelli et al. 2015). Since, in addition, the resulting communities have macroecological properties similar to marine communities (Fung et al. 2013), we expect them to exhibit structural instability of a degree similar to that of marine systems.
We introduce below a method ("conservatism") capable of mitigating the adverse effects of structural instability on management outcomes. In essence, it consists in avoiding management targets that are too different from the current community state to be reliably described by inherently inaccurate management models. Instead, management gradually adapts targets to meet the objective.
It turns out that the answers to above questions Q3 and Q4 are essentially determined by structural instability. Key results of our MSE can be anticipated mathematically from structural instability and general principles of community dynamics alone, as we highlight below by referring to corresponding mathematical considerations developed in Supplementary Material. The good agreement between general theoretical expectations and simulations using a complex model suggests that these expectations are sufficiently robust to be fulfilled in the real world as well. This justifies us in deriving recommendations for practical fisheries management from this study.

Operating model
The operating model representing "reality" in our MSE is the Population Dynamical Matching Model (PDMM; Fung et al. 2013;. It describes structure and dynamics of aquatic food webs with species resolution. Species in the PDMM have different maturation body masses, which determine, through allometric scaling laws (Yodzis and Innes 1992), maximum growth rates and consumptionindependent loss rates (metabolic losses and non-predation mortality). The consumption of a resource species by a consumer species is modelled through a Type II functional response, modified to model prey-switching following van Leeuwen et al. (2013). For any two species in the PDMM, their trophic interaction strength depends on their relative body masses and the match between two kinds of abstract traits that characterise them as consumers and resources. Size preference is parameterized through a population-level predator-prey size ratio window (Rossberg 2012b). For a full description of the model, see Rossberg (2013, Chapter 22) and SM S2.
In the parameterization used here, the model spans approximately five trophic levels and resolves species and their interactions over 17 orders of magnitude in maturation body size (median range 10 −15 to 10 1.8 kg). Following Fung et al. (2013), species with maturation body sizes larger than 1 g are interpreted as fish. When numbering the S species in the model from the largest to the smallest, so the first S F species are the fish, this leads to population dynamics of the general form where t is time. The momentary growth rate g i (. . . ) of the population biomass B i of each species i depends on its direct biological interactions with other species. The parameters F i represent exploitation rates, and for each fish species i the product Y i = F i B i is the yield from that species per unit time. Exploitation rates are proportional to adult fishing mortalities and typically attain numerically similar values (Shephard et al. 2012), hence our use of the symbol F . PDMM model communities are generated by simulating the natural processes of community assembly and turnover by iteratively adding species to the model and removing those that go extinct, until S fluctuates around some equilibrium. This is done with all exploitation rates (F i ) set to zero. Each MSE was performed using 37 different samples from this community, taken after every 10.000 species additions and therefore composed of largely independent sets of species (of 40 communities so sampled, the first 3 were discarded as burn-in). In our parameterization (SM S2), each sample contained around 2000 species, of which between 9 and 38 were fish (mean 22.4).

Management model
As the management model, we used the multispecies extension of Schaefer's (1954) surplus production (or multispecies Lotka-Volterra) model (Pope 1976(Pope , 1979 where the surplus productivities r i and interaction coefficients G ij (1 ≤ i, j ≤ S F ) are constants. This model has frequently been studied for management applications (Bulgakova and Kizner 1986;Collie et al. 2003;Gaichas et al. 2012;ICES 1989;Pope 1979). It has the advantages of compatibility with the abstractions employed in the operating model, Eq. (1), formal simplicity, and a low number of fitting parameters compared to processed-based models. Management models used in practice are often more complex. SM S3 describes a method to calibrate the parameters of the management model (r i , G ij , 1 ≤ i, j ≤ S F ) such that it approximates the dynamics of fish populations in the operating model, Eq. (1), for states similar to the current. The method is based on the assumption that the responses of non-fish species to changes in fish populations are fast compared to these changes. It should therefore work best for fish populations close to equilibrium, where they change slowly.
A mathematical analysis of the calibration method reveals that structural instability of the operating model can lead to inaccurate approximations of the operating model by the management model and to structural instability of the management model itself (SM S4). Our MSE is designed to capture the implications for management outcomes of these effects and of incomplete representation of reality by the management model, because similar issues can arise for management models used in practice. For simplicity, the MSE does not consider measurement errors of population-dynamical parameters or implementation errors, and it assumes perfect knowledge of all B i and and F i .
Management plans are applied in cycles simulating adaptive management. At the beginning of each cycle, the management model is calibrated to approximate the dynamics of the operating model for states similar to the current state. The multispecies harvest control rule (mHCR, see below) corresponding to the plan is then (re)parameterized using the management model and used throughout the management cycle. A rather long period (50 years) is chosen for these cycles to allow the operating model to reach an equilibrium, because, with our simple calibration algorithm, this results in better fit of the management model to the operating model in the next cycle. Conceivable effects of shorter management periods were not investigated.

Multispecies harvest control rules
The three management strategies we consider differ by the formulae used to determine exploitation rates (mHCR). The free parameters in these formulae (F MSY,i , B MSY,i ,Ĝ i ) are chosen such that, if the management model was the correct model of reality, the corresponding equilibrium state of the system would meet the objective of the plan. Corresponding analytic expressions are derived in SM S5 following Pope (1979).
The simplest case is Pressure Target Control, where exploitation rates are kept constant at for each fish species i. State Target Control is implemented by a rule where the F 0,i are defined in the following paragraph, B MSY,i are the targeted stock sizes, and the parameter T (dimension Time) depends on how fast management attempts to reach the target. Here we chose T = 1 yr rather small, to model "fixed escapement" management as recommend by bio-economic analyses with low discount rate. The operation max [0, . . .] replaces any negative exploitation rate by zero. As the (time dependent) neutralising exploitation rate F 0,i we define the exploitation rate that would keep stock i at its current size, provided all other stocks and long-term environmental conditions remain in their current state as well. The values F 0,i can be computed independently of the management model Eq.
(2), in practice during yearly stock-assessments using standard methods, or their multispecies extensions. In our MSE we set F 0,i = g i (B 1 , . . . , B S ), with g i (B 1 , . . . , B S ) as in Eq. (1a). As a result, Eq. (1a) become effectively whenever F i , as given by Eq. (4), is positive. That is, the State Target Control rule is designed to achieve exponential approach of all stocks to their target sizes B MSY,i at a rate 1/T . By modifying Eq. (4) and/or the parameter T , other dynamics for approaching B MSY,i can be obtained. For the type of plan we recommend in Conclusions below, one such variant will briefly be discussed. For the Nash Pressure and Nash State objective, the Self-Optimising Control rule has the form where the constantsĜ i depend on the interaction coefficients G ij of the fitted management model, Eq. (2). On average, approximately 4% of the calculatedĜ i were negative for both Nash Pressure and Nash State and therefore F i set to zero. Self-Optimising Control for the Total Yield objective requires The coefficients G ji entering Eq. (7) are directly given by interaction strengths from the fitted surplus production model, Eq. (2). Because the matrix (G ji ) entering Eq. (7) is the transpose of the interaction matrix G = (G ij ) from Eq.

Single species harvest control rule
To model Single Species Control as currently being implemented in the EU (Task Force on Multiannual Plans 2014), we fixed target exploitation rates as in Eq.
(3), but now computed the targets F MSY,i = F SSC,i using single-species surplus production models, fitted to the operating model's apparent dynamics for that species when disregarding all others. See SM S6.

Addressing structural instability through conservative management plans
Mathematical analyses suggest that, due to structural instability, the parameters of mHCRs can sensitively respond to small errors in the population-dynamical parameters of management models (SM S7) and that inaccuracies in mHCR parameters can lead to large reductions in yields (SM S8). It might therefore be advantageous to pursue management strategies that are conservative in the sense of preferring targets that differ less from current states than other, non-conservative strategies, even if the latter are predicted to be optimal by a naive evaluation of the management model. We achieve such conservatism through modifications of the formulae for mHCR parameters. Technically, we implement a variant of Lavrentiev regularisation of matrix inversion (Engl et al. 2000), which improves the numerical robustness of the calculation, see SM S9. A regularisation parameter controls the degree of conservatism. Where applicable, we performed the MSE both with that parameter set to zero (no conservatism) and to a fixed reasonable non-zero value, given in SM S11.2 (standard conservatism). To obtain upper bounds on conceivable further improvements of outcomes through conservatism, we evaluated in addition plans for the case where the parameter was chosen for each plan and sample community such that the actual total long-term yield attained was maximised (optimal conservatism).

Theoretical maximum sustainable total yield
As a yardstick to compare the outcomes of management plans against, we computed the theoretical maximum long-term total yield (MSY theo ) achievable from each sampled model community. This was determined using an evolutionary optimisation algorithm (Stafford 2008). Specifically, we applied the method of  as implemented in NLopt (http://ab-initio.mit.edu/wiki/index.php/NLopt), modified to improve convergence (SM S2.2).

Comparison of management plans
We compared all possible management plans, i.e., combinations of objectives and strategies, based on the total yield they achieved and the resulting impacts on biodiversity and community size structure. The long-term yield resulting from applying a plan to a sample community was computed as the time-averaged yield over the duration of 8 adaptive management cycles, after a 24 cycle transition phase. Because MSY theo varied substantially among sample communities (CV = 0.40), we quantified the yield from each plan as the average over all sample communities of the percentage of MSY theo attained.
As secondary criteria for comparison, we computed the proportion of fish species conserved (100% − proportion extirpated) and the mean logarithmic maturation body mass of species at the end of each MSE run. The latter addresses the specific concern that maximisation of yields might requires culling large piscivorous species to release pressure on smaller and more productive planktivores. As for yields, both indices were averaged over all sample communities. For comparison, we calculated these indices also for unfished communities after simulating them over a period corresponding to 32 management cycles, because slow residual dynamics of stocks (Fung et al. 2013) can result in natural extirpations.
In the policy context, other criteria will also play a role, for example the similarity of plans to those currently in place or, based on game theoretical considerations, the likelihood with which stakeholders will agree to plans or small modifications of them. Rather than attempting to score plans based on such criteria, key considerations are suggested in Discussion. Table 1 displays for each type of plan the average proportion of MSY theo achieved in the MSE, together with the proportion of surviving species and changes in size structure. Typical standard errors are around 3% for yield scores, 2% for survival rates, and 0.12 for mean log size. More detailed statistical results, including error estimates and statistical significance of differences, are documented in SM S11. There we also report yields for a second, statistically independent MSE based on model communities of half the size used here. Remarkably, differences in yield scores between the two MSE are small, suggesting that our findings depend little on the size of communities. Table 1 shows that exploitation leading to MSY theo would come at the cost of conserving only about 72% of fish species and a change in species size structure corresponding a drop of 52% (= 1 − 10 1.65−1.97 ) in geometric mean maturation body mass. Such impacts are plausible in view of a corresponding survival rate of 51% obtained by Matsuda and Abrams (2006). Yields from Single Species Control were only about half of MSY theo and resulted in survival of 63% of species and a decline in geometric mean size by 45%.

Results
Among the full range of plans considered, a number of general trends can be observed: 1. Independent of objective, State Target Control gave higher long-term yields than Pressure Target Control. This is expected from structural instability (Rossberg 2013, Chapter 24), as explained in SM S8. Remarkably, these yield increases were paralleled by higher species survival rates.
2. For the Total Yield objective, conservatism did substantially increase yields with both Pressure Target Control and State Target Control. The same holds, to a lesser extent, for the Nash Pressure objective. As theoretically expected from structural instability (SM S7 and S9), this was not the case for the Nash State objective. Conservatism increased species survival without exception.
3. With standard conservatism, Self-Optimising Control resulted in yields intermediate between State Target Control and Pressure Target Control, independent of objective.
4. Optimal conservatism led to improvements over standard conservatism by 8.3-16.9% of MSY theo . The pattern of improvements was similar to that for the difference between standard and no conservatism.
5. Independent of objective, yields from Pressure Target Control using standard conservatism improved only slightly over Single Species Control. However, yields from Pressure Target Control with optimal conservatism were notably higher, which indicates potential for some improvements through advances in mathematical methods. Without conservatism, yield score, survival rate and impact on community size structure for Pressure Target Control with Nash Pressure objective were almost identical to Single Species Control.

Biodiversity impacts
In our MSE, strategies that yield more are also beneficial to species survival and community size structure, while for management objectives this relation is reversed. While noteworthy, we caution not to over-interpret these relations, because the management plans could be modified to explicitly incorporate conservation constraints (Matsuda and Abrams 2006). To illustrate possible implications of conservation constraints, Table 2 compares yield, survival rate and community size structure with and without explicit inclusion of such constraints in State Target Control management plans for Nash Pressure and Total Yield objectives. The constraint in this case is to not purposely deplete fish populations to less than 10% of their virgin biomass, see SM S5.2. Interestingly, this modification reduces extirpations and erosion of community size structure substantially, without significantly affecting yields. In view of this robustness of yields, we disregard conservation in the following, acknowledging that conservation while aiming at multispecies MSY requires further study.

Structural Instability
Our MSE demonstrates how structural instability of ecosystems, and so of ecosystem models, can strongly impact management outcomes. Results from our MSE based on a complex food-web model closely follow expectations derived from general mathematical arguments (Rossberg 2013, Chapter 18;Supplementary Material). The mathematical arguments also appear to be robust to many real-world complications not considered in the MSE, as already demonstrated for conservation constraints. Structural instability issues persist when transitioning from stock-to fleet-based management (Rossberg 2012a) or when multiplying yields by market values. Recruitment variability and imperfect control of applied fishing pressures lead to additional uncertainty, which structural instability can enhance (SM S8). The MSE also highlighted options to mitigate these impacts. The most important may be adoption of state-rather than pressure targets in management plans. Another mitigation option is conservatism, i.e. regularization of the inverse problems underlying the computation of multispecies MSY reference points. The mathematical scheme for this proposed in SM S9 can be extended to more detailed management models if required. To determine good degrees of conservatism when direct comparison of long-term management outcomes is not practical, one can compare simulated outcomes using members of a plausible model ensemble (Ianelli et al. 2015), applying conservative reference points derived with another model.
The EU's current approach to implementing its Common Fisheries Policy (CFP, EU 2013), here modelled as Single Species Control, is mathematically similar to aiming at the Nash Pressure objective using Pressure Target Control without systematic conservatism (SM S6). Our analysis shows that this type of management plan is particularly vulnerable to structural instability. Expected yields are therefore comparatively low. These difficulties are illustrated in simulations by Lynam and Mackinson (2014) predicting that a transition to the MSY exploitation rates recommended by the International Council for the Exploration of the Sea (ICES) will lower rather than raise total long-term yields.

Choice of objective
The MSE showed that the consequences of the two major policy choices we considered, i.e. of objective and of strategy to achieve it, are qualitatively independent. They are therefore discussed separately hereafter. In doing so, we shall assume that management models do take structural instability sufficiently into account such that the "standard conservatism" case in Tab. 1 is representative.
Among the three objectives analysed, sensitivity to structural instability increases with increasing expected total yield, i.e. in the expected order Nash State → Nash Pressure → Total Yield objective (theoretical expectations are summarised in SM S10). In view of the high empirical uncertainty in the strengths of multispecies interactions, the Nash State objective might therefore be favoured by stakeholders, despite the comparatively low yields this implies. When combined with State Target Control, it has the added advantages of transparency and simplicity; stakeholders can focus on simply reaching agreement on the targeted ecosystem state. Exploitation rates and corresponding yields required to achieve this state would play a role in negotiations, but not explicitly become part of agreements. Our results suggest this approach would yield about 20-30% higher than Single Species Control.
Despite Nash Pressure and Nash State objective having similar definitions, corre-sponding expected long-term yields clearly differ. The higher yields predicted for Nash Pressure could be reason for stakeholders to favour this objective over Nash State. We therefore caution that we could not find a reason why Nash Pressure should generally yield more. Indeed, SM S10 constructs counter-examples. This pattern should therefore be verified using MSEs that vary our setup. Independent of the higher expected yields, there are two considerations that favour Nash Pressure over Nash State objective. Firstly, Nash Pressure could be more acceptable to stakeholders familiar with Single Species Control because it is conceptually similar. Secondly, the fictitious game underlying the Nash Pressure objective comes closer to the situation of real multi-stakeholder fisheries than that underlying Nash State. Agreeing on this objective or variations of it might therefore be easier.
Improvements in yields when going over from Total Yield to Nash Pressure objective were surprisingly small, if present at all. Even with State Target Control, they were just 16% (standard conservatism). This might be insufficient reason to adopt the Total Yield objective, considering that this might lead to difficult negotiations, as explained in the Introduction. Nevertheless we anticipate that, when either Nash Pressure or Nash State objective is taken as the starting point for negotiations, they will develop towards the Total Yield objective by agreements on modifications of targets that give higher overall yields or incomes.

Choice of strategy
Considering the choice of strategy, out MSE strongly favours State Target Control over Pressure Target Control. This result is expected (SM S10), because structural instability leads to large changes in community structure in response to small changes (or errors) in applied pressures. It is worth noting that, for the same reasons, community structure is expected to be highly sensitive to changes in the physical environment. Within limits, State Target Control counteracts this sensitivity. We warn that MSE that do not use a structurally unstable operating model, e.g. an operating model with non-interacting species, will be unable to reproduce these phenomena! While Pressure Target Control is currently favoured by the EU for most stocks, State Target Control is either legal or de-facto policy, e.g., in the U.S.A., Canada, and New Zealand. Even the EU does not appear to be legally bound to Pressure Target Control. The CFP in its current form requires, according to Article 2, that "management [...] restores and maintains populations of harvested species above levels which can produce the maximum sustainable yield." It does then set deadlines for adjusting exploitation rates accordingly, but these are not the principal objective. Further, the CFP requires formulation of multiannual management plans with "objectives that are consistent with the objectives set out in Article 2 [...]". The plans shall include "quantifiable targets such as fishing mortality rates and/or spawning stock biomass". However, considering the low expected yields from Pressure Target Control, we suggest that multiannual plans setting targets for fishing mortality or exploitation rates may not actually be consistent with the objectives set out in Article 2.
Self-Optimising Control strategies have their own strengths. They are largely insensitive to structural instability, do not require conservatism, and yet achieve long-term yields nearly as large as those from State Target Control. Further, the corresponding mHCRs are independent of productivity rates (r i in Eq. (2)), which lets these strategies automatically adjusts stock sizes to corrected MSY levels if productivities change because of short-term or long-term environmental change.
A specific disadvantage of Self-Optimising Control with Total Yield objective, Transposed Interactions Control, is that the prescribed exploitation rate of a stock depends on the abundances of all stocks on which it has a direct ecological effect: prey, predators, and competitors. This can increase measurement and planning uncertain. (For State Target Control such dependencies enter through neutralising exploitation rates, but there only as corrections to past observations, so uncertainty is lower.) Despite its elegance, Transposed Interactions Control might therefore be preferred only in exceptional cases, e.g. for data-rich single-stakeholder fisheries that have already approached the Total Yield objective and require a scheme to effectively adapt to uncertain environmental change.

Implementation
Most management plans we discussed require for their implementation two models with different skills (Dickey-Collas et al. 2014). The first is what we called the management model. It serves to compute parameters of mHCRs consistent with objectives. Emphasis is on the skill to make long-term MSY projections.
In practice, a second type of model is needed to determine the stock sizes B i required to compute exploitation rates F i from mHCRs and to convert these into catch allowances. Conventional age or size-structured stock-assessment models can be use for this. They have the desired skill to accurately determine current states of stocks based on current and historical data. State Target Control requires in addition that these models can estimate neutralising exploitation rates. To achieve this, they must be capable of forecasting states for given exploitation rates a few years into the future. The key skill is therefore short-term projection.
Strategies to implement these two model types will be different. For long-term projections, interactions with the ecosystem at large and general principles such as conservation and dissipation of energy as it flows through food chains will be important. For short-term projections, optimal use of available data will be key, as with current stock-assessment models.

Conclusions
Because in fact fish stocks do interact, it is not obvious how the frequently invoked MSY policy objective should be translated to real-world fisheries. Answering questions Q1 raised in Introduction, we listed a range of management objectives that would for an isolated stock all reduce to classical MSY. Policy needs to clarify which one is meant! Addressing Q2, we considered a range of harvest control rules and implementation details, and found that management outcomes depend sensitively on their choice.
Answering Q3, the preceding Discussion highlighted trade-offs among management options: the plans expected to generate highest yields require large changes in management practice, accurate data, and/or good collaboration among stakeholders. The EU's current scheme of Single Species Control yields substantially less than most other options (Table 1), which answers Q4.
While there are tradeoffs among all options considered, we propose as a conceivable step towards increased effectiveness, sophistication, and predictability of EU fisheries management a transition to State Target Control with a Nash Pressure objective. In the present study, we evaluated the scheme F = F 0 + (1 − B MSY /B)/T , where F 0 is the neutralising exploitation rate, B is the current stock size, and T is a relaxation time parameter, arbitrarily chosen as T =1yr. If the resulting value for F is negative, it is replaced by zero. When applied to an isolated species, the graph of F vs. B for the corresponding HCR is similar, e.g., to that of the Fishery Decision-Making Framework used by the Canadian government (http://www.dfo-mpo.gc.ca/fm-gp/pechesfisheries/fish-ren-peche/sff-cpd/precaution-eng.htm), except that F is allowed to exceed the predicted value for F MSY . Increases of F beyond F MSY can be necessary because F MSY is uncertain and an uncontrolled stock size larger than the target value B MSY could lead, through multispecies interactions, to unintended and uncontrolled effects on other stocks, a situation that State Target Control aims to avoid.
Variations of this mHCR are conceivable, and many may be similarly effective. For example, following discussion within the Working Group on Multispecies Assessment Methods of ICES (ICES 2014b), we evaluated the modified rule where F SSC is the target exploitation rate of the Single Species Control scheme. While with this modified rule target stock size is approached at the potentially slower rate |F SSC |, the rule has the intriguing property that it would reduce to F = F SSC , i.e. to conventional Single Species Control, if the single-species Schaefer model was accurate. Reality is more complicated and the HCR will deviate from Single Species Control, but these deviations may be small. Applying our MSE to this rule with a Nash Pressure objective, we obtain on average long-term total yield of 75.4% MSY theo (standard conservatism), not significantly different from the result using our original State Target Control rule (p = 0.15, paired t-test). Such results open up the possibility of swiftly transitioning from the current scheme to State Target Control, enabling higher future yields and more predictable ecological outcomes, with minimal changes to currently recommended exploitation rates.
As discussed above, the formulation of EU's CFP is open to both Pressure Target Control and State Target Control. We recommend using of this flexibility to achieve best possible outcomes for society and the marine environment.

Supplementary material
As Supplementary material, we provide a PDF file with mathematical definitions and analyses of management objectives and strategies, details of the methods used to implement these in our MSE, and statistical details of results. These supplementary materials mainly provide two kinds of information. (A) Mathematical formulations and analyses of management objectives and strategies and (B) further details and formal descriptions of the MSE reported in the paper. Sections S1 to S10 incrementally build upon each other. To follow the arguments it is therefore advised to read the material sequentially. Section S11 provides additional technical information on the MSE, and is largely independent of the previous sections.

S1 Management objectives
We provide formal mathematical descriptions of the three major types of objectives considered in the main text, Nash Pressure, Nash State, and Total Yield. Among the possible ways to achieve these objectives, only idealised community steady states with fixed stock biomasses B * i (i = 1, . . . , S F ) shall here be considered. Even when such states are not stable for a given and rigidly fixed set of exploitation rates F i , it will often be possible to stabilise it by letting exploitation rates vary sightly around their long-term averages according to appropriate harvest control rules. Thus, oscillatory states, which are economically undesirable, can be avoided and need not be considered here.
With the long-term yield per unit time interval from species i given by Y i = B * i F i (i = 1, ..., S F ), the corresponding total yield is Y = S F i Y i . To keep the formalism simple, we assume that the relationship between combinations of (averaged) exploitation rates F i and attainable steady states B * i is one-to-one, so that community steady states and the corresponding yields can be parameterized in terms of the corresponding exploitation rates F i . Coverage of ambiguities in the description of states by corresponding pressures would inflate the mathematical formalism without adding much understanding of use for the present work.
Mathematically, the Total Yield objective is attained with a set of exploitation rates As conventional, R ≥0 denotes the set of all non-negative real numbers and R S F ≥0 a set of S F -component vectors of non-negative real numbers. Specifically, F denotes the vector of fishing mortalities F i .
The Nash Pressure objective is attained for a set of fishing mortality rates

S2.1 Model structure and parameterization
The Population-Dynamical Matching Model (PDMM) generates large model communities using an iterative process representing the natural assembly of complex ecological communities. A full description of the model equations, including all parameters and variables, is given in Rossberg (2013, Chapter 22). We describe here only the model's main features. A PDMM model community consists of S = S P +S C species of which S P are producers and S C are consumers. Each species j is characterised by its adult body mass and a 5-component vector of abstract vulnerability traits, and consumers in addition by a 5component vector of abstract foraging traits. The trait vectors locate species in trophic niche space (Nagelkerke and Rossberg 2014;Rossberg et al. 2010). They are thought to be related in an undescribed way to the phenotype of species. The state of the population of each species j is characterised by its population biomass B j alone.
For any conceivable consumer-resource pair, their trophic interaction strength depends on their body mass ratio and the relative position of their traits in trophic niche space. Interaction strengths are large if the producer is in the consumer's preferred size range and if the vulnerability trait vector of the producer is close to the position of the consumer's foraging traits vector. The available trophic niche space is restricted by constraining the maximum length of vulnerability traits vectors.
Starting from a few producer and consumer species, larger communities are assembled by alternately adding new species that differ from a resident species by a small change in adult body mass and traits, and then simulating the population dynamics of the system to a new equilibrium. The species that go extinct in the extended community are removed in each iteration.
Community dynamics is formulated through a set of coupled ordinary differential equations, each representing the dynamics of one species' population. In the most generic form, this system of ODEs can be expressed by Eq. (1) of the main text.
In system (1), the linear growth rates g i (B 1 , . . . , B S ) are non-linear functions of their arguments, including extended type II functional responses, which are nonlinear functions of the biomasses B 1 , . . . , B S . Details are as in Rossberg (2013, Chapter 22), with modifications as given at the end of this section. These type II functional responses complicates the interpretation of model outputs, compared to simpler interaction terms, but have been shown to lead to dynamically stable food-webs with topologies similar to those seen in natural terrestrial and aquatic communities, which are otherwise difficult to build.
As shown by , the size of PDMM communities scales approximately linearly with a parameter X that controls both the probability of direct competition among producers and the maximum length of vulnerability traits vectors. The outcome of, say, doubling X is to halve the probability of direct competition among producers and to double the volume of trophic niche space, which is similar to simply modelling two communities in parallel and then "rewiring" interactions to go criss-cross between the two communities. Results reported in Table 1 of the main text correspond to the choice X = 2, for comparison we also report corresponding results for X = 1.
For X = 2 and with the model definition used by  we found that producers tended to evolve increasingly larger body sizes, which is undesirable for studies of open-ocean communities. We avoided this by going back to a model formulation of , where producer population growth rates are given by with l = 0.5, rather than using the Lotka-Volterra-type model of . To compensate for the resulting slightly lower primary production, we further had to scale up all consumer attack rates by a factor 5/4-see  for the rationale underlying such fine-tuning of attack rates.

S2.2 The evolutionary search algorithm used to maximise total yield
Next we describe the method used to find the theoretically maximal total yield MSY theo for each model community. For the two types of Nash objectives, no attempts were made to find the corresponding solutions in the operating model. Numerical methods to find Nash equilibria seem to be much less developed than algorithms to maximise a single goal function, such as total equilibrium yield Y . Even the latter problem is formidable in our case. From experience testing different approaches, it appears that the search space might best be visualised to have the structure of a mountain landscape such as the Alps. The main maximum is surrounded by several side maxima of similar magnitude, and such local maxima may not be connected to each other by steep ridges. Further off, lower side-maxima can be found. Such goal functions make numerical optimisation difficult. Standard hill-climbing algorithms, for example, almost inevitably get caught in low local maxima corresponding to the Prealps. To find the fishing mortalities achieving the Total Yield objective, Eq. (S1), we therefore used the Improved Stochastic Ranking Evolution Strategy (ISRES) of , as implemented in nlopt v2.3 (http://ab-initio.mit.edu/wiki/index. php/NLopt). The optimisation library nlopt provides implementations of several optimisation algorithms, easily accessible through various programming languages. We describe ISRES here in some detail, because, at one crucial point, we had to modify it.
ISRES is an evolutionary search algorithm that iteratively computes the goal function for a population of parameters combinations x i (1 ≤ i ≤ n, here x i = ln F i and n = S F ), chooses a sub-set of "surviving" good combinations, and then obtains the next generation of the population by randomly mutating the parameters of survivors (there is no recombination step). Together with the parameters x i themselves, the mutation step sizes are inherited and mutated along lineages. Specifically, parameter x i and its mutation step sizes σ i in the next generation are computed from the values x i and σ i of a survivor of the current generation according to the rulê where ξ, ξ i , and ξ i (1 ≤ i ≤ n) are independent random numbers with standard-normal distributions, and the parameters τ and τ are chosen following Schwefel (1993). The last step, Eq. (S6c), with α = 0.2, follows a proposal by Runarsson (2002). It controls the degree of heritability of the "phenotype"σ i of mutation step size along lineages, and so provides for some memory and accumulation of experience with different mutation step sizes as lineages evolve. The problem we encountered with this scheme is that Eq. (S6c) biases mutations of ln σ i towards larger values in the absence of selection pressures. It so undermines the intention of Eq. (S6a) to find the appropriate order of magnitude of σ i by conducting a free search on the ln σ i axis until selection sets in because a good value of σ i is found. Starting with large values of σ i , the rugged goal function of our problem does not provide sufficient cues for good mutation step sizes to counteract this bias, so that step sizes remain large and the algorithm never converges.
We overcame this problem by replacing Eq. (S6c) with a corresponding exponentially declining memory on the log σ i -scale. Specifically, we set For each iteration of the search algorithm, total yield Y was computed from the steady state reached after starting simulations from the virgin community and fishing it at constant exploitation rates F i . To simplify calculations, no attempt was made to stabilise dynamics when oscillatory attractors were reached. Rather, we averaged yield Y in this case over a sufficiently long time interval to obtain the mean long-term yield.

S3 Calibrating the management model
The idea behind our calibration of the management model, i.e., of the multispecies Schaefer model in Eq. (2) of the main text, is to perform first an adiabatic elimination of all non-fish species from the PDMM, and then to choose parameters of the Schaefer model such that it reproduces the dynamics of the PDMM to linear order in the distance from the PDMM's current state.
Technically, this is achieved to a good approximation by first numerically computing the interaction matrix G (Berlow et al. 2004) of the PDMM model community, given by its entries and then dividing it into blocks pertaining to direct interactions within and between the communities of fish (F) and non-fish (N) species. Elimination of the non-fish species is achieved by computing the effective S F × S F interaction matrix for fish The linear growth rates r i entering the Schaefer model are then obtained as The resulting Schaefer model can be written in matrix notation as where • denotes the Hadamard or the entrywise product.

S4 Effects of structural instability of the operating model on the management model
We now sketch some considerations relating structural instability of the operating model to properties of the management model. Mathematically versed readers will note that we gloss over several technical intricacies which could presumably be ironed out through more elaborate technical arguments. For structurally unstable ecological models, the ecological interaction matrix tends to have several eigenvalues λ near zero . Hence, there are S-component vectors e of lengths |e| = 1 such that G e = λe ≈ 0 (S12) because λ ≈ 0. This phenomenon persists even after scaling out differences in the characteristic populations sizes and rates of change for different ecological groups, as achieved by going over to a representation in terms of the competitive overlap matrix , which is why we shall disregard such differences for the present considerations. Dividing G into blocks as in Eq. (S8) and correspondingly writing Eq. (S12) becomes (S14b) Now one needs to distinguish two cases. The first is given by |e F | 2 1 and |e N | 2 ≈ 1, in the second |e F | 2 is of the order of magnitude of one. In the first case, Eq. (S14b) reduces to G NN e N ≈ 0, implying that G NN has eigenvalues near zero. Matrices that have eigenvalues near zero and others further away are said to be ill-conditioned. The inverse of such a matrix, e.g. G −1 NN , depends sensitively on the values of matrix entries. Since G −1 NN enters the effective interaction matrix for fish through Eq. (S9), this can degrade the accuracy of G in describing interactions among fish.
In the second case, elimination of e N from Eq. (S14a) using Eq. (S14b) leads to implying that the effective interaction matrix G itself has eigenvalues near zero. As a result, the Schaefer model Eq. (S11) is then structurally unstable. Thus, structural instability of the operating model can make the management model (S11) an inaccurate description of the interactions among fish or translate into structural instability of the management model itself. For large, complex operating models, as in our case, one can expect both cases to play out to varying degrees. With the more elaborate management models used in practice, this inheritance of structural instability is likely to prevail in similar, albeit less transparent, form.
S5 Solving the management model for MSY objectives and mHCR parameters

S5.1 Standard case
We now derive the parameters entering the mHCR in Section 2.3 of the main text for the objectives that we consider. Setting aside the possibility of extirpations (B * i = 0), equilibrium solutions of Eq. (S11) satisfy We assume that G is non-singular, even though it might be ill-conditioned. Equation (S16) then implies immediately the relations which we shall use in the following. Making use of the continuity of the derivatives in Eq. (S17), one obtains, using standard arguments, necessary conditions for satisfaction of the objectives described in Section S1. For the Total Yield objective these are, with 1 ≤ i ≤ S F , (with T denoting matrix transposition); for the Nash Pressure objective with diag(M) denoting the diagonal matrix with diagonal entries given by those of M; and for the Nash State objective All three conditions can be re-written as where for the Total Yield objectiveĜ = G, for Nash PressureĜ = diag(G −1 ) −1 , and for Nash StateĜ = diag(G).
The Self-Optimising Control rules are obtained as by simply replacing the equilibrium values B * with the current values B of stock sizes in Eq. (S21). State targets B * = B MSY are obtained by putting Eq. (S21) into Eq. (S16) and solving for B * : Pressure targets F = F MSY are obtained by combining Eqs. (S21) and (S23):

S5.2 Modifications to incorporate conservation constraints
Results in Table 2 of the main text were computed using variants of Eq. (S23) and Eq. (S24) that take conservation constraints into account. These are given by the following algorithm. For each species i, a conservation threshold biomasses B cons,i is first defined as 10% of the species' biomass in the unexploited community (Matsuda and Abrams 2006). If B MSY,i , as computed using Eq. (S23), falls below this threshold for one or more species, the species i with the lowest ratio B MSY,i /B cons,i is identified and its target abundance set to B cons,i . Fixing B i = B cons,i in Eq. (S11) and eliminating species i from the system leads to a new Schaefer (or Lotka-Volterra) model of the same form as Eq. (S11) with one fewer species. For this model Eq. (S23) is re-evaluated. This procedure is repeated until B MSY,i > B cons,i for all remaining species. Equation (S24) is then evaluated with B MSY,i replaced by B cons,i where applicable.

S6 Target fishing mortalities for Single Species Control
To derive a HCR for Single Species Control of the ith fish species, we proceeded essentially along the same route as that leading to Eq. (S24), but now adiabatically eliminating all species but species i to obtain the single-species case of the Schaefer model Eq. (S11).
To obtain the non-linear coefficientg i corresponding to the interaction matrix G in the single-species model observe that G −1 , as defined through Eq. (S9), is the upper left S F × S F block of the inverse of G (Bernstein 2005). Correspondingly,g −1 i is the 1 × 1 block on the diagonal of G −1 corresponding to species i, which this is also the ith diagonal element of G −1 . The coefficientr i for the linear growth rate in Eq. (S25) is obtained following the idea of Eq. (S10), for simplicity assuming equilibrium (g i (...) = 0), asr i =g i B i + F i . The combined Schaefer models so obtained for the S F fish species can be summarised in matrix notation as withG = diag(G −1 ) −1 , which is identical toĜ for the Nash Pressure objective. Managers assuming validity of these S F single-species models will compute MSY target exploitation rates as Even though Eq. (S26) disregards multispecies interactions, adaptive management in which the parameters of Eq. (S26) are iteratively adjusted might eventually converge to a community state B * where indeed, with F = F SSC , the equilibrium condition of Eq. (S26),G is met. The community then satisfies the Nash Pressure objective, as can be verified by noting that Eq. (S28) is logically equivalent to Eq. (S21) withĜ =G =G T . Single Species Control for MSY can therefore be understood as a technically simple route to attaining the Nash Pressure objective.

S7 Effects of structural instability of the management model on mHCR parameters
How do structural instability of the management model and parameter uncertainty in G (derived in Sec. S4), affect the accuracy at which the parameters of the mHCR, as derived in Sec. S5, can be computed? To answer this question, we first have to develop further intuition about the structure of the set of the S F eigenvalues of G. Because G is not generally symmetric, its eigenvalues are not generally real numbers-they are complex numbers, represented by points in the complex plane, given by two coordinates: their real and their imaginary parts. Because the entries of G are real numbers, the distribution of its eigenvalues over the complex plane is symmetric with respect to reflection on the real axis. The eigenvalues of asymmetric random matrices tend to be distributed more or less evenly over limited areas of the complex plane, and one can expect this to be similarly the case for G. For ecological model communities built through an assembly process, one finds that this area is convex and comes close to or touches the origin of the complex plane (which corresponds to the eigenvalue zero), but does not fully cover it, thus marking the onset of structural instability . By the sign convention for G used here, this area is located in the right complex half-plane, where the real parts of eigenvalues are positive.
The formula for B MSY , Eq. (S23), can be rewritten as B MSY = H −1 r/2 with The condition of H, and so the impact of inaccuracies in its entries on the accuracy of the reference points B MSY and F MSY (via Eq. (S24)), depends on G and on the definition of Ĝ in terms of G pertinent to the particular objective. The decisive question is how the operation on G defined by Eq. (S29) transforms the area covered by its eigenvalues in the complex plane. This depends on the particular management objective chosen. One can expect that the condition of H becomes poor if the transformed area covers the origin, i.e. if H has eigenvalues with both negative and positive real parts, and that its condition becomes good if the transformed area gets detached from the origin. A first point to note is that, because the sum of all eigenvalues of a matrix equals the sum of its diagonal elements, the respective arithmetic means are equal as well. The "centre of gravity" of the eigenvalues of H in the complex plane is therefore identical to that of G when the diagonal entries of the two matrices are identical, as is the case for Total Yield and Nash State objective (whereĜ = G andĜ = diag(G), respectively).
In case of the Total Yield objective (Ĝ = G), H is simply the symmetric part of G. The eigenvalues of the symmetric part of a matrix are always spread out at least as widely along the real axis as those of the matrix itself (see, e.g., Bhatia 1997, pp. 73-74, presented more accessibly by Raghavan and Barmish 2006), and generally wider. Because typically the set of eigenvalues of G comes close to the origin, one can so expect that the eigenvalues of H cover the origin, implying that H is ill-conditioned.
For the Nash State objective (Ĝ = diag(G)), Eq. (S29) amounts to reducing all offdiagonal elements of G to half of their values. The effect of this is best understood by first considering the impact on eigenvalues of the opposite transformation: gradually adding off-diagonal elements to a diagonal matrix. A simple perturbation analysis (e.g., Bloch et al. 2012), strictly valid only when the off-diagonal entries of G are small compared to the diagonal entries, already gives the general trends, which can be confirmed numerically for larger off-diagonal entries: When exploitative interaction, where G ij and G ji (i = j) have opposite signs, dominate (as for direct predator-prey interactions), the off-diagonal entries lead to a narrowing of the eigenvalue distribution along the real axis. The general expectation, e.g., with statistically unrelated or positively correlated G ij and G ji , however, is that the distribution will be broadened, see, e.g. Proposition 1 of Raghavan and Barmish (2006), where a situation similar to the uncorrelated case is described. Correspondingly, one can expect that the real parts of the eigenvalues of H are distributed more narrowly than those of G and so the condition of H is improved, except when exploitative interactions dominate.
The case of the Nash Pressure objective (Ĝ = diag(G −1 ) −1 ) is similar to that for Nash State, just that now, in addition, the diagonal elements are modified. Applying the blockwise matrix inversion formula (Bernstein 2005) to compute the first diagonal element of G −1 as a 1 × 1 block, one obtains the first diagonal element of H for Nash Pressure from that for Nash State by adding to it −G 1,2...S F (G 2...S F ,2...S F ) −1 G 2...S F ,1 /2, where we specified matrix blocks by corresponding index ranges. Corresponding formulae apply for the other diagonal elements of H. Sign and magnitude of these additions depends on the strength and nature of interactions between species, and also on the condition of the (S F − 1) × (S F − 1) block inversions. This operation tends to broaden the distribution of eigenvalues of H along the real axis, so generally worsening the condition of H for Nash Pressure compared to that for Nash State.
Concluding, these considerations suggest that, for the computation of B MSY and F MSY reference points, structural instability always has a strong impact under the Total Yield objective, generally a small impact under the Nash State objective, and for Nash Pressure an impact that is typically larger than for Nash State but smaller than for Total Yield.

S8 How yields are affected by inaccuracies in mHCR reference points
Above, we discussed in several steps how structural instability of ecosystems (or operating models simulating these) can lead to inaccurate values for the reference points B MSY and F MSY entering mHCR. To build intuition as to what the implications of these inaccuracies are for the resulting long-term yields, we consider here briefly the hypothetical situation where the actual ecosystem was describable by a Schaefer model of the same structure as Eq.
(2) of the main text, with but with true parameters r and G that generally differ from the "empirical" parameters r and G. If the mHCR are successfully implemented and so the respective targets attained, the resulting total yield from State Target Control will be and from Pressure Target Control Similar formulae apply for the individual yield from each stock. Because in Eq. (S30) G enters directly, small errors in B MSY will only have small effect on yields. In Eq. (S31), small errors in F MSY can have large impacts on yields if G is ill-conditioned and so the system structurally unstable (which we expect). Hence Pressure Target Control is expected to result in larger deviations from the expected yield than State Target Control. From other sources of uncertainty in pressures, such as recruitment variability or insufficient control by managers of fishing mortalities, similar adverse effects of structural instability can be expected.
A more detailed analysis would take account of the error structure in reference points resulting from the proposed calculation methods. This analysis is not performed here, because our numerical MSE provides similar information.

S9 Conservatism through matrix regularisation
The manifestation of structural instability discussed in Section S7 implies that the mathematical problem underlying the computation of multispecies MSY reference points falls into the category of ill-posed problems (Hadamard 1923). In practice, such problems are solved using so-called regularisation methods (Engl et al. 2000). In our case we propose to compute mHCR parameters from the parameters r and G of the management model not using Eq. (S23) directly, but instead as and F MSY =Ĝ T B MSY , where κ ≥ 0 is a parameter of dimensions Time −1 Biomass −1 or Time −1 (Biomass/Area) −1 , depending on how biomass is quantified.
With κ = 0, Eq. (S32) reduces to the simpler Eq. (S23) without conservatism, derived above. The term κI in Eq. (S32) shifts the eigenvalues of H away from zero, so improving the condition of the matrix inversion. This method is known as Lavrentiev regularisation (Engl et al. 2000). With this modification, Eq. (S32) implements a form of conservatism: if κ is very large, Eq. (S32) reduces to B MSY = B, i.e. a recommendation to keep stock sizes at their current levels. Equation (S32) applies conservatism only to those components of B that are aligned with eigenvectors of H corresponding to eigenvalues near zero, i.e. where conservatism is most useful because numerical uncertainty is large. This distinguishes Eq. (S32) from alternative conceivable formulae, e.g. B MSY = (1 + κ) −1 [H −1 r/2 + κB], which is also "conservative" for large κ. Conservatism is expected to be beneficial if H is ill-conditioned, and otherwise to be neutral or detrimental in its effect on the ability of managers to meet objectives.
One might expect that the matrix inversion involved in computingĜ = diag(G −1 ) −1 , required for management towards Nash Pressure objectives, could also be made more robust using regularisation methods, but we did not find a regularisation schemes that would improve outcomes among the schemes we experimented with. Except for this, Self-Optimising Control does not require any matrix inversion, so the above considerations play no role if Self-Optimising Control is employed.

S10 Summary of theoretical expectations
Based on the preceding discussions, one can expect that, among the objectives we considered, conservatism will be most important for achieving Total Yield and least important for Nash State, while for achievement of the Nash Pressure objective conservatism has an intermediate role to play (Section S7). On the other hand, it is clear that too strong conservatism will unnecessarily increase the number of adaptive management cycles required to approach management goals. Among the target-oriented management strategies, State Target Control is expected to be more robust than Pressure Target Control (Section S8). For Self-Optimising Control strategies, conservatism is unnecessary (Section S9).
We note that, based on theoretical considerations alone, it is difficult to formulate general expectations of whether total yields will be higher under the Nash Pressure or the Nash State objective. Consider, for example, the case S F = 2, where, based on above calculations and disregarding model uncertainty, the difference between yields under the two objectives is While for perfectly symmetric competition (G 11 = G 22 , G 12 = G 21 , r 1 = r 2 ) this difference evaluates to 4G 11 G 2 12 r 2 1 (4G 2 11 − G 2 12 ) −2 (G 11 + G 12 ) −1 , which is typically positive, it is clear that the sign of Eq. (S33) will flip, e.g., whenever either that of G 12 or G 21 flips. In general, pursuit of either Nash Pressure or Nash State objective could therefore result in higher total yields.

S11 Detailed quantitative outcomes of the MSE
We now present details of the numerical MSEs we conducted based on two sets of communities, generated by using the PDMM with the scaling parameter set to X = 1 and X = 2, respectively (Section S2.1).  Figure S2: Yield time series when applying Pressure Target Control (red), State Target Control (green) and Self-Optimising Control (blue) strategies for the Total Yield objective on the X = 2 sample communities saved after (a) 6 · 10 4 , (b) 10 5 , (c) 1.5 · 10 5 and (d) 2 · 10 5 species added. Figure S2 shows typical time series of simulated yields. They correspond to applying the Pressure Target Control, State Target Control, and Self-Optimising Control management strategies for the Total Yield objective over a period of 32 management cycles to four of the sample communities (X = 2) that we studied. It is apparent that yields can vary substantially between management strategies and between sample communities, and that even the ranking of yields achieved using the three strategies can vary. In order to establish the relative performance of different management plans we therefore tested them on a larger number of sample communities.

S11.2 Standard regularisation and optimal regularisation
The numerical results for plans involving a "standard" degree of conservatism were obtained using a fixed regularisation parameter κ = 5 g −1 m 2 yr −1 in all cases considered and for all sample communities. To put this value into context, we computed the ratios F MSY /B MSY for Single Species Control for all fish species from all virgin model communities, and found that our choice for κ corresponds approximately to the 2/3 quantile of their distribution. This comparison should be useful for choosing a corresponding regularisation parameter in practical applications. In particular, it is noteworthy that κ is not such a "small" parameter here as the theory of matrix regularisation might suggest.
To compare this choice of κ with outcomes for other choices, we conducted, for each plan and each sample community, a series of 9 MSE, using a set of 9 log-equidistant values of the regularisation parameter ranging from κ = 0.5 g −1 m 2 yr −1 to κ = 50 g −1 m 2 yr −1 . To facilitate the comparison, we normalised for each sample community and each plan the average yields to the maximum value across the 9 runs, and from this computed mean and standard deviations across all sample communities. Thus, the value 1 for the mean (and standard deviation 0) associated with a certain regularisation parameter value would indicate that the use of that value produces always (i.e for all sample communities) the highest yield for the plan under consideration. The calculated means and standard deviations for both X = 1 and X = 2 communities are shown in Figure S3 (a) and (b), respectively.
The results of Figure S3 show that our choice of the standard value κ = 5 g −1 m 2 yr −1 generally performs close to the "best" value (i.e. the value corresponding to the largest means shown in Figure S3). While some modest yield improvements can be achieved in some cases by adjusting κ, it should be observed that the standard deviations shown in Figure S3 remain large even for the points with the highest mean (i.e. the "best" κ), indicating high variability among sample communities in terms of the best-performing κ.
This idiosyncratic behaviour appears to be most pronounced for Pressure Target Control and for Nash-type objectives.

S11.3 Detailed quantitative results
In Figures S4 and S5 we display a detailed quantitative comparison of the outcomes of the MSEs. A key to the codes we used for different plans is given in Table S3. Figure S4 displays graphically and numerically mean and standard deviation of the proportion of theoretical MSY attained, as well as the standard error of the mean. In order to account for a weak correlation between results from subsequently sampled communities, the standard errors of the mean in Figure S4 were computed based on the "effective"   Figure S4: The means (bars), standard deviations and the standard errors of the mean (error bars) for total yields (in proportion to theoretical maximum MSYt), calculated across the all X = 1 (left) X = 2 (right) model communities, when applying Self-Optimising Control (blue), Pressure Target Control (red) and State Target Control (green) strategies, compared with results for Single Species Control (grey).  Figure S5: The p-values from two-sided paired t-tests comparing normalised yields (top), survival rates (centre) and mean log sizes (bottom) for each pairing of 23 plans over the 37 sample communities for X = 2, corresponding to Figure S4 (right). Values larger than 0.05 and 0.1 are highlighted in yellow and red, respectively.
number of degrees of freedom following Pyper and Peterman (1998), rather than the actual sample size (total number of sample communities). Because all plans were evaluated based on the same set of sample communities, so that MSE were not statistically independent, we used two-sided paired t-tests to assess the statistical significance of differences between outcomes of different plans. Outcomes of the t-tests for X = 2 for each pair of plans are shown in Figure S5. For the smaller (X = 1) communities, this comparison had slightly lower statistical power, but otherwise produced very similar results.