# Quadratic equation fit
# Written by S. Fushinobu (2019.10.25)
# Thanks to Carly Huitema & Geoff P Horsman
# Edit the parameters below
Quadra.df <- read.csv("/Users/fushi/Documents/R/quadra.csv")
Mtitle <- "Quadratic fitting of spectroscopic titration absorbance change"  # main title of the plot
Ligname <- "Ligand"
Ligunit <- " (μM)"
Vunit <- " dAbs (387-420 nm)"
Econc <- 5.1              # Protein concentration with the same unit of the ligand
Kdinit <- 0.7               # initial values for fitting
Bmaxinit <- 0.25

plot(Quadra.df$conc,Quadra.df$dA,main=Mtitle,xlim=c(0,max(Quadra.df$conc)),ylim=c(0,max(Quadra.df$dA)),xlab=paste(Ligname,Ligunit),ylab=paste(Vunit))

Quadra.nls <- nls(dA ~ (Bmax/2/Econc)*((Kd + Econc + conc) - sqrt((Kd + Econc + conc)^2 - (4*Econc*conc))), data=Quadra.df, start=list(Kd=Kdinit, Bmax=Bmaxinit))

Kd <- unname(coef(Quadra.nls)["Kd"])
Bmax <- unname(coef(Quadra.nls)["Bmax"])
curve((Bmax/2/Econc)*((Kd+Econc+x)-sqrt((Kd+Econc+x)^2-(4*Econc*x))),0,max(Quadra.df$conc),add=TRUE, lty=2)

text(0,max(Quadra.df$dA),paste("Kd = ",signif(Kd,digits=3),"±",signif(summary(Quadra.nls)$parameters[1,2],digits=3),Ligunit,"\nBmax = ",signif(Bmax,digits=3),"±",signif(summary(Quadra.nls)$parameters[2,2],digits=3)),adj=0)

text(max(Quadra.df$conc)/20,0,expression(paste(Delta, italic(A) == "(",italic(B)[max]/2,italic(E),")","{","(",italic(K)[d]+italic(E)+italic(L),")"-"{(",italic(K)[d]+italic(E)+italic(L), ")"^"2" - "4",italic(EL),"}"^"1/2","}"   )), adj=0)