###  Potential for high toxicity of polystyrene nanoplastics to the European Daphnia longispina  ######
###############################################################################
##Authors: Anderson Abel de Souza Machado, Nesar Ghadernezhad, and Justyna Wolinska


######################### Read me: Script guide  ##############################
###############################################################################

# Section 1: Libraries
# Section 2: Data importing
# Section 3: Dose-response fitting
# Section 4: Further exploratory analyses


#########################  Section 1: Libraries  ##############################
###############################################################################


###### Packages required for Dose-response analyses in R ######

#install.packages("drc")
#install.packages("sandwich")
#install.packages("lmtest")
#install.packages("multcomp")

library(drc)  # Required for fitting dose-response curves
library(sandwich) # to obtain better estimates for dose-response parameters
library(lmtest)   # to obtain better estimates for dose-response parameters
library(multcomp)  # to obtain model comparisions

###### Packages required for Data exploratory analyses in R ######

#install.packages("dplyr")
#install.packages("ggplot2")
#install.packages("cowplot")
#install.packages("nlme")
#install.packages("lme4")
#install.packages("MuMIn")
#install.packages("tidyverse")
#packageVersion("ggplot2")
#install.packages("broom")
#install.packages("ggpubr")

library(tidyverse) # for data manipulation instead to combo dplyr+ggplot
library(cowplot) # for multiple ggplot pannels
library(nlme) # for gls and lme function
library(lme4) # for lmer function
library(MuMIn) # for linear model inference
library(ggpubr) # to add equation in ggplot graphs

# Removing undesired objects from the work place 
rm(list=ls())

#################### Section 2:  Getting all data ###########################
#############################################################################

#Polystyrene density (mg/ nm3)
density.ps= 1050*1e-21

#Particle volume (nm3)
v.ps50nm= ((4/3)*pi*(25)^2) 
v.ps100nm= ((4/3)*pi*(50)^2)
# for transformation to cm? use 1e-21

#Area of particle (nm3)
a.ps50nm= (4*pi*(25)^2) 
a.ps100nm= (4*pi*(50)^2)
# for transformation to cm? use 1e-14

#The computation was not correct using the raw values for particle numbers, and ps.area
#So it was necessary to work with log values (remember a/b= log(a)-log(b))
#Try next weekend with the If else statement for log transformation of the control

###Amme3
Amme3= read.csv("Amme_3.time.csv", header=T, na.strings="na", dec = ",",sep =";")%>% 
  mutate(ps.volume= concentration/density.ps)
View(Amme3)
summary(Amme3)

#small size (50 nm)
Amme3S= Amme3[Amme3$nanoplastic %in% c(50, "control"), ]%>%
  filter(concentration != 0.01) %>% mutate(particle.nr= ps.volume/v.ps50nm)%>%
  mutate(ps.area= particle.nr*a.ps50nm)
View(Amme3S)

#large size (100 nm)
Amme3L= Amme3[Amme3$nanoplastic %in% c(100, "control"), ]%>%
  filter(concentration != 0.01)%>% mutate(particle.nr= ps.volume/v.ps100nm)%>%
  mutate(ps.area= particle.nr*a.ps100nm)
View(Amme3L)


###Amme12
Amme12= read.csv("Amme_12.time.csv", header=T, na.strings="na", dec = ",",sep =";")%>% 
  mutate(ps.volume= concentration/density.ps)
View(Amme12)
summary(Amme12)

#small size (50 nm)
Amme12S= Amme12[Amme12$nanoplastic %in% c(50, "control"), ]%>%
  filter(concentration != 0.01)%>% mutate(particle.nr= ps.volume/v.ps50nm)%>%
  mutate(ps.area= particle.nr*a.ps50nm)
View(Amme12S)

#large size (100 nm)
Amme12L= Amme12[Amme12$nanoplastic %in% c(100, "control"), ]%>%
  filter(concentration != 0.01 & concentration != 0.1 & concentration != 1) %>% 
  mutate(particle.nr= ps.volume/v.ps100nm)%>%
  mutate(ps.area= particle.nr*a.ps100nm)
View(Amme12L)


###Amme51
Amme51= read.csv("Amme_51.time.csv", header=T, na.strings="na", dec = ",",sep =";") %>% 
  mutate(ps.volume= concentration/density.ps)
