[R-sig-dyn-mod] a question on ODE.1D, nspec

Hi all:

         I am kind of confused by the ODE.1D function in deSolve.

        I used deSolve to simulate the oral absorption of drug with
gastric emptying.

         I divided the human digestive tract into 100 sections, and
simulate the drug travel through and get absorbed in the same time.
         in the state function, the first 100 elements are the drug amount
in each section, the 101th section is the drug absorbed into blood (with
elemination).
          My question is on "nspec". It was fine when I set "nspec" as 1
and names as mass. But I think there are multiple species in the
simulation. When I set "nspec" as 2 and names as "mass" and "plasma", the
system won't work.
          I am just wondering how "nspec" works inside deSolve. Is it an
important parameter to set?

          Regards!
Dabing




rm(list=ls())
library (deSolve)


# simulate the plug flow of solution dosing with gastric emptying
N       <- 100   # divide the GIT into 100 segments
dx      <- 600/N        # grid size of 6 cm, d, total length of GIT 600 cm
v       <- 600/(4*60)       # velocity of mass movement, cm/min
x       <- seq(dx/2, 600, by = dx)  # m, distance from river

ke       <- 0.693/60      # /min, first-order in vivo elemination rate from
plasma
mass.0 <- 100   # initial dose in mg into GIT


state <- c(100,rep(0, N-1),0)

plug <- function(t, state, paras) {

 mass <- state[1:N]
 plasma <- state[N+1]
 ka <- c(0,rep(0.02,N-1))

 conc <- mass / dx  # conc of drug in mg/cm
 ## BOD dynamics
 Flux <- v * c(0, conc)   # fluxes due to water transport
 Flux[2] <- 0.01*mass[1] +0.2*mass[1]*(tanh(t%%240-150)-tanh(t%%240-180))

 absorption <- ka*mass             # oral absorption

 ## rate of change = plug flow movement  - absorption
 dmass        <- -diff(Flux) - absorption
 dplasma <- sum(absorption)   -ke*plasma


 return(list(c(dmass,dplasma)))
}


times <- seq(0, 1000, by = 1)


systime <- system.time(out <- ode.1D(y = state, times, plug, parms = NULL,
             nspec = 2, names = c("mass","plasma")))


#matplot.1D(out, type = "l", grid = x)
#image(out, grid = x,  xlab = "time, min", ylab = "GIT")
color = topo.colors
drug <- out[,2:(N+1)]
filled.contour(x = times, y = x, drug, color = color, nlevels = 10,
              xlab = "time, min", ylab = "GIT distance, cm",main =
"drug",ylim=c(0,20), xlim=c(0,200))

#matplot.1D(out, type = "l", grid = x,subset = time == 500)
#matplot(out[,1],out[,102],type="l",log="y")
matplot(out[,1],out[,N+2],type="l")
systime

	[[alternative HTML version deleted]]