A systems-modelling approach to predict biological responses to extreme heat

This file contains the code and data needed for demonstrating how a systems modelling approach can be applied in practice and can be found on GitHub. It is the code for reproducing the examples in Box 1 and Box 2 for the manuscript:

Noble et al. 2025. A systems-modelling approach to predict biological responses to extreme heat. Trends in Ecology and Evolution, under review.

It also demonstrates how to couple models from different disciplines to predict the impacts of extreme heat on plant-insect interactions.

Please get in touch with Daniel Noble (daniel.noble@anu.edu.au) and Mike Kearney (m.kearney@unimelb.edu.au) if you have any questions about the code or the models.

1 Introduction

Here, we provide the code for reproducing the systems-modelling approaches described in Boxes 1 & 2.

2 Load packages

Our code is built using the statistical software R. We use a number of packages to run our simulations, including NicheMapR [1,2], microclima [3,4], and plantecophys [5].

We recommend that readers interested in working through this the code download our GitHub repository, which contains all the data required to run the models. This code file can also located in the docs/supp.qmd file. You should then be able to run the code in this document to reproduce the figures and tables without having to setup your repository.

The code below loads the required packages that are needed.

# You will need to install a few packages first, a couple from GitHub as they are not on CRAN. Note that you will need to remove '#' from the code below to install the packages.
# install.packages("pacman")
# devtools::install_github('mrke/NicheMapR')
# devtools::install_github("ilyamaclean/microclima")

# Load required libraries after the installation
pacman::p_load(devtools, deSolve, ggplot2, dplyr, tidyr, purrr, patchwork, equatags, flextable, NicheMapR, terra, microclima, plantecophys, magick, here, ggrepel, lutz, R.matlab, magick, grid, DT, latex2exp, popbio, popdemo)

Throughout we’ll also need some helper functions. These can be added to your working environment by running the code below.

#' @title convert_doy
#' @description Convert day of year (DOY) to date, month-day, ymd, or julian date
#' @param doy Numeric vector of day of year
#' @param year Numeric vector of year
#' @param format_type Character vector of format type. Default is "date"
#' @return Date, month, month-day, ymd, or julian date
#' @examples
#' convert_doy(120, 2024, "date") # "2024-04-29"
#' convert_doy(120, 2024, "month_day") # "April 29" 
#' convert_doy(120, 2024, "ymd") # "2024-04-29"

convert_doy <- function(doy, year, format_type = "date") {
  # Convert DOY to Date
  date <- as.Date(doy - 1, origin = paste0(year, "-01-01"))
  
  # Return different formats based on format_type
  switch(format_type,
         "date" = date,                     # Default: Date format
         "month_day" = format(date, "%B %d"), # "April 29"
         "month" = format(date, "%B"), # "April"
         "ymd" = format(date, "%Y-%m-%d"),    # "2024-04-29"
         "julian" = format(date, "%Y%j"),     # "2024120"
         stop("Invalid format_type. Choose from 'date', 'month_day', 'ymd', or 'julian'.")
  )
}

#' @title convert_time
#' @description Convert decimal time to hh:mm, hh:mm:ss, or 12-hour format
#' @param time_minutes Numeric vector of decimal time in minutes
#' @param format_type Character vector of format type. Default is "hh:mm"
#' @return Time in hh:mm, hh:mm:ss, or 12-hour format
#' @examples
#' convert_time(13.5, "hh:mm") # "13:30"
#' convert_time(13.5, "hh:mm:ss") # "13:30:00"
#' convert_time(13.5, "12-hour") # "1:30 PM"

convert_time_for_plot <- function(time_minutes, format_type = "decimal") {
  # Extract hours and minutes
  hours <- time_minutes %/% 60  
  minutes <- time_minutes %% 60  

  # Format output based on format_type
  switch(format_type,
         "decimal" = hours + minutes / 60,  # e.g., 13.5 for 13:30
         "POSIXct" = as.POSIXct(sprintf("%02d:%02d", hours, minutes), format="%H:%M", tz="UTC"),  # e.g., "13:30:00 UTC"
         "hms" = hms::hms(hours * 3600 + minutes * 60),  # e.g., 13:30:00
         stop("Invalid format_type. Choose from 'decimal', 'POSIXct', or 'hms'.")
  )
}

Note

gFortran can be a problem for NicheMapR so you need to set the paths properly if you are getting errors. Follow details here. The following code will help you set the paths properly.

# But do the following:
# 1. Install gfortran using homebrew in the terminal/bash window
   brew install gcc 
# 2. Take note of the path to gfortran. You'll need this for later
    which gfortran

# 3. In home directory of terminal create a Makevars file
   mkdir ~/.R
   touch ~/.R/Makevars
   open  ~/.R/Makevars

# 4. Add the following to the Makevars file. Note that you need to make sure your version of fortran is correct and replace 124.2.0_1 with your version. You can find version in terminal by typing `gfortran --version`. In addition you need to replace the path to gfortran with the path you found in step 2.
    FC = /opt/homebrew/Cellar/gcc/14.2.0_1/bin/gfortran
   F77 = /opt/homebrew/Cellar/gcc/14.2.0_1/bin/gfortran
 FLIBS = -L/opt/homebrew/Cellar/gcc/14.2.0_1/lib/gcc/11

#5. Restart R and try to install the package again via devtools. It should now work.

3 Box 1 Code: Building a systems modelling approach to capture multi-species exposure to extreme heat

3.1 Getting the microclimate data

We first need to use climate data to calculate the microclimate at our location of interest. Here, we use the micro_silo function to calculate the microclimate at a location near Renmark, South Australia (longitude = 140.799, latitude = -34.159) for the years 2018-2019.

# First, set the location where NCEAP climate data is located. This is a global dataset. We need this to calculate the microclimate. Only need to do once.
    get.global.climate(folder=here("NCEP_dan", "NCEAP" ))

# Location we want
  longlat <- c(140.79970255974192, -34.159176277782564) # Near Renmark, South Australia
    micro <- micro_silo(     loc = longlat,      # Location
                          dstart = "01/01/2018", # Start year
                         dfinish = "31/12/2019", # End year
                        minshade = 40,           # Minimum shade, %
                        maxshade = 80,           # Maximum shade, %
                          Usrhyt = 0.2,          # Height above ground, m 
                      microclima = 1,            # Use mircoclima 
                              PC = -1500,        # Water potential for stomatal closure
                              SP = 10,           # Stability parameter for stomatal closure
                              
                           email = "daniel.noble@anu.edu.au") # Replace with your email

# Save the microclimate object for later in case we need it again as it avoids needing to rerun.                           
  save(micro, file = here("output", "microclim", "microclimate.RDa"))

3.2 Plant biophysical model with stomatal conductance

Now that we have the microclimate data, we can run the biophysical model to calculate the realised environmental conditions for a plant leaf using the ectotherm function in NicheMapR. We first need to set up the parameters for the plant leaf.

   # Load the microclimate object
    load(here("output", "microclim", "microclimate.RDa"))

  # First, lets extract the relevant information from the microclimate object
      PDIFs <- micro$diffuse_frac             # use variable diffuse fraction
        P_a <- get_pressure(micro$elev)       # air pressure, Pa
        P_a <- rep(P_a, nrow(micro$metout))   # create hourly vector
      dates <- micro$dates                    # dates

  # Set leaf functional traits
         live <- 0      # don't simulate behaviour or respiration
         leaf <- 1      # use the stomatal conductance formulation for evaporation
         Ww_g <- 1      # wet weight, g
        shape <- 2      # 0=plate, 1=cylinder, 2=ellipsoid
      shape_b <- 0.0025 # ratio of b axis:a axis for ellipsoid 
      shape_c <- 0.1176 # ratio of c axis:a axis for ellipsoid 
      g_vs_ab <- 0.4    # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
      g_vs_ad <- 0.0    # leaf vapour conductance, adaxial (top of leaf), mol/m2/s
  g_vs_ab_max <- 0.4    # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
  g_vs_ad_max <- 0.0    # leaf vapour conductance, adaxial (top of leaf), mol/m2/s
  epsilon_sub <- 0.95   # emissivity of the substrate (-)
  epsilon_sky <- 1      # emissivity of the sky (-), accounted for already in 'TSKY' output from microclimate model so make it 1
    epsilon_L <- 0.97   # emissivity of the leaf (-)
      alpha_L <- 0.5    # solar absorptivity of the leaf (-)
       fatosk <- 0.5    # radiation configuration factor to sky (-)
       fatosb <- 0.5    # radiation configuration factor to substrate (-)
 conv_enhance <- 1.4    # convective enhancement factor (-)
     pct_cond <- 0      # percent of leaf conducting to the ground (%)
       postur <- 3      # leaf oriented horizontally, not tracking sun
          M_1 <- 0      # make metabolic rate zero

# Set up vapour conductance vectors and simulate stomatal closure at night
                              f_open <- rep(1, nrow(micro$metout))
  f_open[micro$metout[,"ZEN"] == 90] <- 0      # close stomata when the sun is set
                            g_vs_abs <- g_vs_ab * f_open 
                            g_vs_ads <- g_vs_ad * f_open

# Run the heat-budget model for the plant leaf. Note that the microclimate object is used to get the relevant environmental temperatures to make the calculations of leaf temperature (by default, ectotherm expects an object in the environment called 'micro' and uses it to obtain the relevant environmental data, but this can be overridden).
    plant <- ectotherm(  leaf = leaf, 
                         live = live,
                         Ww_g = Ww_g, 
                        shape = shape, 
                      shape_b = shape_b, 
                      shape_c = shape_c, 
                      g_vs_ab = g_vs_abs, 
                      g_vs_ad = g_vs_ads, 
                  g_vs_ab_max = g_vs_ab_max, 
                  g_vs_ad_max = g_vs_ad_max, 
                  epsilon_sub = epsilon_sub, 
                  epsilon_sky = epsilon_sky,
                      epsilon = epsilon_L, 
                    alpha_max = alpha_L, 
                    alpha_min = alpha_L,
                          M_1 = M_1,
                       fatosk = fatosk, 
                       fatosb = fatosb, 
                 conv_enhance = conv_enhance, 
                     pct_cond = pct_cond, 
                       postur = postur,
                       preshr = P_a, 
                         PDIF = PDIFs)

# Save model
  plant$dates <- dates  # Add in dates for clarity
  save(plant, file = here("output", "microclim", "plant_micro.RDa"))

3.3 Plant biophysical model with stomatal conductance from plantecophys

Allow for stomatal conductance to be simulated using the plantecophys package. We can do this by calculating stomatal conductance using the Photosyn function in plantecophys. This function requires vapour pressure deficit (VPD) which we can calculate from the microclimate data using the WETAIR function in NicheMapR.

# Need to extract a few things from the plant object we just created
  # Load the plant object
    load(here("output", "microclim", "plant_micro.RDa"))

# Include simulated stomatal conductance from plantecophys
  UMOLPERJ <- 4.57                                                   # mu_mol photons / J
WETAIR_out <- WETAIR(db = plant$metout[,"TALOC"], rh = plant$metout[,"RHLOC"], bp = P_a)
    VPD_Pa <- WETAIR_out$esat - WETAIR_out$e
 photo_out <- Photosyn(Tleaf = plant$environ[,"TC"],                 # Leaf temperature (C)
                          g1 = 5.2,                                  # Parameter of Ball-Berry model
                          Ca = 400,                                  # Atmospheric CO2 concentration (ppm)
                         VPD = VPD_Pa/1000,                          # Vapour pressure deficit (kPa)
                      vpdmin = 0.5,                                  # lower threshold on VPD
                        PPFD = plant$metout[,"SOLR"] * UMOLPERJ * 1, # mu mol m-2 s-1
                        Patm = P_a/1000                              # Atmospheric pressure (kPa)
                      )
g_vs_abs_pho <- photo_out$GS

Now calculate leaf temperature using these stomatal conductance values from Photosyn. Note that we have re-written the g_vs_abs_pho to capture the stomatal conductance from Photosyn. We just replace this in the plant biophysical model.

# Rerun the models. Note that the parameters as the same as above except for g_vs_ab which is now from plantecophys 
plant2 <- ectotherm( leaf = leaf, 
                     live = live,
                     Ww_g = Ww_g, 
                    shape = shape, 
                  shape_b = shape_b, 
                  shape_c = shape_c, 
                  g_vs_ab = g_vs_abs_pho, 
                  g_vs_ad = g_vs_ads, 
              g_vs_ab_max = g_vs_ab_max, 
              g_vs_ad_max = g_vs_ad_max, 
              epsilon_sub = epsilon_sub, 
              epsilon_sky = epsilon_sky,
                  epsilon = epsilon_L, 
                alpha_max = alpha_L, 
                alpha_min = alpha_L,
                      M_1 = M_1,
                   fatosk = fatosk, 
                   fatosb = fatosb, 
             conv_enhance = conv_enhance, 
                 pct_cond = pct_cond, 
                   postur = postur,
                   preshr = P_a, 
                    PDIF = PDIFs)
# Save model
  plant2$dates <- dates  # Add in dates for clarity
  save(plant2, file = here("output", "microclim", "plant2_micro.RDa"))

