#### Abstract

Multiple fisheries have collapsed as a result of overfishing and strong limitations in our knowledge of system conditions and consequential ecological interactions. Fishery managers need to establish harvesting strategies that balance economic benefits against ecological objectives, including avoiding unintended catastrophic consequences. Our results show that classical assumptions for fisheries management can yield severe instabilities in our quantified views of socioecological tradeoffs, making their ability to inform stakeholder preferences questionable. The complex ecological interactions implied by different parameterizations of such systems yield highly complex and nonlinear dynamic properties with multiple distinct basins of attraction. We show that small changes in our deeply uncertain representations of predator-prey systems can fundamentally shift their dynamics and the validity of candidate management strategies for harvest. Insights from this study highlight the importance of ensuring models capture deep uncertainties, as well as a breadth of financial and ecological criteria when searching for robust management options for resilient fisheries.

#### 1. Introduction

One in four fisheries has collapsed in the latter half of the 20^{th} century [1]. In the northwest Atlantic, most Atlantic cod (*Gadus morhua*) stocks collapsed in the early 1990s in a decline that was considered to be sudden and unexpected at the time [2]. The collapse has since been attributed to the decades-long overexploitation of the system at unsustainable levels [3] and changes in ocean climate conditions [4]. Overlooked changes in environmental conditions and system interactions led to the collapse of the sardine (*Sardinops sagax*) and anchovy (*Engraulis encrasicolus*) fisheries in the northern Benguela ecosystem off the coast of Namibia in the 1970s. As both species are energy-rich prey, their collapse culminated in significant population declines for their predators as well [5]. In the Volga River, the construction of dams has interrupted spawning migrations and reduced habitat sizes, resulting in the collapse of inconnu (*Stenodus leucichthys*), beluga (*Huso huso*), Russian sturgeon (*A. gueldenstaedtii*), and herrings (*Clupea harengus*) [6]. These catastrophic events have been attributed to imprudent human action on marine and freshwater ecosystems and to deep uncertainty in system conditions and poorly understood interactions [7, 8]. Deep uncertainty refers to situations where parameters and relationships describing the system can be complex and difficult to estimate from empirical data, and experts cannot agree on probability density functions to describe them or on the relationships themselves [9, 10].

Such is the case for predator-prey theory in the trophic ecology field. The standard theory of predator-prey interactions has largely been based on the Lotka–Volterra equations that describe a system of two differential equations with a simple relation of proportionality between prey consumption and predator growth. The core assumptions of this model (prey growth and trophic function) have been challenged by multiple authors on both empirical and theoretical bases (e.g., [11–13]). Arditi and Ginzburg [14] argued that the trophic function should account for the intricacies of the predation process at the macroscale (i.e., that the predators need to share the prey available to them). This proposition (the “ratio-dependent” trophic function) sparked a strong debate [15–18] that appeared to come to its conclusion with the general agreement that a range of predator interference levels can be found in nature (i.e., a predator-dependent trophic function with case-specific levels of interference, [19]).

In this study, we consider a fishery management problem where a fleet must develop a harvesting strategy that balances profits with the ecological stability of a predator-prey system. We use the predator-dependent predator-prey system of equations, proposed by Arditi and Akçakaya [20], that includes parameter *m* for the level of predator interference. Accounting for predator interaction (*m*) and time adaptive human prey harvesting (), the following discrete-time forms of the predator-prey system equations are defined:where *x* is prey abundance, *y* is predator abundance, *z* is the fraction of harvested prey, and is the time index in years. *b* is the prey growth rate and *K* is its environmental carrying capacity. *α* represents the rate at which the predator encounters the prey, *h* is the time it needs to consume the prey, and *c* is the rate at which it converts the consumed prey to new predators. *d* is the predator death rate and *m* is the level of predator interference. Parameterizing the level of predator interference in this way, allows us to move between the two contested equation forms (predator dependence versus not) and examine the effects of this type of parametric uncertainty on management tradeoffs. Furthermore, this is the first study to use this specific form of predator-prey relationships and pair it with adaptive human harvesting of prey. Environmental stochasticity is included using “process noise” factors and , modeled as coming from a log-normal distribution, [21–23]. Additional information on the model and its parameterization is provided in Section 2; detailed parameter descriptions, units, and default values are listed in Table 1.

Different system parameterizations can have profound implications on the ways the system behaves. Even in systems without human disruption (), the parameter values and population interactions can affect system dynamics in complex ways, changing, for example, the presence and stability of equilibria between the species (Figure 1). Focusing on the zero isoclines (i.e., the black lines designating the prey and predator population levels that result in either of the species having a zero growth rate), we can identify equilibria (at the intersection of the two zero isoclines). For the prey-dependent model, if the nontrivial (coexistence) equilibrium is stable, it is also a global attractor, which is the specific value that the system tends to evolve toward, irrespective of initial conditions (Figure 1(a)). If the equilibrium is unstable then the global attractor of the system is a stable limit cycle (Figure 1(d)) (i.e., a closed trajectory in the phase space with at least one other trajectory spiraling into it, [24]). For the generalized predator-dependent model, when the nontrivial (coexistence) equilibrium is stable, it is also a global attractor (Figure 1(b)). For this model, if the nontrivial equilibrium is unstable, it can either lead to stable limit cycles (Figure 1(c)) or deterministic extinction (Figure 1(e)). When predator interference is high () and there is sufficient carrying capacity (*K*), there are two nontrivial equilibria, only one of which is stable (Figure 1(f)). The instability of the nontrivial equilibrium can therefore lead to sustained population oscillations or extinction, depending on the model parameters and the initial conditions [24]. As a result, the complex ecological interactions implied by each model formulation are nonlinear and yield highly complex dynamical impacts on the ecosystem, leading to distinct basins of attraction [25, 26]. The complex dynamics of these systems can be further seen in the bifurcation diagrams in Figure S1 of the Supplementary Material, for all respective systems presented in Figure 1. The diagrams are plotted with respect to parameter *m* and demonstrate both the fixed points and the periodic orbits possible under different parameterizations. Specifically, for the systems presented in Figures 1(b)–1(d), decreasing values of *m* change the stability of the equilibrium point, resulting in periodic behavior.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Precise estimates of predator interference, predation, and growth and death rates are difficult to estimate from empirical data, especially for nonartificial environments [19, 20]. As a result, the form of functional responses for a large number of species remains unknown [27]. This is concerning, as the modern ecological paradigm highlights the significance of species interactions to the management of populations, communities, and ecosystems [28–30]. The consequences of such deeply uncertain ecological dynamics on the management objectives can therefore be significant, potentially leading to unattained harvesting profits, or worse, unintentional population collapse.

