# Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke
# $Id$
# 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(..., 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'")
# 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="")
# 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]])
}
}
# 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, 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)
}