Now, lets have a look at the leaf temperature when stomatal closure behaviour is incorporated.

  # Load the plant object
    load(here("output", "microclim", "plant_micro.RDa"))
    load(here("output", "microclim", "plant2_micro.RDa"))
  
  # Extract the microclimate temperatures. Because this was simulated for a plant the 'body tempertaure'is the leaf temperature and is stored as 'TC' in the environ object. Note now that this new biophysical model has stomatal behaviour incorporated in the new plant2 object.
           environ_plant <- as.data.frame(plant$environ)
      environ_plant$dates <- plant$dates  # For clarity add in dates from microclimate that ecto model uses

            environ_plant_plantecophys2 <- as.data.frame(plant2$environ)
      environ_plant_plantecophys2$dates <- plant2$dates  # For clarirty add in dates from microclimate that ecto model uses

  # Re-organise the data for plotting

  environ_plant2 <- environ_plant %>% filter(dates > "2019-01-15" & dates < "2019-01-18") %>%
              select(dates, TA, TC) %>%
              pivot_longer(cols = c(TA, TC), 
                           names_to = "Temps", 
                           values_to = "Temperature") %>% 
              mutate(Temps = ifelse(Temps == "TA", "Air Temperature (2 m) (40% shade)", "Leaf Temperature (Stomata Open) (40% shade)"))

  environ_plant_plantecophys <- environ_plant_plantecophys2 %>% filter(dates > "2019-01-15" & dates < "2019-01-18") %>%
              dplyr::select(dates, TC) %>%
              pivot_longer(cols = c(TC), 
                           names_to = "Temps", 
                           values_to = "Temperature") %>% 
              mutate(Temps = ifelse(Temps == "TC",  "Leaf Temperature (Stomatal Closure) (40% shade)"))

  # Combine the two dataframes environ_plant_plantecophys and environ_plant2
  new_environ_plant2 <- rbind(environ_plant_plantecophys, environ_plant2)

  # Produce element of the fgure in box 1. Note it was adjusted in Biorender.com for final version.
  ggplot(new_environ_plant2) + 
              geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend = TRUE, size = 0.8) +
              labs(y = "Temperature (°C)", x = "Date") + theme_classic() +
              theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),legend.position = "none", legend.text = element_text(size = 14), axis.text = element_text(size = 14), axis.title = element_blank()) + 
              scale_colour_viridis_d()
  # ggsave(here("figs", "fig_plant_temp2_box1.png"), width = 8, height = 5, dpi = 600)

Here, we can see that stomatal closure in the middle of the day causes the leaf temperatures to rise substantially compared to if they were always open and transpiring (i.e., evaporative cooling). This is important because it means that the grasshopper will be exposed to higher temperatures when plant stomata close midday.

3.4 Run biophysical-physiological model for two grasshoppers using the plant microclimate as the environment for grasshoppers

Define the biophysical model parameters for the grasshopper. Note we also include DEB model parameters here, which relate to Box 2, but it’s easier to run the model once with all the parameters defined and we can (and will) extract what we need for Box 1.

# We first need to load the DEB model paraemeters. These are estimated using MatLab. We've already estimated these so we'll simply just load them into R.

     pars <- readMat(here("sims", "results_Warramaba_virgo.mat"))

par.names <- unlist(labels(pars$par))
for(i in 1:length(par.names)){
  assign(par.names[i], unlist(pars$par[i]))
}

     pars <- readMat(here("sims", "stat_Warramaba_virgo.mat"))
par.names <- unlist(labels(pars$stat))
for(i in 1:length(par.names)){
  assign(par.names[i], unlist(pars$stat[i]))
}

# Additional parameters we need
        E_sm <- 1800            # J/cm^2, maximum stomach energy content
       fract <- 1               # scaling factor to change size

# Inital state for DEB model

# Egg Stage. Starting conditions
     E_init <- (E.0 * fract ^ 4) / (3E-9 * fract ^ 3)
   E_H_init <- 0
     V_init <- 3E-9
      stage <- 0

# Life history traits
 clutchsize <- 15               # Eggs per clutch
     stages <- 7                # Number of stages, 0 = egg, 1-5 = nymphs, 6 = pre egg-laying adult, 7 = post-laying adult
   S_instar <- rep(s.1, stages) # Stress at instar n: L_n^2/ L_n-1^2 (-)}\cr}
        L_b <- 0.09             # Structural length at birth (cm)}\cr}
 photostart <- 0                # no photoperiod effect on reproduction
photofinish <- 0                # no photoperiod effect on reproduction
      reset <- 0                # don't reset lifecycle if death occurs
   soilnode <- 3                # depth (soil node) at which eggs laid - 3 = 5 cm
         sG <- 0.001           # Gompertz stress coefficient

# Thermoregulatory parameters
       Ww_g <- 0.2              # g, body mass, not used though because DEB model computes this
      shape <- 1                # cylinder
    shape_b <- 10               # ratio of long to short axis of cylinder
     T_pref <- 40               # preferred body temperature (animal will attempt to regulate as close to this value as possible)
    T_F_min <- 13               # degrees C, lower body temperature for foraging
    T_F_max <- 42               # degrees C, upper body temperature for foraging
    T_B_min <- 13               # degrees C, minimum basking temperature
   T_RB_min <- 13               # degrees C, temperature at which animal will move to a basking site
     CT_max <- 50               # degrees C, critical thermal maximum (animal will die if CT_kill = 1 and this threshold is exceeded)
     CT_min <- 0                # degrees C, critical thermal minimum (used by program to determine depth selected when inactive and burrowing)

# Behavioural traits (made stage-specific later)
      diurn <- 1                # Diurnal activity allowed (1) or not (0)?
    nocturn <- 1                # Nocturnal activity allowed (1) or not (0)?
     crepus <- 1                # Crepuscular activity allowed (1) or not (0)?
     burrow <- 0                # Shelter in burrow allowed (1) or not (0)?
 shade_seek <- 0                # Shade seeking allowed (1) or not (0)?
   pct_cond <- 0                # Percentage of animal surface contacting the substrate (\%)

Run the biophysical-physiological model species A.

# First, we need to load in the two microclimates scenarios. We'll do this in a loop so that we can run the DEB model for each of the two microclimates. So, here we'll just create a few vectors that cell specific objects we need.
        plant_micro <- c("plant_micro.RDa")
  plant_mirco_names <- c("stomata_closed")

# We also need to load in the original microclimate object so that we can get the dates and other information from it. This is the original microclimate object that we used to simulate the plant microclimate. We'll also overide the `mirco` object so that the ectotherm model can use it.
  load(here("output", "microclim", "microclimate.RDa"))
  load(here("output", "microclim", "plant_micro.RDa"))

# Use the plant microclimate data to reset the original microclimate object. Note that we need to adjust both shade and nonshade conditions
            

  # Unshaded microclimate
  micro$metout <- plant$metout
  micro$metout [,"TALOC"] <- plant$environ[,"TC"]

  # Update the IR profiles. 
    micro$metout[,"TSKYC"] <- plant$environ[,"TSKY"]  
  
  # Update the relative humidity.
  micro$metout[,"RHLOC"] <- plant$environ[,"RELHUM"]

  # Shaded microclimate  
  microshadmet <- plant$shadmet
  micro$shadmet [,"TALOC"] <- plant$environ [,"TC"]
 
  # Soil microclimate
       micro$soil <- plant$soil
   micro$shadsoil <- plant$shadsoil
  micro$soilmoist <- plant$soilmoist
  micro$shadmoist <- plant$shadmoist

# Get dates from microclimate object
 dates <- micro$dates  # dates, hourly
dates2 <- micro$dates2 # dates, daily

# Set starting date by truncating microclimate model output to when grasshoppers are epected to oviposit eggs. This is the date when the life cylce of the insect is expected to start.
oviposition_date <- '2018-03-31 23:00:00'
            subs <- which( dates > as.POSIXct(oviposition_date))
           subs2 <- which(dates2 > as.POSIXct(oviposition_date))

# Now run the biophysical model for the grasshopper

  ectoA <- ectotherm( Ww_g = Ww_g,
                    shape = shape,
                  shape_b = shape_b,
                   T_pref = T_pref,
                  T_F_min = T_F_min,
                  T_F_max = T_F_max,
                  T_B_min = T_B_min,
                 T_RB_min = T_RB_min,
                   CT_min = CT_min,
                   CT_max = CT_max,
                    diurn = diurn,
                  nocturn = nocturn,
                   crepus = crepus,
               shade_seek = shade_seek,
                   burrow = burrow,
                 pct_cond = pct_cond,
                   nyears = ceiling(length(subs2)/365),
                   metout = micro$metout[subs,],
                     soil = micro$soil[subs,], 
                  soilpot = micro$soilpot[subs,],
                soilmoist = micro$soilmoist[subs,], 
                    humid = micro$humid[subs,],
                    tcond = micro$tcond[subs,], 
                  shadmet = micro$shadmet[subs,],
                 shadsoil = micro$shadsoil[subs,], 
                  shadpot = micro$shadpot[subs,],
                shadmoist = micro$shadmoist[subs,],
                shadhumid = micro$shadhumid[subs,],
                shadtcond = micro$shadtcond[subs,], 
                 rainfall = micro$RAINFALL[subs2],
                minshades = rep(unique(micro$minshade), length(subs)),
                maxshades = rep(unique(micro$maxshade), length(subs)),
                      DEB = 1, # turns DEB model on
               metab_mode = 1, # turns on hemimetabolous model, abp
                        z = z * fract,
                    del_M = del.M,
                     p_Xm = p.Am * 100 / 24 / kap.X,
                    kap_X = kap.X,
                        v = v / 24,
                      kap = kap,
                      p_M = p.M / 24,
                      E_G = E.G,
                    kap_R = kap.R,
                      k_J = k.J  / 24,
                     E_Hb = E.Hb * fract ^ 3,
                     E_Hp = (E.Hj+1e-5) * fract ^ 3,
                     E_Hj = E.Hj * fract ^ 3,
                      h_a = h.a / (24 ^ 2),
                      s_G = sG,   # Gompertz stress coefficient
                      E_0 = E.0 * fract ^ 4,
                     mu_E = mu.E,
                      d_V = d.V,
                      d_E = d.E,
                    T_REF = T.ref,
                      T_A = T.A,
                     T_AL = T.AL,
                     T_AH = T.AH,
                      T_L = T.L,
                      T_H = T.H,
                     E_sm = E_sm,
                   E_init = E_init,
                 E_H_init = E_H_init,
                   V_init = V_init,
                    stage = stage,
               clutchsize = clutchsize,
                   stages = stages,
                 S_instar = rep(s.1, stages),
                      L_b = L_b, # structural length at birth (cm),
               photostart = photostart,
              photofinish = photofinish,
                    reset = reset,
                 soilnode = soilnode)
# Save the ectotherm model
    ectoA$dates <- dates[subs]  # Add in dates for clarity
  save(ectoA, file = here("output", "microclim", paste0("grasshopperspA_", plant_mirco_names, "_Box1.RDa")))

Run the biophysical-physiological model again, but for species B which uses different parameters to capture the different life histories of the species. Note that species B occupies a slightly different microhabitat which is cooler than species A so we will assume they will be closet to leaf temperatures when stomata are open.

# Using the same starting parameters as species A, so will not re-load the matrix again. 
# Additional parameters we need
        E_sm <- 2000            # J/cm^2, maximum stomach energy content
       fract <- 1               # scaling factor to change size

# Inital state for DEB model. Note that DEB can be turned off and only heat and water budget models will be used so some parameters are not needed if DEB is off.

# Egg Stage. Starting conditions
     E_init <- (E.0 * fract ^ 4) / (3E-9 * fract ^ 3)
   E_H_init <- 0
     V_init <- 3E-9
      stage <- 0

# Life history traits
 clutchsize <- 8                # Eggs per clutch
     stages <- 7                # Number of stages, 0 = egg, 1-5 = nymphs, 6 = pre egg-laying adult, 7 = post-laying adult
   S_instar <- rep(s.1, stages) # Stress at instar n: L_n^2/ L_n-1^2 (-)}\cr}
        L_b <- L.b              # Structural length at birth (cm)}\cr}
 photostart <- 0                # no photoperiod effect on reproduction
photofinish <- 0                # no photoperiod effect on reproduction
      reset <- 0                # don't reset lifecycle if death occurs
   soilnode <- 3                # depth (soil node) at which eggs laid - 3 = 5 cm

# Thermoregulatory parameters
       Ww_g <- 0.2              # g, body mass
      shape <- 1                # cylinder
    shape_b <- 10               # ratio of long to short axis of cylinder
     T_pref <- 38               # preferred body temperature species B
    T_F_min <- 13               # degrees C, lower body temperature for foraging
    T_F_max <- 40               # degrees C, upper body temperature for foraging
    T_B_min <- 13               # degrees C, minimum basking temperature
   T_RB_min <- 13               # degrees C, temperature at which animal will move to a basking site
     CT_max <- 47               # degrees C, critical thermal maximum species B
     CT_min <- 0                # degrees C, critical thermal minimum (used by program to determine depth selected when inactive and burrowing)

# Behavioural traits 
      diurn <- 1                # Diurnal activity allowed (1) or not (0)?
    nocturn <- 1                # Nocturnal activity allowed (1) or not (0)?
     crepus <- 1                # Crepuscular activity allowed (1) or not (0)?
     burrow <- 0                # Shelter in burrow allowed (1) or not (0)?
 shade_seek <- 0                # Shade seeking allowed (1) or not (0)?
   pct_cond <- 0                # Percentage of animal surface contacting the substrate (\%)

   plant_micro <- c("plant2_micro.RDa")
  plant_mirco_names <- c("stomata_open")

# Use the plant microclimate data to reset the original microclimate object. Note that we need to adjust both shade and nonshade conditions
  
  # Unshaded microclimate
             micro$metout <- plant2$metout
  micro$metout [,"TALOC"] <- plant2$environ[,"TC"]
   micro$metout[,"TSKYC"] <- plant2$environ[,"TSKY"]  
   micro$metout[,"RHLOC"] <- plant2$environ[,"RELHUM"]

  # Shaded microclimate  
              microshadmet <- plant2$shadmet
  micro$shadmet [,"TALOC"] <- plant2$environ [,"TC"]
 
  # Soil microclimate
       micro$soil <- plant2$soil
   micro$shadsoil <- plant2$shadsoil
  micro$soilmoist <- plant2$soilmoist
  micro$shadmoist <- plant2$shadmoist