In this study, we quantify and analyze the tradeoffs of managing (harvesting) a two-species fishery governed by a predator-prey relationship and assess how deep uncertainty in population interactions can affect the management tradeoffs. This application expands on the work by [31, 32] and [33] that also sought to identify management tradeoffs resulting from a socioecological system, as well as the implications of deep uncertainty in the system’s parameters. As illustrated in Figure 1, the system investigated in this application exhibits a much richer variety of dynamics, in ways that may potentially alter the topological characteristics of the attained tradeoffs (e.g., if the parameterization shifts in a manner that changes the stability of an equilibrium, Figures 1(b) and 1(c)). Furthermore, the system considered here includes an additional dimension, increasing the complexity of maximizing economic benefits without inadvertently driving either of the species to collapse. The prey can be overconsumed by the predator and harvester; the predator can collapse if there is not enough available prey. Finally, we demonstrate how the concept of multiobjective robustness can be valuable for selecting harvesting policies that avoid triggering potentially catastrophic consequences in harvested fisheries.

This study aims to bridge and complement the two intellectual threads shaping the management of fisheries: economics and biology. As recounted by [34], seminal work by [35, 36] helped ground a biological rationale for the management of fisheries by illuminating the connections between fishing effort, mortality, and dynamics. From the economic perspective, early work by [37, 38] established an operational foundation for the management of fisheries, based on capital theory and investment concepts, that addressed resource conservation by establishing an objective that ought to be pursued by management. Optimal control theory methods [39, 40] complemented these efforts by describing optimal action paths to achieve this objective. In recent decades, progressive discoveries on the importance of complex multispecies relationships, trophic connectivity, and ecosystem interactions have questioned the traditional and, to this day, predominant view of fisheries management: that of a single-species- and single-objective-based control that establishes a “maximum sustainable yield” or an “allowable biological catch” [29, 41, 42]. Work by [43, 44] and others demonstrated that single-species assessments and management controls may indeed produce misleading predictions and destructive impacts on ecosystems with multispecies interactions, and in marine environments, all species have predator-prey relationships with other species in their ecosystem.

On these grounds, ecosystem-based fishery management [29, 45] has been promoted as the new paradigm for fisheries management, advocating for the consideration of multiple species and values. However, explicit incorporation of nonharvesting values has been limited in fisheries management studies. In cases where multispecies interactions are considered [46, 47], the inclusion of nonharvesting values has been hindered by the fact that these commodities are not typically traded in markets. This might lead to the underappreciation or even complete omission of the nonmarket values of environmental and ecosystem goods and services from studies aiming to provide support for socioecological systems management [48]. Such formulations might consequently result in the inappropriate suggestion of strategies that promote resource exploitation and degrade the ecosystem and its provisions [49, 50]. Broadening the set of fishery management objectives to include nonharvesting values could result in fundamentally different management strategies [47, 51, 52]. The latest review of fisheries’ decision-making applications using multiple criteria [53] identified several studies that considered either several species or several objectives (including nonharvesting). However, all of the presented approaches (multiattribute utility, linear and nonlinear goal programming, and weighted goal programming) collapsed these multiple objectives into one, using an a-priori formulation of the preferences of the stakeholders to be included in the model in the form of weights. These weights may not accurately reflect the true stakeholder preferences, especially before exploring the wider set of possible options, and in cases of nonlinear, threshold responses in the objective space [54–56]. Furthermore, specifying specific goals or weights before the search may miss potential solutions that are of interest by unnecessarily limiting the search space [57]. To the best of the authors’ knowledge, there has been one study that applied heuristic global optimization for the identification of management strategies for a fishery without collapsing objectives into one, Mardle et al. [58]. The application used the GENOCOP III genetic algorithm [59], albeit inappropriately, with a search population and number of function evaluations that were too small for the problem at hand. This work aims to expand on the current literature by optimizing state-dependent, adaptive harvesting strategies that explicitly consider a broad range of objectives (including nonharvesting), in a multiobjective optimization framework. We believe this is the first application of stochastic multiobjective control for the identification of robust harvesting strategies for a fishery.

The rest of this manuscript is organized as follows. In Section 2 we first explain the system under study and discuss the presence and stability of equilibria. We then detail how we used multiobjective optimization to identify candidate harvesting strategies for five management objectives, while constraining strategies to avoid predator collapse. Finally, we explain how the robustness of each candidate strategy was assessed using several satisficing criteria. In Section 3, we present the potential values achieved in each objective by each of the candidate management strategies and demonstrate the significant instability of these tradeoffs, when considering uncertainty in parameter values. We then explore how interactions between parameters affect system stability and, consequently, the attainment of management objectives, including avoiding predator collapse. Lastly, we show how alternative preferences in management strategy may affect the system and potentially avoid population collapse in unstable systems.

