Build model

You have some life-history parameters of some important species in your ecosystem.

sub_rakhine <- readRDS("data/sub_rakhine.rds")
sub_rakhine

The variables we will be interested in are the species name, w_inf the asymptotic size, w_mat the maturity size, beta and sigma the parameters of the log-normal predation kernel, a and b the parameter in the length to weight conversion \(w = a\, l^b\), and k_vb the von Bertalanffy growth curve parameter.

species_params <- sub_rakhine %>% 
    select(commonName, w_inf, w_mat, beta, sigma, w_a, w_b, k_vb) %>% 
    rename(a = w_a, b = w_b, species = commonName)
species_params

For later purposes we also set up some fishing gear. We give one gear to each species and set knife_edge selectivity with a knife_edge_size equal to the maturity size.

species_params <- species_params %>% 
    mutate(gear = species,
           sel_func = "knife_edge",
           knife_edge_size = w_mat)

We want to study these species embedded in a set of background species creating a complete balanced ecosystem. We use a trait-based model for the background species.

params <- newTraitParams(
    min_w_pp = 10^-12,
    no_sp = 16, 
    min_w_inf = 1, 
    max_w_inf = 10^6,
    min_egg = 10^-5,
    min_w_mat = 10^-0.6,
    w_pp_cutoff = 2 * 10^-0.6,
    no_w = 200,
    sigma = 2,
    ks = 2.4,
    bmort_prop = 0.4) %>% 
    steady(t_max = 100) %>% 
    markBackground()

Note: Using f0, h, lambda, kappa and the predation kernel to calculate gamma.
Note: Negative plankton abundances
Note: Negative plankton abundance values overwritten with zeros
Steady state was reached before 30 years.

Note how we used the pipe operator %>% to pipe the result of newTraitParams() into steady() and then the result from that into markBackground(). The mizer functions are pipe-friendly in that their first parameter is always the MizerParams object or the MizerSim object. The call to steady() ran the dynamics with fixed reproductive rate to find the steady state abundances. The call to markBackground() just specifies that these are background species. This has the effect for example that they are grey when plotted:

plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

and we will see other consequences below.

We don’t want to fish the background species, so we set their catchability to zero.

params@species_params$catchability[] <- 0

We now add our foreground species with the parameters specified in the data frame species_params and again run to steady state.

params <- addSpecies(params, species_params, initial_effort = 1) %>% 
    steady()

Note: No h provided for some species, so using f0 and k_vb to calculate it.
Note: Using f0, h, lambda, kappa and the predation kernel to calculate gamma.
Steady state was reached before 45 years.

plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

We see that addSpecies() adds the new species at a very low abundance compared to the background species. This is on purpose, because at that low abundance there is very little backreaction from the background species, and therefore we are in a regime where we understand the steady-state size distribution of the species analytically. Therefore mizer was able to find initial values for the size distributions that are close to steady state values and so the call to steady() did not have much work to do.

We make a copy of this model for later purposes.

params2 <- params

Now we can increase the abundance of the foreground species.

params <- rescaleAbundance(params, factor = 16)
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

Note how rescaleAbundance() knew to only rescale the foreground species and leave the background species alone. Note also, that the introduction of the foreground species has led to a bulge in the total community spectrum. However the abundance of the background species was originally chosen so that the community abundance would be close to the Sheldon power-law spectrum. We therefore retune their abundance so that again the community spectrum is as close to the Sheldon power-law as possible.

params <- retuneBackground(params)
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

We see that as we introduce the foreground species, the background species become less important.

The above plot does not yet show a steady state. We still need to call steady() to find the steady state.

params <- steady(params)

Steady state was reached before 90 years.

plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

Clearly there is some tuning to be done. In particular, the abundances of all species need to be increased and tuned to the observed values and the fishing parameters need to be tuned to produce the observed catches.

Now I want to introduce a very convenient gadget that allows you to tune the parameters of the system interactively. This is not yet very well documented, so I will demonstrate when we have our video chat.

params <- tuneParamsCourse(params)