# Get dates from microclimate object
 dates <- micro$dates  # dates, hourly
dates2 <- micro$dates2 # dates, daily

# Set starting date by truncating microclimate model output to when grasshoppers are epected to oviposit eggs. This is the date when the life cylce of the insect is expected to start.
oviposition_date <- '2018-03-31 23:00:00'
            subs <- which( dates > as.POSIXct(oviposition_date))
           subs2 <- which(dates2 > as.POSIXct(oviposition_date))

# Now run the DEB life-cyle model for the grasshopper

  ectoB <- ectotherm( Ww_g = Ww_g,
                    shape = shape,
                  shape_b = shape_b,
                   T_pref = T_pref,
                  T_F_min = T_F_min,
                  T_F_max = T_F_max,
                  T_B_min = T_B_min,
                 T_RB_min = T_RB_min,
                   CT_min = CT_min,
                   CT_max = CT_max,
                    diurn = diurn,
                  nocturn = nocturn,
                   crepus = crepus,
               shade_seek = shade_seek,
                   burrow = burrow,
                 pct_cond = pct_cond,
                   nyears = ceiling(length(subs2)/365),
                   metout = micro$metout[subs,],
                     soil = micro$soil[subs,], 
                  soilpot = micro$soilpot[subs,],
                soilmoist = micro$soilmoist[subs,], 
                    humid = micro$humid[subs,],
                    tcond = micro$tcond[subs,], 
                  shadmet = micro$shadmet[subs,],
                 shadsoil = micro$shadsoil[subs,], 
                  shadpot = micro$shadpot[subs,],
                shadmoist = micro$shadmoist[subs,],
                shadhumid = micro$shadhumid[subs,],
                shadtcond = micro$shadtcond[subs,], 
                 rainfall = micro$RAINFALL[subs2],
                minshades = rep(unique(micro$minshade), length(subs)),
                maxshades = rep(unique(micro$maxshade), length(subs)),
                      DEB = 1, # turns DEB model on
               metab_mode = 1, # turns on hemimetabolous model, abp
                        z = z * fract,
                    del_M = del.M,
                     p_Xm = p.Am * 100 / 24 / kap.X,
                    kap_X = kap.X,
                        v = v / 24,
                      kap = kap,
                      p_M = p.M / 24,
                      E_G = E.G,
                    kap_R = kap.R,
                      k_J = k.J  / 24,
                     E_Hb = E.Hb * fract ^ 3,
                     E_Hp = (E.Hj+1e-5) * fract ^ 3,
                     E_Hj = E.Hj * fract ^ 3,
                      h_a = h.a / (24 ^ 2),
                      s_G = s.G,
                      E_0 = E.0 * fract ^ 4,
                     mu_E = mu.E,
                      d_V = d.V,
                      d_E = d.E,
                    T_REF = T.ref,
                      T_A = T.A,
                     T_AL = T.AL,
                     T_AH = T.AH,
                      T_L = T.L,
                      T_H = T.H,
                     E_sm = E_sm,
                   E_init = E_init,
                 E_H_init = E_H_init,
                   V_init = V_init,
                    stage = stage,
               clutchsize = clutchsize,
                   stages = stages,
                 S_instar = rep(s.1, stages),
                      L_b = L.b,
               photostart = photostart,
              photofinish = photofinish,
                    reset = reset,
                 soilnode = soilnode)
# Save the DEB model
    ectoB$dates <- dates[subs]  # Add in dates for clarity
  save(ectoB, file = here("output", "microclim", paste0("grasshopperB_", plant_mirco_names, "_Box1.RDa")))

Now we can have a look at the body temperature of both species under both microclimate scenarios.

## load the two grasshopper model outputs
load(here("output", "microclim", "grasshopperspA_stomata_closed_Box1.RDa"))
load(here("output", "microclim", "grasshopperB_stomata_open_Box1.RDa"))

# Create dataframes for plotting for both species
           temp_spA_open <- data.frame(ectoA$environ)
  temp_spA_open$datetime <- ectoA$dates

         temp_spB_open <- data.frame(ectoB$environ)
temp_spB_open$datetime <- ectoB$dates

# Merge the two dataframes
  temp_combined <- rbind(
    temp_spA_open %>% mutate(Species = "Species A (Stomata Closed)"),
    temp_spB_open %>% mutate(Species = "Species B (Stomata Open)")
  )

# Now subset to the desired dates
  temp_combined_sub <- temp_combined %>% filter(datetime > as.POSIXct('2019-01-15') & datetime < as.POSIXct('2019-01-18'))

# Plot for the figure in Box 1. Note it was adjusted in Biorender.com for final version.
ggplot(temp_combined_sub, aes(x = datetime, y = TA, color = Species)) +
  geom_line(size = 0.8) +
  labs(x = "Date",
       y = "Body Temperature (°C)") +
  theme_classic() +
  theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),legend.position = "none", legend.text = element_text(size = 14), axis.text = element_text(size = 14), axis.title = element_blank()) +
  scale_colour_manual(values = c("Species A (Stomata Closed)" = "#356920", "Species B (Stomata Open)" = "#664229"))

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

# ggsave(here("figs", "fig_grasshopper_temp_box1.png"), width = 8, height = 5, dpi = 600)

Figure 1

4 Box 2 Code: Scaling up extreme heat effects from impacts on individuals of species to communities

Now that we have predictions about body temperatures we can incorporate thermal load accumulation and use the DEB models to predict life-cycles of the two grasshopper species. The output from DEB models can then be used to calculate vital rates such as growth, development, reproduction, and survival, which we can use to calculate population growth rates (\(r_{max}\)). We can then model population dynamics and interactions between the two species.

Note that in Box 1, we ran the biophysical model with DEB turned on, but we only extracted body temperature information. So, in essence, we already coupled heat and water-budget models with dynamic energy budget models for grasshoppers. Here, we will extract the DEB model and do a bit more with it.

4.1 Capturing impacts of extreme heat on physiological function across the life-cycle

First, we’ll need some functions to extract life-cycle details from the DEB model output and compute survival probabilities.

#' @title alive_data: Subset DEB output to period when alive
#' @description Extracts DEB model output corresponding only to periods when the individual
#' is alive (post-hatch). This is done by subsetting the `debout` component
#' of the ectotherm model output (`ecto`) to include only rows where both 
#' volume (`V`) and stage (`STAGE`) are greater than zero (indicating they are alive).
#' @param ecto A list object (from NicheMapR's `ectotherm` function)
#'   containing components `debout` (developmental output) and `dates`.
#' @return A data frame of post-hatch DEB output, including only periods when
#'   the organism is alive. Adds a `dates` column corresponding to those time points.
#' @examples
#' alive_data <- get_alive(ecto)
#'
get_alive <- function(ecto){
    debout <- data.frame(ecto$debout)
    # Then, we can subset to just the post-hatch period
  debout_post_hatch <- debout[debout$V > 0 & debout$STAGE > 0, ] 
  debout_post_hatch$dates <-  ecto$dates[which(debout$V > 0 & debout$STAGE > 0)]
  return(debout_post_hatch)
}

#' @title get_environ: Extract environmental data from DEB model output
#' @description Extracts and formats the environmental data (`environ`) from a DEB model
#' object, appending the corresponding datetime information.
#' @param ecto A list object (from NicheMapR's `ectotherm` function) containing components `environ` (heat budget calculations) and `dates`
#' @return A data frame of environmental variables with an added `datetime` column.
#' @examples
#' env_data <- get_environ(ecto)
#'
get_environ <- function(ecto){
  environ <- data.frame(ecto$environ)
  environ$datetime <- ecto$dates
  return(environ)
}

#' @title get_environ_alive: Limit environmental data to period when organism is alive
#' @description Filters environmental data to match the period when the organism is alive
#' (based on `dates` from the DEB output). Useful for synchronizing survival
#' and environmental records.
#' @param sub_alive A data frame of post-hatch DEB output (from [get_alive()]),
#'   containing a `dates` column for the alive period.
#' @param sub_env A data frame of environmental data (from [get_environ()]),
#'   containing a `datetime` column.
#' @return A data frame of environmental conditions corresponding only to the
#'   period when the organism is alive.
#' @examples
#' env_alive <- get_environ_alive(sub_alive = alive_data, sub_env = env_data)
#'
get_environ_alive  <- function(sub_alive, sub_env){
  sub_env_alive <- data.frame(sub_env[which(sub_env$datetime %in% sub_alive$dates), ])
  return(sub_env_alive)
}

#' @title update_p_surv: Update survival probabilities with additional mortality sources
#' @description Adjusts survival probabilities in post-hatch DEB output to include additional
#' sources of mortality such as hatch mortality, inactivity mortality, activity-based
#' mortality, and thermal load stress mortality.
#' @param debout_post_hatch A data frame of post-hatch DEB output (from [get_alive()]),
#'   containing a column `P_SURV` with baseline survival probabilities.
#' @param m_h Numeric scalar; proportional mortality due to hatching (0–1).
#' @param m_i Numeric scalar; proportional mortality due to inactivity (0–1).
#' @param m_a Numeric scalar; proportional mortality due to activity (0–1).
#' @param m_tls Numeric vector; proportional mortality due to thermal load stress (0–1)
#'   for each time step (must match `nrow(debout_post_hatch)`).
#' @return A modified data frame with an additional column `P_SURV2` representing
#'   updated survival probabilities accounting for all mortality components.
#' @examples
#' debout_updated <- update_p_surv(debout_post_hatch, m_h = 0.1, m_i = 0.02, m_a = 0.03, m_tls = runif(nrow(debout_post_hatch), 0, 0.1))
#'
update_p_surv <- function(debout_post_hatch, m_h, m_i, m_a, m_tls){
  # Update survival at hatch
  debout_post_hatch$P_SURV2 <- debout_post_hatch$P_SURV * (1 - m_h) 

  # Update survival for additional mortality rates
  for(i in 2:nrow(debout_post_hatch)){
    debout_post_hatch$P_SURV2[i] <- debout_post_hatch$P_SURV2[i-1] * (1 - m_a) * (1 - m_tls[i]) * (1 - m_i)
    # Given TLS mortality is incremented each time step we have one less row, so add row by adding on the last value of P_SURV2
    if(i == nrow(debout_post_hatch)){
      debout_post_hatch$P_SURV2[i] <- debout_post_hatch$P_SURV2[i-1]
    }
  }

  debout_post_hatch$P_SURV2  <- round(debout_post_hatch$P_SURV2, 2)
  return(debout_post_hatch)
}

#' @title get_lifecycle_vitals: Extract life-cycle vitals from DEB model output
#' @description Derives key life-history parameters from DEB model output, including survival,
#' fecundity, longevity, generation time, and the population growth rate (`rmax`).
#' The function builds a stage-structured matrix and estimates the dominant
#' eigenvalue (λ) of the matrix to calculate `rmax`.
#' @param alive_debout A data frame of post-hatch DEB output (from [get_alive()]),
#'   with columns `STAGE`, `E_B`, and updated survival probabilities (`P_SURV2`).
#' @param ecto A DEB model object containing at least a `yearout` component with
#'   a column `FECUNDITY` (representing reproductive output). This is summarised by stage.
#' @return A list containing:
#' \describe{
#'   \item{lambda}{Dominant eigenvalue of the life table (finite rate of increase).}
#'   \item{rmax}{Intrinsic rate of increase (log of λ).}
#'   \item{longevity}{Longevity in years (based on time steps in `alive_debout`).}
#'   \item{generation_time}{Age at first reproduction (years).}
#'   \item{stage_str}{Stage-structured matrix.}
#' }
#' @examples
#' vitals <- get_lifecycle_vitals(alive_debout = alive_data, ecto = ecto)
#' vitals$rmax
#'
get_lifecycle_vitals <- function(alive_debout, ecto, type = STAGE){
  
  # Fecundity data
  life_cycle_open <- data.frame(ecto$yearout)
  
  # Calculate longevity and generation time
        longevity <- nrow(alive_debout) # hours
         E_B_diff <- alive_debout$E_B[2:longevity] - alive_debout$E_B[1:(longevity - 1)]
  age_first_repro <- which(E_B_diff < 0)[1]
  generation_time <- age_first_repro / 24 / 365

  # Create stage-specific survival and fecundity
   stage_str <- alive_debout  %>% 
                group_by({{type}}) %>% 
                summarise(P_SURV = mean(P_SURV2),
                             E_B = mean(E_B),
                       fecundity = ifelse(E_B == 0, 0, life_cycle_open$FECUNDITY/life_cycle_open$CLUTCHES)) # Reproductive buffer energy, can see we only have a single reproductive event.

    # Build the Leslie matrix. First row fecundnity across stages, sub-diagonal survival probabilities
    build_life_table <- matrix(0, nrow = nrow(stage_str), ncol = nrow(stage_str))
    build_life_table[1,] <- stage_str$fecundity

    # Now we need survival probabilities on the sub-diagonal. Because these are just hypothetical situations and only one life cyle realisation to be sure matrices are correct make sure that all survival is > 0 so that it allows for reproduction to occur at least once.
    for(i in 2:nrow(stage_str)){
      build_life_table[i, i - 1] <- stage_str$P_SURV[i - 1]
      if(build_life_table[i, i - 1] == 0){
        build_life_table[i, i - 1] <- 0.01
      }
    }
    build_life_table[nrow(build_life_table), ncol(build_life_table)] <- stage_str$P_SURV[nrow(build_life_table)]  

    if(build_life_table[nrow(build_life_table), ncol(build_life_table)] == 0){
      build_life_table[nrow(build_life_table), ncol(build_life_table)] <- 0.01
    }

   #Estimate rmax and other demographic parameters from eigenvalue approach
   life_table_eigen <- Re(eigen(build_life_table)$values[1])
               rmax_log <- log(life_table_eigen)
            lambda_dem <- popbio::lambda(build_life_table)  # double check
           gen_time_dem <- popbio::generation.time(build_life_table)
    net_reproductive_rate_dem <- popbio::net.reproductive.rate(build_life_table)
  
  return(list(lambda = round(life_table_eigen, 2),
              rmax = round(rmax_log, 2),
         longevity = longevity / 24 / 365,   # days
          net_reproductive_rate_dem = round(net_reproductive_rate_dem, 2),
        gen_time_dem = round(gen_time_dem, 2),
   generation_time = generation_time,
         matrix = build_life_table))
}

