Formulation, General Features and Global Calibration of a Bioenergetically-Constrained Fishery Model

Human exploitation of marine resources is profoundly altering marine ecosystems, while climate change is expected to further impact commercially-harvested fish and other species. Although the global fishery is a highly complex system with many unpredictable aspects, the bioenergetic limits on fish production and the response of fishing effort to profit are both relatively tractable, and are sure to play important roles. Here we describe a generalized, coupled biological-economic model of the global marine fishery that represents both of these aspects in a unified framework, the BiOeconomic mArine Trophic Size-spectrum (BOATS) model. BOATS predicts fish production according to size spectra as a function of net primary production and temperature, and dynamically determines harvest spectra from the biomass density and interactive, prognostic fishing effort. Within this framework, the equilibrium fish biomass is determined by the economic forcings of catchability, ex-vessel price and cost per unit effort, while the peak harvest depends on the ecosystem parameters. Comparison of a large ensemble of idealized simulations with observational databases, focusing on historical biomass and peak harvests, allows us to narrow the range of several uncertain ecosystem parameters, rule out most parameter combinations, and select an optimal ensemble of model variants. Compared to the prior distributions, model variants with lower values of the mortality rate, trophic efficiency, and allometric constant agree better with observations. For most acceptable parameter combinations, natural mortality rates are more strongly affected by temperature than growth rates, suggesting different sensitivities of these processes to climate change. These results highlight the utility of adopting large-scale, aggregated data constraints to reduce model parameter uncertainties and to better predict the response of fisheries to human behaviour and climate change.


Introduction
Global oceanic wild fish harvest grew at a tremendous rate over the 20th century, increasing by approximately a factor of four between 1950 and 1990 [1].As a result, biomass has been depleted relative to its pristine state [2,3], altering ecosystem structures worldwide, with collapses documented in between 7 and 25% of fisheries [4,5].Anthropogenic climate change is altering primary production and temperature distributions worldwide [6], with impacts on the ranges of various marine species [7,8], and on their growth, mortality, reproduction, and recruitment rates [9,10].Meanwhile, rapid and ongoing changes in fishing technology [11], the demand for fish products [12], and regulation frameworks are driving changes in the distribution and intensity of fishing effort.Given future projections of increasing human population, and consequent pressure on food resources, it is critical to develop improved predictive understanding of how fisheries will evolve with interacting environmental and human drivers.
Numerical models of fishery-human interactions that can be coupled to representations of the environment, often called end-to-end models, have been extremely helpful to shed light on these issues [13,14].However, most of these models are designed primarily in order to study internal ecosystem dynamics, and often include complex parameterizations that are challenging to quantify, such as feeding relationships [15].In addition, most of these models are designed to represent the present-day community in a specific region, and therefore lack the generality and flexibility required to represent long-term structural changes such as those caused by new invasive species, climate-driven range shifts [7], or fishing-driven evolution [16].The incomplete view provided by a patchwork of regional models also hampers studies of global food security, given the development of globally interwoven fishing fleets and fish-trading since the 1950's.
Alternative ecological models that are structurally simpler, founded on macroecological principles, and appropriate for a global range of conditions, have recently begun to allow the direct prediction of fish biomass from environmental variables [17][18][19][20][21][22][23].However, most of these models were designed to study fish, rather than fishing.As a result, fish harvests are generally prescribed in these models, rather than treating the fishing effort as an interactive component, which limits the consideration of economic factors such as changes in fishing technology, market conditions, and regulatory regimes.Furthermore, an extensive exploration of parameter uncertainty is challenging for most of these models due to complexity and/or computational cost, and the simulated fish species may not be equivalent to those harvested, making direct comparisons with observed historical harvest records difficult.
Here we describe the economic component and general characteristics of the BiOeconomic mArine Trophic Size-spectrum model (BOATS), a computationally-efficient model which directly couples a size spectrum model of fish biomass with an interactive economic model on a spatially-resolved global grid.The ecological model, which is described in detail in [24], focuses on a small number of robust bioenergetic principles of marine ecosystems that are of direct relevance for fisheries and environmental change, predicting the growth of all commercially-harvested marine consumers, which we refer to as 'fish'.The simulated fish spectrum is coupled directly to a prognostic model of fishing effort, at the same grid scale, that evolves over time in response to the local bioeconomically-determined profit.Because BOATS was designed to include only fundamental processes and relatively few parameters, parameter uncertainty can be addressed in a robust and comprehensive way by using a Monte Carlo approach.As discussed below, this allows us to select global parameter combinations that agree with observations at the scale of the Large Marine Ecosystems (LMEs), rather than being locally-tuned.We show that, after optimization, the ensemble has a high degree of skill, making it suitable for studies of first-order interactions between climate, the marine ecosystem and humans.We also point out the most important implications of the remaining parameter uncertainty.