#### 2. Methods

##### 2.1. System Equilibria and Stability

A generalized predator-dependent predator-prey system of equations has been modified for the purposes of this study to account for human action by means of harvesting:where *z* describes the harvesting effort performed by the fleet. The parameter values describing the system are listed in Table 1 and represent our best current state of knowledge. Since the system in this study represents a stylistic example, these values were not derived from a specific empirical system, but were based on values and ranges that appear in multiple literature sources (e.g., [21–24]). In an unharvested system, for the nontrivial (coexistence) equilibrium to exist, the following equation must hold [24]:

A mathematical proof of how this condition also holds for a system with harvested prey is provided in the Supplementary Material. In an unharvested system, for the nontrivial equilibrium to be stable, the following equation must hold:

Biologically, when (3) holds, predators convert the consumed prey to new predators at a higher rate than the rate of their death *d* and handling time *h* (i.e., their losses in time and energy). When (4) holds, the prey isocline (detailed in the Supplementary Material) decreases as a function of *x*, stabilizing the system [24]. For a system where prey is harvested (i.e., the additional parameter *z* only decreases the prey level), the condition must also hold as the prey isocline can only further decrease as a function of *x*. This is also demonstrated experimentally by our study.

We use the stochastic, closed loop control method direct policy search (DPS) [60], also known as parameterization-simulation-optimization, to identify harvesting policies. This allows for the use of a state-aware control rule that maps the prey population level () to the harvesting effort at the next time step (), instead of optimizing all individual harvesting efforts. The following sections describe this approach in detail, beginning with the management objectives, the policies that were optimized, and the algorithm used.

##### 2.2. Optimization of Harvesting Strategies

The optimization is aimed at determining dynamic harvesting policies that describe how much prey to harvest over time in order to optimize five objectives and meet the specified constraint. The objectives are designed to address financial goals and to ensure that the fish population is maintained at natural levels; this gives rise to tradeoffs between candidate harvesting strategies. We identify a set of “nondominated” solutions, which is comprised of the harvesting strategies that perform better than any other strategy in at least one of the five objectives. The nondominated solutions compose optimal tradeoffs where improvement in any single objective comes at the cost of degraded performance in one or more of the remaining objectives. The objectives and constraint are described in more detail below.

###### 2.2.1. Maximize Harvesting Discounted Profits (Net Present Value)

The net present value of harvesting profits for each realization of environmental stochasticity is given by for *T* years, where *δ* is the discount rate used to convert future benefits to present value, is prey abundance at the *t*th year of the *i*th realization, and is the harvesting effort performed for that prey. The expected harvesting discounted profits are estimated as the average across *N* realizations:where .

###### 2.2.2. Minimize Prey Population Deficit

The prey population deficit for each realization of environmental stochasticity is given by , where *K* is the prey carrying capacity (i.e., the maximum population abundance that can be achieved if the prey is not subjected to predation or harvest). The expected prey population deficit is estimated as the average deficit over all time steps, averaged across *N* realizations:

When policies are re-evaluated or reoptimized in systems with different parameter combinations (as elaborated in Section 2.6), the respective value of *K* is adjusted accordingly, so as to reflect the deficit of population as it relates to the carrying capacity implied by the new set of parameters.

###### 2.2.3. Minimize Longest Duration of Consecutive Low Harvest

Given operational costs, the fishery managers would like to avoid long durations of consecutively low harvest, defined by a minimum harvest limit. The duration of low harvest is given by both holding for all *t* in a realization of *T* years. The expected worst case of consecutive low harvest is defined as the average of the max duration across *N* realizations:where .

###### 2.2.4. Maximize Worst Harvest Instance

Given operational costs, the fishery managers would like to avoid exposure to financial risks. Variance-minimizing strategies have been widely employed in the literature to model behavior under risk; authors have noted, however, that the costs of risks associated with uncertainty often depend on higher order moments, such as skewness and kurtosis (tail events in the distribution) [61]. This was approximated in this analysis by maximizing the worst harvest instance as well as minimizing the variance of harvest in every realization (, explained in the next section). The worst harvest instance is approximated here by the 1st percentile of harvest for every T-year realization. The expected worst harvest instance is calculated as the average 1st percentile across *N* realizations:

###### 2.2.5. Minimize Harvest Variance

More traditionally encountered in the literature is the minimization of the variance of deviations from the expected profits [61, 62]. This objective has been approximated by estimating the variance of the obtained harvest in every *T* year realization and averaging across all *N* realizations:

##### 2.3. Avoid Collapse of Predator Population

Considering the unharvested predator population as a valuable species, the population of which the managers would like to maintain, the optimization is also subject to a constraint:

##### 2.4. Formulation of Harvesting Policy

For the purposes of this study, candidate DPS control rules were used to map the current levels of prey population () to the harvesting effort at the next time step (). The optimization was aimed at identifying the parameters describing the control rules, instead of the harvesting efforts themselves, allowing for state-based feedback control strategies. The control rules were in the form of Gaussian radial basis functions (RBFs), following the formulation by [63]. The optimization problem was formulated as follows:where and *n* is the number of RBFs used in the function mapping; in this study . is the weight of the *i*th RBF and the weights are formulated so as to be positive (i.e., ) and sum to one (i.e., ). and are the center and radius of the *i*th RBF. All . We used two RBFs and one input in this study, resulting in six decision variables that need to be optimized for the control rule mapping current prey population to next harvest. More inputs can be used in the policy formulation, but were omitted from equation (10) for the purposes of simplicity. Note that with the application of these control rules, the harvesting effort will not necessarily be the same in each of the *N* realizations as the harvesting action is informed by the respective level of prey, , which is subject to the environmental stochasticity at each realization. Furthermore, this formulation assumes a perfectly accurate measurement of the prey level at all times, as well as a perfectly accurate execution of the harvesting strategy. These simplifying assumptions have the benefit in this study of demonstrating the severe difficulty of the class of fisheries management problems even in optimistic formulations of available information.

