# Kinetic fitting of irreversible inhibition
# Remaining activity = 100 x exp (-kobs*t)
# kobs = kinact*[I] / (Ki + [I])
# Created by S. Fushinobu (2019.1.17)

# Edit the parameters below
irrinh.df <- read.csv("/Users/fushi/Documents/R/rate-irrinh.csv")
Mtitle <- "Irreversible Inhibition"  # main title of the plot
Inhname <- "Inhibitor"
Inhunit <- "mM"
Timeunit <- "min"
kinactinit <- 0.01               # initial values for fitting
Kiinit <- 1.1

# determine inhibitor concentrations used in experiment
Inh.conc <- unique(irrinh.df$conc)

# generate a blank plot (sub1 - v) and then plot the raw data
plot(irrinh.df$time,irrinh.df$rate, pch="", xlim=c(0,max(irrinh.df$time)),ylim=c(0,100),main=paste(Mtitle,"(time course)"),xlab=paste("Time (",Timeunit,")",sep=""),ylab="Remaining Activity (%)")
for (i in 1:length(Inh.conc)) {points(subset(irrinh.df$time, irrinh.df$conc==Inh.conc[i]), subset(irrinh.df$rate, irrinh.df$conc==Inh.conc[i]), pch=i)}

# perform the fitting
irrinh.nls <- nls(rate ~ 100*exp(-(kinact*conc/(Ki + conc))*time), data=irrinh.df, start=list(kinact=kinactinit, Ki=Kiinit) )

# generate a summary of the fit
summary(irrinh.nls)

# confidence intervals of parameters
confint(irrinh.nls)

# extract coefficients
kinact <- unname(coef(irrinh.nls)["kinact"])
Ki <- unname(coef(irrinh.nls)["Ki"])

# plot lines of best fit
for (i in 1:length(Inh.conc)) {curve(100*exp(-(kinact*Inh.conc[i]/(Ki + Inh.conc[i]))*x),0,max(irrinh.df$time), add=TRUE, lty=i)}

# write legends
labels <- Inh.conc
legend("topright", paste(Inhname,labels,Inhunit), pch=c(1:length(Inh.conc)), bg="white")
text(0,0,paste("kinact = ",signif(kinact,digits=3)," 1/",Timeunit,", Ki = ",signif(Ki,digits=3)," ",Inhunit,sep=""),adj=0)

# make a dataframe for kobs
kobsvec <- c()
for (i in 1:length(Inh.conc)) {
  tc.df <- subset(irrinh.df,conc==Inh.conc[i]);
  exp.nls <- nls(rate ~ 100*exp(-kobs*time), data = tc.df,start=list(kobs=kinact));
  kobsvec[i] <- c(unname(coef(exp.nls)["kobs"]))
}
kobs.df <- data.frame(Inhconc=Inh.conc,kobs=kobsvec)

# generate a replot (Inh - kobs) with kobs data and the calculated curve
dev.new()
plot(kobs.df$Inhconc,kobs.df$kobs,xlim=c(0,max(kobs.df$Inhconc)),ylim=c(0,max(kobs.df$kobs)),main=paste(Mtitle," ([",Inhname,"]-kobs)",sep=""),xlab=paste("[",Inhname,"] (",Inhunit,")",sep=""),ylab=paste("kobs (1/",Timeunit,")",sep=""))
curve(kinact*x/(Ki + x),add=TRUE)
text(max(kobs.df$Inhconc),max(kobs.df$kobs)/20,expression(italic(k)[obs] == over(italic(k)[inact] * "[I]", italic(K)[i] + "[I]")),adj=1)
text(0,max(kobs.df$kobs),paste("kinact = ",signif(kinact,digits=3)," 1/",Timeunit,", Ki = ",signif(Ki,digits=3)," ",Inhunit,sep=""),adj=0)

# generate a Lineweaver-Burk-like replot (1/Inh - 1/kobs) with kobs data and the calculated line
dev.new()
plot(1/kobs.df$Inhconc,1/kobs.df$kobs,xlim=c(-1/Ki,max(1/kobs.df$Inhconc)),ylim=c(0,max(1/kobs.df$kobs)),main=paste(Mtitle," (1/[",Inhname,"]-1/kobs)",sep=""),xlab=paste("1/[",Inhname,"] (1/",Inhunit,")",sep=""),ylab=paste("1/kobs (",Timeunit,")",sep=""))
abline(v = 0, h = 0)
abline(1/kinact,Ki/kinact, lty=2)
text(-1/Ki,max(1/kobs.df$kobs)/10,expression(-1/italic(K)[i]))
text(-1/Ki/3,1/kinact,expression(1/italic(k)[inact]))