Dear R users:
I encountered difficulties in michaelis-menten equation. I found
that when I use right model definiens, I got wrong Km vlaue,
and I got right Km value when i use wrong model definiens.
The value of Vd and Vmax are correct in these two models.
#-----right model definiens--------
PKindex<-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
mm.model <- function(time, y, parms) {
dCpdt <-
-(parms["Vm"]/parms["Vd"])*y[1]/(parms["Km"]+y[1])
list(dCpdt)}
Dose<-300
modfun <- function(time,Vm,Km,Vd) {
out <- lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
rtol=1e-8,atol=1e-8)
out[,2] }
objfun <- function(par) {
out <- modfun(PKindex$time,par[1],par[2],par[3])
sum((PKindex$conc-out)^2) }
fit <- optim(c(10,1,80),objfun, method="Nelder-Mead)
print(fit$par)
[1] 10.0390733 0.1341544 34.9891829 #--Km=0.1341544,wrong value--
#-----wrong model definiens--------
#-----Km should not divided by Vd--
PKindex<-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
mm.model <- function(time, y, parms) {
dCpdt <-
-(parms["Vm"]/parms["Vd"])*y[1]/(parms["Km"]/parms["Vd"]+y[1])
list(dCpdt)}
Dose<-300
modfun <- function(time,Vm,Km,Vd) {
out <- lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
rtol=1e-8,atol=1e-8)
out[,2]
}
objfun <- function(par) {
out <- modfun(PKindex$time,par[1],par[2],par[3])
sum((PKindex$conc-out)^2)}
fit <- optim(c(10,1,80),objfun, method="Nelder-Mead)
print(fit$par)
[1] 10.038821 4.690267 34.989239 #--Km=4.690267,right value--
What did I do wrong, and how to fix it?
Any suggestions would be greatly appreciated.
Thanks in advance!!
"Chun-Ying Lee" <u9370004 at cc.kmu.edu.tw> writes:> Dear R users: > I encountered difficulties in michaelis-menten equation. I found > that when I use right model definiens, I got wrong Km vlaue, > and I got right Km value when i use wrong model definiens. > The value of Vd and Vmax are correct in these two models.How do you know what the correct value is? Are you sure that the other values are right? I'm a bit rusty on MM, but are you sure your "right" model is right? Try doing a dimensional analysis on the ODE. I kind of suspect that Vd is entering in the wrong way. Since you're dealing in concentrations, should it enter at all (except via the conc. at time 0, of course)? Not knowing the context, I can't be quite sure, but generally, I'd expect Vm*Km/(Km+y) to be the reaction rate, so that Vm is the maximum rate, attained when y is zero and Km is the conc. at half-maximum rate. This doesn't look quit like what you have.> #-----right model definiens-------- > PKindex<-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24), > conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89)) > mm.model <- function(time, y, parms) { > dCpdt <- -(parms["Vm"]/parms["Vd"])*y[1]/(parms["Km"]+y[1]) > list(dCpdt)} > Dose<-300 > modfun <- function(time,Vm,Km,Vd) { > out <- lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd), > rtol=1e-8,atol=1e-8) > out[,2] } > objfun <- function(par) { > out <- modfun(PKindex$time,par[1],par[2],par[3]) > sum((PKindex$conc-out)^2) } > fit <- optim(c(10,1,80),objfun, method="Nelder-Mead) > print(fit$par) > [1] 10.0390733 0.1341544 34.9891829 #--Km=0.1341544,wrong value-- > > > #-----wrong model definiens-------- > #-----Km should not divided by Vd-- > PKindex<-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24), > conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89)) > mm.model <- function(time, y, parms) { > dCpdt <- -(parms["Vm"]/parms["Vd"])*y[1]/(parms["Km"]/parms["Vd"]+y[1]) > list(dCpdt)} > Dose<-300 > modfun <- function(time,Vm,Km,Vd) { > out <- lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd), > rtol=1e-8,atol=1e-8) > out[,2] > } > objfun <- function(par) { > out <- modfun(PKindex$time,par[1],par[2],par[3]) > sum((PKindex$conc-out)^2)} > fit <- optim(c(10,1,80),objfun, method="Nelder-Mead) > print(fit$par) > [1] 10.038821 4.690267 34.989239 #--Km=4.690267,right value-- > > What did I do wrong, and how to fix it? > Any suggestions would be greatly appreciated. > Thanks in advance!! > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html-- O__ ---- Peter Dalgaard ??ster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
Chun-Ying Lee wrote:> Dear R users: > I encountered difficulties in michaelis-menten equation. I found > that when I use right model definiens, I got wrong Km vlaue, > and I got right Km value when i use wrong model definiens. > The value of Vd and Vmax are correct in these two models. > > #-----right model definiens-------- > PKindex<-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24), > conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89)) > mm.model <- function(time, y, parms) { > dCpdt <- -(parms["Vm"]/parms["Vd"])*y[1]/(parms["Km"]+y[1]) > list(dCpdt)} > Dose<-300 > modfun <- function(time,Vm,Km,Vd) { > out <- lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd), > rtol=1e-8,atol=1e-8) > out[,2] } > objfun <- function(par) { > out <- modfun(PKindex$time,par[1],par[2],par[3]) > sum((PKindex$conc-out)^2) } > fit <- optim(c(10,1,80),objfun, method="Nelder-Mead) > print(fit$par) > [1] 10.0390733 0.1341544 34.9891829 #--Km=0.1341544,wrong value-- > > > #-----wrong model definiens-------- > #-----Km should not divided by Vd-- > PKindex<-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24), > conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89)) > mm.model <- function(time, y, parms) { > dCpdt <- -(parms["Vm"]/parms["Vd"])*y[1]/(parms["Km"]/parms["Vd"]+y[1]) > list(dCpdt)} > Dose<-300 > modfun <- function(time,Vm,Km,Vd) { > out <- lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd), > rtol=1e-8,atol=1e-8) > out[,2] > } > objfun <- function(par) { > out <- modfun(PKindex$time,par[1],par[2],par[3]) > sum((PKindex$conc-out)^2)} > fit <- optim(c(10,1,80),objfun, method="Nelder-Mead) > print(fit$par) > [1] 10.038821 4.690267 34.989239 #--Km=4.690267,right value-- > > What did I do wrong, and how to fix it? > Any suggestions would be greatly appreciated. > Thanks in advance!! > > >it is not clear to me what you are trying to do: you seem to have a time-concentration-curve in PKindex and you seem to set up a derivative of this time dependency according to some model in dCpdt. AFAIKS this scenario is not directly related to the situation described by the Michaelis-Menten-Equation which relates some "input" concentration with some "product" concentration. If Vm and Km are meant to be the canonical symbols, what is Vd, a volume of distribution? it is impossible to see (at least for me) what exactly you want to achieve. (and in any case, I would prefer "nls" for a least squares fit instead of 'optim'). joerg> ------------------------------------------------------------------------ > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html