Description Usage Details Value References See Also Examples
This is the batch version of the model described by Jähnichen et. al (submitted).
1 |
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:
|
parms |
Vector with the named parameters of the model:
|
inputs |
matrix with time-dependend signals. The two parameters
|
times |
Simulation time and integration interval. |
init |
Vector with start values for |
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 | 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.