Following an introduction to the nature of fisheries catches and their information content, a new development of CMSY, a data-limited stock assessment method for fishes and invertebrates, is presented. This new version, CMSY++, overcomes several of the deficiencies of CMSY, which itself improved upon the “Catch-MSY” method published by S. Martell and R. Froese in 2013. The catch-only application of CMSY++ uses a Bayesian implementation of a modified Schaefer model, which also allows the fitting of abundance indices should such information be available. In the absence of historical catch time series and abundance indices, CMSY++ depends strongly on the provision of appropriate and informative priors for plausible ranges of initial and final stock depletion. An Artificial Neural Network (ANN) now assists in selecting objective priors for relative stock size based on patterns in 400 catch time series used for training. Regarding the cross-validation of the ANN predictions, of the 400 real stocks used in the training of ANN, 94% of final relative biomass (B/k) Bayesian (BSM) estimates were within the approximate 95% confidence limits of the respective CMSY++ estimate. Also, the equilibrium catch-biomass relations of the modified Schaefer model are compared with those of alternative surplus-production and age-structured models, suggesting that the latter two can be strongly biased towards underestimating the biomass required to sustain catches at low abundance. Numerous independent applications demonstrate how CMSY++ can incorporate, in addition to the required catch time series, both abundance data and a wide variety of ancillary information. We stress, however, the caveats and pitfalls of naively using the built-in prior options, which should instead be evaluated case-by-case and ideally be replaced by independent prior knowledge.

data limited stock assessments, Elasmobranchii, finfish, global fisheries, informative priors, shellfish, stock status, Teleostei

National and international organizations, notably the Food and Agricultural Organization of the United Nations (FAO), have been tasked with evaluating the global status of fisheries, including countries or regions without age-structured stock assessments. Thus, their staff resorted to developing graphical typologies of annual catch time series, allowing inference on the status of the underlying stocks (e.g., Caddy and Gulland 1983). One of these typologies (see fig. 7 in Csirke 1984) was quite influential and was reprinted in a major textbook (fig. 11 in Hilborn and Walters 1992) and in a review (Kleisner et al. 2013). Additionally, based on these typologies, other FAO staff created stock-status plots (Grainger and Garcia 1996), now a regularly updated part of the bi-annual State of Fisheries and Aquaculture (SOFIA) (Fig. 1). The detailed methodology for creating FAO’s recent stock status plots and for selecting the stocks they represent could be more transparent (Ye 2011); notably, it is unclear how much they depend on the catch trend typologies mentioned above versus “formal” stock assessment outputs, particularly in the Global South, where these assessments are less frequent. However, as shown in Kleisner et al. (2013), their basic trends can be straightforwardly reproduced by the simplified method detailed in Froese and Kesner-Reyes (2002) and Pauly et al. (2008).

Figure 1. 

Common versions of “catch-only” assessments, which may be verified against real stock assessments: A: The key graph in Grainger and Garcia (1996), for which the status of the world’s marine exploited stocks was inferred from the slope of polynomials fitted to catch time series, as subsequently done for years by FAO; B: Reproduction of fig. 19 in SOFIA (FAO 2020), which summarizes the status of global stocks based on an opaque mix of methods ranging from real stock assessments to subjective interpretation of catch trends in data-poor fisheries (see Ye 2011); C: Stock status plots based on the simplification of the catch-only method in A by Froese and Kesner-Reyes (2002), as implemented by Pauly et al. (2008) to work with the catch “reconstructed” globally by the Sea Around Us. Note that C considers stock rebuilding, but that, as a whole, A, B, and C convey the same message.

Although helpful to get a “big picture” overview of fisheries, the empirical catch-only methods behind the stock status plots (Csirke 1984; Grainger and Garcia 1996; Froese and Kesner-Reyes 2002; Froese et al. 2012) are not sufficient for stock-specific status classification and were never meant to assess single stocks for fisheries management. However, especially the approach of Froese and Kesner-Reyes (2002), which scales catches relative to their maximum and evaluates them relative to their occurrence before or after the maximum, has formed the basis of subsequent approaches (Froese et al. 2012; Kleisner et al. 2013) and was used to derive preliminary priors of stock status in the Catch-MSY method of Martell and Froese (2013) and the CMSY method of Froese et al. (2017).

The purpose of this study was to present recent developments and advances in these methods, such as considering the inverse correlation between productivity and carrying capacity for an examined stock, the application of an Artificial Neural Network (ANN) to predict preliminary stock status from time series of catches, and the use of scatterplot catch and biomass data of hundreds of stocks to derive empirical uncertainty ranges for relative biomass priors predicted from relative catch.

The origins of CMSY++. CMSY++ and its predecessors are based on the first derivative of the logistic curve of population growth (Verhulst 1838), with numbers of individuals replaced by the sum of their body weights (Schaefer 1954, 1957)

$B_{t+1}=\left\{B_t+r\left(1-B_t/k\right)-C_t{e}^{{\epsilon }_t}\right\}{e}^{{\eta }_t}$ [1]

where B_t is the biomass and C_t is the catch in tonnes in year t, r (year^–1) is the intrinsic rate of population growth, and k is the carrying capacity of the environment for this population in tonnes, ε_t represents the normally distributed observation error of catches and η_t the process error, respectively, and are implemented as lognormal error terms. Display of these lognormal error terms is omitted in subsequent equations. Thus, if a reasonable estimate of start biomass and k is available to quantify the unexploited and initial stock size, and if a reasonable estimate of r can be inferred from life-history traits (as done in FishBase; Froese and Pauly 2023), a time series of biomass, based on the time series of catches, can be projected, with maximum sustainable fishing mortality F_MSY = r/2 and minimum biomass that can produce maximum sustainable yield (MSY) as B_MSY = k/2. This approach, known as “stochastic reduction analysis” (Kimura and Tagart 1982; Walters et al. 2006), was applied, e.g., by Martell (unpublished*) to a well-studied stock of lingcod (Ophiodon elongatus Girard, 1854) in British Columbia, and by Christensen (2006) to marine mammals globally.

While examining in depth the method that formed the basis of his Bachelor’s thesis (Martell unpublished), Steven Martell found that in a wide-ranging plot of k vs. r, only a small cluster of k and r pairs were “viable”, i.e., did not lead to a population crash nor suggested that a heavily fished population remained close to carrying capacity, and resulted in a final relative biomass within the expected range based on independent information.

