From d89e3d22eb9dc383897b09e9c5aa1b57f65cdbf0 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 21 Feb 2019 14:34:45 +0100 Subject: Add the logistic model --- R/endpoints.R | 8 ++++++++ R/logistic.solution.R | 4 ++++ R/mkinfit.R | 9 ++++++--- R/mkinmod.R | 20 +++++++++++++------- R/mkinpredict.R | 8 ++++++-- R/transform_odeparms.R | 11 ++++++----- 6 files changed, 43 insertions(+), 17 deletions(-) create mode 100644 R/logistic.solution.R (limited to 'R') 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])) { -- cgit v1.2.1