In ggcostoya/limon: Evolutionary Constraints among Thermal Performance Traits

knitr::opts_chunk$set(echo = TRUE)
devtools::load_all()

The goal of this vignette is to illustrate how survival and reproduction are determined in any simulation run in this project. The explanation will drawn from an example case.

Data

The basic unit we will be working with is an individual's thermal-performance data (TPD). As shown in the Basics vignettes, an individual's TPD is a tibble of temperature (t) and performance (p) values. We will first generate a basic TPD using the function gen_tpd (for further details go to the Generate TPD vignette). From this data, we will also draw a thermal-performance curve (TPC) using the gen_tpc function.

ex_tpts <- tibble(tpt = c("topt", "tb", "skw", "ctmin", "ctmax", "pmax", "pmin"),
               value = c(30, 3, -1, 20, 35, 10, 0.1))
tpd <- gen_tpd(tpts = ex_tpts, samples = 10, error = 0.75)
fit <- fit_tpd(tpd = tpd)
tpc <- gen_tpc(fit = fit)
plot(tpd, pch = 19, cex = 1.25, xlab = "Temperature", ylab = "Performance", ylim = c(0,max(tpd$p,tpc$p, na.rm = T)))
lines(tpc, type = "l", col = "grey", lwd = 2)
lines(tpd, type = "p", pch = 19, cex = 1.25)

Principle

In the context of this project, the principle to determine an individual's survival and reproduction based on its relative performance at a certain environmental temperature. As an example, let's imagine that the environmental temperature on a particular moment in time is of et = 27. From the individual's TPC we can determine the predicted p at the same temperature, in this case:

p_et <- tpc$p[tpc$t == Closest(tpc$t, 27, na.rm =T)]
p_et

Once we have determined p at et, we then determine the individual's maximal performance capacity (pmax), which is p at optimal temperature (topt). We can do so using the function get_tpts. In this case is:

tpts <- get_tpts(tpc, pmin = 0.1) # get tpts from the tpc
pmax <- tpts %>% filter(tpt == "pmax") %>% select(value) %>% as.numeric() # get pmax
pmax

topt <- tpts %>% filter(tpt == "topt") %>% select(value) %>% as.numeric()
ps_data <- tibble(t = c(27, topt), p = c(p_et,pmax), type = c("P at Env. T", "Pmax"))

plot(tpd, pch = 19, cex = 1, col = "grey", xlab = "Temperature", ylab = "Performance", ylim = c(0,max(tpd$p,tpc$p, na.rm = T)))
lines(tpc, type = "l", col = "grey", lwd = 2)
segments(x0 = 27, y0 = -1, x1 = 27, y1 = p_et, lwd = 2, lty = 2, col = "red")
segments(x0 = 27, y0 = p_et, x1 = 0, y1 = p_et, lwd = 2, lty = 2, col = "red")
segments(x0 = topt, y0 = -1, x1 = topt, y1 = pmax, lwd = 2, lty = 2)
segments(x0 = topt, y0 = pmax, x1 = 0, y1 = pmax, lwd = 2, lty = 2)
lines(x = ps_data$t, y = ps_data$p, col = c("red","black"), type = "p", pch = 19, cex = 1.25)
legend(x = max(tpd$t) - 3, y = max(tpd$p, tpc$p, ps_data$p, na.rm = T), unique(ps_data$type),col=c("red", "black"),pch=19, bty = "n", cex = 0.9)

With these two values we can determine the individual's relative performance (p_r) as p_et/pmax.

p_r <- p_et/pmax
p_r

p_r represents at what percentage of pmax at which the individual is performing. The relative performance will be later used to calculate the individual's probability of survival and reproduction, as is explained in the following sections.

Simulate Survival

To simulate the survival or death of a particular individual we use the sim_surv function included in this project, which uses the individual's p_r as the the probability of a binomial distribution of size 1. This approach guarantees a binomial response (0 for death, 1 for survival) being survival as likely as p_r. As an example to illustrate this approach, below we plot the result of drawing 1000 samples from a size 1 and probability p_r binomial distribution. Percentually, the count of samples with value 1 should be very close to p_r.

samples <- rbinom(1000,1,p_r)
binomial <- tibble(count = c(sum(samples), 1000 - sum(samples)), value = c("1", "0"))
ggplot(binomial, aes(x = value, y = count)) + geom_bar(stat = "identity", width = 0.5) + geom_text(aes(label = count), vjust = -0.5)+ theme_classic() + xlab("Value") + ylab("Count")

As an optional feauture, our simulations contemplate the possibility of adding a survival probability scaling parameter (spsp). A constant we can use to multiply with p_r in order to adjust relative performance to more realistic values. Considering the previous example in realistic terms, the survival probability for a particular moment in time would be:

p_r

However high, if the moment in time considered is in day units that is actually an extremely low daily survival probability. If we set spsp to 1.25 that same probability is boosted to more realistic values.