Thus, was born the “Catch-MSY” method of Martell and Froese (2013), where in Monte-Carlo fashion, thousands of potential biomass trajectories were filtered according to the above criteria, and “viable” r–k pairs then formed a triangle in log space (see fig. 1c in Martell and Froese 2013). Reasonable estimates of MSY (with confidence limits) could be derived from the geometric mean of viable r and k values inserted in the equation MSY = r · k/4 where an under-estimation of r was compensated for by an over-estimation of k or vice versa (Schiller 2014).

However, as Martell and Froese (2013) explicitly pointed out, the geometric means of viable r and k values diverged systematically from independent estimates based on full assessments, with r being typically underestimated and k being overestimated. Despite the explicit warning to use only MSY and not the reference points F_MSY = r/2 and B_MSY = k/2 for management purposes, subsequent independent tests of the Catch-MSY method (Rosenberg et al. 2014; Chong et al. 2019; Dowling et al. 2019; Zhou et al. 2018; Free et al. 2020; Pons et al. 2020) often gave low scores to the method because they ignored that warning and reproduced the known biases in r and k.

From Catch-MSY to CMSY. Froese et al. (2017) addressed several shortcomings of the Catch-MSY method, mainly by devising an empirical method that identified the most probable r–k pair not at the center, but near the tip of the triangle of viable pairs. Basically, the 75^th percentile of viable r values was used as the best r estimate and a central value of viable k values within the confidence limits of r was used as the best k estimate. This enhanced method, called CMSY, could produce reasonable predictions of relative biomass (B_t/B_MSY) and exploitation rate (F_t/F_MSY). CMSY also included a number of additional improvements, such as a narrower range of potential k values based on maximum catch and the prior r-range (derived from FishBase; Froese and Pauly 2023), and default ranges (i.e., if not supplied by the user) for priors for relative biomass derived from the empirical rules used by Froese and Kesner-Reyes (2002) and Froese et al. (2012) for the construction of stock-status plots. Also, a Bayesian state-space surplus production model (BSM) was added to provide a second set of results based on abundance indices (catch per unit of effort, stock size index, acoustic or trawl survey trends, etc.) if such information was available.

Moreover, CMSY was formulated to account for the generally observed reduction of recruitment at low population sizes (Ricker 1954; Beverton and Holt 1957; Barrowman and Myers 2000; Froese et al. 2017) replacing Equation 1 by

$B_{t+1}=B_t+\left(4B_t/k\right)r\left(1-B_t/k\right)B_t-C_t\mid B_t/k<0.25$ [2]

where (4B_t/k) creates a linear decline of r if biomass falls below k/4, to account for reduced recruitment and thus productivity at low population size. Half of B_MSY, which is k/4 in the Schaefer model context, is usually chosen as the proxy demarcation of the biomass below which recruitment may be impaired (e.g., ICES 2021, p. 3). This empirical feature makes it unnecessary to require additional parameters for a stock-recruitment relation (as in, e.g., Schnute and Richards 2002), which are typically not known in data-limited stocks. This feature was subsequently adopted in other surplus production modelling software, such as JABBA (Winker et al. 2018) and sraplus (Ovando et al. 2021). It accounts for the “weak” Allee effect (Allee 1938; Walters and Kitchell 2001) or the beginning of depensation, as defined by Perälä et al. (2022), but not for the “strong” Allee effect of Hutchings (2015), which describes an even more drastic decline of recruitment and thus resilience and productivity at very low population size. More importantly, by reducing r, it reduces the reference point for maximum sustainable fishing mortality F_MSY = r/2, a major advance over models that consider F_MSY as a constant reference point even when the rate of recruitment, which is a major component of productivity and thus of r and F_MSY, is known to be impaired. Fig. 2 shows a schematic plot of Equation 2, indicating areas of recruitment overfishing (B/k < 0.25), growth overfishing (B/k < 0.5) and sustainable fishing (B/k > 0.5) (Froese and Proelss 2012). While the equilibrium curve of Equation 2 suggests that, theoretically, a stock can be maintained indefinitely at any B/k level, stocks below B/k = 0.5 are likely to fail to fulfill their natural roles as prey and/or predator (Pauly and Froese 2021; Scotti et al. 2022) and the resulting changes in the ecosystem are likely to make such fisheries unsustainable.

Figure 2. 

Schematic representation of the surplus production model used by CMSY, with indication of impaired recruitment due to small stock size, where F_MSY is reduced linearly with decline in biomass.

A recent study by Bouch et al. (2021) showed that if CMSY was applied with default prior settings and compared against data-rich assessments, 14 of 17 final biomass estimates (82%) differed less than 50% from official biomass estimates, while a data-moderate method (SPiCT) only had 6 (35%) estimates in that range (fig. 4 in Bouch et al. 2021). On the other hand, Ovando et al. (2022) criticized default priors in catch-only methods and highlighted the importance of setting reasonably accurate priors of r–k for providing an accurate estimate of stock status. Sharma et al. (2021) found that CMSY results for the final year improved substantially with the quality of the prior information. Pons et al. (2020) reported that CMSY reduced the bias in estimates of F/F_MSY by more than 50% compared to Catch-MSY based on a simulation study. Zhou et al. (2018) also found that CMSY is more accurate than the original Catch-MSY method and that it performed well in estimating the final relative stock size of 13 Australian stocks in comparison with the official age-structured assessment estimates. As a further confirmation of usefulness, CMSY was fitted with results that were deemed reasonable to a multitude of species and stocks globally (Cheung et al. 2022), in Europe (Froese et al. 2018a), Eastern Mediterranean (Demirel et al. 2020; Saygu et al. 2023), Canada (Schijns 2020; Schijns et al. 2021), China and surroundings (Zhang et al. 2018; Liang et al. 2020; Roa-Ureta 2020; Wang et al. 2020a, 2020b; Kang et al. 2022), Africa (Musinguzi et al. 2020; Palomares et al. 2020a), Mexico (Ferrer et al. 2022), and other countries and regions (Cruz et al. 2020; Palomares et al. 2019, 2020b, 2021) (Fig. 3). Interestingly, CMSY has recently been extended to a bio-economic stock assessment using price data in addition to the time series of catch (Lancker et al. 2023).

