# Kinetic fitting of Ping-Pong mechanism
# Created by S. Fushinobu (2018.11.28)

# Edit the parameters below
pingpong.df <- read.csv("/Users/fushi/Downloads/PP-R/rate-pp.csv")
Mtitle <- "Ping-Pong 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
Vmaxinit <- 2.7
grid <- 10                 # grid number for best fit plot lines
mono <- FALSE              # FALSE for color, TRUE for monochrome

# determine concentrations of Substrates 1 & 2 used in experiment
Sub1.conc <- unique(pingpong.df$conc1)
Sub2.conc <- unique(pingpong.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))
}

# generate a blank plot (sub1 - v) and then plot the raw data
plot(pingpong.df$conc1,pingpong.df$rate, pch="", xlim=c(0,max(pingpong.df$conc1)),ylim=c(0,max(pingpong.df$rate)),main=paste(Mtitle,"(",Sub1name,"- v )"),xlab=paste(Sub1name," (",Sub1unit,")",sep=""),ylab=paste("v",Vunit))
for (i in 1:length(Sub2.conc)) {points(subset(pingpong.df$conc1, pingpong.df$conc2==Sub2.conc[i]), subset(pingpong.df$rate, pingpong.df$conc2==Sub2.conc[i]), col=Sub2.color[i], pch=i)}

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

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

# confidence intervals of parameters
confint(pingpong.nls)

# extract coefficients
KmA <- unname(coef(pingpong.nls)["KmA"])
Vmax <- unname(coef(pingpong.nls)["Vmax"])
KmB <- unname(coef(pingpong.nls)["KmB"])

# use the values directly in the equation to calculate line of best fit
gridw <- max(pingpong.df$conc1)*grid
fit.data <- expand.grid(x=(0:gridw)/grid, Sub2=Sub2.conc)
fit.data$y <- (Vmax*fit.data$x*fit.data$Sub2) / (KmA*fit.data$Sub2 + KmB*fit.data$x + fit.data$x*fit.data$Sub2)

# plot lines of best fit
for (i in 1:length(Sub2.conc)) {lines(subset(fit.data$x, fit.data$Sub2==Sub2.conc[i]), subset(fit.data$y, fit.data$Sub2==Sub2.conc[i]), lty=i, col=Sub2.color[i])}

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

# generate a blank plot (sub2 - v) and then plot the raw data
dev.new()
plot(pingpong.df$conc2,pingpong.df$rate, pch="", xlim=c(0,max(pingpong.df$conc2)),ylim=c(0,max(pingpong.df$rate)),main=paste(Mtitle,"(",Sub2name,"- v )"),xlab=paste(Sub2name," (",Sub2unit,")",sep=""),ylab=paste("v",Vunit))
for (i in 1:length(Sub1.conc)) {points(subset(pingpong.df$conc2, pingpong.df$conc1==Sub1.conc[i]), subset(pingpong.df$rate, pingpong.df$conc1==Sub1.conc[i]), col=Sub1.color[i], pch=i)}

# use the values directly in the equation to calculate line of best fit
gridw <- max(pingpong.df$conc2)*grid
fit.data <- expand.grid(x=(0:gridw)/grid, Sub1=Sub1.conc)
fit.data$y <- (Vmax*fit.data$x*fit.data$Sub1) / (KmB*fit.data$Sub1 + KmA*fit.data$x + fit.data$x*fit.data$Sub1)

# plot lines of best fit
for (i in 1:length(Sub1.conc)) {lines(subset(fit.data$x, fit.data$Sub1==Sub1.conc[i]), subset(fit.data$y, fit.data$Sub1==Sub1.conc[i]), lty=i, col=Sub1.color[i])}

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

#prepare for double-reciplocal (Lineweaver-Burk) plot
dr.df <- data.frame(rconc1=1/pingpong.df$conc1, rrate=1/pingpong.df$rate, rconc2=1/pingpong.df$conc2, conc1=pingpong.df$conc1, conc2=pingpong.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
dev.new()
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, 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")

# 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, 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")