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

[Daniel Reed]

Hi Dabing:

When you set nspec to 2, ode.1D expects an array of length 2*N (i.e., 1:N for the first species, (N+1):(2*N) for the second species). However, your vector is only of length N+1. In other words, ode.1D expects both species to be spatially resolved. I wonder if you just used ode() instead, if that would resolve your issue? Or alternatively would it make sense to spatially resolve your plasma variable?

Cheers,
Daniel

_______________________________
Daniel C. Reed, PhD.
Postdoctoral Fellow,
Dept. of Earth & Environmental Sciences,
University of Michigan,
Ann Arbor, MI, USA.
email: reeddc at umich.edu
web: www.danielreed.org



On Jun 26, 2013, at 3:53 PM, Dabing Chen <dabing.c at gmail.com> wrote:

> 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]]
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models