Figure 3. 

Maps showing the locations of the centroids of the over 2000 stock assessments performed with CMSY (~20%) and CMSY ++ (~80%) in all parts of the world. Generally, the assessments in wealthier countries (USA, Canada, Australia, New Zealand, EU-member states) were complemented with CPUE or other ancillary data, and well-informed priors for the terminal B/k values; such information was often lacking for many countries in low latitudes, but the results still matched what was known of their fisheries. Based on original Sea Around Us data, and previous analyses in Froese et al. (2018a); Zhang et al. (2018); Demirel et al. (2020); Cruz et al. (2020); Schijns (2020); Liang et al. (2020); Roa-Ureta (2020); Wang et al. (2020a,b); Musinguzi et al. (2020); and Palomares et al. (2019, 2020a, 2020b).

Figure 4. 

Plot of MSY prior derived from maximum catch over MSY estimated with BSM for the 400 stocks used for training ANN. The outliers are stocks where catches never exceeded MSY, for which neither CMSY nor the method to drive MSY priors should be used. The dashed 1:1 line indicates identical values whereas the dotted lines indicate deviations of ±50%. Note that CMSY++ estimates of MSY would fall vertically between the MSY priors and the 1:1 line.

However, there remained a stumbling block also identified by Bouch et al. (2021): the use of rigid constraints, either entered by the user or generated by the software as internal heuristics, that the biomass trajectory was forced to accommodate. Setting appropriate initial biomass priors at the start of the time series (B_start/k) is a general challenge in cases where the catch series is short and does not include historical catches that would reflect the initial lightly exploited stock biomass (i.e., B_start/k = 0.9–1.0). Noting this, the ICES benchmark workshop on SPiCT assessments (ICES 2021: WKMSYSPiCT), for example, recommended that priority should be on reconstructing catch time series, but where this is not possible and initial catches are already close to the observed maximum, the B_start/k prior should be set to 0.5 or lower.

The remaining constraints, i.e., the biomass priors in the Bayesian context, were:

• The fraction expressing the biomass depletion (from carrying-capacity) assumed to have already occurred at the start of the time series (B_start/ k);
• The fraction expressing the biomass depletion assumed to have occurred at some intermediate year of the time series (B_int/ k); and
• The fraction expressing the biomass at the end of the catch time series (B_end/ k).

From CMSY to CMSY++. Following the publication of the CMSY method (Froese et al. 2017), the code implementing it underwent a series of improvements, also in response to feedback from its users. Notably an option to consider the degree of technological creep, i.e., an increase in catch per unit of effort (CPUE) not caused by increase in biomass, was introduced based on Palomares and Pauly (2019).

In a workshop held in November 2019 in Thessaloniki, Greece, it became evident that the Bayesian model (BSM) requiring time series of catch and abundance as main input could also be run without abundance data, making it a “catch-only” Bayesian method that could replace the Monte-Carlo approach used in the original CMSY. In other words, CMSY++ and BSM are nested within a single JAGS model and use the same parameterization and catch input, the only difference being that CMSY++ has no input of abundance data. This enables a consistent and continuous transition to fitting abundance indices to as little as two observations should such information become available. In cases where the abundance index is informative about the trend of the abundance trajectory, the user has the option to relax or disable the terminal and intermediate depletion priors. This is particularly relevant for estimating the short-term response to management interventions such as catch quota reductions (Wetzel and Punt 2015).

Also, in Catch-MSY and CMSY, a prior for carrying capacity was derived from the reasoning that a lower limit of k should be larger than the highest observed catch, that an upper limit of k should be 10–100 times higher than the lower limit, and that higher productivity r would suggest a narrower prior range of k. Building on the good correlation between maximum catch and MSY observed in hundreds of stocks (Fig. 4), CMSY++ instead derives a heuristic prior for MSY from maximum catch and then obtains a prior for k from k_prior = 4 MSY_prior/r_prior. Note that this method only works if catches in the time series have exceeded MSY (see Fig. 4 for examples where that requirement was violated).

Another important improvement was the replacement of the rigid uniform r–k prior space with a multivariate lognormal (MVLN) prior that accounts for the negative correlation between k and r within a population, and where lower probabilities are assigned to peripheral r–k pairs further away from the core of the ellipsoid r–k distribution (Fig. 5b). The MVLN prior for r–k is a function of the means and the covariance matrix of log(r_prior) and log(k_prior); see details below.

Figure 5. 

Examples of graphical output of CMSY++, here for European plaice (Pleuronectes platessa) in the eastern English Channel. Panel (A) shows the time series of catch from 1980 to 2011, with the thin blue curve representing smoothed catch and the red circles the smoothed minimum and maximum values. Panel (B) shows as dotted box the prior range for r and k. The dots in light grey indicate potential r–k pairs and the dark grey dots indicate pairs determined as viable by the catch-only CMSY++ analysis. The blue cross indicates the best CMSY++ estimate for r–k, with approximate 95% confidence limits. The red cross indicates the corresponding estimate derived from catch and CPUE by BSM. Panel (C) shows the time series of relative biomass B/k as estimated by CMSY++ (blue curve) and BSM (red curve) with dotted 95% confidence limits. The grey points indicate the available CPUE data. The horizontal lines indicate B_MSY at 0.5 k and B_lim at 0.25 k. The vertical purple line in the lower left corner indicates the B/k prior set by the user to 0.01–0.1. The dotted vertical lines in 2005 and 2011 are the prior B/k ranges set by the Neural Network. Panel (D) compares the density of the light-grey B/k prior set by the user for 1980 with the corresponding dark-grey posterior density estimated by BSM.

The challenge of deriving more rigid biomass priors was also addressed in another major development. An Artificial Neural Network (ANN) now provides the option to predict default relative biomass priors (B/k) from catch relative to prior MSY, based on traits of catch patterns derived from hundreds of test stocks (but see discussion below, stressing that ANN is just a “convenience-add-on”, meant to assist users in selecting the best available prior information).

Other improvements include:

