#################################################### # POPULATION DYNAMICS IN FLUCTUATING ENVIRONEMENTS # #################################################### ######## Working directory and workspace - if needed ###### # setwd("D:/DONNEESPROFESSIONNELLES100GOMAX/Manip Rescue in Flucuating environment/PDFE/Final_PDFE") # save.image("stage_stochastic.RData") #Save workspace # load("stage_stochasticR.RData") #Load previous workspace rm(list=ls()) #Clean workspace ######## PREMABULE ######## #LOAD PACKAGES library("data.table") library("ggplot2") library("cowplot") library("MASS") library("TMB") library("glmmTMB") library("lme4") library("moments") library("survival") library("survminer") library("purrr") library("epiR") library("car") library("ellipse") library("plyr") library("Cairo") library("gridExtra") library("mice") #FUNCTIONS LRtest <- function(full,restricted){ statistic <- 2*(restricted$value - full$value) df <- length(full$par) - length(restricted$par) p.value <- 1-pchisq(statistic,df=df) data.frame(statistic,df,p.value) } #Compute likelihood ratio test between 2 optimized model (from TMB & optim) AIC = function(model){ aic = 2 * length(model$par) + 2 * model$value return(aic) } #Compute the AIC of an optimized model (from TMB & optim) plotylinlog = function(Concentration){ cplot = ifelse(is.na(Concentration), NA, ifelse(Concentration < 1000, Concentration, 1000 * log(Concentration) / log(10) - 2000)) return(cplot) } #Plot population size < 1000 on a linear scale and population sizes > 1000 on a log scale. # Plot commands { #Ticks for the linear / log scale breackticks = c(0., 1000., 2000., 3000., 4000., 5000., 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000., 1301.03, 1477.12, 1602.06, 1698.97, 1778.15, 1845.1, 1903.09, 1954.24, 2000., 2301.03, 2477.12, 2602.06, 2698.97, 2778.15, 2845.1, 2903.09, 2954.24, 3000., 3301.03, 3477.12, 3602.06, 3698.97, 3778.15, 3845.1, 3903.09, 3954.24, 4000., 4301.03, 4477.12, 4602.06, 4698.97, 4778.15, 4845.1, 4903.09, 4954.24) labelticks = c(0., 1000., 10000., 100000., 1.*10^6, 1.*10^7, " ", " ", " ", " ", " \ ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " \ ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " \ ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ") } ############## READ AND CLEAN DATA############# dilution = 0.15;# Dilution level at each transfer #READ, ORDER POPULATION DYNAMICS DATASET AND TRANSFORM VARIABLES data = read.table("population_dynamics.csv", header = TRUE, sep=",") data = as.data.table(data, keep.rownames=FALSE) data$Date = as.Date(data$Date) data$Plate = as.factor(data$Plate) data$Serie = as.factor(data$Serie) data = data[with(data, order(Date,Plate, Well)),] data = data[, Day := as.numeric(data$Date - as.Date("2017-08-10"))] #Calculate the day data = data[, Transfer :=as.numeric(floor(2 * Day/ 7))] #Calculate the transfer data = data[, ID := paste(Genotype, Autocorrelation, Serie, Replicate, sep="_")] #Generate line ID data = data[, Volume := 0.00002 + round(TotalCount/TotalConcentration /0.00023) * 0.00023] #Calculate the volume of cytometer sample Daymax = max(data$Day) TransferMax = max(data$Transfer) #CORRECT FOR CYTOMER CROSS CONTAMINATION AND CALCULATE CONCENTRATION #Our cytometer is a plate reader. Carry over occur between wells and a preliminary experiment show that #it represents in average 0.06% of the cells in the previously readed well. PrevAliveCount = c(0,data$DunaAliveCount[1:length(data$data$DunaAliveCount)-1]) data = cbind(data, PrevAliveCount) #Concentration in the previous well data = data[, DunaAliveCorr := ifelse(is.na(DunaAliveCount), NA, ifelse(Well == 14, DunaAliveCount, ifelse(round(DunaAliveCount - 0.0006 * PrevAliveCount) > 0, round(DunaAliveCount - 0.0006 * PrevAliveCount), 0)))] data = data[ , DunaConc := ifelse(is.na(DunaAliveCorr), NA, DunaAliveCorr / Volume)] #DELETE OUTSIDERS (identified in a Mathematica11 File using raw data from cytometry and spectrometry) data[Plate == 3 & Date == "2017-09-25" & Well == 14, Fcorrfit := NA ] data[Plate == 3 & Date == "2017-09-25" & Well == 14, ODcorrfit := NA ] data[Plate == 5 & Date == "2017-11-02", ODcorrfit := NA ] data[Plate == 5 & Date == "2017-11-02", Fcorrfit := NA ] data[Plate == 5 & Date == "2017-11-02", DunaConc := NA ] data[Plate == 6 & Date == "2017-11-02", ODcorrfit := NA ] data[Plate == 6 & Date == "2017-11-02", Fcorrfit := NA ] data[Plate == 6 & Date == "2017-11-02", DunaConc := NA ] data[Plate == 7 & Date == "2017-10-16" & Well %in% c(18,30,42,54,66,78,79), ODcorrfit := NA ] data[Plate == 7 & Date == "2017-10-16" & Well %in% c(18,30,42,54,66,78,79), Fcorrfit := NA ] data[Plate == 7 & Date == "2017-10-16" & Well %in% c(18,30,42,54,66,78,79), DunaConc := NA ] data[Plate == 8 & Date == "2017-12-21", ODcorrfit := NA ] data[Plate == 8 & Date == "2017-12-21", Fcorrfit := NA ] data[Plate == 8 & Date == "2017-12-21", DunaConc := NA ] data[Plate == 8 & Date == "2017-09-25" & Well %in% c(32, 80), ODcorrfit := NA ] data[Plate == 8 & Date == "2017-09-25" & Well %in% c(32, 80), Fcorrfit := NA ] data[Plate == 8 & Date == "2017-09-25" & Well %in% c(32, 80), DunaConc := NA ] data[Plate == 8 & Date == "2017-09-07" & Well == 77, ODcorrfit := NA ] data[Plate == 8 & Date == "2017-09-07" & Well == 77, Fcorrfit := NA ] data[Plate == 8 & Date == "2017-09-07" & Well == 77, DunaConc := NA ] data[Plate == 12 & Date == "2017-09-04" & Well == 17, ODcorrfit := NA ] data[Plate == 12 & Date == "2017-09-04" & Well == 17, Fcorrfit := NA ] data[Plate == 12 & Date == "2017-09-04" & Well == 17, DunaConc := NA ] data[Plate == 13 & Date == "2017-09-14" & Well %in% c(20, 32, 44, 56, 68, 80), ODcorrfit := NA ] data[Plate == 13 & Date == "2017-09-14" & Well %in% c(20, 32, 44, 56, 68, 80), Fcorrfit := NA ] data[Plate == 13 & Date == "2017-09-14" & Well %in% c(20, 32, 44, 56, 68, 80), DunaConc := NA ] data[Plate == 14 & Date == "2017-09-28" & Well == 42, ODcorrfit := NA ] data[Plate == 14 & Date == "2017-09-28" & Well == 42, Fcorrfit := NA ] data[Plate == 14 & Date == "2017-09-28" & Well == 42, DunaConc := NA ] data[Plate == 16 & Date == "2017-10-12", ODcorrfit := NA ] data[Plate == 16 & Date == "2017-10-12", Fcorrfit := NA ] data[Plate == 16 & Date == "2017-10-12", DunaConc := NA ] data[Plate == 17 & Date == "2017-11-13", ODcorrfit := NA ] data[Plate == 17 & Date == "2017-11-13", Fcorrfit := NA ] data[Plate == 17 & Date == "2017-11-13", DunaConc := NA ] #DELETE CONTAMINATED TIME SERIES (After ITS sequences check) - # series that went extinct before a visible contamination event (Mathematica plot) were kept) data = data[!(Genotype == "B" & Autocorrelation == "cst" & Serie == 1 & Replicate == 1)] data = data[!(Genotype == "B" & Autocorrelation == "cst" & Serie == 2 & Replicate == 5 & Day > 70)] data = data[!(Genotype == "B" & Autocorrelation == "0.5" & Serie == 25 & Transfer > 5)] data = data[!(Genotype == "B" & Autocorrelation == "0.5" & Serie == 36 & Transfer > 2)] data = data[!(Genotype == "B" & Autocorrelation == "0.5" & Serie == 1 & Transfer > 2)] data = data[!(Genotype == "B" & Autocorrelation == "0.5" & Serie == 1 & Replicate == 2)] data = data[!(Genotype == "B" & Autocorrelation == "0.9" & Serie == 20 )] data = data[!(Genotype == "B" & Autocorrelation == "0.9" & Serie == 24 & Transfer > 9)] data = data[!(Genotype == "A" & Autocorrelation == "0" & Serie == 8)] data = data[!(Genotype == "A" & Autocorrelation == "0" & Serie == 36)] data = data[!(Genotype == "A" & Autocorrelation == "0.5" & Serie == 39)] data = data[!(Genotype == "C" & Autocorrelation == "0" & Serie == 14)] #FIT SPECTROMETERS / CYTOMETERS CONCENTRATION MEASURES { #CORRECT ABSORBANCE AND FLUORESCENCE MEASURES WITH BLANK #Find the minimal fluorescence value for each plate, at each transfer minFluo = c(rep(NA, length(data[Transfer %in% c(0),F685]))) # No empty wells measured (no extinction) in the first transfer for (transfer in 1:max(as.numeric(data[,Transfer]))){ for (plate in unique(as.numeric(data[Transfer == transfer,Plate]))){ minFluo =c(minFluo, rep(min(data[Plate == plate & Transfer == transfer, F685]),length(data[Transfer == transfer & Plate == plate,F685]))) } } minFluo = replace(minFluo, is.na(minFluo), mean(minFluo,na.rm = TRUE)) # If no value provided, replace by the average blank value data = cbind(data,minFluo) data = data[ , Fcorr := ifelse(F685 - minFluo > 0, F685 - minFluo, 0)] #Find the minimal absorbance value for each plate, at each transfer (sale process as for Fluo) minOD = c(rep(NA, length(data[Transfer %in% c(0),OD680]))) for (transfer in 1:max(as.numeric(data[,Transfer]))){ for (plate in unique(as.numeric(data[Transfer == transfer,Plate]))){ minOD =c(minOD, rep(min(data[Plate == plate & Transfer == transfer, OD680]),length(data[Transfer == transfer & Plate == plate,OD680]))) } } minOD = replace(minOD, is.na(minOD), mean(minOD,na.rm = TRUE)) data = cbind(data,minOD) data = data[ , ODcorr := ifelse(OD680 - minOD > 0, OD680 - minOD, 0)] # FLUORESCENCE REGRESSION regfluo = lm(DunaConc ~ 0 + Fcorr, data = data[F685 !="NA" & DunaAliveCorr !="NA"]) summary(regfluo) coeff_Fluo_hyp = summary(regfluo)$coefficients[1]; par(mfrow=c(2,2)) plot(regfluo) ggplot() + geom_point(data = data[Fcorr!="NA" & DunaAliveCorr !="NA"], aes(x = Fcorr, y = DunaConc, color = as.numeric(Day))) + geom_line(data = data[Fcorr!="NA" & DunaAliveCorr != "NA"], aes(x = Fcorr, y = Fcorr * summary(regfluo)$coefficients[1])) + scale_y_continuous(trans='log10') + scale_x_continuous(trans='log10') + scale_colour_gradient2(low = rgb(0.8,0.8,0.8), high = "black")+ guides(colour = FALSE) + labs( x = "Fluorescence", y = "Population density (cytometer estimation)") + geom_vline(xintercept = 400) data = data[, Fcorrfit := Fcorr * coeff_Fluo_hyp] # ABSORBANCE REGRESSION regOD = lm(DunaConc ~ 0 + ODcorr , data = data[OD680 !="NA" & DunaAliveCorr > 0]) summary(regOD) coeff_OD_hyp = summary(regOD)$coefficients[1]; par(mfrow=c(2,2)) plot(regOD) ggplot() + geom_point(data = data[OD680!="NA" & DunaAliveCorr != "NA"], aes(x = ODcorr, y = DunaConc, color = factor(Transfer))) + geom_line(data = data[OD680!="NA" & DunaAliveCorr != "NA"], aes(x = ODcorr, y = ODcorr * summary(regOD)$coefficients[1] # + OD680 * summary(regOD)$coefficients[2] # + OD2 * summary(regOD)$coefficients[3] )) + scale_y_continuous(trans='log') + scale_x_continuous(trans='log') + scale_colour_hue(l=60) + guides(colour = FALSE) + labs( x = "OD 680", y = "Population density (cytometer estimation)") + geom_vline(xintercept = 0.01) data = data[, ODcorrfit := ODcorr * coeff_OD_hyp] } # Delete empty wells used for contamination data = data[Autocorrelation != conta] #Compute a first estimation of population size based on fluorescence and cytometer counts - used # for survival analysis. data = data[ , ConcApprox := colMeans(rbind(Fcorrfit, DunaConc), na.rm=TRUE)] #DATA GLYCEROL data_glycerol = read.table("data_cinetique.csv", header = TRUE, sep=",") data_glycerol = as.data.table(data_glycerol) data_glycerol$ID0 = paste(data_glycerol$Traitement, data_glycerol$Temps, sep = "_") # delete bubbles data_glycerol[ID0 == "33_135" & Replicate == 1 & total_extra == "extra"]$GlycerolEstimate = NA data_glycerol[ID0 == "33_135" & Replicate == 2 & total_extra == "total"]$GlycerolEstimate = NA data_glycerol[ID0 == "31_45" & Replicate == 2 & total_extra == "extra"]$GlycerolEstimate = NA ######## SURVIVAL ANALYSIS ######## # A population is considered dead if it remains for 5 transfers under 1000 individuals/mL Deadcontrol = 1000 Tdeadcontrol = 5; #Split data according to lines ID data = data[with(data, order(ID)),] data_split = split(data$ConcApprox,data$ID) time_split = split(data$Day,data$ID) rho_split = split(data$Autocorrelation,data$ID) gen_split = split(data$Genotype,data$ID) #Determine extinction time for each line data_extinction = c(); data_survival = c(); for(id in 1:length(unique(data$ID))){ #For each line data_line = unlist(data_split[id]); #split data by line time_line = unlist(time_split[id]); pos = which(unlist(data_line) Tdeadcontrol){ #For each position, check if N is also under ET for the 4 following transfers list_test = vector("list", min(Tdeadcontrol - 1, length(pos) - Tdeadcontrol + 1)) for (i in 1:min(Tdeadcontrol - 1, length(pos) - Tdeadcontrol + 1)){ list_test[[i]] = pos[diff(pos, lag = i) == i] } if (length(Reduce(intersect, list_test) > 0)){ time_ext = min(Reduce(intersect, list_test)) #if yes, extinction time = first transfer with N < ET } } data_survival = rbind(data_survival, c(as.character(rho_split[id][[1]][1]), as.character(gen_split[id][[1]][1]), time_line[time_ext], ifelse(time_ext == length(time_line), 0, 1))) # timext = ifelse(time_ext == length(time_line), time_ext, time_ext + 1) extinction_line = c(rep(1, time_ext), rep(0, length(time_line) - time_ext)) data_extinction = c(data_extinction, extinction_line) } data$survival_status = data_extinction data_survival = as.data.table(data_survival) colnames(data_survival) = c("Autocorrelation","Genotype","Tdead","census") data_survival$Tdead = as.numeric(data_survival$Tdead) data_survival$census = as.numeric(data_survival$census) # data_survival = data_survival[, Tmax := ifelse(Tdead == Daymax, Inf, Tdead + 1)] # Compute the Kaplan-Meier estimator KM = function(gen, rho){ KMfit = survfit(Surv(time = Tdead, event = as.numeric(census)) ~ Autocorrelation, data=data_survival[Genotype %in% c(gen)& Autocorrelation %in% c(rho)]) return(KMfit) } #COMPARISON BETWEEN TREATMENTS (plot + log rank test) { data_survival$rhoLH = ifelse(data_survival$Autocorrelation %in% c(0,-0.5),"L","H") dataKMtest = data_survival[Genotype %in% c("BC")] KMfit = survfit(Surv(time = Tdead, event = census) ~ rhoLH, data=dataKMtest) ggsurvplot(KMfit, data = dataKMtest, pval = TRUE, pval.method = TRUE) } data_ts = data[survival_status == 1 ] ######## DENSITY - DEPENDANCE ANLYSIS WITH CONSTANT TREATMENTS ########### data_ts = data[survival_status == 1] dataDD = data_ts[Autocorrelation == "cst"] dataDD = dataDD[with(dataDD, order(Well,Plate,Date)),] #Perform the population dynamics analysis unsing TMB #Data NFluo = dataDD$Fcorrfit; NOD = dataDD$ODcorrfit; Ncyto = dataDD$DunaConc; IDlines = as.factor(dataDD$ID) Salinity = dataDD$Serie Genotype = ordered(dataDD$Genotype, levels = c("A","B", "C", "AB", "AC","BC")) timesimul = dataDD$Day; latent_process = log(20000) + 0 * dataDD$S1; data_DD = list(jfluo = NFluo, jOD = NOD, jcyto = Ncyto, ID = IDlines, jtime = timesimul, jSalinity = Salinity, jGenotype = Genotype, dil = dilution) #Parameters parameters = list(r = rep(1, length(unique(Genotype)) * length(unique(Salinity))), log_K = c(rep(log(max(data[Genotype=="A", ConcApprox],na.rm = TRUE)),3), rep(log(max(data[Genotype=="B", ConcApprox],na.rm = TRUE)),3), rep(log(max(data[Genotype=="C", ConcApprox],na.rm = TRUE)),3), rep(log(max(data[Genotype=="AB", ConcApprox],na.rm = TRUE)),3), rep(log(max(data[Genotype=="AC", ConcApprox],na.rm = TRUE)),3), rep(log(max(data[Genotype=="BC", ConcApprox],na.rm = TRUE)),3)), log_beta = 0, log_theta_fluo = -1.5, log_theta_OD = 10, log_size_cyto = 4, jx = latent_process) compile("density_dependence_cst.cpp") dyn.load(dynlib("density_dependence_cst")) #All parameters are estimated no_restri = MakeADFun(data_DD, parameters, random = "jx", smartsearch=TRUE , DLL = "density_dependence_cst") opt_no_restri = optim(no_restri$par, no_restri$fn, no_restri$gr, method = "BFGS", control = list(maxit=10000)) #log_K is not fitted for B lines (lack of data - biologically not coherent results) map_bugB = list(log_K=factor(c(1,2,3,NA,NA,NA,4,5,6,7,8,9,10,11,12,13,14,15))) no_restri_bug_B = MakeADFun(data_DD, parameters, random = "jx", smartsearch=TRUE , DLL = "density_dependence_cst", map = map_bugB) opt_bugB = optim(no_restri_bug_B$par, no_restri_bug_B$fn, no_restri_bug_B$gr, method = "BFGS", control = list(maxit=10000)) #log_K is supposed to be constant accross salinities map_test_K_all = list(log_K=factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6))) restri_all = MakeADFun(data_DD, parameters, random = "jx", smartsearch=TRUE , DLL = "density_dependence_cst", map = map_test_K_all) opt_restri_all = optim(restri_all$par, restri_all$fn, restri_all$gr, method = "BFGS", control = list(maxit=10000)) #log_K is supposed to be constant accross salinities but not estimated for B lines map_test_K_bugB = list(log_K=factor(c(1,1,1,NA,NA,NA,2,2,2,3,3,3,4,4,4,5,5,5))) restri_all_bugB = MakeADFun(data_DD, parameters, random = "jx", smartsearch=TRUE , DLL = "density_dependence_cst", map = map_test_K_bugB) opt_restri_all_bugB = optim(restri_all_bugB$par, restri_all_bugB$fn, restri_all_bugB$gr, method = "BFGS", control = list(maxit=10000)) #log_K is assumed to be the same at 0.8 and 2.4M map_test_K_lowmedium = list(log_K=factor(c(1,1,2,3,3,4,5,5,6,7,7,8,9,9,10,11,11,12))) restri_lowmedium = MakeADFun(data_DD, parameters, random = "jx", smartsearch=TRUE , DLL = "density_dependence_cst", map = map_test_K_lowmedium) opt_restri_lowmedium = optim(restri_lowmedium$par, restri_lowmedium$fn, restri_lowmedium$gr, method = "BFGS", control = list(maxit=10000)) #log_K is assumed to be the same at 0.8 and 2.4M but not estimated for B lines map_test_K_lowmedium_bugB = list(log_K=factor(c(1,1,2,NA,NA,NA,3,3,4,5,5,6,7,7,8,9,9,10))) restri_lowmedium_bugB = MakeADFun(data_DD, parameters, random = "jx", smartsearch=TRUE , DLL = "density_dependence_cst", map = map_test_K_lowmedium_bugB) opt_restri_lowmedium_bugB = optim(restri_lowmedium_bugB$par, restri_lowmedium_bugB$fn, restri_lowmedium_bugB$gr, method = "BFGS", control = list(maxit=10000)) rbind(LRtest(opt_no_restri, opt_restri_lowmedium), LRtest(opt_no_restri, opt_restri_all)) rbind(LRtest(opt_bugB, opt_restri_lowmedium_bugB), LRtest(opt_bugB, opt_restri_all_bugB)) #log_K vary with salinities (Likelihood ratio test significant) but CV is small - supposed constant in # further analysis rep_K_no_restri = sdreport(no_restri_bug_B) rep_K_no_restri_B = sdreport(no_restri) rep_K_bugB = sdreport(restri_all_bugB) logKA = rep_K_bugB[[4]][19] logKB = log(max(data[Genotype == "B"]$ConcApprox,na.rm = TRUE)); logKC = rep_K_bugB[[4]][20] logKAB = rep_K_bugB[[4]][21] logKAC = rep_K_bugB[[4]][22] logKBC = rep_K_bugB[[4]][23] save.image("stage_stochasticR.RData") ####### STOCHASTIC R ANALYSIS - NORMAL + AUTOCORRELATION ############ data_ts = data[survival_status == 1 & Autocorrelation != "conta" ] compile("logistic_data_K_unknownS_normal_autocorr_all_genotype.cpp") dyn.load(dynlib("logistic_data_K_unknownS_normal_autocorr_all_genotype")) #Data datanoS=data_ts[Autocorrelation %in% c(-0.5,0,0.5,0.9)] datanoS = datanoS[with(datanoS, order(ID, Date)),] NFluo = datanoS$Fcorrfit; NOD = datanoS$ODcorrfit; Ncyto = datanoS$DunaConc; rho = datanoS$Autocorrelation rho = ordered(rho, levels = c(-0.5,0,0.5,0.9)) IDlines = as.factor(datanoS$ID) genotypes = datanoS$Genotype genotypes = ordered(genotypes, levels = c("A","B","C","AB","AC","BC")) timesimul = datanoS$Day; latent_process = log(20000) + 0 * datanoS$Day; data_stat = list(jfluo = NFluo, jOD = NOD, jcyto = Ncyto, ID = IDlines, jAutocorrelation = rho, jGenotype = genotypes, jtime = timesimul, dil = dilution) #Parameters parameters = list(mu = rep(0, nlevels(rho) * nlevels(genotypes)), sigma = rep(1, nlevels(rho) * nlevels(genotypes)), rho = rep(0, nlevels(rho) * nlevels(genotypes)), log_K = c(logKA, logKB , logKC, logKAB, logKAC, logKBC), log_theta_fluo = -3, log_theta_OD = 12, log_size_cyto = 2, jx = latent_process ) # State-Space model analysis(TMB) # Carrying capacities could not been estimated with the fluctuating salnities data #(biologically incoherent optimum) { # obj_normal_rho = MakeADFun(data_stat, parameters, random = "jx", smartsearch=TRUE , # DLL = "logistic_data_K_unknownS_normal_autocorr_all_genotype") # opt_normal_rho = optim(obj_normal_rho$par, obj_normal_rho$fn, obj_normal_rho$gr, method = "BFGS", control=list(maxit=1000)) # rep_normal_rho = sdreport(obj_normal_rho) } # Carrying capacities were fixed to their estimations from the constant treatments map_normal_rho1 = list(log_K = factor(c(NA, NA, NA, NA, NA, NA))) # fix K # TMB obj_normal_rho1 = MakeADFun(data_stat, parameters, random = "jx", smartsearch=TRUE , DLL = "logistic_data_K_unknownS_normal_autocorr_all_genotype", map = map_normal_rho1) opt_normal_rho1 = optim(obj_normal_rho1$par, obj_normal_rho1$fn, obj_normal_rho1$gr, method = "BFGS", control=list(maxit=1000)) rep_normal_rho1 = sdreport(obj_normal_rho1) # Extract parameters fixed_normal_rho1 = as.data.table(summary(rep_normal_rho1 , "fixed"),keep.rownames=TRUE) mu_estimated = fixed_normal_rho1[rn == "mu", Estimate]; mu_error = fixed_normal_rho1[rn == "mu", "Std. Error"]; sigma_estimated = fixed_normal_rho1[rn == "sigma", Estimate]; sigma_error = fixed_normal_rho1[rn == "sigma", "Std. Error"]; autocorr_estimated = fixed_normal_rho1[rn == "rho", Estimate]; autocorr_error = fixed_normal_rho1[rn == "rho", "Std. Error"]; rho_estimated = rep(c("-0.5","0","0.5","0.9"), length(unique(genotypes))) genotype_estimated = c(rep(1,length(unique(rho))), rep(6,length(unique(rho))), rep(2,length(unique(rho))), rep(3,length(unique(rho))), rep(4,length(unique(rho))), rep(5,length(unique(rho)))) parameters_estimated_normalrho1 = as.data.table(cbind(genotype_estimated, rho_estimated, mu_estimated,sigma_estimated, autocorr_estimated, mu_error, sigma_error, autocorr_error)) colnames(parameters_estimated_normalrho1) = c("genotype", "rho", "mu","sigma", "Autocorrelation", "mu_error", "sigma_error", "autocorr_error") # Plot parameters estimations mu_plot = ggplot(parameters_estimated_normalrho1[genotype != 6]) + geom_errorbar(aes(x = rho, ymin = mu - mu_error , ymax = mu + mu_error, group = genotype), width = 0.2,position=position_dodge(width=0.3)) + geom_point(aes(x = rho, y = mu, group = genotype), position=position_dodge(width=0.3)) + ggtitle("mu") + theme_bw() sigma_plot = ggplot(parameters_estimated_normalrho1[genotype != 6]) + geom_errorbar(aes(x = rho, ymin = sigma - sigma_error , ymax = sigma + sigma_error, group = genotype), width = 0.3,position=position_dodge(width=0.3)) + geom_point(aes(x = rho, y = sigma, group = genotype), position=position_dodge(width=0.3)) + ggtitle("sigma") + theme_bw() Rho_plot = ggplot(parameters_estimated_normalrho1[genotype != 6]) + geom_errorbar(aes(x = rho, ymin = Autocorrelation - autocorr_error , ymax = Autocorrelation + autocorr_error, group = genotype), width = 0.3,position=position_dodge(width=0.3)) + geom_point(aes(x = rho, y = Autocorrelation, group = genotype, shape = factor(genotype)), position=position_dodge(width=0.3)) + guides(shape = FALSE)+ theme_bw() plot_grid(mu_plot, sigma_plot, Rho_plot, ncol = 1) # extract estimated process: estN_normal_rho1 = as.data.table(cbind(timesimul, IDlines, as.character(genotypes), rho, exp(summary(rep_normal_rho1, "random")[, "Estimate"]))) colnames(estN_normal_rho1) = c("Day", "ID", "gen", "rho", "N") estN_normal_rho1$Day = as.numeric(estN_normal_rho1$Day) N_next = as.numeric(c(estN_normal_rho1$N[2:length(estN_normal_rho1$N)], 0)) DayNext = as.numeric(c(estN_normal_rho1$Day[2:length(estN_normal_rho1$Day)], 0)) log_K_used = ifelse(estN_normal_rho1$gen == "A", logKA, ifelse(estN_normal_rho1$gen == "B", logKB, ifelse(estN_normal_rho1$gen == "C", logKC, ifelse(estN_normal_rho1$gen == "AB", logKAB, ifelse(estN_normal_rho1$gen == "AC", logKAC, logKBC))))) estN_normal_rho1 = as.data.table(cbind(estN_normal_rho1, log_K_used, N_next, DayNext)) estN_normal_rho1$N = as.numeric(estN_normal_rho1$N) deltaT =DayNext - as.numeric(estN_normal_rho1$Day) estN_normal_rho1 = as.data.table(cbind(estN_normal_rho1, deltaT)) Rest = log(estN_normal_rho1$N_next * (exp(estN_normal_rho1$log_K_used) - dilution * estN_normal_rho1$N) / (dilution * estN_normal_rho1$N * (exp(estN_normal_rho1$log_K_used) - estN_normal_rho1$N_next))) / deltaT estN_normal_rho1 = as.data.table(cbind(estN_normal_rho1, Rest)) #Simulate r distribution distri_fit_normal = c() for (geno in 1:6){ for(i in 1:4){ for (r in seq(-3, 2, by=.02)){ dens = dnorm(r, parameters_estimated_normalrho1[(geno - 1) * 4 + i , mu], parameters_estimated_normalrho1[(geno - 1) * 4 + i , sigma]) distri_fit_normal = rbind(distri_fit_normal, c(c("A","B","C","AB","AC","BC")[geno], i, r, dens)) } } } distri_fit_normal = as.data.table(distri_fit_normal) colnames(distri_fit_normal) = c("gen", "i", "r", "dens") save.image("stage_stochasticR.RData") ####### STOCHASTIC R ANALYSIS - GAMMA ############ # Population dynamics analysis using TMB - same as in the normal + autocorrelatino analysis compile("logistic_data_K_unknownS_gamma_rho_all_genotype.cpp") dyn.load(dynlib("logistic_data_K_unknownS_gamma_rho_all_genotype")) data_ts = data_ts[with(data_ts, order(ID)),] datanoS=data_ts[Autocorrelation %in% c(-0.5,0,0.5,0.9), c("ID","Genotype","Day","Autocorrelation","Serie","Replicate","S1","S2", "Fcorrfit","ODcorrfit","DunaConc","ConcApprox")] NFluo = datanoS$Fcorrfit; NOD = datanoS$ODcorrfit; Ncyto = datanoS$DunaConc; rho = datanoS$Autocorrelation rho = ordered(rho, levels = c(-0.5,0,0.5,0.9)) IDlines = as.factor(datanoS$ID) genotypes = datanoS$Genotype genotypes = ordered(genotypes, levels = c("A","B","C","AB","AC","BC")) timesimul = datanoS$Day; latent_process = 0 * datanoS$Day + log(20000); parameters = list(Rmax = mu_estimated+ 2 * sigma_estimated, Shape = rep(4, length(unique(rho))*length(unique(genotypes))), Scale = sigma_estimated / 2, log_K = c(logKA, logKB , logKC, logKAB, logKAC, logKBC), log_theta_fluo = -3, log_theta_OD = 12, log_size_cyto = 2, jx = latent_process ) data_stat = list(jfluo = NFluo, jOD = NOD, jcyto = Ncyto, ID = IDlines, jAutocorrelation = rho, jGenotype = genotypes, jtime = timesimul, dil = dilution) # obj_gamma = MakeADFun(data_stat, parameters, random = c("jx"), smartsearch=TRUE , DLL = "logistic_data_K_unknownS_gamma_rho_all_genotype") # opt_gamma = optim(obj_gamma$par, obj_gamma$fn, obj_gamma$gr, method = "BFGS", control=list(maxit=1000)) # rep_gamma = sdreport(obj_gamma) map_gamma1 = list(log_K = factor(c(NA, NA, NA, NA, NA, NA))) obj_gamma1 = MakeADFun(data_stat, parameters, random = "jx", smartsearch=TRUE , DLL = "logistic_data_K_unknownS_gamma_rho_all_genotype", map = map_gamma1) opt_gamma1 = optim(obj_gamma1$par, obj_gamma1$fn, obj_gamma1$gr, method = "BFGS", control=list(maxit=1000)) rep_gamma1 = sdreport(obj_gamma1) fixed = as.data.table(summary(rep_gamma1),keep.rownames=TRUE) Rmax_estimated = fixed[rn == "Rmax", Estimate]; Shape_estimated = fixed[rn == "Shape", Estimate]; Scale_estimated = fixed[rn == "Scale", Estimate]; Mean_estimated = fixed[rn == "Rmean", Estimate]; Mean_error = fixed[rn == "Rmean", "Std. Error"]; Variance_estimated = fixed[rn == "Variance", Estimate]; Variance_error = fixed[rn == "Variance", "Std. Error"]; Skewness_estimated = fixed[rn == "Skewness", Estimate]; Skewness_error = fixed[rn == "Skewness", "Std. Error"]; rho_estimated = rep(c("-0.5","0","0.5","0.9"), length(unique(genotypes))) genotype_estimated = c(rep(1,length(unique(rho))), rep(6,length(unique(rho))), rep(2,length(unique(rho))), rep(3,length(unique(rho))), rep(4,length(unique(rho))), rep(5,length(unique(rho)))) gen_estimated = c(rep("A",length(unique(rho))), rep("B",length(unique(rho))), rep("C",length(unique(rho))), rep("AB",length(unique(rho))), rep("AC",length(unique(rho))), rep("BC",length(unique(rho)))) parameters_estimated_gamma = as.data.table(cbind(genotype_estimated, gen_estimated, rho_estimated, Mean_estimated,Variance_estimated, Skewness_estimated, Mean_error, Variance_error, Skewness_error, Rmax_estimated, Shape_estimated, Scale_estimated)) colnames(parameters_estimated_gamma) = c("genotype", "gen","rho", "Mean","Variance","Skewness", "Mean_error", "Variance_error", "Skewness_error", "Rmax","Shape","Scale") # extract estimated process: estN_gamma1 = as.data.table(cbind(timesimul, IDlines, as.character(genotypes), datanoS$Autocorrelation, datanoS$S1,datanoS$S2, exp(summary(rep_gamma1, "random")[, "Estimate"]), summary(rep_gamma1, "random")[, "Std. Error"])) colnames(estN_gamma1) = c("Day", "ID", "gen", "rho","S1","S2","N","error") estN_gamma1$Day = as.numeric(estN_gamma1$Day) N_next = as.numeric(c(estN_gamma1$N[2:length(estN_gamma1$N)], 0)) DayNext = as.numeric(c(estN_gamma1$Day[2:length(estN_gamma1$Day)], 0)) log_K_used = ifelse(estN_gamma1$gen == "A", logKA, ifelse(estN_gamma1$gen == "B", logKB, ifelse(estN_gamma1$gen == "C", logKC, ifelse(estN_gamma1$gen == "AB", logKAB, ifelse(estN_gamma1$gen == "AC", logKAC, logKBC))))) estN_gamma1 = as.data.table(cbind(estN_gamma1, log_K_used, N_next, DayNext)) estN_gamma1$N = as.numeric(estN_gamma1$N) deltaT =DayNext - as.numeric(estN_gamma1$Day) estN_gamma1 = as.data.table(cbind(estN_gamma1, deltaT)) Rest = log(estN_gamma1$N_next * (exp(estN_gamma1$log_K_used) - dilution * estN_gamma1$N) / (dilution * estN_gamma1$N * (exp(estN_gamma1$log_K_used) - estN_gamma1$N_next))) / deltaT estN_gamma1 = as.data.table(cbind(estN_gamma1, Rest)) write.csv(estN_gamma1, file = "file:///D:/DONNEESPROFESSIONNELLES100GOMAX/Manip Rescue in Flucuating environment/PDFE/Final_PDFE/estimationR_gamma.csv", row.names = FALSE) save.image("stage_stochasticR.RData") #Generate r distribution for parameter estimation distri_fit_gamma = c() for (geno in 1:6){ for(i in 1:4){ for (r in seq(-3, 2, by=.02)){ dens = dgamma(parameters_estimated_gamma[(geno - 1) * 4 + i , Rmax] - r, shape = parameters_estimated_gamma[(geno - 1) * 4 + i , Shape], scale = parameters_estimated_gamma[(geno - 1) * 4 + i , Scale]) distri_fit_gamma = rbind(distri_fit_gamma, c(c("A","B","C","AB","AC","BC")[geno], i, r, dens)) } } } distri_fit_gamma = as.data.table(distri_fit_gamma) colnames(distri_fit_gamma) = c("gen", "i", "r", "dens") save.image("stage_stochastic.RData") #Save workspace ############### Bivariate TC - 1 genotype ##### # For each gentoypes, paramters of the function r = a + b S1 + c S2 + c S1² + d S2² + f S1 S2 # are estimated, using the population dynamics in all environemental salinity series compile("logistic_data_K.cpp") dyn.load(dynlib("logistic_data_K")) data_ts = data_ts[with(data_ts, order(ID)),] genselect = c("A","B","C","AB","AC","BC") Kselect = c(logKA, logKB, logKC, logKAB, logKAC, logKBC) #Functions { fitdynpop = function(gen, ai, bi, ci, di, ei, fi){ dataline = data_ts[Genotype == genselect[gen]] NFluo = dataline$Fcorrfit; NOD = dataline$ODcorrfit; Ncyto = dataline$DunaConc; rho = dataline$Autocorrelation IDlines = as.factor(dataline$ID) timesimul = dataline$Day; serieS1 = dataline$S1 serieS2 = dataline$S2 latent_process = log(20000) + 0 * dataline$S1; parameters = list(a = ai, b = bi, c = ci, d = di, e = ei, f = fi, log_K = Kselect[gen], log_beta = -1, log_theta_fluo = -3, log_theta_OD = 12, log_size_cyto = 2, jx = latent_process ) data_stat = list(jS1 = serieS1, jS2 = serieS2, jfluo = NFluo, jOD = NOD, jcyto = Ncyto, ID = IDlines, jtime = timesimul, dil = dilution) map1 = list(log_K = factor(NA)) obj1 = MakeADFun(data_stat, parameters, random = "jx", smartsearch=TRUE , DLL = "logistic_data_K", map = map1) opt1 = optim(obj1$par, obj1$fn, obj1$gr, method = "BFGS", control=list(maxit=100)) rep = sdreport(obj1) result = list(opt1, rep) return(result) } tablefitplot = function(gen , repgen){ dataplot=data_ts[Genotype== genselect[gen], c("ID","Day","Autocorrelation","Replicate")] x_evaluated = summary(repgen[[2]], "random")[, "Estimate"] xsd_evaluated = summary(repgen[[2]], "random")[, "Std. Error"] dataplotevaluated = as.data.table(cbind(dataplot, x_evaluated, xsd_evaluated)) return(dataplotevaluated) } } repA = fitdynpop(1, 0.5, 0.2, -0.05, -0.05, -0.05, 0) repB = fitdynpop(2, 0.5, 0.2, -0.05, -0.05, -0.05, 0) repC = fitdynpop(3, 0.5, 0.2, -0.05, -0.05, -0.05, 0) repAB = fitdynpop(4, 0.5, 0.2, -0.05, -0.05, -0.05, 0) repAC = fitdynpop(5, 0.5, 0.2, -0.05, -0.05, -0.05, 0) repBC = fitdynpop(6, 0.5, 0.2, -0.05, -0.05, -0.05, 0) result_acclimation = cbind(repA[[1]]$par,repB[[1]]$par,repC[[1]]$par, repAB[[1]]$par,repAC[[1]]$par,repBC[[1]]$par) result_acclimation = as.data.table(result_acclimation, keep.rownames = TRUE) colnames(result_acclimation) = c("parameter","A","B","C","AB","AC","BC") write.csv(result_acclimation, file = "file:///D:/DONNEESPROFESSIONNELLES100GOMAX/Manip Rescue in Flucuating environment/PDFE/Final_PDFE/result_acclimation.csv", row.names = FALSE) repgen = repC gen = 3 #Extract estimated process { tableGen = tablefitplot(gen, repgen)[Replicate == 1] tableGen = tableGen[, Conc := exp(as.numeric(x_evaluated))] tableGen = tableGen[, ConcPlot := ifelse(Conc < 1000, Conc, 1000 * log(Conc) / log(10) - 2000)] tableGen$Autocorrelation = as.factor(tableGen$Autocorrelation) } #Extract tolerances curves and predict process from tolerance curve { growth_pop = function(S1,S2, sdresult){ r = sdresult[1] + sdresult[2] * S1 + sdresult[3] * S2 + sdresult[4] * S1^2 + sdresult[5] * S2^2 + sdresult[6] * S1 * S2; return(r) } plot_growth = c() for (gen in genselect){ for (s2 in seq(-0.5,4.8,by=0.05)){ for (s1 in seq(-0.5,4.8,by=0.05)) plot_growth = rbind(plot_growth, c(gen,s1,s2,growth_pop(s1, s2, result_acclimation[,..gen]))) } } plot_growth = as.data.table(plot_growth) colnames(plot_growth) = c("gen","S1","S2","R") } save.image("stage_stochastic.RData") #Save workspace ############### Univariate TC - 1 genotype ##### # For each gentoypes, paramters of the function r = a + b S2 + d S2² # are estimated, using the population dynamics in all environemental salinity series compile("logistic_data_K_univariate.cpp") dyn.load(dynlib("logistic_data_K_univariate")) data_ts = data_ts[with(data_ts, order(ID)),] genselect = c("A","B","C","AB","AC","BC") Kselect = c(logKA, logKB, logKC, logKAB, logKAC, logKBC) #Functions { fitdynpopuniv = function(gen, ai, bi, di){ dataline = data_ts[Genotype == genselect[gen]] NFluo = dataline$Fcorrfit; NOD = dataline$ODcorrfit; Ncyto = dataline$DunaConc; rho = dataline$Autocorrelation IDlines = as.factor(dataline$ID) timesimul = dataline$Day; serieS1 = dataline$S1 serieS2 = dataline$S2 latent_process = log(20000) + 0 * dataline$S1; parameters = list(a = ai, b = bi, d = di, log_K = Kselect[gen], log_beta = -1, log_theta_fluo = -3, log_theta_OD = 12, log_size_cyto = 2, jx = latent_process ) data_stat = list(jS1 = serieS1, jS2 = serieS2, jfluo = NFluo, jOD = NOD, jcyto = Ncyto, ID = IDlines, jtime = timesimul, dil = dilution) map1 = list(log_K = factor(NA)) obj1 = MakeADFun(data_stat, parameters, random = "jx", smartsearch=TRUE , DLL = "logistic_data_K_univariate", map = map1) opt1 = optim(obj1$par, obj1$fn, obj1$gr, method = "BFGS", control=list(maxit=100)) rep = sdreport(obj1) result = list(opt1, rep) return(result) } tablefitplot = function(gen , repgen){ dataplot=data_ts[Genotype== genselect[gen], c("ID","Day","Autocorrelation","Replicate")] x_evaluated = summary(repgen[[2]], "random")[, "Estimate"] xsd_evaluated = summary(repgen[[2]], "random")[, "Std. Error"] dataplotevaluated = as.data.table(cbind(dataplot, x_evaluated, xsd_evaluated)) return(dataplotevaluated) } } repAuniv = fitdynpopuniv(1, 0.5, 0.2, -0.05) repBuniv = fitdynpopuniv(2, 0.5, 0.2, -0.05) repCuniv = fitdynpopuniv(3, 0.5, 0.2, -0.05) repABuniv = fitdynpopuniv(4, 0.5, 0.2, -0.05) repACuniv = fitdynpopuniv(5, 0.5, 0.2, -0.05) repBCuniv = fitdynpopuniv(6, 0.5, 0.2, -0.05) rbind(cbind(LRtest(repA[[1]],repAuniv[[1]])[3], "deltaAIC" = AIC(repA[[1]])-AIC(repAuniv[[1]])), cbind(LRtest(repB[[1]],repBuniv[[1]])[3],"deltaAIC" = AIC(repB[[1]])-AIC(repBuniv[[1]])), cbind(LRtest(repC[[1]],repCuniv[[1]])[3],"deltaAIC" = AIC(repC[[1]])-AIC(repCuniv[[1]])), cbind(LRtest(repAB[[1]],repABuniv[[1]])[3],"deltaAIC" = AIC(repAB[[1]])-AIC(repABuniv[[1]])), cbind(LRtest(repAC[[1]],repACuniv[[1]])[3],"deltaAIC" = AIC(repAC[[1]])-AIC(repACuniv[[1]])), cbind( LRtest(repBC[[1]],repBCuniv[[1]])[3],"deltaAIC" = AIC(repBC[[1]])-AIC(repBCuniv[[1]]))) result_no_acclimation = cbind(repAuniv[[1]]$par,repBuniv[[1]]$par,repCuniv[[1]]$par, repABuniv[[1]]$par,repACuniv[[1]]$par,repBCuniv[[1]]$par) result_no_acclimation = as.data.table(result_no_acclimation, keep.rownames = TRUE) colnames(result_no_acclimation) = c("parameter","A","B","C","AB","AC","BC") write.csv(result_no_acclimation, file = "file:///D:/DONNEESPROFESSIONNELLES100GOMAX/Manip Rescue in Flucuating environment/PDFE/Final_PDFE/result_no_acclimation.csv", row.names = FALSE) repgen = repC gen = 3 #Extract estimated process { tableGen = tablefitplot(gen, repgen)[Replicate == 1] tableGen = tableGen[, Conc := exp(as.numeric(x_evaluated))] tableGen = tableGen[, ConcPlot := ifelse(Conc < 1000, Conc, 1000 * log(Conc) / log(10) - 2000)] tableGen$Autocorrelation = as.factor(tableGen$Autocorrelation) } #Extract tolerances curves and predict process from tolerance curve { growth_pop_univ = function(S1,S2, sdresult){ r = sdresult[1] + sdresult[2] * S2 + sdresult[3] * S2^2; return(r) } plot_growth_univ = c() for (gen in genselect){ for (s2 in seq(-0.5,4.8,by=0.05)){ plot_growth_univ = rbind(plot_growth_univ, c(gen,s1,s2,growth_pop(s1, s2, result_no_acclimation[,..gen]))) } } plot_growth_univ = as.data.table(plot_growth_univ) colnames(plot_growth_univ) = c("gen","S1","S2","R") } save.image("stage_stochastic.RData") #Save workspace ############### Simulations : with or without acclimation #### #Extract experimental salinity time series extract_salinity = data[Autocorrelation %in% c(-0.5,0,0.5,0.9) & Genotype == "A" & Replicate == 1] extract_salinity = extract_salinity[with(extract_salinity, order(ID, Autocorrelation, Date))] salinity_simulation = split(extract_salinity$S1, extract_salinity$ID) rho_simulation = split(extract_salinity$Autocorrelation, extract_salinity$ID) #Simulate population dynamics following either the bivariate or the univariate tolerance curve table_simul = c() for(g in 1:6){ x0ref = log(mean(data[Genotype == genselect[g] & Day == 0]$ConcApprox)) K = exp(Kselect[g]); dil = 0.15; #parameters bivariate TC gen = genselect[g] aB = result_acclimation[1, gen, with=FALSE] bB = result_acclimation[2, gen, with=FALSE] cB = result_acclimation[3, gen, with=FALSE] dB = result_acclimation[4, gen, with=FALSE] eB = result_acclimation[5, gen, with=FALSE] fB = result_acclimation[6, gen, with=FALSE] betaB = exp(result_acclimation[7, gen, with=FALSE]) #Parameters univariate TC aU = result_no_acclimation[1, gen, with=FALSE] bU = result_no_acclimation[2, gen, with=FALSE] dU = result_no_acclimation[3, gen, with=FALSE] betaU = exp(result_no_acclimation[4, gen, with=FALSE]) #Loop over lines for (i in 1:length(salinity_simulation)){ x0B = x0ref; x0U = x0ref; S0 = salinity_simulation[[i]][1]; rho_used = rho_simulation[[i]][1] table_simul = rbind(table_simul,c(g,rho_used, i,0, S0, x0ref,x0ref)) #Loop over the time serie for (t in 2:length(salinity_simulation[[i]])){ St = salinity_simulation[[i]][t]; deltat = ifelse(t %% 2 == 0, 3, 4 ) if (x0U > 0) { rU = aU + bU * St + dU * St * St; xtU = log(exp(rU * deltat) * K * dil * exp(x0U) / (K + (exp(rU * deltat) - 1) * dil * exp(x0U))); } else { xtU = 0 } if (x0B > 0) { rB = aB + bB * S0 + cB * St + dB * S0 * S0 + eB * St * St + fB * S0 * St; xtB = log(exp(rB * deltat) * K * dil * exp(x0B) / (K + (exp(rB * deltat) - 1) * dil * exp(x0B))); } else { xtB = 0 } table_simul = rbind(table_simul,c(g, rho_used, i, t - 1, St, xtB, xtU)) x0B = xtB x0U = xtU S0 = St } } } table_simul = as.data.table(table_simul) colnames(table_simul) = c("Genotype","Autocorrelation","Line","Transfer","S","Xbiv","Xuniv") table_simul = table_simul[with(table_simul,order(Line))] table_simul$ID = paste(table_simul$Genotype, table_simul$Autocorrelation,table_simul$Line, sep = "_") table_simul = table_simul[with(table_simul,order(ID))] table_simul$Time = floor(7 * as.numeric(table_simul$Transfer)/2) #Survival analysis # Extinction after 1 transfer with N < 1000 individuals #(threshold < data analysis but no demographoc stochasticity in the simulations) survival_simul = split(table_simul,table_simul$ID) Deadcontrol = 1000 Tdeadcontrol = 0; simul_KM = c(); if(Tdeadcontrol > 1){ for(id in 1:length(unique(table_simul$ID))){ #For each line simul_line = survival_simul[[id]]; #split data by line posUniv = which(simul_line$Xuniv < log(Deadcontrol)) posBiv = which(simul_line$Xbiv < log(Deadcontrol)) time_ext_univ = nrow(simul_line) time_ext_biv = nrow(simul_line) if (length(posUniv) > Tdeadcontrol){ #For each position, check if N is also under ET for the 4 following transfers list_test = vector("list", min(Tdeadcontrol - 1, length(posUniv) - Tdeadcontrol + 1)) for (i in 1:min(Tdeadcontrol - 1, length(posUniv) - Tdeadcontrol + 1)){ list_test[[i]] = posUniv[diff(posUniv, lag = i) == i] } if (length(Reduce(intersect, list_test) > 0)){ time_ext_univ = min(Reduce(intersect, list_test)) #if yes, extinction time = first transfer with N < ET } } if (length(posBiv) > Tdeadcontrol){ #For each position, check if N is also under ET for the 4 following transfers list_test = vector("list", min(Tdeadcontrol - 1, length(posBiv) - Tdeadcontrol + 1)) for (i in 1:min(Tdeadcontrol - 1, length(posBiv) - Tdeadcontrol + 1)){ list_test[[i]] = posBiv[diff(posBiv, lag = i) == i] } if (length(Reduce(intersect, list_test) > 0)){ time_ext_biv = min(Reduce(intersect, list_test)) #if yes, extinction time = first transfer with N < ET } } simul_KM = rbind(simul_KM, c(genselect[simul_line[1]$Genotype], simul_line[1]$Autocorrelation, simul_line[time_ext_univ,Time], ifelse(time_ext_univ == length(simul_line$Time), 0, 1), simul_line[time_ext_biv,Time], ifelse(time_ext_biv == length(simul_line$Time), 0, 1))) } } else { for(id in 1:length(unique(table_simul$ID))){ #For each line simul_line = survival_simul[[id]]; #split data by line posUniv = which(simul_line$Xuniv < log(Deadcontrol)) posBiv = which(simul_line$Xbiv < log(Deadcontrol)) if(length(posBiv) > 0){ time_ext_biv = simul_line[posBiv[1]]$Time census_biv = 1 } else { time_ext_biv = max(simul_line$Time) census_biv = 0 } if(length(posUniv) > 0){ time_ext_univ = simul_line[posUniv[1]]$Time census_univ = 1 } else { time_ext_univ = max(simul_line$Time) census_univ = 0 } simul_KM = rbind(simul_KM, c(genselect[simul_line[1]$Genotype[[1]]], simul_line[1]$Autocorrelation[[1]], time_ext_univ, census_univ, time_ext_biv, census_biv)) } } simul_KM = as.data.table(simul_KM) colnames(simul_KM) = c("Genotype","Autocorrelation","Tdeaduniv","censusuniv","Tdeadbiv","censusbiv") simul_KM$Tdeadbiv = as.numeric(simul_KM$Tdeadbiv) simul_KM$censusbiv = as.numeric(simul_KM$censusbiv) simul_KM$Tdeaduniv = as.numeric(simul_KM$Tdeaduniv) simul_KM$censusuniv = as.numeric(simul_KM$censusuniv) KMfitbiv = survfit(Surv(time = Tdeadbiv, event = as.numeric(censusbiv)) ~ Autocorrelation, data=simul_KM[Genotype != "B"]) KMfituniv = survfit(Surv(time = Tdeaduniv, event = as.numeric(censusuniv)) ~ Autocorrelation, data=simul_KM[Genotype != "B"]) ggsurvplot(KMfitbiv, data = simul_KM, conf.int = TRUE, legend = "none", conf.int.style = "step", xlab = rho, ylab = "", pval = TRUE)+ scale_colour_manual(values = c(rgb(0,0.5,1),rgb(0,2/3,2/3), rgb(1,0.5,0), rgb(2/3,0,0))) ggsurvplot(KMfituniv, data = simul_KM, conf.int = TRUE, legend = "none", conf.int.style = "step", xlab = rho, ylab = "") + scale_colour_manual(values = c(rgb(0,0.5,1),rgb(0,2/3,2/3), rgb(1,0.5,0), rgb(2/3,0,0))) # Population moments in the simulation - univariate { #Compute Population moments { table_stationnarysimul = c() for (genotype in c(1,3,4,5,6)){ for (autocorr in 1:4){ samplestat = table_simul[Genotype == genotype & Autocorrelation == autocorr & Time > 40 & Xuniv > 1] momentteststat = c() for(i in 1:1000){ sample_i_stat = unlist(samplestat$Xuniv[sample(1:length(samplestat$Xuniv), replace=TRUE)]) momentteststat = rbind(momentteststat, c(mean(sample_i_stat), var(sample_i_stat), skewness(sample_i_stat))) } distribrho = c() for(i in 1:1000){ lines_to_take = unique(samplestat$ID)[sample(1:length(unique(samplestat$ID)))] samplerho = samplestat[ID %in% lines_to_take] for (j in 1:length(lines_to_take)){ distribrho = c(distribrho, acf(unlist(samplerho[ID == lines_to_take[j]]$Xuniv), lag.max = 1, plot = FALSE)[[1]][2]) } } table_stationnarysimul = rbind(table_stationnarysimul, c(genselect[genotype], c(-0.5,0,0.5,0.9)[autocorr], mean(momentteststat[,1]), sqrt(var(momentteststat[,1])), quantile(momentteststat[,1], probs = c(0.025, 0.975), na.rm = TRUE), mean(momentteststat[,2]), sqrt(var(momentteststat[,2])), quantile(momentteststat[,2], probs = c(0.025, 0.975), na.rm = TRUE), mean(momentteststat[,3]), sqrt(var(momentteststat[,3])), quantile(momentteststat[,3], probs = c(0.025, 0.975), na.rm = TRUE), mean(distribrho,na.rm = TRUE), sqrt(var(distribrho,na.rm = TRUE)), quantile(distribrho, probs = c(0.025, 0.975), na.rm = TRUE))) } } table_stationnarysimul = as.data.table(table_stationnarysimul) colnames(table_stationnarysimul) = c("genotype", "autocorrelation", "mean", "meansigma","meanlow","meanup", "var", "varsigma","varlow","varup", "skew", "skewsigma","skewlow","skewup", "rho", "rhosigma","rholow","rhoup") } #Linear regression: MICE { mod_mean_simul = list() mod_var_simul = list() mod_skew_simul = list() mod_rho_simul = list() for (i in 1:1000){ sample_N = as.data.table(cbind(mean = rnorm(nrow(table_stationnarysimul), mean = as.numeric(table_stationnarysimul$mean), sd = as.numeric(table_stationnarysimul$meansigma)), var = rnorm(nrow(table_stationnarysimul), mean = as.numeric(table_stationnarysimul$var), sd = as.numeric(table_stationnarysimul$varsigma)), skew = rnorm(nrow(table_stationnarysimul), mean = as.numeric(table_stationnarysimul$skew), sd = as.numeric(table_stationnarysimul$skewsigma)), autocorr = rnorm(nrow(table_stationnarysimul), mean = as.numeric(table_stationnarysimul$rho), sd = as.numeric(table_stationnarysimul$rhosigma)), rho = table_stationnarysimul$autocorrelation, gen = table_stationnarysimul$genotype)) mod_mean_simul[[i]] = lmer(as.numeric(mean) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) mod_var_simul[[i]] = lmer(as.numeric(var) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) mod_skew_simul[[i]] = lmer(as.numeric(skew) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) mod_rho_simul[[i]] = lmer(as.numeric(autocorr) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) } # Combining estimates following Rubin's rule # as.mira takes the list of models and create an object to be used by the mice package pool_mean_simul = pool(as.mira(mod_mean_simul)) pool_var_simul = pool(as.mira(mod_var_simul)) pool_skew_simul = pool(as.mira(mod_skew_simul)) pool_rho_simul = pool(as.mira(mod_rho_simul)) summary(pool_mean_simul) summary(pool_var_simul) summary(pool_skew_simul) summary(pool_rho_simul) } } ####### PLOTS AND STATISTICAL TESTS ############ #Figure 2 : population dynamics analysis { datasize = cbind(summary(rep_gamma1, "random")[, "Estimate"], summary(rep_gamma1, "random")[, "Std. Error"], plotylinlog(cbind(datanoS$Fcorrfit, datanoS$ODcorrfit, datanoS$DunaConc, exp(summary(rep_gamma1, "random")[, "Estimate"]), exp(summary(rep_gamma1, "random")[, "Estimate"] - summary(rep_gamma1, "random")[, "Std. Error"]), exp(summary(rep_gamma1, "random")[, "Estimate"] + summary(rep_gamma1, "random")[, "Std. Error"])))) datasize = as.data.table(datasize) colnames(datasize) = c("log_estimate","error_log","Fluo", "OD", "Cyto","estimate","lower","upper") data_plot_oneline = cbind(datanoS, datasize) #Plot all lines for one genotype { data_plot_dyn = data_plot_oneline[Replicate == 1 & Autocorrelation %in% c(-0.5,0,0.5,0.9)] plot_mean = c() for (gen in c("A","B","C","AB","AC","BC")) { for (rho in c(-0.5, 0, 0.5, 0.9)) { for (day in unique(data_plot_dyn$Day)) { patchx = as.numeric(data_plot_dyn[Genotype == gen & Autocorrelation == rho & Day == day,log_estimate]) SD = ifelse(length(patchx) > 0, sd(patchx), NA) concentration = ifelse(length(patchx) > 0, exp(mean(patchx)), NA) L = exp(mean(patchx) + SD) U = exp(mean(patchx) - SD) plot_mean = rbind(plot_mean, c(gen, rho, day, concentration, U, L,SD )) } } } plot_mean = as.data.table(plot_mean) colnames(plot_mean) = c("Genotype","Autocorr","Day","Concentration","L","U","SD") plot_mean = plot_mean[is.na(plot_mean$Concentration) == FALSE] plot_mean$Concentration = as.numeric(plot_mean$Concentration) plot_mean$L = as.numeric(plot_mean$L) plot_mean$U = as.numeric(plot_mean$U) plot_mean[ , ConcentrationPlot := plotylinlog(Concentration)] plot_mean[ , LPlot := plotylinlog(L)] plot_mean[ , UPlot := plotylinlog(U)] plot_mean$Day = as.numeric(plot_mean$Day) plot_mean = plot_mean[order(plot_mean$Genotype, plot_mean$Autocorr,plot_mean$Day),] popdynplot= function(gen){ ggplot() + geom_line(data = data_plot_dyn[Genotype == gen], aes(x = Day, y = estimate, group = interaction(Autocorrelation, ID), color = Autocorrelation), size = 0.2, alpha = 0.1) + geom_line(data = plot_mean[Genotype == gen], aes(x = Day, y = ConcentrationPlot, group = factor(Autocorr), color = factor(Autocorr)), size = 0.5)+ geom_line(data = plot_mean[Genotype == gen], aes(x = Day, y = UPlot, group = factor(Autocorr), color = factor(Autocorr)), linetype = "dashed", size = 0.5)+ geom_line(data = plot_mean[Genotype == gen], aes(x = Day, y = LPlot, group = factor(Autocorr), color = factor(Autocorr)), linetype = "dashed", size = 0.5)+ geom_segment(aes(x = 0, y = 0, xend = 0, yend = 1000), colour = "black", size = 2) + scale_y_continuous(breaks = breackticks, labels = labelticks,expand = c(0, 0),limits = c(0,plotylinlog(2E6))) + scale_x_continuous(expand = c(0, 0)) + scale_colour_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0), "black")) + labs(x = "Time (days)", y = " Population size", colour = "Autocorrelation") + theme( axis.title.y=element_text(size=10, margin = margin(0,-10,0,-5)), axis.title.x=element_text(size=10), axis.text=element_text(size=8), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid") ) + expand_limits(x = 0, y = 0) } popdynplot("AC") } # plot_hist_N plot_hist_distrib = function(gen){ datadistrib = data_plot_dyn[Genotype == gen & Day > 40] xmax = max(datadistrib$log_estimate) + 1 Mean = mean(xmax - datadistrib$log_estimate) Var = var(xmax - datadistrib$log_estimate) # plot_hist_N = # rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0)) return(ggplot(datadistrib) + geom_freqpoly(aes(x=log_estimate, y=..density.., group = Autocorrelation, color = as.factor(Autocorrelation)), size = 0.5) + # geom_line(data = gammatable, aes(x = x, y = rgammadensity, color = as.factor(rho)),size = 0.8)+ scale_fill_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + scale_color_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0)))+ scale_x_continuous(breaks = log(c(1000,2000,3000,4000,5000,6000,7000,8000,9000, 10000,20000,30000,40000,50000,60000,70000,80000,90000, 100000,200000,300000,400000,500000,600000,700000,800000,900000, 1000000)), labels = c(), expand = c(0, 0), limits = c(log(100),log(2E6)))+ scale_y_continuous(expand = c(0, 0), limits = c(0,0.65)) + guides(fill = FALSE) + guides(color = FALSE) + labs(x = "density")+ theme(axis.title.y=element_blank(), axis.title.x=element_text(size=10), axis.text=element_text(size=8), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) + coord_flip()) } plotC = plot_hist_distrib("C") plotA = plot_hist_distrib("A") #Population moments { table_stationnarystat = c() for (genotype in c("A","C","AB","AC","BC")){ # change for (autocorr in 1:4){ print(paste(genotype,autocorr)) samplestat = data_plot_oneline[Replicate == 1 & Genotype == genotype & Autocorrelation == c(-0.5,0,0.5,0.9)[autocorr] & Day > 40] momentteststat = c() for(i in 1:10000){ num_i_stat = sample(1:nrow(samplestat), replace=TRUE) #bootstrapp # num_i_stat = 1:nrow(samplestat) # without bootstrapp sample_i_stat = rnorm(length(num_i_stat), mean = samplestat[num_i_stat,log_estimate], sd = as.numeric(samplestat[num_i_stat,error_log])) momentteststat = rbind(momentteststat, c(mean(sample_i_stat), var(sample_i_stat), skewness(sample_i_stat))) } distribrho = c() for(i in 1:10000){ lines_to_take = unique(samplestat$ID)[sample(1:length(unique(samplestat$ID)))] #bootstrapp # lines_to_take = unique(samplestat$ID) # without bootstrapp samplerho = samplestat[ID %in% lines_to_take] samplerho$sampleN = rnorm(nrow(samplerho),mean = samplerho$log_estimate, sd = samplerho$error_log) for (j in 1:length(lines_to_take)){ distribrho = c(distribrho, acf(samplerho[ID == lines_to_take[j]]$sampleN, lag.max = 1, plot = FALSE)[[1]][2]) } } save.image("stage_stochastic.RData") #Save workspace table_stationnarystat = rbind(table_stationnarystat, c(genotype, c(-0.5,0,0.5,0.9)[autocorr], mean(momentteststat[,1]), sqrt(var(momentteststat[,1])), quantile(momentteststat[,1], probs = c(0.025, 0.975), na.rm = TRUE), mean(momentteststat[,2]), sqrt(var(momentteststat[,2])), quantile(momentteststat[,2], probs = c(0.025, 0.975), na.rm = TRUE), mean(momentteststat[,3]), sqrt(var(momentteststat[,3])), quantile(momentteststat[,3], probs = c(0.025, 0.975), na.rm = TRUE), mean(distribrho,na.rm = TRUE), sqrt(var(distribrho,na.rm = TRUE)), quantile(distribrho, probs = c(0.025, 0.975), na.rm = TRUE))) } } table_stationnarystat = as.data.table(table_stationnarystat) colnames(table_stationnarystat) = c("genotype", "autocorrelation", "mean", "meansigma","meanlow","meanup", "var", "varsigma","varlow","varup", "skew", "skewsigma","skewlow","skewup", "rho", "rhosigma","rholow","rhoup") } #Linear regression: MICE { mod_mean = list() mod_var = list() mod_skew = list() mod_rho = list() for (i in 1:1000){ sample_N = as.data.table(cbind(mean = rnorm(nrow(table_stationnarystat), mean = as.numeric(table_stationnarystat$mean), sd = as.numeric(table_stationnarystat$meansigma)), var = rnorm(nrow(table_stationnarystat), mean = as.numeric(table_stationnarystat$var), sd = as.numeric(table_stationnarystat$varsigma)), skew = rnorm(nrow(table_stationnarystat), mean = as.numeric(table_stationnarystat$skew), sd = as.numeric(table_stationnarystat$skewsigma)), autocorr = rnorm(nrow(table_stationnarystat), mean = as.numeric(table_stationnarystat$rho), sd = as.numeric(table_stationnarystat$rhosigma)), rho = table_stationnarystat$autocorrelation, gen = table_stationnarystat$genotype)) mod_mean[[i]] = lmer(as.numeric(mean) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) mod_var[[i]] = lmer(as.numeric(var) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) mod_skew[[i]] = lmer(as.numeric(skew) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) mod_rho[[i]] = lmer(as.numeric(autocorr) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_N) } # Combining estimates following Rubin's rule # as.mira takes the list of models and create an object to be used by the mice package pool_mean = pool(as.mira(mod_mean)) pool_var = pool(as.mira(mod_var)) pool_skew = pool(as.mira(mod_skew)) pool_rho = pool(as.mira(mod_rho)) summary(pool_mean) summary(pool_var) summary(pool_skew) summary(pool_rho) } #Plot population moments { meanstat = ggplot(table_stationnarystat[genotype != "B"]) + geom_point(aes(x = as.numeric(autocorrelation) , y = as.numeric(mean), shape = genotype), position=position_dodge(width=0.18), size = 0.7) + geom_errorbar(aes(x = as.numeric(autocorrelation), ymin = as.numeric(meanlow) , ymax = as.numeric(meanup), group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.1) + labs(title = "Mean")+ geom_abline(intercept = summary(pool_mean)[1,1], slope = summary(pool_mean)[2,1])+ scale_shape_manual(values= c(15,3,17,4,16))+ theme(axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) + guides(shape = FALSE) varstat = ggplot(table_stationnarystat[genotype != "B"]) + geom_point(aes(x = as.numeric(autocorrelation) , y = as.numeric(var), shape = genotype), position=position_dodge(width=0.18), size = 0.7) + geom_errorbar(aes(x = as.numeric(autocorrelation), ymin = as.numeric(varlow) , ymax = as.numeric(varup), group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.1) + labs(title ="Variance")+ scale_shape_manual(values= c(15,3,17,4,16))+ geom_abline(intercept = summary(pool_var)[1,1], slope = summary(pool_var)[2,1])+ theme(axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) + guides(shape = FALSE) skewstat = ggplot(table_stationnarystat[genotype != "B"]) + geom_point(aes(x = as.numeric(autocorrelation) , y = as.numeric(skew), shape = genotype), position=position_dodge(width=0.18), size = 0.7) + geom_errorbar(aes(x = as.numeric(autocorrelation), ymin = as.numeric(skewlow) , ymax = as.numeric(skewup), group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.1) + labs(title ="Skewness", x = expression(paste(rho,"(E)")))+ geom_abline(intercept = summary(pool_skew)[1,1], slope = summary(pool_skew)[2,1])+ scale_shape_manual(values= c(15,3,17,4,16))+ theme(axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) + guides(shape = FALSE) rhostat = ggplot(table_stationnarystat[genotype != "B"]) + geom_point(aes(x = as.numeric(autocorrelation) , y = as.numeric(rho), shape = genotype), position=position_dodge(width=0.18), size = 0.7) + geom_errorbar(aes(x = as.numeric(autocorrelation), ymin = as.numeric(rholow) , ymax = as.numeric(rhoup), group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.1) + labs(title ="Autocorrelation", x = expression(paste(rho,"(E)")))+ scale_shape_manual(values= c(15,3,17,4,16))+ geom_abline(intercept = summary(pool_rho)[1,1], slope = summary(pool_rho)[2,1])+ theme(axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) + guides(shape = FALSE) } fig2 = plot_grid(plot_grid(popdynplot("A"), plotA, align = 'h', rel_widths = c(2.3, 1)), plot_grid(meanstat, varstat, skewstat, rhostat,ncol=4,align = 'hv', rel_widths = c(1,1,1,1)), nrow= 2,rel_heights = c(1.7, 1), axis='tblr',align = 'hv') ggsave("Figure2.pdf",fig2, width = 16, height = 11, unit = "cm") } #Figure 3 : growth rate analysis { #from Mathematica meanUnipred = function(p){0.617917} varUnipred = function(p){0.058461} skewUnipred = function(p){-1.73575} rhoUnipred = function(p){0.190902 * p * (4.23829 + 1. * p)} meanApred = function(p){0.578694 + 0.0513912 * p} varApred = function(p){0.141818 - 0.136184 * p + 0.0119903 * p^2} skewApred = function(p){(-5.93986 + 10.3971 * p - 5.15806 * p^2 + 0.263068 * p^3) / (sqrt(0.141818 - 0.136184 * p + 0.0119903 * p^2) * (11.8277 - 11.3578 * p + 1. * p^2))} rhoApred = function(p){(-4.36473 + 9.52358 * p -3.15448 * p^2 - 0.9243 * p^3 + 0.389867 * p^4) / (11.8277 - 11.3578 * p + 1. * p^2)} table_theoretical_moments = c() for (rho in seq(-1,1,0.001)){ table_theoretical_moments = rbind(table_theoretical_moments, c("A",rho, meanApred(rho),varApred(rho),skewApred(rho),rhoApred(rho))) table_theoretical_moments = rbind(table_theoretical_moments, c("uni",rho, meanUnipred(rho),varUnipred(rho),skewUnipred(rho),rhoUnipred(rho))) } table_theoretical_moments = as.data.table(table_theoretical_moments) colnames(table_theoretical_moments) = c("Genotype","Autocorrelation","Mean","Variance","Skewness","Rho") table_theoretical_moments$Mean = as.numeric(table_theoretical_moments$Mean) table_theoretical_moments$Autocorrelation = as.numeric(table_theoretical_moments$Autocorrelation) #Compute rho from the gamma model (realized growth rate) { table_rstat = c() for (gen in c(1,3:6)){ genotype = genselect[gen] for (autocorr in 1:4){ print(paste(gen,autocorr)) samplestat = data_plot_oneline[Replicate == 1 & Genotype == genotype & Autocorrelation == c(-0.5,0,0.5,0.9)[autocorr] & Day > 40] distribrho = c() for(i in 1:10000){ lines_to_take = unique(samplestat$ID)[sample(1:length(unique(samplestat$ID)))] samplerho = samplestat[ID %in% lines_to_take] samplerho$sampleN = rnorm(nrow(samplerho),mean = samplerho$log_estimate, sd = samplerho$error_log) samplerho$sampleN_next = as.numeric(c(samplerho$sampleN[2:length(samplerho$sampleN)], 0)) DayNext = as.numeric(c(samplerho$Day[2:length(samplerho$Day)], 0)) samplerho$deltaT = DayNext - as.numeric(samplerho$Day) samplerho$Rest = log(exp(samplerho$sampleN_next) * (exp(Kselect[gen]) - dilution * exp(samplerho$sampleN)) / (dilution * exp(samplerho$sampleN) * (exp(Kselect[gen]) - exp(samplerho$sampleN_next)))) / samplerho$deltaT for (j in 1:length(lines_to_take)){ distribrho = c(distribrho, acf(samplerho[!is.na(Rest) & ID == lines_to_take[j]]$Rest, lag.max = 1, plot = FALSE)[[1]][2]) } } table_rstat = rbind(table_rstat, c(genotype, c(-0.5,0,0.5,0.9)[autocorr], mean(distribrho,na.rm = TRUE), sqrt(var(distribrho,na.rm = TRUE)), quantile(distribrho, probs = c(0.025, 0.975), na.rm = TRUE))) } } table_rstat = as.data.table(table_rstat) colnames(table_rstat) = c("genotype", "autocorrelation", "rho", "rhosigma","rholow","rhoup") } #Linear regression: MICE { mod_meanR = list() mod_varR = list() mod_skewR = list() mod_rhoR = list() for (i in 1:1000){ sample_R = as.data.table(cbind(mean = rnorm(nrow(parameters_estimated_gamma[gen != "B"]), mean = as.numeric(parameters_estimated_gamma[gen != "B"]$Mean), sd = as.numeric(parameters_estimated_gamma[gen != "B"]$Mean_error)), var = rnorm(nrow(parameters_estimated_gamma[gen != "B"]), mean = as.numeric(parameters_estimated_gamma[gen != "B"]$Variance), sd = as.numeric(parameters_estimated_gamma[gen != "B"]$Variance_error)), skew = rnorm(nrow(parameters_estimated_gamma[gen != "B"]), mean = as.numeric(parameters_estimated_gamma[gen != "B"]$Skewness), sd = as.numeric(parameters_estimated_gamma[genotype != "B"]$Skewness_error)), autocorr = rnorm(nrow(table_rstat), mean = as.numeric(table_rstat$rho), sd = as.numeric(table_rstat$rhosigma)), rho = parameters_estimated_gamma[genotype != "B"]$rho, gen = parameters_estimated_gamma[genotype != "B"]$gen)) mod_meanR[[i]] = lmer(as.numeric(mean) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_R) mod_varR[[i]] = lmer(as.numeric(var) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_R) mod_skewR[[i]] = lmer(as.numeric(skew) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_R) mod_rhoR[[i]] = lmer(as.numeric(autocorr) ~ (1|gen) + as.numeric(rho)*(1|gen), data = sample_R) } # Combining estimates following Rubin's rule # as.mira takes the list of models and create an object to be used by the mice package pool_meanR = pool(as.mira(mod_meanR)) pool_varR = pool(as.mira(mod_varR)) pool_skewR = pool(as.mira(mod_skewR)) pool_rhoR = pool(as.mira(mod_rhoR)) summary(pool_meanR) summary(pool_varR) summary(pool_skewR) summary(pool_rhoR) } #Plot moments Mean_plot = ggplot(parameters_estimated_gamma[genotype != "6"]) + geom_errorbar(aes(x = as.numeric(rho), ymin = Mean - Mean_error , ymax = Mean + Mean_error, group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.18) + geom_line(data = table_theoretical_moments, aes(x = Autocorrelation, y = Mean, linetype = Genotype ), size = 0.5)+ geom_point(aes(x = as.numeric(rho), y = Mean, group = genotype, shape = factor(genotype)), size = 0.8, position=position_dodge(width=0.18)) + guides(shape = FALSE, linetype = FALSE)+ scale_shape_manual(values= c(15,3,17,4,16))+ scale_x_continuous(limits = c(-0.58,0.98))+ labs(x = "", y = "",title="Mean")+ theme(plot.margin=unit(c(0.3,0.3,0,0),"lines"), axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) Variance_plot = ggplot(parameters_estimated_gamma[genotype != "6"]) + geom_errorbar(aes(x = as.numeric(rho), ymin = Variance - Variance_error , ymax = Variance + Variance_error, group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.18) + geom_line(data = table_theoretical_moments, aes(x = Autocorrelation, y = as.numeric(Variance), linetype = Genotype ), size = 0.5)+ geom_point(aes(x = as.numeric(rho), y = Variance, group = genotype, shape = factor(genotype)), size = 0.8, position=position_dodge(width=0.18)) + guides(shape = FALSE, linetype = FALSE)+ scale_shape_manual(values= c(15,3,17,4,16))+ scale_x_continuous(limits = c(-0.58,0.98))+ scale_y_continuous(limits = c(0.03,0.27))+ labs(x = "", y = "",title="Variance")+ theme(plot.margin=unit(c(0.3,0.3,0,0),"lines"), axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) Skewness_plot = ggplot(parameters_estimated_gamma[genotype != "6"]) + geom_errorbar(aes(x = as.numeric(rho), ymin = Skewness - Skewness_error , ymax = Skewness + Skewness_error, group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.18) + geom_line(data = table_theoretical_moments, aes(x = Autocorrelation, y = as.numeric(Skewness), linetype = Genotype ), size = 0.5)+ geom_point(aes(x = as.numeric(rho), y = Skewness, group = genotype, shape = factor(genotype)), size = 0.8, position=position_dodge(width=0.18)) + guides(shape = FALSE, linetype = FALSE)+ scale_shape_manual(values= c(15,3,17,4,16))+ scale_x_continuous(limits = c(-0.58,0.98))+ scale_y_continuous(limits = c(-1.8,-0.45))+ labs(x = "", y = "",title="Skewness")+ theme(plot.margin=unit(c(0.2,0,0.2,0),"lines"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) Rho_plot = ggplot(table_rstat) + geom_errorbar(aes(x = as.numeric(autocorrelation), ymin = as.numeric(rho) - as.numeric(rhosigma) , ymax = as.numeric(rho) + as.numeric(rhosigma), group = genotype), width = 0,position=position_dodge(width=0.18), size = 0.18) + scale_x_continuous(limits = c(-0.58,0.98))+ geom_point(aes(x = as.numeric(autocorrelation), y = as.numeric(rho), group = genotype, shape = factor(genotype)), position=position_dodge(width=0.18), size = 0.8) + geom_line(data = table_theoretical_moments, aes(x = Autocorrelation, y = as.numeric(Rho), linetype = Genotype ), size = 0.5) + guides(shape = FALSE)+ scale_shape_manual(values= c(15,3,17,4,16))+ labs(x = "", title = "Autocorrelation",y="")+ scale_y_continuous(limits = c(-0.55, 0.3))+ theme(plot.margin=unit(c(0.2,0,0.2,0),"lines"), axis.title = element_blank(), plot.title = element_text(size = 10, face = "plain"), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) #Contour plot - genotype c contour_data = as.data.table( rbind(cbind(rho = -0.5,ldply(0.5,ellipse, x = matrix(c(1,-0.5,-0.5,1), nrow=2), scale = c(1.8, 1.8), centre=c(2.4, 2.4))), cbind(rho =0,ldply(0.5,ellipse,x=matrix(c(1,0,0,1), nrow=2), scale = c(1.8, 1.8), centre=c(2.4, 2.4))), cbind(rho =0.5,ldply(0.5,ellipse,x=matrix(c(1,0.5,0.5,1), nrow=2), scale = c(1.8, 1.8), centre=c(2.4, 2.4))), cbind(rho =0.9,ldply(0.5,ellipse,x=matrix(c(1,0.9,0.9,1), nrow=2), scale = c(1.8, 1.8), centre=c(2.4, 2.4))))) CTplot = ggplot(plot_growth[gen == "C"], aes(as.numeric(S2),as.numeric(S1))) + geom_raster(aes(fill = as.numeric(R)), interpolate = TRUE) + geom_polygon(data = contour_data, aes(x = as.numeric(x),y = as.numeric(y), color = factor(rho)), alpha = 0, linetype = 2, size = 1) + geom_contour(aes(z = as.numeric(R)), colour = "black", binwidth = 0.1, size = 0.1) + geom_contour(aes(z = as.numeric(R)), colour = "black", breaks = c(0.54203428), size = 1) + geom_contour(aes(z = as.numeric(R)), colour = "black", binwidth = 1, size = 0.5) + # geom_point(data = data_plot_dyn[Genotype == "C"], aes(x = S2, y = S1, color = Autocorrelation), # size = 0.2, alpha = 0.4)+ labs(x = "Current salinity", y = "Previous salinity")+ scale_fill_gradient2(low ="black", mid = "gray", high = "white", midpoint = 0) + scale_colour_manual(values = c(rgb(0.3,0.5,1),rgb(0.3,0.3,0.3),rgb(1,0.6,0.6),rgb(0.8,0,0))) + scale_x_continuous(limits = c(0,4.8), expand = c(0,0))+ scale_y_continuous(limits = c(0,4.8), expand = c(0,0))+ guides(fill = FALSE, color = FALSE)+ theme(plot.margin=unit(c(0.1,0,0.1,0),"lines"), axis.title.y=element_text(size=10), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) plot_grid(plot_grid(CTplot), plot_grid(Mean_plot, Variance_plot, Skewness_plot, Rho_plot, ncol =2, align = 'vh'), rel_widths = c(1,1),axis = "btlr") plotRmoment = plot_grid(Mean_plot, Variance_plot, Skewness_plot, Rho_plot,ncol =2, align = 'vh') ggsave(filename="Fig3_Rmoments_test.pdf",plot=plotRmoment,width=8,height=8,units="cm") } #Figure 4: Glycerol analysis { estimations_glycerol = c() for (id in unique(data_glycerol$ID0)){ salinity = as.numeric(strsplit(as.character(data_glycerol[ID0 == id][1]$Traitement), split="")[[1]][2]) prevsal = as.numeric(strsplit(as.character(data_glycerol[ID0 == id][1]$Traitement), split="")[[1]][1]) time = data_glycerol[ID0 == id][1]$Temps extra_mean = mean(data_glycerol[ID0 == id & total_extra == "extra"]$GlycerolEstimate, na.rm = TRUE) extra_sigma = sqrt(var(data_glycerol[ID0 == id & total_extra == "extra"]$GlycerolEstimate, na.rm = TRUE))/sqrt(2) total_mean = mean(data_glycerol[ID0 == id & total_extra == "total"]$GlycerolEstimate, na.rm = TRUE) total_sigma = sqrt(var(data_glycerol[ID0 == id & total_extra == "total"]$GlycerolEstimate, na.rm = TRUE))/sqrt(2) total_cells = mean(data_glycerol[ID0 == id & total_extra == "total"]$CytoAlive + data_glycerol[ID0 == id & total_extra == "total"]$CytoDead, na.rm = TRUE) sigma_cells = sqrt(var(data_glycerol[ID0 == id & total_extra == "total"]$CytoAlive + data_glycerol[ID0 == id & total_extra == "total"]$CytoDead, na.rm = TRUE))/sqrt(2) estimations_glycerol = rbind(estimations_glycerol, c(salinity,prevsal, time, total_mean, total_sigma, extra_mean, extra_sigma, total_cells, sigma_cells)) } estimations_glycerol = as.data.table(estimations_glycerol) colnames(estimations_glycerol) = c("salinity","prevsal","time", "total_mean", "total_sigma", "extra_mean", "extra_sigma", "total_cells", "sigma_cells") estimations_glycerol glycerol_intra = c(); for (i in 1:nrow(estimations_glycerol)){ sigma = deltaMethod(object = c("total_mean" = estimations_glycerol[i]$total_mean, "extra_mean" = estimations_glycerol[i]$extra_mean, "total_cells" = estimations_glycerol[i]$total_cells), g = "(total_mean - extra_mean) / total_cells", vcov = matrix(c(estimations_glycerol[i]$total_sigma^2,0,0, 0,estimations_glycerol[i]$extra_sigma^2,0, 0,0,estimations_glycerol[i]$sigma_cells^2),nrow = 3)) glycerol_intra = rbind(glycerol_intra, sigma) } glycerol_intra = as.data.table(glycerol_intra) colnames(glycerol_intra) = c("intra_mean","intra_sigma","intraL","intraU") estimations_glycerol = cbind(estimations_glycerol, glycerol_intra) #Test for differences: wald test { testgly = function(S2,S1,minut){ mu1 = estimations_glycerol[prevsal == S1 & time == minut & salinity == S2[1]]$intra_mean mu2 = estimations_glycerol[prevsal == S1 & time == minut & salinity == S2[2]]$intra_mean sigma1 = estimations_glycerol[prevsal == S1 & time == minut & salinity == S2[1]]$intra_sigma sigma2 = estimations_glycerol[prevsal == S1 & time == minut & salinity == S2[2]]$intra_sigma waldstat = (mu1-mu2)^2 / (sigma1^2 + sigma2^2) p = 1-pchisq(waldstat,df=2) return(p) } table_stat_gycerol = c() for (t in sort(unique(estimations_glycerol[prevsal == 1]$time))){ table_stat_gycerol = rbind(table_stat_gycerol, c(1,1,2,t, testgly(c(1,2),1,t)), c(3,1,2,t, testgly(c(1,2),3,t)), c(3,1,3,t, testgly(c(1,3),3,t)), c(3,2,3,t, testgly(c(2,3),3,t))) } } glycerolplot = ggplot(estimations_glycerol[!(salinity == 3 & prevsal ==1)])+ geom_line(aes(x = time, y = intra_mean, group = interaction(salinity, prevsal), colour = factor(salinity), linetype = factor(prevsal)), size = 0.4)+ geom_errorbar(aes(x = time, ymin = intraL, ymax = intraU, group = interaction(salinity, prevsal), colour = factor(salinity)), size = 0.3)+ geom_point(aes(x = time, y = intra_mean, group = interaction(salinity, prevsal), colour = factor(salinity)), size = 0.5)+ scale_colour_manual(values = c("blue", "purple", "red"))+ scale_linetype_manual(values=c("longdash", "solid"))+ scale_y_continuous(limits = c(0,2.8*10^-7),expand = c(0, 0))+ scale_x_continuous(expand = c(0, 0))+ guides(colour = FALSE)+ theme(plot.margin = margin(5,10,0,0), axis.title.y=element_text(size=10, face = "plain"), axis.title.x=element_text(size=10, face = "plain"), axis.text=element_text(size=8), axis.line.x = element_line(size = 0.4), axis.line.y = element_line(size = 0.4), axis.ticks = element_line(size = 0.4), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid") ) + labs( x = "Time (hours)", y = "Intracellular glycerol") ggsave("Figure4.pdf",glycerolplot, width = 7, height = 7, unit = "cm") } #Figure 5 : survival analysis from data and simulations (univariate vs bivariate TC) { #Plot Kaplan meier curves - all genotypes but B KMfit = survfit(Surv(time = Tdead, event = census) ~ Autocorrelation, data=data_survival[Genotype != "B" & Autocorrelation %in% c(-0.5,0,0.5,0.9)]) F3A = ggsurvplot(KMfit, data = data_survival[Genotype != "B" & Autocorrelation %in% c(-0.5,0,0.5,0.9)], size = 0.5, conf.int = TRUE, legend = "none", conf.int.style = "step", font.x = c(10, "plain", "black"), font.y = c(10, "plain", "black"), font.tickslab = c(8, "plain", "black")) + scale_colour_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + scale_y_continuous( breaks = c(0.25,0.5,1), labels = c("0.25","0.5","1"), trans='log2', expand = c(0, 0.1), limits = c(1/5,1)) #Plot KM curves from simulations with or without acclimation KMfitbiv = survfit(Surv(time = Tdeadbiv, event = as.numeric(censusbiv)) ~ Autocorrelation, data=simul_KM[Genotype != "B"]) KMfituniv = survfit(Surv(time = Tdeaduniv, event = as.numeric(censusuniv)) ~ Autocorrelation, data=simul_KM[Genotype != "B"]) F3B = ggsurvplot(KMfitbiv, data = simul_KM[Genotype != "B"], size = 0.5, conf.int = FALSE, legend = "none", conf.int.style = "step", font.x = c(10, "plain", "black"), font.y = c(10, "plain", "black"), font.tickslab = c(8, "plain", "black")) + scale_colour_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + scale_y_continuous( breaks = c(0.25,0.5,1), labels = c("0.25","0.5","1"), trans='log2', expand = c(0, 0.1), limits = c(1/5,1)) F3C = ggsurvplot(KMfituniv, data = simul_KM[Genotype != "B"], size = 0.5, conf.int = FALSE, legend = "none", conf.int.style = "step", font.x = c(10, "plain", "black"), font.y = c(10, "plain", "black"), font.tickslab = c(8, "plain", "black")) + scale_colour_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + scale_y_continuous( breaks = c(0.25,0.5,1), labels = c("0.25","0.5","1"), trans='log2', expand = c(0, 0.1), limits = c(1/5,1)) Figure5 = arrange_ggsurvplots(list(F3A, F3B,F3C), print = TRUE, title = NA, ncol = 3, nrow = 1) ggsave("Figure5.pdf",Figure5, width = 18,height = 7, units = "cm") } #Supplementary Figure 1: growth rates distribution, tolerance curves and survival analyis in all genotypes { #Compare fit with estimated r plotR = function(table, genotype){ p = ggplot(table[gen == genotype & deltaT >0], aes(x=Rest)) + geom_histogram(aes(group = rho, fill = as.factor(rho), y=..density..), position = position_dodge(width = 0.2), binwidth = 0.2, alpha = 0.5) + geom_line(data = distri_fit_gamma[gen == genotype], aes(x = as.numeric(r), y = as.numeric(dens), color = as.factor(i)), size = 0.2) + scale_fill_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + scale_color_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + scale_x_continuous(limits = c(-3,2)) + guides(fill = FALSE) + guides(color = FALSE) + labs(x = "", y="density") + theme(plot.margin=unit(c(-1,-1,-1,-1),"lines"), plot.title =element_text(size=12, hjust = 0.1, face = "plain"), axis.title.y=element_text(size=10), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) return(p) } plotR(estN_gamma1[deltaT>0], "A") #Compute and plot bivariate tolerance curve plotCT = function(genotype){ CTplot = ggplot(plot_growth[gen == genotype], aes(as.numeric(S2),as.numeric(S1))) + geom_raster(aes(fill = as.numeric(R)), interpolate = TRUE) + geom_contour(aes(z = as.numeric(R)), colour = "gray", binwidth = 0.1, size = 0.02) + geom_contour(aes(z = as.numeric(R)), colour = "black", breaks = c(0), size = 0.2) + geom_contour(aes(z = as.numeric(R)), colour = "black", breaks = c(0.54203428), size = 0.5) + labs(x = "S2", y = "S1")+ scale_fill_gradient2(low ="black", mid = "gray", high = "white", midpoint = 0) + scale_x_continuous(limits = c(0,4.8), expand = c(0,0))+ scale_y_continuous(limits = c(0,4.8), expand = c(0,0))+ guides(fill = FALSE, color = FALSE)+ theme(plot.margin=unit(c(0.1,0,0.1,0),"lines"), axis.title.y=element_text(size=10), axis.title.x=element_blank(), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid")) return(CTplot) } plotCT("C") #Compute and plot KM curves for each genotype plotKMfit = function(genotype){ KMfit = survfit(Surv(time = Tdead, event = census) ~ Autocorrelation, data=data_survival[Genotype == genotype & Autocorrelation %in% c(-0.5,0,0.5,0.9)]) SPlot = ggsurvplot(KMfit, data = data_survival[Genotype == genotype & Autocorrelation %in% c(-0.5,0,0.5,0.9)], conf.int = FALSE, legend = "none", conf.int.style = "step", font.y = c(10, "plain", "black"), font.tickslab = c(8, "plain", "black"), size = 0.5) + scale_colour_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0))) + labs(x="",y="p(survival)")+ scale_x_continuous(limits=c(0,130))+ scale_y_continuous( breaks = c(0.125,0.25,0.5,1), labels = c("1/8","1/4","1/2","1"), trans='log2', expand = c(0, 0.1), limits = c(1/10,1)) return(SPlot) } plotKMfit("AC") supp1_part2 = arrange_ggsurvplots(list(plotKMfit("A"), plotKMfit("B"),plotKMfit("C"),plotKMfit("AB"),plotKMfit("AC"),plotKMfit("BC")), print = TRUE, title = NA, ncol = 1,nrow = 6) supp1_part1 = plot_grid(plotR(estN_gamma1[deltaT>0], "A") + labs(x = ""), plotCT("A")+ labs(x = ""), plotR(estN_gamma1[deltaT>0], "B")+ labs(x = ""), plotCT("B")+ labs(x = ""), plotR(estN_gamma1[deltaT>0], "C")+ labs(x =""),plotCT("C")+ labs(x = ""), plotR(estN_gamma1[deltaT>0], "AB")+ labs(x = ""),plotCT("AB")+ labs(x = ""), plotR(estN_gamma1[deltaT>0], "AC")+ labs(x = ""),plotCT("AC")+ labs(x = ""), plotR(estN_gamma1[deltaT>0], "BC"), plotCT("BC"), ncol = 2,align = 'hv') ggsave("Supp1_part1.pdf",supp1_part1, width = 10,height = 25, units = "cm") ggsave("Supp1_part2.pdf",supp1_part2, width = 5,height = 25, units = "cm") } #Supplementary Figure 2: Popultion dynamics analysis and r and K estiamtion in constant treatments { datasize = cbind(summary(rep_K_bugB, "random")[, "Estimate"], summary(rep_K_bugB, "random")[, "Std. Error"], plotylinlog(cbind(dataDD$Fcorrfit, dataDD$ODcorrfit, dataDD$DunaConc, exp(summary(rep_K_bugB, "random")[, "Estimate"]), exp(summary(rep_K_bugB, "random")[, "Estimate"] - summary(rep_K_bugB, "random")[, "Std. Error"]), exp(summary(rep_K_bugB, "random")[, "Estimate"] + summary(rep_K_bugB, "random")[, "Std. Error"])))) datasize = as.data.table(datasize) colnames(datasize) = c("log_estimate","error_log","Fluo", "OD", "Cyto","estimate","lower","upper") data_plot_cst = cbind(dataDD, datasize) #Plot all lines for one genotype { plot_mean = c() for (gen in c("A","B","C","AB","AC","BC")) { for (serie in c("1","2","3")) { for (day in unique(data_plot_dyn$Day)) { patchx = as.numeric(data_plot_cst[Genotype == gen & Serie == serie & Day == day,log_estimate]) SD = ifelse(length(patchx) > 0, sd(patchx), NA) concentration = ifelse(length(patchx) > 0, exp(mean(patchx)), NA) L = exp(mean(patchx) + SD) U = exp(mean(patchx) - SD) plot_mean = rbind(plot_mean, c(gen, serie, day, concentration, U, L,SD )) } } } plot_mean = as.data.table(plot_mean) colnames(plot_mean) = c("Genotype","Salinity","Day","Concentration","L","U","SD") plot_mean = plot_mean[is.na(plot_mean$Concentration) == FALSE] plot_mean$Concentration = as.numeric(plot_mean$Concentration) plot_mean$L = as.numeric(plot_mean$L) plot_mean$U = as.numeric(plot_mean$U) plot_mean[ , ConcentrationPlot := plotylinlog(Concentration)] plot_mean[ , LPlot := plotylinlog(L)] plot_mean[ , UPlot := plotylinlog(U)] plot_mean$Day = as.numeric(plot_mean$Day) plot_mean = plot_mean[order(plot_mean$Genotype, plot_mean$Salinity,plot_mean$Day),] popdynplotcst= function(gen){ ggplot() + geom_line(data = data_plot_cst[Genotype == gen], aes(x = Day, y = estimate, group = interaction(Serie, ID), color = Serie), size = 0.3, alpha = 0.2) + geom_line(data = plot_mean[Genotype == gen], aes(x = Day, y = ConcentrationPlot, group = factor(Salinity), color = factor(Salinity)), size = 0.5) + geom_line(data = plot_mean[Genotype == gen], aes(x = Day, y = UPlot, group = factor(Salinity), color = factor(Salinity)), linetype = "dashed", size = 0.5)+ geom_line(data = plot_mean[Genotype == gen], aes(x = Day, y = LPlot, group = factor(Salinity), color = factor(Salinity)), linetype = "dashed", size = 0.5)+ geom_segment(aes(x = 0, y = 0, xend = 0, yend = 1000), colour = "black", size = 2) + scale_y_continuous(breaks = breackticks, labels = labelticks,expand = c(0, 0),limits = c(0,plotylinlog(2E6))) + scale_x_continuous(expand = c(0, 0)) + scale_colour_manual(values = c(rgb(0.5,0.8,1),rgb(0,0,0.9),rgb(0,0,0.3), "black")) + labs(x = "Time (days)", y = " Population size", colour = "Autocorrelation") + theme(plot.margin = margin(t = 5, r = 0, b = 0, l = 0, unit = "pt"), axis.title.y=element_text(size=10, margin = margin(t = 0, r = -8, b = 0, l = 0)), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4)) + expand_limits(x = 0, y = 0) } } # plot_hist_N = plot_hist_distrib_cst = function(gen){ datadistrib = data_plot_cst[Genotype == gen & Day > 40] # plot_hist_N = return(ggplot(datadistrib[Day>40]) + geom_freqpoly(aes(x=log_estimate, y=..density.., group = Serie, color = as.factor(Serie)), position = position_dodge(width = .2), binwidth = .3, size = 0.5) + # geom_line(data = gammatable, aes(x = x, y = rgammadensity, color = as.factor(Salinity)),size = 0.8)+ # geom_rect(aes(xmin=log(100), xmax=log(1000), ymin=0, ymax=0.2), fill="white")+ # scale_fill_manual(values = c("#D55E00","#E69F00","#009E73","#0072B2")) + # scale_color_manual(values = c("#D55E00","#E69F00","#009E73","#0072B2"))+ scale_fill_manual(values = c(rgb(0.5,0.8,1),rgb(0,0,0.9),rgb(0,0,0.3))) + scale_color_manual(values = c(rgb(0.5,0.8,1),rgb(0,0,0.9),rgb(0,0,0.3)))+ scale_x_continuous(breaks = log(c(1000,2000,3000,4000,5000,6000,7000,8000,9000, 10000,20000,30000,40000,50000,60000,70000,80000,90000, 100000,200000,300000,400000,500000,600000,700000,800000,900000, 1000000)), labels = c(), expand = c(0, 0), limits = c(log(100),log(2E6)))+ scale_y_continuous(expand = c(0, 0), limits = c(0,1.5)) + guides(fill = FALSE, color = FALSE) + labs(y = "density", x ="") + theme(plot.margin = margin(t = 0, r = -5, b = 0, l = 0, unit = "pt"), axis.title.y=element_text(size=10), axis.title.x=element_blank(), axis.text.x=element_text(size=8), axis.text.y=element_blank(), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4))+ coord_flip())+ expand_limits(x = 0, y = 0) } plot_hist_distrib_cst("A") #Plot r and K { #Plot plot_r_K = as.data.table(summary(rep_K_no_restri_B, "fixed"),keep.rownames=TRUE) strain = c(rep(c(rep(1,3),rep(6,3),rep(2,3),rep(3,3),rep(4,3),rep(5,3)),2),"n","n","n","n") salinity = c(rep(c(0.8, 2.4, 3.2), 12),"n","n","n","n") plot_r_K=cbind(plot_r_K,strain,salinity) colnames(plot_r_K) = c("rn","estimate","error","strain","salinity") plot_r_K=plot_r_K[with(plot_r_K, order(salinity,estimate)),] plot_r = ggplot() + geom_errorbar(data = plot_r_K[rn == "r" & strain != 6], aes(x = salinity, ymin = estimate - error ,ymax = estimate + error, group = strain), width = 0,position=position_dodge(width=0.18), size = 0.18) + geom_point(data = plot_r_K[rn == "r" & strain != 6], aes(x = salinity, y = estimate, group = strain, shape = strain), position=position_dodge(width=0.18), size = 0.8) + scale_shape_manual(values= c(15,3,17,4,16))+ labs(y = "r") + guides(shape = FALSE)+ theme(plot.margin = margin(t = 5, r = 5, b = 0, l = 0, unit = "pt"), axis.title.y=element_text(size=10), axis.title.x=element_blank(), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4)) plot_K=ggplot() + geom_errorbar(data = plot_r_K[rn == "log_K" & strain != 6], aes(x = salinity, ymin = pmax(0, estimate - error),ymax = pmin(16.2, estimate + error), group = strain), width = 0,position=position_dodge(width=0.18), size = 0.18) + geom_point(data = plot_r_K[rn == "log_K" & strain != 6], aes(x = salinity, y = estimate, group = strain, shape = strain), position=position_dodge(width=0.18), size = 0.8) + scale_shape_manual(values= c(15,3,17,4,16))+ scale_y_continuous(limits = c(log(100), 16.2),breaks = c(log(100),log(1000),log(10000),log(100000),log(1000000),log(10000000)), labels = c(100,1000,10000,100000,1000000,10000000),expand = c(0, 0))+ labs(y = "K") + guides(shape = FALSE)+ theme(plot.margin = margin(t = 0, r = 5, b = 0, l = 0, unit = "pt"), axis.title.y=element_text(size=10), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4)) } Supp2 = plot_grid(plot_grid(popdynplotcst("A"),plot_hist_distrib_cst("A"), align = "hv", rel_widths = c(2,1)), plot_grid(plot_r,plot_K,ncol=1,align = 'hv'), ncol = 2, align = 'hv', rel_widths = c(2,1)) ggsave("Supp2.pdf",Supp2, width = 18, height = 8, units= "cm") } #Supplementary Figure 3: Raw measurments analysis (cytometry annd spectrometry) { plotfluo = ggplot() + geom_point(data = data[Fcorr!="NA" & DunaAliveCorr !="NA" & Plate == 1], aes(x = Fcorr, y = DunaConc, color = Day)) + geom_line(data = data[Fcorr!="NA" & DunaAliveCorr != "NA"], aes(x = Fcorr, y = Fcorr * summary(regfluo)$coefficients[1])) + scale_y_log10(limits = c(10,2000000)) + scale_x_continuous(trans='log10') + guides(colour = FALSE) + theme( plot.margin = unit(c(0,0,0,-1),"pt"), axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), axis.text=element_text(size=14), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid") ) + labs( x = "Fluo 685", y = "Population density (cytometer estimation)") + geom_vline(xintercept = 400) plotOD = ggplot() + geom_point(data = data[OD680!="NA" & DunaAliveCorr != "NA" & Plate ==1], aes(x = ODcorr, y = DunaConc, color = Day)) + geom_line(data = data[OD680!="NA" & DunaAliveCorr != "NA"], aes(x = ODcorr, y = ODcorr * summary(regOD)$coefficients[1] # + OD680 * summary(regOD)$coefficients[2] # + OD2 * summary(regOD)$coefficients[3] )) + scale_y_log10(limits = c(10,2000000)) + scale_x_continuous(trans='log10') + theme( plot.margin = unit(c(0,0,0,-1),"pt"), axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), axis.text=element_text(size=14), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line=element_line(linetype = "solid") ) + labs( x = "OD 680", y = "") + geom_vline(xintercept = 0.01) plot_grid(plotfluo, plotOD, ncol = 2, rel_widths = c(1,1) ) } #supplementary Figure 4: Salinity distribution { extract_salinity = data[Autocorrelation %in% c(-0.5,0,0.5,0.9) & Genotype == "A" & Replicate == 1] extract_salinity = extract_salinity[with(extract_salinity, order(ID, Autocorrelation, Date))] salinity_used = split(extract_salinity$S2,extract_salinity$ID) rho_used = split(extract_salinity$Autocorrelation,extract_salinity$ID) plot_S = as.data.table(cbind(rho = as.numeric(as.character(extract_salinity$Autocorrelation)), salinities = as.numeric(extract_salinity$S2))) normal_target = c() for (s in seq(0,5,0.02)){ d = dnorm(s, mean = 2.4, sd = 1, log = FALSE) normal_target = rbind(normal_target, c(s, d)) } normal_target = as.data.table(normal_target) colnames(normal_target) = c("salinity","density") hist_plot_S = ggplot(plot_S) + geom_histogram(aes(x = as.numeric(salinities), y=..density.., group = rho, fill = as.factor(rho)), position = position_dodge(width = .2), binwidth = .2)+ geom_line(data = normal_target, aes(x = salinity, y = density))+ scale_fill_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0)))+ guides(fill = FALSE)+ theme(plot.margin = unit(c(t =0,r=0,b=0,l=0),"pt"), axis.title.y=element_text(size=10), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4)) + labs(x = "")+ coord_flip() Sm05 = ggplot(extract_salinity[Autocorrelation == -0.5 & Day < Daymax & as.numeric(Serie) < 10]) + geom_line(aes(x = Day, y = S2, group = Serie), alpha = 0.7, color = rgb(0,0.43,0.86), size = 0.2)+ theme(plot.margin = unit(c(t =0,r=0,b=0,l=0),"pt"), axis.title.y=element_text(size=10), axis.title.x=element_blank(), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4))+ scale_y_continuous(limits = c(0,5))+ labs(x = "")+ labs(y = "S") S0 = ggplot(extract_salinity[Autocorrelation == 0 & Day < Daymax & as.numeric(Serie) < 10] ) + geom_line(aes(x = Day, y = S2, group = Serie), alpha = 0.7, color = rgb(0,0.57,0.57), size = 0.2)+ theme(plot.margin = unit(c(t =0,r=0,b=0,l=0),"pt"), axis.title.y=element_text(size=10), axis.title.x=element_blank(), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4))+ scale_y_continuous(limits = c(0,5))+ labs(x = "")+ labs(y = "") S05 = ggplot(extract_salinity[Autocorrelation == 0.5 & Day < Daymax & as.numeric(Serie) < 10]) + geom_line(aes(x = Day, y = S2, group = Serie), alpha = 0.7, color = rgb(0.86,0.43,0), size = 0.2)+ theme(plot.margin = unit(c(t =0,r=0,b=0,l=0),"pt"), axis.title.y=element_text(size=10), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4))+ scale_y_continuous(limits = c(0,5))+ labs(x = "days")+ labs(y = "S") S09 = ggplot(extract_salinity[Autocorrelation == 0.9 & Day < Daymax & as.numeric(Serie) < 10]) + geom_line(aes(x = Day, y = S2, group = Serie), alpha = 0.7, color = rgb(0.57,0,0), size = 0.2)+ theme(plot.margin = unit(c(t =0,r=0,b=0,l=0),"pt"), axis.title.y=element_text(size=10), axis.title.x=element_text(size=10), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4))+ scale_y_continuous(limits = c(0,5))+ labs(x = "days")+ labs(y = "") Supp4 = plot_grid(plot_grid(Sm05, S0, S05, S09, ncol =2, align = 'vh'), hist_plot_S, rel_widths = c(2,1.3), axis = "btlr", label_size = 12) ggsave("Supp4.pdf",Supp4, width = 16, height = 8, units= "cm") } #Supplementary Figure 6 : Compare Population sizes estimated with gamma vs normal distribution of r { compareR = as.data.table(cbind(estN_gamma1$gen, estN_gamma1$rho, as.numeric(estN_gamma1$N), as.numeric(estN_normal_rho1$N))) colnames(compareR) = c("gen","rho","normal","gamma") # summary(lm(gamma~normal, data = compareR, na.exclude = TRUE)) compareR$normal = plotylinlog(as.numeric(compareR$normal)) compareR$gamma = plotylinlog(as.numeric(compareR$gamma)) compareR$genotype = ifelse(compareR$gen == "A",1, ifelse(compareR$gen =="B", 6, ifelse(compareR$gen == "C", 2, ifelse(compareR$gen == "AB", 3, ifelse(compareR$gen == "AC", 4,5))))) Ncompare = ggplot(compareR) + geom_point(aes(x = as.numeric(normal), y = as.numeric(gamma), color = factor(rho), shape = factor(genotype)), size = 0.2) + guides(shape = FALSE, color = FALSE) + labs(x = "normal", y = "gamma")+ scale_y_continuous(breaks = breackticks, labels = labelticks,expand = c(0, 0),limits = c(0,plotylinlog(2E6))) + scale_x_continuous(breaks = breackticks, labels = labelticks,expand = c(0, 0),limits = c(0,plotylinlog(2E6))) + scale_color_manual(values = c(rgb(0,0.43,0.86),rgb(0,0.57,0.57),rgb(0.86,0.43,0),rgb(0.57,0,0), "black")) + scale_shape_manual(values= c(15,3,17,4,16,5))+ theme( plot.margin = unit(c(5,5,0,0),"pt"), axis.title.y=element_text(size=10, margin = margin(t = 0, r = -8, b = 0, l = 0)), axis.title.x=element_text(size=10, margin = margin(t = 0, r = 0, b = 5, l = 0)), axis.text=element_text(size=8), legend.position="none", panel.background = element_rect(fill = "white", colour = "White",size = 0, linetype="solid"), axis.line.x=element_line(linetype = "solid", size = 0.4), axis.line.y=element_line(linetype = "solid", size = 0.4), axis.ticks = element_line(size = 0.4) ) + expand_limits(x = 0, y = 0) ggsave("Supp6.pdf",Ncompare, width = 7, height = 7, unit = "cm") } #Plot few single lines { datasize = cbind(summary(rep_gamma1, "random")[, "Estimate"], summary(rep_gamma1, "random")[, "Std. Error"], plotylinlog(cbind(datanoS$Fcorrfit, datanoS$ODcorrfit, datanoS$DunaConc, exp(summary(rep_gamma1, "random")[, "Estimate"]), exp(summary(rep_gamma1, "random")[, "Estimate"] - summary(rep_gamma1, "random")[, "Std. Error"]), exp(summary(rep_gamma1, "random")[, "Estimate"] + summary(rep_gamma1, "random")[, "Std. Error"])))) datasize = as.data.table(datasize) colnames(datasize) = c("log_estimate","error_log","Fluo", "OD", "Cyto","estimate","lower","upper") data_plot_oneline = cbind(datanoS, datasize) id = "C_0_13_1" ggplot() + geom_line(data = data_plot_oneline[ID == id], aes(x = Day, y = estimate), color = "black") + geom_ribbon(data = data_plot_oneline[ID == id],aes(x = Day, ymin = lower, ymax = upper), color = "gray", alpha = 0.1) + geom_point(data = data_plot_oneline[ID == id & NFluo > plotylinlog(4000)], aes(x = Day, y = Fluo), color = gray(0.5)) + geom_point(data = data_plot_oneline[ID == id & NOD > plotylinlog(50000)], aes(x = Day, y = OD), color = gray(0.7)) + geom_point(data = data_plot_oneline[ID == id], aes(x = Day, y = Cyto), color = "black")+ geom_segment(aes(x = 0, y = 0, xend = 0, yend = 1000), colour = "black", size = 2) + scale_y_continuous(limits = c(0,plotylinlog(8e5)), breaks = breackticks, labels = labelticks,expand = c(0, 0)) + scale_x_continuous(limits = c(0,127), expand = c(0, 0)) + labs(x = "days", y = "population size") + expand_limits(x = 0, y = 0) ggplot() + geom_line(data = data_plot_oneline[ID == "A_-0.5_3_1"], aes(x = Day, y = estimate), color = "black") + geom_ribbon(data = data_plot_oneline[ID == "A_-0.5_3_1"],aes(x = Day, ymin = lower, ymax = upper), color = "gray", alpha = 0.1) + geom_point(data = data_plot_oneline[ID == "A_-0.5_3_1"& Fluo > plotylinlog(4000)], aes(x = Day, y = Fluo), color = gray(0.5)) + geom_point(data = data_plot_oneline[ID == "A_-0.5_3_1"& OD > plotylinlog(50000)], aes(x = Day, y = OD), color = gray(0.7)) + geom_point(data = data_plot_oneline[ID == "A_-0.5_3_1"], aes(x = Day, y = Cyto), color = "black") + geom_line(data = data_plot_oneline[ID == "B_0_10_1"], aes(x = Day, y = estimate), color = "black", linetype = 3) + geom_ribbon(data = data_plot_oneline[ID == "B_0_10_1"],aes(x = Day, ymin = lower, ymax = upper), color = "gray", alpha = 0.1, linetype = 3) + geom_point(data = data_plot_oneline[ID == "B_0_10_1"& Fluo > plotylinlog(4000)],aes(x = Day, y = Fluo), color = gray(0.5), shape = 15) + geom_point(data = data_plot_oneline[ID == "B_0_10_1"& OD > plotylinlog(50000)],aes(x = Day, y = OD), color = gray(0.7), shape = 15) + geom_point(data = data_plot_oneline[ID == "B_0_10_1"], aes(x = Day, y = Cyto), color = "black", shape = 15) + geom_line(data = data_plot_oneline[ID == "C_0.9_3_1"], aes(x = Day, y = estimate), color = "black", linetype = 2) + geom_ribbon(data = data_plot_oneline[ID == "C_0.9_3_1"],aes(x = Day, ymin = lower, ymax = upper), color = "gray", alpha = 0.1, linetype = 2) + geom_point(data = data_plot_oneline[ID == "C_0.9_3_1" & Fluo > plotylinlog(4000)], aes(x = Day, y = Fluo), color = gray(0.5), shape = 17) + geom_point(data = data_plot_oneline[ID == "C_0.9_3_1" & OD > plotylinlog(50000)], aes(x = Day, y = OD), color = gray(0.7), shape = 17) + geom_point(data = data_plot_oneline[ID == "C_0.9_3_1"], aes(x = Day, y = Cyto), color = "black", shape = 17) + geom_segment(aes(x = 0, y = 0, xend = 0, yend = 1000), colour = "black", size = 2) + scale_y_continuous(limits = c(0,plotylinlog(8e5)), breaks = breackticks, labels = labelticks,expand = c(0, 0)) + scale_x_continuous(limits = c(0,127), expand = c(0, 0)) + labs(x = "days", y = "population size") + expand_limits(x = 0, y = 0) }