setwd('...')

### Importation des données

dta = read.table("donnees.txt",header=TRUE,sep="\t",dec=",")
summary(dta)

### Visualisation des cinétiques

tps = 0:14

plot(tps,dta$DO[dta$Binome==1],type="b",pch=16,col="blue",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,xlab="temps (min)",
   ylab="Densité optique",ylim=c(0,3),bty="l")
points(tps,dta$DO[dta$Binome==2],type="b",col="red",pch=16,lwd=2,cex=1.25)
points(tps,dta$DO[dta$Binome==3],type="b",col="green",pch=16,lwd=2,cex=1.25)
points(tps,dta$DO[dta$Binome==4],type="b",col="purple",pch=16,lwd=2,cex=1.25)
points(tps,dta$DO[dta$Binome==5],type="b",col="brown",pch=16,lwd=2,cex=1.25)
legend('bottomright',lty=rep(1,5),lwd=2,pch=rep(16,5),cex=1.25,bty="n",legend=paste('Expérience',1:5,sep=" "),
   col=c("blue","red","green","purple","brown"))

### Estimations des vitesses initiales

delta = diff(dta$DO)
vinit = delta[seq(from=1,to=length(delta),by=15)]
dtavinit = data.frame(Binome=1:5,Vinit=vinit,S0=dta$S[dta$Temps=="T0"],E0=dta$E[dta$Temps=="T0"])

### Représentation de Lineweaver et Burk

plot(1/dtavinit$S0,1/dtavinit$Vinit,type="p",pch=16,col="blue",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab=expression(1/s[0]),ylab=expression(1/hat(v)[0]),ylim=c(0,0.65),bty="l",xlim=c(-4,1))

mod = lm(I(1/Vinit)~I(1/S0),data=dtavinit)

abline(coef(mod),lwd=2,col="blue")
abline(h=0,lwd=2,col="darkgray")
abline(v=0,lwd=2,col="darkgray")

a = coef(mod)[1]
b = coef(mod)[2]

text(-a/b,0.05,expression(-over(1,K[m])),cex=1.5,col="red",adj=1)
text(0,1.1*a,expression(over(1,k[3]*e[0])),cex=1.5,col="red",adj=1)

text(-3,0.6,expression(hat(K)[m]==0.283),cex=1.5,col="blue",adj=1)
text(-3,0.6,expression(hat(K)[m]==0.283),cex=1.5,col="blue",adj=1)

### Représentation d'Eadie

plot(dtavinit$Vinit[-1]/dtavinit$S0[-1],dtavinit$Vinit[-1],type="p",pch=16,col="blue",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab=expression(hat(v)[0]/s[0]),ylab=expression(hat(v)[0]),ylim=c(0,3),bty="l",xlim=c(0,10),main="Diagramme d'Eadie")

mod = lm(Vinit~I(Vinit/S0),data=dtavinit[-1,])

abline(coef(mod),lwd=2,col="blue")

a = coef(mod)[1]
b = coef(mod)[2]

text(-a/b,0.1,expression(over(k[3]*e[0],K[m])),cex=1.5,col="red",adj=-1)
text(0,1.1*a,expression(k[3]*e[0]),cex=1.5,col="red",adj=0)

text(6,3,expression(hat(K)[m]==0.256),cex=1.5,col="blue",adj=1)
text(6,2.75,expression(hat(k)[3]==2.09),cex=1.5,col="blue",adj=1)
  
### Validation du modèle de Michaelis-Menten

Km = 0.256
Vmax = 2.09

f = function(x,Vmax,Km,s0,tt) {
   x-s0+Km*log(x/s0)+Vmax*tt
}

tps = seq(0,14,0.1)
cins = matrix(0,nrow=length(tps),ncol=4)
for (i in 1:length(tps)) {
   cins[i,1] = uniroot(f,Vmax=Vmax,Km=Km,s0=1,tt=tps[i],interval=c(0,50))$root
   cins[i,2] = uniroot(f,Vmax=Vmax,Km=Km,s0=3,tt=tps[i],interval=c(0,50))$root
   cins[i,3] = uniroot(f,Vmax=Vmax,Km=Km,s0=10,tt=tps[i],interval=c(0,50))$root
   cins[i,4] = uniroot(f,Vmax=Vmax,Km=Km,s0=30,tt=tps[i],interval=c(0,50))$root
}

