In pwoodworth-jefcoats/Size-Based-Modeling: Mizer With The Effects Of Temperature Included

Load libraries

Mizer 2.0.1 updated at 23/05/2020 or later version is needed to use the background extension.

A. First and best option: download the most updated version of the package directly from GitHub using devtools::install_github("sizespectrum/mizer")

B. Second option: merge upstream master with your branch and upload code from local repository

rm(list = ls())

# option A
# devtools::install_github("sizespectrum/mizer") # run only once
library(mizer)

# option 2 
# library(Jmisc)
# path<-"/Users/nov017/R-projects/Mizer-trait/R"
# sourceAll(path)

Set up multispecies model

All parameters are made up

species_params = data.frame(species = c("speciesA", "speciesB"), w_inf = c(500, 5000), k_vb = c(0.8, 0.3), w_min = c(0.001, 0.001), beta = c(1000,100), sigma = c(3,3))

# available PP to species 
species_params$interaction_resource <- c(1,0.5)

# create the param object
# note: some def values are different from old mizer (e.g. kappa)
params <- newMultispeciesParams(species_params, no_w = 200, kappa = 0.0001) |>
    steady(tol = 0.001)

# run projections:
sim <- project(params, t_max = 500, effort = 0)
plot(sim)

Plankton forcing through time

We want to enable the background resource to change through time via external forcing files used form ISIMIP 3a. We have external code to create these estimates. So here the user will only need to alter how: kappa, lambda, and r_pp change through time. For now, the background resource is only plankton but later this could type of forcing could also be used for the benthic resource or detritus.

Let's start by having a look at how the resource is structured in params. Then we make up a time series of changing plankton parameters through time. For this to work in mizer it will need to have the same timesteps as rest of the model.

# set up new resource forcing "function" - just changes to different timeslot in the n_pp_array array
plankton_state <- new.env(parent = emptyenv())
plankton_state$time <- 1
plankton_forcing <- function(params, t, ...) {
  time_idx <- floor(t) + 1  
  return(other_params(params)$n_pp_array[time_idx,])  
    # plankton_state$time <- plankton_state$time + 1

  }
# set up array of n_pp - currently just initial values
times = 0:500
sizes = names(params@initial_n_pp) 
n_pp_array <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))

# keep constant
# for (i in 1:501) n_pp_array[i,] <- initial_n_pp(params)

# add a slight amount of variation in intercept each timestep
for (i in 1:501) n_pp_array[i,] <- params@initial_n_pp*runif(1, 0.5, 1.5)
other_params(params)$n_pp_array <- n_pp_array 
params <- setResource(params, resource_dynamics = "plankton_forcing")

Run project() using the newly set up parameters

sim2 <- project(params, dt = 0.1, t_max = 500, effort = 0) 
plot(sim2)

Create parameters that must be input by the user

For each species in the model, users must input the minimum and maximum values of its thermal tolerance range. These values are used to scale temperature effects to values ranging from 0 - 1.

species_params(params)$temp_min <- c(15, 10)
species_params(params)$temp_max <- c(25, 20)

Create parameters that are derived from user-input parameters

To scale the effect of temperature on encounter rate and predation mortality to a value ranging from 0 - 1, it is necessary to divide by the maximum possible value for each species. To scale the effect of temperature on metabolism to a value ranging form 0 - 1, it is necessary to subtract the minimum value for each species and then divide by the range. This requires a bit of straightforward arithmetic, and users could do this on their end if they're so inclined. These parameters handle the necessary arithmetic.

# Create parameter for scaling encounter and mortality rates
species_params(params)$encounterpred_scale <- rep(NA, length(species_params(params)$temp_min))

for (indv in seq(1:length(species_params(params)$temp_min))) {

            # Create a vector of all temperatures each species might encounter
      temperature <- seq(species_params(params)$temp_min[indv], species_params(params)$temp_max[indv], by = 0.1)

      # Find the maximum value of the unscaled effect of temperature on encounter and predation rate for each species 
            species_params(params)$encounterpred_scale[indv] <- max((temperature) * (temperature - species_params(params)$temp_min[indv]) * (species_params(params)$temp_max[indv] - temperature))

}

# Determine the minimum, maximum, and range of value for the effect of temperature on metabolism

    min_metab_value <- (exp(25.22 - (0.63/((8.62e-5)*(273 + species_params(params)$temp_min)))))
        max_metab_value <- (exp(25.22 - (0.63/((8.62e-5)*(273 + species_params(params)$temp_max)))))

        species_params(params)$metab_min <- min_metab_value
        species_params(params)$metab_range <- max_metab_value - min_metab_value

