############################################## ############################################## # Results from: # García-Callejas, D., Molowny-Horas, R. and Araújo, M. B. (2017) # Multiple interactions networks: towards more realistic representations of the web of life # http://onlinelibrary.wiley.com/doi/10.1111/oik.04428/full ########### # contents (by line number) # Expanded Food Webs analyses: # - model with NTIs: 53-264 # - exploratory figures: 268-281 # - model without NTIs: 287-482 # - comparison (Fig. 2 of main text): 488-516 # Equal footing networks analyses: # - common parameters: 522-597 # - simulations with varying interaction strengths: 604-787 # - exploratory analyses: 795-829 # - comparison (Fig. 3 of main text and SM Fig. S7): 835-912 # - Supplementary material Fig. S6: 918-951 ########### # NOTE: # this is a self-contained script, that can be used to replicate (approximately) the results of the paper, # the only requirement being that the referenced packages are installed on the local machine. ########### # contact: David García-Callejas (david.garcia.callejas@gmail.com) ########### rm(list=ls(all=TRUE)) gc() assign("last.warning", NULL, envir = baseenv()) wd <- "/home/david/CREAF/FPU/review/" # wd <- "/your/path/here" setwd(wd) ############################################## ############################################## library(deSolve) library(tidyverse) library(rootSolve) library(scatterplot3d) ############################################## ############################################## # Expanded food web analyses # time steps for the ode function n.time.steps <- 2500 # adjust the number of replicates if necessary n.replicates <- 300 time.steps <- seq(0,n.time.steps,1) r_FT_min <- -0.01 r_FT_max <- -0.001 m_FT_min <- 0.0001 m_FT_max <- 0.0015 a_FT_PL_min <- 1e-6 a_FT_PL_max <- 1e-4 a_PL_FT_min <- -1e-2 a_PL_FT_max <- -1e-4 r_PL_0_min <- 0.01 r_PL_0_max <- 0.1 r_PL_nti_min <- 0.15 r_PL_nti_max <- 0.25 m_PL_min <- 0.001 m_PL_max <- 0.009 a_PL_DP_nti_min <- 1e-4 a_PL_DP_nti_max <- 1e-2 a_PL_DP_0_min <- 1e-7 a_PL_DP_0_max <- 1e-5 a_DP_PL_nti_min <- -1e-2 a_DP_PL_nti_max <- -1e-4 a_DP_PL_0_min <- -1e-5 a_DP_PL_0_max <- -1e-7 r_DP_0_min <- 0.4 r_DP_0_max <- 0.6 r_DP_nti_min <- 0.65 r_DP_nti_max <- 0.75 m_DP_min <- 0.005 m_DP_max <- 0.015 r_SV_min <- 0.05 r_SV_max <- 0.15 m_SV_min <- 0.00005 m_SV_max <- 0.00015 r_HM_0_min <- 0.15 r_HM_0_max <- 0.25 r_HM_nti_min <- 0.45 r_HM_nti_max <- 0.55 m_HM_0_min <- 0.0005 m_HM_0_max <- 0.0015 m_HM_nti_min <- 0.00005 m_HM_nti_max <- 0.00015 r_PLe_0_min <- 0.05 r_PLe_0_max <- 0.15 r_PLe_nti_min <- 0.15 r_PLe_nti_max <- 0.25 m_PLe_min <- 0.0002 m_PLe_max <- 0.0004 r_CM_0_min <- 0.15 r_CM_0_max <- 0.25 r_CM_nti_min <- 0.35 r_CM_nti_max <- 0.45 m_CM_min <- 0.00005 m_CM_max <- 0.00015 ################### # PL N_PL_0 <- 2189 # Pérez-Mellado et al. (2006) # DP N_DP_0 <- 150 # SV N_SV_0 <- 3000 # HM N_HM_0 <- 7187 # PLe N_PLe_0 <- 187 # CM N_CM_0 <- 10000 #### #### # rational function (Kéfi et al. 2012) # generalized for X = r, m, a nti.fun <- function(X_nti,X_0,N_0,...){ (X_nti * sum(...) + X_0 * N_0)/(sum(...) + N_0) } # dynamic equations model EFW.model <- function(time.steps, initial.abundances, param.list){ with(as.list(c(initial.abundances, param.list)),{ d_FT <- (r_FT * N_FT) - (m_FT * N_FT^2) + (a_FT_PL * N_FT * N_PL) d_PL <- (nti.fun(r_PL_nti,r_PL_0,sum(N_HM_0,N_PLe_0,N_CM_0),N_HM,N_PLe,N_CM) * N_PL) - (m_PL * N_PL^2) + ( a_PL_FT * N_PL * N_FT + nti.fun(a_PL_DP_nti,a_PL_DP_0,N_HM_0,N_HM) * N_PL * N_DP) d_DP <- (nti.fun(r_DP_nti,r_DP_0,sum(N_HM_0),N_HM) * N_DP) - (m_DP * N_DP^2) + (nti.fun(a_DP_PL_nti,a_DP_PL_0,N_HM_0,N_HM) * N_DP * N_PL) d_SV <- (r_SV * N_SV) - (m_SV * N_SV^2) d_CM <- (nti.fun(r_CM_nti,r_CM_0,N_PL_0,N_PL) * N_CM) - (m_CM * N_CM^2) d_PLe <- (nti.fun(r_PLe_nti,r_PLe_0,N_PL_0,N_PL) * N_PLe) - (m_PLe * N_PLe^2) d_HM <- (nti.fun(r_HM_nti,r_HM_0,sum(N_DP_0,N_PL_0),N_DP,N_PL) * N_HM) - (nti.fun(m_HM_nti,m_HM_0,N_SV_0,N_SV) * N_HM^2) list(c(d_FT, d_PL, d_DP, d_SV, d_CM, d_PLe, d_HM)) }) } initial.abundances <- c(N_FT = 8, N_PL = 5000, N_DP = 200, N_SV = 5000, N_CM = 5000, N_PLe = 200, N_HM = 5000) # results dataframe EFW.abundances <- data.frame(replicate = 1:n.replicates,FT = 0,PL = 0,DP = 0,SV = 0,CM = 0,PLe = 0,HM = 0) # run the model for(i.replicate in 1:n.replicates){ r_FT <- runif(1,r_FT_min,r_FT_max) m_FT <- runif(1,m_FT_min,m_FT_max) a_FT_PL <- runif(1,a_FT_PL_min,a_FT_PL_max) a_PL_FT <- runif(1,a_PL_FT_min,a_PL_FT_max) r_PL_0 <- runif(1,r_PL_0_min,r_PL_0_max) r_PL_nti <- runif(1,r_PL_nti_min,r_PL_nti_max) m_PL <- runif(1,m_PL_min,m_PL_max) a_PL_DP_nti <- runif(1,a_PL_DP_nti_min,a_PL_DP_nti_max) a_PL_DP_0 <- runif(1,a_PL_DP_0_min,a_PL_DP_0_max) a_DP_PL_nti <- runif(1,a_DP_PL_nti_min,a_DP_PL_nti_max) a_DP_PL_0 <- runif(1,a_DP_PL_0_min,a_DP_PL_0_max) r_DP_0 <- runif(1,r_DP_0_min,r_DP_0_max) r_DP_nti <- runif(1,r_DP_nti_min,r_DP_nti_max) m_DP <- runif(1,m_DP_min,m_DP_max) r_SV <- runif(1,r_SV_min,r_SV_max) m_SV <- runif(1,m_SV_min,m_SV_max) r_HM_0 <- runif(1,r_HM_0_min,r_HM_0_max) m_HM_0 <- runif(1,m_HM_0_min,m_HM_0_max) r_HM_nti <- runif(1,r_HM_nti_min,r_HM_nti_max) m_HM_nti <- runif(1,m_HM_nti_min,m_HM_nti_max) r_PLe_0 <- runif(1,r_PLe_0_min,r_PLe_0_max) r_PLe_nti <- runif(1,r_PLe_nti_min,r_PLe_nti_max) m_PLe <- runif(1,m_PLe_min,m_PLe_max) r_CM_0 <- runif(1,r_CM_0_min,r_CM_0_max) r_CM_nti <- runif(1,r_CM_nti_min,r_CM_nti_max) m_CM <- runif(1,m_CM_min,m_CM_max) param.list <- c(r_FT, m_FT, a_FT_PL, a_PL_FT, r_PL_0, r_PL_nti, m_PL, N_PL_0, a_PL_DP_nti, a_PL_DP_0, a_DP_PL_nti, a_DP_PL_0, r_DP_0, r_DP_nti, m_DP, N_DP_0, r_SV, m_SV, N_SV_0, r_HM_0, r_HM_nti, m_HM_0, m_HM_nti, N_HM_0, r_PLe_0, r_PLe_nti, m_PLe, N_PLe_0, r_CM_0, r_CM_nti, m_CM, N_CM_0) EFW.dynamics <- ode(y = initial.abundances,times = time.steps,func = EFW.model,parms = param.list) EFW.abundances[i.replicate,2:8] <- EFW.dynamics[nrow(EFW.dynamics),!dimnames(EFW.dynamics)[[2]] %in% c("time")] } # for i.replicate # most simulations are within reasonable abundance values # PL, DP and FT are the ones more prone to overshoot or extinction tt <- which(EFW.abundances$PL < 10000 & EFW.abundances$PL > 1e-4 & EFW.abundances$DP < 10000 & EFW.abundances$DP > 1e-4 & EFW.abundances$FT < 10000 & EFW.abundances$FT > 1e-4) EFW.valid.abundances <- EFW.abundances[tt,] # store the valid results write.table(x = EFW.valid.abundances,file = paste(wd,"/Aire_EFW_abundances_NTI.csv",sep=""),append = F,sep = ";",dec = ".") ################ ################ # Fig. 2 of main text, groups with NTIs # EFW.valid.gather <- gather(EFW.valid.abundances,key = "species", value = "abundance", -replicate) # # abundances.boxplot <- ggplot(EFW.valid.gather) # abundances.boxplot <- abundances.boxplot + geom_boxplot(aes(x = species, y = abundance)) # x11() # command for opening a plotting window in Linux # abundances.boxplot # # # abundance time series for a single replicate # # in this case, the last one stored at EFW.dynamics # time.series <- gather(as.data.frame(EFW.dynamics),key = "species", value = "abundance",-time) # time.series.plot <- ggplot(time.series) + geom_line(aes(x = time,y = abundance,color = species)) + ylim(0,max(EFW.dynamics[nrow(EFW.dynamics),!dimnames(EFW.dynamics)[[2]] %in% c("time")])) # x11() # time.series.plot ############################ ############################ # EFW without NTIs n.time.steps <- 2500 # adjust the number of replicates if necessary n.replicates <- 300 time.steps <- seq(0,n.time.steps,1) r_FT_min <- -0.01 r_FT_max <- -0.001 m_FT_min <- 0.0001 m_FT_max <- 0.0015 a_FT_PL_min <- 1e-6 a_FT_PL_max <- 1e-4 r_PL_min <- 0.01 r_PL_max <- 0.1 m_PL_min <- 0.001 m_PL_max <- 0.009 a_PL_FT_min <- -1e-2 a_PL_FT_max <- -1e-4 a_PL_DP_min <- 1e-4 a_PL_DP_max <- 1e-2 a_PL_HM_min <- 1e-6 a_PL_HM_max <- 1e-4 a_PL_CM_min <- 1e-7 a_PL_CM_max <- 1e-5 a_PL_PLe_min <- 1e-7 a_PL_PLe_max <- 1e-5 r_DP_min <- 0.4 r_DP_max <- 0.6 m_DP_min <- 0.005 m_DP_max <- 0.015 a_DP_PL_min <- -1e-5 a_DP_PL_max <- -1e-7 a_DP_HM_min <- 1e-6 a_DP_HM_max <- 1e-4 r_SV_min <- 0.05 r_SV_max <- 0.15 m_SV_min <- 0.00005 m_SV_max <- 0.00015 r_HM_min <- 0.15 r_HM_max <- 0.25 m_HM_min <- 0.0005 m_HM_max <- 0.0015 a_HM_DP_min <- -1e-6 a_HM_DP_max <- -1e-8 a_HM_PL_min <- -1e-6 a_HM_PL_max <- -1e-8 r_PLe_min <- 0.05 r_PLe_max <- 0.15 m_PLe_min <- 0.0002 m_PLe_max <- 0.0004 a_PLe_PL_min <- -1e-6 a_PLe_PL_max <- -1e-8 r_CM_min <- 0.15 r_CM_max <- 0.25 m_CM_min <- 0.00005 m_CM_max <- 0.00015 a_CM_PL_min <- -1e-6 a_CM_PL_max <- -1e-8 #### #### # dynamic equations model EFW.model.2 <- function(time.steps, initial.abundances, param.list){ with(as.list(c(initial.abundances, param.list)),{ d_FT <- (r_FT * N_FT) - (m_FT * N_FT^2) + (a_FT_PL * N_FT * N_PL) d_PL <- (r_PL * N_PL) - (m_PL * N_PL^2) + ((a_PL_FT * N_PL * N_FT) + (a_PL_DP * N_PL * N_DP) + (a_PL_HM * N_PL * N_HM) + (a_PL_PLe * N_PL * N_PLe) + (a_PL_CM * N_PL * N_CM)) d_DP <- (r_DP * N_DP) - (m_DP * N_DP^2) + ((a_DP_PL * N_DP * N_PL) + (a_DP_HM * N_DP * N_HM)) d_SV <- (r_SV * N_SV) - (m_SV * N_SV^2) d_CM <- (r_CM * N_CM) - (m_CM * N_CM^2) + (a_CM_PL * N_CM * N_PL) d_PLe <- (r_PLe * N_PLe) - (m_PLe * N_PLe^2) + (a_PLe_PL * N_PLe * N_PL) d_HM <- (r_HM * N_HM) - (m_HM * N_HM^2) + ((a_HM_PL * N_HM * N_PL) + (a_HM_DP * N_HM * N_DP)) list(c(d_FT, d_PL, d_DP, d_SV, d_CM, d_PLe, d_HM)) }) } initial.abundances <- c(N_FT = 8, N_PL = 500, N_DP = 200, N_SV = 5000, N_CM = 5000, N_PLe = 200, N_HM = 5000) EFW.abundances <- data.frame(replicate = 1:n.replicates,FT = 0,PL = 0,DP = 0,SV = 0,CM = 0,PLe = 0,HM = 0) for(i.replicate in 1:n.replicates){ #FT r_FT <- runif(1,r_FT_min,r_FT_max) m_FT <- runif(1,m_FT_min,m_FT_max) a_FT_PL <- runif(1,a_FT_PL_min,a_FT_PL_max) #PL r_PL <- runif(1,r_PL_min,r_PL_max) m_PL <- runif(1,m_PL_min,m_PL_max) a_PL_FT <- runif(1,a_PL_FT_min,a_PL_FT_max) a_PL_DP <- runif(1,a_PL_DP_min,a_PL_DP_max) a_PL_HM <- runif(1,a_PL_HM_min,a_PL_HM_max) a_PL_PLe <- runif(1,a_PL_PLe_min,a_PL_PLe_max) a_PL_CM <- runif(1,a_PL_CM_min,a_PL_CM_max) #DP r_DP <- runif(1,r_DP_min,r_DP_max) m_DP <- runif(1,m_DP_min,m_DP_max) a_DP_PL <- runif(1,a_DP_PL_min,a_DP_PL_max) a_DP_HM <- runif(1,a_DP_HM_min,a_DP_HM_max) #SV r_SV <- runif(1,r_SV_min,r_SV_max) m_SV <- runif(1,m_SV_min,m_SV_max) #HM r_HM <- runif(1,r_HM_min,r_HM_max) m_HM <- runif(1,m_HM_min,m_HM_max) a_HM_PL <- runif(1,a_HM_PL_min,a_HM_PL_max) a_HM_DP <- runif(1,a_HM_DP_min,a_HM_DP_max) #PLe r_PLe <- runif(1,r_PLe_min,r_PLe_max) m_PLe <- runif(1,m_PLe_min,m_PLe_max) a_PLe_PL <- runif(1,a_PLe_PL_min,a_PLe_PL_max) #CM r_CM <- runif(1,r_CM_min,r_CM_max) m_CM <- runif(1,m_CM_min,m_CM_max) a_CM_PL <- runif(1,a_CM_PL_min,a_CM_PL_max) param.list <- c(r_FT, m_FT, a_FT_PL, r_PL, m_PL, a_PL_FT, a_PL_DP, a_PL_HM, a_PL_PLe, a_PL_CM, r_DP, m_DP, a_DP_PL, a_DP_HM, r_SV, m_SV, r_HM, m_HM, a_HM_PL, a_HM_DP, r_PLe, m_PLe, a_PLe_PL, r_CM, m_CM, a_CM_PL) EFW.dynamics <- ode(y = initial.abundances,times = time.steps,func = EFW.model.2,parms = param.list) EFW.abundances[i.replicate,2:8] <- EFW.dynamics[nrow(EFW.dynamics),!dimnames(EFW.dynamics)[[2]] %in% c("time")] } # for i.replicate tt <- which(EFW.abundances$PL < 10000 & EFW.abundances$PL > 1e-4 & EFW.abundances$DP < 10000 & EFW.abundances$DP > 0 & EFW.abundances$FT < 10000 & EFW.abundances$FT > 1e-4) EFW.valid.abundances <- EFW.abundances[tt,] # store the valid results write.table(x = EFW.valid.abundances,file = paste(wd,"/Aire_EFW_abundances_NO_NTI.csv",sep=""),append = F,sep = ";",dec = ".") ############################ ############################ # FIG 2 of main text (before editing) abundances.nti <- read.table(file = paste(wd,"Aire_EFW_abundances_NTI.csv",sep=""),header = T,sep = ";",dec = ".") abundances.NOnti <- read.table(file = paste(wd,"Aire_EFW_abundances_NO_NTI.csv",sep=""),header = T,sep = ";",dec = ".") abundances.nti$interactions <- "trophic + non trophic" abundances.NOnti$interactions <- "trophic" all.abundances <- rbind(abundances.NOnti,abundances.nti) all.abundances <- gather(all.abundances,key = "species", value = "abundance", -c(replicate,interactions)) ############## ############## abundances.boxplot <- ggplot(all.abundances) abundances.boxplot <- abundances.boxplot + geom_boxplot(aes(x = species, y = abundance,fill = interactions)) abundances.boxplot <- abundances.boxplot + scale_fill_manual(values = c("#969696","#636363")) x11() abundances.boxplot ############## ############## # associated wilcoxon tests CM.EFW.wilcox <- wilcox.test(abundances.nti$CM, abundances.NOnti$CM) DP.EFW.wilcox <- wilcox.test(abundances.nti$DP, abundances.NOnti$DP) FT.EFW.wilcox <- wilcox.test(abundances.nti$FT, abundances.NOnti$FT) HM.EFW.wilcox <- wilcox.test(abundances.nti$HM, abundances.NOnti$HM) PL.EFW.wilcox <- wilcox.test(abundances.nti$PL, abundances.NOnti$PL) PLe.EFW.wilcox <- wilcox.test(abundances.nti$PLe, abundances.NOnti$PLe) SV.EFW.wilcox <- wilcox.test(abundances.nti$SV, abundances.NOnti$SV) ############################ ############################ # Equal footing analyses n.time.steps <- 2500 # set to 100 for testing purposes. Actual number was 100000 n.replicates <- 100 # build dynamic system # FT r_FT_min <- -0.005 r_FT_max <- -0.0005 b_FT_min <- 1e-5 b_FT_max <- 1e-4 c_FT_min <- 1e-3 c_FT_max <- 1e-3 # PL r_PL_min <- 0.01 r_PL_max <- 0.1 b_PL_min <- 1e-5 b_PL_max <- 1e-4 c_PL_min <- 1e-3 c_PL_max <- 1e-3 # DP r_DP_min <- 0.4 r_DP_max <- 0.6 b_DP_min <- 1e-5 b_DP_max <- 1e-4 c_DP_min <- 1e-3 c_DP_max <- 1e-3 # SV r_SV_min <- 0.05 r_SV_max <- 0.15 b_SV_min <- 1e-5 b_SV_max <- 1e-4 c_SV_min <- 1e-3 c_SV_max <- 1e-3 # HM r_HM_min <- 0.15 r_HM_max <- 0.25 b_HM_min <- 1e-5 b_HM_max <- 1e-4 c_HM_min <- 1e-3 c_HM_max <- 1e-3 # PLe r_PLe_min <- 0.05 r_PLe_max <- 0.15 b_PLe_min <- 1e-5 b_PLe_max <- 1e-4 c_PLe_min <- 1e-3 c_PLe_max <- 1e-3 # CM r_CM_min <- 0.15 r_CM_max <- 0.25 b_CM_min <- 1e-5 b_CM_max <- 1e-4 c_CM_min <- 1e-3 c_CM_max <- 1e-3 ### interaction coefficients, a_A_B = on A from B # interaction type can be either weak,strongANT (strong antagonism) or strongMUT (strong mutualism) # run the model again for each of the three types and store the results interaction.type <- "strongANT" # interaction.type <- "strongMUT" # interaction.type <- "weak" # FT-PL a_FT_PL_min <- ifelse(interaction.type == "strongANT",1e-4,1e-7) a_FT_PL_max <- ifelse(interaction.type == "strongANT",1e-2,1e-5) # PL-FT a_PL_FT_min <- ifelse(interaction.type == "strongANT",-1e-2,-1e-5) a_PL_FT_max <- ifelse(interaction.type == "strongANT",-1e-4,-1e-7) # PL-DP a_PL_DP_min <- ifelse(interaction.type == "strongANT",1e-4,1e-7) a_PL_DP_max <- ifelse(interaction.type == "strongANT",1e-2,1e-5) # PL-HM a_PL_HM_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_PL_HM_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # PL-CM a_PL_CM_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_PL_CM_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # PL-PLe a_PL_PLe_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_PL_PLe_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # DP-PL a_DP_PL_min <- ifelse(interaction.type == "strongANT",-1e-2,-1e-5) a_DP_PL_max <- ifelse(interaction.type == "strongANT",-1e-4,-1e-7) # DP-HM a_DP_HM_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_DP_HM_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # HM-DP a_HM_DP_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_HM_DP_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # HM-PL a_HM_PL_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_HM_PL_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # HM-SV a_HM_SV_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_HM_SV_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # PLe-PL a_PLe_PL_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_PLe_PL_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # CM-PL a_CM_PL_min <- ifelse(interaction.type == "strongMUT",1e-4,1e-7) a_CM_PL_max <- ifelse(interaction.type == "strongMUT",1e-2,1e-5) # dynamic equations model EF.model <- function(time.steps, initial.abundances, param.list){ abundances.data <- rep(0,length(initial.abundances)) ## calculate the overall effect of the interactions (some terms will be zero) SI <- numeric(length(initial.abundances)) for(i.sp in 1:length(initial.abundances)){ for(j.sp in 1:length(initial.abundances)){ # if(i.sp != j.sp){ actually it's zero SI[i.sp] <- SI[i.sp] + param.list$adjacency.matrix[i.sp,j.sp]*initial.abundances[j.sp] # } } } # self-regulating term SL <- numeric(length(initial.abundances)) for(i.sp in 1:length(initial.abundances)){ for(j.sp in 1:length(initial.abundances)){ SL[i.sp] <- SL[i.sp] + param.list$c[i.sp]*param.list$adjacency.matrix[i.sp,j.sp]*initial.abundances[j.sp] } SL[i.sp] <- (SL[i.sp] + param.list$b[i.sp]) * initial.abundances[i.sp] } ## equations for(i.sp in 1:length(initial.abundances)){ abundances.data[i.sp] <- (param.list$r0[i.sp] + SI[i.sp] - SL[i.sp]) * initial.abundances[i.sp] # abundances.data[i.sp] <- ifelse(abundances.data[i.sp] > param.list$extinction.threshold,abundances.data[i.sp],0) } return(list(abundances.data)) } EF.abundances <- data.frame(replicate = 1:n.replicates, leading.eigenvalue = 0, average.antagonist.strength = 0, average.mutualist.strength = 0, FT = 0,PL = 0,DP = 0,SV = 0,CM = 0,PLe = 0,HM = 0) for(i.replicate in 1:n.replicates){ # FT r_FT <- runif(1,r_FT_min,r_FT_max) b_FT <- runif(1,b_FT_min,b_FT_max) c_FT <- runif(1,c_FT_min,c_FT_max) # PL r_PL <- runif(1,r_PL_min,r_PL_max) b_PL <- runif(1,b_PL_min,b_PL_max) c_PL <- runif(1,c_PL_min,c_PL_max) # DP r_DP <- runif(1,r_DP_min,r_DP_max) b_DP <- runif(1,b_DP_min,b_DP_max) c_DP <- runif(1,c_DP_min,c_DP_max) # SV r_SV <- runif(1,r_SV_min,r_SV_max) b_SV <- runif(1,b_SV_min,b_SV_max) c_SV <- runif(1,c_SV_min,c_SV_max) # HM r_HM <- runif(1,r_HM_min,r_HM_max) b_HM <- runif(1,b_HM_min,b_HM_max) c_HM <- runif(1,c_HM_min,c_HM_max) # PLe r_PLe <- runif(1,r_PLe_min,r_PLe_max) b_PLe <- runif(1,b_PLe_min,b_PLe_max) c_PLe <- runif(1,c_PLe_min,c_PLe_max) # CM r_CM <- runif(1,r_CM_min,r_CM_max) b_CM <- runif(1,b_CM_min,b_CM_max) c_CM <- runif(1,c_CM_min,c_CM_max) ########## # adjacency matrix a_FT_PL <- runif(1,a_FT_PL_min,a_FT_PL_max) a_PL_FT <- runif(1,a_PL_FT_min,a_PL_FT_max) a_PL_DP <- runif(1,a_PL_DP_min,a_PL_DP_max) a_PL_HM <- runif(1,a_PL_HM_min,a_PL_HM_max) a_PL_PLe <- runif(1,a_PL_PLe_min,a_PL_PLe_max) a_PL_CM <- runif(1,a_PL_CM_min,a_PL_CM_max) a_DP_PL <- runif(1,a_DP_PL_min,a_DP_PL_max) a_DP_HM <- runif(1,a_DP_HM_min,a_DP_HM_max) a_HM_PL <- runif(1,a_HM_PL_min,a_HM_PL_max) a_HM_DP <- runif(1,a_HM_DP_min,a_HM_DP_max) a_HM_SV <- runif(1,a_HM_SV_min,a_HM_SV_max) a_PLe_PL <- runif(1,a_PLe_PL_min,a_PLe_PL_max) a_CM_PL <- runif(1,a_CM_PL_min,a_CM_PL_max) adjacency.matrix <- matrix(c(0,a_FT_PL,0,0,0,0,0, a_PL_FT,0,a_PL_DP,0,a_PL_HM,a_PL_PLe,a_PL_CM, 0,a_DP_PL,0,0,a_DP_HM,0,0, 0,0,0,0,0,0,0, 0,a_HM_PL,a_HM_DP,a_HM_SV,0,0,0, 0,a_PLe_PL,0,0,0,0,0, 0,a_CM_PL,0,0,0,0,0),nrow = 7,byrow = T) EF.initial.abundances <- c(N_FT = 8, N_PL = 5000, N_DP = 200, N_SV = 5000, N_CM = 5000, N_PLe = 200, N_HM = 5000) ###### EF.param.list <- list() EF.param.list$r0 <- c(r_FT,r_PL,r_DP,r_SV,r_HM,r_PLe,r_CM) EF.param.list$b <- c(b_FT,b_PL,b_DP,b_SV,b_HM,b_PLe,b_CM) EF.param.list$c <- c(c_FT,c_PL,c_DP,c_SV,c_HM,c_PLe,c_CM) EF.param.list$adjacency.matrix <- adjacency.matrix EF.dynamics <- ode(y = EF.initial.abundances,times = time.steps,func = EF.model,parms = EF.param.list) EF.abundances[i.replicate,c("FT","PL","DP","SV","CM","PLe","HM")] <- EF.dynamics[nrow(EF.dynamics),!dimnames(EF.dynamics)[[2]] %in% c("time")] community.matrix <- jacobian.full(y = EF.dynamics[nrow(EF.dynamics),!dimnames(EF.dynamics)[[2]] %in% c("time")],func = EF.model,parms = EF.param.list) community.matrix[is.na(community.matrix)] <- 0 ### EF.abundances[i.replicate,"average.antagonist.strength"] <- mean(abs(c(a_FT_PL,a_PL_FT,a_PL_DP,a_DP_PL))) EF.abundances[i.replicate,"average.mutualist.strength"] <- mean(a_PL_HM,a_PL_PLe,a_PL_CM,a_DP_HM,a_HM_PL,a_HM_SV,a_HM_DP,a_PLe_PL,a_CM_PL) EF.abundances[i.replicate,"leading.eigenvalue"] <- max(Re(eigen(community.matrix)$values)) } EF.abundances[is.na(EF.abundances)] <- 0 # store the results write.table(x = EF.abundances,file = paste(wd,"Aire_EF_abundances_",interaction.type,".csv",sep=""),append = F,sep = ";",dec = ".") ################# ################# # after running lines 595-780 for the three types of interactions, three files should be in the working directory # These files should be small enough in order to append them to a single file in the script abundances.EF.weak <- read.table(file = paste(wd,"Aire_EF_abundances_weak.csv",sep=""),header = T,sep = ";",dec = ".") abundances.EF.strongANT <- read.table(file = paste(wd,"Aire_EF_abundances_strongANT.csv",sep=""),header = T,sep = ";",dec = ".") abundances.EF.strongMUT <- read.table(file = paste(wd,"Aire_EF_abundances_strongMUT.csv",sep=""),header = T,sep = ";",dec = ".") abundances.EF.weak$interactions <- "weak interactions" abundances.EF.strongMUT$interactions <- "stronger mutualism" abundances.EF.strongANT$interactions <- "stronger antagonism" all.abundances <- rbind(abundances.EF.weak,abundances.EF.strongMUT,abundances.EF.strongANT) ################# ################# # some exploratory results: #how many stable/unstable communities? all.abundances %>% group_by(interactions) %>% summarize(stable = sum(leading.eigenvalue < 0), unstable = sum(leading.eigenvalue >= 0)) # how many with all species present all.abundances %>% group_by(interactions) %>% summarize(all.sp = sum(FT > 0 & PL > 0 & DP > 0 & SV > 0 & HM > 0 & PLe > 0 & CM > 0)) # how many with all species within reasonable limits all.abundances %>% group_by(interactions) %>% summarize(all.sp = sum(FT > 0 & PL > 0 & DP > 0 & SV > 0 & HM > 0 & PLe > 0 & CM > 0 & FT < 1e5 & PL < 1e5 & DP < 1e5 & SV < 1e5 & HM < 1e5 & PLe < 1e5 & CM < 1e5)) # unstable communities, # how many with all species present all.abundances %>% group_by(interactions) %>% filter(leading.eigenvalue >= 0) %>% summarize(all.sp = sum(FT > 0 & PL > 0 & DP > 0 & SV > 0 & HM > 0 & PLe > 0 & CM > 0)) # how many with all species within reasonable limits all.abundances %>% group_by(interactions) %>% filter(leading.eigenvalue >= 0) %>% summarize(all.sp = sum(FT > 0 & PL > 0 & DP > 0 & SV > 0 & HM > 0 & PLe > 0 & CM > 0 & FT < 1e5 & PL < 1e5 & DP < 1e5 & SV < 1e5 & HM < 1e5 & PLe < 1e5 & CM < 1e5)) # set of stable communities stable.communities <- subset(all.abundances, leading.eigenvalue <0) table(stable.communities$interactions) ################# ################# # Fig. 3 of main text all.abundances.cutoff <- all.abundances max.eigen <- 0.06 all.abundances.cutoff <- subset(all.abundances.cutoff, leading.eigenvalue < max.eigen) all.abundances.cutoff$color <- ifelse(all.abundances.cutoff$leading.eigenvalue >= 0,"#252525","#636363") eigen.3d <- scatterplot3d(x = log(all.abundances.cutoff$average.mutualist.strength), y = log(all.abundances.cutoff$average.antagonist.strength), z = all.abundances.cutoff$leading.eigenvalue, color = all.abundances.cutoff$color, # highlight.3d = T, # type = "h", cex.symbols = .7, col.axis = "#969696", scale.y = .5, # lab = c(4,4,2), xlim = c(-17,-3), ylim = c(-16,-4.5), zlim = c(-0.02,0.02), grid = F, xlab = "log(mutualism strength)", ylab = "log(antagonism strength)", zlab = "leading eigenvalue", angle = 35) eigen.3d$plane3d(0,0,0,col="#636363",lty="dashed") # paint again the black dots for better visualization eigen.3d$points3d(log(all.abundances.cutoff$average.mutualist.strength[all.abundances.cutoff$leading.eigenvalue > 0]), log(all.abundances.cutoff$average.antagonist.strength[all.abundances.cutoff$leading.eigenvalue > 0]), all.abundances.cutoff$leading.eigenvalue[all.abundances.cutoff$leading.eigenvalue > 0], col="#252525") ####################### ####################### # Supplementary material Fig. S7, i.e. records with leading.eigenvalue >> 0 all.abund.greater <- subset(all.abundances, leading.eigenvalue > max.eigen & leading.eigenvalue < 1e100) eigen2.3d <- scatterplot3d(x = log(all.abund.greater$average.mutualist.strength), y = log(all.abund.greater$average.antagonist.strength), z = log(all.abund.greater$leading.eigenvalue), # color = all.abundances.cutoff$color, # highlight.3d = T, # type = "h", cex.symbols = .7, col.axis = "#969696", scale.y = .5, # lab = c(4,4,2), xlim = c(-17,-3), ylim = c(-16,-4.5), # zlim = c(-0.02,0.06), grid = F, xlab = "log(mutualism strength)", ylab = "log(antagonism strength)", zlab = "log(leading eigenvalue)", angle = 35) ######################## ######################## # distribution of interaction strengths for the 4 groups (lower panels of Fig. 3) all.abund <- gather(all.abundances,key = "species", value = "abundance", -c(replicate,leading.eigenvalue,average.antagonist.strength,average.mutualist.strength,interactions)) all.abund <- gather(all.abund,key = "interaction.type", value = "strength", -c(replicate,leading.eigenvalue,interactions,species,abundance)) all.abund$interactions <- factor(all.abund$interactions, levels = c("weak interactions", "stronger mutualism", "stronger antagonism")) levels(all.abund$interactions)[levels(all.abund$interactions)=="weak interactions"] <- "a)" levels(all.abund$interactions)[levels(all.abund$interactions)=="stronger mutualism"] <- "b)" levels(all.abund$interactions)[levels(all.abund$interactions)=="stronger antagonism"] <- "c)" line.size <- 0.5 strength.density <- ggplot(all.abund) strength.density <- strength.density + geom_density(aes(x = log(strength), linetype = interaction.type), size = line.size, alpha = .3) strength.density <- strength.density + facet_grid(.~interactions) strength.density <- strength.density + scale_linetype_discrete(guide=FALSE) x11() strength.density ################## ################## # Supplementary material Fig. S6 abundances.EF.weak <- read.table(file = paste(wd,"Aire_EF_abundances_weak.csv",sep=""),header = T,sep = ";",dec = ".") abundances.EF.strongANT <- read.table(file = paste(wd,"Aire_EF_abundances_strongANT.csv",sep=""),header = T,sep = ";",dec = ".") abundances.EF.strongMUT <- read.table(file = paste(wd,"Aire_EF_abundances_strongMUT.csv",sep=""),header = T,sep = ";",dec = ".") abundances.EF.weak$interactions <- "weak interactions" abundances.EF.strongMUT$interactions <- "stronger mutualism" abundances.EF.strongANT$interactions <- "stronger antagonism" all.abundances <- rbind(abundances.EF.weak,abundances.EF.strongMUT,abundances.EF.strongANT) all.abund <- gather(all.abundances,key = "species", value = "abundance", -c(replicate,leading.eigenvalue,average.antagonist.strength,average.mutualist.strength,interactions)) all.abund <- subset(all.abund, interactions != "stronger antagonism") all.abund <- droplevels(all.abund) abundances.boxplot <- ggplot(all.abund) abundances.boxplot <- abundances.boxplot + geom_boxplot(aes(x = species, y = abundance,fill = interactions)) abundances.boxplot <- abundances.boxplot + ylim(0,6000) x11() abundances.boxplot ################## ################## # associated wilcoxon tests abund.strongmut <- subset(all.abundances, interactions == "stronger mutualism") abund.weak <- subset(all.abundances, interactions == "weak interactions") CM.EF.wilcox <- wilcox.test(abund.strongmut$CM, abund.weak$CM) DP.EF.wilcox <- wilcox.test(abund.strongmut$DP, abund.weak$DP) FT.EF.wilcox <- wilcox.test(abund.strongmut$FT, abund.weak$FT) HM.EF.wilcox <- wilcox.test(abund.strongmut$HM, abund.weak$HM) PL.EF.wilcox <- wilcox.test(abund.strongmut$PL, abund.weak$PL) PLe.EF.wilcox <- wilcox.test(abund.strongmut$PLe, abund.weak$PLe) SV.EF.wilcox <- wilcox.test(abund.strongmut$SV, abund.weak$SV)