Re: [R] Michaelis-menten equation

2005-07-20 Thread Peter Dalgaard
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

2005-07-20 Thread joerg van den hoff
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

2005-07-19 Thread Chun-Ying Lee
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

2005-07-19 Thread Peter Dalgaard
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

2005-07-19 Thread joerg van den hoff
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

2005-07-19 Thread Peter Dalgaard
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

2005-07-19 Thread Chun-Ying Lee
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