#' @title damage_accum
#' @description Calculates a cumulative percentage increase through time,
#' where each interval’s increment depends on the temperature.
#' @param time Numeric vector of time points (e.g., minutes).
#' @param temp Numeric vector of temperatures (°C), same length as `time`.
#' @param intercept Numeric scalar controlling the denominator (base-10 intercept).
#' @param slope Numeric scalar controlling how temperature affects the rate.
#' @param threshold Numeric scalar between 0 and 1 indicating the damage threshold.For example, if LT50 is used as the threshold, set `threshold = 0.5`. Used to convert heat injury to mortality probabilities
#'
#' @return A numeric vector of cumulative percentages of mortality (0–100).
damage_accum <- function(time, temp, intercept, slope, threshold = 0.5) {
    
    output = data.frame(damage_i = NA, cum_damage = NA, mortality_prob = NA)
  
  # Calculate the damage in a given increment
  for (i in 1:(length(time) - 1)) {
    # Calculate the accumulated damage increment for each time step
    acc = (100 * (time[i + 1] - time[i])) / (10^(slope * max(temp[i + 1], temp[i]) + intercept))
    output[i, "damage_i"] = acc
  }

  # Sum the increments cumulatively to get accumulated heat damage through time
  output$cum_damage[1] <- output$damage_i[1]
  for (j in 2:nrow(output)) {
    output$cum_damage[j] <- output$cum_damage[j - 1] + output$damage_i[j]
  }

  # Convert HI to survival probability at each time step given a threshold value used (50%, or 80% usually)
  output$mortality_prob = (output$cum_damage * threshold) / 100
    
  # Now check that mortality cum damage and probabilities do not exceed 1 (100%)
  output$mortality_prob[output$mortality_prob > 1] <- 1
  output$cum_damage[output$cum_damage > 100] <- 100
  output$temp <- temp[1:(length(time) - 1)]
  output$time_min <- time[1:(length(time) - 1)]

  return(output)
}


# par(mfrow = c(3, 1))
# plot(temp ~ time, data = damage, type = "l", xlim = c(109000, 120780))
# plot(cum_damage_repair ~ time, data = damage, type = "l", xlim = c(109000, 120780))
# plot(mortality_prob ~ time, data = damage, type = "l", xlim = c(109000, 120780))

#' @title damage_accum
#' @description Calculates cumulative heat damage through time, where each
#' interval’s increment depends on temperature, and converts this to
#' mortality probability via an exponential hazard model. The threshold
#' is defined such that when cumulative damage reaches 100 * threshold,
#' mortality is 50% (e.g., threshold = 0.5 for LT50).
#'
#' @param time Numeric vector of time points (e.g., minutes).
#' @param temp Numeric vector of temperatures (°C), same length as `time`.
#' @param intercept Numeric scalar controlling the denominator (base-10 intercept).
#' @param slope Numeric scalar controlling how temperature affects the rate.
#' @param threshold Numeric scalar between 0 and 1 indicating the damage
#'   threshold. For example, if LT50 is used as the threshold, set
#'   `threshold = 0.5`. When cum_damage reaches 100 * threshold, mortality = 0.5.
#'
#' @return A data.frame with columns:
#'   - damage_i: damage increment per interval (in arbitrary "% damage" units)
#'   - cum_damage: cumulative damage (same units as damage_i)
#'   - hazard: cumulative hazard H
#'   - mortality_prob: mortality probability (0–1)
#'   - temp: interval temperature (°C)
#'   - time_min: interval end time (minutes)
damage_accum <- function(time, temp, intercept, slope, threshold = 0.5) {

  if (length(time) != length(temp)) {
    stop("'time' and 'temp' must have the same length.")
  }
  if (length(time) < 2L) {
    stop("Need at least two time points to compute damage increments.")
  }

  # interval lengths
  dt <- diff(time)

  # use max temp over each interval
  temp_step <- pmax(temp[-1], temp[-length(temp)])

  # damage increment per interval (arbitrary "% damage" units)
  damage_i <- (100 * dt) / (10^(slope * temp_step + intercept))

  # cumulative damage (same units as damage_i)
  cum_damage <- cumsum(damage_i)

  # map to hazard so that cum_damage = 100 * threshold -> mortality = 0.5
  D50   <- 100 * threshold             # "damage" at LT50
  haz   <- log(2) * (cum_damage / D50) # cumulative hazard
  surv  <- exp(-haz)
  mort  <- 1 - surv

  # numerical safety
  mort[mort < 0] <- 0
  mort[mort > 1] <- 1

  out <- data.frame(
    damage_i       = damage_i,
    cum_damage     = cum_damage,
    hazard         = haz,
    mortality_prob = mort,
    temp           = temp_step,
    time_min       = time[-1]
  )

  return(out)
}

#' @title damage_repair_accum
#' @description Calculates cumulative heat damage through time with
#' temperature-dependent repair (Sharpe–Schoolfield form), then converts
#' this to mortality probability via an exponential hazard model.
#' The threshold is defined such that when cumulative damage reaches
#' 100 * threshold, mortality is 50% (e.g., threshold = 0.5 for LT50).
#'
#' @param time Numeric vector of time points (e.g., minutes).
#' @param temp Numeric vector of temperatures (°C), same length as `time`.
#' @param intercept Numeric scalar controlling the denominator (base-10 intercept)
#'   for damage accumulation.
#' @param slope Numeric scalar controlling how temperature affects the rate
#'   of damage accumulation.
#' @param threshold Numeric scalar between 0 and 1 indicating the damage
#'   threshold (e.g., 0.5 for LT50). When cum_damage_repair reaches
#'   100 * threshold, mortality = 0.5.
#' @param TA  Arrhenius activation parameter for repair (Kelvin).
#' @param TAL Low-temperature inactivation parameter for repair (Kelvin).
#' @param TAH High-temperature inactivation parameter for repair (Kelvin).
#' @param TL  Low-temperature inactivation midpoint (Kelvin).
#' @param TH  High-temperature inactivation midpoint (Kelvin).
#' @param TREF Reference temperature for repair rate (Kelvin).
#' @param kdot_ref Reference repair rate at `TREF` (per minute).
#'
#' @return A data.frame with columns:
#'   - damage_i: damage increment per interval (no repair)
#'   - repair_inc: repair amount per interval
#'   - net_damage_inc: net damage increment after repair
#'   - cum_damage_repair: cumulative net damage through time
#'   - hazard: cumulative hazard H
#'   - mortality_prob: mortality probability (0–1)
#'   - temp: interval temperature (°C)
#'   - time_min: interval end time (minutes)
damage_repair_accum <- function(time, temp,
                                intercept, slope,
                                threshold = 0.5,
                                TA, TAL, TAH, TL, TH, TREF, kdot_ref) {

  if (!is.numeric(temp)) stop("'temp' must be numeric (°C).")
  if (length(time) != length(temp)) {
    stop("'time' and 'temp' must have the same length.")
  }
  if (length(time) < 2L) {
    stop("Need at least two time points to compute damage increments.")
  }

  # check scalar parameters
  for (nm in c("TA", "TAL", "TAH", "TL", "TH", "TREF", "kdot_ref")) {
    val <- get(nm, inherits = FALSE)
    if (!is.numeric(val) || length(val) != 1L) {
      stop(sprintf("'%s' must be a numeric scalar.", nm))
    }
  }

  # Kelvin temps for repair
  TK <- temp + 273.15

  # Step 1: baseline damage increments (without repair)
  damage <- damage_accum(
    time      = time,
    temp      = temp,
    intercept = intercept,
    slope     = slope,
    threshold = threshold
  )

  n <- nrow(damage)

  # prepare columns for repair; overwrite hazard/mortality later
  damage$repair_inc        <- NA_real_
  damage$net_damage_inc    <- NA_real_
  damage$cum_damage_repair <- NA_real_

  # interval lengths and Kelvin temps per interval
  dt      <- diff(time)                          # minutes
  TK_step <- pmax(TK[-1], TK[-length(TK)])       # K, max over each interval

  # Step 2: compute repair and net damage per interval
  for (i in seq_len(n)) {
    T_K  <- TK_step[i]
    dt_i <- dt[i]

    # Sharpe–Schoolfield repair rate (per minute)
    num <- kdot_ref * exp(TA * (1 / TREF - 1 / T_K))
    den <- 1 + exp(TAL * (1 / T_K - 1 / TL)) +
              exp(TAH * (1 / TH  - 1 / T_K))
    repair_rate <- num / den

    # amount repaired during this interval (same units as damage_i)
    repair_inc <- repair_rate * dt_i

    # net damage increment after repair, truncated at zero
    damage_inc <- damage$damage_i[i]
    net_inc    <- max(damage_inc - repair_inc, 0)

    damage$repair_inc[i]     <- repair_inc
    damage$net_damage_inc[i] <- net_inc

    # cumulative net damage
    if (i == 1L) {
      damage$cum_damage_repair[i] <- net_inc
    } else {
      damage$cum_damage_repair[i] <- damage$cum_damage_repair[i - 1L] + net_inc
    }
  }

  # Step 3: map cumulative net damage to hazard and mortality
  D50   <- 100 * threshold
  haz   <- log(2) * (damage$cum_damage_repair / D50)
  surv  <- exp(-haz)
  mort  <- 1 - surv

  mort[mort < 0] <- 0
  mort[mort > 1] <- 1

  damage$hazard         <- haz
  damage$mortality_prob <- mort

  # ensure temp/time columns correspond to intervals
  damage$temp     <- pmax(temp[-1], temp[-length(temp)])
  damage$time_min <- time[-1]

  return(damage)
}

#' @title damage_repair_accum_state
#' @description Calculates cumulative heat damage through time with
#' temperature-dependent repair (Sharpe–Schoolfield form), tracking
#' the damage state over time and computing instantaneous hazard and
#' mortality probability at each time step.
#' The threshold is defined such that when cumulative damage reaches
#' 100 * threshold, mortality is 50% (e.g., threshold = 0.5 for LT50).
#' @param time Numeric vector of time points (e.g., minutes).
#' @param temp Numeric vector of temperatures (°C), same length as `time`.
#' @param intercept Numeric scalar controlling the denominator (base-10 intercept)
#'   for damage accumulation.
#' @param slope Numeric scalar controlling how temperature affects the rate
#'   of damage accumulation.
#' @param threshold Numeric scalar between 0 and 1 indicating the damage
#'   threshold (e.g., 0.5 for LT50). When cum_damage_repair reaches
#'   100 * threshold, mortality = 0.5.
#' @param TA  Arrhenius activation parameter for repair (Kelvin).
#' @param TAL Low-temperature inactivation parameter for repair (Kelvin).
#' @param TAH High-temperature inactivation parameter for repair (Kelvin).
#' @param TL  Low-temperature inactivation midpoint (Kelvin).     
#' @param TH  High-temperature inactivation midpoint (Kelvin).
#' @param TREF Reference temperature for repair rate (Kelvin).
#' @param kdot_ref Reference repair rate at `TREF` (per minute).
#' @return A data.frame with columns:
#'   - damage_i: damage increment per interval (no repair)
#'   - repair_inc: repair amount per interval     
#'  - net_damage_inc: net damage increment after repair
#'  - cum_damage_repair: cumulative net damage through time
#'  - hazard: instantaneous hazard h from current damage state
#'  - mortality_prob: mortality probability (0–1)
#'  - temp: interval temperature (°C)
#'  - time_min: interval end time (minutes)
#' 
damage_repair_accum_state <- function(time, temp,
                                      intercept, slope,
                                      threshold = 0.5,
                                      TA, TAL, TAH, TL, TH, TREF, kdot_ref) {

  if (!is.numeric(temp)) stop("'temp' must be numeric (°C).")
  if (length(time) != length(temp)) {
    stop("'time' and 'temp' must have the same length.")
  }
  if (length(time) < 2L) {
    stop("Need at least two time points to compute damage increments.")
  }

  # Kelvin temps
  TK <- temp + 273.15

  # interval lengths and step temperatures (max over interval)
  dt      <- diff(time)                         # minutes
  temp_st <- pmax(temp[-1], temp[-length(temp)])
  TK_st   <- pmax(TK[-1],   TK[-length(TK)])

  n <- length(dt)

  # damage rate per interval (without repair), "% damage per minute"
  # from your original damage_accum:
  #   damage_i = (100 * dt) / 10^(slope * T + intercept)
  # so:
  damage_rate <- 100 / (10^(slope * temp_st + intercept))

  # Sharpe–Schoolfield repair rate per interval (per minute)
  repair_rate <- vapply(seq_len(n), function(i) {
    T_K <- TK_st[i]
    num <- kdot_ref * exp(TA  * (1 / TREF - 1 / T_K))
    den <- 1 + exp(TAL * (1 / T_K - 1 / TL)) +
              exp(TAH * (1 / TH  - 1 / T_K))
    num / den
  }, numeric(1))

  # Raw increments for bookkeeping
  damage_i   <- damage_rate * dt
  repair_inc <- repair_rate * dt

  # State variables
  net_damage_inc    <- numeric(n)  # can be negative
  cum_damage_state  <- numeric(n)  # damage STATE with repair
  hazard_inst       <- numeric(n)  # instantaneous hazard from state
  mortality_prob    <- numeric(n)

  D50 <- 100 * threshold
  D_prev <- 0

  for (i in seq_len(n)) {
    # net increment in damage state this interval
    net_inc <- (damage_rate[i] - repair_rate[i]) * dt[i]  # can be < 0
    D_new   <- max(D_prev + net_inc, 0)                   # repair can undo damage

    # instantaneous hazard given CURRENT damage state
    h_i <- log(2) * (D_new / D50)

    # interpret this as "mortality if tested now" via your mapping
    mort_i <- 1 - exp(-h_i)

    net_damage_inc[i]   <- net_inc
    cum_damage_state[i] <- D_new
    hazard_inst[i]      <- h_i
    mortality_prob[i]   <- mort_i

    D_prev <- D_new
  }

  out <- data.frame(
    damage_i          = damage_i,
    repair_inc        = repair_inc,
    net_damage_inc    = net_damage_inc,
    cum_damage_repair = cum_damage_state,  # current damage state
    hazard            = hazard_inst,       # from state, not cumulative
    mortality_prob    = mortality_prob,    # can go up or down
    temp              = temp_st,
    time_min          = time[-1]
  )

  return(out)
}

