1
#' Function to calculate endpoints for further use from kinetic models fitted
2
#' with mkinfit
3
#'
4
#' This function calculates DT50 and DT90 values as well as formation fractions
5
#' from kinetic models fitted with mkinfit. If the SFORB model was specified
6
#' for one of the parents or metabolites, the Eigenvalues are returned. These
7
#' are equivalent to the rate constants of the DFOP model, but with the
8
#' advantage that the SFORB model can also be used for metabolites.
9
#'
10
#' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
11
#' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
12
#'
13
#' @param fit An object of class [mkinfit], [nlme.mmkin] or
14
#'  [saem.mmkin]. Or another object that has list components
15
#'  mkinmod containing an [mkinmod] degradation model, and two numeric vectors,
16
#'  bparms.optim and bparms.fixed, that contain parameter values
17
#'  for that model.
18
#' @importFrom stats optimize
19
#' @return A list with a matrix of dissipation times named distimes,
20
#'   and, if applicable, a vector of formation fractions named ff
21
#'   and, if the SFORB model was in use, a vector of eigenvalues
22
#'   of these SFORB models, equivalent to DFOP rate constants
23
#' @note The function is used internally by [summary.mkinfit],
24
#'   [summary.nlme.mmkin] and [summary.saem.mmkin].
25
#' @author Johannes Ranke
26
#' @examples
27
#'
28
#'   fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
29
#'   endpoints(fit)
30
#'   \dontrun{
31
#'     fit_2 <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE)
32
#'     endpoints(fit_2)
33
#'     fit_3 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE)
34
#'     endpoints(fit_3)
35
#'   }
36
#'
37
#' @export
38
endpoints <- function(fit) {
39 2
  ep <- list()
40 2
  mkinmod <- fit$mkinmod
41 2
  degparms <- c(fit$bparms.optim, fit$bparms.fixed)
42 2
  obs_vars <- names(mkinmod$spec)
43 2
  ep$ff <- vector()
44 2
  ep$SFORB <- vector()
45 2
  ep$distimes <- data.frame(
46 2
    DT50 = rep(NA, length(obs_vars)),
47 2
    DT90 = rep(NA, length(obs_vars)),
48 2
    row.names = obs_vars)
49 2
  for (obs_var in obs_vars) {
50 2
    type = names(mkinmod$map[[obs_var]])[1]
51

52
    # Get formation fractions if directly fitted, and calculate remaining fraction to sink
53 2
    f_names = grep(paste("^f", obs_var, sep = "_"), names(degparms), value=TRUE)
54 2
    if (length(f_names) > 0) {
55 2
      f_values = degparms[f_names]
56 2
      f_to_sink = 1 - sum(f_values)
57 2
      names(f_to_sink) = ifelse(type == "SFORB",
58 2
        paste(obs_var, "free", "sink", sep = "_"),
59 2
        paste(obs_var, "sink", sep = "_"))
60 2
      for (f_name in f_names) {
61 2
        ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
62
      }
63 2
      ep$ff = append(ep$ff, f_to_sink)
64
    }
65

66
    # Get the rest
67 2
    if (type == "SFO") {
68 2
      k_names = grep(paste("^k", obs_var, sep="_"), names(degparms), value=TRUE)
69 2
      k_tot = sum(degparms[k_names])
70 2
      DT50 = log(2)/k_tot
71 2
      DT90 = log(10)/k_tot
72 2
      if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) {
73 2
        for (k_name in k_names)
74
        {
75 2
          ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
76
        }
77
      }
78
    }
79 2
    if (type == "FOMC") {
80 2
      alpha = degparms["alpha"]
81 2
      beta = degparms["beta"]
82 2
      DT50 = beta * (2^(1/alpha) - 1)
83 2
      DT90 = beta * (10^(1/alpha) - 1)
84 2
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
85 2
      ep$distimes[obs_var, c("DT50back")] = DT50_back
86
    }
87 2
    if (type == "IORE") {
88 2
      k_names = grep(paste("^k__iore", obs_var, sep="_"), names(degparms), value=TRUE)
89 2
      k_tot = sum(degparms[k_names])
90
      # From the NAFTA kinetics guidance, p. 5
91 2
      n = degparms[paste("N", obs_var, sep = "_")]
92 2
      k = k_tot
93
      # Use the initial concentration of the parent compound
94 2
      source_name = mkinmod$map[[1]][[1]]
95 2
      c0 = degparms[paste(source_name, "0", sep = "_")]
96 2
      alpha = 1 / (n - 1)
97 2
      beta = (c0^(1 - n))/(k * (n - 1))
98 2
      DT50 = beta * (2^(1/alpha) - 1)
99 2
      DT90 = beta * (10^(1/alpha) - 1)
100 2
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
101 2
      ep$distimes[obs_var, c("DT50back")] = DT50_back
102 2
      if (mkinmod$use_of_ff == "min") {
103 0
        for (k_name in k_names)
104
        {
105 0
          ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
106
        }
107
      }
108
    }
109 2
    if (type == "DFOP") {
110 2
      k1 = degparms["k1"]
111 2
      k2 = degparms["k2"]
112 2
      g = degparms["g"]
113 2
      f <- function(log_t, x) {
114 2
        t <- exp(log_t)
115 2
        fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
116 2
        (fraction - (1 - x/100))^2
117
      }
118 2
      DT50_k1 = log(2)/k1
119 2
      DT50_k2 = log(2)/k2
120 2
      DT90_k1 = log(10)/k1
121 2
      DT90_k2 = log(10)/k2
122

123 2
      DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum),
124 2
                  silent = TRUE)