View(Amme51)
summary(Amme51)

#small size (50 nm)
Amme51S= Amme51[Amme51$nanoplastic %in% c(50, "control"), ] %>%
  filter(concentration != 0.01)%>% mutate(particle.nr= ps.volume/v.ps50nm)%>%
  mutate(ps.area= particle.nr*a.ps50nm)
View(Amme51S)

#large size (100 nm)
Amme51L= Amme51[Amme51$nanoplastic %in% c(100, "control"), ]%>%
  filter(concentration != 0.01) %>% mutate(particle.nr= ps.volume/v.ps100nm)%>%
  mutate(ps.area= particle.nr*a.ps100nm)
View(Amme51L)

###Fluorometry data
se = function(x) { sd(x, na.rm=T) / sqrt(length(na.omit(x)))} # Creating a function for computing standard error
azrange = list(y=c(1000,100000),x=c(0,100))
ayrange = list(y=c(0,150),x=c(0,100))


Fluo= read.csv("Fluo.csv", header=T, na.strings="na", dec = ",",sep =";") %>% 
  group_by(Concentration, NP) %>% 
  summarise_all(funs(mean(., na.rm=T),sd(., na.rm=T))) %>% 
  ungroup() %>% 
  as.data.frame
View(Fluo)
summary(Fluo)

Blank_dif= 1564.4-2056.9


Difference= read.csv("Fluo.csv", header=T, na.strings="na", dec = ",",sep =";") %>% 
  mutate(FluoLoss= 100* (1-(Emmission0-Emmission96)/Emmission0))%>% 
  filter(Concentration != 0) %>% group_by(Concentration, NP) %>% 
  summarise_all(funs(mean(., na.rm=T),sd(., na.rm=T))) %>% 
  ungroup() %>% 
  as.data.frame
  

View(Difference)
summary(Difference)

scientific_10 = function(x) {gsub("e", " x 10^", scientific_format()(x))}

p1= ggplot(Fluo,  aes(Concentration, y = Emmission0_mean)) + scale_color_manual(values=c("#7EC8E3","#0000FF"))+
  geom_point(size=2, aes(colour = NP))+ 
  stat_smooth(formula = y ~ I(log(x+1)), method="lm", se= T, size=0.5, alpha= 0.15, linetype ="longdash", aes(colour = NP)) +
  facet_grid(~NP, scales="fixed") + 
  scale_y_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
    labels = scales::trans_format("log10", scales::math_format(10^.x))) +
  coord_cartesian(ylim = c(1000,100000))+
  labs(x = expression("Polystyrene beads "(mg.L^-1)), y = expression(Fluorescence~at~T[" 0 h"]~(AU))) +
  theme_classic() + theme(axis.line.y=element_blank(),axis.line.x=element_blank())+ 
  annotate("segment",x=-Inf,xend=-Inf, y=azrange$y[1],yend=azrange$y[2]) +
  annotate("segment",y=1000,yend=1000, x=azrange$x[1],xend=azrange$x[2]) + 
  theme(legend.position = "none")

Difference=Fluo%>% filter (Concentration == 0 | Concentration > 1) %>% mutate(Percentage=100* (1-(Emmission0_mean-Emmission96_mean-Blank_dif)/Emmission0_mean))

Fluo_loss= lm (Difference$Percentage ~ Difference$Concentration*Difference$NP)
summary(Fluo_loss)  # The model is significant for effect of concentration and not significant for NM size


p2= ggplot(Difference,  aes(Concentration, y = Percentage)) + scale_color_manual(values=c("#7EC8E3","#0000FF"))+
  geom_point(size=2, aes(colour = NP))+ 
  stat_smooth(method="lm", se= F, size=0.5, alpha= 0.15, linetype ="longdash", aes(colour = NP)) +
  facet_grid(~NP, scales="fixed") + scale_x_continuous(tran= "log10")+
  coord_cartesian(ylim = c(0,150), xlim = c(1,100))+
  labs(x = expression("Polystyrene beads "(mg.L^-1)), y = expression(Fluorescence~Loss~at~T[" 96 h"]~("%"~T[" 0 h"]))) + 
  theme_classic() + theme(legend.position = "none", axis.line.y=element_blank(),axis.line.x=element_blank())+ 
  annotate("segment",x=0,xend=0, y=ayrange$y[1],yend=ayrange$y[2]) +
  annotate("segment",y=-Inf,yend=-Inf, x=ayrange$x[1],xend=ayrange$x[2]) 

