[R] deSolve question

[Karline]

In-Sun, 

It is very simple. You define your state variables in the following order:
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)

and your rates of change in another order:
dy = c(dAar, dAve, dAlu, dAli, dAbr, dAh, dApa, dAsp, dAgi, dAk,
dAm, dAad, dAsk)

So, The solver uses the rate of change of Aar to update Agi etc... and you
get nonsense in the end.

I hope correcting this will solve your problem. 

I have seen this multiple times; it is definitely the most common type of
mistake in using R for solving differential equations.

Btw, why use daspk ? It is really meant for solving DAEs, not ODEs. 
lsoda or lsode or vode might be a better choice.

Cheers,
Karline


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))
> 
> 
> 
> 
> -- 
> Dr In-Sun Nam Knutsson
> Research Associate
> The Centre for Applied Pharmacokinetic Research (CAPKR)
> School of Pharmacy and Pharmaceutical Sciences
> University of Manchester
> Stopford Building
> Oxford Road
> Manchester
> U.K.
> Phone: +44 161 275 2355
> Email: In-sunnam.Knutsson at manchester.ac.uk
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: http://www.nabble.com/deSolve-question-tp23985008p23994033.html
Sent from the R help mailing list archive at Nabble.com.