Model description
Overview.The guiding principles during the development of the BiOeconomic mArine Trophic Size-spectrum (BOATS) model were: 1.To maintain general applicability to a broad range of marine ecosystems and fisheries; 2. To include sufficient ecological complexity to represent the fundamental dynamics of commercial fish biomass production, but remain simple enough to be easily interpretable; 3. To minimize the number of uncertain parameters; 4. To include a dynamical representation of fishing activity in the same framework as the ecosystem; 5. To enable long-term projections of the coupled natural-social system with robust and simple external forcings.
Fig 1 gives a schematic representation of the model.Biomass production is determined by the energy supplied from photosynthesis [25], the life history of individuals according to temperature-dependent growth and mortality rates [26,27], and stock-productivity-dependent recruitment [28].Feeding relationships are not explicit, a simplifying assumption motivated by the diversity of marine life and predatory strategies, and the heterogeneity of species interrelationships in different regions.Life history is resolved in terms of growth from juveniles to reproductively-mature individuals, through a continuous spectrum of size.The model is deterministic, discretized on a 2-dimensional spatial grid (representing the vertically-integrated ecosystem), and steps forward in time from a set of initial conditions.The model can include an arbitrary number of size spectra, each of which represents a 'superorganism' that exists in all grid cells at all times.There is no simulation of horizontal movement by fish or fishing boats, so that each cell is independent (that is, no interaction with other cells).Although this ignores the life history aspects of larval transport and migration, such as important environmental contrasts between spawning grounds and adult feeding grounds [29], these were determined to be difficult problems when considering all types of commercial species, and resolving movement would incur a large computational burden [21].Fishing effort is interactively determined in each grid cell, and evolves over time as a function of the harvest in that grid cell, dependent on globally-constant economic forcings.In the following, we refer to globally-constant scalars that determine the behavior of a model but do not vary during or between simulations as 'parameters', and to a model with a particular set of parameters as a 'model variant'.Prescribed quantities that influence the model behaviour and may vary during or between simulations are referred to as 'forcings' (which may be globally homogenous, or locally variable), while model state variables are referred to as 'variables'.Key assumptions are summarized in Table 1.
Below, we elaborate on details of the model.First, we present a brief overview of the ecosystem module formulation and equations, and refer the reader to [24] for a full description.Second, we describe the dynamic determination of fishing effort.Third, we discuss the allocation of effort over the size spectra, representing the selectivity of fishing gear.
Ecological module.BOATS resolves an arbitrary number of size spectra that are defined by their asymptotic sizes and, together, are defined as the sum of all consumers (including both finfish and invertebrates) that were commercially harvested before the year 2006.Here we discretize the commercial species among three spectra, which we refer to as small, medium, and large groups, referring to their asymptotic size.The evolution of biomass in each group spectrum is represented using the McKendrickvon Foerster model [30,31], where m is the individual fish mass, t is time, f k (gwB m −2 g −1 , where gwB are grams of wet biomass) is the biomass spectrum of group k, γ k (g s −1 ) the group growth rate, and Λ k (s −1 ) the group natural mortality rate.Fish biomass growth rates are limited by a production spectrum that estimates the energy (here assumed to be equivalent to biomass) that can be provided to any given size through trophic transfer from primary production ξ P,k (g s −1 ) [18,25,32], and are also limited by allometric estimates of the rates at which individual fish can grow ξ VB,k (g s −1 ) [28,33,34].Underlying the simulated biomass spectra are the calculated flows of biomass energy from photosynthesis to fish via the trophic transfer, which we express as where π (gwB m −2 g −1 s −1 ) is the spectrum of energy from primary production available to each individual of size m that depends on the trophic scaling τ as the exponent of a power law dependence, and ϕ C,k is the fraction of primary production allocated to a commercial group k [18,25,32].Furthermore, we assume that the flow of biomass energy from small to large fish, which occurs as individuals grow according to empirical allometric rules, cannot exceed a maximum physiological rate where b (unitless) is the allometric scaling constant, A (g s −1 ) is the allometric growth rate, and k a (s −1 ) is the mass specific investment in activity [28,33,34].The allometric growth rate A is the product of a temperature-independent rate A 0 (g 1−b s −1 ) and a van't Hoff-Arrhenius (exponential) temperature dependence a A (T) which is set by the growth activation energy of metabolism ω a,A (eV) parameter [25].
The realized input energy to growth and reproduction in each size class, ξ I (g s −1 ), is the minimum of ξ P and ξ VB .This implies that, given an excess of food, individual fish will grow as fast as biologically possible, given their temperature-dependent metabolism.If, instead, the production of food at a given location is insufficient to fully feed the existing biomass, this production will be equally partitioned among all individuals and will limit growth rates.When there is no harvest applied, primary production, and thus ξ P , generally limits growth over the biomass spectrum.
The natural mortality rate Λ k is based on an empirical relationship that depends on the individual and asymptotic mass, and is first order in fish biomass, according to where h is the mortality allometric scaling and m 1,k is the asymptotic mass of group k [26,27].
The component of the mortality rate that is independent of mass, λ, is proportional to e z 1 A 0 a l ðTÞ, where z 1 is the mortality rate parameter and a λ (T) the temperature dependence of mortality, which, as with growth, takes the van't Hoff-Arrhenius functional form and is set by the mortality activation energy of metabolism parameter ω a,λ (eV).
Recruitment forms the bottom boundary condition to Eq (1), and is written as where R e,k is the potential generation of new recruits from egg production, itself a function of 1. the primary production and the fish biomass that is able to reproduce and 2. the fraction of larvae that survive to the recruit size m 0 [34], whereas R P,k is the potential flux of primary production input to juvenile fish.The model is forced with two-dimensional, global, monthly climatological fields of net primary production and temperature (S1 and S2 Figs).A schematic of BOATS, illustrating the main modules, components, and processes, is presented in The determination of fishing effort.Harvest depends upon fish biomass and the fishing effort.Here, we interpret effort as the combined inputs used in the fishing process [35], essentially fuel, labour and the construction and maintenance of fishing fleets, and we express it in energy units of W m −2 (1 W = 1 J s −1 ).This choice facilitates comparison with the global effort database of [1,36].Effort of this type is often described as 'nominal effort', and is defined here as being independent of fishing technology.
The allocation of fishing effort, that is, the dynamics of fleets of fishing vessels, involves biological, social, and economic aspects that depend on numerous variables that operate on a variety of timescales [37,38].In order to provide a first step in prognostic modeling of these interactions in a simple yet robust way, we couple the Gordon-Schaeffer open-access fisheries economics model [39,40] to our ecological model [24] in a discrete, temporally-resolved framework.This open-access model assumes an absence of regulation, so that fishermen tend to individually seek the greatest total harvest [39,41], which ultimately leads to overharvest and to the tragedy of the commons [42,43].Although this is not an economically or biologically optimal outcome, it is common, given the difficulty in establishing property rights over living marine resources.Although management has made great advances in effectively counteracting this tendency in some fisheries, particularly since the 1990s [44,45], the open-access result has historically dominated, and remains a cause of overfishing in many open-ocean and coastal fisheries [46].Subsidies can also shift behaviour from the open-access solution, although wherever subsidies are capacity-building they often exacerbate the open-access problem [47].
In open-access fisheries, fishermen seek to maximize their individual harvests, and so are driven by average instead of marginal profit [39].Effort in a fish group k, E k (t) (W m −2 ), therefore evolves in time as a linear function of the average profit (profit per unit effort), as where κ e (W 2 m −2 $ −1 ) is the fleet dynamics parameter [48] that determines the timescale at which a change in average profit affects effort, and revenue and cost are in units of $ m −2 s −1 .
The fleet dynamics parameter κ e is set so that the adjustment of effort to changes in profit is approximately 10 years, since fleets adjust to changes in biomass on roughly decadal timescales [49][50][51].Although the entry of effort into fisheries often occurs at a faster rate than exit, particularly when subsidies are applied, for simplicity we assume here that both rates are the same.We assume that harvest is linear in both biomass and effort, and so for a given level of harvesting technology, a fixed amount of effort always catches the same fraction of biomass [40].The ecological model discretizes biomass along a size spectrum, and so we define an analagous harvest spectrum.The harvest h k (gwB m −2 s −1 g −1 ) in the mass range dm is that which results from effort E k applied to group k during time interval dt, which we define as where q k (m 2 W −1 s −1 ) is the catchability and σ k (unitless) is the size-dependent selectivity of harvest on group k (discussed further in the following subsection).
The catchability q k of fish in group k is the fraction of selectable biomass harvested for a unit amount of effort [52].The catchability encapsulates the ease with which fish can be harvested given the inherent characteristics of fish, and the level of applied technology, including both embodied and disembodied aspects of technology [11,53].Thus, the adoption of improved vessels and fishing gear, sonar, communications, and navigation technologies, as well as improved skipper knowledge and efficiency practices, would all cause an increase of q k .In contrast, q k does not include technological gains that reduce the cost of fishing per unit effort.Nor does q k include technological gains that might enhance the market value of fish or provide access to new markets, which would increase the ex-vessel price.
The revenue spectrum r k ($ m −2 s −1 g −1 ) is the ex-vessel price p k ($ gwB −1 , the price that fishermen are paid at port) of the harvest multiplied by the harvest spectrum h k of Eq (7), and so r k dtdm gives the revenue associated with the harvest of group k for fish of mass range dm at time t during dt, Studies have shown that when numerous species are considered together, a scale that conceptually matches our approach, the price is not strongly related to the trophic level [54].Further, given the unpredictability of fish prices as a function of size, which relates to complex factors such as transportation, processing and societal preferences, we presently assume that the exvessel price is globally constant.The revenue component of the change in effort expressed in Eq ( 6) is therefore calculated by integrating over mass, so that Note that, although the harvest is resolved as a size spectrum, there is only a single scalar value of effort for a given fish group k at a given location and time, since it is assumed that the fishing effort does not discriminate among size within each group, other than as determined by the gear type through the selectivity function.
We assume that cost is proportional to effort [37], such that where C k is the cost ($ m −2 s −1 ) associated with the effort exerted on fishing group k, and c k ($ W −1 s −1 ) is the cost per unit effort.We estimate the global cost per unit effort forcing by assuming that the real world was close to an open-access equilibrium between 1990 and 2006, so that total global revenue was equal to the total global cost, which is a reasonable first-order assumption based on the fishing cost database of Lam et al., [55] (see Table 5 of that article).Thus, we divide the global revenue from the SAUP harvest database by the global effort as reported in [1].We presently assume that cost has no spatial dependence, for simplicity.Although this ignores the potentially-important role of the transit distance between fishing grounds and ports, it avoids a great deal of complexity, and could be revisited in future work.
Linking effort to harvest: harvest selectivity.The size selectivity curve σ k (m) plays a fundamental role in determining the size spectrum of fish biomass and harvest [56][57][58][59].A variety of functional forms have been used to describe the selectivity, all of which avoid the smallest sizes, and can be generalized as either dome-shaped or sigmoidal.Dome-shaped curves are used for selective gears such as gillnets, driftnets, and longlines, which are designed to avoid larger organisms that belong to non-targeted species.Sigmoidal curves are used for towed gears such as trawls and seines, as well as traps and dredges, that do not discriminate against large organisms [56].Given that between 65% [55] and 75% [60,61] of global harvest is associated with sigmoidal selectivity gear types (S2 Table ), and because the modeled effort is by construction perfectly selective by group, perfectly avoiding the harvest of larger organism types, a sigmoidal form is most appropriate.We adapt the formulation of [34], where m Θ,k is the threshold mass of fish in group k, and c σ is the selectivity slope, which we assume takes the same constant value over all groups (S3 Fig) .We interpret the selectivity to be sigmoidal in length due to the physical constraints of nets [57], and therefore convert to mass using the mass-length relationship m ¼ d 1 l d 2 , where m and l are fish mass and length, respectively, and δ 1 and δ 2 are constants [62].
The threshold mass of fish for a group is defined in terms of the maturity mass . By means of the parameter d m Θ ,k , we adjust the threshold mass of the large group to 25% of the maturity mass, that of the middle group to 50% of the maturity mass, and leave the threshold mass of the small group equal to the maturity mass.This could be thought of as representing proportionally greater bycatch of the juveniles of larger fish species.The parameter e m Θ ,k is used solely to test uncertainty in the threshold mass.Table 2 lists the economic model forcings and variables.Variables Eq (10) [37] doi:10.1371/journal.pone.0169763.t002

