[R-sig-dyn-mod] PBTK model - unexpected behavior

[Haselman, Jon]

Hello,

(resending same inquiry with code)

Below is a very simple model that simulates growth of a larval amphibian over a 10 day period with calculation of the volume of distribution of toxicant. The apparent volume of distribution constant is 1 for simplicity sake and to test whether the model will simulate the true Vd equal to the body weight. However, the resulting data suggests the AVD constant is not staying constant because the Vd increases by a factor approximately equal to the time step in relation to the body weight. I apologize if I am overlooking something simple, but could someone please explain why this is happening and a possible solution. The body weight data turns out exactly as we would expect, but if the rest of our larger, more complex PBTK model is to work correctly, we need to see the output data for terms such as Vd ('ordinary' variables?) being more in line with what we would expect.

I understand that Vd in this case is being treated as a state variable that will respond to a 'rate' - the same as growth rate is driving the body weight term. If Vd is not defined as a state variable, but is changing over time as a function of body weight and driving other equations in the model, how do we extract the behavior/data of that term in the model for verification purposes? It appears that only state variables can be 'tracked' and output by the ode function, but in this case when we track Vd as a state variable, it is inappropriately changing because it perhaps should not be a state variable.

Your help is MUCH appreciated.

Thank you,

Jon Haselman
Toxicologist


library(deSolve)

dynmodel <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {

    #-------Calculate Body Weight and volume of distribution

    dBW <- RBW                                               # BW=0.354 g to start
    VdTOX <- AVDTOXC*BW                           # ml


    list(c(dBW, VdTOX))
  })
}
parms <- c(RBW=0.00237, AVDTOXC=1)
state <- c(BW=0.354, VdTOX=0.354)
times <- seq(0, 240, by = 1)

# Model output
out <- ode(y = state, times = times, func = dynmodel, parms = parms)
plot(out)




	[[alternative HTML version deleted]]