-
Notifications
You must be signed in to change notification settings - Fork 1
/
synergyTheory.R
262 lines (226 loc) · 11 KB
/
synergyTheory.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
# Copyright: (C) 2017-2018 Sachs Undergraduate Research Apprentice Program
# This program and its accompanying materials are distributed
# under the terms of the GNU General Public License v3.
# Filename: synergyTheory.R
# Purpose: Concerns radiogenic mouse Harderian gland tumorigenesis.
# Contains relevant synergy theory models, information coefficient
# calculations and useful objects. It is part of the
# source code for the NASAmouseHG project.
# Contact: Rainer K. Sachs
# Website: https://github.com/sachsURAP/NASAmouseHG
# Mod history: 18 Jun 2018
# Details: See dataAndInfo.R for further licensing, attribution,
# references, and abbreviation information.
source("dataAndInfo.R") # Load in the data. Remark: dose is in units of cGy;
# LET usually in keV/micron; prevalence Prev always < 1
# (i.e. not in %, which would mean prevalence < 100 but is strongly deprecated).
library(deSolve) # Solving differential equations.
#========================= MISC. OBJECTS & VARIABLES ==========================#
# In next line phi controls how fast NTE build up from zero; not really needed
# during calibration since phi * Dose >> 1 at every observed Dose !=0. phi is
# needed for later synergy calculations.
# d_0 = 1 / phi = 5 x 10-4 cGy = 5 x 10^-6 Gy.
phi <- 2000 # even larger phi should give the same final results,
# but might cause extra problems with R.
#================================ DER MODELS ==================================#
#================ PHOTON MODEL =================#
# Linear model fit on beta_decay_data dataset.
# We will never recalculate this unless new data comes in but here it is
# just in case.
# beta_decay_lm <- lm(HG ~ dose + I(dose ^ 2), data = beta_decay_data)
# summary(beta_decay_lm, correlation = TRUE)
#=============== HZE/NTE MODEL =================#
# (HZE = high charge and energy;
# NTE = non-targeted effects are included in addition to TE)
HZE_data <- ion_data[13:47, ] # Includes 1-ion data iff Z > 3
# Uses 3 adjustable parameters.
HZE_nte_model <- nls( # Calibrating parameters in a model that modifies the hazard function NTE models in 17Cuc.
Prev ~ .0275 + (1 - exp ( - (aa1 * LET * dose * exp( - aa2 * LET)
+ (1 - exp( - phi * dose)) * kk1))),
data = HZE_data,
weights = NWeight,
start = list(aa1 = .00009, aa2 = .001, kk1 = .06)) # Use extra argument trace = TRUE if you want to watch convergence.
summary(HZE_nte_model, correlation = TRUE) # Parameter values & accuracy.
# If a paper uses dose in Gy some care is required in the preceding and following lines to rescale from cGy.
vcov(HZE_nte_model) # Variance-covariance matrix.
HZE_nte_model_coef <- coef(HZE_nte_model) # Calibrated central values of the 3 parameters.
# The DER, = 0 at dose 0.
calibrated_nte_hazard_func <- function(dose, LET, coef) { # Calibrated hazard function.
return(coef[1] * LET * dose * exp( - coef[2] * LET)
+ (1 - exp( - phi * dose)) * coef[3])
}
calibrated_HZE_nte_der <- function(dose, LET, coef = HZE_nte_model_coef) { # Calibrated HZE NTE DER.
return(1 - exp( - calibrated_nte_hazard_func(dose, LET, coef)))
}
#================ HZE/TE MODEL =================#
# (TE = targeted effects only). This chunk runs and gives good results.
# We will not use it in the minor paper, only for later papers.
HZE_te_model <- nls( # Calibrating parameters in a TE only model.
Prev ~ .0275 + (1 - exp ( - (aate1 * LET * dose * exp( - aate2 * LET)))),
data = HZE_data,
weights = NWeight,
start = list(aate1 = .00009, aate2 = .01))
summary(HZE_te_model, correlation = TRUE) # Parameter values & accuracy.
vcov(HZE_te_model) # Variance-covariance matrix.
HZE_te_model_coef <- coef(HZE_te_model) # Calibrated central values of the 2 parameters.
# The DER, = 0 at dose 0.
calibrated_te_hazard_func <- function(dose, LET, coef) { # Calibrated hazard function.
return(coef[1] * LET * dose * exp( - coef[2] * LET))
}
calibrated_HZE_te_der <- function(dose, LET, coef = HZE_te_model_coef) {
return(1 - exp( - calibrated_te_hazard_func(dose, LET, coef))) # Calibrated HZE TE one-ion DER.
}
#==== LIGHT ION, LOW Z (<= 3), LOW LET MODEL ===#
low_LET_data = ion_data[1:12, ] # Swift protons and alpha particles.
low_LET_model <- nls(
Prev ~ .0275 + 1 - exp( - alpha_low * dose), # alpha is used throughout radioiology for dose coefficients.
data = low_LET_data,
weights = NWeight,
start = list(alpha_low = .005))
summary(low_LET_model, correlation = TRUE)
low_LET_model_coef <- coef(low_LET_model) # Calibrated central values of the parameter.
# Calibrated Low LET model. Use L = 0, but maybe later will use small L > 0.
calibrated_low_LET_der <- function(dose, LET, alph_low = low_LET_model_coef[1]) {
return(1 - exp( - alph_low * dose))
}
# Slope dE/dd of the low LET, low Z model.
low_LET_slope <- function(dose, LET) {
low_LET_model_coef * exp( - low_LET_model_coef * dose)
}
#=========================== INFORMATION CRITERION ============================#
info_crit_table <- cbind(AIC(HZE_te_model, HZE_nte_model),
BIC(HZE_te_model, HZE_nte_model))
print(info_crit_table)
#=========== BASELINE NO-SYNERGY/ANTAGONISM MIXTURE DER FUNCTIONS =============#
#======== SIMPLE EFFECT ADDITIVITY (SEA) =======#
#' @description Applies Simple Effect Additivity to get a baseline mixture DER.
#'
#' @param dose Numeric vector corresponding to the sum dose in cGy.
#' @param LET Numeric vector of all LET values, must be length n.
#' @param ratios Numeric vector of all dose ratios, must be length n.
#' @param lowLET Boolean of whether a low LET DER should be included in the mixture DER.
#' @param n Number of DERs, optional argument used to check parameter validity.
#'
#' @details Corresponding elements of ratios, LET should be associated with the
#' same DER.
#'
#' @return Numeric vector representing the estimated Harderian Gland
#' prevalence from a SEA mixture DER constructed from the given DER
#' parameters.
#'
#' @examples
#' calculate_SEA(.01 * 0:40, c(70, 195), c(1/2, 1/2), n = 2)
#' calculate_SEA(.01 * 0:70, c(0.4, 195), c(4/7, 3/7))
calculate_SEA <- function(dose, LET, ratios, lowLET = FALSE, n = NULL) {
if (!is.null(n) && (n != length(ratios) | n != length(LET))) {
stop("Length of arguments do not match.")
} else if (sum(ratios) != 1) {
stop("Sum of ratios do not add up to one.")
} # End error handling
total <- 0
i <- 1
if (lowLET == TRUE) {
# First elements of ratios and LET should be the low-LET DER
total <- total + calibrated_low_LET_der(dose * ratios[i], LET[i])
i <- i + 1
}
while (i < length(ratios) + 1) { # Iterate over HZE ions in the mixture.
total <- total + calibrated_HZE_nte_der(dose * ratios[i], LET[i])
i <- i + 1
}
return(total)
}
#====== INCREMENTAL EFFECT ADDITIVITY (IEA) ====#
#' @description Applies IEA to get a baseline no-synergy/antagonism mixture DER.
#'
#' @param dose Numeric vector corresponding to the sum dose in cGy.
#' @param LET Numeric vector of all LET values, must be length n.
#' @param ratios Numeric vector of all dose ratios, must be length n.
#' @param model String value corresponding to the model to be used, either "NTE" or "TE".
#' @param coef Named list of numeric vectors containing coefficients for one-ion DERs.
#' @param ders Named list of functions containing relevant one-ion DER models.
#' @param calculate_dI Named vector of functions to calculate dI depending on
#' the selected model.
#' @param phi Numeric value, used in NTE models.
#'
#' @details Corresponding elements of ratios, LET should be associated with the
#' same DER.
#'
#' @return Numeric vector representing the estimated Harderian Gland
#' prevalence from an IEA mixture DER constructed from the given
#' one-ion DERs parameters.
#'
#' @examples
#' calculate_id(.01 * 0:40, c(70, 195), c(1/2, 1/2))
#' calculate_id(.01 * 0:70, c(.4, 195), c(4/7, 3/7), model = "TE")
calculate_id <- function(dose, LET, ratios, model = "NTE",
coef = list(NTE = HZE_nte_model_coef,
TE = HZE_te_model_coef,
lowLET = low_LET_model_coef),
ders = list(NTE = calibrated_HZE_nte_der,
TE = calibrated_HZE_te_der,
lowLET = calibrated_low_LET_der),
calculate_dI = c(NTE = .calculate_dI_nte,
TE = .calculate_dI_te),
phi = 2000) {
dE <- function(yini, state, pars) { # Constructing an ODE from the DERS.
with(as.list(c(state, pars)), {
# Screen out low LET values.
lowLET_total <- lowLET_ratio <- 0
remove <- c()
for (i in 1:length(LET)) {
if (LET[i] <= 3) { # Ion is low-LET.
lowLET_total <- lowLET_total + LET[i]
lowLET_ratio <- lowLET_ratio + ratios[i]
remove <- unique(c(remove, LET[i]))
ratios[i] <- 0
}
}
LET <- c(setdiff(LET, remove))
ratios <- ratios[! ratios == 0]
# Begin calculating dI values.
aa <- u <- dI <- vector(length = length(LET))
if (length(LET) > 0) {
for (i in 1:length(LET)) {
aa[i] <- pars[1] * LET[i] * exp( - pars[2] * LET[i])
u[i] <- uniroot(function(dose) HZE_der(dose, LET[i], pars) - I,
interval = c(0, 20000),
extendInt = "yes",
tol = 10 ^ - 10)$root
dI[i] <- ratios[i] * calc_dI(aa[i], u[i], pars[3])
}
}
if (lowLET_ratio > 0) {
# If low-LET DER is present then include it at the end of the dI vector.
u[length(LET) + 1] <- uniroot(function(dose)
low_der(dose, LET = lowLET_total,
alph_low = coef[["lowLET"]]) - I,
interval = c(0, 20000),
extendInt = "yes",
tol = 10 ^ - 10)$root
dI[length(LET) + 1] <- lowLET_ratio * low_LET_slope(u[length(LET) + 1],
LET = lowLET_total)
}
return(list(sum(dI)))
})
}
p <- list(pars = coef[[model]],
HZE_der = ders[[model]],
low_der = ders[["lowLET"]],
calc_dI = calculate_dI[[model]])
return(ode(c(I = 0), times = dose, dE, parms = p, method = "radau"))
}
#============= dI HIDDEN FUNCTIONS =============#
.calculate_dI_nte <- function(aa, u, kk1) {
return((aa + exp( - phi * u) * kk1 * phi) *
exp( - (aa * u + (1 -exp( - phi * u)) * kk1)))
}
.calculate_dI_te <- function(aa, u, pars = NULL) {
return(aa * exp(- aa * u))
}
#============================ DEVELOPER FUNCTIONS =============================#
test_runtime <- function(f, ...) { # Naive runtime check
start_time <- Sys.time()
f(...)
Sys.time() - start_time
}