Assessment of parameter uncertainty and calibration
In order to evaluate uncertain model parameter values, we compare model simulations with two primary sources of global data.The first is stock assessment data as compiled in the RAM Legacy database [63].These data provide an estimate of fish biomass, which is extremely useful; however, for most ecosystems, the stock assessments do not include the entirety of harvested species, limiting their comparability to BOATS.Therefore, we also make use of global catch data, described in greater detail below, as reconstructed by the Sea Around Us Project (SAUP) [60,61].
Because the harvest depends on economic as well as ecological forcings, there are multiple economic and ecological combinations that can produce the same harvest in a given region.On its own, this could lead to an intractable problem, with too many degrees of freedom.However, we can take advantage of the fact that there is a maximum rate of fish production that an ecosystem can achieve at any given location.Because of the historical increase of fishing intensity, most ecosystems are presently either close to this peak value (fully exploited), or have passed it and gone into decline (over-exploited) or collapse [64].We identify the observed maximum multi-year harvest for each LME (which we denote as l H obs LMEmax , where l = 1, . .., 66 is an index for the LMEs) as an estimate of the peak fish that can be produced by the ecosystem when subjected to a transient increase of fishing intensity.Although most of these estimates will individually include uncertainties on the order of a factor of two, arising from uncertainty in catch reconstructions [65], the large number of LMEs and their global distribution should provide sufficient sampling to compensate these errors.
To evaluate the full uncertainty range of model parameters, including interactions between the parameters, we used a Monte Carlo approach [66].This involved first determining a probability distribution for each parameter of interest, and then randomly generating over 10,000 values for each parameter from each distribution (Fig 2).In order to compare with the observed harvest peaks, l H obs LMEmax , we ran a transient simulation for each parameter combination, gradually increasing catchability over 200 years.In the following, we describe the simulations, the data, and the comparison between them.
Transient simulations with increasing catchability.The parameters considered in the Monte Carlo simulations are presented in Table 3, and further detailed in S1 Table.We consider the parameters that have the strongest influence on model biomass, as detailed in the full description of the ecological module [24].For each parameter, we reviewed the literature for statistically-estimated ranges of values and their probability distributions.When such a statistical range was not available, we used the literature to inform the parameter choices assuming a uniform probability distribution.
For each parameter set, we first spin up the model with a 100-year simulation at a low and constant level of catchability (1 × 10 −5 m 2 W −1 s −1 ), which is designed to result in negligible harvests that are consistent with that set of parameters.To span an appropriate range of catchability values, following the first 100 years of constant catchability, we increase catchability by 5% per year for 200 years.This rate of increase ensures that each simulation will produce a peak harvest under all parameter combinations, and is broadly consistent with estimates of technological increase [53].
Observational data.Under intense harvest, the harvest to biomass ratio reflects the ability of a standing stock to produce harvestable biomass each year.We obtained harvest and biomass data from the RAM Legacy stock assessment database [63], and used them to calculate harvest to biomass ratios (H:B hereafter).Stock assessments, by their nature, target individual species in particular geographic regions.Because BOATS represents all commercial species, its results can only be compared with stock assessments when they represent a significant fraction  of the total ecosystem.Analysis of the RAM Legacy stock assessment database demonstrated that there were a total of 8 LMEs where at least 40% of the of the harvest from the SAUP database were represented: Baltic Sea, Barents Sea, Patagonian Shelf, Benguela Current, North Sea, Okhotsk Sea, Gulf of Mexico, and East Bering Sea.For each of these LMEs, we calculated average harvest and biomass from the years that were common to all of the stocks represented.We calculate peak harvest using the Sea Around Us Project (SAUP) database [60,61] aggregated over functional groups at the scale of Large Marine Ecosystems (LMEs).To estimate the peak harvest at each LME, l H obs LMEmax , we take the mean harvest of the 10 years with largest harvest.To compare observed LME peak harvests to those simulated by BOATS, i;l H sim LMEmax where i is an index for the Monte Carlo simulations, we convert both to area-specific LME harvests [67], and ignore a number of high latitude Large Marine Ecosystems (LMEs) due to well-established biases in chlorophyll estimates [68,69] which carry through to estimates of net primary production.We therefore neglect the following LMEs: Antarctica, Kara Sea, Laptev Sea, East Siberian Sea, Canadian Eastern Arctic-West Greenland, Hudson Bay Complex, Beaufort Sea, Canadian High Arctic-North Greenland, Central Arctic, and Northern Bering-Chukchi Seas.Since the Black Sea is an inland sea ecosystem, and not comparable to the other LMEs, we ignore it as well [67].The harvest in these 11 LMEs is a negligible fraction (< 1%) of the total harvest over all LMEs, and so ignoring them has a negligible impact on our analysis.S3 Table lists all of the LMEs.
The fish biomass modeled in BOATS includes three groups defined by asymptotic size: small (S), medium (M), and large (L) groups (S4 Table ).One of our aims is to compare this simple size-based community structure with the SAUP observations.Although most of the SAUP functional groups clearly denote size, we define an association when this is not the case (eg.large shrimp).For cases where the size group is vague (eg.cephalopods) or there are mixed sizes in a functional group (eg.Small to medium flatfishes) we denote these as other (O).We assume the other (O) fraction of harvest follows the same size distribution as the identified groups, such that where G represents a size group (S, M, or L) and the subscript a represents the adjusted harvest.This adjustment is usually small since the other (O) fraction of harvest is generally small.Analysis of Monte Carlo suite of simulations.The Monte Carlo simulations (Table 3) reveal a wide range of biomass and harvest results.It is notable that, despite the 4-order of magnitude range in the i H sim Globalmax across the full Monte Carlo suite, the mean of the i H sim Globalmax is less than a factor of 2 different than the reconstructed value of 64 Mt y −1 [60,61], suggesting that our priors were not unreasonable.
We then select 100 optimal parameter combinations, representing approximately 1% of the total, by discarding any that violate acceptable ranges among four types of constraints: 1. Global harvest integrated over all LMEs for the year of maximum harvest, H obs Globalmax .There are significant uncertainties in the results reported in the SAUP harvest database [70], particularly in the corrections for illegal, unreported and under-reported harvest.We therefore consider a generous estimate of the uncertainty, allowing the i H sim Globalmax to range between 70 and 150 × 10 12 gwB y −1 .We find that 9.1% of the simulations satisfy this constraint, while 84.5% and 6.4% have too small and too large a peak harvest, respectively.
2. Size structure of catch.Founded on analysis of the SAUP database, we select simulations where the i H sim Globalmax for the medium size group is at least 30% that of the small group harvest, and those where the large harvest is at least 10% but less than 80% of the small harvest.This guarantees that the modeled harvest is not unrealistically dominated by a single group.We find that 29.1% of the simulations satisfy the group constraints, which, when taken together with the total harvest constraint, leaves us with 6.1% of the total simulations.