##### 2.5. Optimization Algorithm

The multiobjective formulation was solved for realizations of randomly generated environmental stochasticity, with each realization spanning years. Initial prey and predator population values were assumed to be and . The parameter values for the assumed state of the world (SOW) are listed in Table 1. The default SOW is assumed to have a single stable equilibrium that is a global attractor. The parameters of the formulated policy were optimized using the Borg multiobjective evolutionary algorithm (MOEA). MOEAs are heuristic optimization algorithms that follow an iterative search procedure to adapt and evolve a population of possible solutions toward an optimized set. To do so, MOEAs apply various probabilistic operators for mutation, mating, archiving, and selection [64].

The Borg MOEA has been designed for the optimization of a wide array of many-objective, multimodal problems [65]. The Borg MOEA is a stochastic population-based search algorithm. Multiple diagnostic studies have demonstrated its capacity to perform consistently well across multiple challenging applications [64]. Its success in identifying high quality Pareto-approximate solution sets has been attributed to its use of *ϵ*-dominance archiving, *ϵ*-progress, and its autoadaptive exploitation of multiple search operators based on their success in generating high quality solutions [65]. The numerical precision required for evaluating each objective, *ϵ*, is specified by the decision maker, creating multidimensional *ϵ*-boxes for ranking and archiving solutions as the algorithm completes its search [65]. Additionally, *ϵ*-values provide a control on the resolution of the Pareto-approximate set (e.g., [66]). The standard means of specifying *ϵ*-values is to consider the significance of precision in each objective. The *ϵ*-values for the five objectives were set as follows: : *ϵ* = 5, : *ϵ* = 0.005, : *ϵ* = 1, : *ϵ* = 1, and : *ϵ* = 5. The Borg MOEA was implemented using its default parameter values, as recommended by [65]. Since Borg is a stochastic optimization algorithm, its search is affected by the random seed that initializes the population and impacts its stochastic operators. The optimization problem was hence solved with 20 random seed runs, each of which exploited 3,000 function evaluations.

##### 2.6. Robustness Analysis

As illustrated in Figure 1, the parameters describing the predator-prey system can fundamentally shift its dynamics and, by extension, the harvesting strategy one should choose to employ. Given the limitations of our knowledge of the parameters describing the modeled system [19, 20, 27], decision makers may consider identifying potential policies that continue to perform satisfactorily when operated under a broad range of alternative system characteristics, also referred to as SOWs. Such policies are termed “robust” and can be identified using various metrics available in the decision analysis literature [67–69]. The domain criterion satisficing measure [70] quantifies the fraction of potential SOWs, in which a solution meets a desired performance (e.g., a set of criteria), and has been widely used in the robustness literature. When compared against other satisficing- and regret-based measures, this metric has been found to identify solutions in line with the stakeholders’ performance criteria [67]. We generated 4,000 alternative SOWs using a Latin Hypercube Sample across the ranges of the deeply uncertain parameters listed in Table 1, assuming uniform parameter distributions. The prey population at each SOW was initialized at the respective sampled *K-*simulated unharvested population trajectories of both prey and predator for all SOWs can be found in Figure S3 in the Supplementary Material. The purpose of the exploratory parameter ranges is to encompass the best available nominal estimates, while sampling broadly enough in their feasible ranges to discover consequential impacts. The intent of this approach is to shift focus from predicting system conditions to, instead, discovering scenarios that are consequential to the decision makers [71–73]. Our multivariate satisficing robustness measure quantifies the percentage (%) of the SOWs in which harvest management is possible (i.e., there is no deterministic extinction) that meet the following requirements:(1)Net Present Value (NPV) (2)Prey population deficit (Prey Deficit) (3)Longest duration of low harvest (Low Harvest Duration) (4)Worst harvest instance (Worst Harvest) (5)Harvest variance (6)Duration of predator population collapse (Predator Collapse)

These performance criteria should be elicited by relevant stakeholders in a real world application. In this exploratory work, they are set so as to reflect possible performance levels expected by decision makers managing this system. Criteria 1, 3, 4, and 5 are each met by at least 75% of the solutions in the assumed SOW. Criterion 2 was set based on critical fish biomass levels for sustainable yields, as reported in the literature [74]. Criterion 6 was met by all solutions in the assumed SOW as it was defined as a constraint.

Even though more robust solutions may be discovered if the optimization is conducted over a large ensemble of deeply uncertain SOWs, said robustness might also result in increased losses (regrets) in the assumed SOW as well as in SOWs that are never encountered [31]. Instead, in this study, we aim to establish a baseline of performance in the assumed SOW, representing our best current state of knowledge, and then re-evaluate this performance across a subjective and wide enough range of SOWs. This allows us to both minimize regret in the assumed SOW and also highlight the existence of areas in the parametric space where management plans might fail, despite their highly adaptive and optimistic design, due to their ignorance of a fundamental shift in system dynamics. Lastly, this intentionally broad ensemble of parameter combinations produces SOWs where extinction of one or both species is unavoidable (detailed in the following section). Such SOWs were omitted from the robustness analysis of the candidate solutions, so as to only evaluate them in contexts where the choice of strategy actually matters and affects the outcomes.

