## Code for model in Appendix 2 article entitled "Phytophagous insect oviposition shifts in response to probability of flower abortion owing to presence of basal fruits" ## Authors: Shivani Jadeja and Brigitte Tenhumberg ## Journal Ecology and Evolution # We developed an Individual-Based Model (IBM) to support our prediction that with increasing number of basal fruits fewer yucca moth larvae will emerge from distal fruits because of fewer ovipositions in those flowers with increasing number of basal fruits. # Libraries --------------------------------------------------------------- library(ggplot2) library(Hmisc) library(gridExtra) # Model description/assumptions ------------------------------------------- # all inflorescences start flowering at time step 1 # time in flowering season is correlated with flower position # each flower is available for pollination/ovipostion for one time step # if a flower receives at least one oviposition, it is considered pollinated. # we assume all eggs laid develop into larvae, if flower is retained # we assume once flowers are retained, all fruits survive to maturation # we assume probability of survival of flowers, eggs or fruits if any will only shift the intercept and not the pattern, hence we ignore such scaling parameters # function for probability of oviposition/retention ----------------------- prob_func = function(E, frbelow, n = 1, b = 1, fl = 30){ exp(-(E/n)^b) * (1 - frbelow/(fl-1))} # function for zero-truncated poisson distribution (from https://stat.ethz.ch/pipermail/r-help/2005-May/070684.html) rtpois <- function(N, lambda) qpois(runif(N, dpois(0, lambda), 1), lambda) # Parameters -------------------------------------------------------------- infl = 100 # number of inflorescences fl = 30 # number of flowers on an inflorescence meanmoths = 5 # mean number of moths that arrive at an inflorescence (for a poison distribution) novi = 12 # value of coefficient n for logistic function for probability of oviposition given number of prior ovipositions bovi = 3 # value of coefficient b for logistic function for probability of oviposition given number of prior ovipositions nret = 20 # value of coefficient n for logistic function for probability of flower retention given number of prior ovipositions bret = 5 # value of coefficient b for logistic function for probability of flower retention given number of prior ovipositions # Figures for functions used in model ------------------------------------- ovicurves = data.frame(expand.grid(seq(0, length = 30), 0:(fl-1))) colnames(ovicurves) = c("E", "frbelow") ovicurves$oprob = prob_func(ovicurves$E, 0, novi, bovi, fl) p1a = ggplot(data = ovicurves, aes(x = E+1, y = oprob)) + geom_line(aes(group = frbelow,color = rev(frbelow) )) + xlab("Prior ovipositions") + ylab("Oviposition probability")+ ggtitle("(a)") + theme_bw() + theme(panel.grid = element_blank(), panel.background = element_blank(),legend.position = "none", plot.title = element_text(size = 14, face = "plain")) ovicurves = data.frame(expand.grid(seq(0, length = 30), 0:(fl-1))) colnames(ovicurves) = c("E", "frbelow") ovicurves$oprob = prob_func(ovicurves$E, ovicurves$frbelow, novi, bovi, fl) p1b = ggplot(data = ovicurves, aes(x = E+1, y = oprob)) + geom_line(aes(group = frbelow,color = frbelow) ) + xlab("Prior ovipositions") + ylab("Oviposition probability")+ ggtitle("(b)") + theme_bw() + theme(panel.grid = element_blank(), panel.background = element_blank(), plot.title = element_text(size = 14, face = "plain"),axis.title.y=element_blank(),axis.text.y=element_blank())+ scale_color_continuous(name="Basal fruits") ggsave("FigS2_1ab.tiff",grid.arrange(p1a, p1b, ncol = 2, widths = c(4.5,6)), width = 5, height = 3, dpi = 600) ovimeans = data.frame(frbelow = 0:(fl-1)) ovimeans$mean = 6 p3a = ggplot(data = ovimeans, aes(x = frbelow, y = mean)) + geom_line() + xlab("Number of basal fruits") + ylab("Mean ovipositions")+ ylim(c(0,6)) + ggtitle("(a)") + theme_bw() + theme(panel.grid = element_blank(), panel.background = element_blank(), plot.title = element_text(size = 14, face = "plain")) ovimeans = data.frame(frbelow = 0:(fl-1)) ovimeans$mean = 6*(1 - ovimeans$frbelow/(fl-1)) p3b = ggplot(data = ovimeans, aes(x = frbelow, y = mean)) + geom_line() + ylim(c(0,6))+ xlab("Number of basal fruits") + ylab("Mean ovipositions")+ ggtitle("(b)") + theme_bw() + theme(panel.grid = element_blank(), panel.background = element_blank(), legend.position = c(0.9, 0.8), plot.title = element_text(size = 14, face = "plain"),axis.title.y=element_blank(),axis.text.y=element_blank()) ggsave("FigS2Xab.tiff",grid.arrange(p3a, p3b, ncol = 2, widths = c(5,4.5)), width = 5, height = 3, dpi = 600) ##### retcurves = data.frame(expand.grid(seq(0, length = 30), 0:(fl-1))) colnames(retcurves) = c("E", "frbelow") retcurves$retprob = prob_func(retcurves$E, retcurves$frbelow, nret, bret, fl) p2 = ggplot(data = retcurves, aes(x = E+1, y = retprob, col = frbelow)) + geom_line(aes(group = frbelow)) + xlab("Number of prior ovipositions") + ylab("Retention probability") + theme_bw() + theme(panel.grid = element_blank(), panel.background = element_blank(), plot.title = element_text(size = 14, face = "plain"))+ scale_color_continuous(name="Basal fruits") ggsave("FigS2_2.tiff",p2, width = 4, height = 3, dpi = 600) # Simulation function ----------------------------------------------------- simfunc = function(infl = 10, fl = 30, meanmoths = 5,novi = 12, bovi = 3, nret = 20, bret = 5, probdep = TRUE, numdep = TRUE, label = NULL){ # object for recording information for each flower position for each inflorescence outovi = array(rep(0,infl*fl), dim=c(fl,infl)) # columns are infl and rows are flower positions with row 1 being the bottom first flower # IBM for(f in 1:fl){ # each time step one flower opens starting with the bottom most flower which is stored in row 1 for(i in 1:infl){ # loop over all inflorescences sumovi = 0 # object to store ovipositons received by flower, if any basal = ifelse(f == 1, 0, length(outovi[1:(f-1),i][outovi[1:(f-1),i]>0])) # determining number of fruits below M = rpois(1,meanmoths) # number of moths that arrive at flower if(M != 0){ # proceed if non-zero moths arrive at flower then proceed, else flower receives zero ovipositions and aborts for(m in 1:M){ # loop over the number of moths that arrive at the flower if(rbinom(1,1, prob = prob_func(sumovi,ifelse(probdep == TRUE,basal,0),novi, bovi, fl)) == 1){ # probablity moth m oviposits sumovi = sumovi + rtpois(1,ifelse(numdep == TRUE, 6*(1 - basal/(fl-1)),6)) # adds moth m's ovipositions to the object where the flower's total ovipositions are recorded if(sumovi > 0 & rbinom(1,1, prob = prob_func(sumovi,basal, nret, bret, fl)) == 1){ # if flower received at least one oviposition, and plant retained flower based on its probability of retention. outovi[f,i] = sumovi}}}}}} # records total ovipositions from flower f # output for figure to match figure from empirical data ibmdat = data.frame(infl = NA, fpos = NA, ovi = NA, frbelow = NA) for(i in 1:infl){ # loop over all inflorescneces for(f in (round(fl*2/3,0)+1):fl){ # loop over all fruits from TOP third flowers on the inflorescences if(outovi[f,i] > 0){ ibmdat = rbind(ibmdat, c(i,f, outovi[f,i], length(outovi[1:(f-1),i][outovi[1:(f-1),i]>0])))}}} # if fruit was developed, extract number of oviposition and fruits below for that flower ibmdat = ibmdat[complete.cases(ibmdat),] # removing first row with NA for(i in 1:dim(ibmdat)[1]){ # adding a column to get frequency of points ibmdat$fruits[i] = length(which(ibmdat$ovi == ibmdat$ovi[i] & ibmdat$frbelow == ibmdat$frbelow[i])) } # calculating mean, median, and upper and lower quantiles for simulated values of number of basal fruits and number of ovipositions predout = data.frame(frbelow = NA, mean = NA, median = NA, lower = NA, upper = NA) # to store meams, medians and 1st and 4th quantiles of number of ovipositions at each value of basal fruits for(n in unique(ibmdat$frbelow)){ predout = rbind(predout, c(n, mean(ibmdat$ovi[ibmdat$frbelow == n]), median(ibmdat$ovi[ibmdat$frbelow == n]), quantile(ibmdat$ovi[ibmdat$frbelow == n],0.025), quantile(ibmdat$ovi[ibmdat$frbelow == n],0.975)))} predout = predout[complete.cases(predout),] # removing first row with NA ggplot(data = ibmdat, aes(x = frbelow)) + geom_point(aes(y = ovi, size = fruits), col = "violetred1") + geom_point(data = predout, aes(x = frbelow, y = median)) + geom_errorbar(data = predout, aes(x = frbelow, ymin = lower, ymax = upper), width = 0.2) + xlab("Number of basal fruits") + ylab("Number of eggs") + theme_bw() + ggtitle(label) + theme(panel.grid = element_blank(), panel.background = element_blank(), plot.title = element_text(size = 14, face = "plain"))+ scale_size_continuous(name="Frequency") } # running function for four types of simulations a = simfunc(infl = 10000, fl = 30, meanmoths = 5, novi = 12, bovi = 3, nret = 20, bret = 5, probdep = FALSE, # oviposition probability dependent on basal fruits (TRUE/FALSE) numdep = FALSE, # number of ovipositions dependent on basal fruits (TRUE/FALSE) label = expression(paste("(a) ",P[ovi],"(E) and ", ~lambda, " = 6"))) b = simfunc(infl = 10000, fl = 30, meanmoths = 5, novi = 12, bovi = 3, nret = 20, bret = 5, probdep = TRUE, numdep = FALSE, label = expression(paste("(b) ",P[ovi],"(E,B) and ", ~lambda, " = 6"))) c = simfunc(infl = 10000, fl = 30, meanmoths = 5, novi = 12, bovi = 3, nret = 20, bret = 5, probdep = FALSE, numdep = TRUE, label = expression(paste("(c) ",P[ovi],"(E) and ", ~lambda, "(B)"))) d = simfunc(infl = 10000, fl = 30, meanmoths = 5, novi = 12, bovi = 3, nret = 20, bret = 5, probdep = TRUE, numdep = TRUE, label = expression(paste("(d) ",P[ovi],"(E,B) and ", ~lambda, "(B)"))) ggsave("FigS2_3abcd.tiff", grid.arrange(a, b, c, d, ncol = 2), width =6, height = 4, dpi = 600)