The harvest to biomass (H:B
) ratio of well-assessed ecosystems.In our analysis of the RAM legacy stock assessment database, we found a wide range of H:B among the 8 analyzed LMEs, but in all cases, it was less than 0.4 y −1 .We apply this constraint by calculating the modeled H:B at peak harvest ( i;l H sim LMEmax ) at each of the 8 LMEs.74.2% of the Monte Carlo simulations satisfy this H:B constraint when applied to the 8 LMEs.Taken together with the prior two constraints, we arrive at 3.9% of the simulations.Although there is great potential to use improved stock assessment data in the future, for the present work we only apply the upper limit on the feasible H:B of these 8 LMEs.
4. The relative productivity of LMEs.We rank the remaining simulations by calculating the correlation between the modeled peak harvest of each LME ( i;l H sim LMEmax ) and the reconstructed peak LME harvest from the SAUP ( l H obs LMEmax ).Because the LMEs span a wide range of NPP and temperature, analyzing the simulated distribution of harvest among the LMEs provides a powerful means by which to evaluate model sensitivity to these two critical input variables.We consider the Pearson product-moment correlation coefficient r (plotted as r 2 ), of the harvest per unit area in each LME, but also the log10 transformed harvest per unit area [67] (not shown), as well as Spearman's rank correlation r s .We find a wide range of correlation values from among the members of the Monte Carlo suite, but there are clear patterns of improving correlation with increasing biomass within the range of acceptable harvest which is consistent over the three correlation techniques considered (Fig 3).From among the simulations with acceptable total harvest, the r 2 ranges from close to zero to 0.59, whereas Spearman's r s ranges from −0.42 to 0.76.We select the best-correlated 1% of the total simulations, corresponding to an r 2 larger than 0.45.Table 3 summarizes the results of the analysis, and presents the mean of the optimized (column Mean Opt.) and non-optimized (column Mean N.O.) parameter values, as well as their standard deviations (columns SD Opt. and SD N.O.).Fig 5 shows the transient time series of total LME fish biomass, harvest, and effort for the optimized set.Total LME harvest steadily increases until it reaches a peak, beyond which a sufficient portion of the global ocean becomes recruitment-limited and harvest begins to fall.All simulations show a monotonic decrease of biomass as catchability increases, since the selectable biomass is progressively reduced towards a limit, the critical biomass, that is discussed in the next section.This occurs in an increasing number of grid cells as greater catchability allows fishing to become profitable.Effort follows the same general pattern as harvest, but consistently lags it.The timescale of biomass and effort change is set mostly by the timescale of technology change, yet effort also responds with the timescale of adjustment to revenue changes (the fleet dynamics parameter κ e ) which delays the change of effort relative to biomass and harvest changes (Eq (13)).We overlay the optimized ensemble of simulated harvests in Fig 5B with two versions of the reconstructed SAUP harvest (black line [60,61], red line [71]), assuming that year 1950 of the SAUP harvest corresponds to year 100 of the 300-year transient simulation (when catchability begins to increase).Both reconstructions show the same general trend as the simulated transient simulations over the historical period, and the red line (which attempts to correct for illegal and unreported catch) is close to the median of the simulations.Although the highest simulated values are possibly over-estimates of the global potential fish production, it is useful to include these variants in the ensemble given that their parameter combinations may represent aspects of the ecosystems realistically, even though they allow too much harvest.

