R/SでODE - Clinical PK/PD

R で ODE

ex1 <-
function(t, y, p) {
  dy1 <- -p["ka"] * y[1] # amount
  dy2 <-  p["ka"] * y[1] / p["V"] - p["CL"] / p["V"] * y[2] # concentration

  list(c(dy1, dy2)) # list で返す
}

ka <- 0.7; CL <- 0.1; V <- 1
parms <- c(ka=ka, CL=CL, V=V)

TIME <- seq(0, 24, by=1)
Dose <- 1000

require(odesolve)

my.atol <- c(1e-6, 1e-10)
out <-
lsoda(c(Dose, 0), TIME, ex1, parms, rtol=1e-5, atol=my.atol)

out

S-PLUS で ODE

ex1 <-
function(x, y, CL, V, ka) {
  dy1 <- -ka * y[1] # amount
  dy2 <-  ka * y[1] / V - CL * y[2] # concentration
  c(dy1, dy2)
}

TIME <- seq(1, 24, by=1)
CONC <- double(length(TIME))

CL <- 0.1; V <- 1; ka <- 0.7
Dose <- 1000
out <- ivp.ab(TIME[1], c(0, c(Dose, 0)), ex1, aux=list(CL=CL, V=V, ka=ka))
CONC[1] <- out$values[3]

for (i in seq(2, length(TIME))) {
  out <- ivp.ab(TIME[i], restart=out)
  CONC[i] <- out$values[3]
}

CONC