While the stock-recruit relationship is constant in `openMSE`

, it is possible to change the effective stock-recruit relationship by changing the mean \(\mu\) of the recruitment deviations through `OM@cpars$Perr_y`

. This article describes how reference points implicitly change with \(\mu\). The internal code used to calculate annual reference points in `MSEtool::runMSE`

is not altered when the mean of `OM@cpars$Perr_y`

is altered, so it is up to the user to decide whether reference points (see code below) need to be re-calculated.

## Stock-recruit parameters

With the Beverton-Holt stock recruit relationship, the recruitment in year \(y\) is:

\[ R_y = \dfrac{\alpha SB_y}{1 + \beta SB_y} \delta_y\]

where \(\delta_y\) are recruitment deviates. For stationary productivity around the recruitment predicted by the S-R curve, we sample \(\delta_y\) from a lognormal distribution with \(E[\delta_y] = 1\) and \(V[\delta_y] = \exp(\sigma^2-1)\). These values are passed to the operating model in `OM@cpars$Perr_y`

. Parameters \(\alpha\) and \(\beta\), along with \(\phi_0\), can be used to calculate the corresponding unfished recruitment \(R_0\) and steepness \(h\):

\[ h = \dfrac{\alpha \phi_0}{4 + \alpha \phi_0}\]

\[ R_0 = \dfrac{1}{\beta}\left(\alpha - \dfrac{1}{\phi_0}\right)\]

We define the unfished biomass as \(SB_0 = R_0 \phi_0\). For easier reading, we have dropped the superscript \(\textrm{SR}\) associated with \(R_0\), \(h\), and \(\phi_0\).

If we want to change the productivity via recruitment deviations, we can re-scale the recruitment deviations. Set \(\delta'_y = \mu \delta_y\), with \(\mu < 1\) for a less productive system. In effect, we will have created a different stock-recruit relationship:

\[ R_y = \dfrac{\alpha SB_y}{1 + \beta SB_y} \delta'_y = \dfrac{\alpha' SB_y}{1 + \beta SB_y} \delta_y\]

