diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinfit.R | 66 | ||||
| -rw-r--r-- | R/mkinmod.R | 76 |
2 files changed, 142 insertions, 0 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R new file mode 100644 index 00000000..9651fd66 --- /dev/null +++ b/R/mkinfit.R @@ -0,0 +1,66 @@ +mkinfit <- function(mkinmod, observed,
+ parms.ini = rep(0.1, length(mkinmod$parms)),
+ state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
+ fixed_parms = rep(FALSE, length(mkinmod$parms)),
+ fixed_initials = c(FALSE, rep(TRUE, length(mkinmod$diffs) - 1)),
+ plot = NULL,
+ err = NULL, weight = "none", scaleVar = FALSE,
+ ...)
+{
+ # Name the parameters if they are not named yet
+ if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms
+ # Create a function calculating the differentials specified by the model
+ mkindiff <- function(t, state, parms) {
+ diffs <- vector()
+ for (box in names(mkinmod$diffs))
+ {
+ diffname <- paste("d", box, sep="_")
+ diffs[diffname] <- with(as.list(c(state, parms)),
+ eval(parse(text=mkinmod$diffs[[box]])))
+ }
+ return(list(c(diffs)))
+ }
+
+ # Name the inital parameter values if they are not named yet
+ if(is.null(names(state.ini))) names(state.ini) <- names(mkinmod$diffs)
+
+ # TODO: Collect parameters to be optimised
+ parms.optim <- parms.ini[!fixed_parms]
+ parms.fixed <- parms.ini[fixed_parms]
+
+ state.ini.optim <- state.ini[!fixed_initials]
+ state.ini.optim.boxnames <- names(state.ini.optim)
+ names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
+ state.ini.fixed <- state.ini[fixed_initials]
+
+ # Define the model cost function
+ cost <- function(P)
+ {
+ if(length(state.ini.optim) > 0) {
+ odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+ names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
+ } else odeini <- state.ini.fixed
+
+ odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
+ # Solve the ODE
+ out <- ode(
+ y = odeini,
+ times = unique(observed$time),
+ func = mkindiff,
+ parms = odeparms)
+
+ # Output transformation for models with ghost compartments like SFORB
+ out_transformed <- data.frame(time = out[,"time"])
+ for (var in names(mkinmod$map)) {
+ if(length(mkinmod$map[[var]]) == 1) {
+ out_transformed[var] <- out[, var]
+ } else {
+ out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
+ }
+ }
+
+ return(modCost(out_transformed, observed, y = "value",
+ err = err, weight = weight, scaleVar = scaleVar))
+ }
+ modFit(cost, c(state.ini.optim, parms.optim), ...)
+}
diff --git a/R/mkinmod.R b/R/mkinmod.R new file mode 100644 index 00000000..7e64f879 --- /dev/null +++ b/R/mkinmod.R @@ -0,0 +1,76 @@ +mkinmod <- function(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))
+{
+ if(!is.list(spec)) stop("spec must be a list of model specifications for each observed variable")
+
+ # The returned model will be a list of two character vectors, containing
+ # differential equations and parameter names
+ parms <- vector()
+ diffs <- vector()
+ map <- list()
+
+ # Establish list of differential equations
+ for (varname in names(spec))
+ {
+ if(!spec[[varname]]$type %in% c("SFO", "SFORB")) stop("Available types are SFO and SFORB only")
+ new_parms <- vector()
+
+ # New (sub)compartments (boxes) needed for the model type
+ new_boxes <- switch(spec[[varname]]$type,
+ SFO = varname,
+ SFORB = paste(varname, c("free", "bound"), sep="_")
+ )
+ map[[varname]] <- new_boxes
+
+ # Start a new differential equation for each new box
+ new_diffs <- paste("d_", new_boxes, " =", sep="")
+
+ # Construct terms for transfer to sink and add if appropriate
+ if(spec[[varname]]$sink) {
+ k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")
+ sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])
+ # Only add sink term to first (or only) box for SFO and SFORB models
+ if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
+ new_diffs[[1]] <- paste(new_diffs[[1]], sink_term)
+ }
+ new_parms <- k_compound_sink
+ }
+
+ # Add reversible binding if appropriate
+ if(spec[[varname]]$type == "SFORB") {
+ k_free_bound <- paste("k", varname, "free", "bound", sep="_")
+ k_bound_free <- paste("k", varname, "bound", "free", sep="_")
+ reversible_binding_terms <- c(
+ paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]),
+ paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]]))
+ new_diffs <- paste(new_diffs, reversible_binding_terms)
+ new_parms <- c(new_parms, k_free_bound, k_bound_free)
+ }
+
+ # Add observed variable to model
+ parms <- c(parms, new_parms)
+ names(new_diffs) <- new_boxes
+ diffs <- c(diffs, new_diffs)
+ }
+
+ # Transfer between compartments
+ for (varname in names(spec)) {
+ to <- spec[[varname]]$to
+ if(!is.na(to)) {
+ origin_box <- switch(spec[[varname]]$type,
+ SFO = varname,
+ SFORB = paste(varname, "free", sep="_"))
+ for (target in to) {
+ target_box <- switch(spec[[target]]$type,
+ SFO = target,
+ SFORB = paste(target, "free", sep="_"))
+ }
+ k_from_to <- paste("k", origin_box, target_box, sep="_")
+ diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", k_from_to, "*", origin_box)
+ diffs[[target_box]] <- paste(diffs[[target_box]], "+", k_from_to, "*", origin_box)
+ parms <- c(parms, k_from_to)
+ }
+ }
+ model <- list(diffs = diffs, parms = parms, map = map)
+ class(model) <- "mkinmod"
+ return(model)
+}
|