[R-sig-dyn-mod] Turbidostat

Uta,

thanks for the nice example. I would do it like follows: declare the 
correction factors (i.e. time dependent parameters) as parameters for 
the modCost function, so that they can be fitted like ordinary model 
parameters. Note also that in the example below we now pass the "signal" 
as an optional parameter explicitly to the turbidostat function.

A few additional notes:

1) The parameters may be not simultaneously identifiable because of 
strong correlation. Ypou may consider to use other start values, e.g.

parms <- c(r=1, korr_s=1)

and/or to adapt the interpolation method (e.g. "constant").

2) Why does the simulation start at time step 1 and not at zero like the 
data?

3) Is there any evidence why this inhibition occurs at the specified 
times only? Is there any chance to describe this a little bit more 
process oriented?

Thomas P.



###########################################################
library(deSolve)
library(FME)

maxTime <- 15
korr_s <- 0.3  #not fitted but badly assumed - just as an example
signal <- approxfun(1:maxTime, c(rep(korr_s,3),rep(1,10),rep(korr_s,2)))

turbidostat <- function(time, E0, parms, signal){
     korr <- signal(time)
     with(as.list(c(E0, parms)), {
       dE <- r*korr*E
     list(c(dE))
   })
}

E0 <- c(E=0.1)
parms <- c(r=0.5)
times <- seq(1,15,0.1)

eventdat <- data.frame(var = c("E","E","E","E","E"), time = c(7,9,11,13,15),
   value=c(rep(0.5,5)),method=c("rep","rep","rep","rep","rep"))


#####observed data
Obs <- (data.frame(time = c(0,1,2,3,7,9,11,13,15),
                    E = c(0.1,0.11,0.2,0.3,1.8,1.4,1.45,1.46,1.0)))

out <- ode(E0,times,turbidostat,parms, events=list(data=eventdat), 
signal=signal)

plot(Obs$time,Obs$E,pch=16,cex=2.0)
lines(out, type = "l", lwd = 2)

Cost <- function(parms){
   pks <- parms["korr_s"]
   signal <- approxfun(1:maxTime, c(rep(pks, 3),rep(1, 10),rep(pks, 2)))
   model <- ode(E0, times, turbidostat, parms,
     events=list(data=eventdat), signal=signal)
   modCost(model = model, obs = Obs, x = "time")
}

parms <- c(r=0.5, korr_s=0.3)

Cost(parms)

Fit <- modFit(f=Cost,p=parms)
deviance(Fit)
Fit$par

out <- ode(E0,times,turbidostat,Fit$par, events=list(data=eventdat), 
signal=signal)

plot(Obs$time,Obs$E,pch=16,cex=2.0)
lines(out, type = "l", lwd = 2)

## But note that there is *very strong* correlation
## between these two parameters:
summary(Fit)

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20101208/c3ed8b52/attachment.pl>

Thomas,

thanks again.

2) Why does the simulation start at time step 1 and not at zero like
the data?

It was my mistake.

3) Is there any evidence why this inhibition occurs at the specified
times only? Is there any chance to describe this a little bit more
process oriented?

Toxic substances were given during the initial phase (first four time
steps) and before the last time step. During the "middle time steps",
algae could grow without inhibition.
Unfortunately, I do not have a chance for a mechanistic explanation
since we have only the data for the algae. I did not present all of
them. Three experiments were done made with three different levels of
toxic substances. (I know the levels by number, e.g., Uran = 0; 0.2 and
0.5mu mol, but do not know the absorption by the algae).

The three data sets thus differ mainly in their initial phase and in the
last phase. The "mid phase" is determined by the last concentration of
the initial phase (time step four in the example).

I see. So you may consider to handle korr_s as an ordinary parameter of 
the model that can be fitted with the usual procedure. To model the 
presence of the inhibitor one may use either (1) a signal set to 0 resp. 
1 or (2) a state variable for the toxic substance that is added at 
specified times (by help of a signal or an event) and that is removed by 
dilution of the culture -- and optionally something like a 1st order 
decay or an uptake by the cells.

Thomas P.