plot_grid(p1, p2,labels = "AUTO", ncol = 1)
  

#################### Section 3: Dose-response fitting ########################
##############################################################################

###Amme3
# Nesar and I tried several models (LL3 to LL6) and various CDR 4-5abc in previous versions it seemed that the LL3 was the one best fitting.
# That holds true even when trying to fit Cedergreen models that account for non-monotonic effects
# Here I run the LL3 as it was essentially similar to Cedergreen with 4 parameters I opted to use LL3

#CONCENTRATION
#small mobility
amme3.S.m = drm(mobility~ concentration, time, weights = number, data = Amme3S,
                 fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme3.S.m)
plot(amme3.S.m, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE)

ED(amme3.S.m, c(5, 10, 50), interval = "delta")
ED(amme3.S.m, c(10), interval = "delta")
modelFit(amme3.S.m)

#large mobility
amme3.L.m = drm(mobility ~ concentration, time, weights = number, data = Amme3L,
                 fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme3.L.m)

plot(amme3.L.m, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme3.L.m, c(5, 10, 50), interval = "delta")
ED(amme3.L.m, c(10), interval = "delta")
modelFit(amme3.L.m)


#### Amme12

#small mobility
amme12.S.m = drm(mobility~ concentration, time, weights = number, data = Amme12S,
                 fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme12.S.m)
plot(amme12.S.m, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE)

ED(amme12.S.m, c(5, 10, 50), interval = "delta")
ED(amme12.S.m, c(10), interval = "delta")
modelFit(amme12.S.m)


#large mobility
amme12.L.m = drm(mobility ~ concentration, time, weights = number, data = Amme12L,
                 fct =  LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme12.L.m)
plot(amme12.L.m, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme12.L.m, c(5, 10, 50), interval = "delta")
ED(amme12.L.m, c(10), interval = "delta")
modelFit(amme12.L.m)


#### Amme51

#small mobility
amme51.S.m = drm(mobility~ concentration, time, weights = number, data = Amme51S,
                 fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme51.S.m)
plot(amme51.S.m, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE)

ED(amme51.S.m, c(5, 10, 50), interval = "delta")
ED(amme51.S.m, c(10), interval = "delta")
modelFit(amme51.S.m)

#large mobility
amme51.L.m = drm(mobility ~ concentration, time, weights = number, data = Amme51L,
                 fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme51.L.m)
plot(amme51.L.m, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme51.L.m, c(5, 10, 50), interval = "delta")
ED(amme51.L.m, c(10), interval = "delta")
modelFit(amme51.L.m)

EC10.48h= data.frame(Clone = c("Amme 12", "Amme 3", "Amme 51","Amme 12", "Amme 3", "Amme 51"), EC10mass = c(7.48479, 4.0936, 0.00066303, 69.1255,60.779,25.524), Error.mass =c(10.06901, 1.0856,0.00132480, 164.1892,NA,25.467), Size = c(50,50,50,100,100,100))%>%
  mutate(EC10pn= ifelse(Size== 100, (EC10mass/density.ps)/v.ps100nm, (EC10mass/density.ps)/v.ps50nm)) %>%
  mutate(Error.pn= ifelse(Size== 100, (Error.mass/density.ps)/v.ps100nm, (Error.mass/density.ps)/v.ps50nm)) %>%
  mutate(EC10pa= ifelse(Size== 100, EC10pn/a.ps100nm, EC10pn/a.ps50nm)) %>%
  mutate(Error.pa= ifelse(Size== 100, Error.pn/a.ps100nm, Error.pn/a.ps50nm))

View(EC10.48h)

options(scipen = -1)
################  Compiling All Effect Concentrations  ########################
### Only mobility is considered here as it is the only OECD standard endpoint #

controls = Amme_raw %>% 
  filter(nanoplastic == "control")

control_long = controls %>%
  mutate(nanoplastic = "100") %>% 
  bind_rows(controls %>%
              mutate(nanoplastic = "50"))

