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@manchester.ac.uk
[[alternative HTML version deleted]]
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))
>
>
>
>
>
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[ ... long example deleted, see original post ...] Hi In-Sun, as far as I see the *first* of your problems is your high demand on precision. Double precision allows 16 digits, so 1e-10 may lead internally to values that are close to the maximum number of digits in double precision arithmetics, because you set relative *and* absolute tolerance of all states to such a small value. Reduce atol and rtol to reasonable values (e.g. the default of 1e-6) and it should work (it does on my machine). Another possibility would be to increase maxsteps, but this will not help if the requested precision is higher than would be possible with double precision. If you try the following: ODE <- as.data.frame(daspk(y = y, times = times, func = Fun_ODE, parms = pars, atol = 1e-6, rtol = 1e-6)) ... then the simulation goes through in a technical sense, however state variables reach very high values which, again, are beyond what is possible with double precision arithmetics -- is this what you mean with "the numbers are everywhere"? Are you sure that your model formulation is correct? Please check your equations, especially parentheses and signs. Thomas Petzoldt
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.