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

[Jinsong Zhao]

Hi there,

I got the solution. It's my fault.

The function revised to the following form, then I got the same result.

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

Best,
Jinsong

On 2020/3/5 22:16, Jinsong Zhao wrote:
> 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