diff options
Diffstat (limited to 'R/mkinmod.R')
| -rw-r--r-- | R/mkinmod.R | 168 |
1 files changed, 61 insertions, 107 deletions
diff --git a/R/mkinmod.R b/R/mkinmod.R index 31c778cb..b49c813b 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,7 +1,7 @@ # $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $
-# Copyright (C) 2010 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.de
+# Copyright (C) 2010-2012 Johannes Ranke
+# Contact: jranke@uni-bremen.de
# This file is part of the R package mkin
@@ -33,21 +33,27 @@ mkinmod <- function(...) 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",
+ 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
# Establish list of differential equations
+ # as well as map from observed compartments to 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()
+ if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {
+ stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
+ "which is assumed to be the source compartment"))
+ }
# New (sub)compartments (boxes) needed for the model type
new_boxes <- switch(spec[[varname]]$type,
@@ -62,124 +68,75 @@ mkinmod <- function(...) # Start a new differential equation for each new box
new_diffs <- paste("d_", new_boxes, " =", sep="")
+ names(new_diffs) <- new_boxes
+ diffs <- c(diffs, new_diffs)
+ }
+
+ # 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
+ 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)
+ }
- # Turn on sink if not specified otherwise
- if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
# 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]])
- new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
- new_parms <- c("alpha", "beta")
- ff <- vector()
+ origin_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(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]])
- new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
- new_parms <- c("k1", "k2", "g")
- ff <- vector()
+ 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)
+ parms <- c(parms, "k1", "k2", "g")
}
# 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]])
- new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
- new_parms <- c("k1", "k2", "tb")
- ff <- vector()
+ origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)
+ parms <- c(parms, "k1", "k2", "tb")
}
- # Construct terms for transfer to sink and add if appropriate
-
- 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")) {
- 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 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 if appropriate
+ # Add reversible binding for SFORB models
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)
+ 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)
+ 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)
}
- # 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) {
+ # Transfer between compartments
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
+ # Name of box from which transfer takes place
+ origin_box <- box_1
+
+ # Add transfer terms to listed compartments
for (target in to) {
target_box <- switch(spec[[target]]$type,
SFO = target,
SFORB = paste(target, "free", sep="_"))
- if(spec[[varname]]$type %in% c("SFO", "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)
- }
- if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) {
- fraction_to_target = paste("f_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 {
- fraction_really_to_target = paste(fraction_left, " * ",
- fraction_to_target, sep="")
- fraction_left = paste(fraction_left, " * ",
- fraction_not_to_target, sep="")
- }
- ff[target_box] = fraction_really_to_target
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- ff[target_box], "*", nonlinear_term)
- parms <- c(parms, fraction_to_target)
- }
+ 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)
}
}
}
@@ -193,26 +150,23 @@ mkinmod <- function(...) for (from in boxes) {
for (to in boxes) {
if (from == to) {
- 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")
+ k.candidate = paste("k", from, sep="_")
+ m[from,to] = ifelse(k.candidate %in% model$parms,
+ paste("-", k.candidate), "0")
} else {
+ 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)
- k.effective = intersect(model$parms, k.candidate)
- m[to, from] = ifelse(length(k.effective) > 0,
- k.effective, "0")
+ 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)
+ return(model)
}
|