# Load each grasshopper DEB model and extract life-cycle details
load(here("output", "microclim", "grasshopperspA_stomata_closed_Box1.RDa")) #ectoA
load(here("output", "microclim", "grasshopperB_stomata_open_Box1.RDa"))     #ectoB

# Define Common Parameters Applied to Both Species

# Baseline mortality rates to add to DEB. TLS and differences in hatch mortality will be added to these baseline rates differently for each species
  m_i <- 0     # mortality per hour inactive
  m_a <- 1e-4  # mortality per hour active

## ------------------## 
# Species A
## ------------------## 

# First, find out dates and times when animal is alive. We need these to extract the body temperatures AND the relevant aspects of the lifecyle

# Get the DEB output while alive
  debout_post_hatch_A <- get_alive(ectoA)

# Get the environ data while alive, which we use for TLS as we need the body temperature every hour
          environ_A <- get_environ(ectoA)
    environ_A_alive <- get_environ_alive(debout_post_hatch_A, environ_A)
    environ_A_alive <- environ_A_alive %>% mutate(time_min = c(0:(n()-1))*60) # time in minutes

# Thermal Load sensitivity parameters for species A. Note that CTmax and z are values derived from a regression between temperture (yaxis) and log time to 50% mortality (xaxis). See Arnold et al. 2025 for more details. We need to invert these values to get the intercept and slope for the damage accumulation function.
 ctmaxA <- 50            # CTmax1h, 50
     zA <- -1.8         # z (slope), negative indicates decrease in survival with time which it must do as animals don't get better the longer they are at stressful temps - try -0.9
alphaA  <- -ctmaxA/zA    # Need these for the log time accumulation function
  betaA <- 1 / zA        # Need these for the log time accumulation function

# Define temperatures and times
      tbA <- environ_A_alive$TC-3       # body temps, Celsius
    timeA <- environ_A_alive$time_min # time in minutes

# Using body temperatures through time, estimate the accumulation of physiological damage towards a threshold of 50% mortality under no repair. Note that threshold (50% = 0.5) is used to scale % HI to % mortality
damageA <- damage_accum(time = timeA, 
                       temp = tbA, 
                  intercept = alphaA, 
                      slope = betaA)

# Set Arrhenius model and repair rate parameters        
TA <- 14065 # Arrhenius temperature, K
TL <- 10.5 + 273.15 # Arrhenius temperature lower threshold, K
TH <- 28.5 + 273.15 # Arrhenius temperature upper threshold, K
TAL <- 50000 # Arrhenius temperature lower, K
TAH <- 100000 # Arrhenius temperature upper, K
TREF <- 28 + 273.15 # reference temperature, K
kdot_ref <- 0.005 # repair rate % per minute at ref temp 28 deg C

damageA_repair <- damage_repair_accum_state(time = timeA, 
                       temp = tbA, 
                  intercept = alphaA, 
                      slope = betaA,
                         TA = TA, 
                        TAL = TAL, 
                        TAH = TAH,
                         TL = TL, 
                         TH = TH, 
                       TREF = TREF, 
                   kdot_ref = kdot_ref)
#plot(damageA_repair$time_min, damageA_repair$mortality_prob, type = "l")
# Now that we have the probability of mortality from thermal load stress, we can update the survival probabilities in the DEB model output. Note that the DEB model already has some baseline mortality rates from the Gompertz ageing model, so we are adding additional sources of mortality.

m_h <- 0.05   # mortality at hatch

debout_post_hatch_A <- update_p_surv(debout_post_hatch_A, 
                                    m_h = m_h,    # mortality at hatch
                                    m_a = m_a,    # mortality for active hour
                                    m_i = m_i,    # mortality for inactive hour
                                    m_tls = damageA_repair$mortality_prob)    # Heat load mortality

# Now we can extract life-cycle details. Our life cycle model assumes all animals progress to the next stage that survive. Given we have multiple life stages we create a stage-structured matrix (or could be an age-structuted model – for example, based on 'day') to capture survival and fecundity at each stage, and quantify the probability of transitioning through each stage. From this, we can estimate population growth rate (rmax) and other life-cycle details. This could be changed to an age structure here by setting type = DAY and should be more comparable across species with different development times.

lifecycle_A <- get_lifecycle_vitals(debout_post_hatch_A, ectoA, type = STAGE) 

## ------------------## 
# Species B
## ------------------## 

debout_post_hatch_B <- get_alive(ectoB)
          environ_B <- get_environ(ectoB)
    environ_B_alive <- get_environ_alive(debout_post_hatch_B, environ_B) %>% mutate(time_min = c(0:(n()-1))*60) # time in minutes
   
# Thermal Load sensitivity parameters for species A. Note that CTmax and z are values derived from a regression between temperture (yaxis) and log time to 50% mortality (xaxis). See Arnold et al. 2025 for more details. We need to invert these values to get the intercept and slope for the damage accumulation function.
 ctmaxB <- 47              # CTmax1h
     zB <- -1            # z (slope), negative indicates decrease in survival with time
alphaB  <- -ctmaxB/zB      # Need these for the log time accumulation function
  betaB <- 1 / zB          # Need these for the log time accumulation function

# Define temperatures and times
      tbB <- environ_B_alive$TC-3 # body temps, Celsius
    timeB <- environ_B_alive$time_min # time in minutes

# Using body temperatures through time, estimate the accumulation of physiological damage towards a threshold of 50% mortality under no repair
damageB <- damage_accum(time = timeB, 
                        temp = tbB, 
                   intercept = alphaB, 
                       slope = betaB)

damageB_repair <- damage_repair_accum_state(time = timeB, 
                       temp = tbB, 
                  intercept = alphaB, 
                      slope = betaB,
                         TA = TA, 
                        TAL = TAL, 
                        TAH = TAH,
                         TL = TL, 
                         TH = TH, 
                       TREF = TREF, 
                   kdot_ref = kdot_ref)
# plot(damageB_repair$time_min, damageB_repair$mortality_prob, type = "l")
# Now that we have the probability of mortality from thermal load stress, we can update the survival probabilities in the DEB model output

  m_h <- 0.05   # mortality at hatch

debout_post_hatch_B <- update_p_surv(debout_post_hatch_B, 
                                    m_h = m_h,    # mortality at hatch
                                    m_a = m_a,    # mortality for active hour
                                    m_i = m_i,    # mortality for inactive hour
                                    m_tls = damageB_repair$mortality_prob)    # Heat load mortality

# Now we can extract life-cycle details
lifecycle_B <- get_lifecycle_vitals(debout_post_hatch_B, ectoB)

### Create a TLS figure for box 2

# Species A

multiplier <- 50  # to scale the two y axes
p1_tlA <- ggplot() +
annotate("rect",
           xmin = 1800*60, xmax = 2500*60,
           ymin = -Inf, ymax = Inf,
           fill = "red", alpha = 0.2) +
  labs(x = "Time (hours)", y = "Body Temperature (°C)") +
geom_line(aes(x = timeA, y = tbA), color = "#356920", size = 0.5) + ylim(0, 50) +
geom_line(aes(x = timeA, y = c(0, damageA_repair$mortality_prob) * multiplier), color = "red") +
scale_y_continuous(
    name = "Probability of Mortality (Thermal Load Stress)", 
    sec.axis = sec_axis(~ . / multiplier, name = "Body temperature (°C)")) +
  theme_classic() +
  theme( axis.text.x = element_blank(), legend.position = "none", legend.text = element_text(size = 14), axis.text = element_text(size = 14), axis.title = element_blank())  

# Species B
p2_tlB <- ggplot() +
annotate("rect",
           xmin = 1800*60, xmax = 2500*60,
           ymin = -Inf, ymax = Inf,
           fill = "red", alpha = 0.2) +
  labs(x = "Time (hours)", y = "Body Temperature (°C)") +
geom_line(aes(x = timeB, y = tbB), color = "#664229", size = 0.5) + ylim(0, 50) +
geom_line(aes(x = timeB, y = c(0, damageB_repair$mortality_prob) * multiplier), color = "red") + 
scale_y_continuous(
    name = "Probability of Mortality (Thermal Load Stress)", 
    sec.axis = sec_axis(~ . / multiplier, name = "Body temperature (°C)")) +
  theme_classic() +
  theme( axis.text.x = element_blank(), legend.position = "none", legend.text = element_text(size = 14), axis.text = element_text(size = 14), axis.title = element_blank())  

# Combine plots
ggsave(p1_tlA | p2_tlB + plot_layout(ncol = 1), width = 11.876543, height = 4.839506, dpi = 600, filename = here("figs", "fig_grasshopper_TLS_box2.png"))

# Note this is a trimed down version for the plot for the figure in main MS. It was edited in BioRender to add axis labels and panel labels.
p1_tlA | p2_tlB + plot_layout(ncol = 1)

4.2 Modelling population dynamics of the two grasshopper species including density dependence and interactions

We now have key vital rate information for each species that incorperates physiologically explicit processes governing growth, survival and reproduction. We can now use this information to model population dynamics of the two species including density dependence and interspecific interactions.

One way (of many) we could model the population dynamics of the two species is to use a discrete-time Ricker model that incorporates both intra and interspecific competition. The general form of the model is as follows:

\[ N_{t+1}^i = N_t^i \cdot \lambda_i^{\left(1 - \frac{N_t^i + \sum_{j \neq i} \alpha_{ij} N_t^j}{K_i}\right)} \tag{1}\]

where \(N_t^i\) is the population size of species \(i\) at time \(t\), \(\lambda_i\) is the intrinsic growth rate of species \(i\), \(K_i\) is the carrying capacity of species \(i\), and \(\alpha_{ij}\) is the competition coefficient that describes how much species \(j\) affects species \(i\). When \(\alpha_{ij}\) = 1, species \(j\) competes just as strongly as species \(i\) does with itself. When \(\alpha_{ij}\) < 1, then species \(j\) has a weaker effect on species \(i\) then conspecifics and \(\alpha_{ij}\) > 1, indicates that species \(j\) has a stronger effect on species \(i\) than conspeecifics. The sum runs over all species \(j\) that are not equal to \(i\), which means that the model can be extended to include any number of species within a community.

As with any model there are assumptions that need to be considered. In this case:

1. it again assumes discrete non-overlapping generations.
2. it assumes \(\lambda\), \(K\) and now \(\alpha_{ij}\) are constant over time, ignoring any effects of environmental stochasticity, changing interaction strengths between species or evolution. Again, these assumptions can be relaxed. For example, \(\alpha_{ij}\) can of course vary depending on the environment and could be viewed as a function in of itself. As such, we could, for example, couple our mircoclimate and biophysical model predictions to estimate how \(\alpha_{ij}\) varies as a function of temperature, water or even resource availability based on experimental data describing how \(\alpha_{ij}\) varies as across different environments.
3. it assumes that population growth is limited by the total effective density based on both intra and interspecific effects. It also assumes these act additively on the log scale.
4. it assumes a closed system with no emmigration or immigration.

For now, we will assume that the population is closed and that \(K\) is constant over time. This model allows for asymmetric interaction (i.e., \(\alpha_{ij} \neq \alpha_{ji}\)), and so, we’ll give each species its own competition coefficient. We will also assume that \(\lambda\) varies slightly over time and add this to our model to demonstrate how you can relax assumptions of \(\lambda\), \(K\) and \(\alpha_{ij}\) being constant over time.

We’ll need some functions to run the population dynamics model.

#' @title Ricker two species function
#' @description A simple Ricker model function to simulate population dynamics for multiple species
#' @param N0   Initial population size of species i
#' @param r_i    Intrinsic growth rate for species i
#' @param sd_r_i Standard deviation of the growth rate for species i
#' @param K_i    Carrying capacity for species i
#' @param alpha Competition coefficient between species 
#' @param N0_j   Initial population size of species j
#' @return     Next years population size
ricker_two <- function(N0_i, r_i, sd_r_i, K_i, N0_j, alpha){
  N0_i * exp(rnorm(1, r_i, sd_r_i)*(1 - ((N0_i + (alpha * N0_j))/K_i)))
}

## Create a function to do a stochastic simulation of the two species Ricker model
#' @title Two species ricker model simulation function
#' @description A simple two species Ricker model function to simulate population dynamics for two species
#' @param N0_i     Initial population size of species i
#' @param r_i      Intrinsic growth rate for species i
#' @param sd_r_i   Standard deviation of the growth rate for species i
#' @param K_i      Carrying capacity for species i
#' @param N0_j     Initial population size of species j
#' @param r_j      Intrinsic growth rate for species j
#' @param sd_r_j   Standard deviation of the growth rate for species j
#' @param K_j      Carrying capacity for species j
#' @param alpha_ij Competition coefficient for species j on species i
#' @param alpha_ji Competition coefficient for species i on species j
#' @param time     Time step
#' @param nsims    Number of replicates
#' @return         A matrix of population sizes at each time step for nsims

