Development of a new control rule for managing anthropogenic removals of protected, endangered or threatened species in marine ecosystems

Operating models

An operating model is a “digital twin” of the socio-ecosystem to be managed. It is a mathematical (and hence idealized) representation of the phenomenon under study that encapsulates important knowledge about the phenomenon (Punt et al., 2016), in this case population dynamics of PETS. In the case of cetaceans, commonly used operating models include a population dynamics model that can be age-aggregated (Wade, 1998; Wade, 2002) or age-disaggregated (Wade, 2002; Brandon et al., 2017). The interested reader is referred to Punt (2016) for a review of operating models developed for cetaceans; to Tuck (2011) and Tinker et al. (2022) for seabird operating models; and to Chaloupka & Balazs (2007) and Snover (2008) for a discussion on a turtle operating model. These operating models are variants of SPM. In the remainder, the operating model is age-disaggregated (that is age-structured), it assumes a Pella-Tomlinson functional form for density-dependent fecundity (see Algorithm 2 in Genu et al., 2021 for a full description of the model) and no Allee effect (see e.g. Haider et al., 2017).

Control rules for managing removals

There are two main control rules currently used for managing marine mammal removals: the Potential Biological Removal (PBR; Wade, 1998) and the Removals Limit Algorithm (Cooke, 1995).

Potential Biological Removal

The calculation of the PBR is straightforward (Wade, 1998): ${\displaystyle \mathrm{PBR}=N_{\min }\frac{1}{2}{r}_{\max }{F}_{\mathrm{R}},}$ N_min being an estimate of minimum population size, $\frac{1}{2}{r}_{\max }$ one half of the maximum theoretical productivity rate of the stock and F_R a recovery factor chosen between 0.1 and 1. The computation of PBR does not require data on removals, although they are necessary to gauge the level of removals against the PBR limit (that is, to carry out an assessment). PBR is extensively used in the United States where it is enshrined in the US Marine Mammal Protection Act (MMPA). It has also been used in several cases, including dolphins, for example in New Zealand (Slooten & Dawson, 2008) or in the EU (e.g., Baltic harbour porpoise Berggren et al., 2002; grey seal Halichoerus grypus Taylor et al., 2022; common dolphin Delphinus delphis ICES, 2020).

Removals limit algorithm

The other harvest control rule currently in use is known as the Removals Limit Algorithm (RLA; Winship, 2009; Hammond, Paradinas & Smout, 2019). It also aims at setting an upper limit to anthropogenic mortality of PETS, and was applied to the population of Harbour porpoise in the North Sea (Hammond, Paradinas & Smout, 2019; Genu et al., 2021; Taylor et al., 2022). The RLA requires both a time-series of abundance/biomass estimates (whereas PBR only requires one such estimate) and a time-series of removals (whereas PBR requires none for its computation). RLA is a variant of the Catch Limit Algorithm (CLA), formulated by Cooke (1994) for baleen whales (Internation Whaling Commission, 2012). The RLA is an ad hoc algorithm that uses historical series of captures and estimates of abundance (Hammond, Paradinas & Smout, 2019): (1)${\displaystyle N_{t+1}=N_t+r{N}_t\left(1-{\left(\frac{N_t}{N_0}\right)}^2\right)-R_t,}$

where N_t and R_t are respectively the abundance/biomass and removals at time t. The computation of the RLA control rule for setting a removals limit (as a fraction of the best available abundance estimate) is: (2)${\displaystyle \mathrm{removals limit}=r\times \max \left(0,D_T-\mathrm{IPL}\right),}$

where T stands for the current year, D_T is the current depletion (that is, $D_T=\frac{N_T}{K}$, K being the carrying capacity) and IPL (Internal Protection Level) is the depletion threshold below which the limit is set to 0 (Punt, 1993). Both r and D_T are estimated from the model defined by Eq. (1) in a Bayesian framework (Hammond, Paradinas & Smout, 2019; Genu et al., 2021) and removals limit is computed from the joint posterior distribution of (r, D_T). A point estimate is used in practice by selecting a quantile of the posterior distribution of the removal limit to account for uncertainty.

Development of a stochastic SPM

Population models for PETS are often based on SPM (e.g., Snover, 2008; Punt, 2016) which are standard models of population dynamics in situations of strong uncertainty and low information. SPMs seek to encompass important population processes governing the dynamics of abundance change over time (Hilborn & Walters, 1992; Ver Hoef, 2006): ${\displaystyle \mathrm{next}\mathrm{abundance}=\mathrm{previous}\mathrm{abundance}+\mathrm{recruitment}+\mathrm{growth}-\mathrm{natural}\mathrm{mortality}-\mathrm{anthropogenic removals}}$

