## Size spectrum dynamics

In previous tutorials we have concentrated on the steady state of the mizer model, where for each size class and each species, the rate at which individuals grow into the size class balances the rate at which individuals grow out of the size class or die, thus keeping the size spectrum constant. In this tutorial we explore the dynamic that takes place when this balance is changed.

Size-spectrum dynamics is described by the beautiful partial differential equation

\frac{\partial N(w)}{\partial t} + \frac{\partial g(w) N(w)}{\partial w} = -\mu(w) N(w)

together with the boundary condition

N(w_{min}) = \frac{R_{dd}}{g(w_{min})},

where N(w) is the number density at size w, g(w) is the growth rate and \mu(w) is the death rate of individuals of size w, w_{min} is the egg size and R_{dd} is the birth rate. Luckily it is easy to describe in words what these equations are saying.

Important

Size spectrum dynamics is very intuitive: The rate at which the number of individuals in a size class changes is the difference between the rate at which individuals are entering the size class and the rate at which they are leaving the size class. Individuals can enter a size class by growing in or, in the case of the smallest size class, by being born into it. They can leave by growing out or by dying.

What makes these seemingly obvious dynamics interesting is how the growth rate and the death rate are determined in terms of the abundance of prey and predators and the feeding preferences and physiological parameters of the individuals. We have discussed a bit of that in previous tutorials and will discuss it much more in upcoming tutorials. We will discuss the birth rate R_{dd} below in the section on how reproduction is modelled. But first we want to look at the results of simulating the size spectrum dynamics.

You will be doing your own work for this tutorial in the accompanying worksheet file “worksheet5-dynamics-of-spectra.Rmd”. You will find this file in the worksheet repository for this week that you already used in previous tutorials.

As always, we start by making sure we have the latest version of the mizerExperimental package and load it as well as the tidyverse.

remotes::install_github("sizespectrum/mizerExperimental")
library(mizerExperimental)
library(tidyverse)

### Projections

In the previous tutorial, in the section on trophic cascades, we already simulated the size-spectrum dynamics to find the new steady state. But we only looked at the final outcome once the dynamics had settled down to the new steady state. We reproduce the code here:

# Create trait-based model
mp <- newTraitParams() |>
# run to steady state with constant reproduction rate
# turn of reproduction and instead keep egg abundance constant
setRateFunction("RDI", "constantEggRDI") |>
setRateFunction("RDD", "noRDD")
The resource carrying capacity has been commented and therefore will not be recalculated from the resource parameters.
Convergence was achieved in 6 years.
The resource carrying capacity has been commented and therefore will not be recalculated from the resource parameters.
# We make a copy of the model
mp_lessRes <- mp
# and set the resource interaction to 0.8 for species 8 to 11
species_params(mp_lessRes)\$interaction_resource[8:11] <- 0.8

# We run the dynamics until we reach steady state
mp_lessRes_steady <- projectToSteady(mp_lessRes)
Convergence was achieved in 24 years.
# We compare the steady states
ylim = c(1e-8, NA), wlim = c(1e-3, NA)) But we can also save and then display the spectra of all the species at intermediate times. This is what the project() function does. It projects the current state forward in time and saves the result of that simulation in a larger object, a MizerSim object, which contains the resulting time series of size spectra. Let’s use it to project forward by 24 years.
sim_lessRes <- project(mp_lessRes, t_max = 24)
We can now use this MizerSim object in the animateSpectra() function to create an animation showing the change in the size spectra over time.
animateSpectra(sim_lessRes, total = TRUE, power = 2,
ylim = c(1e-8, NA), wlim = c(1e-3, NA))