sim_ricker_two_sp <- function(N0_i, r_i, sd_r_i, K_i, N0_j, r_j, sd_r_j, K_j, alpha_ij, alpha_ji, time, nsims){
  # Create a vector to store the population size at each time step
  N_array_i <- array(NA, dim = c(time, nsims)) 
  N_array_j <- array(NA, dim = c(time, nsims))
  
  for(j in 1:nsims){
    # Create arrays to store population sizes and initilise with initial population sizes
        N_array_i[1,j] <- N0_i
        N_array_j[1,j] <- N0_j
    
    for(i in 2:time){
      # Calculate the population size in the next time step for species i
        N_array_i[i,j] <- ricker_two(N_array_i[i-1, j], r_i, sd_r_i, K_i, N_array_j[i-1, j], alpha_ji)

      # Calculate the population size in the next time step for species j
        N_array_j[i,j] <- ricker_two(N_array_j[i-1, j], r_j, sd_r_j, K_j, N_array_i[i-1, j], alpha_ij)
    }
  }
  # Return the population size at each time step for species i and j
  return(list(N_i = N_array_i, N_j = N_array_j))
}

set.seed(123) # Set seed for reproducibility

# Set up the parameters for the Ricker model
      N0_A <- 1000     # Initial population size, species A
      N0_B <- 500      # Initial population size, species B
    sd_r_A <- 0.15     # Standard deviation of the growth rate species A
    sd_r_B <- 0.15     # Standard deviation of the growth rate species B
       K_A <- 10000    # Carrying capacity species A
       K_B <- 10000    # Carrying capacity species B
r_A <- lifecycle_A$rmax # Intrinsic growth rate species A and B open watering scenerio 
r_B <- lifecycle_B$rmax # Intrinsic growth rate species A and B closed watering scenerio
  alpha_AB <- 0.8      # Competition coefficient. Assumed not to vary with environment
  alpha_BA <- 1.2      # Competition coefficient. Assumed not to vary with environment
    time <- 5         # Time steps
   nsims <- 1000       # Number of simulations

# Lets run the simulations

N_heat <- sim_ricker_two_sp( N0_i = N0_A, 
                              r_i = r_A, 
                           sd_r_i = sd_r_A, 
                              K_i = K_A, 
                             N0_j = N0_B, 
                              r_j = r_B, 
                           sd_r_j = sd_r_B, 
                              K_j = K_B, 
                         alpha_ij = alpha_AB, alpha_ji = alpha_BA, time, nsims)


# Create a data frame to store the results for each species for each scenerio
sims_pop <- plyr::ldply(lapply(N_heat, function(x){ 
                            pivot_longer(as.data.frame(x), 
                            cols = everything(), 
                        names_to = "sim", 
                       values_to = "N")  %>% 
               mutate(time = rep(1:time, each = nsims))}))

Code to plot the results and export part of the box figure below:

  # Create a summary of the simulations to make plotting easier
   sims_heat  <- sims_pop  %>% group_by(.id, time) %>%
  summarise(N_m = mean(N), sd = sd(N)) 

# Plot the results
p1 = ggplot(sims_heat, aes(x = time, y = N_m, colour = .id, fill = .id)) +
geom_ribbon(aes(ymin = N_m - sd, ymax = N_m + sd), alpha = 0.3, colour = NA) +
  geom_line() +
  labs(x = "Generation (time)", y = "Population size (N)") + theme_classic()  + scale_colour_manual(values = c("#356920", "#664229")) + scale_fill_manual(values = c("#356920", "#664229")) + theme(legend.position = "none", axis.text = element_text(size = 14), axis.title = element_blank())
  
# Save the plot
# ggsave(here("figs", "fig-ricker_twosp_box2.png"), p1, width = 6, height = 5,  dpi = 600)

More complex models that incorporate a multitude of species and their interactions can be built using the same principles [e.g., 6,7,8]. Information on interaction strength between species can be used to understand how species coexistence is expected to change under extreme heat scenarios.

5 Code for “Hot Solutions for a changing world”: Using a systems modelling approach to help develop solutions to extreme heat events

A systems modelling approach can also be used to design possible solutions by understanding better how ‘thermal interventions’ affect heat experienced by organisms and how it influences growth, survival and reproduction. For example, we can manipulate rainfall at a particular time to simulate artificial watering which might be a possible intervention used during extreme heat. The added water can be used to make calculations of heat exchange, temperatures and transpiration under extreme heat events. We can also optimise the amount of water to balance costs and benefits of watering.

5.0.1 Simulating the effect of watering on plant and insect body temperatures

5.0.1.1 Setting the scene

Let’s assume our plant species is of interest for conservation or agricultural reasons and we want to understand how watering during a heatwave might help it survive (Fig. 2). The challenge with watering is that often we are drought stressed in combination with the extreme heat [9]. As such, it’s ideal to reduce water use because this saves money and water. So, the questions become:

1. How much water do we need to protect the plant from extreme heat while also not overwatering?
2. In addition, how does this affect pest insects that might be residing on the plant?

Figure 2— Evaluating the impacts of ‘thermal interventions’ on temperatures experienced by plants and animals. Under extreme heat we might add water to the system to help keep plants cool, but such intervention may also have unexpected consequences on insects. Using a systems approach we can evaluate the impact of watering on plant and insect body temperatures. This can help us understand how to best manage the system to reduce the impacts of extreme heat on plants and insects.

Below, we load in our original microclimate object. We can now use this to simulate the effect of watering on the microclimate. First, lets just look at the rainfall that occurred over the period of interest (Fig. 3). The rainfall data is stored in the microclimate object and is provided as a daily value (rather than hourly as for temperature).

# Lets reload our microclimate object so we can use it again
  load(here("output", "microclim", "microclimate.RDa"))

# Get the dates from the microclimate object. We'll take it here because heat-budget objects are not saved with the dates
  dates <- micro$dates

# Lets plot out the actual rainfall data for the location. This is the rainfall data that was used in the microclimate model. 
            rainfall_actual <- data.frame(micro$RAINFALL)
      rainfall_actual$dates <- as.POSIXct(micro$dates2)
  
# Plot the rainfall data each day through time

  p1 <- ggplot(rainfall_actual) + 
              geom_line(aes(y = micro.RAINFALL, x = dates), col = "blue") +
              geom_hline(yintercept = 10, colour = "red", lty = 2) +
              labs(y = "Rainfall (mm)", x = "Date") + theme_classic() +
              theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + 
              scale_colour_viridis_d()

# Now a plot of ran over the heat wave even just to zoom in

   p2 <-  rainfall_actual %>% filter(dates > "2019-01-15" & dates < "2019-01-31")  %>% 
  ggplot() + 
              geom_line(aes(y = micro.RAINFALL, x = dates), col = "blue") +
              geom_hline(yintercept = 10, colour = "red", lty = 2) +
              labs(y = "Rainfall (mm)", x = "Date") + theme_classic() +
              theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + 
              scale_colour_viridis_d()

# Combine the two plots

p1 + p2 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A") +
  theme(plot.tag = element_text(size = 12, face = "bold"))

Figure 3— Rainfall events over time for A) over years and B) only during the heat wave event from 2019-01-15 to 2019-01-31.

It’s clear that we had very little water during the heatwave event. If this was an agriculturally important plant species, we might than as the question “how would the temperature of the plants change if we added water to the system?”. Remember, stomatal closure during the heat wave event meant that the leaf temperatures were much higher than the air temperature. Spraying water will cool the plants through evaportive cooling, but by how much can we expect it to do that? Well, we can use the microclimate model to simulate this. Lets set up our scenerios.

# Lets assume we are planning on adding 10mm or 20 mm of water to the system to see how much it makes a difference on plant leaf temperaures. 
   rainfall_2 <- 2  # mm water to add
   rainfall_5 <- 5  # mm water to add
  rainfall_10 <- 10 # mm water to add
  rainfall_20 <- 20 # mm water to add

# Specify the dates we are going to add water. 
    date_waters <- as.character(seq(as.POSIXct("2019-01-01"), as.POSIXct("2019-01-31"), by = "days"))

# Lets create a new vector of rainfall and then add it to the microclimate object. This is done by creating a new vector of rainfall and then adding it to the microclimate object.
      RAINFALL <- micro$RAINFALL                        # extract rainfall from previous sim
     RAINFALL2 <- RAINFALL                              # copy rainfall into a new vector
  RAINFALL2[as.character(micro$dates2) %in% date_waters] <- rainfall_2 # add 2mm of water
     RAINFALL3 <- RAINFALL                              # copy rainfall into a new vector
  RAINFALL3[as.character(micro$dates2) %in% date_waters] <- rainfall_5 # add 5mm of water
     RAINFALL4 <- RAINFALL                              # copy rainfall into a new vector
  RAINFALL4[as.character(micro$dates2) %in% date_waters] <- rainfall_10 # add 5mm of water
       RAINFALL5 <- RAINFALL                              # copy rainfall into a new vector
  RAINFALL5[as.character(micro$dates2) %in% date_waters] <- rainfall_20 # add 5mm of water

Now we have created two new rainfall vectors simulating different amounts of water added to the system during the heatwave event. Lets have a look at the new rainfall data.

# Add in rainfall manipulations to columns
rainfall_actual$rainfall_2 <- RAINFALL2
rainfall_actual$rainfall_5 <- RAINFALL3
rainfall_actual$rainfall_10 <- RAINFALL4
rainfall_actual$rainfall_20 <- RAINFALL5

# Re-arrange data for plotting
rainfall_actual2 <- rainfall_actual  %>% 
  select(dates, micro.RAINFALL, rainfall_2, rainfall_5, rainfall_10, rainfall_20) %>%
  pivot_longer(cols = c(micro.RAINFALL, rainfall_2, rainfall_5, rainfall_10, rainfall_20), 
               names_to = "Rainfall", 
               values_to = "Rainfall_mm") %>% 
  mutate(Rainfall = ifelse(Rainfall == "micro.RAINFALL", "Original Rainfall", 
                           ifelse(Rainfall == "rainfall_10", "Watering (10 mm)", 
                                  ifelse(Rainfall == "rainfall_20", "Watering (20 mm)", 
                                         ifelse(Rainfall == "rainfall_2", "Watering (2 mm)", 
                                                ifelse(Rainfall == "rainfall_5", "Watering (5 mm)", NA))))),
         Rainfall = factor(Rainfall, levels = c("Watering (2 mm)", "Watering (5 mm)", "Watering (10 mm)", "Watering (20 mm)", "Original Rainfall")))
# Plot the rainfall data each day through time
 water <-  ggplot(rainfall_actual2) + 
              geom_vline(xintercept = as.POSIXct("2019-01-01"), linetype = "dashed", color = "red") +
              geom_vline(xintercept = as.POSIXct("2019-01-31"), linetype = "dashed", color = "red") +
              geom_line(aes(y = Rainfall_mm, x = dates, col = Rainfall), show.legend = TRUE) +
              labs(y = "Rainfall (mm)", x = "Date") + theme_classic() +
              theme(                 legend.position = "inside",
                 legend.position.inside = c(0.25,0.8),
                 axis.title = element_text(size =20),
                 axis.text.x = element_blank(),
                 legend.title = element_text(size = 20),
                 legend.text = element_text(size = 20)) + 
              scale_colour_viridis_d()
 water

Figure 4— Rainfall events over time with two watering scenarios (10 mm and 20mm of water applied) during a heat wave event 2019-01-15 to 2019-01-31 (red dashed lines).

5.0.1.2 Calculating the different mircoclimates using the watering scenerios

We can now use the micro_silo function to run the microclimate model with the new rainfall data. We simply need to update the rainfall data in the rainhour argument. The rainhourly argument is set to 1 to indicate that we want to use hourly rainfall data.

# Create hourly values, but water at midnight for one hour (00:00-01:00) (midnight is used to avoid changes too rapid for the integrator to solve).
 rainhour_2 <- unlist(lapply(RAINFALL2, function(x) c(rep(0,23), x)))
 rainhour_5 <- unlist(lapply(RAINFALL3, function(x) c(rep(0,23), x)))
rainhour_10 <- unlist(lapply(RAINFALL4, function(x) c(rep(0,23), x)))
rainhour_20 <- unlist(lapply(RAINFALL5, function(x) c(rep(0,23), x)))

rain_list <- list(rainhour_2, rainhour_5, rainhour_10, rainhour_20)
names(rain_list) <- c("rainhour_2", "rainhour_5", "rainhour_10", "rainhour_20")

# Now lets run the microclimate model with the new rainfall data. We will do this for both 10mm and 20mm of water added to the system.

#Set up the location and parameters for the microclimate model again. 
    longlat = c(140.79970255974192, -34.159176277782564) 
     dstart = "01/01/2018" # Start year
    dfinish = "31/12/2019" # End year
   minshade = 40           # Minimum shade, %
   maxshade = 80           # Maximum shade, %
     Usrhyt = 0.2          # Height above ground, m 
 microclima = 1            # Use mircoclima 
         PC = -1500        # Water potential for stomatal closure
         SP = 10

for(i in 1:length(rain_list)){
  # run the microclimate model for different watering regimes and save the output.
         rain <- micro_silo(loc = longlat, 
                         dstart = dstart,
                        dfinish = dfinish, 
                         Usrhyt = Usrhyt, 
                     microclima = microclima,
                       minshade = minshade,
                       maxshade = maxshade, 
                             PC = PC, 
                             SP = SP, 
                     rainhourly = 1, 
                       rainhour = rain_list[[i]],
                          email = "daniel.noble@edu.au")

# Save the microclimate object
save(rain, file = here("output", "microclim", paste0(names(rain_list)[i], "_shade.RDa")))
}

