# Kinetic fitting of bi-substrate reaction to determine Ping-Pong or sequential mechanisms
# Created by S. Fushinobu (2018.11.28)

# Edit the parameters below
bisub.df <- read.csv("/Users/fushi/Downloads/PP-R/rate-pp.csv")
Mtitle <- "Bi-substrate sequential Mechanism"  # main title of the plot
Sub2name <- "Substrate 2"
Sub2unit <- "mM"
Sub1name <- "Substrate 1"
Sub1unit <- "mM"
Vunit <- " (arb-unit)"
KmAinit <- 2.5               # initial values for fitting
KmBinit <- 0.5
KiAinit <- 1.0
KiBinit <- 1.0
Vmaxinit <- 2.7
grid <- 10                 # grid number for best fit plot lines
mono <- FALSE              # FALSE for color, TRUE for monochrome
AorB <- TRUE              # FALSE for plotting KiB*KmA, TRUE for plotting KiA*KmB

# determine concentrations of Substrates 1 & 2 used in experiment
Sub1.conc <- unique(bisub.df$conc1)
Sub2.conc <- unique(bisub.df$conc2)

# create colors for each Substrates 1 & 2 concentration
if (mono) {
Sub1.color <- rep(1,length(Sub1.conc))
Sub2.color <- rep(1,length(Sub2.conc))
} else {
Sub1.color <- rainbow(length(Sub1.conc))
Sub2.color <- rainbow(length(Sub2.conc))
}

# perform the fitting with KiA
bisub1.nls <- nls(rate ~  (Vmax*conc1*conc2) / (KmA*conc2 + KmB*conc1 + conc1*conc2 + KiA*KmB), data=bisub.df, start=list(KmA=KmAinit, Vmax=Vmaxinit, KmB=KmBinit, KiA=KiAinit) )
summary(bisub1.nls)

# perform the fitting with KiB
bisub2.nls <- nls(rate ~  (Vmax*conc1*conc2) / (KmA*conc2 + KmB*conc1 + conc1*conc2 + KiB*KmA), data=bisub.df, start=list(KmA=KmAinit, Vmax=Vmaxinit, KmB=KmBinit, KiB=KiBinit) )
summary(bisub2.nls)

# extract coefficients
if (AorB){
Vmax <- unname(coef(bisub1.nls)["Vmax"])
KmA <- unname(coef(bisub1.nls)["KmA"])
KmB <- unname(coef(bisub1.nls)["KmB"])
KiAB <- unname(coef(bisub1.nls)["KiA"])
KmAB <- unname(coef(bisub1.nls)["KmB"])
Kilabel <- Sub1name
} else{
Vmax <- unname(coef(bisub2.nls)["Vmax"])
KmA <- unname(coef(bisub2.nls)["KmA"])
KmB <- unname(coef(bisub2.nls)["KmB"])
KiAB <- unname(coef(bisub2.nls)["KiB"])
KmAB <- unname(coef(bisub1.nls)["KmA"])
Kilabel <- Sub2name
}

#prepare for double-reciplocal (Lineweaver-Burk) plot
dr.df <- data.frame(rconc1=1/bisub.df$conc1, rrate=1/bisub.df$rate, rconc2=1/bisub.df$conc2, conc1=bisub.df$conc1, conc2=bisub.df$conc2)
Sub1.rconc <- unique(dr.df$rconc1)
Sub2.rconc <- unique(dr.df$rconc2)

# generate a blank plot (1/sub1 - 1/v) and then plot the raw data
plot(dr.df$rconc1,dr.df$rrate, pch="", xlim=c(0,max(dr.df$rconc1)),ylim=c(0,max(dr.df$rrate)),main=paste(Mtitle," (1/",Sub1name," - 1/v)",sep=""),xlab=paste("1 /",Sub1name),ylab="1/v")
for (i in 1:length(Sub2.rconc)) {points(subset(dr.df$rconc1, dr.df$rconc2==Sub2.rconc[i]), subset(dr.df$rrate, dr.df$rconc2==Sub2.rconc[i]), col=Sub2.color[i], pch=i)}

#draw lines
for (i in 1:length(Sub2.rconc)) {  abline(KmB/Vmax*Sub2.rconc[i] + 1/Vmax, KmA/Vmax + KiAB*KmAB/Vmax, col=Sub2.color[i], lty=i) }
abline(v = 0, h = 0)

# write legends
labels <- Sub2.conc
legend("topleft", paste(Sub2name,labels,Sub2unit), col=Sub2.color, pch=c(1:length(Sub2.conc)), bg="white")
text(max(dr.df$rconc1)/2,0,paste("Km(",Sub1name,") = ",signif(KmA,digits=3)," ",Sub1unit,", Km(",Sub2name,") = ",signif(KmB,digits=3)," ",Sub2unit," ,Ki(",Kilabel,") = ",signif(KiAB, digits=3),sep=""),cex=0.8)

# generate a blank plot (1/sub2 - 1/v) and then plot the raw data
dev.new()
plot(dr.df$rconc2,dr.df$rrate, pch="", xlim=c(0,max(dr.df$rconc2)),ylim=c(0,max(dr.df$rrate)),main=paste(Mtitle," (1/",Sub2name," - 1/v)",sep=""),xlab=paste("1 /",Sub2name),ylab="1/v")
for (i in 1:length(Sub1.rconc)) {points(subset(dr.df$rconc2, dr.df$rconc1==Sub1.rconc[i]), subset(dr.df$rrate, dr.df$rconc1==Sub1.rconc[i]), col=Sub1.color[i], pch=i)}

#draw lines
for (i in 1:length(Sub1.rconc)) {  abline(KmA/Vmax*Sub1.rconc[i] + 1/Vmax, KmB/Vmax + KiAB*KmAB/Vmax, col=Sub1.color[i], lty=i) }
abline(v = 0, h = 0)

# write legends
labels <- Sub1.conc
legend("topleft", paste(Sub1name,labels,Sub1unit), col=Sub1.color, pch=c(1:length(Sub1.conc)), bg="white")
text(max(dr.df$rconc2)/2,0,paste("Km(",Sub1name,") = ",signif(KmA,digits=3)," ",Sub1unit,", Km(",Sub2name,") = ",signif(KmB,digits=3)," ",Sub2unit," ,Ki(",Kilabel,") = ",signif(KiAB, digits=3),sep=""),cex=0.8)