We continue our investigation of size-spectrum dynamics in simple examples. Where in the previous tutorial notebook we had studied the size-spectrum dynamics of a single species in a fixed power-law background, here we will study a species in a multi-species background, where the background species interact with the foreground species and vice versa.

Again in this tutorial we are not yet interested in a realistic model, but in a model that helps us to get a feel for the multispecies effects that we can expect. We start looking at the real world only in tomorrow’s notebook.

A trait-based multi-species model

First we will create a generic model consisting of several interacting species. We use another wrapper function newTraitParams() that sets up a MizerParams object choosing a lot of defaults for us.

params_background <- newTraitParams(
    lambda = 2.12, # Exponent of community spectrum
    no_sp = 16,    # Number of species
    min_w_inf = 1, # Asymptotic size of smallest species
    max_w_inf = 10^6,  # Asymptotic size of largest species
    min_w = 10^-5,     # Egg size
    min_w_mat = 0.1,   # Maturity size of smallest species
    min_w_pp = 10^-12, # Size of smallest plankton
    bmort_prop = 0.2   # Proportion of mortality that is not due to predation
    )
The search volume exponent q has been set to 0.786666666666667 so as to lead to the desired community exponent lambda = 2.12.
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
plotSpectra(params_background, power = 2, total = TRUE)

We only had to specify a small number of parameters and the function makes default choices for all the others.

Note how the total community, aggregated over all species and the plankton, follows a Sheldon spectrum with the Sheldon exponent lambda of our choice.

These species will form our background, modelling all the species about which we do not have much information and don’t really care about. We tell mizer about this with

params_background <- markBackground(params_background)

This will for example have the effect that these species are plotted in grey. We also don’t want to fish them, so we set their catchability to zero.

params_background@species_params$catchability[] <- 0

Adding a foreground species

Let’s look at hake. Let’s assume we know the following parameters for hake:

species_params <- tribble(
    ~species, ~w_inf, ~w_mat, ~beta, ~sigma, ~a,      ~b,    ~k_vb, ~l50, ~l25,
    "Hake",   10470,  250,    11,    1.1,    0.00667, 3.035, 0.178, 16.6, 16
)
species_params$sel_func <- "sigmoid_length"

Here w_inf is 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\), k_vb the von Bertalanffy growth curve parameter and l50 and l25 the parameters of the sigmoid_length() selectivity function.

In practice you would probably not type them into R as above but create a spreadsheet with the parameters and then read that in with read_csv().

We can now add Hake to our existing model with the addSpecies() function.