• The rigid, uniform prior B/ k ranges of CMSY were replaced with beta-distributions which have the desirable property for biomass depletion priors of being bounded by 0 and 1 (Winker et al. 2018) with increasing skewness as either end of the spectrum is approached; this resembles the uneven distribution of relative biomass around the equilibrium curve in Fig. 6.
• To correct a remaining bias of r– k pairs towards high k and low r, the lower right focus of the ellipse containing “viable” r– k pairs was used as “best” estimate, with confidence intervals derived from pairs where r was larger than the 25 ^th percentile of “viable” r.

In summary, the purpose of this study is to present the history and latest developments of the catch-only CMSY++ method, and to compare its predictions with those of a regular surplus production model, which has time series of abundance information as additional input, everything else being equal.

Figure 6. 

Scatterplot of relative biomass (B_t/k) over relative catch (C_t/MSY), both estimated with BSM, with 18 341 points for 400 stocks. The blue curve is the equilibrium biomass prediction from Equation 2. The vertical blue line indicates the range that contains 90% of the (B_t/k) points for catches above MSY. The red dashed lines indicate approximate 95% confidence ranges for prior B_t/k.

Description of the stocks used for preliminary testing and training. A data set with times series of catch and abundance for 400 different stocks was assembled to train the Artificial Neural Network and to understand the correlation between the MSY prior derived from maximum catch and MSY estimated by BSM from catch and abundance (Fig. 4). The stocks stem from 11 large marine ecosystems covering 87 marine ecoregions worldwide, including 10 Arctic stocks, 101 North Pacific stocks, 181 North Atlantic stocks, 14 South Pacific and South Atlantic stocks, 36 tropical stocks, 5 stocks from South Africa, 26 stocks from Australia, and 27 wide-ranging stocks. About three quarters of the species are demersal and one quarter is pelagic. There are 321 teleost, 22 elasmobranch, and 57 invertebrate stocks. Invertebrates are mainly represented by crustaceans (shrimps and lobsters) and mollusks (bivalves and cephalopods). Resilience categories range from very low to medium. Recent biomass was severely depleted in 106 (27%) of the stocks. Twenty-five stocks (6%) had recent biomass close to the unexploited level. The longest time series started in 1876 (with the actual years analyzed starting in 1940) and the shortest in 2005. For use in training, some of the time series were shortened to exclude recent periods where e.g., declining catch was not caused by low or declining biomass and thus the assumed relation between catch and biomass was broken (see also discussion of caveats below). Catch and abundance data were derived from official stock assessment reports such as summarized in ICES Stock Assessment Graphs (https://standardgraphsicesdk/stockListaspx) for the Northeast Atlantic and in NOAA Stock SMART (https://wwwstnmfsnoaagov/stocksmart?app=browse-by-stock) for North America. More details on the stocks and the sources are available in the files Train_Catch_9.csv, Train_ID_9.csv and Out_Train_ID_9.csv, which are available from https://oceanrep.geomar.de/53324/.

Simulated data. Simulated catch and CPUE data (24 stocks) were created so that the simulated parameter values and stock status estimates were “true” known quantities for performance evaluation. For convenience, k was set to 1000 and r was set at 0.06, 0.25, 0.5, and 1.0 to represent species with very low, low, medium, and high resilience, respectively. For a simulated time series horizon of 50 years, biomass patterns of continuously high, continuously low, high to low, low to high, low–high–low, and high–low–high were created. The desired patterns were produced by inserting high or low catches into Equation 2 and calculating the biomass in subsequent years. The simulated data and the CMSY++ results are available from https://oceanrep.geomar.de/53324/ [files SimSpecCPUE_4_NA.csv, SimCatchCPUE_4.csv, Out_July082021_SimSpecCPUE_4_NA.csv, CMSY++16_Sim_8.R].

Derivation of MSY_prior and multivariate-lognormal r–k distribution. Similar to other well-established surplus productions models such as SPiCT (Pedersen and Berg 2017) or JABBA (Winker et al. 2018), CMSY++ assumes lognormal prior distributions for r, k, and MSY, thus avoiding negative values in ranges of uncertainty for these parameters. Building on the good empirical relation between maximum catch and plausible MSY ranges observed among the 400 stocks used for testing (Fig. 4), a prior for MSY was attained as follows: if the time series of catch was more or less flat or ascending such that the maximum catch occurred in its last 5 years, then the mean of the three highest catches was taken as prior for MSY. This was done because a flat or ascending time series without recent decline was deemed unlikely to have exceeded MSY by much, if at all. In contrast, if catches were declining after a peak, that peak was likely overshooting MSY. Therefore, under such conditions, ¾ of the mean of the 5 highest catches in the time series was taken as prior for MSY. The mean of three or five catches was chosen to reduce the misleading impact of single, extraordinary high catches. Note that this procedure assumes that MSY is equal to or smaller than the highest catch values, i.e., this approach and CMSY in general are not suitable for lightly exploited stocks where catch never approached MSY (see fig. 3 in Martell and Froese 2013), or where such catches were not included in the time series.

Log(r_prior) is derived from life-history traits and log(MSY_prior) is derived from maximum catch, i.e., these methods of derivation are uncorrelated, there is no circularity in the derivation, and the priors can thus be drawn from lognormal distributions without violating statistical assumptions about independence. These priors then provide the solution for k_prior = 4 MSY_prior/r_prior. Note that if no variability were assumed for MSY_prior, this would result in a fixed log(k)–log(r) correlation of –1 (equation 5 in Froese et al. 2017). In reality, correlations vary between zero and –1, e.g., between –0.44 and –0.98 in the BSM results for the 400 training stocks (see Out_Train_ID_9.csv available from https://oceanrep.geomar.de/53324/).

The Schaefer model can be expressed as a function of r and MSY, without k (Equation 3); however, this arrangement does not change the dynamics of the model and the new term for surplus production seems less intuitive than the original one (Equation 1).

$B_{t+1}=B_t+r·B_t-\frac{{\left(r{B}_t\right)}^2}{4\mathrm{MSY}}-C_t$ [3]

To retain the original form of the CMSY base model (Equation 2) with parameters r and k, the within-stock correlation between r and k was accounted for in a MVLN distribution implemented by: (1) drawing a large sample (n = 10 000) of independent random deviates of log(r̃) and log($\widetilde{MSY}$) from their prior distributions; (2) computing the corresponding log(k̃) = log(4) + log($\widetilde{MSY}$) – log(r̃) and (3) computing the means and the covariance of log(r̃) and log(k̃), which are (4) then passed on as covariance matrix for the r–k ~ MVLN prior in the CMSY++ and BSM model formulations (see bsm() function in CMSY++16R code, which is available from https://oceanrep.geomar.de/53324/). The biomass dynamic in Equation 3 was implemented as a Bayesian state-space model that accounted for random variability in population dynamics (process error) and catch (observation error) (see Equation 1). This way, biomass over time was modelled as a sequence of random variables. This avoided the model to be completely driven by priors, which occurs when random variables are linked through a deterministic function (Borel’s paradox: Schweder and Hjort 1996).

Application of an Artificial Neural Network in prediction of B/k priors. A feed-forward Artificial Neural Network (ANN) (Fritsch et al. 2019) was chosen for classifying stock status as being above or below the MSY level to accommodate Equation 4. ANN input consists of characteristics of the catch time series such as overall shape, difference between minimum and maximum catch, and slope in the first and final years. The network was trained with time series of 400 stocks (see detailed description above), which were selected to reflect the interplay of their catch and abundance data as described by Equation 2, so that the ANN could detect and learn typical patterns that allowed for the prediction of relative biomass (B/k) priors from relative catch (C/MSY). The time series of B/k were estimated with BSM from catch and abundance data and were treated as “true” for the purpose of the training. Specifically, ANN was set to predict whether B was above or below B_MSY for the start year, an intermediate year, and the final year of the time series. ANN was designed to have one classification output neuron that simulated a variable A_t with a binary probability distribution with values of either 1 (B/k above 0.5) or –1 (B/k below 0.5) The equilibrium B/k prior for reference year t, with catch value C_t, was then derived as

$(B/k{)}_{t\text{ prior }}=\frac{1+A_t\sqrt{1-C_t/{\mathrm{MSY}}_{\text{prior }}}}{2}$ [4]

Note that Equation 4 only gives real number solutions if C_t < = MSY_prior. Therefore, its application was restricted to cases where C_t < 0.99 MSY_prior. The optimal ANN topology was found using the growing strategy (Bishop 1995), while performing 20-fold cross-validations with 95% to 5% random separation into training and test data sets and was made up of one hidden layer with 71 neurons. As an alternative model, we also tested a Long Short Memory “end-to-end” model (Hochreiter and Schmidhuber 1997) that accepted the entire catch time series as input but did not produce better predictions. During the cross-validation process, the ANN accuracy was assessed based on the calculation of correct classifications of relative biomass (B/k) being above or below the MSY threshold for the start, intermediate and end years. Using more classes or continuous output ended in lower accuracy and more frequent over-fitting, i.e., in a lower generalization capacity. Moreover, using ANNs instead of other models was justified by the importance of processing the time series data as a whole, i.e., by automatically modelling inter-sample correlation, which is a primary driver of time series classification accuracy (Coro et al. 2021). After selecting and training the optimal model, we used simulated stocks to assess the ANN prediction accuracy of biomass being above or below B_MSY for the start, intermediate, and end year of the time series, while using the “true” values to calculate the percentage of the correct predictions. Also, the “true” B/k value in the last year was compared with the respective CMSY++ estimate, with approximate 95% confidence limits.

How uncertainty of B/k priors was established. Equation 4 describes how a point estimate of relative equilibrium biomass (B/k) was derived from catch relative to MSY. Catch and biomass are rarely in equilibrium in real world stocks and the width and shape of uncertainty vary with the position of the equilibrium point estimate in B_t/k–C_t/MSY space (see distribution of points around the equilibrium curve in Fig. 6). As a pragmatic solution, ranges of uncertainty were derived as follows: (1) the 5^th and 95^th percentile of B_t/k values and the median of C_t/MSY values were determined for all points where catch exceeded MSY (C_t/MSY > 1); (2) B_t/k values that bracketed most of the variability for the case of close to zero catches (C_t/MSY = 0) were chosen for high and low biomass; (3) the B_t/k values in (2) were treated as intercepts of linear regression lines that connected them to the 5^th and 95^th B_t/k percentiles determined in (1); and (4) for catches larger than 1.21 MSY (the median of catches above MSY), a fixed range of uncertainty in B_t/k was used. The resulting equations for prior B_t/k ranges as a function of C_t/MSY and their being above or below B_MSY are shown in Table 1.

For example, for a catch of 0.5·MSY and a biomass below 0.5 · k, the dashed red lines in Fig. 6 propose a prior B/k range of 0.11–0.42; for catches at or above MSY, the parallel red lines propose a prior B/k range of 0.26–0.72. Remember that the pattern-based ANN B/k priors thus derived are only a “convenience-add-on” meant to assist CMSY++ users in evaluating and objectively selecting the best possible prior ranges for the analysis (see also discussion of biomass priors below).

Derivation of equilibrium curves. The equilibrium curve for the interplay between relative biomass (B/k) and relative catch (C/MSY) for the modified Schaefer model shown in Figs 5 and 7 was derived from

Figure 7. 

Scatterplot of 4805 observations of abundance relative to maximum abundance for 94 stocks where maximum abundance was deemed close to unexploited (B/k) and catch relative to a prior for MSY derived from maximum catch, i.e., no modelling was involved in generating the data. The upper blue curve represents the modified Schaefer model (mSchaefer) used by CMSY++. The middle black curve represents the Fox model. The lower red curve with approximate 95% confidence limits represents 14 stocks assessed with the Stock Synthesis model (SS3). The short green bold line indicates the median of relative population size = 0.497 for available points from 0.95 to 1.05 relative catch levels.

$\frac{C}{\mathrm{MSY}}=\left(4\frac{B}{k}-{\left(2\frac{B}{k}\right)}^2\right)·\mathrm{RC}$ [5]

where RC stands for recruitment correction with RC = 4 B/k if B/k < 0.25 and RC = 1 otherwise (same as in Equation 2)

The equilibrium curve for the Fox (1970) model shown in Fig. 7 was derived from

$\frac{C}{\mathrm{MSY}}=e\frac{B}{k}\left(1-\log \left(e·\frac{B}{k}\right)\right)$ [6]

where e stands for Euler’s number 2.718.

For comparative purposes, equilibrium yield curves were extracted from 14 Stock Synthesis models, which had been used for quota advice for tunas, billfishes, hakes, monkfish, snapper, herring, and sardine by national or Regional Fisheries Management Organizations including NOAA, ICES, ICCAT, and IOTC. All Stock Synthesis models had been fitted assuming a Beverton–Holt stock-recruitment function. In Stock Synthesis, the equilibrium curves are computed internally based on the age-structured equilibrium dynamics (cf. Winker et al. 2020) and can be extracted using the R package “r4ss” (Taylor et al. 2021). The median of the 14 equilibrium curves with approximate 95% confidence limits (2.5^th and 97.5^th percentiles) are illustrated in Fig. 7, where the spawning stock biomass ratio SSB/SSB₀ corresponds to B/k and SSB_MSY/SSB₀ corresponds to B_MSY/k where yield (surplus production) is at its maximum (MSY).

All data and code used in this study are available from https://oceanrep.geomar.de/53324/ and https://github.com/SISTA16/cmsyPlusPlus.

Cross-validation of the ANN predictions. The task of ANN was to predict from properties of the time series of catches whether relative biomass (B/k) was above or below the MSY threshold (B_MSY/k) in a given year. The percentages of correct classifications are presented in Table 2 and are 99%–100% for stocks that were included in the training, as expected. For stocks that were excluded from the training, cross-validation accuracy ranged from 68% accuracy for the first year in the time series to 91% for the intermediate and final year (Table 2). Of the 400 real stocks used in the training of ANN, 377 (94%) BSM estimates of final B/k were within the approximate 95% confidence limits of the respective CMSY++ estimate. Of the 23 mismatches, 13 were cases where BSM estimated relative biomass of less than 0.1 B/k, and the approximate lower confidence limit of the CMSY++ estimate was also below 0.1 B/k, albeit above the BSM estimate, i.e., both methods predicted the stock as severely depleted. Detailed results are available in the file Out_Train_ID_9.csv that can be retrieved from https://oceanrep.geomar.de/53324/.

Performance of ANN and CMSY++ against simulated data. The results of applying CMSY++ with B/k priors predicted by ANN to simulated stocks are given in Table 3. ANN made correct predictions of biomass being above or below B_MSY in 10 of 21 applicable cases of start biomass (48%), 16 of 24 cases of intermediate biomass (67%), and 12 of 23 applicable cases of end biomass (52%). The “true” value of B/k was contained in the approximate 95% confidence limits of the CMSY++ estimate in 12 of 24 cases (50%). Of the 12 cases where the ANN predictions for end biomass were correct, 10 CMSY++ predictions (83%) contained the true value in their approximate 95% confidence limits. Of the 11 applicable cases where ANN predictions for relative end biomass were wrong, 9 CMSY++ predictions (82%) were also wrong. Note that ANN was on purpose not trained on the artificial, often unrealistic and sometimes extreme catch patterns of the simulated stocks, which were designed to test the limits of the method; see Discussion below.

ANN and CMSY++ performance against real and simulated stocks. CMSY++ performed well (68%–91% correct ANN predictions in cross-validation, see Table 2; 94% correspondence of CMSY++ predictions for final biomass with BSM results) against 400 real world stocks, which had been selected because their interplay between catch and biomass largely followed Equation 2. We believe that it is more appropriate and informative to compare data-limited results (here from CMSY++ with catch data) against results obtained with the same accepted model, but with substantially more data (here from BSM with time series of catch and abundance), rather than against very different models with very different data requirements and assumptions. Exploring in depth the differences in results obtained from surplus production models (such as BSM) versus data-rich age-structured models is beyond the scope of this study; but see the discussion around Fig. 7 below for the existence of and possible reasons for some of those differences.

Not surprisingly, ANN predictions for biomass being above or below B_MSY were less satisfactory for simulated stocks (only 48%–67% correct predictions, Table 3), because several of the simulated stocks had catch and abundance patterns which purposely were not present in the real-world data on which ANN had been trained. For example, 50 years ago many stocks were still underexploited (confirm Fig. 1) and ANN thus correlated low catches in that period with large stock size. Instead, half of the stocks in the simulation were set to start with low catches and low biomass, thus causing about half of the wrong ANN classifications (5 A/B of 11 wrong classifications for the start year in Table 3). Similarly, the simulations used high catches in the first year to bring down high biomass and to force the desired depletion patterns, whereas in the real world, high catches at or above MSY and high biomass rarely occur together (see Figs 5 and 7). Other rarely found patterns in the real world pertain to the simulations that included very lightly exploited stocks where catch never exceeded MSY or biomass never fell below B_MSY and thus were never fully exploited (the High–High scenarios in Table 3). While ANN and CMSY++ underestimated final biomass in all four High–High scenarios, the “true” B/k values from the simulations were still contained within the confidence limits of the predictions. In other words, with the training applied to ANN, the results that will be obtained from simulation exercises fully depend on the resemblance of the simulated scenarios to real world stocks. Instead, the simulations were used to explore the behavior of CMSY++ in extreme scenarios to better understand its limitations.

In addition, priors are as important as data in a Bayesian context, especially in data-limited applications, and it should not come as a surprise that wrong input (here: wrong prior information about the likely B/k range) led to wrong results. In contrast, the simulations suggest that if the final B/k prior range is broadly set correctly, then there is a high probability that CMSY++ will give reasonable predictions of stock status (see fig. 4 in Bouch et al. 2021 and results by Sharma et al. 2021 for independent confirmation). The main reason for the introduction of ANN was to make the derivation of the default B/k prior ranges more objective, with unknown rules developed and applied by the neural network. However, it is important to stress again that these default B/k priors are just a convenient add-on to CMSY++, and that the ANN predictions are to be replaced by evidence-based prior knowledge of stock status wherever possible (Table 4).

Addressing some common misconceptions. In medicine, asking a patient (or others who know that person) specific questions about their medical history, a process called anamnesis, is an essential part of formulating a diagnosis and developing a plan for recovery and wellbeing. The similarities to the process of stock assessment and management are obvious. Yet, one of the most common criticisms of CMSY is its strong dependence on such anamnesis or “anecdotal” knowledge about past and present fishing pressure or stock status. Some have even suggested the dependence of CMSY on B/k priors is so strong that the analysis might as well be skipped, and the priors be used directly for stock status classification. This would be analogous to using the anamnesis directly and forego its subsequent verification in the full diagnostic examination, surely not a serious proposal in a medical context.

Figure 5 shows an example of a misleading prior being strongly corrected by the data, here for European plaice (Pleuronectes platessa Linnaeus, 1758) in the eastern English Channel. Because abundance of that stock was very low in 1980 (Fig. 5a; ICES 2020), the user had set the B/k prior for that year to 0.01–0.1, with a central value around 0.05 (lower left corner of Fig. 5c). However, the B/k posteriors of both CMSY++ and BSM overlap only marginally with that prior and instead estimate a central value close to 0.2, i.e., about 4 times higher than suggested (Figs 5c–5d). In other words, the prior user perception of stock status was not compatible with the available data, and the change from uniform to beta-distribution of the B/k prior allowed for the substantial correction of the prior knowledge.

Another misunderstanding of the Bayesian approach is the use of very wide, uninformative priors with the explicit purpose of reducing their influence on the results, e.g., by providing a uniform B/k prior range for final biomass of 0.1–0.9 k. Such prior informs the analysis that, with equal probability, the stock may be nearly collapsed or nearly unexploited. We are not aware of a single stock where such statement would be true. In other words, the objectivity that an uninformative prior is supposed to bring to the analysis is in reality the feeding of knowingly erroneous input to the model. Instead, in recognition of the importance of a realistic prior for a realistic analysis, real effort must be invested to determine the best possible prior information. We stress again that the built-in B/k prior predictions by ANN are a not-required add-on of CMSY++, to be replaced by independent B/k prior knowledge whenever possible (Table 4). To that end, stock assessments whose results will be used by managers should be done case-by-case, and not in a batch mode, so that due attention is paid to selecting the appropriate priors.

Setting appropriate initial biomass priors at the start of the time series (B_start/k) is not specific to CMSY++, but a general challenge for parameterizing surplus production models in cases where the catch series is fairly short and does not include historical catches that would reflect the initial lightly exploited stock biomass (i.e., B_start/k = 0.9–1.0).

There may also be concerns that deriving priors from the time series of catch data violates the requirement of Bayesian prior beliefs to be established before the data are considered. We agree with this principle and Table 4 gives examples on how independent beliefs about stock status can be translated into numerical prior ranges. Only for cases where such information is not available should the priors proposed by CMSY++ be used, which are based on ANN having looked at patterns in the catch data and comparing them with catch patterns of 400 stocks for which the biomass was known.

The CMSY user guide (available from https://oceanrep.geomar.de/id/eprint/52147/) provides a table with suggested B/k ranges according to the perception or “narrative” about the depletion of the stock. This approach is expanded upon in Table 4, giving examples of typical terms used to describe depletion levels. It is hard to imagine a fisher, fish processor or fisheries manager being unable to correctly assign an important stock to one of these broad categories. Also, with the exception of stocks that are nearly unexploited or collapsed, the proposed B/k ranges in Table 4 span 40% of the possible range and include, respectively, total collapse to near sustainability (0.01–0.4 k), outside of safe biological limits to securely sustainable level (0.2–0.6 k), and below B_MSY to near k (0.4–0.8 k). In other words, these ranges do not pre-empt the CMSY++ analysis or unduly predetermine the results. It makes a big difference for managers whether a stock is likely to be recruitment-impaired (0.2 k) or safely above the MSY threshold (0.6 k). If a stock seems to fall between two of these categories, intermediate prior B/k ranges can of course be used. If length frequency data are available, these can be used to obtain independent and objective B/k priors (Froese et al. 2018b, 2019; Musinguzi et al. 2020).

Caveats to using CMSY++. We have argued above that in the absence of abundance data, the catch-only implementation of CMSY++ performs reasonably well in classifying the stock status and thus in assisting the prioritization of management interventions (cf. Sharma et al. 2021). However, CMSY++ does strongly depend on accurate catch data and if, for example, a reduction in catch is enforced by management, a catch-only approach has clear limitations in monitoring its effectiveness to meet the rebuilding objectives (Wetzel and Punt 2015). This is because validating that the model can correctly predict the stock’s response will ultimately require more observations (Kell et al. 2021). If, for example, a fishery is completely closed by law, the zero-catch signal does not allow any inferences about stock size development. Ideally, rebuilding of stocks should be aligned with data collection programs designed to monitor the progress and provide a feedback control to stakeholders and to BSM, which is included in the CMSY++ package and automatically activated if abundance data are provided. This enables a seamless transition from the CMSY catch-only approach to regular surplus-production modelling (BSM).

As indicated above, catch patterns can be used to make empirical predictions about relative stock sizes at selected points (e.g., start, end, intermediate) in the time series. In CMSY, that was done by a list of empirical if-then rules; in CMSY++, this is now done by an Artificial Neural Network (ANN). However, in both cases empirical predictions of relative biomass from catch patterns will only work if the interplay between catch and biomass is mainly driven by Equation 2. That will be the case if a more or less constant fishing effort is applied or if management follows a harvest control rule and sets catch limits based on relative stock size. It will, however, not work if strong changes in catch occur for other reasons, such as drastic variation in demand, or a species being newly protected from fishing, or declining carrying capacity because of warming waters, or increasing carrying capacity because the stock is released from predation mortality because a main predator has collapsed or disappeared. The presence of such circumstances has to be considered by the local experts, and the default biomass priors may then have to be corrected accordingly. To help with a better understanding of cases where CMSY++ works well with its default settings and cases where expert knowledge is required to get meaningful results, a number of selected stock assessments is presented and discussed in the Suppl. material 1.

More generally, especially in depleted stocks, a minor overestimation in stock size, such as estimating final B/k as 0.2 instead of 0.1, will lead to a substantial underestimation of fishing mortality. In addition, especially in large species where several age classes contribute to the catch, the contribution of early year classes may already be substantial although they are not yet fully selected by the gear. This reduces their specific F and the overall estimate of fishing mortality if compared with official assessments, which typically base their estimate of F only on fully selected age classes (see examples in the Suppl. material 1). Therefore, management based on surplus production models such as CMSY++ or BSM should focus on predicted biomass and not on predicted fishing mortality, which may be underestimated for fully selected age classes.

Comparing the modified Schaefer model to other models. In the course of searching for relative biomass priors for CMSY++, we realized that the equilibrium biomass predicted by the modified Schaefer model provided a very reasonable fit for the widely scattered catch and biomass data of the 400 stock assessments that we examined (Figs 6 and 7). In contrast, problematic implementations of surplus-production models, such as those of Fox (1970) and Pella and Tomlinson (1969) became apparent (Fig. 7). As Pauly and Froese (2021) pointed out, these models are based on unjustified departures from the ecological foundation of the original model proposed by Graham (1935), operationalized by Schaefer (1957), and modified to account for reduced recruitment at low stock sizes by Froese et al. (2017) (equation 2, fig. 2).

The critique of Pauly and Froese (2021) centered on the low estimates of the minimum biomass required to produce MSY, fixed at 37% of carrying capacity in the Fox model, with some assessments using target values as low as 30% in Pella and Tomlinson (1969) or age-structured models (MRAG 2016). Recall that MSY is instead generated at 50% of carrying capacity in the original Schaefer model, based on the widely corroborated logistic curve of population growth (Verhulst 1838).

An apparently overlooked objective comparison of different models is how well their predictions fit observed catch and abundance data across many stocks. To avoid any confounding model assumption effects, a model-independent approach was used to generate points of the ratios Catch/MSY and B/k in Fig. 7, which is a novel feature that allows an objective comparison of the assumptions and prediction of the various models and should be explored in more depth in subsequent studies. Instead of using one of the above-mentioned assessment models to predict relative catch and relative biomass, MSY values were derived based on the maximum catch as described for MSY prior generation above and B/k observations were approximated as observed CPUE relative to the highest CPUE apparently close to representing carrying capacity. A subset of 94 of the 400 real stocks analyzed in this study met these requirements. While one could argue about the appropriateness of the methods to standardize the biomass and catch values shown in Fig. 7, the point here is that they did not favor one model over the other.

Similar to Fig. 6, Fig. 7 shows that the equilibrium yield curve of the modified Schaefer model traces the highest density of points reasonably well, whereas the Fox model traces only the lowest biomass values for a given catch. Equilibrium yield curves, which were extracted from age-structured assessments using Stock Synthesis (SS3; Methot and Wetzel 2013) with an integrated Beverton and Holt stock-recruitment function fell completely outside of the range of points at low stock sizes (i.e., B/k < 0.25). In other words, the Fox and SS3 models systematically overestimated the average productivity (i.e., the equilibrium yield curve) at low stock sizes, which can have severe implications for rebuilding of stocks (cf. Hutchings 2015; Perälä et al. 2022). Mangel et al. (2013) also criticize data-rich models for their strong dependence on typically unknown inputs such as the steepness parameter of the Beverton and Holt stock recruitment relation (Mace and Doonan 1988) or the rate of natural mortality, the arbitrary settings of which may then pre-determine key outputs, such as the biomass required to produce MSY.

The model comparison presented in Fig. 7 is preliminary and qualitative and should be repeated with a larger number of stocks, considering that some of the low catches at lower stock sizes may have been mandated by harvest control rules. However, Fig. 7 clearly points at a potentially large and highly risk prone bias associated with models that predict maximum sustainable catches at low stock sizes and overestimate the population growth potential at dangerously low stock levels.

Summary. CMSY++ has developed into a versatile integrative method that can incorporate, in addition to the required catch time series, abundance data and a wide variety of ancillary information (e.g., Froese et al. 2018b) in a rigorous Bayesian context that tends to reduce the dependency on prior information while remaining robust and thus usable in data-limited situations such as common in many parts of the world.

The majority of the independent tests of CMSY used the default priors and thus did not really test the CMSY method per se. With CMSY++, such tests would reproduce the 9%–32% failure rate of the Artificial Neural Network, with even higher percentages if applied to stocks whose catches were reduced for reasons external to the dynamics of the fish population in question, such as changes in market demand or environmentally-mediated productivity. Similarly, any failure rate can be produced with simulated stocks that deviate substantially from the 400 stocks used in training ANN. Instead, in order to be more realistic, tests should assume that local experts are able to provide priors that are not wider than about 40% of the maximum possible range and that include the “true” value.

CMSY++, either applied as a data-poor (catch-only) or preferentially as a data-moderate (catch and CPUE) method, allows the assessment of stocks for which at least catch data are known. That is especially important for data-poor areas that have been generally excluded from Ecosystem Based Fisheries Management (EBFM) programs (Link 2010) because of ignorance regarding stock status (Townsend et al. 2019). In addition to the large number of data-deficient areas, stocks that have low commercial interest have generally been overlooked in assessments and the conservation status of marine megafauna is unknown in many areas, such as the Mediterranean and the Black Sea (Stergiou and Tsikliras 2006). CMSY++ and related methods (e.g., Winker et al. 2018; Froese et al. 2018b, 2020) can provide assessments for these important ecosystem components. Such standardized and comparable stock assessments applied on a global scale (Fig. 3) will contribute to a much-needed better understanding of the world’s fisheries and ecosystems (Pauly et al. 2018).

Rainer Froese acknowledges support by the German Federal Nature Conservation Agency (BfN); Maria-Lourdes D Palomares and Daniel Pauly acknowledge support from the Sea Around Us, itself funded by a number of philanthropic foundations, notably the Minderoo Foundation, which underwrote the work leading to the update of the reconstructed catches to 2018 and the majority of the CMSY++ assessments in the global map shown as Fig. 3. We also thank Nicolas Bailly and Elizabeth Bato David for their assistance in creating this map. Athanassios C Tsikliras was partly supported by the European Union’s Horizon 2020 Research and Innovation Program (H2020-BG-10-2020-2), grant number No 101000302—EcoScope (Ecocentric management for sustainable fisheries and healthy marine ecosystems). The article reflects only the authors’ view and the Commission is not responsible for any use that may be made of the information it contains.

All data, scripts, and supplementary material are available from https://oceanrep.geomar.de/53324/ and https://github.com/SISTA16/cmsyPlusPlus.