Re: [R] Michaelis-menten equation
Chun-Ying Lee [EMAIL PROTECTED] writes: Hi, We used known Vm and Km to simulate the data set (time, Cp) without adding random error in there. Yes, the line looks like very close to a straight line. But why can't we obtain the correct values with fitting process? We used optim first and then followed by using nls to fit the model. Thanks. regards, Check your simulation. There is no way that curve is consistent with Km in the middle of the y range! ---Chun-ying Lee Hmm, sorry, no. I'm talking through a hole in my head there. Vm*y/(Km+y) makes OK sense. Vm is what you get for large y - passing from 1st order to 0th order kinetics. However, looking at the data plot(PKindex) abline(lm(conc~time,data=PKindex)) shows that they are pretty much on a straight line, i.e. you are in the domain of 0-order kinetics. So why are you expecting the rate of decrease to have changed by roughly 3/4 (from 2/3*Vm/Vd at y=2*Km to 1/2*Vm/Vd at y=Km when you reach 4.67)?? __ R-help@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 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@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
Re: [R] Michaelis-menten equation
I believe the following is correct: 1. first of all, as peter daalgaard already pointed out, your data Cp(t) are following a straight line very closely, i.e. 0.-order kinetics 2. for your diff. eq. this means that you are permanently in the range cp Km so that dCp/dt = - Vm/Vd = const. =: -b and, therefore, Cp = a - b*t 3. you can't get any reliable information concerning Km from the fit. the solution of the diff. eq (according to Maxima), namely t + const. = -(Km*log(Cp) + Cp)/(Vm/Vd) tells you the same: Km*log(cp) Cp in your data. 4. in any case (even in the range Km ~= Cp) you can't determine Vm _and_ Vd separately according to your diff. eq., you only get the ratio b = Vm/Vd. this does make sense: what you are measuring is the decreasing plasma concentration, you don't have any information concerning the relevant volume fraction, i.e. the Vd, in you data. therefore any variation in the effective max. velocity can either be ascribed to a variation of Vm or to a modified Vd. in other words: you should view Vm* = Vm/Vd as your effective Vm. 5. x-PKindex[,1] y-PKindex[,2] res - lm ( y~x ) yields a=8.561, b=Vm*= 0.279 summary(abs(residuals(r)/y)*100) Min. 1st Qu. MedianMean 3rd Qu.Max. 0.00863 0.09340 0.13700 0.26500 0.22900 1.37000 i.e. 0.265 percent deviation on average between data and fit. I believe that is the max. information you can get from your data ((the result for b is accidently is not so far away from the ratio of your results Vm = 10.04, Vd = 34.99 which actually must be completely unstable (double both parameters and nothing happens to the fit). 6. you state you simulated the data with km=4.69? using the above Vm and Vd, the resulting data are (not unexpectedly) quite different from those you used as input to the fit. maybe you made an error somewhere upstream? 7. in conclusion: don't try to fit Vm _and_ Vd separately, check whether your simulated data are correct, keep in mind that if kmCp, you can't fit Km (at least not reliably). Chun-Ying Lee wrote: Hi, We are doing a pharmaockinetic modeling. This model is described as intravenous injection of a certain drug into the body. Then the drug molecule will be eliminated (or decayed) from the body. We here used a MM eq. to describe the elimination process and the changes of the drug conc.. So the diff. eq. would be: dCp/dt = -Vm/Vd * Cp/(Km+Cp). Vd is a volume of distribution. We used lsoda to solve the diff. eq. first and fit the diff. eq. with optim first (Nelder-Mead simplex) and followed by using nls to take over the fitting process of optim. However, we can not obtain the correct value for Km if we used the above model. The correct Km can be obtained only when we modeled the diff eq. with dCp/dt= -Vm/Vd * Cp/(Km/vd + Cp). Now we lost. The data were from simulation with known Vm and Km. Any idea? Thanks. regards, --- Chun-ying Lee 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@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 __ R-help@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 __ R-help@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
[R] Michaelis-menten equation
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!! __ R-help@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
Re: [R] Michaelis-menten equation
Chun-Ying Lee [EMAIL PROTECTED] 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@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 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@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
Re: [R] Michaelis-menten equation
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@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 __ R-help@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
Re: [R] Michaelis-menten equation
Peter Dalgaard [EMAIL PROTECTED] writes: Chun-Ying Lee [EMAIL PROTECTED] 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. Hmm, sorry, no. I'm talking through a hole in my head there. Vm*y/(Km+y) makes OK sense. Vm is what you get for large y - passing from 1st order to 0th order kinetics. However, looking at the data plot(PKindex) abline(lm(conc~time,data=PKindex)) shows that they are pretty much on a straight line, i.e. you are in the domain of 0-order kinetics. So why are you expecting the rate of decrease to have changed by roughly 3/4 (from 2/3*Vm/Vd at y=2*Km to 1/2*Vm/Vd at y=Km when you reach 4.67)?? #-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@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 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@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 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@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
Re: [R] Michaelis-menten equation
Hi, We used known Vm and Km to simulate the data set (time, Cp) without adding random error in there. Yes, the line looks like very close to a straight line. But why can't we obtain the correct values with fitting process? We used optim first and then followed by using nls to fit the model. Thanks. regards, ---Chun-ying Lee Hmm, sorry, no. I'm talking through a hole in my head there. Vm*y/(Km+y) makes OK sense. Vm is what you get for large y - passing from 1st order to 0th order kinetics. However, looking at the data plot(PKindex) abline(lm(conc~time,data=PKindex)) shows that they are pretty much on a straight line, i.e. you are in the domain of 0-order kinetics. So why are you expecting the rate of decrease to have changed by roughly 3/4 (from 2/3*Vm/Vd at y=2*Km to 1/2*Vm/Vd at y=Km when you reach 4.67)?? __ R-help@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