# Time dependence fitting
# Exponential (psendo-first-order) kinetics starting from low substrate concentration
# Multiple curves version
# Written by S. Fushinobu (2020.1.20)
# Edit the parameters below
ExpLin.df <- read.csv("/Users/fushi/Documents/R/Pseudo1st.csv")
Mtitle <- "Psendo-first-order (exponential) fit"  # main title of the plot
Xaxisname <- "Time (h)"
Yaxisname <- "Abs (400 nm)"
krateinit <- 0.5              # initial values for fitting.
Pmaxinit <- 0.3               # Put the maximum value of prodct here.
offsetinit <- 1
mono <- FALSE              # FALSE for color, TRUE for monochrome

ncurve <- ncol(ExpLin.df) -1
if (mono) {
plot.color <- rep(1,ncurve)
} else {
plot.color <- rainbow(ncurve)
}

#plot blank
plot(0,0, pch="", xlim=c(0,max(ExpLin.df[,1], na.rm=TRUE)),ylim=c(0,max(ExpLin.df[,-1], na.rm=TRUE)), main=Mtitle,xlab=Xaxisname,ylab=Yaxisname)

# write legends
labels <- colnames(ExpLin.df)[-1]
legend("topleft", legend = labels, pch=c(1:length(labels)), lty=c(1:length(labels)), col=plot.color[1:length(labels)], bg="white")

for (i in 1:ncurve){
ithdata.df <- ExpLin.df[,1]
ithdata.df <- cbind(ithdata.df,ExpLin.df[i+1])
colnames(ithdata.df) <- c("time","Pconc")
ithdata.df <- subset(ithdata.df, Pconc != "NA")
points(ithdata.df$time,ithdata.df$Pconc, col=plot.color[i], pch=i)

Expfit2.nls <- nls(Pconc ~ (Ak/krate) * (1 - exp(-krate*time)+offset), data=ithdata.df, start=list(Ak=Pmaxinit*krateinit, krate=krateinit, offset=offsetinit))
summary(Expfit2.nls)$parameters
Ak <- unname(coef(Expfit2.nls)["Ak"])
krate <- unname(coef(Expfit2.nls)["krate"])
offset <- unname(coef(Expfit2.nls)["offset"])
Pmax <- Ak / krate
curve((Ak/krate)*(1 - exp(-krate*x)+offset),0,max(ithdata.df$time),add=TRUE,col=plot.color[i],lty=1)

mtext(paste(colnames(ExpLin.df)[i+1],":Ak=",signif(Ak,digits=3),"± ",signif(summary(Expfit2.nls)$parameters[1,2],digits=3)),side=1,line=-i,adj=0)

#linearMod <- lm(Pconc ~ time, data=ithdata.df)      #For linear data
#summary(linearMod)$coefficients
#Aval <- unname(coef(linearMod)["time"])
#Bval <- unname(coef(linearMod)["(Intercept)"])
#abline(linearMod,col=plot.color[i],lty=2)

#mtext(paste(colnames(ExpLin.df)[i+1],": slope=",signif(Aval,digits=3),"± ",signif(summary(linearMod)$coefficients[2,2],digits=3)),side=1,line=-i,adj=0)

}