#### 3. Results and Discussion

This section is organized as follows, with Figure 2 summarizing the motivation and main findings of each section. In Section 3.1, we first present and discuss the objective values attained in the assumed SOW, presented in Figure 3, highlighting the strong tradeoffs between NPV and Prey Deficit objectives. In Section 3.2, we explore the impacts of deep uncertainty on the inferred tradeoffs for four other potential SOWs. Figures 4 and 5 demonstrate severe instabilities in these objective values when strategies are re-evaluated in different SOWs. Further, Figure 5 shows that when compared with solutions identified assuming perfect knowledge of the SOW, there are significant losses in the populations of the two species, despite the highly adaptive harvesting policies. In Section 3.3, we explore how system stability is generally affected by parametric uncertainty, by looking at the SOWs leading to predator population collapse. In Figure 6, we show that there is a multidimensional threshold surface dividing sustained harvesting and recovery from fishery collapse, where there are no management tradeoffs to be attained. Finally, in Section 3.4, we posit that the concept of multiobjective robustness can be a valuable driver for selecting management strategies that meet necessary system performance and avoid system collapse. We present the robustness values of each candidate management strategy in Figure 7, and the implications of choosing two of them for the two fish populations in Figure 8. Figure 2 summarizes the methodological approach of this study, with the motivation and main findings of each section.

**(a)**

**(b)**

**(c)**

##### 3.1. Fishery Management Tradeoffs in Assumed State of the World

The general challenge considered by fishery management is how to best exploit the system through harvesting while balancing multiple societal objectives with respect to environmental sustainability and profit. Figure 3 presents a parallel axis plot of the objective values achieved by each optimized solution in the assumed SOW (i.e., the set of model parameters describing the system, to which the policies were optimized). The performance on each objective is represented by a vertical axis. The points where each line crosses a vertical axis represent the average performance value for that objective across 100 realizations of environmental stochasticity. An upward shift in one of the vertical axes indicates increased preference in the equivalent objective performance. All lines have been shaded according to their performance on the NPV objective. The plot is oriented such that the ideal solution would be a dark horizontal line crossing the top of each vertical axis. Diagonal, intersecting lines indicate pairwise tradeoffs between two objectives, where improving the performance in one objective is only possible with reduced performance in the other. One should note that the identified solutions carry no stakeholder preference (or weight) towards the objectives but have been instead identified by searching the space of possibilities as widely as possible. Presenting objective performance and tradeoffs in such a format allows and facilitates an a-posteriori elicitation and negotiation of preferences by the decision makers [76]. For example, in systems or decision-making contexts, where species conservation is valued more than harvest profits, one can express such weighting by imposing upward moving limits on the respective axes. This preference elicitation process should be iterative, allowing for stakeholders with diverging preferences to evaluate the performance and appropriateness of the candidate solutions for the system at hand.

As indicated by the intersecting lines and the switch in the color gradient, there appears to be strong tradeoffs between the NPV and Prey Deficit objectives, as well as between the Prey Deficit and the Low Harvest Duration objectives. NPV-maximizing policies do so by harvesting large parts of the available prey population, depleting it from its natural levels by 58.2% on average, across all realizations. Similarly, policies that minimize the Low Harvest Duration to zero (i.e., harvesting at least 5% of the available prey at all times) also severely deplete the prey population. In contrast, policies that minimize the Prey Deficit to 5%, achieve very low values in the NPV objective (indicated by their very light color), as well as in the Low Harvest Duration objective. Looking at the two objectives incorporated to minimize financial risk, Worst Harvest and Harvest Variance, one may note that some solutions ranking poorly with regards to their Worst Harvest perform very well in minimizing Harvest Variance. These solutions also achieve low NPVs (as indicated by their light color), suggesting that the low variance is achieved by simply consistently harvesting very little at every time step. This observation is consistent with past robust optimization studies (e.g., [77, 78]) that have found variance-minimizing objectives penalizing outcomes both above and below the mean (i.e., both higher and lower profits when only lower profits are of concern).

The identified solutions have been optimized under the assumption that the abstracted harvesting agent has access to an accurate reading of the prey population at each timestep (i.e., “error-free information”). Furthermore, the solutions and tradeoffs that arise as a result of this formulation only take into account uncertainty in the form of environmental stochasticity (), which is traditionally assumed to be well characterized. These assumptions are commonly employed in the literature [52, 79, 80] and are highlighted here to emphasize the highly optimistic nature of such optimized strategies, even when they account for standard sources of environmental uncertainty. In addition, the harvesting policies assume no incidental by-catch of the predator. Our intent is to illustrate that even in a highly optimistic harvested predator-prey management context, with a well-informed and highly adaptive harvesting agent, deep uncertainties that are traditionally neglected pose severe challenges.

##### 3.2. Tradeoff Implications of Deep Uncertainty

Figure 3 presents the objective values achieved by each policy under a SOW representing our best knowledge of the system, a predator-dependent system with a stable global attractor, the dynamics of which are presented in Figure 1(b). As previously elaborated, precise estimates of growth and death rates, predation, and interference are often difficult to estimate from empirical data [19, 20]. Given this uncertainty in the parameter estimates, Figure 1(c) presents the dynamics of a predator-dependent system with a now unstable equilibrium, resulting from slight shifts in our best estimates of the system’s parameter values. Figures 1(a)–1(f) illustrate the wide variety of dynamic behavior that can occur as a result of uncertainty in the system parameterization, in absence of any human disturbance (harvesting). Having identified the Pareto-approximate set of solutions for the assumed SOW (Figure 3), Figure 4 presents the same set of optimized solutions with their performance in a three-dimensional objective space. The solutions in the assumed SOW are presented in blue, and, in light red, light orange, light green, and brown, the plot shows how their objective values change when re-evaluated in four other SOWs. The parameter values for the four SOWs presented in Figure 4 are provided in Table 2.