The critical fish biomass
In considering the fundamental dynamics that arise from the BOATS formulation, we arrive at a useful observation regarding the equilibrium solution of Eq (6) for the evolution of effort.We begin with the equations for revenue and cost, Eqs ( 9) and (10), respectively, and apply them to Eq (6) to find that Further assuming that price does not depend on mass, we can write the equilibrium solution of Eq (13) as This states that the selectable biomass (the quantity on the left hand side) of group k is, at steady state, determined solely by the ratio of economic parameters, c k /p k q k .This ratio can also be interpreted as a critical threshold, which we refer to as F crit k , which determines the biomass density that is required at any given grid point in order for harvest to be profitable.If the equilibrium biomass ignoring harvest is less than F crit k , then harvest cannot be profitable and there will be no harvest at equilibrium.In this case, the biomass will be controlled by primary production and temperature, by means of Eq (1).Alternatively, if the selectable group biomass exceeds F crit k , harvest will occur, driving the selectable biomass to asymptotically approach F crit k , while the unselectable portion of the biomass will continue to be driven by primary production and temperature.
The exclusive dependence of F crit k on economic parameters is an important consequence of the open access assumption used in the present version of the BOATS model.Although it would not be expected to hold true in the real world, given for example the movement of fish, multi-species fisheries, and the lack of steady-state, it is a useful simplifying concept for explaining the first-order tendencies of fish biomass in the absence of effective management.