125 2
      DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum),
126 2
                  silent = TRUE)
127 0
      if (inherits(DT50, "try-error")) DT50 = NA
128 0
      if (inherits(DT90, "try-error")) DT90 = NA
129 2
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
130

131 2
      ep$distimes[obs_var, c("DT50back")] = DT50_back
132 2
      ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
133 2
      ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
134
    }
135 2
    if (type == "HS") {
136 2
      k1 = degparms["k1"]
137 2
      k2 = degparms["k2"]
138 2
      tb = degparms["tb"]
139 2
      DTx <- function(x) {
140 2
        DTx.a <- (log(100/(100 - x)))/k1
141 2
        DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
142 2
        if (DTx.a < tb) DTx <- DTx.a
143 2
        else DTx <- DTx.b
144 2
        return(DTx)
145
      }
146 2
      DT50 <- DTx(50)
147 2
      DT90 <- DTx(90)
148 2
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
149 2
      DT50_k1 = log(2)/k1
150 2
      DT50_k2 = log(2)/k2
151 2
      ep$distimes[obs_var, c("DT50back")] = DT50_back
152 2
      ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
153 2
      ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
154
    }
155 2
    if (type == "SFORB") {
156
      # FOCUS kinetics (2006), p. 60 f
157 2
      k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(degparms), value=TRUE)
158 2
      k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
159 2
      k_1output = sum(degparms[k_out_names])
160 2
      k_12 = degparms[paste("k", obs_var, "free", "bound", sep="_")]
161 2
      k_21 = degparms[paste("k", obs_var, "bound", "free", sep="_")]
162

163 2
      sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
164 2
      b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
165 2
      b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
166

167 2
      DT50_b1 = log(2)/b1
168 2
      DT50_b2 = log(2)/b2
169 2
      DT90_b1 = log(10)/b1
170 2
      DT90_b2 = log(10)/b2
171

172 2
      SFORB_fraction = function(t) {
173 2
        ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
174 2
        ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
175
      }
176

177 2
      f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2
178 2
      log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum,
179 2
                      silent = TRUE)
180 2
      f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2
181 2
      log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum,
182 2
                      silent = TRUE)
183

184 2
      DT50 = if (inherits(log_DT50, "try-error")) NA
185 2
             else exp(log_DT50)
186 2
      DT90 = if (inherits(log_DT90, "try-error")) NA
187 2
             else exp(log_DT90)
188

189 2
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
190

191 2
      for (k_out_name in k_out_names)
192
      {
193 2
        ep$ff[[sub("k_", "", k_out_name)]] = degparms[[k_out_name]] / k_1output
194
      }
195

196
      # Return the eigenvalues for comparison with DFOP rate constants
197 2
      ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
198 2
      ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
199

200 2
      ep$distimes[obs_var, c("DT50back")] = DT50_back
201 2
      ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1
202 2
      ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2
203
    }
204 2
    if (type == "logistic") {
205
      # FOCUS kinetics (2014) p. 67
206 2
      kmax = degparms["kmax"]
207 2
      k0 = degparms["k0"]
208 2
      r = degparms["r"]
209 2
      DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
210 2
      DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))
211

212 2
      DT50_k0 = log(2)/k0
213 2
      DT50_kmax = log(2)/kmax
214 2
      ep$distimes[obs_var, c("DT50_k0")] = DT50_k0
215 2
      ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax
216
    }
217 2
    ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
218
  }
219 2
  if (length(ep$ff) == 0) ep$ff <- NULL
220 2
  if (length(ep$SFORB) == 0) ep$SFORB <- NULL
221 2
  return(ep)
222
}

Read our documentation on viewing source code .

Loading