Create an ocean_temp parameter

# Create temperature array and fill it
times <- 0:500
species <- species_params(params)$species
ocean_temp_array <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))

temp_inc <- 0
for (i in 1:501) {
  ocean_temp_array[i,] <- c(17 + temp_inc, 17 + temp_inc)
  temp_inc <- temp_inc + 0.01
}

other_params(params)$ocean_temp <- ocean_temp_array

Create the new rate functions

Temperature will be added to the functions to determine encounter rate, predation mortality, and energy available for growth and reproduction.

To scale encounter rate and predation mortality with temperature, we're essentially taking a temperature-dependent proportion of the value calculated by the mizerEncounter and mizerPredMort functions. If species are at their thermal optimum, we take the full value. Elsewhere in their thermal range, we take a proportion that goes to zero at the limits of species' thermal tolerance.

# Calculate the temperature scaling factor for the encounter rate, 
# predation mortality rate and resource mortality rate
scaled_temp_effect <- function(t) {
    # Using t+1 to avoid calling ocean_temp[0,] at the first time step
    temp_at_t <- other_params(params)$ocean_temp[t + 1,]

    # Calculate unscaled temperature effect using a generic polynomial rate equation
    unscaled_temp_effect <- 
        temp_at_t * (temp_at_t - species_params(params)$temp_min) * 
        (species_params(params)$temp_max - temp_at_t)

    # Scale using new parameter
    scaled_temp_effect <- unscaled_temp_effect / species_params(params)$encounterpred_scale

    # Set temperature effect to 0 if temperatures are outside thermal tolerance limits
    above_max <- temp_at_t > species_params(params)$temp_max
    below_min <- temp_at_t < species_params(params)$temp_min

    scaled_temp_effect[above_max | below_min] = 0

    scaled_temp_effect
}

therMizerEncounter <- function(params, t, ...) {

      # Calculate maximum possible encounter rate
      max_encounter <- mizerEncounter(params, t, ...)

      # Apply temperature effect
      # return(sweep(max_encounter, 1, scaled_temp_effect, '*', check.margin = FALSE))
      return(max_encounter * scaled_temp_effect(t))

}

therMizerPredRate <- function(params, t, ...) {

      # Calculate maximum possible encounter rate
      max_predrate <- mizerPredRate(params, t, ...)

      # Apply temperature effect
      # return(sweep(max_predrate, 1, scaled_temp_effect, '*', check.margin = FALSE))
      return(max_predrate * scaled_temp_effect(t))

}

To calculate the effect of temperature on metabolism, we use an Arrhenius function to scale the cost of metabolism. When species are at their thermal maximum, the cost of metabolism is at its maximum. When species are at their thermal minimum, the cost of metabolism is at its minimum

therMizerEReproAndGrowth <- function(params, t, encounter, feeding_level, ...) {

    # Using t+1 to avoid calling ocean_temp[0,] at the first time step
    temp_at_t <- other_params(params)$ocean_temp[t + 1,]

    # Arrhenius equation
    unscaled_temp_effect <- (exp(25.22 - (0.63/((8.62e-5)*(273 + temp_at_t)))))

    # Arrhenius equation scaled to a value between 0 and 1
        temp_effect_metabolism <- (unscaled_temp_effect - species_params(params)$metab_min) / species_params(params)$metab_range

        # Set temperature effect to 0 if temperatures are outside thermal tolerance limits
    above_max <- temp_at_t > species_params(params)$temp_max
    below_min <- temp_at_t < species_params(params)$temp_min

    temp_effect_metabolism[above_max | below_min] = 0

        # Apply scaled Arrhenius value to metabolism
    sweep((1 - feeding_level) * encounter, 1,
               species_params(params)$alpha, "*", check.margin = FALSE) - 
      metab(params)*temp_effect_metabolism  

}

Update the functions

params <- setRateFunction(params, "Encounter", "therMizerEncounter")
params <- setRateFunction(params, "PredRate", "therMizerPredRate")
params <- setRateFunction(params, "EReproAndGrowth", "therMizerEReproAndGrowth")

Run project() using the newly set up parameters

sim3 <- project(params, t_max = 500, effort = 0) 
plot(sim3)

pwoodworth-jefcoats/Size-Based-Modeling documentation built on Sept. 15, 2023, 8:13 a.m.