[R] deSolve question

[spencerg]

Dear In-Sun: 


      I have not seen a reply, so I will offer some suggestions, hoping 
I can help without understanding all the details. 


           1.  Have you run that code with "options(warn=2)"?  It 
produced over 50 warnings for me, and "options(warn=2)" will convert the 
warnings to errors, making it easier for you to find the exact problem 
lines of code.  With this, have you used the "debug" function to trace 
through your code line by line to provide more detail about what the 
code is doing?  I use that routinely. 


           2.  Have you considered the "PKfit" package?  I don't know if 
it includes your model, but it includes code for many pharmacokinetic 
models.  If your interest in "deSolve" is as a means to solve this 
problem, you might consider "PKfit". 


           3.  Have you tried writing directly to the authors?  Names 
and email addresses are available in "help(pac=deSolve)". 


      Hope this helps. 
      Spencer Graves


insun nam wrote:
> Dear All,
>
> I like to simulate a physiologically based pharmacokinetics model using R
> but am having a problem with the daspk routine.
>
> The same problem has been implemented in Berkeley madonna and Winbugs so
> that I know that it is working. However, with daspk it is not, and the
> numbers are everywhere!
>
> Please see the following and let me know if I am missing something...
>
> Thanks a lot in advance,
> In-Sun
>
> #---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
> library("deSolve")
>
> y <- c(Agi = 0,Alu = 0, Abr  = 0, Ah = 0, Ali = 0,  Ak = 0,  Am = 0,  Ask =
> 0,  Aad  = 0,  Apa  = 0, Asp  = 0, Aar  = 0,  Ave  = 0)
> times = seq(0, 100, length=100)
>
> pars <- c(
> dose = 80 * 0.26,
> doseduration = 10,
> Vmax = 1.44,
> Km = 0.96,
> s = 1.33,
> fp = 0.236,
> Kpfgi=0.324,
> Kpflu = 1.092,
> Kpfbr= 0.155 ,
> Kpfh=0.767,
> Kpfli = 0.551,
> Kpfk=0.537,
> Kpfm=0.339,
> Kpfsk=0.784,
> Kpfad=0.465,
> Kpfpa=0.595,
> Kpfsp=0.410,
> Qar = 51.9,
> Qve = 51.9,
> Qgi = 12.3,
> Qlu = 51.9,
> Qbr = 3.2,
> Qh = 6.4,
> Qli = 16.5,
> Qk = 12.8,
> Qm = 7.6,
> Qsk = 5.0,
> Qad = 0.4,
> Qpa = 1.0,
> Qsp = 1.0,
> Var = 7.0,
> Vve = 14.1,
> Vgi = 12.4,
> Vlu = 1.3,
> Vbr = 1.3,
> Vh = 1.2,
> Vli = 12.4,
> Vk = 2.2,
> Vm = 140.0,
> Vsk = 49.0,
> Vad = 11.2,
> Vpa = 1.0,
> Vsp = 1.0
> )
>
> Fun_ODE <- function(t,y, pars){
> with (as.list(c(y, pars)), {
> It <- dose/doseduration
> Car <- Aar/Var
> Cve <- Ave/Vve
> Clu <- Alu/Vlu
> Cli <- Ali/Vli
> Cbr <- Abr/Vbr
> Ch <- Ah/Vh
> Cpa <- Apa/Vpa
> Csp <- Asp/Vsp
> Cgi <- Agi/Vgi
> Ck <- Ak/Vk
> Cm <- Am/Vm
> Cad <- Aad/Vad
> Csk <- Ask/Vsk
>
> Kpbbr <- s*fp*Kpfbr
> Kpbli <- s*fp*Kpfli
> Kpbh <- s*fp*Kpfh
> Kpbpa <- s*fp*Kpfpa
> Kpbsp <- s*fp*Kpfsp
> Kpbgi <- s*fp*Kpfgi
> Kpbk <- s*fp*Kpfk
> Kpbm <- s*fp*Kpfm
> Kpbad <- s*fp*Kpfad
> Kpbsk <- s*fp*Kpfsk
> Kpblu <- s*fp*Kpflu
>
> dAar <- (Clu/Kpblu - Car)*Qar
> dAve <- if (t < 10) It + Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli +
> Ck*Qk/Kpbk + Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
>
>         else Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli + Ck*Qk/Kpbk +
> Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
> dAlu <- (Cve-Clu/Kpblu)*Qlu
>
> dAli <- ((Qli - Qgi- Qpa-Qsp)*Car + Cgi*Qgi/Kpbgi + Csp*Qsp/Kpbsp +
> Cpa*Qpa/Kpbpa - Cli*Qli/Kpbli) - Vmax*Cli/Kpfli/(Km + Cli/Kpfli)
> dAbr <- (Car - Cbr/Kpbbr)*Qbr
> dAh <- (Car - Ch/Kpbh)*Qh
> dApa <- (Car - Cpa/Kpbpa)*Qpa
> dAsp <- (Car - Csp/Kpbsp)*Qsp
> dAgi <- (Car - Cgi/Kpbgi)*Qgi
> dAk <- (Car - Ck/Kpbk)*Qk
> dAm <- (Car - Cm/Kpbm)*Qm
> dAad <- (Car - Cad/Kpbad)*Qad
> dAsk <- (Car - Csk/Kpbsk)*Qsk
>
> return(list(dy = c(dAar, dAve, dAlu, dAli, dAbr, dAh, dApa, dAsp, dAgi, dAk,
> dAm, dAad, dAsk),
>             Car = Car, Cve=Cve, Clu=Clu, Cli=Cli, Cbr=Cbr, Ch=Ch, Cpa=Cpa,
> Csp=Csp, Cgi=Cgi, Ck=Ck, Cm=Cm, Cad=Cad, Csk=Csk))
> })
> }
>
> ODE <- as.data.frame(daspk(y = y, times = times, func = Fun_ODE,
> parms = pars, atol = 1e-10, rtol = 1e-10))
>
>
>
>
>