[R-sig-dyn-mod] Please help one Law of Mass Action Model

[Jinsong Zhao]

Hi there,

I am trying to use deSolve package to reproduce the Law of Mass Action 
Model that demonstrated in Hannon & Ruth (2014), Chapter 6. The original 
model is run under STELLA. The model equation is listed as:

H(t) = H(t - dt) + (-ΔH) * dt
INIT H = 200
OUTFLOWS:
ΔH = ΔH2O*H_PER_H2O
H2O(t) = H2O(t - dt) + (ΔH2O) * dt
INIT H2O = 0
INFLOWS:
ΔH2O = K*H*O
O(t) = O(t - dt) + (-ΔO) * dt
INIT O = 100
OUTFLOWS:
ΔO = O_PER_H2O*ΔH2O
H_PER_H2O = 2
K = 0.005
O_PER_H2O = 1

The computer models of the book are also available at
www.iseesystems.com/modelingdynamicbiologicalsystems.

And the code I written in R is:

mass <- function(time, state, params) {
    with(as.list(c(state, params)), {
            dH2O <- K3 * H * O
            dH <- -K1 * H2O
            dO <- -K2 * H2O
            list(c(dH, dO, dH2O))
    })
}

params <- c(K1 = 2,
             K2 = 1,
             K3 = 0.005)
state <- c(H = 200,
            O = 100,
            H2O = 0)
time <- seq(0,5)
out <- ode(state, time, mass, params)
plot(out)

However, I can not get the same or similar results.

I would really appreciate if you would like to give me any suggestion or 
hint.

Best regards,
Jinsong