General characteristics of the optimized bio-economic model
In order to illustrate the basic characteristics of the coupled model, we use the parameter combination that gives the highest correlation (r 2 ) between modeled and observed harvests at the scale of LMEs.While we initially focus on results with this single set of model parameters for clarity, it is preferable to use an ensemble of optimal models, in order to better span the uncertainty.
The predictive ability of this model variant can be assessed by comparing the correlation of the SAUP harvest with the primary production used to force the model (Fig 6A ), the simulated unharvested LME biomass by LME (Fig 6B ), and the modeled harvest (Fig 6C).We find that NPP alone explains 32% of the variability in the SAUP harvest, similar to prior works [72], while unharvested model biomass explains 54%, and modeled harvest explains 59%.All correlations have a p-value<0.01.Thus, the inclusion of the interactive economics provides an improvement to the prediction of harvest, even though it did not introduce additional parameters into the optimization procedure.
Fig 7 illustrates some aspects of the model by focusing on a single model grid point in the East Bering Sea LME (64˚N, 165˚W).A detailed description of the ecological model at this site is provided in [24], and showed that the slope of modeled biomass spectra is consistent with observations for a wide range of temperature and NPP values.Here, harvest alters the slope and intercept of biomass spectra (Fig 7A ), truncating the large ends of the spectra, consistent with other modeling studies [73].The truncation occurs because larger fish must grow from smaller fish, and so larger fish are affected both by direct harvesting and the reduction in   growth from smaller fish.Moreover, the loss of large, spawning individuals results in a reduced production of eggs and so a reduction in recruitment (Eq (5)) and juvenile biomass.Harvests peak near the threshold size m Θ (Fig 7B ), and decline both at smaller sizes (because of sharply decreasing selectivities), and at larger sizes (because of rapidly decreasing biomasses and growth rates).
It is worth noting that, at low harvest and high biomass, the input energy from the transfer of primary production to fish ξ P (Eq (2)) is partitioned across a larger number of individual, and is generally the limiting factor for growth.On the other hand, at high harvest and low biomass, the same energy from primary production is partitioned to a smaller number of individuals, increasing their potential growth rates, to the point at which growth becomes limited by the maximum physiological rate, expressed by ξ VB .Thus, all else being equal, increasing harvest and decreasing biomass determines a transition from productivity-limited to physiologically-limited growth in the model.14)).With increasing catchability or increasing NPP, we eventually surpass the F crit k limit for a profitable fishery, and arrive at positive harvest.For a given level of NPP, increasing catchability initially results in rising harvest, until the peak is reached, after which further increases of catchability deplete the ability of the population to produce biomass due to stock-limitation of recruitment [74].As a result, there is a peak equilibrium harvest for every value of NPP.For high catchability, both harvest and the fraction of pristine biomass approach zero (0.01 gwB m −2 y −1 and 0.001 contours at the right of Fig 8A and 8B, respectively).