cinv = Vmax*cins/(Km+cins)

library(caTools)

resp = matrix(0,nrow=length(tps),ncol=4)
for (i in 2:length(tps)) {
   resp[i,1] = dta$DO[(dta$Binome==5)&(dta$Temps=="T0")]+
      dta$DO[(dta$Binome==5)&(dta$Temps=="T14")]*trapz(tps[1:i],cinv[1:i,1])/trapz(tps[1:length(tps)],cinv[1:length(tps),1])
   resp[i,2] = dta$DO[(dta$Binome==4)&(dta$Temps=="T0")]+
      dta$DO[(dta$Binome==4)&(dta$Temps=="T14")]*trapz(tps[1:i],cinv[1:i,2])/trapz(tps[1:length(tps)],cinv[1:length(tps),2])
   resp[i,3] = dta$DO[(dta$Binome==3)&(dta$Temps=="T0")]+
      dta$DO[(dta$Binome==3)&(dta$Temps=="T14")]*trapz(tps[1:i],cinv[1:i,3])/trapz(tps[1:length(tps)],cinv[1:length(tps),3])
   resp[i,4] = dta$DO[(dta$Binome==2)&(dta$Temps=="T0")]+
      dta$DO[(dta$Binome==2)&(dta$Temps=="T14")]*trapz(tps[1:i],cinv[1:i,4])/trapz(tps[1:length(tps)],cinv[1:length(tps),4])
}

plot(tps,resp[,1],type="l",col="brown",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab="Temps (min)",ylab=expression(p[t]),ylim=c(0,4),bty="l",xlim=c(0,14),main="Modèle de catalyse")
lines(tps,resp[,2],col="purple",lwd=2)
lines(tps,resp[,3],col="green",lwd=2)
lines(tps,resp[,4],col="red",lwd=2)
legend('bottomright',lty=rep(1,4),lwd=2,cex=1.25,bty="n",legend=paste('S0=',c(1,3,10,30),sep=" "),
   col=c("brown","purple","green","red"))

savePlot("pic/modelp.pdf",type="pdf")

plot(tps,resp[,1],type="l",col="brown",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab="Temps (min)",ylab=expression(p[t]),ylim=c(0,4),bty="l",xlim=c(0,14),main="Modèle de catalyse - S0 = 1")
points(0:14,dta$DO[(dta$Binome==5)],type="p",col="brown",pch=16,cex=1.25)

savePlot("pic/modelp1.pdf",type="pdf")

plot(tps,resp[,2],type="l",col="purple",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab="Temps (min)",ylab=expression(p[t]),ylim=c(0,4),bty="l",xlim=c(0,14),main="Modèle de catalyse - S0 = 3")
points(0:14,dta$DO[(dta$Binome==4)],col="purple",pch=16,cex=1.25)

savePlot("pic/modelp3.pdf",type="pdf")

plot(tps,resp[,3],type="l",col="green",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab="Temps (min)",ylab=expression(p[t]),ylim=c(0,4),bty="l",xlim=c(0,14),main="Modèle de catalyse - S0 = 10")
points(0:14,dta$DO[(dta$Binome==3)],col="green",pch=16,cex=1.25)

savePlot("pic/modelp10.pdf",type="pdf")

plot(tps,resp[,4],type="l",col="red",lwd=2,cex=1.25,cex.lab=1.25,cex.axis=1.25,
   xlab="Temps (min)",ylab=expression(p[t]),ylim=c(0,4),bty="l",xlim=c(0,14),main="Modèle de catalyse - S0 = 30")
points(0:14,dta$DO[(dta$Binome==2)],col="red",pch=16,cex=1.25)

savePlot("pic/modelp30.pdf",type="pdf")