params <- addSpecies(params_background,
                     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 52.5 years.
plotSpectra(params, power = 2, total = TRUE)

Increasing abundance of foreground species

Initially this has added hake with a low abundance, so that there is little effect of the presence of the hake on the background species. We can now slowly increase the abundance of the hake, making sure we always run to steady state.

First we double the abundance of the foreground species with respect to the background species.

params <- rescaleAbundance(params, factor = 2) %>% 
    steady()
Simulation run in steady() did not converge after 105 years. Residual relative rate of change = 0.0111441404629417
plotSpectra(params, power = 2, total = TRUE)

Now we double it again, which we can do simply by running the previous R chunk again.

After doing this 4 times we obtain

The multi-species effects are now clearly visible.

Exploring with the tuning gadget

Instead of exploring further by hand, we now switch to a convenient shiny gadget, which provides sliders for changing parameters and a set of plots that update in real time to reflect the effects of the parameter changes. I will demonstrate this on the projector.

params_new <- tuneParams(params)
LS0tCnRpdGxlOiAiQSBmaXJzdCBtdWx0aS1zcGVjaWVzIHNpemUtc3BlY3RydW0gbW9kZWwiCmF1dGhvcjogIkd1c3RhdiBEZWxpdXMiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogVFJVRQotLS0KCldlIGNvbnRpbnVlIG91ciBpbnZlc3RpZ2F0aW9uIG9mIHNpemUtc3BlY3RydW0gZHluYW1pY3MgaW4gc2ltcGxlIGV4YW1wbGVzLgpXaGVyZSBpbiB0aGUgcHJldmlvdXMgdHV0b3JpYWwgbm90ZWJvb2sgd2UgaGFkIHN0dWRpZWQgdGhlIHNpemUtc3BlY3RydW0KZHluYW1pY3Mgb2YgYSBzaW5nbGUgc3BlY2llcyBpbiBhIGZpeGVkIHBvd2VyLWxhdyBiYWNrZ3JvdW5kLCBoZXJlIHdlIHdpbGwgc3R1ZHkKYSBzcGVjaWVzIGluIGEgbXVsdGktc3BlY2llcyBiYWNrZ3JvdW5kLCB3aGVyZSB0aGUgYmFja2dyb3VuZCBzcGVjaWVzIGludGVyYWN0CndpdGggdGhlIGZvcmVncm91bmQgc3BlY2llcyBhbmQgdmljZSB2ZXJzYS4KCkFnYWluIGluIHRoaXMgdHV0b3JpYWwgd2UgYXJlIG5vdCB5ZXQgaW50ZXJlc3RlZCBpbiBhIHJlYWxpc3RpYyBtb2RlbCwgYnV0IGluCmEgbW9kZWwgdGhhdCBoZWxwcyB1cyB0byBnZXQgYSBmZWVsIGZvciB0aGUgbXVsdGlzcGVjaWVzIGVmZmVjdHMgdGhhdCB3ZSBjYW4KZXhwZWN0LiBXZSBzdGFydCBsb29raW5nIGF0IHRoZSByZWFsIHdvcmxkIG9ubHkgaW4gdG9tb3Jyb3cncyBub3RlYm9vay4KCgojIEEgdHJhaXQtYmFzZWQgbXVsdGktc3BlY2llcyBtb2RlbAoKRmlyc3Qgd2Ugd2lsbCBjcmVhdGUgYSBnZW5lcmljIG1vZGVsIGNvbnNpc3Rpbmcgb2Ygc2V2ZXJhbCBpbnRlcmFjdGluZyBzcGVjaWVzLgpXZSB1c2UgYW5vdGhlciB3cmFwcGVyIGZ1bmN0aW9uIGBuZXdUcmFpdFBhcmFtcygpYCB0aGF0IHNldHMgdXAgYSBNaXplclBhcmFtcwpvYmplY3QgY2hvb3NpbmcgYSBsb3Qgb2YgZGVmYXVsdHMgZm9yIHVzLgpgYGB7cn0KcGFyYW1zX2JhY2tncm91bmQgPC0gbmV3VHJhaXRQYXJhbXMoCiAgICBsYW1iZGEgPSAyLjEyLCAjIEV4cG9uZW50IG9mIGNvbW11bml0eSBzcGVjdHJ1bQogICAgbm9fc3AgPSAxNiwgICAgIyBOdW1iZXIgb2Ygc3BlY2llcwogICAgbWluX3dfaW5mID0gMSwgIyBBc3ltcHRvdGljIHNpemUgb2Ygc21hbGxlc3Qgc3BlY2llcwogICAgbWF4X3dfaW5mID0gMTBeNiwgICMgQXN5bXB0b3RpYyBzaXplIG9mIGxhcmdlc3Qgc3BlY2llcwogICAgbWluX3cgPSAxMF4tNSwgICAgICMgRWdnIHNpemUKICAgIG1pbl93X21hdCA9IDAuMSwgICAjIE1hdHVyaXR5IHNpemUgb2Ygc21hbGxlc3Qgc3BlY2llcwogICAgbWluX3dfcHAgPSAxMF4tMTIsICMgU2l6ZSBvZiBzbWFsbGVzdCBwbGFua3RvbgogICAgYm1vcnRfcHJvcCA9IDAuMiAgICMgUHJvcG9ydGlvbiBvZiBtb3J0YWxpdHkgdGhhdCBpcyBub3QgZHVlIHRvIHByZWRhdGlvbgogICAgKQoKcGxvdFNwZWN0cmEocGFyYW1zX2JhY2tncm91bmQsIHBvd2VyID0gMiwgdG90YWwgPSBUUlVFKQpgYGAKCldlIG9ubHkgaGFkIHRvIHNwZWNpZnkgYSBzbWFsbCBudW1iZXIgb2YgcGFyYW1ldGVycyBhbmQgdGhlIGZ1bmN0aW9uIG1ha2VzCmRlZmF1bHQgY2hvaWNlcyBmb3IgYWxsIHRoZSBvdGhlcnMuCgpOb3RlIGhvdyB0aGUgdG90YWwgY29tbXVuaXR5LCBhZ2dyZWdhdGVkIG92ZXIgYWxsIHNwZWNpZXMgYW5kIHRoZSBwbGFua3RvbiwKZm9sbG93cyBhIFNoZWxkb24gc3BlY3RydW0gd2l0aCB0aGUgU2hlbGRvbiBleHBvbmVudCBgbGFtYmRhYCBvZiBvdXIgY2hvaWNlLgoKVGhlc2Ugc3BlY2llcyB3aWxsIGZvcm0gb3VyIGJhY2tncm91bmQsIG1vZGVsbGluZyBhbGwgdGhlIHNwZWNpZXMgYWJvdXQgd2hpY2ggd2UKZG8gbm90IGhhdmUgbXVjaCBpbmZvcm1hdGlvbiBhbmQgZG9uJ3QgcmVhbGx5IGNhcmUgYWJvdXQuIFdlIHRlbGwgbWl6ZXIgYWJvdXQKdGhpcyB3aXRoCmBgYHtyfQpwYXJhbXNfYmFja2dyb3VuZCA8LSBtYXJrQmFja2dyb3VuZChwYXJhbXNfYmFja2dyb3VuZCkKYGBgClRoaXMgd2lsbCBmb3IgZXhhbXBsZSBoYXZlIHRoZSBlZmZlY3QgdGhhdCB0aGVzZSBzcGVjaWVzIGFyZSBwbG90dGVkIGluIGdyZXkuCldlIGFsc28gZG9uJ3Qgd2FudCB0byBmaXNoIHRoZW0sIHNvIHdlIHNldCB0aGVpciBjYXRjaGFiaWxpdHkgdG8gemVyby4KYGBge3J9CnBhcmFtc19iYWNrZ3JvdW5kQHNwZWNpZXNfcGFyYW1zJGNhdGNoYWJpbGl0eVtdIDwtIDAKYGBgCgoKCiMgQWRkaW5nIGEgZm9yZWdyb3VuZCBzcGVjaWVzCgpMZXQncyBsb29rIGF0IGhha2UuIExldCdzIGFzc3VtZSB3ZSBrbm93IHRoZSBmb2xsb3dpbmcgcGFyYW1ldGVycyBmb3IgaGFrZToKYGBge3J9CnNwZWNpZXNfcGFyYW1zIDwtIHRyaWJibGUoCiAgICB+c3BlY2llcywgfndfaW5mLCB+d19tYXQsIH5iZXRhLCB+c2lnbWEsIH5hLCAgICAgIH5iLCAgICB+a192Yiwgfmw1MCwgfmwyNSwKICAgICJIYWtlIiwgICAxMDQ3MCwgIDI1MCwgICAgMTEsICAgIDEuMSwgICAgMC4wMDY2NywgMy4wMzUsIDAuMTc4LCAxNi42LCAxNgopCnNwZWNpZXNfcGFyYW1zJHNlbF9mdW5jIDwtICJzaWdtb2lkX2xlbmd0aCIKYGBgCkhlcmUgYHdfaW5mYCBpcyB0aGUgYXN5bXB0b3RpYyBzaXplLCBgd19tYXRgIHRoZSBtYXR1cml0eSBzaXplLCBgYmV0YWAgYW5kCmBzaWdtYWAgdGhlIHBhcmFtZXRlcnMgb2YgdGhlIGxvZy1ub3JtYWwgcHJlZGF0aW9uIGtlcm5lbCwgYGFgIGFuZCBgYmAgdGhlCnBhcmFtZXRlciBpbiB0aGUgbGVuZ3RoIHRvIHdlaWdodCBjb252ZXJzaW9uICR3ID0gYVwsIGxeYiQsIGBrX3ZiYCB0aGUgdm9uCkJlcnRhbGFuZmZ5IGdyb3d0aCBjdXJ2ZSBwYXJhbWV0ZXIgYW5kIGBsNTBgIGFuZCBgbDI1YCB0aGUgcGFyYW1ldGVycyBvZiB0aGUKYHNpZ21vaWRfbGVuZ3RoKClgIHNlbGVjdGl2aXR5IGZ1bmN0aW9uLgoKSW4gcHJhY3RpY2UgeW91IHdvdWxkIHByb2JhYmx5IG5vdCB0eXBlIHRoZW0gaW50byBSIGFzIGFib3ZlIGJ1dCBjcmVhdGUgYSAKc3ByZWFkc2hlZXQgd2l0aCB0aGUgcGFyYW1ldGVycyBhbmQgdGhlbiByZWFkIHRoYXQgaW4gd2l0aCBgcmVhZF9jc3YoKWAuCgpXZSBjYW4gbm93IGFkZCBIYWtlIHRvIG91ciBleGlzdGluZyBtb2RlbCB3aXRoIHRoZSBgYWRkU3BlY2llcygpYCBmdW5jdGlvbi4KYGBge3J9CnBhcmFtcyA8LSBhZGRTcGVjaWVzKHBhcmFtc19iYWNrZ3JvdW5kLAogICAgICAgICAgICAgICAgICAgICBzcGVjaWVzX3BhcmFtcywKICAgICAgICAgICAgICAgICAgICAgaW5pdGlhbF9lZmZvcnQgPSAxKSAlPiUgCiAgICBzdGVhZHkoKQpwbG90U3BlY3RyYShwYXJhbXMsIHBvd2VyID0gMiwgdG90YWwgPSBUUlVFKQpgYGAKCiMgSW5jcmVhc2luZyBhYnVuZGFuY2Ugb2YgZm9yZWdyb3VuZCBzcGVjaWVzCkluaXRpYWxseSB0aGlzIGhhcyBhZGRlZCBoYWtlIHdpdGggYSBsb3cgYWJ1bmRhbmNlLCBzbyB0aGF0IHRoZXJlIGlzIGxpdHRsZQplZmZlY3Qgb2YgdGhlIHByZXNlbmNlIG9mIHRoZSBoYWtlIG9uIHRoZSBiYWNrZ3JvdW5kIHNwZWNpZXMuIFdlIGNhbiBub3cKc2xvd2x5IGluY3JlYXNlIHRoZSBhYnVuZGFuY2Ugb2YgdGhlIGhha2UsIG1ha2luZyBzdXJlIHdlIGFsd2F5cyBydW4gdG8gc3RlYWR5CnN0YXRlLgoKRmlyc3Qgd2UgZG91YmxlIHRoZSBhYnVuZGFuY2Ugb2YgdGhlIGZvcmVncm91bmQgc3BlY2llcyB3aXRoIHJlc3BlY3QgdG8gdGhlCmJhY2tncm91bmQgc3BlY2llcy4KYGBge3J9CnBhcmFtcyA8LSByZXNjYWxlQWJ1bmRhbmNlKHBhcmFtcywgZmFjdG9yID0gMikgJT4lIAogICAgc3RlYWR5KCkKcGxvdFNwZWN0cmEocGFyYW1zLCBwb3dlciA9IDIsIHRvdGFsID0gVFJVRSkKYGBgCgpOb3cgd2UgZG91YmxlIGl0IGFnYWluLCB3aGljaCB3ZSBjYW4gZG8gc2ltcGx5IGJ5IHJ1bm5pbmcgdGhlIHByZXZpb3VzIFIgCmNodW5rIGFnYWluLgoKQWZ0ZXIgZG9pbmcgdGhpcyA0IHRpbWVzIHdlIG9idGFpbgpgYGB7ciBlY2hvPUZBTFNFfQpwbG90U3BlY3RyYShwYXJhbXMsIHBvd2VyID0gMiwgdG90YWwgPSBUUlVFKQpgYGAKClRoZSBtdWx0aS1zcGVjaWVzIGVmZmVjdHMgYXJlIG5vdyBjbGVhcmx5IHZpc2libGUuCgojIEV4cGxvcmluZyB3aXRoIHRoZSB0dW5pbmcgZ2FkZ2V0Ckluc3RlYWQgb2YgZXhwbG9yaW5nIGZ1cnRoZXIgYnkgaGFuZCwgd2Ugbm93IHN3aXRjaCB0byBhIGNvbnZlbmllbnQgc2hpbnkKZ2FkZ2V0LCB3aGljaCBwcm92aWRlcyBzbGlkZXJzIGZvciBjaGFuZ2luZyBwYXJhbWV0ZXJzIGFuZCBhIHNldCBvZiBwbG90cyB0aGF0CnVwZGF0ZSBpbiByZWFsIHRpbWUgdG8gcmVmbGVjdCB0aGUgZWZmZWN0cyBvZiB0aGUgcGFyYW1ldGVyIGNoYW5nZXMuIEkgd2lsbApkZW1vbnN0cmF0ZSB0aGlzIG9uIHRoZSBwcm9qZWN0b3IuCgpgYGB7ciBldmFsPUZBTFNFfQpwYXJhbXNfbmV3IDwtIHR1bmVQYXJhbXMocGFyYW1zKQpgYGAK