Version 3.2v3.2

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)