Now we have the new microclimate objects. The original, 2, 5, 10 and 20mm watering events. We have access to temperatures in the air and soil as well as high shade (80%) and low shade (40%). Let’s how feed this into our plant heat-budget model to predict the leaf temperature of the plant.

5.0.1.3 Calculating the plant leaf temperature using the watering scenerios

# Now, run the plant heat-budget model using each of the three different microclimate objects.

 objects <- c("microclimate.RDa", "rainhour_2_shade.RDa", "rainhour_5_shade.RDa", "rainhour_10_shade.RDa", "rainhour_20_shade.RDa")
   names <- c("micro", rep("rain", 4)) # Name of the microclimate objects on load. Original mirciclimate object named "micro" and the rest "rain" because the object saved at each run of loop is "rain"

# Set leaf functional traits
         live <- 0      # don't simulate behaviour or respiration
         leaf <- 1      # use the stomatal conductance formulation for evaporation
         Ww_g <- 1      # wet weight, g
        shape <- 2      # 0=plate, 1=cylinder, 2=ellipsoid
      shape_b <- 0.0025 # ratio of b axis:a axis for ellipsoid 
      shape_c <- 0.1176 # ratio of c axis:a axis for ellipsoid 
      g_vs_ab <- 0.4    # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
      g_vs_ad <- 0.0    # leaf vapour conductance, adaxial (top of leaf), mol/m2/s
  g_vs_ab_max <- 0.4    # leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s
  g_vs_ad_max <- 0.0    # leaf vapour conductance, adaxial (top of leaf), mol/m2/s
  epsilon_sub <- 0.95   # emissivity of the substrate (-)
  epsilon_sky <- 1      # emissivity of the sky (-), accounted for already in 'TSKY' output from microclimate model so make it 1
    epsilon_L <- 0.97   # emissivity of the leaf (-)
      alpha_L <- 0.5    # solar absorptivity of the leaf (-)
       fatosk <- 0.5    # radiation configuration factor to sky (-)
       fatosb <- 0.5    # radiation configuration factor to substrate (-)
 conv_enhance <- 1.4    # convective enhancement factor (-)
     pct_cond <- 0      # percent of leaf conducting to the ground (%)
       postur <- 3      # leaf oriented horizontally, not tracking sun
          M_1 <- 0      # make metabolic rate zero

# Load in the plant heat budget
  load(here("output", "microclim", "plant_micro.RDa"))

# Loop through the objects and run the plant heat-budget model
for(i in 1:length(objects)){

# Load the microclimate object
   load(here("output", "microclim", objects[i]))
      micro <- get(names[i]) # get the name of the object and set it as micro as ectotherm uses object names micro

# First, lets extract the relevant information from the microclimate object to get plant stomatal conductance sorted.
      PDIFs <- micro$diffuse_frac             # use variable diffuse fraction
        P_a <- get_pressure(micro$elev)       # air pressure, Pa
        P_a <- rep(P_a, nrow(micro$metout))   # create hourly vector
      dates <- micro$dates            
 
# Set up vapour conductance vectors and simulate stomatal closure at night
                              f_open <- rep(1, nrow(micro$metout))
  f_open[micro$metout[,"ZEN"] == 90] <- 0      # close stomata when the sun is set
                            g_vs_abs <- g_vs_ab * f_open 
                            g_vs_ads <- g_vs_ad * f_open

# Include simulated stomatal conductance from plantecophys
  UMOLPERJ <- 4.57                                                   # mu_mol photons / J
WETAIR_out <- WETAIR(db = micro$metout[,"TALOC"], rh = micro$metout[,"RHLOC"], bp = P_a)
    VPD_Pa <- WETAIR_out$esat - WETAIR_out$e
 photo_out <- Photosyn(Tleaf = plant$environ[,"TC"],                 # Leaf temperature (C)
                          g1 = 5.2,                                  # Parameter of Ball-Berry model
                          Ca = 400,                                  # Atmospheric CO2 concentration (ppm)
                         VPD = VPD_Pa/1000,                          # Vapour pressure deficit (kPa)
                      vpdmin = 0.5,                                  # lower threshold on VPD
                        PPFD = micro$metout[,"SOLR"] * UMOLPERJ * 1, # mu mol m-2 s-1
                        Patm = P_a/1000                              # Atmospheric pressure (kPa)
                      )
g_vs_abs_pho <- photo_out$GS

# Run the heat-budget model for the plant leaf. Note that the microclimate object is used to get the relevant environmental temperatures to make the calculations of leaf temperature
    plant <- ectotherm(  leaf = leaf, 
                         live = live,
                         Ww_g = Ww_g, 
                        shape = shape, 
                      shape_b = shape_b, 
                      shape_c = shape_c, 
                      g_vs_ab = g_vs_abs_pho, 
                      g_vs_ad = g_vs_ads, 
                  g_vs_ab_max = g_vs_ab_max, 
                  g_vs_ad_max = g_vs_ad_max, 
                  epsilon_sub = epsilon_sub, 
                  epsilon_sky = epsilon_sky,
                      epsilon = epsilon_L, 
                    alpha_max = alpha_L, 
                    alpha_min = alpha_L,
                          M_1 = M_1,
                       fatosk = fatosk, 
                       fatosb = fatosb, 
                 conv_enhance = conv_enhance, 
                     pct_cond = pct_cond, 
                       postur = postur,
                       preshr = P_a, 
                         PDIF = PDIFs)
# Save the plant object
  save(plant, file = here("output", "microclim", paste0("plant_", gsub(".RDa", "", objects[i]), "_micro.RDa")))

}

5.0.1.4 Impacts of watering scenerios on leaf temperature

We’ve now used the heat budget model applying each of the three different microclimate objects to calculate the leaf temperature of the plant. We can now plot this out to see how the watering events affect leaf temperature.

# Extract the microclimate temperatures. Because this was simulated for a plant the 'body temperature' is the leaf temperature and is stored as 'TC' in the environ object

objects <- c("plant_microclimate_micro.RDa", "plant_rainhour_2_shade_micro.RDa", "plant_rainhour_5_shade_micro.RDa", "plant_rainhour_10_shade_micro.RDa", "plant_rainhour_20_shade_micro.RDa")
listnames <- c("envron_plant", "environ_plant_2", "environ_plant_5", "environ_plant_10", "environ_plant_20")

environ_plant <- list()
for(i in 1:length(objects)){
# Load the plant object
        load(here("output", "microclim", objects[i]))
        environ_plant[[i]] <- as.data.frame(plant$environ)
}
names(environ_plant) <- listnames

# Now, lets bind these data.frames together and re-organise the data for plotting
  environ_plant2 <- plyr::ldply(environ_plant) 

# Extract and stores dates
      environ_plant2$dates <- rep(dates, 5) 

# Now create a column that names each scenario
    environ_plant2 <- environ_plant2  %>% 
                     mutate(.id = ifelse(.id == "envron_plant", "Original Rainfall", 
                                        ifelse(.id == "environ_plant_2", "Watering (2 mm)", 
                                              ifelse(.id == "environ_plant_5", "Watering (5 mm)", 
                                                    ifelse(.id == "environ_plant_10", "Watering (10 mm)",  
                                                          ifelse(.id == "environ_plant_20", "Watering (20 mm)", NA))))))  %>% 
                            mutate(.id = factor(.id, levels = c("Original Rainfall", "Watering (2 mm)", "Watering (5 mm)", "Watering (10 mm)", "Watering (20 mm)")))

# Now lets plot the leaf temperature of the plant over time for each of the three different watering scenarios
plant_temps_water <- environ_plant2  %>% filter(dates > "2019-01-15" & dates < "2019-01-20")  %>% 
  ggplot() + 
              geom_line(aes(y = TC, x = dates, col = .id), show.legend = TRUE, size = 1) +
              labs(y = "Leaf Temperature (°C)", x = "Date") + theme_classic() +
              theme(legend.position = "none") + 
              scale_colour_viridis_d() + theme(legend.title = element_blank(),
                                                legend.text = element_blank(),
                                                axis.title = element_text(size =24),
                                                axis.text.x = element_blank(), 
                                                axis.text.y = element_text(size = 20)) 
# Export the figure for Fig 2 in paper
ggsave(plant_temps_water, filename = here("figs", "fig-plant_temp2.png"), width = 8.580247, height = 6.111111, dpi = 600)

Table 1— Mean leaf temperature differences between the original rainfall and the watering scenarios (2, 5, 10 and 20 mm of water added) during the heat wave event (2019-01-15 to 2019-01-31) between the hours 11:00 and 16:00 hr.

Watering Scenario	Mean Leaf Temperature Difference (°C)
Original Rainfall	0.00
Watering (2 mm)	-2.80
Watering (5 mm)	-3.63
Watering (10 mm)	-3.82
Watering (20 mm)	-3.87

?@fig-plant_temp3 shows how much watering at the various levels affects the leaf temperature of the plant during the day (between 11:00-14:00 hr). We can see that watering even just 2mm can reduce the leaf temperature by on average 2.8 \(^\circ\text{C}\) (Tab. 1). Interestingly, we can see that after 5mm of water, the leaf temperature doesn’t change much, suggesting that we can get away with less water. A good thing in a very arid environment!

5.0.2 Impacts of watering scenarios on grasshopper life cycle

5.0.2.1 Setting the scene

Now that we have the leaf temperature of the plant under the various watering scenarios, we can use this to understand how the different scenarios might affect grasshopper body temperature, life cycle and population dynamics like we did in Box 2. This is a nice demonstration of re-applying everything we’ve learnt to a new problem – a problem potentially involving the need to control and understand the probability of pest outbreaks.

We can do this by using the ectotherm model using the leaf temperature under the different watering scenarios as the microclimate for the grasshopper. Given the short period of the heat wave and watering regimes we don’t expect major differences but it’s better to check.

5.0.2.2 Use the different watering scenarios to look at how changes in temperature during the watering regime period affect the grasshopper life cycle

As we did before, we’ll update the microclimate experienced by the grasshopper to coincide with the leaf temperature and environmental conditions of the plant (assuming grasshoppers are within the bush in 40% shade). This will then calculate what temperatures the grasshopper experiences using a heat-budget model to predict the grasshopper body temperature. Then, we will use the heat-budget calculations to run the Dynamic Energy Budget model for the grasshopper to predict whether this impacts life-cycle. Again, these simulations are all done at once using the ectotherm function.

# We first need to load the DEB model paraemeters. These are estimated using MatLab. We've already estimated these so we'll simply just load them into R.

     pars <- readMat(here("sims", "results_Warramaba_virgo.mat"))

par.names <- unlist(labels(pars$par))
for(i in 1:length(par.names)){
  assign(par.names[i], unlist(pars$par[i]))
}

     pars <- readMat(here("sims", "stat_Warramaba_virgo.mat"))
par.names <- unlist(labels(pars$stat))
for(i in 1:length(par.names)){
  assign(par.names[i], unlist(pars$stat[i]))
}

# Additional parameters we need
        E_sm <- 2000            # J/cm^2, maximum stomach energy content
       fract <- 1               # scaling factor to change size

# Inital state for DEB model

# Egg Stage. Starting conditions
     E_init <- (E.0 * fract ^ 4) / (3E-9 * fract ^ 3)
   E_H_init <- 0
     V_init <- 3E-9
      stage <- 0

# Life history traits
 clutchsize <- 8                # Eggs per clutch
     stages <- 7                # Number of stages, 0 = egg, 1-5 = nymphs, 6 = pre egg-laying adult, 7 = post-laying adult
   S_instar <- rep(s.1, stages) # Stress at instar n: L_n^2/ L_n-1^2 (-)}\cr}
        L_b <- L.b              # Structural length at birth (cm)}\cr}
 photostart <- 0                # no photoperiod effect on reproduction
photofinish <- 0                # no photoperiod effect on reproduction
      reset <- 0                # don't reset lifecycle if death occurs
   soilnode <- 3                # depth (soil node) at which eggs laid - 3 = 5 cm

# Thermoregulatory parameters
       Ww_g <- 0.2              # g, body mass, not used though because DEB model computes this
      shape <- 1                # cylinder
    shape_b <- 10               # ratio of long to short axis of cylinder
     T_pref <- 38               # preferred body temperature (animal will attempt to regulate as close to this value as possible)
    T_F_min <- 13               # degrees C, lower body temperature for foraging
    T_F_max <- 40               # degrees C, upper body temperature for foraging
    T_B_min <- 13               # degrees C, minimum basking temperature
   T_RB_min <- 13               # degrees C, temperature at which animal will move to a basking site
     CT_max <- 47               # degrees C, critical thermal maximum (animal will die if CT_kill = 1 and this threshold is exceeded)
     CT_min <- 0                # degrees C, critical thermal minimum (used by program to determine depth selected when inactive and burrowing)


# Behavioural traits (made stage-specific later)
      diurn <- 1                # Diurnal activity allowed (1) or not (0)?
    nocturn <- 1                # Nocturnal activity allowed (1) or not (0)?
     crepus <- 1                # Crepuscular activity allowed (1) or not (0)?
     burrow <- 0                # Shelter in burrow allowed (1) or not (0)?
 shade_seek <- 1                # Shade seeking allowed (1) or not (0)?
   pct_cond <- 0                # Percentage of animal surface contacting the substrate (\%)

#' @title Ricker model function
#' @description A simple Ricker model function to simulate population dynamics
#' @param N0   Initial population size
#' @param r    Intrinsic growth rate
#' @param sd_r Standard deviation of the growth rate
#' @param K    Carrying capacity
#' @return     Next years population size
ricker <- function(N0, r, sd_r, K, time){
  N0 * exp(rnorm(1, r, sd_r)*(1 - (N0/K)))
}