Insights on uncertain parameters
We use the two-sample Kolmogorov-Smirnov test (Table 3, column KS p-value) to compare the distributions of the optimized and non-optimized parameters to determine if they are generated from the same distribution ( [75], p. 428-441).For 7 parameters, we find no difference in the distributions that generate the optimized and non-optimized simulations.However, for the remaining 6 parameters, there is a clear statistical difference (that is, small p-values).Prominent differences are in the mortality parameters z 1 (Fig 4C ) and h, that have significantly lower means and standard deviations.The difference between the optimized and non-optimized mortality rate z 1 is particularly strong.This indicates that, in our model, only the lower range of possible mortality rates [26] can reproduce the present-day global peak harvest distributions.The optimized mean value of the allometric scaling parameter, b, which we assume follows a uniform distribution ranging from 0.61 to 0.78, is 0.65, that is, indistinguishable from the 2/3 value that is widely used in the literature [76], but somewhat smaller than the 0.75-0.79value of other studies [25,77].We also find a preferred trophic efficiency of 0.15, statistically higher than that of the sampling distribution, and encouragingly close to the value of 0.15 that is widely applied.The activation energy of metabolism for mortality, ω a,λ , is not statistically different from the sampling distribution, but that for growth ω a,A is, with a statistically significant lower mean.This lends support to the idea that growth and mortality rates are not subject to the same temperature dependence [78], and suggests that mortality may be more sensitive to temperature variations than growth.
The ensemble optimization reveals that acceptable parameter combinations have significant correlations between some model parameters.Some of these correlations, including the two strongest ones, are shown in    suggesting that higher mortalities can be balanced by increased energy transfer to large size classes.Finally, we highlight a positive correlation (r = 0.41, p-value<0.01)between the activation energies of growth and mortality (Fig 9C ), which control the temperature dependence of these rates.This covariation indicates that, in order to maintain harvests in agreement with observations, temperature-driven increases in mortality need to be balanced by concurrent increases in growth rates.As noted above, the temperature dependence of mortality generally exceeds the temperature dependence of growth, with 80% of the optimized ensemble lying above the 1:1 line.Globalmax , to the total biomass over the LMEs) varies by approximately a factor of 5, indicating that between 5 and 25% of the standing stock of biomass is harvested every year at the peak.This range indicates that, within the limits allowed by the observational constraints, similar harvest can be sustained by substantially different levels of biomass.
Motivated by the range in model H:B, we select three members from the optimized ensemble that represent low, medium, and high values of this quantity, to analyze the spatial patterns of harvest and biomass at the catchabilities corresponding to the peak harvest over the LMEs.The low and high variants are selected based on their H:B and their similar global peak harvest (calculated to be 214 and 193 × 10 12 gwB y −1 , respectively, when including the LMEs as well as the high seas), whereas the medium member is the parameter set that provides the best correlation with SAUP data (described above).When integrated over the LMEs, the total harvest for the low, medium, and high variants are 122, 102, and 86 × 10 12 gwB y −1 , respectively, corresponding to total biomasses of 1.6, 1. Harvest in the low H:B member is dominated by high latitude, high biomass regions with slow fish production, whereas harvest in the high H:B member is mostly in low latitude, low biomass regions with rapid fish production.The spatial concentration of harvest also differs dramatically, with over 60% of harvest taking place at locations with high harvest rates (above 3 gwB m −2 y −1 ) in the low H:B member, whereas harvest is much more widely distributed in the high H:B member, with only about 25% of harvest taking place at location with high harvest rates.Because of important differences in the model sensitivity to temperature (Fig 9C ), contrasting projections of climate change impacts in low versus high latitude ecosystems [79], as well as potentially different resilience of low and high biomass conditions to overfishing, it will be important to further constrain the uncertainty in H:B.
The significant variability of H:B among the 100-member ensemble arises in part from the generous ranges of observational uncertainty applied in the parameter optimization.Expanding and strengthening the observational constraints, for example from stock assessments aggregated at the LME scale, will further narrow the uncertainty in model parameters, and reduce the ensemble range of H:B.Future work should therefore strive to find stronger constraints on harvest and biomass, such as by focusing on a subset of the best-studied LMEs.

Conclusion
The BOATS model represents the first-order features of global fish biomass and economic harvest using fundamental concepts of macroecology, life-history, and resource economics in a unified spatially-resolved framework.The economic module of BOATS is directly coupled to the biophysical model by means of the fishing effort, which responds dynamically to net profits at the grid scale.Given the complex nature of real-world ecosystems and of the social and economic processes that determine fishing effort and harvest, the ability of BOATS to explain a significant amount of the global spatial variability in harvest is promising, and makes it suitable for exploratory global studies of fisheries dynamics that include both human and environmental forcings.
A distinctive feature of this work is the use of a Monte Carlo approach to parse the uncertainty in the model parameters and to determine which parameter sets generate the best representations of present-day harvest from the Sea Around Us Project database.We find that although most parameter combinations among the probability distributions of our priors are unable to produce the observed fish harvest, the central values of the distributions are associated with a global harvest that is reasonably close to observations.By discarding unacceptable parameter combinations using a series of catch and biomass criteria, we arrive at an optimal ensemble of about 100 parameter combinations that perform well, while still retaining a substantial degree of parameter uncertainty.These optimal parameter combinations tend to include rates of natural mortality among the lower range of the sampling distribution, a result that deserves further attention given that natural mortality rates are difficult to constrain from observations.We also find that the temperature dependence of mortality is most often greater than the temperature dependence of fish growth, which has important implications in a warming ocean.Realistic global harvests can be achieved by parameter combinations that result in different biomass distributions, and range from solutions dominated by slowly-growing abundant biomass in high latitudes and coastal regions, to solutions dominated by fast-growing low biomasses in low-latitudes and open-ocean waters.These remaining uncertainties allow a potential range of responses of the global fishery to human and climate-driven perturbations, which could be further narrowed down by considering improved biomass-based observational constraints.

Fig 1 ,
and the full set of biological model parameters are presented in S1 Table.