where \(\alpha' = \mu \alpha\), \(E[\delta'_y] = \mu\) and \(V[\delta'_y] = \exp(\sigma^2-1)/\mu^2\), and \(\beta\) is unchanged. In other words, by setting the mean of `OM@cpars$Perr_y`

to \(\mu\), we create a stock-recruit relationship defined by \(\alpha'\) and recruitment variance of \(V[\delta'_y]\).

## Implied reference points

Annual unfished and MSY reference points can be re-calculated accordingly with \(\alpha'\), \(\beta\), and \(\phi_0\). The implications on the reference points can be seen in “new” parameters \(R_0'\) and \(h'\) corresponding to “new” \(\alpha'\):

\[ h' = \dfrac{\alpha' \phi_0}{4 + \alpha' \phi_0} =\dfrac{\mu h}{1 + h (\mu - 1)}\]

\[ R_0' = \dfrac{1}{\beta}\left(\alpha' - \dfrac{1}{\phi_0}\right) = R_0\dfrac{h(4\mu + 1) - 1}{5h-1}\] Then, \(SB_0' = R_0' \phi_0\).

In the simplest case with constant \(\phi_0\), we expect \(h' < h\) and \(R_0' < R_0\) when \(\mu < 1\), and also expect \(SB_0\), \(F_{\textrm{MSY}}\), \(\textrm{MSY}\), and \(SB_\textrm{MSY}\) to decrease.

The new \(SPR_{\textrm{crash}} = (\alpha'\phi_0)^{-1} = (\mu \alpha \phi_0)^{-1}\) increases (\(F_{\textrm{crash}}\) decreases) when \(\mu < 1\).

## Examples

### Contour plots

Let’s explore the ratio of the new and old values of \(R_0\) and \(h\) as a function of \(\mu\) and the original steepness, with \(\beta = 1\) and \(\phi_0 = 1\):

Interestingly, there can exist a combination of low steepness and low \(\mu\) such that the population crashes, as implied by \(R_0' \le 0\) in red.

### Forward projections

Let’s explore the dynamics in a forward projecting model in `MSEtool`

. We’ll create an operating model with 4 simulations with varying \(\mu\) but constant stock-recruit parameters among simulations. Then we’ll project the model using the `NFref`

(no fishing) MP.

`library(MSEtool)`

`## Loading required package: snowfall`

`## Loading required package: snow`

```
OM <- MSEtool::testOM
# Generate four simulations with same alpha, beta, and phi_0
OM@nsim <- 4
OM@Linfsd <- OM@Msd <- OM@Ksd <- c(0, 0)
OM@Linf <- mean(OM@Linf) |> rep(2)
OM@K <- mean(OM@K) |> rep(2)
OM@t0 <- mean(OM@t0) |> rep(2)
OM@M <- mean(OM@M) |> rep(2)
OM@L50 <- mean(OM@L50) |> rep(2)
OM@L50_95 <- mean(OM@L50_95) |> rep(2)
OM@h <- mean(OM@h) |> rep(2)
OM@Prob_staying <- OM@Frac_area_1 <- OM@Size_area_1 <- c(0.5, 0.5)
# With four simulations, let's set the projection Perr_y to these values:
mu_vec <- c(0.01, 0.75, 1, 1.5)
OM@cpars$Perr_y <- matrix(1, OM@nsim, OM@maxage + OM@nyears + OM@proyears)
OM@cpars$Perr_y[, OM@maxage + OM@nyears + 1:OM@proyears] <- mu_vec
Hist <- Simulate(OM, silent = TRUE)
MSE <- Project(Hist, MPs = "NFref", silent = TRUE)
```

Let’s see what happens with the `NFref`

MP:

Despite the same stock-recruit parameters, the spawning biomass all diverge with \(\mu\):

This terminal SB should be the new \(SB_0\), so let’s create a function to calculate these values analytically instead of projecting:

```
recalc_ref_pt <- function(Hist, mu = rep(1, Hist@OM@nsim)) {
StockPars <- Hist@SampPars$Stock
FleetPars <- Hist@SampPars$Fleet
R0 <- StockPars$R0
hs <- StockPars$hs
if(StockPars$SRrel[1] == 1) {
R0_new <- R0 * (hs * (4 * mu + 1) - 1)/(5 * hs - 1)
hs_new <- mu * hs /(1 + hs * (mu - 1))
} else {
R0_new <- 1 + 0.8 * R0 * log(mu)/log(5 * hs)
hs_new <- hs * mu^0.8
}
MSY_ref_pt <- lapply(1:Hist@OM@nsim, function(x) { # Internal function for calculating unfished and MSY reference points
sapply(Hist@OM@nyears + 1:Hist@OM@proyears, function(y) {
MSEtool:::optMSY_eq(x,
M_ageArray=StockPars$M_ageArray,
Wt_age=StockPars$Wt_age,
Mat_age=StockPars$Mat_age,
Fec_age=StockPars$Fec_Age,
V=FleetPars$V_real,
maxage=StockPars$maxage,
R0=R0_new,
SRrel=StockPars$SRrel,
hs=hs_new,
SSBpR=StockPars$SSBpR,
yr.ind=y,
plusgroup=StockPars$plusgroup)
})
})
crash_ref_pt <- local({ # Internal function for calculating SPR crash
boundsF <- c(1E-3, 3)
F_search <- exp(seq(log(min(boundsF)), log(max(boundsF)), length.out = 50))
lapply(1:Hist@OM@nsim, function(x) {
sapply(Hist@OM@nyears + 1:Hist@OM@proyears, function(y) {
Ref_search <- MSEtool:::Ref_int_cpp(F_search,
M_at_Age = StockPars$M_ageArray[x, , y],
Wt_at_Age = StockPars$Wt_age[x, , y],
Mat_at_Age = StockPars$Mat_age[x, , y],
Fec_at_Age = StockPars$Fec_Age[x, , y],
V_at_Age = FleetPars$V_real[x, , y],
maxage = StockPars$maxage,
plusgroup = StockPars$plusgroup)
SPR_search <- Ref_search[2, ]
RPS <- Ref_search[3, ]
if(StockPars$SRrel[x] == 1) {
CR <- 4 * StockPars$hs[x]/(1 - StockPars$hs[x])
} else if(StockPars$SRrel[x] == 2) {
CR <- (5 * StockPars$hs[x])^1.25
}
alpha <- mu[x] * CR/StockPars$SSBpR[x, 1] # New alpha
if(min(RPS) >= alpha) { # Unfished RPS is steeper than alpha
SPRcrash <- min(1, RPS[1]/alpha) # Should be 1
Fcrash <- 0
} else {
SPRcrash <- MSEtool:::LinInterp_cpp(RPS, SPR_search, xlev = alpha)
Fcrash <- MSEtool:::LinInterp_cpp(RPS, F_search, xlev = alpha)
}
c(SPRcrash = SPRcrash, Fcrash = Fcrash)
})
})
})
ref_pt <- abind::abind(simplify2array(MSY_ref_pt), simplify2array(crash_ref_pt), along = 1)
dimnames(ref_pt)[[1]] <- c("MSY", "FMSY", "SBMSY", "SBMSY_SB0", "BMSY_B0", "BMSY", "VBMSY", "VBMSY_VB0", "RMSY", "SB0", "B0", "R0", "h",
"N0", "SN0", "SPRcrash", "Fcrash")
return(ref_pt)
}
```

Let’s compare the ratio of the calculated \(SB_0\) and the projected terminal SSB:

```
new_ref_pt <- recalc_ref_pt(Hist, mu = mu_vec)
SSB0_prime <- new_ref_pt["SB0", OM@proyears, ]
plot(mu_vec, SSB0_prime/terminal_SB, xlab = expression(mu), ylab = "SB0/Terminal SB")
```

Looks like our new reference points match the values from the projections. Note with \(\mu = 0.01\), the population crashes so \(SB_0 = 0\).

Now let’s map out \(h'/h\) and \(R_0'/R_0\) as a function of \(\mu\):

If \(\mu < 0.08\), our population crashes because \(R_0' < 0\), and \(SPR_{\textrm{crash}} = 1\):

```
SPRcrash <- new_ref_pt["SPRcrash", OM@proyears, ]
plot(mu_vec, SPRcrash, xlab = expression(mu), ylim = c(0, 1))
```

## Ricker function

With a Ricker function where \(R_y = \alpha SB_y \exp(-\beta SB_y) \delta_y\), steepness is \(h = 0.2(\alpha\phi_0)^{0.8}\) and \(R_0 = \log(\alpha\phi_0)/(\beta \phi_0)\).

The new SR relationship is \(R_y = \alpha' SB_y \exp(-\beta SB_y) \delta_y\), with

\[h' = 0.2 (\alpha' \phi_0)^{0.8} = \mu^{0.8} h\]

\[R_0' = 1 + \dfrac{0.8 R_0\log(\mu)}{\log(5h)}\]