From 1718d434efae26de02754c6622c43f4dc9e624b9 Mon Sep 17 00:00:00 2001 From: jranke Date: Thu, 15 Mar 2012 15:54:14 +0000 Subject: Update kinfit and mkin to the latest version published on BerliOS. git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@17 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/mkinmod.R | 154 +++++++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 132 insertions(+), 22 deletions(-) (limited to 'R/mkinmod.R') diff --git a/R/mkinmod.R b/R/mkinmod.R index 643e855e..31c778cb 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,3 +1,23 @@ +# $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $ + +# Copyright (C) 2010 Johannes Ranke +# Contact: mkin-devel@lists.berlios.de + +# This file is part of the R package mkin + +# mkin 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 + mkinmod <- function(...) { spec <- list(...) @@ -5,22 +25,36 @@ mkinmod <- function(...) # The returned model will be a list of character vectors, containing # differential equations, parameter names and a mapping from model variables - # to observed 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 # 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", "SFORB")) stop("Available types are SFO, FOMC and SFORB only") + 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 @@ -29,25 +63,64 @@ mkinmod <- function(...) # Start a new differential equation for each new box new_diffs <- paste("d_", new_boxes, " =", sep="") - # Construct terms for transfer to sink and add if appropriate + # 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() + } + + # 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() + } + + # 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() + } + + # 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 models + # 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 } - if(spec[[varname]]$type == "FOMC") { - if(match(varname, obs_vars) != 1) { - stop("Type FOMC is only allowed for the first compartment, which is assumed to be the source compartment") - } - # From p. 53 of the FOCUS kinetics report - fomc_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]]) - new_diffs[[1]] <- paste(new_diffs[[1]], "-", fomc_term) - new_parms <- c("alpha", "beta") - } } # Add reversible binding if appropriate @@ -66,7 +139,6 @@ mkinmod <- function(...) names(new_diffs) <- new_boxes diffs <- c(diffs, new_diffs) } - # Transfer between compartments for (varname in obs_vars) { to <- spec[[varname]]$to @@ -74,6 +146,8 @@ mkinmod <- function(...) 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) { @@ -82,27 +156,63 @@ mkinmod <- function(...) 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) + 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 == "FOMC") { + 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="") + 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="") + fraction_really_to_target = paste(fraction_left, " * ", + fraction_to_target, sep="") + fraction_left = paste(fraction_left, " * ", + fraction_not_to_target, sep="") } - diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_really_to_target, "*", fomc_term) + 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) } } } } 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) { + 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) } -- cgit v1.2.1