## Change in stock-recruit relationship

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) { 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/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)[] <- 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)}$