Amme_raw = bind_rows(list(Amme3, Amme12, Amme51))
Amme = Amme_raw %>% filter(concentration != 0.01) %>% 
  filter( !(clone== 12 & nanoplastic == 100 & concentration %in% c(0.1, 1) )) %>%  
  filter( nanoplastic != "control") %>%
  bind_rows(control_long) %>%
  group_by(clone, nanoplastic)%>%  
  nest() %>% 
  mutate(model.s = map(data, function(df) drm(survival ~ concentration, time, weights = number, data = df,
    fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial"))) %>%
  mutate(model.m = map(data, function(df) drm(mobility ~ concentration, time, weights = number, data = df,
                                              fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")))

View(Amme)

Amme %>% mutate(res.m = map(model.m, .f = function(x) {ED(x, c(50), interval = "delta")})) %>% mutate(res.s = map(model.s, .f = function(x) {ED(x, c(50), interval = "delta")}))
      
summary(Amme)
#It was taking tricky to use ggplot on the models this way, so I tried another way

models= list("Amme3.100.Mobility"=amme3.L.m, "Amme3.50.Mobility"=amme3.S.m, "Amme12.100.Mobility"=amme12.L.m, "Amme12.50.Mobility"=amme12.S.m, "Amme51.100.Mobility"=amme51.L.m, "Amme51.50.Mobility"=amme51.S.m)

Amme_models= models %>% map_dfr( function(x) { ED(x, c(50), interval = "delta") %>%
    as_tibble()  }, .id = "model")  %>% as.data.frame() %>% add_column(time= rep(c(96,72,48,24), 6))
  
Amme_EC50= cbind(rownames(Amme_models), data.frame(Amme_models, row.names=NULL)) %>% separate(model, sep="\\.", into= c("Clone","Size", "Endpoint"))
View(Amme_EC50)

#Transforming EC50s from mass to log of particle numbers or log of surface area (nm3) 
Amme_EC50.pn= Amme_EC50 %>% mutate (Estimate.v= Estimate / density.ps) %>% 
  mutate (Estimate.log.pn= if_else(Size== 100, log10(Estimate.v)-log10(v.ps100nm), log10(Estimate.v)-log10(v.ps50nm), missing = NULL))%>% 
  mutate (Estimate.log.pa= if_else(Size== 100, Estimate.log.pn+log10(a.ps100nm), Estimate.log.pn+log10(a.ps50nm), missing = NULL))
View(Amme_EC50.pn)  
  
# Computing Species EC50 & EC10 mass-based

EC50sp50= exp(mean(log(c(8.9,8.7,0.2))))
EC50sp100= exp(mean(log(c(90.6,70.7,32.7))))

EC10sp50= exp(mean(log(c(7.48,4.09,0.0007))))
EC10sp100= exp(mean(log(c(69.13,60.78,25.52))))


EC50sp50
EC50sp100
EC10sp50
EC10sp100

###################    Section 4: Graphs for Manuscript      ##################
###############################################################################

par(mfrow=c(1,2))
plot(amme12.S.m, col = TRUE, pch = 20:25,legend=FALSE, bty= "n",
     xlab =expression("Polystyrene beads "["50 nm "]*"(mg.L"^"-1"*")"), 
     ylab = "Fraction of Mobile Daphnia", broken = TRUE)
title(main="A", adj=0, line=0 )
legend("bottomleft", legend=c("24 h", "48 h", "72 h", "96 h"),
       col=c("blue", "green","red", "black"), lty=4:1, bty= "n", cex=0.6)
plot(amme12.L.m, col = TRUE, pch = 20:25,legend= FALSE,bty= "n", xlab = expression("Polystyrene beads "["100 nm "]*"(mg.L"^"-1"*")"), ylab = " ", broken = TRUE)
title(main="B", adj=0, line=0 )

ED(amme12.S.m, c(50), interval = "delta")
ED(amme12.L.m, c(50), interval = "delta")
ED(amme3.S.m, c(50), interval = "delta")
ED(amme3.L.m, c(50), interval = "delta")
ED(amme51.S.m, c(50), interval = "delta")
ED(amme51.L.m, c(50), interval = "delta")


par(mfrow=c(3,2), mai = c(0.6,0.6, 0.2, 0.1))
plot(amme12.S.m, col = c("#050A30","#000C66","#0000FF", "#7EC8E3"), pch = c(15,17,19,21),legend=FALSE, bty= "n",
     xlab = " ", cex.lab=1.5, ylab = "Mobile Amme 12", broken = TRUE)
legend(0.01, 0.7, pch = c(21,19,17,15),legend=c("24 h", "48 h", "72 h", "96 h"),
       col=c("#7EC8E3",  "#0000FF","#000C66", "#050A30"), lty=4:1, bty= "n", cex=0.85, y.intersp= 0.4, x.intersp= 0.4)
title(main="A", adj=0, line=0 )
text(0.22, 0.05, cex=1, expression("EC50"["48h"]== 8.9 %+-% "5.0 "(mg.L^-1)))

plot(amme12.L.m, col = c("#050A30", "#000C66","#0000FF", "#7EC8E3"), pch = c(15,17,19,21),legend= FALSE,bty= "n", xlab =  " ", ylab = " ", broken = TRUE)
title(main="B", adj=0, line=0 )
text(1.4, 0.06, cex=1, expression("EC50"["48h"]== 90.6 %+-% "57.4 "(mg.L^-1)))

plot(amme3.S.m, col = c("#050A30","#000C66","#0000FF", "#7EC8E3"), pch = c(15,17,19,21),legend= FALSE,bty= "n", xlab = " ", cex.lab=1.5, ylab = "Mobile Amme 3", broken = TRUE)
title(main="C", adj=0, line=0 )
text(0.22, 0.05, cex=1, expression("EC50"["48h"]== 8.7 %+-% "1.2 "(mg.L^-1)))

plot(amme3.L.m, col = c("#050A30","#000C66","#0000FF", "#7EC8E3"), pch = c(15,17,19,21),legend= FALSE,bty= "n", xlab = " ", ylab = " ", broken = TRUE)
title(main="D", adj=0, line=0 )
text(0.22, 0.05, cex=1, expression("EC50"["48h"]%~~% 74.7~(mg.L^-1)))

plot(amme51.S.m, col = c("#050A30","#000C66","#0000FF", "#7EC8E3"), pch = c(15,17,19,21),legend= FALSE,bty= "n", cex.lab=1.5, xlab = expression("Polystyrene beads "["50 nm "]*"(mg.L"^"-1"*")"), ylab = "Mobile Amme 51", broken = TRUE)
title(main="F", adj=0, line=0 )
text(0.22, 0.05, cex=1, expression("EC50"["48h"]== 0.2 %+-% "0.1 "(mg.L^-1)))

plot(amme51.L.m, col = c("#050A30","#000C66","#0000FF", "#7EC8E3"), pch = c(15,17,19,21),legend= FALSE,bty= "n", cex.lab=1.5, xlab = expression("Polystyrene beads "["100 nm "]*"(mg.L"^"-1"*")"), ylab = " ", broken = TRUE)
title(main="E", adj=0, line=0 )
text(0.25, 0.05, cex=1, expression("EC50"["48h"]== 32.7 %+-% "41.6 "(mg.L^-1)))
par(mfrow=c(1,1),mai = c(1,1,1,1))

#this figure to be saved as 8 x 6  inches

# Figure 2: EC50s across time and particle types
axrange = list(y=c(0.1,100),x=c(20,90))
apnrange = list(y=c(13,17),x=c(13.5,16.5))
aparange = list(y=c(17,21),x=c(17.5,20.5))


par(mfrow=c(1,3))
p4= ggplot(Amme_EC50,  aes(time, y = Estimate)) + scale_color_manual(values=c("#7EC8E3","#0000FF"))+
  geom_point(size=2, aes(colour = Size))+ 
  stat_smooth(formula = y ~ 0  + I(1/x)+ I((x-1)/x), method="lm", se= T, size=0.5, alpha= 0.15, linetype ="longdash", aes(colour = Size)) +
  facet_grid(~Clone, scales="fixed") + scale_y_continuous( tran= "log10")+
  coord_cartesian(ylim = c(0.05,500))+
  labs(x = "Expoure time (h)", y = expression(EC[50]~(mg.L^-1))) +
  theme_classic() +
  theme(axis.line.y=element_blank(),axis.line.x=element_blank(),legend.position="none")+ 
  annotate("segment",x=-Inf,xend=-Inf, y=axrange$y[1],yend=axrange$y[2]) +
  annotate("segment",y=-Inf,yend=-Inf, x=axrange$x[1],xend=axrange$x[2])
 
p5= ggplot(Amme_EC50.pn,  aes(time, y = Estimate.log.pn)) + scale_color_manual(values=c("#7EC8E3","#0000FF"))+
  geom_point(size=2, aes(colour = Size))+ 
  stat_smooth(formula = y ~ 0  + I(1/x)+ I((x-1)/x), method="lm", se= T, size=0.5, alpha= 0.15, linetype ="longdash", aes(colour = Size)) +
  facet_grid(~Clone, scales="fixed")+ theme_classic()+
  labs(x = "Expoure time (h)", y = expression(EC[50]~(log~count~.L^-1)))+
  theme(axis.line.y=element_blank(),axis.line.x=element_blank(), legend.position="none")+ 
  annotate("segment",x=20,xend=20, y=apnrange$y[1],yend=apnrange$y[2])


p6= ggplot(Amme_EC50.pn,  aes(time, y = Estimate.log.pa)) + scale_color_manual(values=c("#7EC8E3","#0000FF"))+
  geom_point(size=2, aes(colour = Size))+ 
  stat_smooth(formula = y ~ 0  + I(1/x)+ I((x-1)/x), method="lm", se= T, size=0.5, alpha= 0.15, linetype ="longdash", aes(colour = Size)) +
  facet_grid(~Clone, scales="fixed")+ theme_classic()+
  labs(x = "Expoure time (h)", y = expression(EC[50]~(log~area~nm^2~.L^-1)))+
  theme(axis.line.y=element_blank(),axis.line.x=element_blank()) + 
  annotate("segment",x=20,xend=20, y=aparange$y[1],yend=aparange$y[2])

plot_grid(p4, p5, p6, labels = "AUTO", ncol = 1)

#this figure to be saved as 8 x 5  inches

ECs_Time_Clone= lm(Amme_EC50$Estimate~ 0  + I(1/Amme_EC50$time)+ I((Amme_EC50$time-1)/Amme_EC50$time)*Amme_EC50$Size*Amme_EC50$Clone) 

summary(ECs_Time_Clone)

#annotate("segment",y=0.05,yend= 0.05, x=axrange$x[1],xend=axrange$x[2])
#  geom_errorbar(aes(ymin=Estimate-Std..Error,   ymax=Estimate-Std..Error))
# formula that best fit Amme51 50nm is  formula = y ~ 0  + I(1/x)+ I((x-1)/x)

###################  Section 5: Further Exploratoy Analyses  ##################
###############################################################################
# models for mortality. This is not standard endpoint and should not be over-interpreted. 
#Amme 3 small survival
amme3.S.s = drm(survival~ concentration, time, weights = number, data = Amme3S,
                fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme3.S.s)
plot(amme3.S.s, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Survival", broken = TRUE)

ED(amme3.S.s, c(5, 10, 50), interval = "delta")
ED(amme3.S.s, c(50), interval = "delta")
modelFit(amme3.S.s)

#Amme 3 large survival
amme3.L.s = drm(survival ~ concentration, time, weights = number, data = Amme3L,
                fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme3.L.s)
plot(amme3.L.s, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Survival", broken = TRUE, legend=FALSE)
ED(amme3.L.s, c(5, 10, 50), interval = "delta")
ED(amme3.L.s, c(50), interval = "delta")
modelFit(amme3.L.s)

#Amme 12 small survival
amme12.S.s = drm(survival~ concentration, time, weights = number, data = Amme12S,
                 fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme12.S.s)
plot(amme12.S.s, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Survival", broken = TRUE)

ED(amme12.S.s, c(5, 10, 50), interval = "delta")
ED(amme12.S.s, c(50), interval = "delta")
modelFit(amme12.S.s)

#Amme 12 large survival
amme12.L.s = drm(survival ~ concentration, time, weights = number, data = Amme12L,
                 fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme12.L.s)
plot(amme12.L.s, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Survival", broken = TRUE, legend=FALSE)
ED(amme12.L.s, c(5, 10, 50), interval = "delta")
ED(amme12.L.s, c(50), interval = "delta")
modelFit(amme12.L.s)

#Amme 51 small survival
amme51.S.s = drm(survival~ concentration, time, weights = number, data = Amme51S,
                 fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme51.S.s)
plot(amme51.S.s, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Survival", broken = TRUE)

ED(amme51.S.s, c(5, 10, 50), interval = "delta")
ED(amme51.S.s, c(50), interval = "delta")
modelFit(amme51.S.s)


#Amme 51 large survival
amme51.L.s = drm(survival ~ concentration, time, weights = number, data = Amme51L,
                 fct = LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme51.L.s)
plot(amme51.L.s, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Survival", broken = TRUE, legend=FALSE)
ED(amme51.L.s, c(5, 10, 50), interval = "delta")
ED(amme51.L.s, c(50), interval = "delta")
modelFit(amme51.L.s)

#Comparison mortality
models_survival= list("Amme3.100.Survival"=amme3.L.s, "Amme3.50.Survival"=amme3.S.s,"Amme12.100.Survival"=amme12.L.s, "Amme12.50.Survival"=amme12.S.s, "Amme51.100.Survival"=amme51.L.s, "Amme51.50.Survival"=amme51.S.s)

Amme_models_survival= models_survival %>% map_dfr( function(x) { ED(x, c(50), interval = "delta") %>% as_tibble()  }, .id = "model") %>%
  as.data.frame() %>% add_column(time= rep(c(96,72,48,24), 6))

Amme_EC50_survival= cbind(rownames(Amme_models_survival), data.frame(Amme_models_survival, row.names=NULL)) %>% separate(model, sep="\\.", into= c("Clone","Size", "Endpoint"))
View(Amme_EC50_survival)


#ggplot for large
ggplot(Amme12L, aes(concentration, y = survival)) +
  geom_point(size=3, aes(colour = time))+
  stat_smooth(method = "lm", se = TRUE)+
  theme(text=element_text(size=16))+
  facet_wrap(~time, scales="fixed")  + scale_x_continuous(trans = "log10")+
  labs(x = "nanoplastic concentration(log (mg/L)", y = "survival")


#ggplot for small 
ggplot(Amme3S, aes(log(concentration), y = survival)) +
  geom_point(size=3, aes(colour = time))+
  stat_smooth(method = "lm", se = TRUE)+
  theme(text=element_text(size=16))+
  facet_wrap(~time, scales="fixed")+
  labs(x = "PS-NP concentration", y = "survival")

############# Double checking Mobility & PARTICLE NUMBER
#small mobility
amme3.S.m.pn = drm(mobility~ log10(particle.nr+1), time, weights = number, data = Amme3S,
                   fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme3.S.m.pn)
plot(amme3.S.m.pn, col = TRUE, pch = 20:25, xlab = "Particle Number", ylab = "Mobility", broken = TRUE)

ED(amme3.S.m.pn, c(5, 10, 50), interval = "delta")
ED(amme3.S.m.pn, c(50), interval = "delta")
modelFit(amme3.S.m.pn)

#large mobility
amme3.L.m.pn = drm(mobility ~ log10(particle.nr+1), time, weights = number, data = Amme3L,
                   fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme3.L.m.pn)

plot(amme3.L.m.pn, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme3.L.m.pn, c(5, 10, 50), interval = "delta")
ED(amme3.L.m.pn, c(50), interval = "delta")
modelFit(amme3.L.m.pn)


#### Amme12

#small mobility
amme12.S.m.pn = drm(mobility~ log10(particle.nr+1), time, weights = number, data = Amme12S,
                    fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme12.S.m.pn)
plot(amme12.S.m.pn, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE)

ED(amme12.S.m.pn, c(5, 10, 50), interval = "delta")
ED(amme12.S.m.pn, c(50), interval = "delta")
modelFit(amme12.S.m.pn)


#large mobility
amme12.L.m.pn = drm(mobility ~ log10(particle.nr+1), time, weights = number, data = Amme12L,
                    fct =  LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme12.L.m.pn)
plot(amme12.L.m.pn, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme12.L.m.pn, c(5, 10, 50), interval = "delta")
ED(amme12.L.m.pn, c(50), interval = "delta")
modelFit(amme12.L.m.pn)


#### Amme51

#small mobility
amme51.S.m.pn = drm(mobility~ log10(particle.nr+1), time, weights = number, data = Amme51S,
                    fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme51.S.m.pn)
plot(amme51.S.m.pn, col = TRUE, pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE)

ED(amme51.S.m.pn, c(5, 10, 50), interval = "delta")
ED(amme51.S.m.pn, c(50), interval = "delta")
modelFit(amme51.S.m.pn)

#large mobility
amme51.L.m.pn = drm(mobility ~ log10(particle.nr+1), time, weights = number, data = Amme51L,
                    fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme51.L.m.pn)
plot(amme51.L.m.pn, col = TRUE, type= "average", pch = 20:25, xlab = "Concentration (mg/L)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme51.L.m.pn, c(5, 10, 50), interval = "delta")
ED(amme51.L.m.pn, c(50), interval = "delta")
modelFit(amme51.L.m.pn)

#PARTICLE AREA
#small mobility
amme3.S.m.pa = drm(mobility~ log10(ps.area+1), time, weights = number, data = Amme3S,
                   fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme3.S.m.pa)
plot(amme3.S.m.pa, col = TRUE, pch = 20:25, xlab = "Particle Area (log nm3)", ylab = "Mobility", broken = TRUE)

ED(amme3.S.m.pa, c(5, 10, 50), interval = "delta")
ED(amme3.S.m.pa, c(50), interval = "delta")
modelFit(amme3.S.m.pa)

#large mobility
amme3.L.m.pa = drm(mobility ~ log10(ps.area+1), time, weights = number, data = Amme3L,
                   fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme3.L.m.pa)

plot(amme3.L.m.pa, col = TRUE, type= "average", pch = 20:25, xlab = "Particle Area (log nm3)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme3.L.m.pa, c(5, 10, 50), interval = "delta")
ED(amme3.L.m.pa, c(50), interval = "delta")
modelFit(amme3.L.m.pa)


#### Amme12

#small mobility NOT working
amme12.S.m.pa = drm(mobility~ log10(ps.area+1), time, weights = number, data = Amme12S,
                    fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme12.S.m.pa)
plot(amme12.S.m.pa, col = TRUE, pch = 20:25, xlab = "Particle Area (log nm3)", ylab = "Mobility", broken = TRUE)

ED(amme12.S.m.pa, c(5, 10, 50), interval = "delta")
ED(amme12.S.m.pa, c(50), interval = "delta")
modelFit(amme12.S.m.pa)


#large mobility
amme12.L.m.pa = drm(mobility ~ log10(ps.area+1), time, weights = number, data = Amme12L,
                    fct =  LL.3(names = c("slope", "max_value", "LC_50")), type = "binomial")
summary(amme12.L.m.pa)
plot(amme12.L.m.pa, col = TRUE, type= "average", pch = 20:25, xlab = "Particle Area (log nm3)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme12.L.m.pa, c(5, 10, 50), interval = "delta")
ED(amme12.L.m.pa, c(50), interval = "delta")
modelFit(amme12.L.m.pa)


#### Amme51

#small mobility
amme51.S.m.pa = drm(mobility~ log10(ps.area+1), time, weights = number, data = Amme51S,
                    fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme51.S.m.pa)
plot(amme51.S.m.pa, col = TRUE, pch = 20:25, xlab = "Particle Area (log nm3)", ylab = "Mobility", broken = TRUE)

ED(amme51.S.m.pa, c(5, 10, 50), interval = "delta")
ED(amme51.S.m.pa, c(50), interval = "delta")
modelFit(amme51.S.m.pa)

#large mobility  NOT WORKING
amme51.L.m.pa = drm(mobility ~ log10(ps.area+1), time, weights = number, data = Amme51L,
                    fct = LL.3(names = c("slope", "max_value", "EC_50")), type = "binomial")
summary(amme51.L.m.pa)
plot(amme51.L.m.pa, col = TRUE, type= "average", pch = 20:25, xlab = "Particle Area (log nm3)", ylab = "Mobility", broken = TRUE, legend=FALSE)
ED(amme51.L.m.pa, c(5, 10, 50), interval = "delta")
ED(amme51.L.m.pa, c(50), interval = "delta")
modelFit(amme51.L.m.pa)