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