In particular, density-dependence in paramount to take into account (Punt, 2016): a convenient functional form is that described by Pella & Tomlinson (1969), assuming a first-order Markovian process on abundance: (3)${\displaystyle N_{t+1}=N_t+r^{\ast }\left(\frac{z+1}{z}\right)N_t\left(1-{\left(\frac{N_t}{K}\right)}^z\right)-R_t.}$

Setting z = 1 gives the well-known Schaefer model in fisheries (Schaefer, 1991), and N_t and R_t are respectively the abundance/biomass and removals at time t, K being the carrying capacity and r^∗ the growth rate at the Maximum Net Productivity Level (MNPL; r^∗ is also known as the Maximum Sustainable Yield Rate).

Incorporating environmental variability (the so-called process noise ɛ_t) in a multiplicative way in Eq. (3) yields (Polacheck, Hilborn & Punt, 1993; Millar & Meyer, 2000): (4)${\displaystyle N_{t+1}=\left\{N_t+r^{\ast }\left(\frac{z+1}{z}\right)N_t\left(1-{\left(\frac{N_t}{K}\right)}^z\right)-R_t\right\}{\varepsilon }_t.}$ This shift to a stochastic framework, assumed to be more realistic for describing changes in abundance, implies that the {N_t} states are random variables, similar to ɛ_t. At time t, the noise ɛ_t is classically assumed to be unbiased and homoskedastic (and therefore independent on t and past random variables (removals, abundances, noises): ${\displaystyle \mathbb{E}\left[{\varepsilon }_t\right]=1,}$ ${\displaystyle \mathbb{V}\left[{\varepsilon }_t\right]={\sigma }^2.}$

Assuming a simple relationship between removals R_t and abundance N_t: (5)${\displaystyle R_t=\rho N_t}$ where $\rho \in \left]\right.0,1\left[\right.$ is a time-invariant removal rate (Bousquet, Duchesne & Rivest, 2008; Bordet & Rivest, 2014), the removal process becomes: (6)${\displaystyle R_{t+1}=\left\{R_t+\frac{z+1}{z}{r}^{\ast }{R}_t\left(1-{\left(\frac{R_t}{K\rho }\right)}^z\right)-\rho R_t\right\}{\varepsilon }_t.}$

The set of parameters in Eq. (6) is defined as $\theta =\left\{K,\sigma ,\rho ,r^{\ast }\right\}$. Parameter z is usually fixed rather than estimated (Fletcher, 1978; Best & Punt, 2020). For example, setting z = 2.39 corresponds to a MNPL of 60% of K that is customarily assumed for marine mammals (Punt, 2016; Kanaji et al., 2021).

Equation (5) is a simplifying assumption that allows to link the abundance and removal processes. Stochasticity is introduced in Eq. (6) which allows to write a likelihood function (see below and Supplementary Information Appendix 1) for estimating a removal rate from data. In particular, it allows a time-series of past removals to be analyzed as the response variable where most current approaches treat these removal data as covariates (e.g., Kanaji et al., 2021). In effect, this choice means that removals are used as an index of abundance.

Reparametrization

Following the reparametrization of Bordet & Rivest (2014), let ${\displaystyle Z_t=\frac{N_t}{K}{\left(\frac{r^{\ast }\left(z+1\right)}{z-\rho z+r^{\ast }\left(z+1\right)}\right)}^{\frac{1}{z}}}$ which is well defined for z − ρz + r^∗(1 + z) > 0, that is ${\displaystyle \rho <1+r^{\ast }\frac{z+1}{z}.}$ With z > 0 this assumption is not restrictive since 0 < ρ < 1. Eqs. (4) and (6) can be re-arranged: ${\displaystyle Z_{t+1}=\left(1-\rho +r^{\ast }\frac{z+1}{z}\right)Z_t\left(1-Z_t^z\right){\varepsilon }_t,}$ ${\displaystyle R_{t+1}=\left(1-\rho +r^{\ast }\frac{z+1}{z}\right)R_t\left(1-{\left\{D\left(\theta \right)R_t\right\}}^z\right){\varepsilon }_t}$ with Z_t = D(θ)R_t and ${\displaystyle D\left(\theta \right)=\frac{1}{K\rho }{\left(\frac{r^{\ast }\left(z+1\right)}{z-\rho z+r^{\ast }\left(z+1\right)}\right)}^{\frac{1}{z}}.}$ Note that positive removals imply $R_t\le \frac{1}{D\left(\theta \right)}$. To simplify notations and to adopt a more conventional Markovian writing, we rewrite: (7)${\displaystyle R_t=g\left(R_{t-1},\theta \right){\varepsilon }_t}$ where (8)${\displaystyle g\left(R_{t-1},\theta \right)=\left(1-\rho +r^{\ast }\frac{z+1}{z}\right)R_{t-1}\left(1-{\left\{D\left(\theta \right)R_{t-1}\right\}}^z\right).}$ Note that g is neither linear nor log-linear. Appendix 1 in the Supplementary Information details the likelihoods for estimating θ from data. In the remainder, we denote ℓ({R_t}|θ) the likelihood associated with removals.

Initial conditions

For practical reasons, the initial depletion D₀ at the start of the time-series of removals, instead of K, is estimated. This choice is convenient as $D_0=\frac{N_0}{K}$. Initial abundance N₀ is typically unknown for PETS and the first abundance estimate available may not even match the start of the observed removals time-series. In that case, information on initial depletion D₀ may be elicited from expert knowledge (see e.g., Bockting, Radev & Brükner 2023) or historical data, and given a prior distribution: it may be easier to elicit a prior on depletion (a quantity expected to be bounded between 0 and 1) than on K directly. In practice, K is deduced from D₀ and the first observed abundance estimate that is available (e.g., Punt, 2019). The set of parameters to estimate is now: θ = {r^∗, σ, ρ, D₀}. The stochastic SPM defined above (Eqs. 6 & 8) was implemented in the probabilistic programming language Stan (see Appendix 2; Carpenter et al., 2017).

Anthropogenic removals threshold (ART)

The stochastic SPM developed above meshes removal and abundance processes by assuming a crude proportionality between the two, which is simplistic but aligns with stakeholder perceptions of variations in PETS by-catch for example (M Authier, pers. obs.). It is not a realistic model as removals depends strongly on fishing effort, which in turn is influenced by a myriad of factors. However, as Rademeyer, Plagányi & Butterworth (2007) noted, “[e]stimators based on simple population models have often been shown to perform as well or better than those based on more complex ones”. We are thus less interested in biological realism than in the performance of this model for management.³ The stochastic SPM assumes stationarity (that is time-invariance) in ρ, the removal rate, which is at odds with the very purpose of managing removals. If management is meant to be effective, it will take action to precisely change the level of removals. By definition, management aims at changing ρ over time so the model is clearly wrong once management is implemented. Before management is implemented, stationarity is an assumption as there is typically little knowledge or data are too noisy to test it. To obviate this issue while retaining a parsinomious model with four parameters, we used a statistical estimation paradigm based on the notion of weighted (or local) likelihood (Tibshirani & Hastie, 1987; Wang, 2006; Agostinelli & Greco, 2013). This approach allows us to progressively down-weigh data extending the furthest back in time (i.e., to consider that the policy relevance of statistical information provided by each datum in a sample depends on how far back in time the datum was collected). Applied only to the likelihood of removals, it is formalized as follows (see Appendix 1): the likelihood ℓ({R_t}|θ) is replaced by

(9)${\displaystyle {\ell }^{\mathit{w}\left(\eta \right)}\left(\left\{R_t\right\}|\theta \right)=\prod _{t=1}^T{\left(f\left(R_t|R_{t-1},\theta \right)\right)}^{w_t\left(\eta \right)}}$

where the weights w_t(η) have the sense of kernels, i.e., they provide a kernel-based representation of the score function ${\displaystyle s\left(\theta \right)={\nabla }_{\theta }\ell \left(\left\{R_t\right\}|\theta \right).}$ In choosing the weights w_t(η), such that the above representation has good convergence properties, the following properties must hold (Agostinelli & Greco, 2013): w_t(.) should be a bounded differentiable non-negative function of t that may depend on a parameter η which can be consistently estimated by $\hat{\eta }$, such that ${\displaystyle \begin{array}{ccc}\underset{R_t}{sup}|w\left(t,\hat{\eta }\right)-c|\hfill & \underset{p}{\overset{t\to \infty }{\to }}\hfill & 0\hfill \\ \text{almost surely}\hfill & \hfill & \hfill \end{array}}$ where c is a positive constant. The following choice (Gaussian kernel): ${\displaystyle w_t\left(\eta \right)=exp\left(-\frac{{\left(T-t\right)}^2}{2{\eta }^2}\right)}$ obeys these requirements (with c = 1), provided the scale parameter η be estimated by some rule that favors the information provided by the sample of most recent observations in the estimation of θ, trading bias for precision (Hu & Zidek, 2002).

In this study, η was chosen so that data older than 50 years contribute less than 0.05 to the likelihood during estimation (w_t = 0.05 for (t − 50)). This choice (η = 20.4) is arbitrary but was found to work well in practice (see Results). In addition, a sample size of 50 appears adequate for our purposes (although smaller sample size could work but this point requires additional simulations).

Capitalizing from the stochastic SPM (Eq. 9), we propose two candidate control rules (as a fraction of the best available abundance estimate) which we call Anthropogenic Removals Threshold (ART) to emphasize that the quantity derived from these control rules represents a threshold beyond which conservation objectives run a high risk of not being met.

The first candidate is simply the posterior mean of the quantity: (10)${\displaystyle {\mathrm{candidate}}_1=\rho \times F_{\mathrm{R}}}$

where F_R is a recovery factor chosen between 0.1 and 1 (as in the PBR control rule). This rule adapts the historical removal rate and does not directly rely on an estimate of carrying capacity or population growth rate as the RLA control rule (although both a carrying capacity and a population growth rate are in θ). The second candidate is an elaboration that takes into account any trend in abundance, and more specifically any evidence of a decline, to negatively feedback on the removals limit: (11)${\displaystyle {\mathrm{candidate}}_2=\rho \times F_{\mathrm{R}}\times \min \left(1,exp\left(\beta \right)\right)}$

where β is a trend in abundance (on a logarithmic scale) estimated using the regularized approach detailed in Authier et al. (2020). Briefly, β is the slope of a regression line through the abundance estimates (scaled by the first estimate and then log-transformed) and estimated using a weakly-informative prior (namely the ‘skeptical’ prior of Cook, Fúquene & Pericchi, 2011) that favours the hypothesis of no trend over time (see Authier et al., 2020 for an in-depth investigation of this choice in the context of trend estimation). This candidate rule operationalizes the principle on non-deterioration whereby populations or species in need of restoration (that is, that are depleted) should not be allowed to deteriorate further. If no trend or a positive trend in abundance is evidenced, candidate₂ is equivalent to candidate₁. Both candidate₁ and candidate₂ can be computed for the same data necessary to compute RLA, and their posterior mean approximated by the average over MCMC samples of the posterior distribution after estimating θ from data (see Appendices 2 and 3).

(‘Harvest’) control rules

We tested five control rules for managing anthropogenic removals of PETS:

• the fixed percentage rule of ASCOBANS (2000): $\mathrm{ASCOBANS}=0.017\times N_T^{\mathrm{obs}}$;

• the PBR rule of Wade (1998) with N_min defined as the 20% quantile of a log-normally distributed abundance estimate $N_T^{\mathrm{obs}}$: $\mathrm{PBR}=N_{\min }\frac{1}{2}{r}_{\max }{F}_{\mathrm{R}}$;

• the RLA rule, with $\mathrm{RLA}=N_T^{\mathrm{obs}}\times \mathrm{removals limit}=N_T^{\mathrm{obs}}\times r\times \max \left(0,D_T-\mathrm{IPL}\right)$;

• the candidate₁ ART, with ${\mathrm{ART}}_1=N_T^{\mathrm{obs}}\times {\mathrm{candidate}}_1=N_T^{\mathrm{obs}}\times \rho \times F_{\mathrm{R}}$; and

• the candidate₂ ART, with ${\mathrm{ART}}_2=N_T^{\mathrm{obs}}\times {\mathrm{candidate}}_2=N_T^{\mathrm{obs}}\times \rho \times F_{\mathrm{R}}\times \min \left(1,exp\left(\beta \right)\right)$.

To set a removals limit, the fixed percentage rule of ASCOBANS (2000) only requires a point estimate of current abundance $N_T^{\mathrm{obs}}.\mathrm{PBR}$ requires a recent point estimate of abundance and its associated coefficient of variation to compute an estimate of minimum population size N_min. The rules RLA and ART_i∈[1:2] each require a time-series of abundance estimates with their uncertainties as well as a time-series of removals’ estimates. All rules, save for the fixed percentage one (which was included as a negative control), need tuning. For PBR, ART₁ and ART₂, this process means the testing of different values of F_R to identify the minimum one that allows to reach a given CO. For RLA, tuning is achieved by testing different quantiles of the posterior distribution of Eq. (2) (Genu et al., 2021). That quantile tuning was not carried out with ART₁ or ART₂ stemmed from the typically tight posterior concentration observed when estimating ρ during the development of the stochastic SPM (Ouzoulias, 2022). In contrast, posterior concentration does not occur because the log-normal likelihood assumed for Eq. (1) when estimating removals limit is down-weighted by a fixed factor $\frac{1}{16}$ to limit the speed at which the management procedure responds to feedback (Cooke, 1999).

Scenarios

We evaluated control rules on three scenarios: a base case scenario whereby unbiased but noisy data are assumed to be available and collected; and two so-called robustness trials. In the first robustness trial, estimates of abundance were assumed to be observed with a systematic bias resulting in an overestimation by a factor 2. In the second, removals’ estimates were assumed to be biased downward, resulting in an underestimation of true removals by a factor 2. The two robustness trials were found to be the most challenging ones in a previous investigation (Genu et al., 2021).

Performance metrics

The effectiveness of control rules to meet requirements of management was gauged on 6 performance metrics. A period of 6 years was chosen to match the reporting cycle for the MSFD. Accordingly, it was assumed during simulations that abundance estimates would be updated every 6 years. The first performance metric allows tuning by selecting either a value for F_R or a quantile for each control rule in order to meet a CO. The second allowed a comparison with the fixed percentage rule of ASCOBANS. The third and fourth metrics informed on the variability in removals limits for a given control rule, with greater stability with low cv and positive autocorrelation. The fifth indicated any systematic changes in removals limits over time: it is the estimated slope of a linear regression of the computed removals limit against time with time rescaled to the unit interval (the model was fitted using the default settings of function stats::lm in R version 4.2.2) (R Core Team, 2022). This performance metrics reveal how the removals limit is changing on average over time (increasing, decreasing, no change) as it is updated every 6 years in light of monitoring data. The sixth performance metric informed on the expected delay in population recovery, where recovery was defined as the time at which a CO was sustained over 6 consecutive years (to match the MSFD reporting cycle). The performance metrics were:

1. the level of population (true) depletion at the end of each simulation;

2. the average removals limit (average, as a percent of current abundance estimate) across simulations;

3. the coefficient of variation (cv) of removals limits computed within each simulation and averaged over simulations;

4. the 1-lagged autocorrelation (autocorrelation) in removals limits computed within each simulation and averaged over simulations;

5. the temporal trend (trend, in %) in removals limits within a simulation and averaged over simulations; and

6. the average delay in recovery (delay, in years) compared to a hypothetical complete ban on removals.

Implementation

To carry out a (desk) MSE, we used a cross-factorial design with two CO, five control rules, three scenarios, 10 F_R and 100 simulations. The total computational budget was 30, 000 runs which were sent to the computing cluster Curta hosted at Bordeaux University, France. Computer time for this study was provided by the computing facilities of the MCIA (“Mésocentre de Calcul Intensif Aquitain”). All simulations were carried out using the functionalities of the dedicated R package RLA v.0.2.0 (Genu et al., 2021) available at https://gitlab.univ-lr.fr/pelaverse/RLA. The different functions used are mentioned in Table 3. As simulations rely on random numbers, seeds were controlled for and recycled across control rules within scenarios in order to match and compare results (Rademeyer, Plagányi & Butterworth, 2007).

Computation of PBR does not require any model fitting in contrast to RLA and all flavours of ART which have an estimation step for quantities of interest to compute removals limits. With a time horizon of 100 years and an assumed frequency of 6 years (as per MSFD cycle), $\lceil \frac{100}{6}\rceil =17$ model fittings are required within one run. Estimation of θ (Appendices 2 and 4) was carried out in a Bayesian framework using programming language Stan (Carpenter et al., 2017) called from R v.4.2.2 (R Core Team, 2022) with library Rstan (Stan Development Team, 2023). Four chains were initialized and run for a total of 2,000 iterations, discarding the first 1,000 as warm-up and retaining every 1 iterations (no thinning). Control parameters for the No-U-Turns algorithm were kept to their default values. Parameter convergence was assessed using the $\hat{R}$ statistics (Vehtari et al., 2019) and assumed if $\hat{R}<1.025$.