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))