p_r*1.25

Simulate Reproduction

To simulate reproduction, we have developed the function sim_repr, which uses the individual's p_r, multiplied by an always used reproduction probability scaling parameter (rpsp), as the probability parameter of a poisson distribution ($\lambda$) from which we sample from to the determine the amount of offspring an produced in a particular moment in time. Using a Poisson distribution ensures that the sampling output is a 0 or positive integer, to illustrate, below we plot the probability mass function of a Poisson distribution of $\lambda$ = pr. For our p_r the most likely outcomes are 0 or 1, but it is possible for an individual to have >= 2 offspring at a particular moment in time.

x <- seq(0,8, by = 1)
plot(x, dpois(x, lambda=p_r), type='b', pch = 19, xlab = "Number of Offspring Generated", ylab = "Outcome Probability", ylim = c(0,1), lwd = 2)

In this case rpsp becomes essential since, realistically, reproduction is highly unlikey over a short period of time like a day. Setting rpsp = 0.05 results in a much more realistic distribution of potential offspring generated being 0 offspring, by far, the most likely outcome.

plot(x, dpois(x, lambda=p_r), type='b', pch = 19, xlab = "Number of Offspring Generated", ylab = "Outcome Probability", ylim = c(0,1), lwd = 2)
lines(x,dpois(x, lambda = p_r*0.05), type = "b", pch = 19, col = "red", lwd = 2)
legend(x = 6, y = 1, legend = c("lambda = p_r", "lambda = p_r * rpsp"), col = c("black", "red"), lty = c(1,1), cex = 0.75, box.lty = 0, lwd = 4)

Simulate Survival and Reproduction throughout a Generation

In order to make our simulations more realistic, we have developed an complementary functions to the both the sim_surv and sim_repr that allow computing both survival and reproduction throughout an entire generation. First, using the sim_gen_tseq we can generate some temperature data for an example 50 time steps (day equivalent) generation (For more details on the simulation of temperature data check the Simulate Temperature Data vignette).

t_seq <- sim_gen_tseq(mean_te = 29, sd_te = 1.5, n_days = 50)
plot(t_seq, type = "l", lwd = 2, xlab = "Day", ylab = "Temperature")

Using this temperature sequence we can simulate the survival of a individual throughout the entire generation using the sim_gen_surv function. sim_gen_surv applies the previously presented sim_surv function multiple times but whenever the individual dies (surv = 0) all future survival is kept at 0. Below we show how the function works:

generation_survival <- sim_surv_gen(tpc = tpc, tseq = t_seq, spsp = 1.25, pmin = 0.1)
generation_survival

par(mar =c(3,4,3,4))
plot(t_seq, type = "l", lwd = 2, xlab = "Day", ylab = "Temperature", col = "grey")
par(new = TRUE)
plot(generation_survival, type = "b", lwd = 2, xaxt = "n", yaxt = "n", xlab = "", ylab = "", pch = 19, col = "royalblue")
mtext("Survival", side = 4, line = 3)
axis(4, ylim=c(0,1), at = c(0,1), col="black",col.axis="black",las=1)
legend("top", inset = c(0,-0.1), bty = "n", legend = c("Temperature", "Survival"), col = c("grey", "royalblue"), lty = c(1,1), cex = 0.75, lwd = 4, horiz = TRUE, bg = "transparent", xpd = TRUE)

In those time steps in which the individual is alive it can also reproduce. Using the sim_gen_repr we can determine the amount of offspring produced throughout the generation. Below we show how the function works, increasing rpsp up to the totally unrealistic 0.5 since otherwise, as in realistic circumstances, reproduction would be extremely unrealistic.

generation_reproduction <- sim_repr_gen(tpc = tpc, tseq = t_seq, surv = generation_survival, rpsp = 0.5, pmin = 0.1)
generation_reproduction

gen_repr_plot <- c(generation_reproduction, rep(NA,length(generation_survival)-length(generation_reproduction)))

par(mar =c(3,4,3,4))
plot(t_seq, type = "l", lwd = 2, xlab = "Day", ylab = "Temperature", col = "grey")
par(new = TRUE)
plot(generation_survival, type = "b", lwd = 2, xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = "lightblue", pch = 19)
par(new = TRUE)
plot(gen_repr_plot, type = "b", lwd = 2, xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = "red", pch = 19, ylim = c(0,max(generation_reproduction)))
mtext("Reproduction", side = 4, line = 3, col = "black")
axis(4, ylim=c(0,max(generation_reproduction)), col="black",col.axis="black",las=1)
legend("top", inset = c(0,-0.1), bty = "n", legend = c("Temperature", "Survival", "Reproduction"), col = c("grey", "lightblue","red"), lty = c(1,1,1), cex = 0.75, lwd = 4, horiz = TRUE, bg = "transparent", xpd = TRUE)