[R-sig-dyn-mod] Simple model, odd results...

[Jinsong Zhao]

Hi there,

I am newbie to dynamic system modeling. So I start from simple model. I 
have a simple mode with diagram like the following:

  ELc
   ^
   |           EPc
  Enzc  <-------------+
   |                  |
   v                  |
  Dc           Uc     |           Re
---->  DOC  ---->  Biomc  -----> Rm
                                  Rg

The state variable is Enzc, DOC and Biomc. And I have the following model,
### Model equations
C_Only <- function(time, state, parameters) {
     with(as.list(c(state,parameters)), {
          # rate of change
          Dc <- Kd * Enzc
          dEnzc <- EPc - ELc
          dBiomc <- Uc - EPc - Re - Rm - Rg
          dDOC <- Dc - Uc
          EPc <- Ke * Biomc
          ELc <- Kl * Enzc
          Re <- EPc * (1 - SUE) / SUE
          Rm <- Km * Biomc
          Rg <- (Uc - EPc - Re - Rm) * (1 - SUE)
          # return the rate of change
          list(c(Dc, dDOC, Uc, dEnzc, EPc, ELc, dBiomc, Re, Rm, Rg))
     })
}

The parameters are following,
### Model parameters
parameters <- c(Kd = 1,
                 DOC = 0.0,
                 Ke = 0.05,
                 Kl = 0.05,
                 Km = 0.022,
                 SUE = 0.5)

The initial value for state variables are following,
### State variables
state <- c(Dc = 0,
            DOC = 0,
            Uc = 0,
            Enzc = 0.7,
            EPc = 0,
            ELc = 0,
            Biomc = 29,
            Re = 0,
            Rm = 0,
            Rg = 0)

### Model integration
library(deSolve)
times <- seq(0, 2000, by = 1)
out <- ode(y = state, times = times, func = C_Only, parms = parameters)

I plot the Biomc using,

plot(out[,c(1,8)], type = "l")

you may see Biomc changed periodically. (In fact, it is not correct, 
since, all variable list above should be positive).

If I set Uc in state to 0, this pattern seems always there. However, 
when set Uc to any positive value, all state variables are out of control.

Thus, my question are:

(1) based on the diagram, is the model correct?
(2) how to set correct initial values?

I have try another scheme. I eliminate all algebra equation in the 
model, and get:

C_Only <- function(time, state, parameters) {
     with(as.list(c(state,parameters)), {
          # rate of change
          dEnzc <- Ke * Biomc - Kl * Enzc
          dBiomc <- (Uc - Km * Biomc) * SUE - Ke * Biomc
          dDOC <- Kd * Enzc - Uc
          # return the rate of change
          list(c(dDOC, Uc, dEnzc, dBiomc))
     })
}

I think this equals to the above one. I give the same parameters, and 
the following initial values (and I think they are same):

### State variables
state <- c(DOC = 0,
            Uc = 0,
            Enzc = 0.7,
            Biomc = 29)

However, the results are different. I don't understand what's the reason.

Any suggestions or comments will be really appreciated. Thanks in advance.

Best wishes,
Jinsong