diff options
Diffstat (limited to 'CakeModel.R')
| -rw-r--r-- | CakeModel.R | 261 |
1 files changed, 261 insertions, 0 deletions
diff --git a/CakeModel.R b/CakeModel.R new file mode 100644 index 0000000..a93d7af --- /dev/null +++ b/CakeModel.R @@ -0,0 +1,261 @@ +# Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke +# $Id: CakeModel.R 216 2011-07-05 14:35:03Z nelr $ + +# The CAKE R modules are based on mkin +# Portions Johannes Ranke, 2010 +# Contact: mkin-devel@lists.berlios.de + +# This version has been modified to parameterise SFO as k and flow fractions +# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta +# Tessella Project Reference: 6245 + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>.” + +CakeModel <- function(...) +{ + spec <- list(...) + obs_vars <- names(spec) + + # The returned model will be a list of character vectors, containing + # differential equations, parameter names and a mapping from model variables + # to observed variables. If possible, a matrix representation of the + # differential equations is included + parms <- vector() + diffs <- vector() + map <- list() + if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) { + mat = FALSE + if(!is.null(spec[[1]]$to)) { + message <- paste("Only constant formation fractions over time are implemented.", + "Depending on the reason for the time dependence of degradation, this may be unrealistic", + sep="\n") + warning(message) + } else message <- "ok" + } else mat = TRUE + + mat = FALSE # XYZZY Not implemented yet assuming it should be + nlt = list() + + # Establish list of differential equations + for (varname in obs_vars) + { + if(is.null(spec[[varname]]$type)) stop( + "Every argument to mkinmod must be a list containing a type component") + if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop( + "Available types are SFO, FOMC, DFOP, HS and SFORB only") + new_parms <- vector() + + # New (sub)compartments (boxes) needed for the model type + new_boxes <- switch(spec[[varname]]$type, + SFO = varname, + FOMC = varname, + DFOP = varname, + HS = varname, + SFORB = paste(varname, c("free", "bound"), sep="_") + ) + map[[varname]] <- new_boxes + names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes)) + + # Start a new differential equation for each new box + new_diffs <- paste("d_", new_boxes, " =", sep="") + + # Turn on sink if not specified otherwise + if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE + + #@@@ADD SFO k HERE !!!!!!!!!!!!! + # Construct and add SFO term and add SFO parameters if needed + if(spec[[varname]]$type == "SFO") { + # From p. 53 of the FOCUS kinetics report + k_term <- paste("k", new_boxes[[1]], sep="_") + nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ") + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- k_term + ff <- vector() + } + + # Construct and add FOMC term and add FOMC parameters if needed + if(spec[[varname]]$type == "FOMC") { + if(match(varname, obs_vars) != 1) { + stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment") + } +# if(spec[[varname]]$sink == FALSE) { +# stop("Turning off the sink for the FOMC model is not implemented") +# } + # From p. 53 of the FOCUS kinetics report + nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]]) + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c("alpha", "beta") + ff <- vector() + } + + # Construct and add DFOP term and add DFOP parameters if needed + if(spec[[varname]]$type == "DFOP") { + if(match(varname, obs_vars) != 1) { + stop("Type DFOP is only possible for the first compartment, which is assumed to be the source compartment") + } +# if(spec[[varname]]$sink == FALSE) { +# stop("Turning off the sink for the DFOP model is not implemented") +# } + # From p. 57 of the FOCUS kinetics report + nonlinear_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]]) + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c("k1", "k2", "g") + ff <- vector() + } + + # Construct and add HS term and add HS parameters if needed + if(spec[[varname]]$type == "HS") { + if(match(varname, obs_vars) != 1) { + stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment") + } +# if(spec[[varname]]$sink == FALSE) { +# stop("Turning off the sink for the HS model is not implemented") +# } + # From p. 55 of the FOCUS kinetics report +# nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]]) + nonlinear_term <- paste("((k1 - k2) / ( 1 + exp( 10*(time - tb) ) ) + k2)", "*", new_boxes[[1]]) + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c("k1", "k2", "tb") + ff <- vector() + } + + # Construct terms for transfer to sink and add if appropriate + + #@@@@ REMOVE THIS ????? + if(spec[[varname]]$sink) { + # Add first-order sink term to first (or only) box for SFO and SFORB +# if(spec[[varname]]$type %in% c("SFO", "SFORB")) { + if(spec[[varname]]$type == "SFORB") { + k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_") + sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]]) + 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 obs_vars) { + to <- spec[[varname]]$to + if(!is.null(to)) { + origin_box <- switch(spec[[varname]]$type, + SFO = varname, + FOMC = varname, + DFOP = varname, + HS = varname, + SFORB = paste(varname, "free", sep="_")) + fraction_left <- NULL + for (target in to) { + target_box <- switch(spec[[target]]$type, + SFO = target, + SFORB = paste(target, "free", sep="_")) + # SFO is no longer special + #if(spec[[varname]]$type %in% c("SFO", "SFORB")) { + if(spec[[varname]]$type == "SFORB") { + 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) + } + # Handle SFO like the others +# if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) { + if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "SFO")) { + if ( length(to)==1 && !spec[[varname]]$sink ) { + # There is only one output, so no need for any flow fractions. Just add the whole flow from the parent + diffs[[target_box]] <- paste(diffs[[target_box]], "+", spec[[varname]]$nlt) + } else { + fraction_to_target = paste("f", varname, "to", target, sep="_") + fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", sep="") + if(is.null(fraction_left)) { + fraction_really_to_target = fraction_to_target + fraction_left = fraction_not_to_target + } else { + # If this is the last output and there is no sink, it gets what's left + if ( target==tail(to,1) && !spec[[varname]]$sink ) { + fraction_really_to_target = fraction_left + } else { +# (1-fa)fb version +# fraction_really_to_target = paste(fraction_left, " * ", fraction_to_target, sep="") +# fraction_left = paste(fraction_left, " * ", fraction_not_to_target, sep="") +# fb version + fraction_really_to_target = fraction_to_target + fraction_left = paste(fraction_left, " - ", fraction_to_target, sep="") + } + } + ff[target_box] = fraction_really_to_target + diffs[[target_box]] <- paste(diffs[[target_box]], "+", ff[target_box], "*", spec[[varname]]$nlt) + # Add the flow fraction parameter (if it exists) + if ( target!=tail(to,1) || spec[[varname]]$sink ) { + parms <- c(parms, fraction_to_target) + } + } + } + } + } + } + model <- list(diffs = diffs, parms = parms, map = map) + + # Create coefficient matrix if appropriate + if (mat) { + boxes <- names(diffs) + n <- length(boxes) + m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes)) + for (from in boxes) { + for (to in boxes) { + if (from == to) { + #@@@@ OMIT NEXT LINE? !!!!! + k.candidate = paste("k", from, c(boxes, "sink"), sep="_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + k.effective = intersect(model$parms, k.candidate) + m[from,to] = ifelse(length(k.effective) > 0, + paste("-", k.effective, collapse = " "), "0") + } else { + k.candidate = paste("k", from, to, sep="_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + k.effective = intersect(model$parms, k.candidate) + m[to, from] = ifelse(length(k.effective) > 0, + k.effective, "0") + } + } + } + model$coefmat <- m + } + + if (exists("ff")) model$ff = ff + class(model) <- "mkinmod" + invisible(model) +} |