This is the online supplementary material 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.
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.
1Introduction
In this supplementary material, we develop five hypothetical case studies to demonstrate how a systems modelling approach can be applied in practice. Our case studies focus on plant-insect interactions, and how these interactions are influenced by extreme heat events.
We provide theory and code with each of our case studies in the hopes that it stimulates easier application and development of our ideas. Models and code will of course always need to be modified to suite the specific system in question and should be validated against empirical data. As such, we view this tutorial as a starting point for helping readers understand how to implement a systems modelling approach to predict the impact of extreme heat on species assemblages.
2Load packages
Our case studies all use 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 tutorial download our GitHub repository, which contains all the code and data required to run the models. To locate the document that generated this tutorial navigate to 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 for this tutorial.
Code
# 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 installationpacman::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)
Throughout we’ll also use some helper functions. These can be added to your working environment by running the code below.
Code
#' @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_typeswitch(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_typeswitch(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:00stop("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.
Code
# 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.
3Predicting relevant microclimates to organisms: A walk through microclimate modelling
Broad-scale climate data is extremely coarse relative to the types of conditions experienced by organisms in their natural environment. Using biophysical models we can use coarse climate, topological, rainfall, elevation and spatial information to make predictions on what the microclimate a species is likely to experience, which can help us predict its body temperature (exposure). In combination with an understanding of the species’ biology we can simulate realistic scenarios that are relevant to the species in question.
As a first step, it’s important to be able to down-scale climatic data to microclimate conditions relevant to organisms. To model the microclimate, we can use the NicheMapR and the microclima packages. This section will demonstrate some ways that we can achieve this and make predictions for locations relevant to organisms of interest.
There are various climate datasets that can be used to predict the microclimate but for simplicity we’ll use weather data from 1889 to yesterday from https://www.longpaddock.qld.gov.au/silo/.
The code below shows how to predict various microclimates at a location close to Renmark, South Australia (Calperum Station). For example, it will compute temperatures at depths below ground (e.g., 0 and 50cm at or below the ground) and also above ground (20 cm and 2 m air temperature above ground). These might be relevant places for our organisms of interest, but we could be more thorough and simulate at other heights and locations if required.
Code
# 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, # Locationdstart ="01/01/2018", # Start yeardfinish ="31/12/2019", # End yearminshade =40, # Minimum shade, %maxshade =80, # Maximum shade, %Usrhyt =0.2, # Height above ground, m microclima =1, # Use mircoclima PC =-1500, # Water potential for stomatal closureSP =10, # Stability parameter for stomatal closureemail ="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"))
We can now plot the below ground temperatures to see the thermal profiles at different depths. Organisms can exploit this variation, especially if they can burrow (Fig. 1). We can see that surface temperatures are much more heterogeneous than those at 50 cm below ground. This is because the surface is more exposed to solar radiation and can heat up more quickly. Note that the temperatures assume the area is 40% exposed to the sun.
Code
# Load the microclimate objectload(here("output", "microclim", "microclimate.RDa"))# Explore microclimate data. For ease of plotting we'll extract soil temperatures at the 40% shade levels across the full set of dates we requested. soil_temp <-data.frame(micro$soil) # Note this is under 40% shade soil_temp$dates <- micro$dates # Add in the dates# Re-arrange for plotting soil_temp2 <- soil_temp %>%select(D0cm, D50cm, dates) %>%pivot_longer(cols =c(D0cm, D50cm), names_to ="Depth", values_to ="Temperature")# Plot soil tempsggplot(soil_temp2) +geom_line(aes(y = Temperature, x = dates, col = Depth), show.legend =TRUE) +labs(y ="Soil Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1)) +scale_colour_viridis_d()
Figure 1— Average soil temperature at a location near Renmark, South Australia. Soil temperatures are simulated from the surface (0 cm) to 50 cm below ground (and in minimum shade, 40%) and shown across time (hours) between 01/01/2018 to 31/12/2019.
If we want to use air temperatures because our organism does not burrow then we can extract air temperatures at different heights above ground (we have simulated 20 cm and 2 m) (Fig. 2). Of course, we could calculate temperatures for other heights as well, but assume this is the space most relevant to our system (Fig. 4), which we will detail below. We can also simulate microclimates under different scenarios in new calls to the micro_silo function but in our case we have managed to keep it simple and do everything in a single run.
Information on air temperatures at our two heights is also stored in our microclimate object from the simulation we ran above. To obtain our air temperatures, these are stored in a different place in our microclimate object. Again, we have access to these calculations under two shade levels (40% and 80%), but we’ll just use our lowest shade level for now.
Important
The microclimate object has a lot of relevant information that we can use in our simulations. Keep this in mind while exploring the documentation. We’ll come back to this later on.
Code
# Get the above ground temperatures at 20 cm and 2 m air_temp <-data.frame(micro$metout) air_temp$dates <- micro$dates
Now we can have a look at what the air temperature would be 20 cm and 2 m above ground (Fig. 2).
Figure 2— Average air temperature at a location near Renmark, South Australia. Air temperatures at 20 cm and 2 m above ground (and in minimum shade, 40%) are shown for a short times between 2019-01-15 and 2019-01-31.
We can immediately see from Fig. 2 that air temperatures closer to the ground are hotter than those higher up because of increased radiative heat transfer from the ground to air in close proximity and convective cooling higher in the air column.
Interestingly, between 2019-01-15 to 2019-01-31 we had a major heatwave event. Lets have a look at the temperatures during this period.
Code
air_temp2 %>%filter(dates >"2019-01-15"& dates <"2019-01-31") %>%ggplot(aes(y = Temperature, x = dates, col = Height)) +geom_line() +geom_hline(yintercept =40, colour ="red", lty =2) +labs(y ="Air Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1)) +scale_colour_viridis_d()
Figure 3— Average air temperature at a location near Renmark, South Australia. Air temperatures at 20 cm and 2 m above ground (and in minimum shade, 40%) are shown for a short times between 2019-01-15 and 2019-01-31.
Important
It is important to recognize that these are predictions. As with any model, they need field validation [6]. Nonetheless, they provide an important tool in our systems modelling framework to help us understand how heatwaves are likely to affect species interactions.
Now that we’ve explored the various temperatures predicted in our microclimate model, we can use this information to predict the body temperatures of organisms. We can also use this information to predict how these organisms will respond to extreme heat events.
4Developing a systems modelling approach to predict biological responses to extreme heat
4.1Case Study 1: Simulating plant and insect interactions under heatwaves
4.1.1 Setting the scene
Case study one explores how feedback between plants and insects can happen in response to extreme heat events. Here, we are interested in how the behaviour and body temperature of grasshoppers is influenced by changes in plant transpiration which can impact microhabitat temperatures of grasshoppers (Fig. 4). While this is a simple situation it’s a useful starting point to understand how a systems modelling approach can be built up to predict the impact of heatwaves on species interactions.
Figure 4— Grasshopper and plant interactions under heatwaves. Grasshoppers feeding on leaves of a bush in Calperum Station, South Australia. Complex leaf structure provides shade and cooler microclimate for grasshoppers. Leaves exposed to hot temperatures close their stomata to conserve water reducing transpiration and photosynthesis. Leaves that are not exposed to direct sunlight can keep their stomata open and continue photosynthesis while also transpiring and providing cooler microhabitats that grasshoppers can exploit through thermoregulatory behaviours. Importantly, stomatal behaviour is mediated by drought conditions. More stomata will be closed during drought to conserve water and this will impact the cooling effect of the plant on the grasshopper.
4.1.2 Predicting temperatures experienced by plants that form the grasshopper microclimate
To predict the impact of heatwaves on plant-insect interactions, we first need to predict how the microclimate that each organism experiences translates into body temperatures. Taking a systems approach means that we can account for the effects of other organisms on our focal species. To do this, we need to consider the biology of the species that co-exist.
In this case, the plants that the grasshopper consumes will influence the microclimate the grasshopper experiences (Fig. 4). Plants will also respond to the thermal environmental, and this will directly and indirectly influence the grasshopper. For example, plants will change their transpiration through the day and this can impact local temperatures experienced by grasshoppers (Fig. 4).
The microclimate model that we ran in Section 1 can now be used to predict how the plant and grasshopper will respond to the heatwave. Remember, we just used weather data to calculate air temperatures expected at 0.2 and 2m in both 40% and 80% shade. We can use this microclimate data to estimate leaf temperatures for a plant and then use the information on plant temperatures to feedback into the calculations of grasshopper body temperature.
Before we can estimate the temperatures for the grasshopper, we need to use a biophysical model to predict what the plant is doing. We can then use this ‘new’ microclimate data to simulate the growth of insects under different temperature conditions, changing what would happen if plants respond differently Fig. 4.
Let’s use our microclimate simulations to predict the plant leaf and air temperatures at the relevant heights to grasshoppers.
First, the traits of the leaf must be defined. An ellipsoid shape is used in this simulation. Specification of the mass and the ratios of the linear dimensions of the shape (here, semi-minor axes relative to the semi-major axis) together with the density (default value 1 g cm\({^3}\)) allows the geometric functions in the ectotherm model to compute the absolute dimensions and surface areas.
Evaporative exchange is driven by the stomatal conductance parameters, which can be specified for the abaxial and adaxial sides of the leaf as well as the overall base leaf conductance through the cuticle. The pct_wet parameter is set to 0.
Radiative exchange parameters can be set including the emissivities of the leaf and the substrate, the solar absorptivities of the leaf, the configuration factors to the sky and ground, the percentage of the leaf conducting to the ground and the orientation to the sun.
Code
# Load the microclimate objectload(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 modelsave(plant, file =here("output", "microclim", "plant_micro.RDa"))
Code
# Load the plant objectload(here("output", "microclim", "plant_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 environ_plant <-as.data.frame(plant$environ) environ_plant$dates <- micro$dates # For clarity 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-31") %>%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)"))# Plot it outggplot(environ_plant2) +geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend =TRUE) +labs(y ="Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1), legend.position ="top") +scale_colour_viridis_d()
Figure 5— Leaf temperature of a plant at a location near Renmark, South Australia. Leaf temperature is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) between 2019-01-15 and 2019-01-31.
We will now use the plantecophys package to calculate the leaf temperature and transpiration of the plant. This will allow us to understand how the plant is responding to the heatwave and how this might influence the grasshopper. The plantecophys package allows us to get more realistic stomatal behaviour via the Photosyn function. This code runs Photosyn assuming the leaf temperature calculated with the ectotherm model to start with.
Code
# Include simulated stomatal conductance from plantecophys UMOLPERJ <-4.57# mu_mol photons / JWETAIR_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 modelCa =400, # Atmospheric CO2 concentration (ppm)VPD = VPD_Pa/1000, # Vapour pressure deficit (kPa)vpdmin =0.5, # lower threshold on VPDPPFD = plant$metout[,"SOLR"] * UMOLPERJ *1, # mu mol m-2 s-1Patm = P_a/1000# Atmospheric pressure (kPa) )g_vs_abs_pho <- photo_out$GS
Now recalculate 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.
Code
# Rerun the models. Note that the paramaters as the same as 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 modelsave(plant2, file =here("output", "microclim", "plant2_micro.RDa"))
Now, lets have a look at the leaf temperature when stomatal closure behaviour is incorporated.
Code
# Load the plant objectload(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_plantecophys2 <-as.data.frame(plant2$environ) environ_plant_plantecophys2$dates <- micro$dates # For clarirty add in dates from microclimate that ecto model uses# Re-organise the data for plotting environ_plant_plantecophys <- environ_plant_plantecophys2 %>%filter(dates >"2019-01-15"& dates <"2019-01-31") %>% 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)ggplot(new_environ_plant2) +geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend =TRUE) +labs(y ="Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1),legend.position ="top", legend.text =element_text(size =6)) +scale_colour_viridis_d()
Figure 6— Leaf temperature of a plant at a location near Renmark, South Australia. Leaf temperature is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) between 2019-01-15 and 2019-01-31 for two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure.
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.
4.1.3 Coupling plant and animal models: Incorporating cooling effects of plant transpiration to estimate grasshopper body temperature
Now, lets run the ectotherm model for a grasshopper, but instead of using the microclimate object as it is, lets assume that the plant environment forms the microclimate of the grasshopper. This allows us to take in the microclimate and use it to calculate the body temperature of the grasshopper, and eventually understand how this affects growth, development etc within a dynamic energy budget model. The ectotherm model also allows us to incorporate relevant behavioural plasticity in response to extreme heat, which we’ll do.
Code
# First, lets fix up the microclimate object so that we can use it in the ectotherm model. We need to add in the leaf temperature as the microclimate for the grasshopper. There are a couple of ways we can conceptualise this. For example, we could assume that animals are in shade and therefore infrared radiative heat (`TSKYC`) is the temperature the grasshopper experiences from their environment. Alternatively, we could assume that when sitting on a leaf the grasshopper is in direct contact with the leaf and therefore the leaf temperature is the microclimate experienced by the grasshopper. We therefore need to treat the temperatures of the leaf as the environmental temperatures experienced by the grasshopper. Then, we feed this new microclimate into the grasshopper biophysical model to predict its body temperature. # Re-load the original micro object. We'll modify this nowload(here("output", "microclim", "microclimate.RDa")) # Get some key mircoclimate details of the plants from the model that considers the behaviour of the stomata and update the mircoclimate. We'll do this for the 40% shade condition (metout).# First, assume that the air temp at the location is close to or equivilant to the temp of the plant leaf. Update miroclimate accordingly that will be used for grasshopper. micro$metout[,"TALOC"] <- environ_plant_plantecophys2[,"TC"]# Second, lets update the relative humidity. micro$metout[,"RHLOC"] <- environ_plant_plantecophys2[,"RELHUM"]# Third, lets also update the IR profiles. micro$metout[,"TSKYC"] <- environ_plant_plantecophys2[,"TSKY"] # Parameters for the grasshopper 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 (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 <-12.5# degrees C, minimum basking temperatureT_RB_min <-12.5# 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. This tells the model what the grasshopper can do in response to heat. diurn <-1# diurnal activity allowed (1) or not (0)? nocturn <-0# nocturnal activity allowed (1) or not (0)? crepus <-1# crepuscular activity allowed (1) or not (0)?shade_seek <-1# shade seeking allowed (1) or not (0)? burrow <-0# shelter in burrow allowed (1) or not (0)? pct_cond <-0# Percentage of animal surface contacting the substrate# Now that we have set the parameters we can run the modelgrasshopper <-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)# Save the grasshopper objectsave(grasshopper, file =here("output", "microclim", "grasshopper_micro.RDa"))
Now that we have run the ectotherm model we can plot the body temperature of the grasshopper over time along side the air temperature. This will allow us to see how the grasshopper responds to the heatwave. We can also see activity patterns of the grasshopper which predict how we expect it to behave in relation to its body temperature and time of day. At night the animal is inactive while during suitable periods of the day where body temperatures are suitable it will forage (Fig. 7).
Code
# Load grasshopper objectload(here("output", "microclim", "grasshopper_micro.RDa"))# Get the resulting heat budget calculations for the grasshopper grasshopper_environ <-as.data.frame(grasshopper$environ)grasshopper_environ$dates <- micro$dates# Filter to heat wave datesgrasshopper_environ2 <- grasshopper_environ %>%filter(dates >"2019-01-15"& dates <"2019-01-31") %>%select(dates, TA, TC, ACT) %>%pivot_longer( cols =c(TA, TC), names_to ="Temps", values_to ="Temperature") %>%mutate(Temps =ifelse(Temps =="TA", "Air Temperature", "Body Temperature"))ggplot(grasshopper_environ2) +geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend =TRUE) +geom_point(aes(y = ACT, x = dates, fill =factor(ACT)), shape =21, alpha =0.3 ,size =3) +labs(y ="Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1)) +scale_colour_viridis_d() +scale_fill_manual(values =c("0"="blue", "1"="green", "2"="orange"),labels =c("0"="Inactive", "1"="Basking", "2"="Foraging")) +labs(fill ="Behaviour")
Figure 7— Body temperature of a grasshopper at a location near Renmark, South Australia. Body temperature of the grasshopper is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) for a short times between 2019-01-15 and 2019-01-31.
4.2Case Study 2: Understanding the effects of leaf mirochabitat on insect life-cyles: Linking mircoclimates to dynamic energy budget and population ecology models
4.2.1 Setting the scene
Our next case study will explore how the microclimate produced by a plant leaf can influence the life-cycle of an insect. In this case, we will use a dynamic energy budget model to understand how the microclimate of a plant leaf can influence the growth and development of a grasshopper. Here, we will introduce a Dynamic Energy Budget (DEB) model which now takes as input the mircoclimate produced by the plant to predcit body temperatures, development, growth and reproduction of the grasshopper (see Fig. 8 for a simplified DEB model). The DEB model also allows us to incorperate the heat wave event.
Dynamic Energy Budget (DEB) models are a powerful tool for understanding how organisms allocate energy to growth, reproduction and maintenance. To parameterise a DEB model it requires a set of parameters that are derived from empirical data on growth and metabolism. While DEB models are parameter rich, key parameters have been estimated for > 3000 species already [7,8]. New analytical approaches using phylogenetic imputation methods can be used to predict DEB parameters for over 1.3 million species providing immense opportunities to model how heat will affect life-cycles of multi-species assemblages (see https://deb.bolding-bruggeman.com, [9]).
Figure 8— A simple adult stage Dynamic Energy Budget (DEB) model for the grasshopper showing how food gets assimilated into reserve energy stores. The grasshopper than aloocates some fraction of the energy in reserves to maintaining body stucrture (somatic maintenance) and growth (i.e., structure) and to developing reproductive tissue (reproductive maintenance) and allocating to reproduction.
4.2.2 Grasshopper development and reproduction under heatwaves when plants don’t respond
To setup our DEB model we first need to set the DEB model paarameters. Grasshoppers are slightly more complicated because of the differing life-stages but we can use an existing DEB model for the grasshopper Warramaba virgo, and use a model developed by [10]. Thankfully, [10] have already estimated the DEB parameters for us so we’ll use these. The DEB model parameters are available in the sims folder of the repository. We’ll load these parameters into R and then set the initial conditions for the DEB model.
Code
# 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 in1: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 in1: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 reproductionphotofinish <-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 (\%)
4.2.3 Microclimate data for the DEB model
Recall from Case Study 1 that we have already simulated the microclimate for the grasshopper assuming: 1) plant stoma remain open during the day or 2) close during the day. As we seen in Fig. 17, this had major implications for temperatures experienced in the plant microhabitat that grasshoppers inhabit. Under these two different scenerios, the key questions now become:
How do the different microhabitats created by plant stomatal behaviour influence the development, growth and reproduction of the grasshopper?
How do different vital rates impact population dynamics?
Answers to these questions are useful to see how a systems approach can be implemented. Lets first simulate the life-cyle of grasshoppers under the two different microclimate scenarios.
Code
# First, we need to load in the two microclimates scenerios. 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", "plant2_micro.RDa") plant_mirco_names <-c("stomata_open", "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"))# Now we can loop through the two microclimates and run the DEB model for each of themfor(i in1:length(plant_micro)){# Load the microclimate objectload(here("output", "microclim", plant_micro[i]))# Object is labelled plant2 for the plant2_micro.RDa objectif(plant_micro[i] =="plant2_micro.RDa"){ plant <- plant2 } else { plant <- plant }# 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, hourlydates2 <- 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 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 = micro$minshade[subs2],maxshades = micro$maxshade[subs2],DEB =1, # turns DEB model onmetab_mode =1, # turns on hemimetabolous model, abpz = 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 modelsave(ecto, file =here("output", "microclim", paste0("grasshopper_", plant_mirco_names[i], "_DEB.RDa")))}
Now that we have run the DEB model we can have a look at whether the life-cycles of the grasshoppers are different under the two different microclimate scenarios.
Code
# Load in the DEB model for scenerio 1load(here("output", "microclim", "grasshopper_stomata_open_DEB.RDa")) 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# Create the life cycle plotp1 <-ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col ='lightgrey', linewidth =0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col ="forestgreen") +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 ="Predicted Grasshopper life-cycle under stomata open", x ="Date") +scale_y_continuous(name ="Wet mass (mg)", sec.axis =sec_axis(~ . / scale_factor, name ="Body temperature (°C)") )+theme_classic()# Now for scenerio 2load(here("output", "microclim", "grasshopper_stomata_closed_DEB.RDa")) debout <-data.frame(ecto$debout) debout$datetime <- datetime# Create the life cycle plotp2 <-ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col ='lightgrey', size =0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col ="forestgreen") +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 ="Predicted Grasshopper life-cycle under stomata closed", x ="Date") +scale_y_continuous(name ="Wet mass (mg)", sec.axis =sec_axis(~ . / scale_factor, name ="Body temperature (°C)") )+theme_classic()p1/p2 +plot_layout(ncol =1) +plot_annotation(tag_levels ="A") +theme(plot.tag =element_text(size =12, face ="bold"))
Figure 9— Life-cycle of a grasshopper at a location near Renmark, South Australia. Life-cycle of the grasshopper is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure. 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’).The grey lines indicate the body temperature of the grasshopper (°C) (right axis).
We can see from Fig. 9 that the life-cyle of the grasshopper is different when we consider the behaviour of the plant stomata. Notably, when stomata are closed and the leaf temperatures are higher we expect the grasshopper to move through stages faster and reproduce more (Tab. 1). We can have a look at the different traits expected under both scenarios.
Code
# Traits we want to extract from the DEB model traits <-c("FECUNDITY", "CLUTCHES", "LifeSpan", "rmax") # Load in the DEB model for scenerio 1load(here("output", "microclim", "grasshopper_stomata_open_DEB.RDa")) life_cycle_open <-data.frame(ecto$yearout)[,traits]# Load in the DEB model for scenerio 2load(here("output", "microclim", "grasshopper_stomata_closed_DEB.RDa"))life_cycle_closed <-data.frame(ecto$yearout)[,traits]# Finalise the table.table_final <-rbind(life_cycle_open, life_cycle_closed)table_final <- table_final %>%mutate(Scenario =c("Stomata Open", "Stomata Closed")) %>%select(Scenario, everything())# Create a table of the life-cycle traitsflextable::flextable(table_final) %>% flextable::set_header_labels(values =c("FECUNDITY"="Fecundity", "CLUTCHES"="Clutches", "LifeSpan"="Life Span","Lamda Max"="rmax")) %>% flextable::set_table_properties(layout ="autofit")
Table 1— Predicted life-cycle traits of a grasshopper at a location near Renmark, South Australia. Life-cycle traits of the grasshopper are shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure.
Scenario
Fecundity
Clutches
Life Span
rmax
Stomata Open
24
3
1.112329
0.4188939
Stomata Closed
24
3
1.090868
1.0360276
4.2.4 Population dynamics of grasshoppers under two plant response scenerios
As we have seen in the previous section, the life-cycle of the grasshopper is different when we consider the behaviour of the plant stomata and this has noticable effects on vital rates. We can use this information to parameterise a simple stochastic population models to understand how population dynamics of the grasshopper might change.
We can use a simple Ricker model [11], which we’ll use to simulate the population dynamics of the grasshopper. The Ricker model is a simple model that describes how populations grow over time. The model is given by the equation:
where \(N_t\) is the population size at time \(t\), \(r\) is the intrinsic growth rate, and \(K\) is the carrying capacity of the environment. The model assumes that populations grow exponentially when they are small, but that growth slows down as the population approaches the carrying capacity.
We can also re-parameterise to use a finite rate model given that the grasshoppers are an annual species.
While Equation 2 is a useful model it does make some simplifying assumptions. First, it assumes discrete non-overlapping generations. This should be a sensible assumption for insect and probably many plant species. Second, \(\lambda\) and \(K\) are assumed constant over time, ignoring any effects of environmental stochasticity or evolution. However, this assumption can easily be relaxed. Third, it assumes a closed system with no emmigration or immigration. However, again, we can modify the model to include these factors. For now, we will assume that the population is closed and that \(\lambda\) and \(K\) are constant over time.
Recall from our previous section that we already estimated the intrinsic growth rate. Using this information, for the stomata-closed scenario (\(\lambda\) = 2.82), the population could tolerate up to 64.51% mortality per generation and still remain stable.
In contrast, for the stomata-open scenario (\(\lambda\) = 1.52), the population is increasing more slowly. The population can tolerate around 34.22% mortality per generation and still remain stable.
Assuming we had information on the carying capacity of the enviroment and population size we can then use the \(r_{max}\) in the Ricker model to simulate the population dynamics of the grasshopper.
Code
#' @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 sizericker <-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 nrepssim_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 stepfor(j in1:nreps){ N_array[1,j] <- N0for(i in2:time){ N_array[i,j] <-ricker(N_array[i-1, j], r, sd_r, K) } }return(N_array)}
Code
set.seed(123) # Set seed for reproducibility# Set up the parameters for the Ricker model N0 <-1000# Initial population sizer_closed <-1.0360276# Intrinsic growth (DEB model) with stoma closed r_open <-0.4188939# Intrinsic growth (DEB model) with stoma open sd_r <-0.15# Standard deviation of the growth rate K <-10000# Carrying capacity time <-10# Time steps nsims <-1000# Number of simulations# Lets run the simulationsN_closed <-sim_ricker(N0, r_closed, sd_r, K, time, nsims) N_open <-sim_ricker(N0, r_open, sd_r, K, time, nsims)# Create a data frame to store the resultssims_closed <-pivot_longer(as.data.frame(N_closed), cols =everything(), names_to ="sim", values_to ="N") %>%mutate(time =rep(1:time, each = nsims))sims_open <-pivot_longer(as.data.frame(N_open), cols =everything(), names_to ="sim", values_to ="N") %>%mutate(time =rep(1:time, each = nsims))
Code
open_stoma <-image_read(here("figs", "Stomata_open.png")) open <-rasterGrob(open_stoma)closed_stoma <-image_read(here("figs", "Stomata_closed.png")) closed <-rasterGrob(closed_stoma)# Plot the resultsp1 =ggplot(sims_closed, aes(x = time, y = N, colour = sim)) +geom_line() +labs(title ="Stomata Closure", x ="Generation (time)", y ="Population size (N)") +theme_classic() +theme(legend.position ="none") +scale_colour_viridis_d() +theme(axis.text =element_text(size =12), axis.title =element_text(size =14)) +annotation_custom(closed, xmin =5, xmax =9, ymin =500, ymax =8500) +ylim(0, 12000)p2 =ggplot(sims_open, aes(x = time, y = N, colour = sim)) +geom_line() +labs(title ="Stomata Remain Open", x ="Generation (time)", y ="Population size (N)") +theme_classic() +theme(legend.position ="none") +scale_colour_viridis_d() +theme(axis.text =element_text(size =12), axis.title =element_text(size =14)) +annotation_custom(open, xmin =1, xmax =5, ymin =7000, ymax =11000) +ylim(0, 12000)p1 +theme(plot.tag =element_text(size =12, face ="bold")) + p2 +theme(plot.tag =element_text(size =12, face ="bold")) +plot_layout(ncol =2) +plot_annotation(tag_levels ="A")
Figure 10— Population dynamics of a grasshopper at a location near Renmark, South Australia. Population dynamics of the grasshopper is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure.
We can see from Fig. 10 that when stomata are able to respond by closing grasshoppers are expected to grow rapidly until they reach carrying capacity (K = 10,000, and the point they all converge on), whereas if we have stomata remaining open we expect the population to grow much more slowly over time as microhabiats become less suitable for growth and reproduction. This happens because higher body temperatures that result from when stomata close better allows grasshoppers to reach optimal physiological performance.
4.3Case Study 3: Heat wave impacts through small communities
4.3.1 Setting the scene
A systems modelling approach naturally involves attempting to model multi-species assemblies / communities. This is challenging with large communities because of the complexity of the interactions between species. However, for simpler communities such an approach of connecting models across fields is achievable and potentially powerful in making predictions about how communities will respond to extreme heat events.
In this case study, we expand on the impacts of extreme heat on a single species of grasshopper and show how this can be extended to capture interactions between three species: a plant, a grasshopper and a second competing grasshopper species that occupies the same mircohabitat. We extend our Ricker model to capture the impacts that each grasshopper species has on the other. One way to model this is to explicity translate how each species’ interactions impact the population size of the other.
In keeping with our notation, we can extend this model to a two-species Ricker model to include the interactions between the two grasshopper species as follows:
where \(N_t^A\) and \(N_t^B\) are the population sizes of species A and B at time \(t\), \(\lambda_A\) and \(\lambda_B\) are the intrinsic growth rates of species A and B, \(K_A\) and \(K_B\) are the carrying capacities of species A and B, and \(\alpha\) is the competition coefficient that describes how much species B affects species A, whereas \(\beta\) is the competition coefficient that describes how much species A affects species B.
In Equation 4, the competition coefficient \(\alpha\) is a measure of the strength of the interaction between the two species. If \(\alpha\) is greater than 1, it means that species B has a stronger negative effect on species A than species A has on itself. If \(\alpha\) is less than 1, it means that species A has a stronger negative effect on itself than species B has on it. If \(\alpha\) is equal to 1, it means that the two species have equal effects on each other. The opposite is true for \(\beta\) in Equation 5.
A more general form of this model is given by: \[
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{6}\]
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.
Just like Equation 2, Equation 6 is also making similar and a few additional assumptions that wre worth noting.
it again assumes discrete non-overlapping generations.
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.
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.
it assumes a closed system with no emmigration or immigration.
while we assume symmetric interactions, it allows for asymmetric interaction (i.e., \(\alpha_{ij} \neq \alpha_{ji}\))
For now, we will assume that the population is closed and that \(K\) is constant over time. While this model allows for asymmetric interaction (i.e., \(\alpha_{ij} \neq \alpha_{ji}\)), we will assume that the interactions are symmetric (i.e., \(\alpha_{ij} = \alpha_{ji}\)). 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.
4.3.3 Simulating effects of extreme heat on two species
We can simulate the life-cyle of a second grasshopper species that occupies the same habitat but varies in physiology. We’ve already done this for one species Warramaba virgo (Species B; Fig. 11) in Case Study 2.
Figure 11— Simulating two species competitive interactions.
The second species (Species A; Fig. 11) has a different physiology compared to the first species but is usually more abundant than Species B making it highly competitive. We’ll now simulate its life-cycle like we did for Species B in Case Study 2. Then, link this to a two-sepcies population dynamic model. These species are closely related but do have some differences in physiology that include: 1) higher thermal tolerance, 2) larger clutch sizes, and 3) different preferred body temperatures; and 4) smaller size. Otherwise, we will assume the same DEB model parameters as Species B.
Code
# 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 in1: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 in1: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 <- L.b # Structural length at birth (cm)}\cr} photostart <-0# no photoperiod effect on reproductionphotofinish <-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 <-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 <-1# Shade seeking allowed (1) or not (0)? pct_cond <-0# Percentage of animal surface contacting the substrate (\%)
Code
# First, we need to load in the two microclimates scenerios. 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", "plant2_micro.RDa") plant_mirco_names <-c("stomata_open", "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"))# Now we can loop through the two microclimates and run the DEB model for each of themfor(i in1:length(plant_micro)){# Load the microclimate objectload(here("output", "microclim", plant_micro[i]))# Object is labelled plant2 for the plant2_micro.RDa objectif(plant_micro[i] =="plant2_micro.RDa"){ plant <- plant2 } else { plant <- plant }# 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, hourlydates2 <- 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 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 = micro$minshade[subs2],maxshades = micro$maxshade[subs2],DEB =1, # turns DEB model onmetab_mode =1, # turns on hemimetabolous model, abpz = 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 =0.0779, # Gompertz stress coefficient, from DebberE_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 =0.09, # structural length at birth (cm),photostart = photostart,photofinish = photofinish,reset = reset,soilnode = soilnode)# Save the DEB modelsave(ecto, file =here("output", "microclim", paste0("grasshopperspA_", plant_mirco_names[i], "_DEB.RDa")))}
We can have a look at the life-cycle of the new species of grasshopper, Species A, under both scenarios (Fig. 12).
Code
# Load in the DEB model for scenerio 1load(here("output", "microclim", "grasshopperspA_stomata_open_DEB.RDa")) 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# Create the life cycle plotp1 <-ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col ='lightgrey', linewidth =0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col ="forestgreen") +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 ="Predicted Species A Grasshopper life-cycle under stomata open", x ="Date") +scale_y_continuous(name ="Wet mass (mg)", sec.axis =sec_axis(~ . / scale_factor, name ="Body temperature (°C)") )+theme_classic()# Now for scenerio 2load(here("output", "microclim", "grasshopperspA_stomata_closed_DEB.RDa")) debout <-data.frame(ecto$debout) debout$datetime <- datetime# Create the life cycle plotp2 <-ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col ='lightgrey', size =0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col ="forestgreen") +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 ="Predicted Species A Grasshopper life-cycle under stomata closed", x ="Date") +scale_y_continuous(name ="Wet mass (mg)", sec.axis =sec_axis(~ . / scale_factor, name ="Body temperature (°C)") )+theme_classic()p1/p2 +plot_layout(ncol =1) +plot_annotation(tag_levels ="A") +theme(plot.tag =element_text(size =12, face ="bold"))
Figure 12— Life-cycle of a grasshopper (Species A) at a location near Renmark, South Australia. Life-cycle of the grasshopper is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure. 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’).The grey lines indicate the body temperature of the grasshopper (°C) (right axis).
We can now put the summaries of all our life-cycle predictions together for both species of grasshopper and both scenerios (Tab. 2). We’ll use the intrinsic growth rates in our population models.
Table 2— Life-cycle of two grasshopper species at a location near Renmark, South Australia. Life-cycle of the grasshoppers is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure.
4.3.4 Two-species population dynamics
Modelling the dynamics between two species requires some knowledge about how competition impacts the population size of each species. Normally, this is achieved by estimating competition coefficients from field or laboratory data. Experiments involve determining how the presence of one species impacts growth of the other relative to a control state where only the single species is present. Given that our goal is simply to demonstrate how models from different fields can be coupled, we’ll assume that the competition coefficients are known and equal so that \(\alpha_{i,j}\) = 0.2.
Lets re-write our functions to use the two-species model now.
Code
#' @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_ij Competition coefficient for species j on species i#' @param N0_j Initial population size of species j#' @return Next years population sizericker_two <-function(N0_i, r_i, sd_r_i, K_i, N0_j, alpha_ij){ N0_i *exp(rnorm(1, r_i, sd_r_i)*(1- ((N0_i + (alpha_ij * 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 time Time step#' @param nsims Number of replicates#' @return A matrix of population sizes at each time step for nsimssim_ricker_two_sp <-function(N0_i, r_i, sd_r_i, K_i, N0_j, r_j, sd_r_j, K_j, alpha_ij, 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 in1: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_jfor(i in2: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_ij)# 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 jreturn(list(N_i = N_array_i, N_j = N_array_j))}
We can now use these functions to simulate the population dynamics for both species under the two different scenarios. For simplity we have assumed the competition coefficient is the same for both species and does not vary with environment. Under more realistic conditions it’s very likely that the competition coefficient would vary with environment and we could use our microclimate model to estimate this.
Code
set.seed(123) # Set seed for reproducibility# Set up the parameters for the Ricker model N0_i <-1000# Initial population size, species A N0_j <-500# Initial population size, species B sd_r_i <-0.15# Standard deviation of the growth rate species A sd_r_j <-0.15# Standard deviation of the growth rate species B K_i <-10000# Carrying capacity species A K_j <-10000# Carrying capacity species Br_i_j_open <- life_cycle[c(1,3),"Intrinsic growth rate (r)"] # Intrinsic growth rate species A and B open watering scenerio r_i_j_closed <- life_cycle[c(2,4),"Intrinsic growth rate (r)"] # Intrinsic growth rate species A and B closed watering scenerio alpha_ij <-0.8# Competition coefficient. Assumed not to vary with environment time <-10# Time steps nsims <-1000# Number of simulations# Lets run the simulationsN_closed <-sim_ricker_two_sp(N0_i, r_i = r_i_j_closed[1], # Species i growth rate sd_r_i, K_i, N0_j, r_j = r_i_j_closed[2], # Species j growth rate sd_r_j, K_j, alpha_ij, time, nsims)N_open <-sim_ricker_two_sp(N0_i, r_i = r_i_j_open[1], # Species i growth rate sd_r_i, K_i, N0_j, r_j = r_i_j_open[2], # Species j growth rate sd_r_j, K_j, alpha_ij, time, nsims)# Create a data frame to store the results for each species for each sceneriosims_closed <- plyr::ldply(lapply(N_closed, function(x){ pivot_longer(as.data.frame(x), cols =everything(), names_to ="sim", values_to ="N") %>%mutate(time =rep(1:time, each = nsims))}))sims_open <- plyr::ldply(lapply(N_open, function(x){ pivot_longer(as.data.frame(x), cols =everything(), names_to ="sim", values_to ="N") %>%mutate(time =rep(1:time, each = nsims))}))
Now we can plot the results of the simulations for both species under both scenarios (Fig. 13).
Figure 13— Simulated population dynamics for two species of grasshopper under two different watering scenarios. A) Stomata open, B) Stomata closed. The lines show the mean population size across 1000 simulations and the shaded area shows the standard deviation.
More complex models that incorporate a multitude of species and their interactions can be built using the same principles [e.g., 12,13,14]. Information on interaction strength between species can be used to understand how species coexistence is expected to change under extreme heat scenarios.
5Using 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.1Case Study 4: Simulating the effect of watering on plant and insect body temperatures
5.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. 14). The challenge with watering is that often we are drought stressed in combination with the extreme heat [15]. As such, it’s ideal to reduce water use because this saves money and water. So, the questions become:
How much water do we need to protect the plant from extreme heat while also not overwatering?
In addition, how does this affect pest insects that might be residing on the plant?
Figure 14— 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. 15). The rainfall data is stored in the microclimate object and is provided as a daily value (rather than hourly as for temperature).
Code
# Lets reload our microclimate object so we can use it againload(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 plotsp1 + p2 +plot_layout(ncol =1) +plot_annotation(tag_levels ="A") +theme(plot.tag =element_text(size =12, face ="bold"))
Figure 15— 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 (Fig. 6). 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.
Code
# 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.
Code
# Add in rainfall manipulations to columnsrainfall_actual$rainfall_2 <- RAINFALL2rainfall_actual$rainfall_5 <- RAINFALL3rainfall_actual$rainfall_10 <- RAINFALL4rainfall_actual$rainfall_20 <- RAINFALL5# Re-arrange data for plottingrainfall_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 16— 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.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.
Code
# 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 =10for(i in1: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 objectsave(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.1.3 Calculating the plant leaf temperature using the watering scenerios
Code
# 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 budgetload(here("output", "microclim", "plant_micro.RDa"))# Loop through the objects and run the plant heat-budget modelfor(i in1:length(objects)){# Load the microclimate objectload(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 / JWETAIR_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 modelCa =400, # Atmospheric CO2 concentration (ppm)VPD = VPD_Pa/1000, # Vapour pressure deficit (kPa)vpdmin =0.5, # lower threshold on VPDPPFD = micro$metout[,"SOLR"] * UMOLPERJ *1, # mu mol m-2 s-1Patm = 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 objectsave(plant, file =here("output", "microclim", paste0("plant_", gsub(".RDa", "", objects[i]), "_micro.RDa")))}
5.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.
Code
# 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 objectobjects <-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 in1:length(objects)){# Load the plant objectload(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 scenariosplant_temps_water <- environ_plant2 %>%filter(dates >"2019-01-15"& dates <"2019-01-31") %>%ggplot() +geom_line(aes(y = TC, x = dates, col = .id), show.legend =TRUE) +labs(y ="Leaf Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1), legend.position ="top") +scale_colour_viridis_d() +theme(legend.title =element_blank(),legend.text =element_text(size =12),axis.title =element_text(size =18),axis.text =element_text(size =12)) plant_temps_water
Figure 17— Leaf temperature of a plant at a location near Renmark, South Australia. Leaf temperature is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) between 2019-01-15 and 2019-01-31 for three different watering scenarios (original rainfall, 2, 5, 10 and 20 mm of water added).
Table 3— 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. 17 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. 3). 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.2Case Study 5: Impacts of watering scenerios on grasshopper life cycle
5.2.1 Setting the scene
Now that we have the leaf temperature of the plant under the various watering scenerios, we can use this to understand how the different scenerios might affect grasshopper body temperature, life cycle and population dynamics like we did in case study 2. This is a nice demonstrations 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 scenerios 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.2.2 Use the different watering scenerios 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 tempeartures the grasshopper experiences using a heat-budget model to predict the grasshopper body temperature. Then, we will use the heat-budget calculations to run teh 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.
Code
# 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 themfor(i in1:length(plant_micro)){# Load the microclimate objectload(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, hourlydates2 <- 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 = micro$minshade[subs2],maxshades = micro$maxshade[subs2],DEB =1, # turns DEB model onmetab_mode =1, # turns on hemimetabolous model, abpz = 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 modelsave(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. 18).
Code
# 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 listplot_list <-list()for(i in1:length(debs)){# Load the DEB modelload(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 18— 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 4— 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.
Code
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 <-10# Time steps nsims <-1000# Number of simulations# Lets run the simulationssims <-list()for(i in1:length(r)){ sims[[i]] <-sim_ricker(N0, r[i], sd_r, K, time, nsims)names(sims)[i] <- names[i]}# Re-organise the data for plottingdata_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 sceneriodata_organised2 <- plyr::ldply(data_organised, .id ="scenario") %>%mutate(scenerio_sim =interaction(scenario, sim)) plot_water_sims <-ggplot(data_organised2, aes(x = time, y = N, col = scenario)) +geom_line(aes(group = scenerio_sim)) +labs(x ="Time (Generations)", y ="Population size (N)") +theme_classic() +scale_colour_viridis_d() +theme(axis.title =element_text(size =20),axis.text =element_text(size =20),legend.position ="inside",legend.position.inside =c(0.7,0.5),legend.text =element_text(size =20),legend.title =element_blank())plot_water_sims
Figure 19— 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.
De Frenne, P. et al. (2024) Ten practical guidelines for microclimate research in terrestrial ecosystems. Methods Ecol. Evol.
7.
Marques, G.M. et al. (2018) The AmP project: Comparing species on the basis of dynamic energy budget parameters. PLoS Comput. Biol. 14, e1006100
8.
Kooijman, S.A.L.M. et al. (2021) Multidimensional scaling for animal traits in the context of dynamic energy budget theory. Conserv. Physiol. 9, coab086
9.
Bruggeman, J. et al. (2009) PhyloPars: Estimation of missing parameter values using phylogeny. Nucleic Acids Research 37, W179–W184
10.
Llandres, A.L. et al. (2015) A dynamic energy budget for the whole life-cycle of holometabolous insects. Ecol. Monogr. 85, 353–371
11.
Ricker, W.E. (1954) Stock and recruitment. J. Fish. Res. Board Can. 11, 559–623
12.
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
13.
Bimler, M.D. et al. (2023) Estimating interaction strengths for diverse horizontal systems using performance data. Methods Ecol. Evol. 14, 968–980
14.
Bimler, M.D. et al. (2024) Plant interaction networks reveal the limits of our understanding of diversity maintenance. Ecol. Lett. 27, e14376
15.
Dai, A. (2013) Increasing drought under global warming in observations and models. Nature Climate Change 3, 52–58
Source Code
---title: "A systems-modelling approach to predict biological responses to extreme heat"date: "`r Sys.Date()`"bibliography: ../bib/refs.bibcsl: ../bib/trends-in-ecology-and-evolution.cslformat: html: toc: true toc-location: left toc-depth: 3 toc-title: "**Table of Contents**" output-file: "index.html" theme: cosmo embed-resources: true code-fold: show code-tools: true number-sections: true fontsize: "12" max-width: "10" code-overflow: wrapcrossref: fig-title: Figure # (default is "Figure") tbl-title: Table # (default is "Table") title-delim: — # (default is ":") fig-prefix: Fig. # (default is "Figure") tbl-prefix: Tab. # (default is "Table")editor_options: chunk_output_type: console---This is the online supplementary material 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.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.# **Introduction**In this supplementary material, we develop five hypothetical case studies to demonstrate how a systems modelling approach can be applied in practice. Our case studies focus on plant-insect interactions, and how these interactions are influenced by extreme heat events. We provide theory and code with each of our case studies in the hopes that it stimulates easier application and development of our ideas. Models and code will of course always need to be modified to suite the specific system in question and should be validated against empirical data. As such, we view this tutorial as a starting point for helping readers understand how to implement a systems modelling approach to predict the impact of extreme heat on species assemblages.# **Load packages**Our case studies all use the statistical software [R](https://www.r-project.org/). We use a number of packages to run our simulations, including [`NicheMapR`](https://mrke.github.io)[@Kearney2009b; @Kearney2013], [`microclima`](https://github.com/ilyamaclean/microclima)[@Maclean2019; @Kearney2020], and [`plantecophys`](https://github.com/RemkoDuursma/plantecophys)[@Duursma2015]. We recommend that readers interested in working through this tutorial download our [GitHub repository](https://github.com/daniel1noble/thermal_tol_interactions), which contains all the code and data required to run the models. To locate the document that generated this tutorial navigate to 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 for this tutorial.```{r, libraries}#| label: libraries#| echo: true#| warning: false#| message: false# 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 installationpacman::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)``````{r, func}#| label: func#| echo: false#| warning: false#| message: false#| source(here("R", "func.R"))```Throughout we'll also use some helper functions. These can be added to your working environment by running the code below.```{r, functions}#| label: functions#| echo: true#| warning: false#| message: false#| eval: false#' @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'.") )}```:::{.column-margin}:::{.callout-note}gFortran can be a problem for NicheMapR so you need to set the paths properly if you are getting errors. Follow details [here](https://stackoverflow.com/questions/69639782/installing-gfortran-on-macbook-with-apple-m1-chip-for-use-in-r). The following code will help you set the paths properly.```{r, nichmapr}#| label: nichmapr#| echo: true#| eval: false#| results: hide#| code-fold: true# 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.```::::::# **Predicting relevant microclimates to organisms: A walk through microclimate modelling** {#section1}Broad-scale climate data is extremely coarse relative to the types of conditions experienced by organisms in their natural environment. Using biophysical models we can use coarse climate, topological, rainfall, elevation and spatial information to make predictions on what the microclimate a species is likely to experience, which can help us predict its body temperature (*exposure*). In combination with an understanding of the species' biology we can simulate realistic scenarios that are relevant to the species in question.As a first step, it's important to be able to down-scale climatic data to microclimate conditions relevant to organisms. To model the microclimate, we can use the [`NicheMapR`](https://mrke.github.io) and the [`microclima`](https://github.com/ilyamaclean/microclima) packages. This section will demonstrate some ways that we can achieve this and make predictions for locations relevant to organisms of interest. There are various climate datasets that can be used to predict the microclimate but for simplicity we'll use weather data from 1889 to yesterday from https://www.longpaddock.qld.gov.au/silo/. The code below shows how to predict various microclimates at a location close to Renmark, South Australia (Calperum Station). For example, it will compute temperatures at depths below ground (e.g., 0 and 50cm at or below the ground) and also above ground (20 cm and 2 m air temperature above ground). These might be relevant places for our organisms of interest, but we could be more thorough and simulate at other heights and locations if required. ```{r, microclimate_soil}#| label: microclimate_soil#| echo: true#| warning: false#| message: false#| error: false#| eval: false # 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"))```We can now plot the below ground temperatures to see the thermal profiles at different depths. Organisms can exploit this variation, especially if they can burrow (@fig-soil_temp). We can see that surface temperatures are much more heterogeneous than those at 50 cm below ground. This is because the surface is more exposed to solar radiation and can heat up more quickly. Note that the temperatures assume the area is 40% exposed to the sun.```{r, fig-soil_temp}#| label: fig-soil_temp#| fig-cap: "Average soil temperature at a location near Renmark, South Australia. Soil temperatures are simulated from the surface (0 cm) to 50 cm below ground (and in minimum shade, 40%) and shown across time (hours) between 01/01/2018 to 31/12/2019. "#| code-fold: true# Load the microclimate object load(here("output", "microclim", "microclimate.RDa"))# Explore microclimate data. For ease of plotting we'll extract soil temperatures at the 40% shade levels across the full set of dates we requested. soil_temp <- data.frame(micro$soil) # Note this is under 40% shade soil_temp$dates <- micro$dates # Add in the dates# Re-arrange for plotting soil_temp2 <- soil_temp %>% select(D0cm, D50cm, dates) %>% pivot_longer(cols = c(D0cm, D50cm), names_to = "Depth", values_to = "Temperature")# Plot soil temps ggplot(soil_temp2) + geom_line(aes(y = Temperature, x = dates, col = Depth), show.legend = TRUE) + labs(y = "Soil Temperature (°C)", x = "Date") + theme_classic() + theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + scale_colour_viridis_d()```If we want to use air temperatures because our organism does not burrow then we can extract air temperatures at different heights above ground (we have simulated 20 cm and 2 m) (@fig-Air_temp). Of course, we could calculate temperatures for other heights as well, but assume this is the space most relevant to our system (@fig-case1), which we will detail below. We can also simulate microclimates under different scenarios in new calls to the `micro_silo` function but in our case we have managed to keep it simple and do everything in a single run. Information on air temperatures at our two heights is also stored in our `microclimate` object from the simulation we ran above. To obtain our air temperatures, these are stored in a different place in our `microclimate` object. Again, we have access to these calculations under two shade levels (40% and 80%), but we'll just use our lowest shade level for now. :::{.column-margin}:::{.callout-important}The microclimate object has a lot of relevant information that we can use in our simulations. Keep this in mind while exploring the [documentation](https://rdrr.io/github/mrke/NicheMapR/man/). We'll come back to this later on.::::::```{r, microclimate_air}#| label: microclimate_air#| echo: true#| warning: false#| message: false# Get the above ground temperatures at 20 cm and 2 m air_temp <- data.frame(micro$metout) air_temp$dates <- micro$dates```Now we can have a look at what the air temperature would be 20 cm and 2 m above ground (@fig-Air_temp). ```{r, fig-Air_temp}#| label: fig-Air_temp#| fig-cap: "Average air temperature at a location near Renmark, South Australia. Air temperatures at 20 cm and 2 m above ground (and in minimum shade, 40%) are shown for a short times between 2019-01-15 and 2019-01-31."#| code-fold: true air_temp2 <- air_temp %>% select(dates, TALOC, TAREF) %>% pivot_longer(cols = c(TALOC, TAREF), names_to = "Height", values_to = "Temperature") %>% mutate(Height = factor(ifelse(Height == "TALOC", "20 cm", "2 m"), levels = c("20 cm", "2 m")))# Plot soil temps air_temp2 %>% ggplot(aes(y = Temperature, x = dates, col = Height)) + geom_line() + labs(y = "Air Temperature (°C)", x = "Date") + theme_classic() + theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + scale_colour_viridis_d()```We can immediately see from @fig-Air_temp that air temperatures closer to the ground are hotter than those higher up because of increased radiative heat transfer from the ground to air in close proximity and convective cooling higher in the air column.Interestingly, between 2019-01-15 to 2019-01-31 we had a major heatwave event. Lets have a look at the temperatures during this period.```{r, fig-Air_temp_heatwave}#| label: fig-Air_temp_heatwave#| fig-cap: "Average air temperature at a location near Renmark, South Australia. Air temperatures at 20 cm and 2 m above ground (and in minimum shade, 40%) are shown for a short times between 2019-01-15 and 2019-01-31."#| code-fold: true air_temp2 %>% filter(dates > "2019-01-15" & dates < "2019-01-31") %>% ggplot(aes(y = Temperature, x = dates, col = Height)) + geom_line() + geom_hline(yintercept = 40, colour = "red", lty = 2) + labs(y = "Air Temperature (°C)", x = "Date") + theme_classic() + theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + scale_colour_viridis_d()```:::{.column-margin}:::{.callout-important}It is important to recognize that these are predictions. As with any model, they need field validation [@De_Frenne2024]. Nonetheless, they provide an important tool in our systems modelling framework to help us understand how heatwaves are likely to affect species interactions.:::::: Now that we've explored the various temperatures predicted in our microclimate model, we can use this information to predict the body temperatures of organisms. We can also use this information to predict how these organisms will respond to extreme heat events.# **Developing a systems modelling approach to predict biological responses to extreme heat**## **Case Study 1: Simulating plant and insect interactions under heatwaves** {#case1}### Setting the sceneCase study one explores how feedback between plants and insects can happen in response to extreme heat events. Here, we are interested in how the behaviour and body temperature of grasshoppers is influenced by changes in plant transpiration which can impact microhabitat temperatures of grasshoppers (@fig-case1). While this is a simple situation it's a useful starting point to understand how a systems modelling approach can be built up to predict the impact of heatwaves on species interactions. ```{r, fig-case1}#| label: fig-case1#| echo: false#| fig-cap: "Grasshopper and plant interactions under heatwaves. Grasshoppers feeding on leaves of a bush in Calperum Station, South Australia. Complex leaf structure provides shade and cooler microclimate for grasshoppers. Leaves exposed to hot temperatures close their stomata to conserve water reducing transpiration and photosynthesis. Leaves that are not exposed to direct sunlight can keep their stomata open and continue photosynthesis while also transpiring and providing cooler microhabitats that grasshoppers can exploit through thermoregulatory behaviours. Importantly, stomatal behaviour is mediated by drought conditions. More stomata will be closed during drought to conserve water and this will impact the cooling effect of the plant on the grasshopper."fig_case1 <- image_read(here("figs", "figS1.png"))fig_case1```### Predicting temperatures experienced by plants that form the grasshopper microclimate {#step1}To predict the impact of heatwaves on plant-insect interactions, we first need to predict how the microclimate that each organism experiences translates into body temperatures. Taking a systems approach means that we can account for the effects of other organisms on our focal species. To do this, we need to consider the biology of the species that co-exist. In this case, the plants that the grasshopper consumes will influence the microclimate the grasshopper experiences (@fig-case1). Plants will also respond to the thermal environmental, and this will directly and indirectly influence the grasshopper. For example, plants will change their transpiration through the day and this can impact local temperatures experienced by grasshoppers (@fig-case1). The microclimate model that we ran in [Section 1](#section1) can now be used to predict how the plant and grasshopper will respond to the heatwave. **Remember, we just used weather data to calculate air temperatures expected at 0.2 and 2m in both 40% and 80% shade. We can use this microclimate data to estimate leaf temperatures for a plant and then use the information on plant temperatures to feedback into the calculations of grasshopper body temperature.**Before we can estimate the temperatures for the grasshopper, we need to use a biophysical model to predict what the plant is doing. We can then use this 'new' microclimate data to simulate the growth of insects under different temperature conditions, changing what would happen if plants respond differently [@fig-case1].Let's use our microclimate simulations to predict the plant leaf and air temperatures at the relevant heights to grasshoppers.First, the traits of the leaf must be defined. An ellipsoid shape is used in this simulation. Specification of the mass and the ratios of the linear dimensions of the shape (here, semi-minor axes relative to the semi-major axis) together with the density (default value 1 g cm${^3}$) allows the geometric functions in the ectotherm model to compute the absolute dimensions and surface areas. Evaporative exchange is driven by the stomatal conductance parameters, which can be specified for the abaxial and adaxial sides of the leaf as well as the overall base leaf conductance through the cuticle. The `pct_wet` parameter is set to 0.Radiative exchange parameters can be set including the emissivities of the leaf and the substrate, the solar absorptivities of the leaf, the configuration factors to the sky and ground, the percentage of the leaf conducting to the ground and the orientation to the sun.```{r, microclimate_plant}#| label: microclimate_plant#| echo: true#| warning: false#| message: false # 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 save(plant, file = here("output", "microclim", "plant_micro.RDa"))``````{r, fig-plant_temp}#| label: fig-plant_temp#| fig-cap: "Leaf temperature of a plant at a location near Renmark, South Australia. Leaf temperature is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) between 2019-01-15 and 2019-01-31."#| code-fold: true # Load the plant object load(here("output", "microclim", "plant_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 environ_plant <- as.data.frame(plant$environ) environ_plant$dates <- micro$dates # For clarity 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-31") %>% 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)")) # Plot it out ggplot(environ_plant2) + geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend = TRUE) + labs(y = "Temperature (°C)", x = "Date") + theme_classic() + theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), legend.position = "top") + scale_colour_viridis_d()```We will now use the `plantecophys` package to calculate the leaf temperature and transpiration of the plant. This will allow us to understand how the plant is responding to the heatwave and how this might influence the grasshopper. The `plantecophys` package allows us to get more realistic stomatal behaviour via the `Photosyn` function. This code runs `Photosyn` assuming the leaf temperature calculated with the ectotherm model to start with.```{r}#| label: plant_stomatal#| echo: true#| warning: false#| message: false# Include simulated stomatal conductance from plantecophys UMOLPERJ <-4.57# mu_mol photons / JWETAIR_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 modelCa =400, # Atmospheric CO2 concentration (ppm)VPD = VPD_Pa/1000, # Vapour pressure deficit (kPa)vpdmin =0.5, # lower threshold on VPDPPFD = plant$metout[,"SOLR"] * UMOLPERJ *1, # mu mol m-2 s-1Patm = P_a/1000# Atmospheric pressure (kPa) )g_vs_abs_pho <- photo_out$GS```Now recalculate 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.```{r, plant_stomatal2}#| label: plant_stomatal2#| echo: true#| warning: false#| message: false# Rerun the models. Note that the paramaters as the same as 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 save(plant2, file = here("output", "microclim", "plant2_micro.RDa"))```Now, lets have a look at the leaf temperature when stomatal closure behaviour is incorporated.```{r, fig-plant_temp2}#| label: fig-plant_temp2#| fig-cap: "Leaf temperature of a plant at a location near Renmark, South Australia. Leaf temperature is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) between 2019-01-15 and 2019-01-31 for two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure."#| code-fold: true # 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_plantecophys2 <- as.data.frame(plant2$environ) environ_plant_plantecophys2$dates <- micro$dates # For clarirty add in dates from microclimate that ecto model uses # Re-organise the data for plotting environ_plant_plantecophys <- environ_plant_plantecophys2 %>% filter(dates > "2019-01-15" & dates < "2019-01-31") %>% 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) ggplot(new_environ_plant2) + geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend = TRUE) + labs(y = "Temperature (°C)", x = "Date") + theme_classic() + theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),legend.position = "top", legend.text = element_text(size = 6)) + scale_colour_viridis_d()```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. ### Coupling plant and animal models: Incorporating cooling effects of plant transpiration to estimate grasshopper body temperature {#step2}Now, lets run the ectotherm model for a grasshopper, but instead of using the microclimate object as it is, lets assume that the plant environment forms the microclimate of the grasshopper. This allows us to take in the microclimate and use it to calculate the body temperature of the grasshopper, and eventually understand how this affects growth, development etc within a dynamic energy budget model. The ectotherm model also allows us to incorporate relevant behavioural plasticity in response to extreme heat, which we'll do. ```{r, ecto_grasshopper}#| label: ecto_grasshopper#| echo: true#| warning: false#| message: false# First, lets fix up the microclimate object so that we can use it in the ectotherm model. We need to add in the leaf temperature as the microclimate for the grasshopper. There are a couple of ways we can conceptualise this. For example, we could assume that animals are in shade and therefore infrared radiative heat (`TSKYC`) is the temperature the grasshopper experiences from their environment. Alternatively, we could assume that when sitting on a leaf the grasshopper is in direct contact with the leaf and therefore the leaf temperature is the microclimate experienced by the grasshopper. We therefore need to treat the temperatures of the leaf as the environmental temperatures experienced by the grasshopper. Then, we feed this new microclimate into the grasshopper biophysical model to predict its body temperature. # Re-load the original micro object. We'll modify this now load(here("output", "microclim", "microclimate.RDa")) # Get some key mircoclimate details of the plants from the model that considers the behaviour of the stomata and update the mircoclimate. We'll do this for the 40% shade condition (metout). # First, assume that the air temp at the location is close to or equivilant to the temp of the plant leaf. Update miroclimate accordingly that will be used for grasshopper. micro$metout[,"TALOC"] <- environ_plant_plantecophys2[,"TC"] # Second, lets update the relative humidity. micro$metout[,"RHLOC"] <- environ_plant_plantecophys2[,"RELHUM"] # Third, lets also update the IR profiles. micro$metout[,"TSKYC"] <- environ_plant_plantecophys2[,"TSKY"] # Parameters for the grasshopper 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 (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 <- 12.5 # degrees C, minimum basking temperatureT_RB_min <- 12.5 # 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. This tells the model what the grasshopper can do in response to heat. diurn <- 1 # diurnal activity allowed (1) or not (0)? nocturn <- 0 # nocturnal activity allowed (1) or not (0)? crepus <- 1 # crepuscular activity allowed (1) or not (0)?shade_seek <- 1 # shade seeking allowed (1) or not (0)? burrow <- 0 # shelter in burrow allowed (1) or not (0)? pct_cond <- 0 # Percentage of animal surface contacting the substrate# Now that we have set the parameters we can run the modelgrasshopper <- 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)# Save the grasshopper objectsave(grasshopper, file = here("output", "microclim", "grasshopper_micro.RDa"))```Now that we have run the ectotherm model we can plot the body temperature of the grasshopper over time along side the air temperature. This will allow us to see how the grasshopper responds to the heatwave. We can also see activity patterns of the grasshopper which predict how we expect it to behave in relation to its body temperature and time of day. At night the animal is inactive while during suitable periods of the day where body temperatures are suitable it will forage (@fig-grasshopper_temp).```{r, fig-grasshopper_temp}#| label: fig-grasshopper_temp#| code-fold: true#| fig-cap: "Body temperature of a grasshopper at a location near Renmark, South Australia. Body temperature of the grasshopper is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) for a short times between 2019-01-15 and 2019-01-31. "# Load grasshopper object load(here("output", "microclim", "grasshopper_micro.RDa"))# Get the resulting heat budget calculations for the grasshopper grasshopper_environ <- as.data.frame(grasshopper$environ)grasshopper_environ$dates <- micro$dates# Filter to heat wave datesgrasshopper_environ2 <- grasshopper_environ %>% filter(dates > "2019-01-15" & dates < "2019-01-31") %>% select(dates, TA, TC, ACT) %>% pivot_longer( cols = c(TA, TC), names_to = "Temps", values_to = "Temperature") %>% mutate(Temps = ifelse(Temps == "TA", "Air Temperature", "Body Temperature")) ggplot(grasshopper_environ2) + geom_line(aes(y = Temperature, x = dates, col = Temps), show.legend = TRUE) + geom_point(aes(y = ACT, x = dates, fill = factor(ACT)), shape = 21, alpha = 0.3 ,size = 3) + labs(y = "Temperature (°C)", x = "Date") + theme_classic() + theme( axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + scale_colour_viridis_d() + scale_fill_manual(values = c("0" = "blue", "1" = "green", "2" = "orange"), labels = c("0" = "Inactive", "1" = "Basking", "2" = "Foraging")) + labs(fill = "Behaviour")```## **Case Study 2: Understanding the effects of leaf mirochabitat on insect life-cyles: Linking mircoclimates to dynamic energy budget and population ecology models** {#case2}### Setting the sceneOur next case study will explore how the microclimate produced by a plant leaf can influence the life-cycle of an insect. In this case, we will use a dynamic energy budget model to understand how the microclimate of a plant leaf can influence the growth and development of a grasshopper. Here, we will introduce a Dynamic Energy Budget (DEB) model which now takes as input the mircoclimate produced by the plant to predcit body temperatures, development, growth and reproduction of the grasshopper (see @fig-case2 for a simplified DEB model). The DEB model also allows us to incorperate the heat wave event.Dynamic Energy Budget (DEB) models are a powerful tool for understanding how organisms allocate energy to growth, reproduction and maintenance. To parameterise a DEB model it requires a set of parameters that are derived from empirical data on growth and metabolism. While DEB models are parameter rich, key parameters have been estimated for > 3000 species already [@Marques2018; @Kooijman2021]. New analytical approaches using phylogenetic imputation methods can be used to predict DEB parameters for over 1.3 million species providing immense opportunities to model how heat will affect life-cycles of multi-species assemblages (see https://deb.bolding-bruggeman.com, @Bruggeman2009). ```{r, fig-case2}#| label: fig-case2#| echo: false#| fig-cap: "A simple adult stage Dynamic Energy Budget (DEB) model for the grasshopper showing how food gets assimilated into reserve energy stores. The grasshopper than aloocates some fraction of the energy in reserves to maintaining body stucrture (somatic maintenance) and growth (i.e., structure) and to developing reproductive tissue (reproductive maintenance) and allocating to reproduction."fig_case2 <- image_read(here("figs", "FigS3.png"))fig_case2```### Grasshopper development and reproduction under heatwaves when plants don't respondTo setup our DEB model we first need to set the DEB model paarameters. Grasshoppers are slightly more complicated because of the differing life-stages but we can use an existing DEB model for the grasshopper *Warramaba virgo*, and use a model developed by @Llandres2015. Thankfully, @Llandres2015 have already estimated the DEB parameters for us so we'll use these. The DEB model parameters are available in the `sims` folder of the repository. We'll load these parameters into R and then set the initial conditions for the DEB model.```{r, deb_grasshopper_params}#| label: deb_grasshopper_params#| echo: true#| warning: false#| message: false# 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 reproductionphotofinish <- 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 (\%)```### Microclimate data for the DEB modelRecall from [Case Study 1](#case1) that we have already simulated the microclimate for the grasshopper assuming: 1) plant stoma remain open during the day or 2) close during the day. As we seen in @fig-plant_temp3, this had major implications for temperatures experienced in the plant microhabitat that grasshoppers inhabit. Under these two different scenerios, the key questions now become: 1. **How do the different microhabitats created by plant stomatal behaviour influence the development, growth and reproduction of the grasshopper?**2. **How do different vital rates impact population dynamics?**Answers to these questions are useful to see how a systems approach can be implemented. Lets first simulate the life-cyle of grasshoppers under the two different microclimate scenarios.```{r, deb_grasshopper}#| label: deb_grasshopper#| echo: true#| warning: false#| message: false# First, we need to load in the two microclimates scenerios. 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", "plant2_micro.RDa") plant_mirco_names <- c("stomata_open", "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"))# Now we can loop through the two microclimates and run the DEB model for each of themfor(i in 1:length(plant_micro)){# Load the microclimate object load(here("output", "microclim", plant_micro[i]))# Object is labelled plant2 for the plant2_micro.RDa object if(plant_micro[i] == "plant2_micro.RDa"){ plant <- plant2 } else { plant <- plant }# 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, hourlydates2 <- 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 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 = micro$minshade[subs2], maxshades = micro$maxshade[subs2], 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 run the DEB model we can have a look at whether the life-cycles of the grasshoppers are different under the two different microclimate scenarios.```{r, fig-grasshopper_lifecycle}#| label: fig-grasshopper_lifecycle#| fig-cap: "Life-cycle of a grasshopper at a location near Renmark, South Australia. Life-cycle of the grasshopper is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure. 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’).The grey lines indicate the body temperature of the grasshopper (°C) (right axis)."#| code-fold: true#| echo: true#| warning: false#| message: false# Load in the DEB model for scenerio 1 load(here("output", "microclim", "grasshopper_stomata_open_DEB.RDa")) 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# Create the life cycle plotp1 <- ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col = 'lightgrey', linewidth = 0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col = "forestgreen") +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 = "Predicted Grasshopper life-cycle under stomata open", x = "Date") + scale_y_continuous( name = "Wet mass (mg)", sec.axis = sec_axis(~ . / scale_factor, name = "Body temperature (°C)") )+ theme_classic()# Now for scenerio 2load(here("output", "microclim", "grasshopper_stomata_closed_DEB.RDa")) debout <- data.frame(ecto$debout) debout$datetime <- datetime# Create the life cycle plotp2 <- ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col = 'lightgrey', size = 0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col = "forestgreen") +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 = "Predicted Grasshopper life-cycle under stomata closed", x = "Date") + scale_y_continuous( name = "Wet mass (mg)", sec.axis = sec_axis(~ . / scale_factor, name = "Body temperature (°C)") )+ theme_classic()p1/p2 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A") + theme(plot.tag = element_text(size = 12, face = "bold"))```We can see from @fig-grasshopper_lifecycle that the life-cyle of the grasshopper is different when we consider the behaviour of the plant stomata. Notably, when stomata are closed and the leaf temperatures are higher we expect the grasshopper to move through stages faster and reproduce more (@tbl-grasshopper_lifecycle). We can have a look at the different traits expected under both scenarios. ```{r, tbl-grasshopper_lifecycle}#| label: tbl-grasshopper_lifecycle#| echo: true#| message: false#| warning: false#| code-fold: true#| tbl-cap: "Predicted life-cycle traits of a grasshopper at a location near Renmark, South Australia. Life-cycle traits of the grasshopper are shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure."# Traits we want to extract from the DEB model traits <- c("FECUNDITY", "CLUTCHES", "LifeSpan", "rmax") # Load in the DEB model for scenerio 1load(here("output", "microclim", "grasshopper_stomata_open_DEB.RDa")) life_cycle_open <- data.frame(ecto$yearout)[,traits]# Load in the DEB model for scenerio 2load(here("output", "microclim", "grasshopper_stomata_closed_DEB.RDa"))life_cycle_closed <- data.frame(ecto$yearout)[,traits]# Finalise the table.table_final <- rbind(life_cycle_open, life_cycle_closed)table_final <- table_final %>% mutate(Scenario = c("Stomata Open", "Stomata Closed")) %>% select(Scenario, everything())# Create a table of the life-cycle traitsflextable::flextable(table_final) %>% flextable::set_header_labels(values = c("FECUNDITY" = "Fecundity", "CLUTCHES" = "Clutches", "LifeSpan" = "Life Span", "Lamda Max" = "rmax")) %>% flextable::set_table_properties(layout = "autofit")```### Population dynamics of grasshoppers under two plant response sceneriosAs we have seen in the previous section, the life-cycle of the grasshopper is different when we consider the behaviour of the plant stomata and this has noticable effects on vital rates. We can use this information to parameterise a simple stochastic population models to understand how population dynamics of the grasshopper might change.We can use a simple Ricker model [@Ricker1954a], which we'll use to simulate the population dynamics of the grasshopper. The Ricker model is a simple model that describes how populations grow over time. The model is given by the equation:$$N_{t+1} = N_t \cdot e^{r(1 - \frac{N_t}{K})}$$ {#eq-ricker_model}where $N_t$ is the population size at time $t$, $r$ is the intrinsic growth rate, and $K$ is the carrying capacity of the environment. The model assumes that populations grow exponentially when they are small, but that growth slows down as the population approaches the carrying capacity.We can also re-parameterise to use a finite rate model given that the grasshoppers are an annual species.$$N_{t+1} = N_t \cdot \lambda^{(1 - \frac{N_t}{K})}$$ {#eq-ricker_model2}where $\lambda$ is the finite rate of increase. $$\lambda = e^{r_{max}}$$ {#eq-ricker_lamda}While @eq-ricker_model2 is a useful model it does make some simplifying assumptions. First, it assumes discrete non-overlapping generations. This should be a sensible assumption for insect and probably many plant species. Second, $\lambda$ and $K$ are assumed constant over time, ignoring any effects of environmental stochasticity or evolution. However, this assumption can easily be relaxed. Third, it assumes a closed system with no emmigration or immigration. However, again, we can modify the model to include these factors. For now, we will assume that the population is closed and that $\lambda$ and $K$ are constant over time.Recall from our previous section that we already estimated the intrinsic growth rate. Using this information, for the stomata-closed scenario ($\lambda$ = `r round(exp(1.0360276), 2)`), the population could tolerate up to `r round((1-(1/exp(1.0360276)))*100, 2)`% mortality per generation and still remain stable.In contrast, for the stomata-open scenario ($\lambda$ = `r round(exp(0.4188939),2)`), the population is increasing more slowly. The population can tolerate around `r abs(round((1-(1/exp(0.4188939)))*100, 2))`% mortality per generation and still remain stable.Assuming we had information on the carying capacity of the enviroment and population size we can then use the $r_{max}$ in the Ricker model to simulate the population dynamics of the grasshopper. ```{r, ricker_model}#| label: ricker_model#| echo: true#| message: false#| warning: false#' @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 sizericker <- 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 nrepssim_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)}``````{r, ricker_simulation}#| label: ricker_simulation#| echo: true#| message: false#| warning: falseset.seed(123) # Set seed for reproducibility# Set up the parameters for the Ricker model N0 <- 1000 # Initial population sizer_closed <- 1.0360276 # Intrinsic growth (DEB model) with stoma closed r_open <- 0.4188939 # Intrinsic growth (DEB model) with stoma open sd_r <- 0.15 # Standard deviation of the growth rate K <- 10000 # Carrying capacity time <- 10 # Time steps nsims <- 1000 # Number of simulations# Lets run the simulationsN_closed <- sim_ricker(N0, r_closed, sd_r, K, time, nsims) N_open <- sim_ricker(N0, r_open, sd_r, K, time, nsims)# Create a data frame to store the resultssims_closed <- pivot_longer(as.data.frame(N_closed), cols = everything(), names_to = "sim", values_to = "N") %>% mutate(time = rep(1:time, each = nsims))sims_open <- pivot_longer(as.data.frame(N_open), cols = everything(), names_to = "sim", values_to = "N") %>% mutate(time = rep(1:time, each = nsims))``````{r, fig-ricker_simulation}#| label: fig-ricker_simulation#| fig-cap: "Population dynamics of a grasshopper at a location near Renmark, South Australia. Population dynamics of the grasshopper is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure."#| code-fold: true#| echo: true open_stoma <- image_read(here("figs", "Stomata_open.png")) open <- rasterGrob(open_stoma)closed_stoma <- image_read(here("figs", "Stomata_closed.png")) closed <- rasterGrob(closed_stoma)# Plot the resultsp1 = ggplot(sims_closed, aes(x = time, y = N, colour = sim)) + geom_line() + labs(title = "Stomata Closure", x = "Generation (time)", y = "Population size (N)") + theme_classic() + theme(legend.position = "none") + scale_colour_viridis_d() + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 14)) + annotation_custom(closed, xmin = 5, xmax = 9, ymin = 500, ymax = 8500) + ylim(0, 12000)p2 = ggplot(sims_open, aes(x = time, y = N, colour = sim)) + geom_line() + labs(title = "Stomata Remain Open", x = "Generation (time)", y = "Population size (N)") + theme_classic() + theme(legend.position = "none") + scale_colour_viridis_d() + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 14)) + annotation_custom(open, xmin = 1, xmax = 5, ymin = 7000, ymax = 11000) + ylim(0, 12000)p1 + theme(plot.tag = element_text(size = 12, face = "bold")) + p2 + theme(plot.tag = element_text(size = 12, face = "bold")) + plot_layout(ncol = 2) + plot_annotation(tag_levels = "A") ```We can see from @fig-ricker_simulation that when stomata are able to respond by closing grasshoppers are expected to grow rapidly until they reach carrying capacity (K = 10,000, and the point they all converge on), whereas if we have stomata remaining open we expect the population to grow much more slowly over time as microhabiats become less suitable for growth and reproduction. This happens because higher body temperatures that result from when stomata close better allows grasshoppers to reach optimal physiological performance.## **Case Study 3: Heat wave impacts through small communities** {#case3}### Setting the sceneA systems modelling approach naturally involves attempting to model multi-species assemblies / communities. This is challenging with large communities because of the complexity of the interactions between species. However, for simpler communities such an approach of connecting models across fields is achievable and potentially powerful in making predictions about how communities will respond to extreme heat events.In this case study, we expand on the impacts of extreme heat on a single species of grasshopper and show how this can be extended to capture interactions between three species: a plant, a grasshopper and a second competing grasshopper species that occupies the same mircohabitat. We extend our Ricker model to capture the impacts that each grasshopper species has on the other. One way to model this is to explicity translate how each species' interactions impact the population size of the other.### Some TheoryRecall our original Ricker model (@eq-ricker_model2):$$N_{t+1} = N_t \cdot \lambda^{\left(1 - \frac{N_t}{K}\right)}$$In keeping with our notation, we can extend this model to a two-species Ricker model to include the interactions between the two grasshopper species as follows:$$N_{t+1}^A = N_t^A \cdot \lambda_{A}^{\left(1 - \frac{N_t^A + \alpha N_t^B}{K_A}\right)} $$ {#eq-A}$$N_{t+1}^B = N_t^B \cdot \lambda_{B}^{\left(1 - \frac{N_t^B + \beta N_t^A}{K_B}\right)} $$ {#eq-B}where $N_t^A$ and $N_t^B$ are the population sizes of species A and B at time $t$, $\lambda_A$ and $\lambda_B$ are the intrinsic growth rates of species A and B, $K_A$ and $K_B$ are the carrying capacities of species A and B, and $\alpha$ is the competition coefficient that describes how much species B affects species A, whereas $\beta$ is the competition coefficient that describes how much species A affects species B.In @eq-A, the competition coefficient $\alpha$ is a measure of the strength of the interaction between the two species. If $\alpha$ is greater than 1, it means that species B has a stronger negative effect on species A than species A has on itself. If $\alpha$ is less than 1, it means that species A has a stronger negative effect on itself than species B has on it. If $\alpha$ is equal to 1, it means that the two species have equal effects on each other. The opposite is true for $\beta$ in @eq-B.A more general form of this model is given by:$$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)}$$ {#eq-general}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.Just like @eq-ricker_model2, @eq-general is also making similar and a few additional assumptions that wre worth noting. 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. 5. while we assume symmetric interactions, it allows for asymmetric interaction (i.e., $\alpha_{ij} \neq \alpha_{ji}$) For now, we will assume that the population is closed and that $K$ is constant over time. While this model allows for asymmetric interaction (i.e., $\alpha_{ij} \neq \alpha_{ji}$), we will assume that the interactions are symmetric (i.e., $\alpha_{ij} = \alpha_{ji}$). 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.### Simulating effects of extreme heat on two speciesWe can simulate the life-cyle of a second grasshopper species that occupies the same habitat but varies in physiology. We've already done this for one species *Warramaba virgo* (Species B; @fig-two_species_grasshopper) in [Case Study 2](#case2). ```{r, fig-two_species_grasshopper}#| label: fig-two_species_grasshopper#| echo: true#| message: false#| warning: false#| fig-cap: "Simulating two species competitive interactions."fig_casetwospecies <- image_read(here("figs", "FigS4.png"))fig_casetwospecies```The second species (Species A; @fig-two_species_grasshopper) has a different physiology compared to the first species but is usually more abundant than Species B making it highly competitive. We'll now simulate its life-cycle like we did for Species B in [Case Study 2](#case2). Then, link this to a two-sepcies population dynamic model. These species are closely related but do have some differences in physiology that include: 1) higher thermal tolerance, 2) larger clutch sizes, and 3) different preferred body temperatures; and 4) smaller size. Otherwise, we will assume the same DEB model parameters as Species B.```{r, speciesA_grasshopper}#| label: speciesA_grasshopper#| echo: true#| message: false#| warning: false# 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 <- L.b # Structural length at birth (cm)}\cr} photostart <- 0 # no photoperiod effect on reproductionphotofinish <- 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 <- 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 <- 1 # Shade seeking allowed (1) or not (0)? pct_cond <- 0 # Percentage of animal surface contacting the substrate (\%)``````{r, deb_grasshopperA}#| label: deb_grasshopperA#| echo: true#| warning: false#| message: false# First, we need to load in the two microclimates scenerios. 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", "plant2_micro.RDa") plant_mirco_names <- c("stomata_open", "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"))# Now we can loop through the two microclimates and run the DEB model for each of themfor(i in 1:length(plant_micro)){# Load the microclimate object load(here("output", "microclim", plant_micro[i]))# Object is labelled plant2 for the plant2_micro.RDa object if(plant_micro[i] == "plant2_micro.RDa"){ plant <- plant2 } else { plant <- plant }# 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, hourlydates2 <- 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 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 = micro$minshade[subs2], maxshades = micro$maxshade[subs2], 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 = 0.0779, # Gompertz stress coefficient, from Debber 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 = 0.09, # structural length at birth (cm), photostart = photostart, photofinish = photofinish, reset = reset, soilnode = soilnode)# Save the DEB model save(ecto, file = here("output", "microclim", paste0("grasshopperspA_", plant_mirco_names[i], "_DEB.RDa")))}```We can have a look at the life-cycle of the new species of grasshopper, Species A, under both scenarios (@fig-grasshopperspA_lifecycle).```{r, fig-grasshopperspA_lifecycle}#| label: fig-grasshopperspA_lifecycle#| fig-cap: "Life-cycle of a grasshopper (Species A) at a location near Renmark, South Australia. Life-cycle of the grasshopper is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure. 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’).The grey lines indicate the body temperature of the grasshopper (°C) (right axis)."#| code-fold: true#| echo: true#| warning: false#| message: false# Load in the DEB model for scenerio 1 load(here("output", "microclim", "grasshopperspA_stomata_open_DEB.RDa")) 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# Create the life cycle plotp1 <- ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col = 'lightgrey', linewidth = 0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col = "forestgreen") +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 = "Predicted Species A Grasshopper life-cycle under stomata open", x = "Date") + scale_y_continuous( name = "Wet mass (mg)", sec.axis = sec_axis(~ . / scale_factor, name = "Body temperature (°C)") )+ theme_classic()# Now for scenerio 2load(here("output", "microclim", "grasshopperspA_stomata_closed_DEB.RDa")) debout <- data.frame(ecto$debout) debout$datetime <- datetime# Create the life cycle plotp2 <- ggplot() +geom_line(data = ecto$environ, aes(x = datetime, y = TC*scale_factor), col = 'lightgrey', size = 0.5) +ylim(-10, 500) +geom_line(data = ecto$environ, aes(x = datetime, y = SHADE), col = "forestgreen") +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 = "Predicted Species A Grasshopper life-cycle under stomata closed", x = "Date") + scale_y_continuous( name = "Wet mass (mg)", sec.axis = sec_axis(~ . / scale_factor, name = "Body temperature (°C)") )+ theme_classic()p1/p2 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A") + theme(plot.tag = element_text(size = 12, face = "bold"))```We can now put the summaries of all our life-cycle predictions together for both species of grasshopper and both scenerios (@tbl-grasshopperssp_lifecycle). We'll use the intrinsic growth rates in our population models.```{r, tbl-grasshopperssp_lifecycle}#| label: tbl-grasshopperssp_lifecycle#| echo: false#| message: false#| warning: false#| tbl-cap: "Life-cycle of two grasshopper species at a location near Renmark, South Australia. Life-cycle of the grasshoppers is shown in under two scenarios, one where stomata are open all day, and one where stomatal closure is driven by photosynthetically active radiation intensity and vapour pressure."# Now lets extract the life-cycle parameters from the DEB model and put them into a table debs1 <- c(paste0("grasshopperspA_", plant_mirco_names, "_DEB.RDa"), paste0("grasshopper_", plant_mirco_names, "_DEB.RDa")) scenerios <- c("Stomata Open (Species A)", "Stomata Closed (Species A)", "Stomata Open (Species B)", "Stomata Closed (Species B)")# Create a data.frame to store the life-cycle parameters life_cycle <- data.frame()# Loop through the DEB models and extract the life-cycle parametersfor(i in 1:length(debs1)){# Load the DEB model load(here("output", "microclim", debs1[i])) yearout <- data.frame(ecto$yearout) life_cycle <- rbind(life_cycle, yearout[1,])}# Lets add a column for the different scenerios and re-arrange the data so it's the first column life_cycle <- life_cycle %>% mutate(scenario = scenerios) %>% select(scenario, everything())# Remove irrelevant columns life_cycle <- life_cycle %>% select(-c("ANNUALACT", "Length")) colnames(life_cycle) <- c("Watering Scenario", "Development time (days)", "Birth day (day of year)", "Mass at birth (g)", "Months to maturity", "Months to reproducton", "Length (mm) at first reproduction", "Total fecundity (# eggs)", "Total clutches", "Minimum reserve density (J/cm3)", "Food eaten in last year (kg)", "Total food eaten across life (kg)", "Minimum body temperature (°C)", "Maximum body temperature (°C)", "Maximum level of desiccation", "Maximum lifespan (days)", "Generation time (years)", "Net reproductive rate (R0)", "Intrinsic growth rate (r)" )life_cycle %>% mutate(across(where(is.numeric), sprintf, fmt = '%.2f')) %>%datatable(., extensions = "Buttons", rownames = FALSE)```### Two-species population dynamics Modelling the dynamics between two species requires some knowledge about how competition impacts the population size of each species. Normally, this is achieved by estimating competition coefficients from field or laboratory data. Experiments involve determining how the presence of one species impacts growth of the other relative to a control state where only the single species is present. Given that our goal is simply to demonstrate how models from different fields can be coupled, we'll assume that the competition coefficients are known and equal so that $\alpha_{i,j}$ = 0.2. Lets re-write our functions to use the two-species model now. ```{r, ricker_two_species_model}#| label: ricker_two_species_model#| echo: true#| message: false#| warning: false#' @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_ij Competition coefficient for species j on species i#' @param N0_j Initial population size of species j#' @return Next years population sizericker_two <- function(N0_i, r_i, sd_r_i, K_i, N0_j, alpha_ij){ N0_i * exp(rnorm(1, r_i, sd_r_i)*(1 - ((N0_i + (alpha_ij * 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 time Time step#' @param nsims Number of replicates#' @return A matrix of population sizes at each time step for nsimssim_ricker_two_sp <- function(N0_i, r_i, sd_r_i, K_i, N0_j, r_j, sd_r_j, K_j, alpha_ij, 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_ij) # 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))}```We can now use these functions to simulate the population dynamics for both species under the two different scenarios. For simplity we have assumed the competition coefficient is the same for both species and does not vary with environment. Under more realistic conditions it's very likely that the competition coefficient would vary with environment and we could use our microclimate model to estimate this.```{r, ricker_simulation_twosp}#| label: ricker_simulation_twosp#| echo: true#| message: false#| warning: falseset.seed(123) # Set seed for reproducibility# Set up the parameters for the Ricker model N0_i <- 1000 # Initial population size, species A N0_j <- 500 # Initial population size, species B sd_r_i <- 0.15 # Standard deviation of the growth rate species A sd_r_j <- 0.15 # Standard deviation of the growth rate species B K_i <- 10000 # Carrying capacity species A K_j <- 10000 # Carrying capacity species Br_i_j_open <- life_cycle[c(1,3),"Intrinsic growth rate (r)"] # Intrinsic growth rate species A and B open watering scenerio r_i_j_closed <- life_cycle[c(2,4),"Intrinsic growth rate (r)"] # Intrinsic growth rate species A and B closed watering scenerio alpha_ij <- 0.8 # Competition coefficient. Assumed not to vary with environment time <- 10 # Time steps nsims <- 1000 # Number of simulations# Lets run the simulationsN_closed <- sim_ricker_two_sp(N0_i, r_i = r_i_j_closed[1], # Species i growth rate sd_r_i, K_i, N0_j, r_j = r_i_j_closed[2], # Species j growth rate sd_r_j, K_j, alpha_ij, time, nsims)N_open <- sim_ricker_two_sp(N0_i, r_i = r_i_j_open[1], # Species i growth rate sd_r_i, K_i, N0_j, r_j = r_i_j_open[2], # Species j growth rate sd_r_j, K_j, alpha_ij, time, nsims)# Create a data frame to store the results for each species for each sceneriosims_closed <- plyr::ldply(lapply(N_closed, function(x){ pivot_longer(as.data.frame(x), cols = everything(), names_to = "sim", values_to = "N") %>% mutate(time = rep(1:time, each = nsims))}))sims_open <- plyr::ldply(lapply(N_open, function(x){ pivot_longer(as.data.frame(x), cols = everything(), names_to = "sim", values_to = "N") %>% mutate(time = rep(1:time, each = nsims))}))```Now we can plot the results of the simulations for both species under both scenarios (@fig-ricker_twosp).```{r, fig-ricker_twosp}#| label: fig-ricker_twosp#| echo: true#| code-fold: true#| message: false#| warning: false#| fig-cap: "Simulated population dynamics for two species of grasshopper under two different watering scenarios. A) Stomata open, B) Stomata closed. The lines show the mean population size across 1000 simulations and the shaded area shows the standard deviation." open_stoma <- image_read(here("figs", "Stomata_open.png")) open <- rasterGrob(open_stoma)closed_stoma <- image_read(here("figs", "Stomata_closed.png")) closed <- rasterGrob(closed_stoma) # Create a summary of the simulations to make plotting easier sims_closed_sum <- sims_closed %>% group_by(.id, time) %>% summarise(N_m = mean(N), sd = sd(N)) sims_open_sum <- sims_open %>% group_by(.id, time) %>% summarise(N_m = mean(N), sd = sd(N)) # Plot the resultsp1 = ggplot(sims_closed_sum, 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(title = "Stomata Closure", x = "Generation (time)", y = "Population size (N)") + theme_classic() + theme(legend.position = "right") + scale_colour_viridis_d() + scale_fill_viridis_d() + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 14), legend.position = "none") + annotation_custom(closed, xmin = 1, xmax = 5, ymin = 7500, ymax = 12000) + ylim(0, 12000) + scale_colour_viridis_d( labels = c("Species A", "Species B"), name = "Species") +scale_fill_viridis_d( labels = c("Species A", "Species B"), name = "Species")p2 = ggplot(sims_open_sum, 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(title = "Stomata Remain Open", x = "Generation (time)", y = "Population size (N)") + theme_classic() + theme(legend.position = "right") + scale_colour_viridis_d() + scale_fill_viridis_d() + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 14)) + annotation_custom(open, xmin = 1, xmax = 5, ymin = 7500, ymax = 12000) + ylim(0, 12000) + scale_colour_viridis_d( labels = c("Species A", "Species B"), name = "Species") +scale_fill_viridis_d( labels = c("Species A", "Species B"), name = "Species")p1 + theme(plot.tag = element_text(size = 12, face = "bold")) + p2 + theme(plot.tag = element_text(size = 12, face = "bold")) + plot_layout(ncol = 2) + plot_annotation(tag_levels = "A") ```More complex models that incorporate a multitude of species and their interactions can be built using the same principles [e.g., @Bimler2018; @Bimler2023; @Bimler2024]. Information on interaction strength between species can be used to understand how species coexistence is expected to change under extreme heat scenarios.# **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.## **Case Study 4: Simulating the effect of watering on plant and insect body temperatures** {#case4}### Setting the sceneLet'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-case4). The challenge with watering is that often we are drought stressed in combination with the extreme heat [@Dai2013]. 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?**```{r, fig-case4}#| label: fig-case4#| echo: false#| fig-cap: "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."fig_case2 <- image_read(here("figs", "figS2.png"))fig_case2```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-water). The rainfall data is stored in the `microclimate` object and is provided as a daily value (rather than hourly as for temperature). ```{r, fig-water}#| label: fig-water#| fig-cap: "Rainfall events over time for A) over years and B) only during the heat wave event from 2019-01-15 to 2019-01-31."#| code-fold: true#| echo: true# 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 plotsp1 + p2 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A") + theme(plot.tag = element_text(size = 12, face = "bold"))```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 (@fig-plant_temp2). 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. ```{r, watering}#| label: watering#| echo: true#| warning: false#| message: false# 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.```{r, fig-water2}#| label: fig-water2#| fig-cap: "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)."#| code-fold: true#| echo: true # Add in rainfall manipulations to columnsrainfall_actual$rainfall_2 <- RAINFALL2rainfall_actual$rainfall_5 <- RAINFALL3rainfall_actual$rainfall_10 <- RAINFALL4rainfall_actual$rainfall_20 <- RAINFALL5# Re-arrange data for plottingrainfall_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``````{r}#| label: fig-water2.2#| echo: false#| warning: false#| message: false#| include: falseggsave(water, filename =here("figs", "fig-water2.png"), width =7.419753, height =5.938272, dpi =600)```### 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.```{r, water_micoclimate}#| label: water_micoclimate#| echo: true#| eval: false#| warning: false#| message: false# 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 = 10for(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 objectsave(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.### Calculating the plant leaf temperature using the watering scenerios```{r, water_plant}#| label: water_plant#| echo: true#| eval: false#| warning: false#| message: false# 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 modelfor(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 / JWETAIR_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")))}```### Impacts of watering scenerios on leaf temperatureWe'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.```{r}#| label: fig-plant_temp3#| fig-cap: "Leaf temperature of a plant at a location near Renmark, South Australia. Leaf temperature is shown in relation to air temperature at 2 m above ground (and in minimum shade, 40%) between 2019-01-15 and 2019-01-31 for three different watering scenarios (original rainfall, 2, 5, 10 and 20 mm of water added)."#| code-fold: true# 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 objectobjects <-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 in1:length(objects)){# Load the plant objectload(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 scenariosplant_temps_water <- environ_plant2 %>%filter(dates >"2019-01-15"& dates <"2019-01-31") %>%ggplot() +geom_line(aes(y = TC, x = dates, col = .id), show.legend =TRUE) +labs(y ="Leaf Temperature (°C)", x ="Date") +theme_classic() +theme( axis.text.x =element_text(angle =45, vjust =1, hjust =1), legend.position ="top") +scale_colour_viridis_d() +theme(legend.title =element_blank(),legend.text =element_text(size =12),axis.title =element_text(size =18),axis.text =element_text(size =12)) plant_temps_water``````{r}#| label: fig-plant_temp3.1#| echo: false#| warning: false#| message: false#| include: falsepw1<- plant_temps_water +theme(axis.text.x =element_blank(),legend.text =element_text(size =16),axis.title =element_text(size =24)) +annotate("text", x =as.POSIXct("2019-01-15"), y =52, label =TeX("$^\\circ$C Decrease"), size =6, hjust =0)+annotate("text", x =as.POSIXct("2019-01-19"), y =52, label =TeX("-2.80$^\\circ$C"), size =6, hjust =0) +annotate("text", x =as.POSIXct("2019-01-22"), y =52, label =TeX("-3.63$^\\circ$C"), size =6, hjust =0) +annotate("text", x =as.POSIXct("2019-01-25"), y =52, label =TeX("-3.82$^\\circ$C"), size =6, hjust =0) +annotate("text", x =as.POSIXct("2019-01-28"), y =52, label =TeX("-3.87$^\\circ$C"), size =6, hjust =0)ggsave(pw1, filename =here("figs", "fig-plant_temp3.png"), width =12.432099, height =6.876543, dpi =600)``````{r, tbl-plant_temp3}#| label: tbl-plant_temp3#| tbl-cap: "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."#| echo: false# Have a look at the temperature reductions tbl_temps <- environ_plant2 %>% filter(dates > "2019-01-15" & dates < "2019-01-31" & TIME > 11 & TIME < 16) %>% mutate(.id = ordered(.id, levels = c("Original Rainfall", "Watering (2 mm)", "Watering (5 mm)", "Watering (10 mm)", "Watering (20 mm)"))) %>% group_by(.id) %>% summarise(mean_temp = mean(TC, na.rm = TRUE)) %>% mutate(temp_diff = round(mean_temp - mean_temp[.id == "Original Rainfall"],2)) %>% select(.id, temp_diff) %>% data.frame() %>% rename(c("Watering Scenario" = ".id", "Mean Leaf Temperature Difference (°C)" = "temp_diff")) flextable(tbl_temps) %>% flextable::align(align = "center") %>% flextable::bold(part = "header") %>% flextable::align(align="center") %>% flextable::autofit()```@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 `r abs(tbl_temps[2,2])` $^\circ\text{C}$ (@tbl-plant_temp3). 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!## **Case Study 5: Impacts of watering scenerios on grasshopper life cycle** {#case5}### Setting the sceneNow that we have the leaf temperature of the plant under the various watering scenerios, we can use this to understand how the different scenerios might affect grasshopper body temperature, life cycle and population dynamics like we did in [case study 2](#case2). This is a nice demonstrations 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 scenerios 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.### Use the different watering scenerios 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 tempeartures the grasshopper experiences using a heat-budget model to predict the grasshopper body temperature. Then, we will use the heat-budget calculations to run teh 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.```{r, deb_grasshopper_water}#| label: deb_grasshopper_water#| echo: true#| warning: false#| message: false# 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 themfor(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, hourlydates2 <- 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 = micro$minshade[subs2], maxshades = micro$maxshade[subs2], 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-grasshopper_lifecycle_water). ```{r, fig-grasshopper_lifecycle_water}#| label: fig-grasshopper_lifecycle_water#| fig-cap: "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')."#| code-fold: true#| echo: true#| warning: false#| message: false# 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 listplot_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"))```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.```{r, tbl-grasshopper_lifecycle_watering}#| label: tbl-grasshopper_lifecycle_watering#| tbl-cap: "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)."#| echo: false#| code-fold: true#| warning: false#| message: false# Now lets extract the life-cycle parameters from the DEB model and put them into a table debs <- paste0("grasshopper_", plant_mirco_names, "_DEB.RDa") scenerios <- c("Original Rainfall", "Watering (2 mm)", "Watering (5 mm)", "Watering (10 mm)", "Watering (20 mm)")# Create a data.frame to store the life-cycle parameters life_cycle <- data.frame()# Loop through the DEB models and extract the life-cycle parametersfor(i in 1:length(debs)){# Load the DEB model load(here("output", "microclim", debs[i])) yearout <- data.frame(ecto$yearout) life_cycle <- rbind(life_cycle, yearout[1,])}# Lets add a column for the different scenerios and re-arrange the data so it's the first column life_cycle <- life_cycle %>% mutate(scenario = scenerios) %>% select(scenario, everything())# Remove irrelevant columns life_cycle <- life_cycle %>% select(-c("ANNUALACT", "Length")) colnames(life_cycle) <- c("Watering Scenario", "Development time (days)", "Birth day (day of year)", "Mass at birth (g)", "Months to maturity", "Months to reproducton", "Length (mm) at first reproduction", "Total fecundity (# eggs)", "Total clutches", "Minimum reserve density (J/cm3)", "Food eaten in last year (kg)", "Total food eaten across life (kg)", "Minimum body temperature (°C)", "Maximum body temperature (°C)", "Maximum level of desiccation", "Maximum lifespan (days)", "Generation time (years)", "Net reproductive rate (R0)", "Intrinsic growth rate (r)" )life_cycle %>% mutate(across(where(is.numeric), sprintf, fmt = '%.2f')) %>%datatable(., extensions = "Buttons", rownames = FALSE)```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](#case2). ```{r}#| label: fig-grasshopper_growth_rate_water#| fig-cap: "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."#| echo: true#| warning: false#| message: false#| code-fold: trueset.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 <-10# Time steps nsims <-1000# Number of simulations# Lets run the simulationssims <-list()for(i in1:length(r)){ sims[[i]] <-sim_ricker(N0, r[i], sd_r, K, time, nsims)names(sims)[i] <- names[i]}# Re-organise the data for plottingdata_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 sceneriodata_organised2 <- plyr::ldply(data_organised, .id ="scenario") %>%mutate(scenerio_sim =interaction(scenario, sim)) plot_water_sims <-ggplot(data_organised2, aes(x = time, y = N, col = scenario)) +geom_line(aes(group = scenerio_sim)) +labs(x ="Time (Generations)", y ="Population size (N)") +theme_classic() +scale_colour_viridis_d() +theme(axis.title =element_text(size =20),axis.text =element_text(size =20),legend.position ="inside",legend.position.inside =c(0.7,0.5),legend.text =element_text(size =20),legend.title =element_blank())plot_water_sims``````{r}#| label: fig-grasshopper_pop_growth_out#| echo: false#| warning: false#| message: false#| include: falseggsave(plot_water_sims, filename =here("figs", "fig-grasshopper_growth_rate_water.png"), width =7, height =7, dpi =600)```# **References**
