########################################################################################
## 6/10/2006 - Gil Voge [gvoge][@][ac-grenoble][fr]
## 4/12/2006 - Remplacement de mean() par median() dans Eisenthal-Cornish-Bowden
## Les données d'apres L. Michaelis and ML Menten. Die Kinetik der Invertinwirkung. Biochemische Zeitschrift, 49 : 333-369, 1913
## Positionner le répertoire de travail avec getwd() et setwd() ou dans le menu "Divers"
## Appeler le fichier avec source("Enzymo")
## Introduire les valeurs expérimentales de la concentration initiale en substrat s <- c(val1, val2, val3 ... valn) et de la vitesse initiale v <- c(val1, val2, ..., valn) ou toute autre méthode ...
## Attention, le point de coordonnées s=0, v=0 génère une erreur de division par zéro
## Appeler la fonction Enzymo() et passer les paramètres s, v, "Nom de l'expérimentateur"
## En cas d'omission d'un des paramètres, les valeurs ci-dessous leurs sont attribuées par défaut :
#########################################################################################
Enzymo <- function (s=c(0.0052, 0.0104, 0.0208, 0.0416, 0.0833, 0.1670, 0.3330), v=c(0.866, 1.466, 2.114, 2.666, 3.236, 3.636, 3.636), nom="Leonor et Maud")

{

	
	# Partition de la fenêtre graphique
	par_defaut<-par()
	on.exit(par(par_defaut))
	par(mfrow=c(2,2))


	
	
	## Points expérimentaux
	plot(1:length(s),1:length(s)+1, type = "n", xlab = "", ylab = "", axes=FALSE)
	title(main=paste("Valeurs exp\n", nom))
	text(2, (length(s)+1):1, paste(s[1:length(s)]))
	text(length(s)-2, (length(s)+1):1, paste(v[1:length(s)]))
	text(1, length(s)+1, paste("s = "))
	text(length(s)-3, length(s)+1, paste("v = "))
		
	## Lineweaver-Burk
	is <- 1/s; iv <- 1/v
	lb <- lm (iv~is)
	a <- summary(lb)$coefficients[2]
	b <- summary(lb)$coefficients[1]
	r2 <- summary(lb)$r.squared
	VmLB<- 1/b; KmLB<- a*VmLB
	plot (x=is, y=iv, xlab="1/[S]", ylab="1/V", xlim=c(-max(is),max(is)), ylim=c(0, max(iv)), main="Lineweaver & Burk\n (doubles inverses)", axes=FALSE)
	lines (abline(lb))
	axis(1, pos=0)
	axis(2, pos=0)
	legend(-max(is),max(iv),c(as.expression(substitute(hat(Y)==pente*X+b,list(pente=round(a,3),b=round(b,3)))),as.expression(substitute(italic(R)^2==r2,list(r2=round(r2,4)))), as.expression(substitute(Km==K,list(K=round(KmLB,4)))), as.expression(substitute(Vm==V,list(V=round(VmLB,4))))), bty="n", xjust=0,yjust=1)
	mtext(expression(-1/K[m]), side=1, line=0, at=-1/KmLB, col="red")
	
	## Eisenthal-Cornish-Bowden
	n <- length(s)
	a <- numeric(n)
	K <- numeric(2*n)
	V <- numeric(2*n)
	plot(s,v, xlim=c(-max(s),max(s)), ylim=c(0, 2*max(v)), main="Eisenthal Cornish Bowden", type="n", xlab="Concentration initiale en substrat [s]", ylab="Vitesse initiale (v)")
	for (i in 1:n) {
		a[i] <- v[i]/s[i]
		abline(v[i], a[i])
	}
	k <-1
	for (i in 1:(n-1)){
		for (j in (i+1):n){
			K[k] <- (v[j]-v[i])/(a[i]-a[j])
			V[k] <- a[i]*K[k]+v[i]
			k <- k+1
		}
	}
	KmCB <- median (K)
	VmCB <- median (V)
	abline(h=VmCB, col="red", lty=2)
	abline(v=KmCB, col="red", lty=2)
	legend(KmCB, min(v)/0.8, bty="n", as.expression(substitute(Km==K,list(K=round(KmCB,4)))))
	legend(-max(s), VmCB/0.75, bty="n", as.expression(substitute(Vm==V,list(V=round(VmCB,4)))))
	
	
	## Ajustement non lineaire
	## start fournit les valeurs initiales pour Km et Vm
	modnl <- nls(v~(Vm*s)/(Km+s), start=list(Km=s[1], Vm=max(v)))
	modnl
	KmNL<-summary(modnl)$parameters[1]
	VmNL <-summary(modnl)$parameters[2]
	plot(s, v, xlab="Concentration initiale en substrat [s]", ylab="Vitesse initiale (v)", main="Ajustement non lineaire")
	x <- seq(from=0, to=max(s), by=max(s)/500, lengh=500)
	y <- (VmNL*x)/(KmNL+x)
	lines(x, y, col="red")
	legend(max(s)/1.5,max(v)/1.2,c(as.expression(substitute(Km==K,list(K=round(KmNL,4)))), as.expression(substitute(Vm==V,list(V=round(VmNL,4))))), bty="n")
}