# Manuscript Title: Dynamic regulation of partner abundance mediates response of reef coral # symbioses to environmental change # Authors: R Cunning, N Vaughan, P Gillette, T Capo, JL Maté, AC Baker # Description: This R Script fits the cost-benefit model structure # to empirical data by adjusting scaling parameters #-------------------------------------------------------------------------------------------------- # Import data: observed symbiont to host cell ratios and standard deviations # and environmental metadata data <- read.csv("Pdam_dataset.csv") metadata <- read.csv("Pdam_metadata.csv") means <- aggregate(data$logTotal, FUN=mean, by=list(Experiment=data$Experiment, Time=data$Time, Sym=data$Sym)) means <- means[with(means, order(Sym, Experiment, Time)), ] sds <- aggregate(data$logTotal, FUN=sd, by=list(Experiment=data$Experiment, Time=data$Time, Sym=data$Sym)) sds <- sds[with(sds, order(Sym, Experiment, Time)), ] optz.obs <- data.frame(Experiment=rep(metadata$Experiment, 2), Time=rep(metadata$Time, 2), Sym=c(rep("C", 5), rep("D", 5))) optz.obs <- optz.obs[with(optz.obs, order(Sym, Experiment, Time)), ] optz.obs$lower <- 10^(means$x - sds$x) optz.obs$mean <- 10^means$x optz.obs$upper <- 10^(means$x + sds$x) optz.obs # Define match.numeric function for use in fitting model match.numeric <- function(x, table) { are.equal <- function(x, y) isTRUE(all.equal(x, y)) match.one <- function(x, table) match(TRUE, vapply(table, are.equal, logical(1L), x=x)) vapply(x, match.one, integer(1L), table) } #-------------------------------------------------------------------------------------------------- # Define cost-benefit model # Define environmental range t.max <- 32 t.min <- 15 i.max <- 30 i.min <- 1 # Define variable simulation space i.step <- 0.1 t.step <- 0.1 z.step <- 0.001 # Environmental "windows" where data were collected (means ± se from Pdam_metadata.csv) i.windows <- c(3.6, seq(5, 12.6, i.step), seq(13, 14.2, i.step)) t.windows <- c(seq(22.2, 23.2, t.step), seq(23.9, 24.7, t.step), seq(25.5, 26.5, t.step), seq(27.9, 28.4, t.step)) z <- seq(0, 0.2, by=z.step) # Define symbiont types (biotic factors) cladeC <- data.frame(t.opt=28, i.opt=10, v=1, f=1) cladeD <- data.frame(t.opt=28, i.opt=16, v=0.4778, f=0.5) clades <- list(C=cladeC, D=cladeD) # Define fast model structure CostBenModFast <- function(clade, params) { i.opt <- (clade$i.opt - i.min) / (i.max - i.min) t.opt <- (clade$t.opt - t.min) / (t.max - t.min) i.sub <- function(i.amb) ifelse(i.amb > i.opt, 0, (i.amb - i.opt)^2) i.super <- function(i.amb) ifelse(i.amb > i.opt, (i.amb - i.opt)^2, 0) t.sub <- function(t.amb) ifelse(t.amb > t.opt, 0, (t.amb - t.opt)^2) t.super <- function(t.amb) ifelse(t.amb > t.opt, (t.amb - t.opt)^2, 0) S <- outer(i.amb, t.amb, function(i.amb, t.amb) exp(params$b7 * (clade$f * (params$b2 * t.sub(t.amb) + params$b3 * t.super(t.amb)) + params$b8 * (i.amb / i.opt)^2))) P <- params$b1 * clade$f / (outer(i.amb, t.amb, function(i.amb, t.amb) 1 + (clade$f * (params$b2 * t.sub(t.amb) + params$b3 * t.super(t.amb)) + i.super(i.amb)))) B.net <- (aperm(vapply(1:length(t.amb), FUN=function(x) outer(P, z)[, x, ] * exp(outer(-(params$b4 * clade$v * (1 + params$b5 * i.sub(i.amb))), z)), FUN.VALUE=array(0, dim=c(length(i.amb), length(z)))), c(1, 3, 2))) - ((outer(rep(1, length(i.amb)), params$b6 * clade$f * outer(t.amb, z))) + (clade$f * outer(S, z) * outer(S, z, function(S, z) exp(params$b9 * (z - params$b10 / S))))) result <- list(optz=apply(B.net, c(1, 2), function(x) z[which.max(x)])) return(result) } # Set up fast model simulation space i.amb <- (i.windows - i.min) / (i.max - i.min) t.amb <- (t.windows - t.min) / (t.max - t.min) # Create data frame to populate with modeled optimal densities optz.mod <- cbind(rbind(metadata, metadata), Sym=c(rep("C", 5), rep("D", 5))) optz.mod$i.lower <- round(optz.mod$i.mean - optz.mod$i.se, 1) optz.mod$i.upper <- round(optz.mod$i.mean + optz.mod$i.se, 1) optz.mod$t.lower <- round(optz.mod$t.mean - optz.mod$t.se, 1) optz.mod$t.upper <- round(optz.mod$t.mean + optz.mod$t.se, 1) optz.mod$i.lower.s <- match.numeric((optz.mod$i.lower - i.min) / (i.max - i.min), i.amb) optz.mod$i.upper.s <- match.numeric((optz.mod$i.upper - i.min) / (i.max - i.min), i.amb) optz.mod$t.lower.s <- match.numeric((optz.mod$t.lower - t.min) / (t.max - t.min), t.amb) optz.mod$t.upper.s <- match.numeric((optz.mod$t.upper - t.min) / (t.max - t.min), t.amb) optz.mod <- optz.mod[with(optz.mod, order(Sym, Experiment, Time)), ] #-------------------------------------------------------------------------------------------------- # Define function to evaluate scaling parameters Evaluate <- function(Values) { params <- data.frame(b1 = Values[1], b2 = Values[2], b3 = Values[3], b4 = Values[4], b5 = Values[5], b6 = Values[6], b7 = Values[7], b8 = Values[8], b9 = Values[9], b10 = Values[10]) # Scaling parameter starting values model <- lapply(clades, CostBenModFast, params=params) modz <- apply(optz.mod, 1, function (y) mean(model[[as.character(y[[7]])]]$optz[y[12]:y[13], y[14]:y[15]])) total.dens <- sum(dnorm(log10(modz), means$x, sds$x)) return(total.dens) } #start<-c(39.59615, 0.5767399, 10, 4.717651, 46.45819, 16.046361, 25, 0.023545237, 28.43473, 0.2793839) #start<-c(39.58814, 0.5805317, 10, 4.608033, 46.69023, 16.053237, 25, 0.023291565, 28.55235, 0.2652132) #start<-c(39.61344, 0.5803846, 10, 4.631671, 46.96926, 16.065083, 24.960921, 0.023395223, 28.37715, 0.2808066) # Define the starting values, ranges, step sizes, and search momentum for each scaling parameter # used in the gradient search Inputs <- data.frame(start=c(50, 1.8, 10, 6, 55, 7, 10, 0.08, 20, 0.32), min=c(25, 0.2, 10, 1, 5, 1, 2, 0.001, 5, 0.05), max=c(75, 4, 10, 20, 100, 20, 25, 0.5, 50, 1), step=c(1, 0.1, 0.1, 0.5, 1, 0.2, 0.2, 0.01, 0.5, 0.05), momentum=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) # Set the number of search iterations Iterations <- 1000 # Once the fitting reaches half the total iterations the step size starts to shrink to # help the algorithm converge. Decay.Start <- Iterations / 2 # Trace functions provide a record of the goodness of fit during the fitting search # and allow for confirmation of convergence. Trace <- matrix(nrow=(Iterations + 1), ncol=length(Inputs[, 1])) Trace.Eval <- vector(length=(Iterations * length(Inputs[, 1]))) Trace.Eval[1] <- Evaluate(Inputs[, 1]) Opt.Estim <- c(Trace.Eval[1], 50, 1.8, 10, 6, 55, 7, 10, 0.08, 20, 0.32) for (p in 1:length(Inputs[, 1])) Trace[1, p] <- Inputs[p, 1] # Set counters to zero to start fitting routine # ID references the specific parameter being updated ID <- 1 Iteration <- 0 Run <- TRUE # Run the fitting algorithm while (Run) { # Evaluate the step size decay value if (Iteration < Decay.Start) { decay <- 1 } else { decay <- (Iterations - Iteration) / (Iterations - Decay.Start) } # Calculate model fit for one step plus momentum in each direction of the current parameter. vals1 <- Inputs[, 1] vals2 <- Inputs[, 1] vals1[ID] <- max(Inputs[ID, 2], min(Inputs[ID, 3], (Inputs[ID, 1] - Inputs[ID, 4] * decay + Inputs[ID, 5]))) Estim1 <- Evaluate(vals1) vals2[ID] <- max(Inputs[ID, 2], min(Inputs[ID, 3], (Inputs[ID, 1] + Inputs[ID, 4] * decay + Inputs[ID, 5]))) Estim2 <- Evaluate(vals2) # Determine the best step direction corr <- Estim1 - Estim2 if (corr > 0) { va <- -0.1 * Inputs[ID, 4] } else if (corr < 0) { va <- 0.1 * Inputs[ID, 4] } else { va <- 0 } # Update parameter momentum Inputs[ID, 5] <- 0.8 * Inputs[ID, 5] + va * decay # Calculate a weighted mean of the two bounded estimates to determine if they straddle an optimum vals3 <- Inputs[, 1] vals3[ID] <- (vals1[ID]*Estim1+vals2[ID]*Estim2)/(Estim1+Estim2) Estim3 <- Evaluate(vals3) # Pick the best fit and update the estimated value for the current parameter Best.Estim <- max(Estim1, Estim2, Estim3) if(Estim1==Best.Estim){Inputs[ID, 1] <- vals1[ID] }else if(Estim2==Best.Estim){Inputs[ID, 1] <- vals2[ID] }else if(Estim3==Best.Estim){Inputs[ID, 1] <- vals3[ID]} # Iterate search trace vectors and check for new optimum estimate Trace[(Iteration + 2), ID] <- Inputs[ID, 1] Trace.Eval[(Iteration * length(Inputs[, 1]) + ID + 1)] <- Best.Estim if(Best.Estim>Opt.Estim[1]) { Opt.Estim[1]<-Best.Estim for(i in 1:length(Inputs[,1])) { Opt.Estim[i+1]<-Inputs[i,1] } } # Increment iterations and cycle parameter ID if (ID == length(vals1)) { Iteration <- Iteration + 1 } if (Iteration >= Iterations) { Outputs <- Inputs Run <- FALSE } else { if (ID < length(Inputs[, 1])) { ID <- ID + 1 } else { ID <- 1 } } } # Write parameter value trace (Trace), fit trace (Trace.Eval), and fitted parameter estimates (Opt.Estim) to file write.table(Trace, file="Trace.csv", row.names=FALSE, col.names=FALSE, append=FALSE, qmethod="double", dec=".", sep=",") write.table(Trace.Eval, file="TraceEval.csv", row.names=FALSE, col.names=FALSE, append=FALSE, qmethod="double", dec=".", sep=",") write.table(Opt.Estim, file="OptEstim.csv", row.names=FALSE, col.names=FALSE, append=FALSE, qmethod="double", dec=".", sep=",")