# 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 runs the cost-benefit model with fitted scaline parameters # and produces all of the Figures in the manuscript #---------------------------------------------------------- # Import data: observed symbiont to host cell ratios 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 #-------------------------------------------------------------------------------------------------- # 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 i.full <- seq(i.min, i.max, by=i.step) # full range t.full <- seq(t.min, t.max, by=t.step) # full range 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 model structure and output CostBenMod <- 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) light.capture <- exp(outer(-(params$b4 * clade$v * (1 + params$b5 * i.sub(i.amb))), z)) L <- 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))) 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.gross <- 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)) C.maint <- outer(rep(1, length(i.amb)), params$b6 * clade$f * outer(t.amb, z)) C.stress <- clade$f * outer(S, z) * outer(S, z, function(S, z) exp(params$b9 * (z - params$b10 / S))) C.gross <- C.maint + C.stress B.net <- B.gross - C.gross B.max <- apply(B.net, c(1, 2), max) optz <- apply(B.net, c(1, 2), function(x) z[which.max(x)]) result <- list(optz=optz, light.capture=light.capture, L=L, S=S, P=P, B.gross=B.gross, C.maint=C.maint, C.stress=C.stress, C.gross=C.gross, B.net=B.net, B.max=B.max) return(result) } # Set up full model simulation space i.amb <- (i.full - i.min) / (i.max - i.min) t.amb <- (t.full - t.min) / (t.max - t.min) # Define match.numeric function for use 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) } # 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 scaling parameters and run model # Scaling parameter values used here are from 1000 iterations of # the gradient seach routine implemented in accompanying R script: Pdam_fitModel.R finalparams <- data.frame(b1 = 39.288, b2 = 0.577, b3 = 10, b4 = 4.666, b5 = 46.657, b6 = 16.131, b7 = 24.911, b8 = 0.0234, b9 = 28.377, b10 = 0.281) model <- lapply(clades, CostBenMod, params=finalparams) modz <- apply(optz.mod, 1, function (y) { mean(model[[as.character(y[[7]])]]$optz[y[12]:y[13], y[14]:y[15]]) }) # Compare observed and modeled mean densities df <- data.frame(cbind(optz.obs[, c(1, 2, 3, 5)], modz)); df #-------------------------------------------------------------------------------------------- # FIGURES #-------------------------------------------------------------------------------------------- t <- seq(t.min, t.max, by=t.step) i <- seq(i.min, i.max, by=i.step) z <- seq(0, 0.2, by=z.step) #-------------------------------------------------------------------------------------------- #Figure 1 graphics.off() par(mfrow=c(1, 3)) #symbiont self shading (Fig. 2A) par(mar=c(3, 3, 3, 1)) plot(1, type="n", mgp=c(1, 1, 0), cex.lab=1.2, xlim=c(0, 0.2), xaxt="n", xlab="Sym. abundance", xaxs="i", ylim=c(1, 30), yaxt="n", ylab="Irradiance", yaxs="i") .filled.contour(t(model$D$light.capture), x=z, y=i, levels=seq(0, 1, 0.05), col=gray.colors(20, 0.2, 1)) title("B. Self shading", adj=0, cex.main=1.5) #symbiont productivity limitation (Fig. 2B) par(mar=c(3, 2.5, 3, 1.5)) plot(1, type="n", mgp=c(1, 1, 0), cex.lab=1.2, xlim=c(15, 32), xaxt="n", xlab="Temperature", xaxs="i", ylim=c(1, 30), yaxt="n", ylab="Irradiance", yaxs="i") .filled.contour(t(model$D$L), x=t, y=i, col=gray.colors(12, 1, 0), levels=seq(1, 1.55, 0.05)) title("C. Productivity limitation", adj=0, cex.main=1.5) range(t(model$D$L)) #ROS production (Fig. 2C) par(mar=c(3, 2, 3, 2)) plot(1, type="n", mgp=c(1, 1, 0), cex.lab=1.2, xpd=T, xlim=c(15, 32), xaxt="n", xlab="Temperature", xaxs="i", ylim=c(1, 30), yaxt="n", ylab="Irradiance", yaxs="i") .filled.contour(log(t(model$D$S)), x=t, y=i, levels=seq(1, 10, 0.5), col=gray.colors(19, 1, 0)) title("D. Oxidative stress", adj=0, cex.main=1.5) #-------------------------------------------------------------------------------------------- # FIGURE 2 graphics.off() dl <- split(data, interaction(data$Experiment, data$Sym, data$Time), drop=T) pdfs <- lapply(dl, function(x) density(x$logTotal, adjust=1.5)) par(mfrow=c(2, 2)) # Experiment 1: C corals par(mar=(c(2, 3, 2, 1))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(-1.5, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=F) axis(side=2, at=c(0, 2), labels=c(0, 2), las=1) lines(pdfs$"1.C.field", col="gray40") lines(pdfs$"1.C.lab", col="orange") boxplot(dl$"1.C.field"$logTotal, range=0, horizontal=T, add=T, at=-1, xaxt="n", frame=F) points(dl$"1.C.field"$logTotal, rep(-1, length(dl$"1.C.field"$logTotal)), pch=18, col="gray40") boxplot(dl$"1.C.lab"$logTotal, range=0, horizontal=T, add=T, at=-0.5, xaxt="n", frame=F) points(dl$"1.C.lab"$logTotal, rep(-0.5, length(dl$"1.C.field"$logTotal)), pch=18, col="orange") text(c(-2.35, -2.35), c(-0.5, -1), c("lab", "field"), pos=4) text(-2.35, 2.9, "Exp. 1: C corals", pos=4, font=2, xpd=T) # Experiment 2: C corals par(mar=(c(2, 1, 2, 3))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(-1.5, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=F) axis(side=2, at=c(0, 2), labels=F) lines(pdfs$"2.C.feb", col="DodgerBlue3") lines(pdfs$"2.C.apr", col="DarkGreen") lines(pdfs$"2.C.jun", col="firebrick") boxplot(dl$"2.C.feb"$logTotal, range=0, horizontal=T, add=T, at=-1.25, xaxt="n", frame=F) points(dl$"2.C.feb"$logTotal, rep(-1.25, length(dl$"2.C.feb"$logTotal)), pch=18, col="DodgerBlue3") boxplot(dl$"2.C.apr"$logTotal, range=0, horizontal=T, add=T, at=-0.75, xaxt="n", frame=F) points(dl$"2.C.apr"$logTotal, rep(-0.75, length(dl$"2.C.apr"$logTotal)), pch=18, col="DarkGreen") boxplot(dl$"2.C.jun"$logTotal, range=0, horizontal=T, add=T, at=-0.25, xaxt="n", frame=F) points(dl$"2.C.jun"$logTotal, rep(-0.25, length(dl$"2.C.jun"$logTotal)), pch=18, col="firebrick") text(c(-2.35, -2.35, -2.35), c(-1.25, -0.75, -0.25), c("Feb.", "Apr.", "Jun."), pos=4) text(-2.35, 2.9, "Exp. 2: C corals", pos=4, font=2, xpd=T) # Experiment 1: D corals par(mar=(c(3, 3, 1, 1))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(-1.5, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T), mgp=c(3, 0.5, 0)) axis(side=2, at=c(0, 2), labels=c(0, 2), las=1) lines(pdfs$"1.D.field", col="gray40") lines(pdfs$"1.D.lab", col="orange") boxplot(dl$"1.D.field"$logTotal, range=0, horizontal=T, add=T, at=-1, xaxt="n", frame=F) points(dl$"1.D.field"$logTotal, rep(-1, length(dl$"1.D.field"$logTotal)), pch=18, col="gray40") boxplot(dl$"1.D.lab"$logTotal, range=0, horizontal=T, add=T, at=-0.5, xaxt="n", frame=F) points(dl$"1.D.lab"$logTotal, rep(-0.5, length(dl$"1.D.field"$logTotal)), pch=18, col="orange") text(c(-2.35, -2.35), c(-0.5, -1), c("lab", "field"), pos=4) text(-2.35, 2.9, "Exp. 1: D corals", pos=4, font=2, xpd=T) # Experiment 2: D corals par(mar=(c(3, 1, 1, 3))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(-1.5, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T), mgp=c(3, 0.5, 0)) axis(side=2, at=c(0, 2), labels=F) lines(pdfs$"2.D.feb", col="DodgerBlue3") lines(pdfs$"2.D.apr", col="DarkGreen") lines(pdfs$"2.D.jun", col="firebrick") boxplot(dl$"2.D.feb"$logTotal, range=0, horizontal=T, add=T, at=-1.25, xaxt="n", frame=F) points(dl$"2.D.feb"$logTotal, rep(-1.25, length(dl$"2.D.feb"$logTotal)), pch=18, col="DodgerBlue3") boxplot(dl$"2.D.apr"$logTotal, range=0, horizontal=T, add=T, at=-0.75, xaxt="n", frame=F) points(dl$"2.D.apr"$logTotal, rep(-0.75, length(dl$"2.D.apr"$logTotal)), pch=18, col="DarkGreen") boxplot(dl$"2.D.jun"$logTotal, range=0, horizontal=T, add=T, at=-0.25, xaxt="n", frame=F) points(dl$"2.D.jun"$logTotal, rep(-0.25, length(dl$"2.D.jun"$logTotal)), pch=18, col="firebrick") text(c(-2.35, -2.35, -2.35), c(-1.25, -0.75, -0.25), c("Feb.", "Apr.", "Jun."), pos=4) text(-2.35, 2.9, "Exp. 2: D corals", pos=4, font=2, xpd=T) mtext("Symbiont abundance (S/H cell ratio)", side=1, line=-1, outer=T) mtext("Relative occurrence (kernel density estimation)", side=2, line=-1, outer=T) #-------------------------------------------------------------------------------------------- #Figure 3: Costs and benefits at selected levels of temperature and irradiance par(oma=c(4, 4, 2, 4)) par(mfrow=c(3, 3)) bottom=F left=T t.Series=c(23, 26, 29) # Set temperature values i.Series=c(15, 10, 5) # Set irradiance values t.Series <- match.numeric((t.Series - t.min) / (t.max - t.min), t.amb) i.Series <- match.numeric((i.Series - i.min) / (i.max - i.min), i.amb) # Set clade here clade <- 1 # 1=clade C, 2=clade D # Run following code for each clade vect <- vector(length=length(z)) for(j in 1:3) { if(j==3) {bottom=T} for(k in 1:3) { if(k==1) {left=T} par(yaxt="n", xaxt="n") par(mar=c(0, 0, 0, 0), lwd=1) for(l in 1:length(z)) {vect[l] <- model[[clade]]$B.gross[i.Series[j], t.Series[k], l]} plot(vect ~ z, pch=NA, ann=FALSE, xlim=c(0, 0.2), ylim=c(0, 1.5), type="n") #Add standardized axis to the plot window par(yaxt="s", xaxt="s") par(tcl=-0.25) if(left!=F)axis(side=2, lwd=0, lwd.ticks=1, at=seq(0, 2, by=0.5), labels=c("0", "0.5", "1", "1.5", "2")) if(bottom!=F)axis(side=1, lwd=0, lwd.ticks=1, at=seq(0, 0.2, by=0.05), labels=c("0", "", "0.1", "", "0.2")) #Add benefit and cost curves lines(vect ~ z, col="green", lty=1) for(l in 1:length(z)) {vect[l] = model[[clade]]$C.gross[i.Series[j], t.Series[k], l]} lines(vect ~ z, col="red", lty=1) for(l in 1:length(z)) {vect[l] <- model[[clade]]$B.net[i.Series[j], t.Series[k], l]} lines(vect ~ z, col="blue", lty=1) #Add line and label for maximum net benefit X1 <- seq(-1, model[[clade]]$optz[i.Series[j], t.Series[k]], by=0.0001) Y1 <- c(1:length(X1)) for(i in 1:length(Y1)) {Y1[i] <- model[[clade]]$B.max[i.Series[j], t.Series[k]]} par(tcl=-0.25) Y <- round(model[[clade]]$B.max[i.Series[j], t.Series[k]], 3) axis(side=2, lwd=0, lwd.ticks=2, at=Y, labels=FALSE, cex.axis=0.8, font=2) lines(X1, Y1, lty=2, col="black") text(-0.01, model[[clade]]$B.max[i.Series[j], t.Series[k]] + 0.1, Y, pos=4, font=2) #Add line and label for optimal density Y2 <- seq(-1, model[[clade]]$B.max[i.Series[j], t.Series[k]], by=0.0001) X2 <- c(1:length(Y2)) for(i in 1:length(X2)) {X2[i] <- model[[clade]]$optz[i.Series[j], t.Series[k]]} par(tcl=-0.25) X <- round(model[[clade]]$optz[i.Series[j], t.Series[k]], 3) axis(side=1, lwd=0, lwd.ticks=2, at=X, labels=FALSE, cex.axis=0.8, font=2) lines(X2, Y2, lty=2, col="black") text(model[[clade]]$optz[i.Series[j], t.Series[k]] + 0.0025, 0, X, pos=4, font=2) #Add point at optimal density, maximum net benefit points(model[[clade]]$optz[i.Series[j], t.Series[k]], model[[clade]]$B.max[i.Series[j], t.Series[k]]) } } #-------------------------------------------------------------------------------------------- #Figure 4: Optimal density contour plot t <- seq(t.min, t.max, by=t.step) i <- seq(i.min, i.max, by=i.step) z <- seq(0, 0.2, by=z.step) # Set clade here clade <- 1 #1=clade C, 2=clade D # Run following code for each clade filled.contour(t(model[[clade]]$optz), x=t, y=i, xlim=c(18, 32), ylim=c(1, 30), xlab="Temperature", ylab="Irradiance", main="Optimal Density", col=gray.colors(17, start=1, end=0.1), levels=seq(0, 0.17, by=0.01), plot.axes={axis(1, seq(18, 32, by=1)) axis(2, seq(0, 30, by=5)) rect(metadata[1, ]$t.mean - metadata[1, ]$t.se, metadata[1, ]$i.mean - metadata[1, ]$i.se, metadata[1, ]$t.mean + metadata[1, ]$t.se, metadata[1, ]$i.mean + metadata[1, ]$i.se, border="black", col=rgb(1, 1, 1, 0.3137)) #Experiment 1 initial rect(metadata[2, ]$t.mean - metadata[2, ]$t.se, metadata[2, ]$i.mean - metadata[2, ]$i.se, metadata[2, ]$t.mean + metadata[2, ]$t.se, metadata[2, ]$i.mean + metadata[2, ]$i.se, border="orange", col=rgb(255, 165, 0, 80, maxColorValue=255)) #Experiment 1 final rect(metadata[3, ]$t.mean - metadata[3, ]$t.se, metadata[3, ]$i.mean - metadata[3, ]$i.se, metadata[3, ]$t.mean + metadata[3, ]$t.se, metadata[3, ]$i.mean + metadata[3, ]$i.se, border="blue", col=rgb(0, 0, 1, 0.3137)) #Experiment 2, February rect(metadata[4, ]$t.mean - metadata[4, ]$t.se, metadata[4, ]$i.mean - metadata[4, ]$i.se, metadata[4, ]$t.mean + metadata[4, ]$t.se, metadata[4, ]$i.mean + metadata[4, ]$i.se, border="forestgreen", col=rgb(34, 139, 34, 80, maxColorValue=255)) #Experiment 2, April rect(metadata[5, ]$t.mean - metadata[5, ]$t.se, metadata[5, ]$i.mean - metadata[5, ]$i.se, metadata[5, ]$t.mean + metadata[5, ]$t.se, metadata[5, ]$i.mean + metadata[5, ]$i.se, border="red", col=rgb(1, 0, 0, 0.3137))}) range(t(model[[2]]$optz)) #-------------------------------------------------------------------------------------------- # Figure 5: Observed densities with model predictions modranges <- apply(optz.mod, 1, function (y) { range(model[[as.character(y[[7]])]]$optz[y[12]:y[13], y[14]:y[15]]) }) graphics.off() par(mfrow=c(2, 2)) par(mar=(c(2, 3, 2, 1))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(0, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=F) axis(side=2, at=c(0, 2), labels=c(0, 2), las=1) rect(log10(modranges[, 1][1]), 0, log10(modranges[, 1][2]), 2.5, border=NA, col=rgb(102, 102, 102, 128, maxColorValue=255)) rect(log10(modranges[, 2][1]), 0, log10(modranges[, 2][2]), 2.5, border=NA, col=rgb(255, 165, 0, 128, maxColorValue=255)) lines(pdfs$"1.C.field", col="gray40", lwd=2) lines(pdfs$"1.C.lab", col="orange", lwd=2) par(mar=(c(2, 1, 2, 3))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(0, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=F) axis(side=2, at=c(0, 2), labels=F) rect(log10(modranges[, 4][1]), 0, log10(modranges[, 4][2]), 2.5, border=NA, col=rgb(16, 78, 139, 128, maxColorValue=255)) rect(log10(modranges[, 3][1]), 0, log10(modranges[, 3][2]), 2.5, border=NA, col=rgb(0, 100, 0, 128, maxColorValue=255)) rect(log10(modranges[, 5][1]), 0, log10(modranges[, 5][2]), 2.5, border=NA, col=rgb(178, 34, 34, 128, maxColorValue=255)) lines(pdfs$"2.C.feb", col="DodgerBlue4", lwd=2) lines(pdfs$"2.C.apr", col="DarkGreen", lwd=2) lines(pdfs$"2.C.jun", col="firebrick", lwd=2) par(mar=(c(3, 3, 1, 1))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(0, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T), mgp=c(3, 0.5, 0)) axis(side=2, at=c(0, 2), labels=c(0, 2), las=1) rect(log10(modranges[, 6][1]), 0, log10(modranges[, 6][2]), 2.5, border=NA, col=rgb(102, 102, 102, 128, maxColorValue=255)) rect(log10(modranges[, 7][1]), 0, log10(modranges[, 7][2]), 2.5, border=NA, col=rgb(255, 165, 0, 128, maxColorValue=255)) lines(pdfs$"1.D.field", col="gray40", lwd=2) lines(pdfs$"1.D.lab", col="orange", lwd=2) par(mar=(c(3, 1, 1, 3))) plot(1, type="n", xlim=c(log10(0.006), log10(0.7)), xlab="", xaxt="n", ylim=c(0, 2.5), ylab="", yaxt="n") axis(side=1, at=log10(axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T, nint=5)), labels=axisTicks(c(-3, 0), axp=c(0.001, 0.1, 3), log=T), mgp=c(3, 0.5, 0)) axis(side=2, at=c(0, 2), labels=F) rect(log10(modranges[, 9][1]), 0, log10(modranges[, 9][2]), 2.5, border=NA, col=rgb(16, 78, 139, 128, maxColorValue=255)) rect(log10(modranges[, 8][1]), 0, log10(modranges[, 8][2]), 2.5, border=NA, col=rgb(0, 100, 0, 128, maxColorValue=255)) rect(log10(modranges[, 10][1]), 0, log10(modranges[, 10][2]), 2.5, border=NA, col=rgb(178, 34, 34, 128, maxColorValue=255)) lines(pdfs$"2.D.feb", col="DodgerBlue4", lwd=2) lines(pdfs$"2.D.apr", col="DarkGreen", lwd=2) lines(pdfs$"2.D.jun", col="firebrick", lwd=2)