mcyst_batch_mc: Microcystin Production Model (Batch Version for Monte Carlo...
In simecolModels: Model collection for the simecol package

Description Usage Details Value References See Also Examples

This is the batch version of the model described by Jähnichen et. al (submitted).

1	mcyst_batch_mc()

To see all details, please have a look into the implementation and the original publications.

S4 object according to the odeModel specification. The object contains the following slots:

`main`	The differential equations for cell abundance and Microcystin concentration: cells abundance of cells (1/L), mcyst microcystin concentration (fg/L).
`parms`	Vector with the named parameters of the model: mu Growth rate of cells (1/d), K Carrying Capacity (cells /L), dMmin minimal value of MCYST depletion rate (1/d), dMmax maximal value of MCYST depletion rate (1/d), dMrho autocorrelation coefficient for MCYST decay rate, pmin minimal value of MCYST production coefficient (fg/cell), pmax minimal value of MCYST production coefficient (fg/cell), prho autocorrelation coefficient for MCYST production.
`inputs`	matrix with time-dependend signals. The two parameters `dM` and `p` are considered as stochastic, time dependend signals in the Monte Carlo version of the MCYST production model. x time (d), p MCYST production coefficient (fg/cell), dM MCYST depletion rate (1/d)
`times`	Simulation time and integration interval.
`init`	Vector with start values for `cells` and `mcyst`.
`solver`	Character string with the integration method.

Jähnichen, S., Ihle, T. and Petzoldt, T. (2008). Variability of microcystin cell quota: A small model explains dynamics and equilibria. Limnologica - Ecology and Management of Inland Waters, 38, 339-349. URL http://dx.doi.org/10.1016/j.limno.2008.05.003

sim, parms, init, times.

  mcystmod <- mcyst_batch_mc()
  ## sim initializes automatically and draws new random numbers for the signals
  plot(sim(mcystmod))

  ## derive a scenario

  parms(mcystmod)["prho"] <- 0.7
  parms(mcystmod)["dMmin"] <- 0.02
  parms(mcystmod)["dMmax"] <- 0.02
  
  ## Monte Carlo Simulation
  for (i in 1:10) { # number of runs, increase this to > 1000 for real simulations
     out <- out(sim(mcystmod))
    if (i==1) allout <- out else allout <- rbind(allout, out)
  }
  
  mu_cells <- diff(log(allout$cells))/diff(allout$time) # spezific growth rate
  mu_mcyst <- diff(log(allout$mcyst))/diff(allout$time) # spezific MCYST-production
  q_mcyst <- allout$mcyst/allout$cells
  
  mu <- parms(mcystmod)["mu"]
  K  <- parms(mcystmod)["K"]
  
  deriv <- mu * allout$cells*(1-allout$cells/K)
  par(mfrow=c(3,3))
  
  plot(allout$time, allout$cells, col="green", type="p",
    xlab="time", ylab="Cells", ylim=c(0, max(out$cells)))
  plot(allout$time, allout$mcyst*1e-6, col="red", type="p",
    xlab="time", ylab="MCYST", ylim=c(0, max(out$mcyst*1e-6)))
  plot(allout$time, q_mcyst, col="red",  type="p",
    xlab="time", ylab="Q_MCYST")
  plot(mu_cells, mu_mcyst,col="red",type="p",
    xlab="mu_Zellen", ylab="mu_MCYST",ylim=c(-0,0.5))
  
  plot(allout$time[-1],mu_cells,col="blue",type="p",
    xlab="time", ylab="mu_cells")
  plot(allout$time[-1],mu_mcyst,col="blue",type="p",
    xlab="time", ylab="mu_MCYST", ylim=c(0,0.5))
  plot(allout$time[-1],mu_mcyst/mu_cells,col="blue",type="p",
    xlab="time", ylab="mu_MCYST/mu_cells", ylim=c(0.8,1.2))
  plot(allout$time,deriv, col="blue",type="p",
    xlab="time", ylab="dx.dt" )
  
  ## derive some numeric results
  results1 <- data.frame(allout$time,
                         allout$cells,
                         allout$mcyst, q_mcyst)
  results2 <- data.frame(mu_cells, mu_mcyst, mu_mcyst/mu_cells)
  
  qmax  <- aggregate(allout$mcyst/allout$cells,list(allout$time),max)
  qmin  <- aggregate(allout$mcyst/allout$cells,list(allout$time),min)
  qmean <- aggregate(allout$mcyst/allout$cells,list(allout$time),mean)
  qsd   <- aggregate(allout$mcyst/allout$cells,list(allout$time),sd)
  
  mumean<-aggregate(diff(log(allout$cells))/
                    diff(allout$time),list(allout$time[-1]),mean)