When policies are re-evaluated in a SOW with deterministic extinction (dynamics presented in Figure 1(e)), the Pareto front collapses (brown points). Population collapse is deterministic in this SOW and occurs even without any harvest. Deterministic extinction has been the focus of several studies [13, 24, 81], particularly in the context of the prey- versus ratio-dependent predator-prey theory. Deterministic extinction behavior was routinely observed [24], but could not be described by the classic prey-dependent model. The behavior is explained by predators becoming more efficient (higher *α* and *c* and lower *h* values) and dying out after exhausting the prey.

The policies were also re-evaluated in a parametrically nearby SOW with a global attractor (light orange) and in a parametrically distant SOW with a global attractor (light green), both of which exhibit dynamic behavior akin to that in the assumed SOW (presented in Figure 1(b)). Both re-evaluations examine situations where the decision makers find themselves managing a system exhibiting similar dynamic behavior to that assumed, but having different parameter values from their best estimates. In the nearby SOW (light orange), a subset of candidate harvesting policies cause significant losses in the prey population numbers (increased values in prey population deficit) and predator population collapse (Figure 4). In the distant SOW (light green), the conditions for prey growth are more favorable: higher *K* allows for larger prey population and higher *m* stabilizes the system. In this SOW, no significant losses are exhibited in either of the two populations and the adaptive control policies that had been identified appear to outperform the expected objective values achieved in the assumed SOW. However, this is a stable SOW favorable to higher growth and a decision maker might be inclined to inquire into what objective values could have been achieved having had perfect information about the SOW being managed. In other words, if the optimization was performed, having the accurate reading of the parameters describing the system each time, how would the performance of those solutions differ?

Points in dark orange and dark green in Figure 4 present the objective values achieved by the solutions identified when re-optimizing to the equivalent nearby and distant SOWs. In the nearby stable SOW (light and dark orange), the significant losses in prey population numbers are reduced by the solutions identified through optimization with perfect information. This is more clearly seen in the equivalent parallel axis plot (Figure 5(a)). Here, the solutions identified in the assumed SOW and re-evaluated in the one nearby (presented in grey) are contrasted with those achieved in the nearby SOW, having had perfect information during the optimization (presented in shades of orange). The regrets are significant in the Prey Deficit objective, where the value of the worst-performing policy could have been reduced from 97% to 75% with perfect information about the SOW parameterization. More importantly, many of the policies re-evaluated in this SOW fail to meet the Predator Collapse constraint, with the worst-performing among them failing 29.34% of the time. Conversely, all the policies identified with perfect information about the SOW meet that constraint. Looking at the distant stable SOW (light and dark green in Figure 4), even though the state-aware control policies manage to avoid significant prey losses, the regrets in this case come in the form of lower values in the NPV and Worst Harvest objectives. This is despite the fact that the harvesting agent is highly adaptive and has a perfect reading of the prey population at each timestep.

When solutions are re-evaluated in a nearby SOW without a global attractor (light red in Figure 4), multiple basins of attraction exist and changes in initial conditions can lead to different equilibria (as depicted in Figure 1(f)). In this SOW, the predators are more efficient (higher *α* and *c* values), but also exhibit high interference () which has a stabilizing effect on their population [24]. In such a system, the decision maker harvesting the prey competes with more efficient predators and the solutions re-evaluated in this SOW end up significantly depleting the prey (light red in Figure 4), as well as collapsing the predator population in most instances (grey lines in Figure 5(c)). Even if the system dynamics are less favorable in this SOW, predator population collapse could have been entirely avoided by all solutions if the optimization had perfect information about the SOW parameters (red lines in Figure 5(c)).

##### 3.3. System Stability Implications of Deep Uncertainty

Figures 4 and 5 present how the initial perceptions of tradeoffs in the original SOW would be misleading if the fishery was actually described by the parameterizations of the four other distinct SOWs, selected here for the purpose of demonstration. As knowledge limits may yield broad parameter ranges, it is a decision relevant concern to assess how consequential performance changes may be for the tradeoff solutions across a broader sample of the deeply uncertain parametric space. We therefore have explored the impacts of 4,000 SOWs. Figure 6 presents the sampled SOWs in the *α*, *b*, and *m* parametric space, with all other parameter values kept constant. The color of each point indicates the percentage of the original set of tradeoff solutions (%) that do not lead to inadvertent predator collapse within each sampled SOW. Certain parameter combinations cause deterministic extinction, as populations collapse even in the absence of harvest; these SOWs are indicated by the smaller dark red points in Figure 6. Even though the model used in this study is in a discrete-time form, the continuous-time model (equation (2) in Section 2) can be used to study the underlying dynamics of the system. The necessary condition for a stable global attractor equilibrium, as derived for this generalized system with harvest (equation (4) in Section 2) and applied to the parametric space, is represented by the shaded surface shown in Figure 6.

