From 2fc3f680f8621f9b78d5c10184ee8e9e49206f82 Mon Sep 17 00:00:00 2001 From: jranke Date: Mon, 16 Apr 2012 22:11:48 +0000 Subject: Overhaul of mkinmod in order to make it possible to choose between different model specifications - a variant with minimal use of formation fractions, and a variant with formation fractions everywhere. git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@26 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/mkinmod.R | 160 ++++++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 108 insertions(+), 52 deletions(-) (limited to 'R') diff --git a/R/mkinmod.R b/R/mkinmod.R index 8095932..54b3b5f 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -18,10 +18,12 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see #}}} -mkinmod <- function(...) +mkinmod <- function(..., use_of_ff = "min") { 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'") # The returned model will be a list of character vectors, containing#{{{ # differential equations, parameter names and a mapping from model variables @@ -51,7 +53,7 @@ mkinmod <- function(...) { # Check the type component of the compartment specification#{{{ if(is.null(spec[[varname]]$type)) stop( - "Every argument to mkinmod must be a list containing a type component") + "Every part of the model specification 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(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) { @@ -77,46 +79,66 @@ mkinmod <- function(...) # Create content of differential equations and build parameter list#{{{ for (varname in obs_vars) { - # Add first-order term to first (or only) box for SFO and SFORB#{{{ + # 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 - if(spec[[varname]]$type %in% c("SFO", "SFORB")) { - k_compound <- paste("k", box_1, sep="_") - origin_term <- paste(k_compound, "*", box_1) - parms <- c(parms, k_compound) + if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB decline term + if (use_of_ff == "min") { # Minimum use of formation fractions + # Turn on sink if this is not explicitly excluded by the user by specifying sink=FALSE + if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE + + if(spec[[varname]]$sink) { + # If sink is required, add first-order sink term + k_compound_sink <- paste("k", box_1, "sink", sep="_") + parms <- c(parms, k_compound_sink) + decline_term <- paste(k_compound_sink, "*", box_1) + } else { # otherwise no decline term needed here + decline_term = "" + } + } else { + k_compound <- paste("k", box_1, sep="_") + parms <- c(parms, k_compound) + decline_term <- paste(k_compound, "*", box_1) + } }#}}} - # Construct and add FOMC term and add FOMC parameters if needed#{{{ - if(spec[[varname]]$type == "FOMC") { + if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term # From p. 53 of the FOCUS kinetics report - origin_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1) + decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1) parms <- c(parms, "alpha", "beta") } #}}} - # Construct and add DFOP term and add DFOP parameters if needed#{{{ - if(spec[[varname]]$type == "DFOP") { + if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term # From p. 57 of the FOCUS kinetics report - origin_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1) + decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1) parms <- c(parms, "k1", "k2", "g") } #}}} - # Construct and add HS term and add HS parameters if needed#{{{ - if(spec[[varname]]$type == "HS") { + if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term # From p. 55 of the FOCUS kinetics report - origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1) + decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1) parms <- c(parms, "k1", "k2", "tb") } #}}} # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{ - diffs[[box_1]] <- paste(diffs[[box_1]], "-", origin_term)#}}} - # Add reversible binding for SFORB models#{{{ - if(spec[[varname]]$type == "SFORB") { + diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}} + if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms box_2 = map[[varname]][[2]] - k_free_bound <- paste("k", varname, "free", "bound", sep="_") - k_bound_free <- paste("k", varname, "bound", "free", sep="_") - reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", - k_bound_free, "*", box_2) + if (use_of_ff == "min") { # Minimum use of formation fractions + k_free_bound <- paste("k", varname, "free", "bound", sep="_") + k_bound_free <- paste("k", varname, "bound", "free", sep="_") + parms <- c(parms, k_free_bound, k_bound_free) + reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", + k_bound_free, "*", box_2) + reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", + k_bound_free, "*", box_2) + } else { # Use formation fractions also for the free compartment + f_free_bound <- paste("f", varname, "free", "bound", sep="_") + k_bound_free <- paste("k", varname, "bound", "free", sep="_") + parms <- c(parms, f_free_bound, k_bound_free) + reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2) + reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-", + k_bound_free, "*", box_2) + } diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1) - reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", - k_bound_free, "*", box_2) diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2) - parms <- c(parms, k_free_bound, k_bound_free) } #}}} + # Transfer between compartments#{{{ to <- spec[[varname]]$to if(!is.null(to)) { @@ -128,45 +150,79 @@ mkinmod <- function(...) target_box <- switch(spec[[target]]$type, SFO = target, SFORB = paste(target, "free", sep="_")) - fraction_to_target = paste("f", origin_box, "to", target, sep="_") - diffs[[target_box]] <- paste(diffs[[target_box]], "+", - fraction_to_target, "*", origin_term) - parms <- c(parms, fraction_to_target) + if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", + "SFORB")) { + k_from_to <- paste("k", origin_box, target_box, sep="_") + parms <- c(parms, k_from_to) + diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", + k_from_to, "*", origin_box) + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + k_from_to, "*", origin_box) + } else { + fraction_to_target = paste("f", origin_box, "to", target, sep="_") + parms <- c(parms, fraction_to_target) + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + fraction_to_target, "*", decline_term) + } } - }#}}} + } #}}} }#}}} - model <- list(diffs = diffs, parms = parms, map = map) + 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)) - 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 (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="_") + 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") } - if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB - m[from,to] = paste("- k", from, "free", sep="_") + } + } # }}} + } 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="_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + m[to, from] = ifelse(f.candidate %in% model$parms, + paste(f.candidate, " * k_", from, sep=""), + ifelse(k.candidate %in% model$parms, k.candidate, "0")) } - 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="_") - k.candidate = sub("free.*bound", "free_bound", k.candidate) - k.candidate = sub("bound.*free", "bound_free", k.candidate) - 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 }#}}} -- cgit v1.2.1