Add the logistic model

6 files changed, 43 insertions, 17 deletions

diff --git a/R/endpoints.R b/R/endpoints.R
index ac1e3e7c..80450185 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -159,6 +159,14 @@ endpoints <- function(fit) {
       ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1

       ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2

     }

+    if (type == "logistic") {

+      # FOCUS kinetics (2014) p. 67

+      kmax = parms.all["kmax"]

+      k0 = parms.all["k0"]

+      r = parms.all["r"]

+      DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))

+      DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))

+    }

     ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)

   }

   return(ep)

diff --git a/R/logistic.solution.R b/R/logistic.solution.R
new file mode 100644
index 00000000..a3bddab3
--- /dev/null
+++ b/R/logistic.solution.R
@@ -0,0 +1,4 @@
+logistic.solution <- function(t, parent.0, kmax, k0, r)
+{
+	parent = parent.0 * (kmax / (kmax - k0 + k0 * exp (r * t))) ^(kmax/r)
+}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index b27f67b4..40413125 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2018 Johannes Ranke

+# Copyright (C) 2010-2019 Johannes Ranke

 # Portions of this code are copyright (C) 2013 Eurofins Regulatory AG

 # Contact: jranke@uni-bremen.de

 # The summary function is an adapted and extended version of summary.modFit

@@ -48,7 +48,7 @@ mkinfit <- function(mkinmod, observed,
 {

   # Check mkinmod and generate a model for the variable whith the highest value

   # if a suitable string is given

-  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE")

+  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")

   if (class(mkinmod) != "mkinmod") {

     presumed_parent_name = observed[which.max(observed$value), "name"]

     if (mkinmod[[1]] %in% parent_models_available) {

@@ -153,6 +153,9 @@ mkinfit <- function(mkinmod, observed,
     if (parmname == "k2") parms.ini[parmname] = 0.01

     if (parmname == "tb") parms.ini[parmname] = 5

     if (parmname == "g") parms.ini[parmname] = 0.5

+    if (parmname == "kmax") parms.ini[parmname] = 0.1

+    if (parmname == "k0") parms.ini[parmname] = 0.0001

+    if (parmname == "r") parms.ini[parmname] = 0.2

   }

   # Default values for formation fractions in case they are present

   for (box in mod_vars) {

@@ -376,7 +379,7 @@ mkinfit <- function(mkinmod, observed,
     lower[index_k] <- 0

     index_k__iore <- grep("^k__iore_", names(lower))

     lower[index_k__iore] <- 0

-    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))

+    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb", "r"), names(lower))

     lower[other_rate_parms] <- 0

   }

 

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 491b3d0a..2805ef54 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2015 Johannes Ranke {{{

+# Copyright (C) 2010-2015,2019 Johannes Ranke {{{

 # Contact: jranke@uni-bremen.de

 

 # This file is part of the R package mkin

@@ -42,8 +42,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
   parms <- vector()

   # }}}

 

-  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{

-  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS")) {

+  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{

+  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) {

     mat = FALSE

   } else mat = TRUE

   #}}}

@@ -57,10 +57,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
     # Check the type component of the compartment specification {{{

     if(is.null(spec[[varname]]$type)) stop(

       "Every part of the model specification must be a list containing a type component")

-    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB")) stop(

-      "Available types are SFO, FOMC, IORE, DFOP, HS and SFORB only")

-    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,",

+    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop(

+      "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only")

+    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "logistic") & match(varname, obs_vars) != 1) {

+        stop(paste("Types FOMC, DFOP, HS and logistic are only implemented for the first compartment,",

                    "which is assumed to be the source compartment"))

     }

     #}}}

@@ -71,6 +71,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
       IORE = varname,

       DFOP = varname,

       HS = varname,

+      logistic = varname,

       SFORB = paste(varname, c("free", "bound"), sep = "_")

     )

     map[[varname]] <- new_boxes

@@ -141,6 +142,11 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
       decline_term <- paste(HS_decline, "*", box_1)

       parms <- c(parms, "k1", "k2", "tb")

     } #}}}

+    if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term

+      # From p. 67 of the FOCUS kinetics report (2014)

+      decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1)

+      parms <- c(parms, "kmax", "k0", "r")

+    } #}}}

     # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{

     diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}

     if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms

diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 8e0823a8..c36d724a 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2016,2018 Johannes Ranke
+# Copyright (C) 2010-2016,2018,2019 Johannes Ranke
 # Some lines in this code are copyright (C) 2013 Eurofins Regulatory AG
 # Contact: jranke@uni-bremen.de
 
@@ -83,7 +83,11 @@ mkinpredict.mkinmod <- function(x,
           evalparse(parent.name),
           evalparse(paste("k", parent.name, "bound", sep="_")),
           evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
-          evalparse(paste("k", parent.name, "sink", sep="_")))
+          evalparse(paste("k", parent.name, "sink", sep="_"))),
+      logistic = logistic.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("kmax"), evalparse("k0"),
+          evalparse("r"))
     )
     out <- data.frame(outtimes, o)
     names(out) <- c("time", sub("_free", "", parent.name))
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index c871c52a..f69f4ebd 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2014 Johannes Ranke

+# Copyright (C) 2010-2014,2019 Johannes Ranke

 # Contact: jranke@uni-bremen.de

 

 # This file is part of the R package mkin

@@ -71,8 +71,9 @@ transform_odeparms <- function(parms, mkinmod,
   }

 

   # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2

-  # and HS parameter tb if transformation of rates is requested

-  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

+  # and HS parameter tb as well as logistic model parameters kmax, k0 and r if

+  # transformation of rates is requested

+  for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {

     if (!is.na(parms[pname])) {

       if (transform_rates) {

         transparms[paste0("log_", pname)] <- log(parms[pname])

@@ -151,8 +152,8 @@ backtransform_odeparms <- function(transparms, mkinmod,
     }

   }

 

-  # Transform parameters also for FOMC, DFOP and HS models

-  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

+  # Transform parameters also for FOMC, DFOP, HS and logistic models

+  for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {

     if (transform_rates) {

       pname_trans = paste0("log_", pname)

       if (!is.na(transparms[pname_trans])) {