Hi,
I am doing a project on the simulation of glucose metabolism based on a
pharmacokinetic modeling in which we have 4 differential equations. I did
this in R by using the odesolve package. It works very well, but I have two
questions:
Here is the odemodel function
_________________________________________________
Ogtt.Odemodel <- function(t, y, p) {
absx <- c(-60, -45, -30, -15, 0, 5, 10, 15, 20, 25, 30, 45, 60, 75, 90,
105, 120, 140, 160, 180, 210, 240, 270, 301)
absx <- absx + 60
absy <- c(0, 0, 0, 0, 0, 156.1225165, 266.8509934, 432.5794702,
609.8807946, 759.0364238, 824.7649004, 636.0562915, 482.7450332,
395.1870861, 329.1953641, 591.2408942, 214.6274834, 180.316, 123.3885,
111.0497, 88.83972, 53.30383, 21.32153, 0)
sx <- c(-60, -37.5, -22.5, -7.5, 0, 2.5, 7.5, 12.5, 17.5, 22.5, 27.5,
37.5, 52.5, 67.5, 82.5, 97.5, 112.5, 130, 150, 170, 195, 225, 255, 285, 315)
sx <- sx + 60
sy <- c(81.32, 81.32, 81.32, 81.32, 81.32, 218.24, 329.41, 448.82,
527.06, 588.82, 568.24, 613.53, 564.12, 415.88, 354.12, 325.29, 267.65,
222.35, 115.29, 98.31, 89.81, 85.57, 81.32, 81.32, 81.32)
abs <- approxfun(absx, absy, method = "linear")
s <- approxfun(sx, sy, method = "linear")
Rtb <- 81.32
ib <- ((1-p["f"])*Rtb*p["kg"]*(p["hgo0"]/p["vg"])/(p["k21"] + p["k01"]
+ p["kl"] - p["k21"]*p["k12"]/(p["k12"] + p["k02"]))/(p["ki"]*p["vi"]))
# G1(t)-y(1); G2(t)-y(2);I(t)-y(3); X(t)-y(4)
yd1 <- (abs(t) + p["hgo0"] - (p["kl"] + p["fx"]*y[4])*y[1]*p["vg"] +
p["k12"]*y[2])/p["vg"] - (p["k21"] + p["k01"])*y[1]
yd2 <- p["k21"]*y[1]*p["vg"] - (p["k12"] + p["k02"] +
(1-p["fx"])*y[4])*y[2]
yd3 <- (1-p["f"])*s(t)*p["kg"]*y[1]/p["vi"] - p["ki"]*y[3]
yd4 <- -p["p2"]*y[4] + p["p3"]*(y[3]-ib)
list(c(yd1, yd2, yd3, yd4))
}
__________________________________________________
First is that The computing time for each individual simulation is about 6
sec if I set up time from 0 - 360 min and calculate on each minute. In
order to reduce computing time, I want to calculate for every two minutes,
but I don't know how to do it. I try to make t to be t <- seq(0, 360, by =
2), however, the computing time is still about 6 sec instead of 3. I also
found that there is a hmin parameter in odesolve. But again, it doesn't
work.
The second question is that I found that some time, the last two or three
points of calculation is NA.
Thanks.
Guan
[[alternative HTML version deleted]]
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
