The simplest version of the size spectrum model is the community model. In the community model, individuals are only characterized by their size and are represented by a single group representing an across species average. Reproduction is not considered and the recruitment flux \(R\) is set to be constant. The resource spectrum only extends to the start of the community spectrum. Standard metabolism is turned off.
In this section we describe how a community model can be set up and projected through time. We then use a community model to illustrate the idea of a ’trophic cascade’. Due to the relative simplicity of this type of model they are useful for gently introducing some of the concepts behind the mizer package. Consequently, this section should hopefully serve as an introduction to using mizer and some of the main classes and methods.
Setting up a community model
The first stage in implementing a model using mizer is to create an object of class MizerParams. This class contains the model parameters including the life-history parameters of the species in the model, the stock-recruitment relationships, the fishing selectivity functions and the parameters of the plankton spectrum.
To avoid having to make a MizerParams object directly, the set_community_model() wrapper function, has been provided that conveniently creates a MizerParams object specifically for a community model. The documentation for the function can be seen by entering:
As can be seen in the help page, the function can take many arguments. We can ignore most of these for the moment as they almost all come with default values.
The arguments that you should pay attention to are: z0 (the level of background mortality), alpha (the assimilation efficiency of the community), f0 (the average feeding level of the community which is used to calculate \(\gamma,\) as described in the equation for gamma and recruitment (the level of constant recruitment).
Although default values for these parameters are provided, you are encouraged to explore how changing the values affects the simulated community. For example, the default value of z0 is 0.1. Increasing this value effectively ‘shortens’ the length of the community spectrum. The value of the constant recruitment should be set so that the community spectrum is a continuation of the plankton spectrum. This can be done with trial and error. A reasonable value for alpha is 0.2.
The set_community_model() function is called by passing in the arguments by name. Any parameter that is not passed in is set to the default value. For example, the following line sets up the parameters with z0 = 0.1, f0 = 0.7, alpha = 0.2 and recruitment = 4e7. All other parameters will have their default value:
Calling the function creates and returns an object of type MizerParams. We can check this using the class() function.
## [1] "MizerParams"
## attr(,"package")
## [1] "mizer"If you are going through this vignette for the first time, it is likely that you have no idea what a MizerParams object actually is. In R-speak it is an S4 object, which means it is an object made up of slots. Slots are essentially containers that store the object data. The names of these slots can be seen by calling the slotNames() method on the object:
## [1] "w" "dw" "w_full"
## [4] "dw_full" "psi" "initial_n"
## [7] "intake_max" "search_vol" "activity"
## [10] "std_metab" "ft_pred_kernel_e" "ft_pred_kernel_p"
## [13] "mu_b" "rr_pp" "cc_pp"
## [16] "sc" "initial_n_pp" "species_params"
## [19] "interaction" "srr" "selectivity"
## [22] "catchability" "n" "p"
## [25] "lambda" "q" "f0"
## [28] "kappa" "A" "linecolour"
## [31] "linetype"As you can see, the params object is made up of lots of slots. A quick description can be found in the MizerParams help page (run class ? MizerParams. The slots of an object are accessed by using the @ operator. For example, to access the w slot (which contains a vector of the size bins in the model) you would use
Rather than picking through the slots to find out what is in a particular MizerParams object, a summary of the object can be seen by calling the summary() method on it:
## An object of class "MizerParams"
## Community size spectrum:
## minimum size: 0.001
## maximum size: 1e+06
## no. size bins: 100
## Background size spectrum:
## minimum size: 1e-10
## maximum size: 1e+06
## no. size bins: 177
## Species details:
## species w_inf w_mat beta sigma
## Community Community 9e+05 NA 100 2
## Fishing gear details:
## Gear Target species
## Community CommunityIn the summary you can see that the size range of the community spectrum has been set from \(0.001\) to \(10^{6}\) and these are spread over \(100\) size bins. Similar information is available for the plankton spectrum. Additionally, the community is made up of only one species, called Community, which has an asymptotic size of \(9\times 10^{5}\) and a preferred predator prey mass ratio of \(100\). The w_mat parameter has been set to NA as it is not used when running a community model. These values have all been set by default using the set_community_model() function. If you want to set different values for these, you will need to call the set_community_model() function and pass in the desired argument values.
Running the community model
By using the set_community_model() method we now have a MizerParams object that contains all the information we need about the model community. We can use this to perform a simulation and project the community through time. In the mizer package, projections are performed using the project() method. You can see the help page for project() for more details and it is described fully in the section on running a simulation. We will ignore the details for the moment and just use project() to run some simple projections. The arguments for project() that we need to be concerned with are effort, which determines the fishing effort (and therefore fishing mortality) through time, and t_max, which is the length of the simulation. Initial population abundances are set automatically by the get_initial_n() function. It is possible to set your own initial abundances but we will not do this here.
To run a projection for 50 time steps, with no fishing effort (i.e. we want to model an unexploited community) we run:
The resulting object, sim, is of type MizerSim. This class holds the results of the simulation, including the community and plankton abundances at size through time, as well as the original model parameters. It is explained in detail in the section on running a simulation.
After running the projection, it is possible to explore the results using a range of plots and analyses. These are described fully in the section on exploring the simulation results.
To quickly look at the results of the projection you can call the generic plot() method. This plots the feeding level, predation mortality, fishing mortality and abundance by size in the last time step of the simulation, and the biomass through time. Each of the plots can be show individually if desired.
Example plot of the community model.
In the above plot there are several things going on that are worth talking about. Looking at the total biomass of the community against time, you can see that the biomass quickly reaches a stable equilibrium. The other panels show what is happening at the last time step in the simulation, which in this case is when the community is at equilibrium. Fishing mortality is 0 because we set the effort argument to 0 when running the simulation. The predation mortality rate is clearly a function of size, with the smallest sizes experiencing the highest levels of predation. The feeding level describes how satiated an individual is, with 0 being unfed, and 1 being full satiated. The feeding level at size will be strongly affected by the values of the f0 and alpha arguments passed to the set_community_model() function. The plankton and community spectra are shown in the bottom panel of the plot (the plotted plankton spectrum has been truncated to make for a better plot, but really extends all the way back to \(10^{-10}\) g). You can see that the community spectrum forms a continuum with the plankton spectrum. This is strongly affected by the level of fixed recruitment (the recruitment argument passed to set_community_model()).
Note the hump in the biomass at the largest end of the community spectrum. This is because the size spectrum model can be broadly described as ‘big things eating little things’. Given this, what is eating the very biggest things? Without fishing pressure, the mortality of the largest individuals is only from the background mortality (determined by the z0 argument) and the mortality from predation is almost 0. This is difficult to see in the plot due to the predation mortality being so high for the smaller individuals.
We can see this more clearly by extracting the predation mortality information from the MizerSim object, sim, that we created above. This is easily done by using the getM2() function (see the help page for more details). There are several methods that can be used for extracting information from a MizerSim object, e.g. getFeedingLevel() and getFMort(). For more information see the section on exploring the simulation results. Here we just call getM2() using the sim object:
The resulting m2 object is an array that contains the predation mortality at time by species by size. Here we only have one species so the species dimension is dropped, leaving us with a two dimensional array of time by size. We projected the model for \(50\) time steps but the length of the time dimension is \(51\) as the initial population is also included as a time step.
To pull out the predation mortality at size in the final time step we use:
If you plot this predation mortality on a log-log scale you can see how the predation mortality declines to almost zero for the largest sizes.
plot(x = sim@params@w, y = m2[51,], log = "xy", type = "l",
xlab = "Size", ylab = "Predation mortality")
Predation mortality without fishing in the community model (note the log scales for both axes).
Example of a trophic cascade with the community model
It is possible to use the community model to simulate a trophic cascade. To do this we need to perform two simulations, one with fishing and one without.
This means we need to consider how fishing is handled in mizer. The set_community_model() function automatically sets the fishing selectivity to have a knife-edge shape, with only individuals larger than 1 kg selected (the size at the knife-edge can be changed by setting the knife_edge_size argument). Although it is possible to change the selectivity function, here we will use the default knife-edge selectivity. We set up the parameter object exactly as before, but now we are explicitly setting the size at which individuals are selected by the fishing gear.
params_knife <- set_community_model(z0 = 0.1, recruitment = 4e7,
alpha = 0.2, f0 = 0.7, knife_edge_size = 1000)First we perform a simulation without fishing in the same way we did above by setting the effort argument to 0:
Now we want to simulate again, this time with fishing. In the simulations, fishing mortality is calculated as the product of the fishing selectivity, effort and catchability (see the section on fishing gears for more details). By default catchability is set to 1. This means that a fishing effort of 1 will result in a fishing mortality of 1 for fully selected sizes. Here we run a simulation with fishing effort set to 1 for the duration of the simulation:
You can compare the difference between these scenarios by using the plot() method as before. Of particular interest is the fishing mortality at size. The knife-edge selectivity at 1000 g can be clearly seen and an effort of 1 has resulted in a fishing mortality of 1 for the fully selected sizes.
Summary plot for the community model when fishing with knife-edge selectivity at size = 1000 g.
To explore the presence of a trophic cascade, we are interested in looking at the relative change in abundance when the community is fished compared to when it is not fished. To do this we need to get the abundances at size from the simulation objects. The abundances are stored in the n slot of the MizerSim objects. The n slot contains a three dimensional array with dimensions time x species x size. Here we have 51 time steps (50 from the simulation plus one which stores the initial population), 1 species and 100 sizes:
## [1] 51 1 100We want the abundances in the final time step, and we can use these to calculate the relative abundances:
This can then be plotted using basic R plotting commands. The sizes are stored in the params@w slot (a slot of a slot!).
plot(x = sim0@params@w, y = relative_abundance, log = "x", type = "n",
xlab = "Size (g)", ylab = "Relative abundance")
lines(x = sim0@params@w, y = relative_abundance)
lines(x = c(min(sim0@params@w), max(sim0@params@w)), y = c(1,1), lty = 2)
Relative abundances from the unfished (dashed line) and fished (solid line) community model.
The impact of fishing on species larger than 1000g can be clearly seen. The fishing pressure lowers the abundance of large fish (the decrease in relative abundance at 1000 g). This then relieves the predation pressure on their smaller prey (the preferred predator-prey size ratio is given by the \(\beta\) parameter, which is set to 100 by default), leading to an increase in their abundance. This in turn increases the predation mortality on their smaller prey, which reduces their abundance and so on.
The impact of changing \(\sigma\)
As described above, the \(\sigma\) parameter determines the width of the predator prey size preference. Here we take a look at how changing the value of \(\sigma\) can affect the dynamics of the community. In the examples above, \(\sigma\) is set in the set_community_model() function by default to a value of \(2\). When projected through time, the community abundances converge to a stable equilibrium. What happens if we reduce the value of \(\sigma\), for example by setting it to 1.0? We can do this by passing in the new value of \(\sigma\) into set_community_model().
We want to project this new model through time using the project() method. Note that we have introduced a new argument: \(dt\). This is the step size of the solver. It does not have anything to do with the biology in the model. It only affects the internal engine of project() that performs the projection.
As you can see in the underlying model equations in the model description section, the model runs in continuous time. Therefore, to project it forward, project() must solve the system of equations using numerical methods. The quality of these methods is strongly affected by \(dt\). The default value of \(dt\) is 0.1, which will be fine for most of the projections we run in this Vignette. Here it is necessary to reduce the value to 0.01 to avoid introducing any artefacts into the projected values. Decreasing \(dt\) increases the time it takes to run a projection. Here we project the new parameters object for 50 time steps without fishing:
Let’s take a look at how the abundances change through time. We can do this with the plotBiomass() method:
Biomass of the community model when sigma is reduced to a value of 1.0.
The plot above shows that abundances of the community no longer converge to a stable equilibrium and the dynamics appear to be chaotic. The ecological significance of the change in dynamics, and of the ability of simple community models to show chaotic behaviour, is still being debated. It can be argued that the size of the oscillations are too large to be ‘true’. Additionally, when a trait-based model is implemented, the magnitude of the oscillations are much smaller. The next section is about the trait based model.