Fig 2 .
Fig 2. Schematic diagram of the Monte Carlo simulation procedure and selection of ensemble members.doi:10.1371/journal.pone.0169763.g002 Fig 3 displays the maximum global harvest (which we denote as i H obs Globalmax ) and biomass integrated over all LMEs.Biomass ranges over 5 orders of magnitude, has a mean and standard deviation of 0.22 and 0.63 Gt wet biomass, respectively, and is strongly positively skewed (that is, there are relatively more randomly-generated parameter combinations that produce low biomass (S4 Fig).The smallest biomass values are mostly attributable to a low trophic scaling (Fig 4A) or a high mortality constant (Fig 4C).Meanwhile, the i H sim Globalmax values range over more than 4 orders of magnitude (the smallest 3% of the runs were not included in Fig 3, to improve readability).The mean and standard deviation of harvest are 40 and 79 Mt y −1 , respectively, and as with biomass, the harvest distribution is strongly positively skewed (S4 Fig).

Fig 3 .Fig 4 .
Fig 3. Pearson r 2 of Monte Carlo simulations.Each point represents the global biomass and harvest ( i H sim Globalmax ) integrated over all LMEs for the year at which the global harvest is at a maximum, while the color represents the Pearson r 2 calculated for the modeled area-specific individual LME harvests ( i;l H sim LMEmax ) compared to the peak area-specific individual LME harvest from the SAUP ( l H obs LMEmax ).The grey line represents a reference global harvest of 100 × 10 12 gwB y −1 .doi:10.1371/journal.pone.0169763.g003

Fig 5 .
Fig 5. Total (A) biomass, (B) harvest, and (C) effort over all LMEs for the optimized Monte Carlo transient simulations.Line colors represent, for each optimized ensemble member i, the Pearson r 2 of the area-specific BOATS harvest at each LME l ( i;l H sim LMEmax ) and the peak harvest of each LME based on the SAUP ( l H obs LMEmax ).Black [60, 61] and red [71] lines in (B) correspond to two versions of the reconstructed SAUP harvest, assuming that the year 1950 of the SAUP harvest corresponds to the time in our transient simulations when catchability begins to increase.doi:10.1371/journal.pone.0169763.g005

Fig 6 .
Fig 6.SAUP reconstructed harvest as described by primary production, modelled unharvested biomass, and modelled harvest at the LME scale.Correlations are for the variant from the Monte Carlo simulation that best represents the SAUP database.The circle color represents the average temperature of the LME, whereas the number represents the LME number (S3 Table).All quantities are plotted in log10 space to improve readability.r 2 is the squared Pearson product-moment correlation coefficient, and r s is Spearman's rank correlation.All 6 correlations have a p-value<0.01.(A) NPP vs. SAUP harvest.(B) Modelled unharvested biomass vs. SAUP harvest.(C) Modelled harvest at peak harvest state vs. SAUP harvest.The black line is the identity line (1:1 line). doi:10.1371/journal.pone.0169763.g006

Fig 7 .
Fig 7. Steady state biomass (A) and harvest (B) size spectra.Black solid, dashed, and dot-dashed curves represent unharvested small, medium, and large groups, respectively, whereas grey curves represent harvested spectra.Thick curves represent total biomass and harvest spectra.Curves are equilibrium solutions of simulations for a site in the East Bering Sea LME (64˚N, 165˚W) forced with annual average net primary production and temperature using a timestep of 15 days.doi:10.1371/journal.pone.0169763.g007

Fig 8
illustrates the impact of NPP (ranging from 500 to 3000 mg C m −2 d −1 ) and catchability (ranging from 10 −5 to 10 −3 m 2 W −1 s −1 ) on biomass and harvest.Each harvest (Fig 8A) and fraction of prisitine biomass remaining (Fig 8B) value is the equilibrium result of a 1000-year simulation using a 15-day timestep and constant forcing of annually-averaged NPP and temperature.For low catchability and low NPP (and hence low unharvested biomass), harvest is zero (0.01 gwB m −2 y −1 contour in Fig 8A) because the fishery is not profitable (Eq (

Fig 9 .
The mortality constant, z 1 , and the allometric scaling of mortality, h, are negatively correlated in the ensemble (Pearson correlation coefficient r = −0.82,p-value<0.01)and compensate each other so that a high mortality constant is balanced by a weak size-dependence (low h) of mortality (Fig 9A).We also find a positive correlation (r = 0.68, p-value<0.01)between the mortality constant and the trophic scaling (Fig 9B)

Fig 8 .
Fig 8. Model sensitivity to NPP and catchability q. (A) Total harvest and (B) Fraction of pristine biomass remaining as functions of NPP and q. doi:10.1371/journal.pone.0169763.g008

Fig 9 .
Fig 9. Harvest to biomass ratio (H:B) at peak harvest in terms of three sets of parameter combinations.Circles are members of the optimized set of model parameters.Circle color represents the H: B (y −1 ) calculated as the total harvest ( i H sim Globalmax ) over the total biomass at the year of global peak harvest integrated over the LMEs, while the circle size is proportional to the peak harvest.Horizontal and vertical axes correspond to parameter values from the Monte Carlo suite.Thick orange, red, and peach circles correspond to the low, medium, and high H:B members that are used in Fig 10.The black line in (C) is the identity line (1:1 line). doi:10.1371/journal.pone.0169763.g009

Fig 9
also highlights the impact of several fundamental model parameters on simulated H: B. Low mortality constants, high allometric growth exponents, and high growth rate constants tend to occur with high H:B (r between z 1 , b, and A 0 , respectively, with H:B are −0.49,0.65, and 0.30, p-value<0.01).Across the 105 members of the optimized Monte Carlo set, the H:B calculated for the global peak harvest (calculated by taking the ratio of total harvest over the LMEs, i H sim 5, and 0.2 × 10 15 gwB.Maps (Fig 10), and cumulative harvest distributions (S5 Fig), of the peak harvest and biomass state demonstrate fundamental differences in spatial patterns between the members.

Fig 10 .
Fig 10.Two-dimensional global maps of harvest and biomass for three model variants with varying global harvest to biomass ratio (H:B).Maps show the annual average of the year of peak harvest.Harvest for low (A), medium (B), and high (C) H:B.Biomass for low (D), medium (E), and high (F) H:B.doi:10.1371/journal.pone.0169763.g010

Table 1 . Key model assumptions. Essential Present version (could be changed) Ecosystem
Catchability is globally constant, but is varied over time doi:10.1371/journal.pone.0169763.t001