The underlying dynamics of the system can shift such that the coexistence of the two species is deterministically impossible. In other words, a multidimensional threshold surface exists that separates sustained harvesting and recovery from fishery collapse, where there are no management tradeoffs to be attained. Even though this exploratory analysis sampled the parametric space uniformly, its intent has not been to assign equal probabilities to all outcomes, including the deterministic extinction cases, but rather to explore broadly to discover said threshold surface in the parametric space. With regards to the parametric plausibility of such a region, it could be the side effect of permanent change in one or more critical components in the environment of a species. Such changes may be progressive (e.g., climate change [4] or habitat loss [6]) and lead to permanent changes in the growth or death rates of a species, such that it inevitably moves towards extinction [82]. When collapse is not prescribed by the underlying dynamics, it is evident that it can be avoided by some of the optimized strategies, and, as a result, value preference in the objective space becomes the most decisive consideration. When feasible, we hypothesize that the concept of multiobjective robustness [76, 83], achieved through compromise, can be a valuable driver in selecting a strategy that can avoid the collapse of the two species.

##### 3.4. System Behavior as a Result of Decision Preferences

Exploiting knowledge of the deterministic extinction threshold, we shift focus to reevaluating the candidate harvesting strategies in SOWs, where management tradeoffs exist, and assess their robustness. To do so, we re-evaluate each of the Pareto-approximate solutions presented in Figure 3 in the alternative sampled SOWs and compute their robustness using the domain criterion satisficing measure. This measure includes minimum performance requirements defined by six satisficing criteria for the management problems’ objectives and predator population collapse constraint (detailed in Section 2). We then calculate the percentage of SOWs where each solution meets each of the criteria, as well as the percentage of SOWs where each solution meets all criteria. Extinction of the two species is deterministic in some of the SOWs, and thus the domain criterion satisficing measure was only applied in SOWs in which human action (and choice thereof) matters. Figure 7 presents the percentage of SOWs where each of the optimized solutions meets each of the six criteria (placed on a vertical axis). Each line is shaded according to the percentage of SOWs where it meets all specified criteria. As with Figure 3, diagonal, intersecting lines indicate pairwise tradeoffs in the robustness of the equivalent criteria. For example, there appears to be a strong tradeoff between meeting the NPV criterion in many SOWs and meeting the Prey Deficit criterion in many SOWs. Furthermore, solutions most robust in either of those two criteria (top-most lines crossing the vertical axes) fail to meet at least one of the other specified criteria (indicated by their white color).

Two solutions are highlighted in this figure, and , illustrating two different robustness definitions. The solution meets criterion 1 () in the most SOWs (leftmost vertical axis in Figure 7). NPV-maximizing criteria (or similar metrics of discounted profits) are very common in the bioeconomics literature (e.g., [34, 84–86]), typically considered as the only system objective. The solution meets all of the specified satisficing criteria in the most SOWs, meaning that it exhibits the highest multiobjective robustness (indicated by the line color and color bar on the right of Figure 7). This solution is more in line with newer paradigms in fisheries management, with experts calling for multiobjective perspectives and values [29, 45]. The concept of multiobjective robustness builds on that perspective, by identifying policies that are successful in meeting a specified level of performance in all objectives. The solution most robust in NPV performs poorly in the prey population deficit and the harvest variance criteria. The solution most robust across all criteria achieves robustness by compromising individual robustness for the NPV and low harvest duration criteria, as well as performance in the objective space (see Supplementary Material Figure S2).

With regards to the entire set of optimized solutions and their application in other SOWs, Figure 7 highlights that assumptions on their stability are tenuous even in SOWs where species coexistence and harvesting are possible (i.e., without deterministic extinction). Looking at criterion 2 in particular, the most robust of the solutions in meeting that criterion (topmost line crossing the axis) only does so in fewer than 50% of the considered SOWs without deterministic extinction and, in doing so, it harvests at very low rates and fails to ever meet the other criteria (as indicated by its color). Solution (most robust in criterion 1) harvests at very high rates and fails to meet the prey deficit criterion. Finally, solution (meeting all criteria in the most SOWs) only performs acceptably in fewer than 50% of the sampled SOWs. In other words, out of the original set of optimized solutions, all of which adapting their harvest given a perfect reading of the prey population at each timestep, most fail to meet all the criteria in other SOWs, and the most robust among them only does so in less than half of the sampled SOWs without deterministic extinction. At the same time, when looking at criteria 1 and 2, the equivalent objectives of which are conflicting (maximizing NPV and minimizing prey population deficit), increasing robustness preference for any of the two inevitably reduces robustness in the other.

Linking our observations for this socioecological system to the broader concept of resilience, we draw from its definition as proposed by [87]: *“the capacity of a system to absorb disturbance and reorganize while undergoing change so as to retain essentially the same function, structure, identity, and feedbacks.”* While robustness is a related term, it is defined as the ability of the strategies to maintain an acceptable performance (as measured by some criteria) and it more closely pertains to the decision space of the system, as shaped by the decision-makers’ value preferences. As such, the concept can be used to link complex system dynamics, persistence, and transformation to performance measures (albeit in a more narrow sense-that of feedback and control) [88]. With the application of the concept of multiobjective robustness in this paper, the identified strategy is also more resilient across the sampled SOWs.

This is more explicitly demonstrated in Figure 8, where we present the trajectories of the predator-prey system as a result of the two harvesting strategies in the assumed SOW and two alternative SOWs. We include five trajectories in each of the subplots, aiming to capture imperfect readings of the initial population levels, each starting near the natural (unharvested) equilibrium of the system. Figure 8(a) presents the trajectories of the two populations as they result from the implementation of the two policies in the assumed SOW. Solution appears to harvest at significantly higher levels (indicated by the color of each line segment) and reduce the prey population, resulting in a reduction in the predator population as well. Solution harvests at significantly lower levels and reduces both populations as well, albeit to levels much closer to the natural equilibrium. When applied in a nearby SOW with a global attractor (Figure 8(b)), the harvesting actions behave similarly, shifting the prey isocline, and as a result, the equilibrium moves to lower population levels for the two species (the dynamics of the unharvested system are presented in Figure 1(b)). Figure 8(c) shows the two policies applied in a SOW, where the coexistence equilibrium (trajectory starting point) is not a global attractor. Here, the harvesting can shift the system to a different basin of attraction so as to lead to the collapse of the two populations. This is avoided by the strategy selected for being robust across the multiple objectives ().