## Create a function to do a stochastic simulation of the Ricker model
#' @title Ricker model simulation function
#' @description A simple Ricker model function to simulate population dynamics
#' @param N0    Initial population size
#' @param r     Intrinsic growth rate
#' @param sd_r  Standard deviation of the growth rate
#' @param K     Carrying capacity
#' @param time  Time step
#' @param nreps Number of replicates
#' @return      A matrix of population sizes at each time step for nreps

sim_ricker <- function(N0, r, sd_r, K, time, nreps){
  N_array <- array(NA, dim = c(time, nreps)) # Create a vector to store the population size at each time step
  for(j in 1:nreps){
          N_array[1,j] <- N0
    
    for(i in 2:time){
      N_array[i,j] <- ricker(N_array[i-1, j], r, sd_r, K)
    }
  }
  return(N_array)
}

# First, we need to load in the  microclimates scenerios. We'll do this in a loop so that we can run the DEB model for each of the scenerios. So, here we'll just create a few vectors that cell specific objects we need.
        plant_micro <- c("plant_microclimate_micro.RDa", "plant_rainhour_2_shade_micro.RDa", "plant_rainhour_5_shade_micro.RDa", "plant_rainhour_10_shade_micro.RDa", "plant_rainhour_20_shade_micro.RDa")
  plant_mirco_names <- c("envron_plant", "environ_plant_2", "environ_plant_5", "environ_plant_10", "environ_plant_20")

# We also need to load in the original microclimate object so that we can get the dates and other information from it. This is the original microclimate object that we used to simulate the plant microclimate. We'll also overide the `mirco` object so that the ectotherm model can use it.
  load(here("output", "microclim", "microclimate.RDa"))

# Now we can loop through the two microclimates and run the DEB model for each of them
for(i in 1:length(plant_micro)){

# Load the microclimate object
  load(here("output", "microclim", plant_micro[i]))

# Use the plant microclimate data to reset the original microclimate object. Note that we need to adjust both shade and nonshade conditions 
  
  # Unshaded microclimate
  micro$metout <- plant$metout
  micro$metout [,"TALOC"] <- plant$environ[,"TC"]

  # Shaded microclimate  
  microshadmet <- plant$shadmet
  micro$shadmet [,"TALOC"] <- plant$environ [,"TC"]
 
  # Soil microclimate
       micro$soil <- plant$soil
   micro$shadsoil <- plant$shadsoil
  micro$soilmoist <- plant$soilmoist
  micro$shadmoist <- plant$shadmoist


# Get dates from microclimate object
 dates <- micro$dates  # dates, hourly
dates2 <- micro$dates2 # dates, daily

# Set starting date by truncating microclimate model output to when grasshoppers are epected to oviposit eggs. This is the date when the life cylce of the insect is expected to start.
oviposition_date <- '2018-03-31 23:00:00'
            subs <- which( dates > as.POSIXct(oviposition_date))
           subs2 <- which(dates2 > as.POSIXct(oviposition_date))

# Now run the DEB life-cyle model for the grasshopper. See Case Study 2 for more details on the parameters used in the DEB model.

  ecto <- ectotherm( Ww_g = Ww_g,
                    shape = shape,
                  shape_b = shape_b,
                   T_pref = T_pref,
                  T_F_min = T_F_min,
                  T_F_max = T_F_max,
                  T_B_min = T_B_min,
                 T_RB_min = T_RB_min,
                   CT_min = CT_min,
                   CT_max = CT_max,
                    diurn = diurn,
                  nocturn = nocturn,
                   crepus = crepus,
               shade_seek = shade_seek,
                   burrow = burrow,
                 pct_cond = pct_cond,
                   nyears = ceiling(length(subs2)/365),
                   metout = micro$metout[subs,],
                     soil = micro$soil[subs,], 
                  soilpot = micro$soilpot[subs,],
                soilmoist = micro$soilmoist[subs,], 
                    humid = micro$humid[subs,],
                    tcond = micro$tcond[subs,], 
                  shadmet = micro$shadmet[subs,],
                 shadsoil = micro$shadsoil[subs,], 
                  shadpot = micro$shadpot[subs,],
                shadmoist = micro$shadmoist[subs,],
                shadhumid = micro$shadhumid[subs,],
                shadtcond = micro$shadtcond[subs,], 
                 rainfall = micro$RAINFALL[subs2],
                minshades = rep(unique(micro$minshade), length(subs)),
                maxshades = rep(unique(micro$maxshade), length(subs)),
                      DEB = 1, # turns DEB model on
               metab_mode = 1, # turns on hemimetabolous model, abp
                        z = z * fract,
                    del_M = del.M,
                     p_Xm = p.Am * 100 / 24 / kap.X,
                    kap_X = kap.X,
                        v = v / 24,
                      kap = kap,
                      p_M = p.M / 24,
                      E_G = E.G,
                    kap_R = kap.R,
                      k_J = k.J  / 24,
                     E_Hb = E.Hb * fract ^ 3,
                     E_Hp = (E.Hj+1e-5) * fract ^ 3,
                     E_Hj = E.Hj * fract ^ 3,
                      h_a = h.a / (24 ^ 2),
                      s_G = s.G,
                      E_0 = E.0 * fract ^ 4,
                     mu_E = mu.E,
                      d_V = d.V,
                      d_E = d.E,
                    T_REF = T.ref,
                      T_A = T.A,
                     T_AL = T.AL,
                     T_AH = T.AH,
                      T_L = T.L,
                      T_H = T.H,
                     E_sm = E_sm,
                   E_init = E_init,
                 E_H_init = E_H_init,
                   V_init = V_init,
                    stage = stage,
               clutchsize = clutchsize,
                   stages = stages,
                 S_instar = rep(s.1, stages),
                      L_b = L.b,
               photostart = photostart,
              photofinish = photofinish,
                    reset = reset,
                 soilnode = soilnode)
# Save the DEB model
  save(ecto, file = here("output", "microclim", paste0("grasshopper_", plant_mirco_names[i], "_DEB.RDa")))
}

Now that we have the DEB model for each of the watering scenerios, we can plot out the grasshopper life cycles for each scenerio (Fig. 5).

# Scenerios
   plant_mirco_names <- c("envron_plant", "environ_plant_2", "environ_plant_5", "environ_plant_10", "environ_plant_20")
   debs <- paste0("grasshopper_", plant_mirco_names, "_DEB.RDa")
   scenerios <- c("Original Rainfall", "Watering (2 mm)", "Watering (5 mm)", "Watering (10 mm)", "Watering (20 mm)")

# Now, lets load in the DEB model for each of the scenerios and plot them out. We'll save to a list
plot_list <- list()

for(i in 1:length(debs)){
# Load the DEB model
  load(here("output", "microclim", debs[i]))
  datetime <- dates[subs]
              sub <- which(datetime > as.POSIXct('2018-03-31') & datetime < as.POSIXct('2019-04-30'))
           debout <- data.frame(ecto$debout)
  debout$datetime <- datetime
    scale_factor  <-  10


  p <- ggplot() +
ylim(-5, 400) +
geom_line(data = debout, aes(x = datetime, y = WETMASS * 1000), col = "black") +
geom_line(data = debout %>% filter(STAGE == 0), aes(x = datetime, y = WETMASS * 1000), col = "darkgrey", linewidth =2) +
geom_line(data = debout %>% filter(STAGE == 1), aes(x = datetime, y = WETMASS * 1000), col = "lightgreen", linewidth =2) +
geom_line(data = debout %>% filter(STAGE == 2), aes(x = datetime, y = WETMASS * 1000), col = "green", linewidth =2) +
geom_line(data = debout %>% filter(STAGE == 3), aes(x = datetime, y = WETMASS * 1000), col = "tan", linewidth =2) +
geom_line(data = debout %>% filter(STAGE == 4), aes(x = datetime, y = WETMASS * 1000), col = "brown", linewidth =2) +
geom_line(data = debout %>% filter(STAGE == 5), aes(x = datetime, y = WETMASS * 1000), col = "orange", linewidth =2) +
geom_line(data = debout %>% filter(STAGE == 6), aes(x = datetime, y = WETMASS * 1000), col = "pink", linewidth =2) +
labs(title = scenerios[i], x = "Date") + 
scale_y_continuous(
    name = "Wet mass (mg)")+ theme_classic()

plot_list[[i]] <- p 

}

(plot_list[[1]] + plot_list[[2]]) / (plot_list[[3]] + plot_list[[4]]) / (plot_list[[5]])  + plot_annotation(tag_levels = "A") +
  theme(plot.tag = element_text(size = 12, face = "bold"))

Figure 5— Life-cycle of a grasshopper at a location near Renmark, South Australia. Life-cycle of the grasshopper is shown in under 5 different scenarios: A) Original rainfall and B) 2mm, C) 5mm, D) 10 mm and E) 20mm of watering to plants during a heat wave event. Note that the life-cycle is shown as the wet mass of the grasshopper (mg) over time (days) with different colours corresponding to different life stages. The life stages are: 0 = egg (‘grey’), 1 = first instar (‘forest green’), 2 = second instar (‘lime green’), 3 = third instar (‘tan’), 4 = fourth instar (‘brown’), 5 = fifth instar (‘orange’), 6 = adult (‘pink’).

Just like before we can get some key parameters from the DEB model to understand how the different watering scenerios affect the grasshopper life cycle.

Table 2— Predicted life-cycle parameters of a grasshopper at a location near Renmark, South Australia under different watering scenerios (2, 5, 10 and 20 mm of water added) during a heat wave event. The life-cycle parameters are: A) Time to first instar (days), B) Time to adult (days), C) Time to oviposition (days), D) Time to first egg hatch (days), E) Time to first instar hatch (days), F) Time to first adult hatch (days).

We can then look at the expected population growth rate of the grasshopper species under the different watering scenerios as we did in case study 2.

set.seed(123) # Set seed for reproducibility

# Set up the parameters for the Ricker model
        N0 <- 1000                                     # Initial population size
         r <- life_cycle[,"Intrinsic growth rate (r)"] # Intrinsic growth (DEB model) 
     names <- scenerios                                # Scenerios
      sd_r <- 0.15                                     # Standard deviation of the growth rate
         K <- 10000                                    # Carrying capacity
      time <- 5                                     # Time steps
     nsims <- 1000     # Number of simulations

# Lets run the simulations
sims <- list()
for(i in 1:length(r)){
  sims[[i]] <- sim_ricker(N0, r[i], sd_r, K, time, nsims)
  names(sims)[i] <- names[i]
}

# Re-organise the data for plotting
data_organised  <-  lapply(sims, function(x) {
  pivot_longer(as.data.frame(x), 
                            cols = everything(), 
                        names_to = "sim", 
                       values_to = "N")  %>% 
               mutate(time = rep(1:time, each = nsims))
})

# Plot the data for each scenerio
data_organised2  <- plyr::ldply(data_organised, .id = "scenario")  %>% mutate(scenerio_sim = interaction(scenario, sim)) 

# Create a summary of the simulations to make plotting easier
   data_organised2_final  <- data_organised2  %>% group_by(scenario, time) %>%
  summarise(N_m = mean(N), sd = sd(N)) 

# Now plot the simulations
plot_water_sims <- ggplot(data_organised2_final, aes(x = time, y = N_m, color = scenario)) +
geom_ribbon(aes(ymin = N_m - sd, ymax = N_m + sd, fill = scenario), alpha = 0.3, colour = NA) +
  geom_line() +
  labs(x = "Time (Generations)", y = "Population size (N)") +
    theme_classic() + scale_colour_viridis_d() + scale_fill_viridis_d() +
    theme(axis.title = element_text(size = 20),
          axis.text = element_text(size = 20),
           legend.position = "inside",
      legend.position.inside = c(0.25,0.8),
          legend.text = element_text(size = 20),
          legend.title = element_blank())


plot_water_sims

Figure 6— Expected population growth rate of a grasshopper at a location near Renmark, South Australia under different watering scenerios (2, 5, 10 and 20 mm of water added). Models assume that the watering scenerio is implemented each generation leading to a similar intrinsic population growth rate. Watering scenerios implemented at different times of the year are likely to have different effects on the population growth rate, but could be captured by simulating lice cyles for different years.

6 References

1.

Kearney, M. and Porter, W. (2009) Mechanistic niche modelling: Combining physiological and spatial data to predict species’ ranges. Ecol. Lett. 12, 334–350

2.

Kearney, M.R. et al. (2013) Balancing heat, water and nutrients under environmental change: A thermodynamic niche framework. Funct. Ecol. 27, 950–966

3.

Maclean, I.M.D. et al. (2019) Microclima: An r package for modelling meso‐ and microclimate. Methods Ecol. Evol. 10, 280–290

4.

Kearney, M.R. et al. (2020) A method for computing hourly, historical, terrain‐corrected microclimate anywhere on earth. Methods Ecol. Evol. 11, 38–43

5.

Duursma, R.A. (2015) Plantecophys–an R package for analysing and modelling leaf gas exchange data. PLoS One 10, e0143346

6.

Bimler, M.D. et al. (2018) Accurate predictions of coexistence in natural systems require the inclusion of facilitative interactions and environmental dependency. J. Ecol. 106, 1839–1852

7.

Bimler, M.D. et al. (2023) Estimating interaction strengths for diverse horizontal systems using performance data. Methods Ecol. Evol. 14, 968–980

8.

Bimler, M.D. et al. (2024) Plant interaction networks reveal the limits of our understanding of diversity maintenance. Ecol. Lett. 27, e14376

9.

Dai, A. (2013) Increasing drought under global warming in observations and models. Nature Climate Change 3, 52–58