diff options
Diffstat (limited to 'CakeModel.R')
| -rw-r--r-- | CakeModel.R | 287 |
1 files changed, 101 insertions, 186 deletions
diff --git a/CakeModel.R b/CakeModel.R index ce6bc01..3a32779 100644 --- a/CakeModel.R +++ b/CakeModel.R @@ -1,15 +1,13 @@ # Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke -# $Id$ -# The CAKE R modules are based on mkin +# Some of 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 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414 -# This program is free software: you can redistribute it and/or modify +# The CAKE R modules are 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. @@ -20,98 +18,80 @@ # 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/>.” +# along with this program. If not, see <http://www.gnu.org/licenses/>. CakeModel <- function(..., use_of_ff = "max") { spec <- list(...) obs_vars <- names(spec) - if (!use_of_ff %in% c("min", "max")) - stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'") + if (!use_of_ff %in% c("min", "max")) stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'") + # 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") + if(is.null(spec[[varname]]$type)) stop("Every argument to CakeModel must be a list containing a type component") + if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) stop("Available types are SFO, FOMC, DFOP, HS and IORE 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="_") - ) + if (spec[[varname]]$type == "DFOP" && varname != "Parent"){ + # If this is a DFOP metabolite, we need to form a system of two differential equations + # (see FOCUS page 137) + new_boxes <- paste(varname, c("1", "2"), sep="_") + } else{ + new_boxes <- varname + } + 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="") - # Get the name of the box(es) we are working on for the decline term(s) - box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB - # 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") { if (use_of_ff == "min") { # Minimum use of formation fractions if(spec[[varname]]$sink) { - # From p. 53 of the FOCUS kinetics report - k_term <- paste("k", new_boxes[[1]], "sink", 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() - decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) - } else { # otherwise no decline term needed here - decline_term = "0" - } - } else { - 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() - decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) - } - } + # From p. 53 of the FOCUS kinetics report + k_term <- paste("k", new_boxes[[1]], "sink", 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() + decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) + } else { # otherwise no decline term needed here + decline_term = "0" + } + } else { + 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() + decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) + } + } # 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 @@ -122,18 +102,31 @@ CakeModel <- function(..., use_of_ff = "max") # 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") + k1_term = paste("k1", varname, sep="_") + k2_term = paste("k2", varname, sep="_") + g_term = paste("g", varname, sep="_") + + new_parms <- c(k1_term, k2_term, g_term) ff <- vector() + + if (varname == "Parent"){ + # From p. 57 of the FOCUS kinetics report + # Looks like this: paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]]) + nonlinear_term <- paste("((", k1_term, " * ", g_term, " * exp(-", k1_term, " * time) + ", k2_term, " * (1 - ", g_term, ") * exp(-", k2_term, " * time)) / (", g_term, " * exp(-", k1_term, " * time) + (1 - ", g_term, ") * exp(-", k2_term, " * time))) *", new_boxes[[1]]) + + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + } else { + # From p. 137 of the FOCUS kinetics report - two sub-compartments with g applied to initial conditions (in ODE pre-processing) and in formation terms. + first_nonlinear_term <- paste(k1_term, new_boxes[[1]], sep=" * ") + second_nonlinear_term <- paste(k2_term, new_boxes[[2]], sep=" * ") + + overall_term_for_formation <- paste("(", g_term, " * ", first_nonlinear_term, " + (1 - ", g_term, ") *", second_nonlinear_term, ")") + spec[[varname]]$nlt <- overall_term_for_formation + + new_diffs[[1]] <- paste(new_diffs[[1]], "-", first_nonlinear_term) + new_diffs[[2]] <- paste(new_diffs[[2]], "-", second_nonlinear_term) + } } # Construct and add HS term and add HS parameters if needed @@ -141,82 +134,52 @@ CakeModel <- function(..., use_of_ff = "max") 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]]) + # From p. 55 of the FOCUS kinetics report: nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]]) + # Replaced this with the smoothed version on the next line to ease convergence. + nonlinear_term <- paste("((k1 - k2) / ( 1 + exp( 300*(time - tb)/(tb + 0.00001) ) ) + 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 - } + # Construct and add IORE term and add IORE parameters if needed + if(spec[[varname]]$type == "IORE") { + k_term <- paste("k", new_boxes[[1]], sep="_") + nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ") + nonlinear_term <- paste(nonlinear_term, "^N") + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c(k_term, "N") + ff <- vector() } - - # 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, + target_boxes <- 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")) { + DFOP = paste(target, c("1", "2"), sep="_")) + + if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "SFO", "IORE")) { 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) + formation_term <- 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 @@ -225,86 +188,38 @@ CakeModel <- function(..., use_of_ff = "max") 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) + + ff[target] = fraction_really_to_target + formation_term <- paste(ff[target], "*", 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) } } + + if (spec[[target]]$type == "DFOP") { + g_term <- paste("g", target, sep="_") + + first_formation_term <- paste(g_term, "*", formation_term) + second_formation_term <- paste("(1 -", g_term, ") *", formation_term) + + diffs[[target_boxes[[1]]]] <- paste(diffs[[target_boxes[[1]]]], "+", first_formation_term) + diffs[[target_boxes[[2]]]] <- paste(diffs[[target_boxes[[2]]]], "+", second_formation_term) + } else{ + diffs[[target_boxes[[1]]]] <- paste(diffs[[target_boxes[[1]]]], "+", formation_term) + } } } } } + model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff) - # Create coefficient matrix if appropriate#{{{ - if (mat) { - boxes <- names(diffs) - n <- length(boxes) - m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes)) - - if (use_of_ff == "min") { # Minimum use of formation fractions - for (from in boxes) { - for (to in boxes) { - if (from == to) { # diagonal elements - 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 { # off-diagonal elements - k.candidate = paste("k", from, to, sep="_") - if (sub("_free$", "", from) == sub("_bound$", "", to)) { - k.candidate = paste("k", sub("_free$", "_free_bound", from), sep="_") - } - if (sub("_bound$", "", from) == sub("_free$", "", to)) { - k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep="_") - } - k.effective = intersect(model$parms, k.candidate) - m[to, from] = ifelse(length(k.effective) > 0, - k.effective, "0") - } - } - } - } else { # Use formation fractions where possible - for (from in boxes) { - for (to in boxes) { - if (from == to) { # diagonal elements - k.candidate = paste("k", from, sep="_") - m[from,to] = ifelse(k.candidate %in% model$parms, - paste("-", k.candidate), "0") - if(grepl("_free", from)) { # add transfer to bound compartment for SFORB - m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_")) - } - if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB - m[from,to] = paste("- k", from, "free", sep="_") - } - m[from,to] = m[from,to] - } else { # off-diagonal elements - f.candidate = paste("f", from, "to", to, sep="_") - k.candidate = paste("k", from, to, sep="_") - # SFORB with maximum use of formation fractions not implemented, see above - m[to, from] = ifelse(f.candidate %in% model$parms, - paste(f.candidate, " * k_", from, sep=""), - ifelse(k.candidate %in% model$parms, k.candidate, "0")) - } - } - } - } - model$coefmat <- m - } - if (exists("ff")) model$ff = ff class(model) <- "mkinmod" invisible(model) |