#### 4. Conclusions

Progressive discoveries in the ecological literature have highlighted the importance of multispecies relationships and trophic connectivity in fisheries management, resulting in a new proposed paradigm, that of ecosystem-based fishery management [29, 45]. Within this new direction, multiple species and objective values need to be taken into account. This simple exploratory experiment was designed with the specific aim of complementing the current fisheries management literature with the inclusion of a multiobjective perspective for the management of a two-species fishery. The formulation of the system in such a manner allows us to demonstrate the potentially catastrophic consequences of deep uncertainty, as manifested in our parameterized mathematical representations of predator-prey systems that are coupled with human actions (harvesting) and human preferences (multiobjective tradeoffs). For this purpose, the isoclines, equilibria, and conditions for stability were analytically derived for the harvested system and used to discover regions of the deeply uncertain parametric space that lead to system instability and fundamentally impact the estimated tradeoffs. The impacts of deep uncertainty on such a system are significant as distinct basins of attraction can be present and also shift with only marginal changes in assumed mathematical parameterizations. Human action, manifested in the form of harvest, can move the populations into basins of attraction where the attractor is a collapsed system.

As a result, harvesting strategies able to navigate the dynamic complexities of such a system need to be identified. We demonstrate that multiobjective robustness can be valuable as a driver for identifying fishery harvest policies that can navigate deep uncertainties in system parameters and relationships. It is important to note here that the concept of multiobjective robustness is not treated as entirely equivalent to resilience, but rather as an alternative principle by which one can operationalize the identification of such policies in a systematic manner. The principle is applied here with the specific aim of achieving acceptable performance, as measured by explicit criteria, that are defined subjectively and specifically for the system at hand. The broader concept of resilience includes a wider array of aspects and spatial and temporal scales, as well as system boundaries [88], that multiobjective robustness (as applied here) does not claim to consider. Nevertheless, management policies that are not robust fail to meet their decision-relevant criteria and, by extension, lead to a system that is also not resilient. For these reasons, we believe that the findings have broader implications for socioecological management in general, where balancing conflicting economic and ecological objectives is an essential yet complex task that is further impeded by severe uncertainties in determining tipping points and the consequences of crossing them [89].

The system used in the study is specifically formulated at a reduced complexity relative to those found in the mathematical biology literature, where multispecies systems are typical, or in resource economics, where operational costs and more complex financial accounting take place. However, the more complex models in the respective fields are not typically paired with each other and the impacts of deep uncertainty in the ecological system are not manifested through to the decision space, to assess its social implications. The socioecological model employed in this study pairs the two perspectives to demonstrate significant implications for the decision maker even when management is operating in the optimistic context of adaptive, perfectly informed, and perfectly implemented policies. Furthermore, the harvesting policies assume no incidental by-catch nor intentional harvest of the predator, as the aim has been to highlight the implications of deep uncertainty for such a system.

Be that as it may, this simplicity in design bears certain limitations that future work should aim to address. First, the information used by the state-aware control rules is limited to the current levels of the prey fish population. In the interest of avoiding crossing tipping points critical to the recovery and sustainability of a fishery, additional information can be incorporated, for example, predator population levels. This would be particularly valuable for system parameterizations that exhibit multistability (for example, in Figure 1(f)), where remaining with the basin of attraction of the preferred equilibrium is a vital concern. In a real application, this would likely be accompanied by additional costs borne by sensing and measurements. Secondly, the control rules mapping system state to action could take different forms to more realistically represent the harvesting actions performed by fisheries, for example, to avoid significantly shifting the prescribed harvest between sequential time steps, or accidentally harvesting other species (by-catch). Lastly, the study assumes a perfect reading of the prey population at each time step as well as a perfect execution of the prescribed harvesting effort. These are common assumptions in the literature [52, 79, 80], but result in highly optimistic quantifications of the tradeoffs and should be addressed in future applications. Inclusion of a learning procedure about system parameters and dynamics, in the form of a probabilistic approach with sequentially updated estimates [90, 91], would also be a valuable next step.

#### Data Availability

The predator-prey harvesting optimization code, re-evaluation code, robustness analysis, and identified solutions are available on Github at https://github.com/antonia-had/Generalized_fish_game. The optimization and Pareto-sorting can be replicated using the software code available for the Borg MOEA (http://borgmoea.org/), pareto.py (https://github.com/matthewjwoodruff/pareto.py), and the MOEA framework (http://moeaframework.org/).

#### Disclosure

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the funding entities.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Acknowledgments

This study was partially supported by the National Science Foundation (NSF) through the Network for Sustainable Climate Risk Management (SCRiM) under NSF cooperative agreement GEO-1240507 and the Penn State Center for Climate Risk Management.

#### Supplementary Materials

S1 Appendix: derivation of condition (3). S2 Appendix: derivation of the prey isocline. S1 Figure: bifurcation diagrams for systems presented in Figure 1, with regard to predator interference parameter *m*. S2 Figure: population trajectories for prey and predator populations under all sampled SOWs. S3 Figure: parallel axis plot of the objective values achieved by each optimized solution in the assumed SOW.* (Supplementary Materials)*