from “Interpretation of thyroid-relevant bioactivity data for comparison to in vivo exposures: A prioritization approach for putative chemical inhibitors of in vitro deiodinase activity”
Kimberly T Truong, John F Wambaugh, Dustin F Kapraun, Sarah E Davidson-Fritz, Stephanie Eytcheson, Richard S Judson, and Katie Paul Friedman
Toxicological Sciences, 2025
In this vignette, we show how to use the full human gestational model from httk to perform in vitro-in vivo extrapolation (IVIVE) on a relevant in vitro bioactivity dataset. Briefly, this demonstrates how to take thyroid-related in vitro bioactivity data from ToxCast, perform IVIVE using the 3-compartment steady state model and the full human pregnancy model within httk, and compare it with the latest exposure estimates from ExpoCast as a method of prioritization. This is detailed in Truong et al. 2025 (in review).
Truong et al. 2025 describes a prioritization workflow for thyroid-related bioactivity data related to deiodinase inhibition by applying the full human gestational model to compute administered equivalent doses (AEDs) and bioactivity-to-exposure ratios (BERs) by integrating SEEM3 ExpoCast predictions. See Figure 1 for an overview of the workflow.
This vignette reproduces Figures 9 and 10 from Truong et al. 2025.
This vignette was created with httk v2.6.0 which includes new model functionality (see Simulate full-term human pregnancy model).
rm(list = ls())Several chunks of code in this vignette have “eval = execute.vignette” and we set execute.vignette to be FALSE, by default. We do this to include all the code that was necessary to produce the figures but was not run when the vignette was built since some of these chunks take a considerable amount of time and the checks on CRAN limit how long we can spend building the package. You can still reproduce the end figures knitting the .Rmd file, since all the necessary and minimal data objects have been stored in httk. If you want to rerun the whole workflow on your own, you must rerun all code from top to bottom, either by changing execute.vignette to TRUE or clicking on the “green arrow” in RStudio.
execute.vignette <- FALSEWe need to load in the processed in vitro bioactivity data.
This is stored within httk and can be accessed by
httk::thyroid.ac50s. From the R command prompt,
?thyroid.ac50s for more details.
ac50.dt <- httk::thyroid.ac50slibrary(httk)
library(data.table)
library(dplyr)
library(openxlsx)
library(ggpubr)
library(cowplot)
library(ggrepel)
library(latex2exp)
library(viridis)
library(scales)
library(ggstar)Only 50 chemicals have experimental toxicokinetic (TK) parameters sufficient for httk.
# how many chems had experimental TK params?
human.data.pre <- subset(get_cheminfo(info = 'all', species = 'Human',
model = '3compartmentss'),
DTXSID %in% ac50.dt$dsstox_substance_id)
nrow(human.data.pre)We need to supplement these with in silico predictions of TK parameters.
# predict Clint/fup values for missing chems
# this loads new chem.physical_and_invitro.data table into the global env
load_dawson2021() # most data for ToxCast chems
load_sipes2017() # most data for pharma compounds
load_pradeep2020() # best ML method for in silico prediction human.httk.dt <- subset(ac50.dt, dsstox_substance_id %in% get_cheminfo(info = 'DTXSID',
species = 'Human',
model = '3compartmentss',
physchem.exclude = FALSE,
median.only = TRUE))
missing.dtxsids <- setdiff(ac50.dt$dsstox_substance_id, human.httk.dt$dsstox_substance_id)
View(ac50.dt[dsstox_substance_id %in% missing.dtxsids])There are two ways to compute an oral AED in httk:
1. compute a steady-state plasma concentration via
calc_mc_css and then: \[
AED = \frac{AC_{50}}{CSS_{50}}
\] 2. use calc_mc_oral_equiv passing the AC50 to the
conc argument and specifying
which.quantile = c(0.5).
for (i in 1:nrow(human.httk.dt)) {
message('Calculating for chemical: ', human.httk.dt$chnm[i], '\n')
set.seed(123)
out <- calc_mc_css(dtxsid = human.httk.dt$dsstox_substance_id[i],
species = 'Human',
daily.dose = 1,
which.quantile = c(0.5, 0.95),
model = '3compartmentss',
output.units = 'uM',
parameterize.args.list=list(physchem.exclude=FALSE),
suppress.messages = TRUE)
human.httk.dt$mc.css50[i] <- out["50%"]
set.seed(123)
oralequiv.out <- calc_mc_oral_equiv(conc = human.httk.dt$min_ac50[i],
dtxsid = human.httk.dt$dsstox_substance_id[i],
species = 'Human',
which.quantile = c(0.5, 0.95),
input.units = 'uM',
output.units = 'mgpkgpday',
model = '3compartmentss',
restrictive.clearance = TRUE,
parameterize.args.list=list(physchem.exclude=FALSE),
suppress.messages = TRUE)
human.httk.dt$oral_equiv50[i] <- oralequiv.out["50%"]
}
human.httk.dt[, calc.aed50 := min_ac50/mc.css50]
human.httk.dt[, fold_diff50 := log10(oral_equiv50/calc.aed50)]
# let's use the aed values from the division
human.httk.dt2 <- human.httk.dt[, c('dsstox_substance_id', 'chnm',
'DIO1', 'DIO2', 'DIO3', 'IYD', 'min_ac50',
'mc.css50', 'calc.aed50')]
SEEM <- httk::truong25.seem3
# integrate SEEM3 predictions with merge
human.httk.dt2.seem3 <- merge.data.table(human.httk.dt2,
SEEM[, c("dsstox_substance_id", "seem3", "seem3.u95")],
by = c("dsstox_substance_id"),
all.x = TRUE)
human.httk.dt2.seem3[, plasmaBER50 := calc.aed50/seem3.u95]
human.httk.dt2.seem3[, log.pBER50 := log10(plasmaBER50)]
ivive.moe.tb <- human.httk.dt2.seem3
setnames(ivive.moe.tb, "dsstox_substance_id", "dtxsid")for (i in 1:nrow(ivive.moe.tb)) {
ivive.moe.tb$half.life[i] <- calc_half_life(dtxsid = ivive.moe.tb$dtxsid[i], # t1/2 in hours
physchem.exclude = FALSE)
ivive.moe.tb$half.life.days[i] <- ivive.moe.tb$half.life[i]/24 # in days
}
# tissues where DIO enzymes live
impacted_tissues <- c('Cliver', 'Cthyroid', 'Cfliver', 'Cfbrain', 'Cfthyroid',
'Cplacenta', 'Cconceptus', 'Cplasma')
targets <- c('DIO1', 'DIO2', 'DIO3', 'IYD')
# max conc in every tissue for every chemical
Cmax.tissues <- data.frame(matrix(NA, ncol = length(impacted_tissues) + 1,
nrow = nrow(ivive.moe.tb),
dimnames = list(NULL, c('dtxsid', impacted_tissues))))
for(i in 1:nrow(ivive.moe.tb)) {
sol.out <- solve_full_pregnancy(dtxsid = ivive.moe.tb$dtxsid[i],
daily.dose = 1,
doses.per.day = 1,
time.course = seq(0, 40*7, 1/24), # view solution by hour
track.vars = impacted_tissues,
physchem.exclude = FALSE)
# get max of every column re: compt
Cmax.tissues[i, impacted_tissues] <- apply(sol.out[, impacted_tissues], 2, max, na.rm = T)
Cmax.tissues[i, 'dtxsid'] <- ivive.moe.tb$dtxsid[i]
}tissue.aeds <- list()
# tissues specific for each target
tissue_list <- list(DIO1 = c('Cliver', 'Cfliver', 'Cthyroid', 'Cfthyroid', 'Cconceptus'),
DIO2 = c('Cthyroid', 'Cfthyroid', 'Cfbrain', 'Cconceptus'),
DIO3 = c('Cthyroid', 'Cfthyroid', 'Cfbrain', 'Cplacenta', 'Cconceptus'),
IYD = c('Cthyroid', 'Cfthyroid', 'Cconceptus'))
for (i in names(tissue_list)) {
active.chem.ind <- which(!is.na(ivive.moe.tb[, ..i]))
n.chems <- length(active.chem.ind)
empty.mat <- matrix(NA, nrow = n.chems, ncol = length(tissue_list[[i]]) + 1)
df.target <- data.frame(empty.mat)
colnames(df.target) <- c('dtxsid', tissue_list[[i]])
df.target[, 2: (length(tissue_list[[i]])+1) ] <- mapply(function(x,y) x/y, ivive.moe.tb[active.chem.ind, ..i], Cmax.tissues[active.chem.ind, tissue_list[[i]]])
df.target$dtxsid <- ivive.moe.tb$dtxsid[active.chem.ind]
tissue.aeds[[i]] <- df.target
}impacted_tissues <- c('Cliver', 'Cthyroid', 'Cfliver', 'Cfbrain', 'Cfthyroid',
'Cplacenta', 'Cconceptus', 'Cplasma')
get_ber_data <- function(name, df) {
df$exposure <- ivive.moe.tb$seem3.u95[match(df$dtxsid, ivive.moe.tb$dtxsid)]
df.m <- df %>% select(-exposure) %>%
reshape2::melt(id.vars = c('dtxsid'),
variable.name = 'compt', value.name = 'aed')
df.m$exposure <- ivive.moe.tb$seem3.u95[match(df.m$dtxsid, ivive.moe.tb$dtxsid)]
df.m$ber = df.m$aed/df.m$exposure
df.m2 <- df.m %>%
select(-exposure) %>%
reshape2::melt(id.vars = c('dtxsid', 'compt', 'ber'))
# add melted ExpoCast data as rows for each unique chem
expocast.dat <- df %>% select(dtxsid, exposure) %>%
reshape2::melt(id.vars = c('dtxsid'))
df.m2 <- bind_rows(df.m2, expocast.dat)
df.m2$value <- log10(df.m2$value)
df.m2$ber <- log10(df.m2$ber)
df.m2$chnm <- ivive.moe.tb$chnm[match(df.m2$dtxsid, ivive.moe.tb$dtxsid)]
df.m2$enzyme <- name
return(df.m2)
}
new.list <- Map(get_ber_data, names(tissue.aeds), tissue.aeds)
ggdata = bind_rows(new.list)
setDT(ggdata)
# order chemicals within each enzyme group by the most sensitive tissue/lifestage BER
ggdata[, min_ber := min(ber, na.rm = TRUE), keyby = .(enzyme, dtxsid)]
ggdata <- ggdata[order(min_ber), .SD, keyby = .(enzyme)]
# exposure values are dependent on dtxsid not httk compt
ggdata[variable == 'exposure', compt := 'exposure']
# make a deep copy or else it will be modified by reference
tissue.bers <- copy(ggdata)
# round all numerical values including BERs to 2 sigfigs
num.cols <- c("value", "ber", "min_ber")
tissue.bers[, c(num.cols) := lapply(.SD, signif, digits = 2), .SDcols = num.cols]min_val <- tissue.bers[min_ber <= 3, min(value), by = enzyme][, min(V1)]
max_val <- tissue.bers[min_ber <= 3, max(value), by = enzyme][, max(V1)]
ylims <- c(min_val, max_val)
# generate a set of 3 viridis colors for lifestage
# yellow represents neither maternal or fetal, i.e., conceptus
my_viridis_colors <- viridis(n = 3)
listgg <- list()
for(k in names(tissue.aeds)) {
aeds <- tissue.bers[enzyme == k & ber <= 3]
chems <- tissue.bers[enzyme == k & ber <= 3, dtxsid]
exposures <- tissue.bers[enzyme == k & variable == "exposure" & dtxsid %in% chems]
listgg[[k]] <- ggplot(data = rbind(aeds,exposures),
mapping = aes(x = reorder(chnm, -min_ber),
y = value)) +
geom_star(aes(starshape = compt, fill = compt, color = compt),
alpha = 0.65, size = 1.25, # default: 1.5
starstroke = 0.25, # default: 0.5
show.legend = T) +
scale_y_continuous(breaks = seq(ylims[1], ylims[2], 1),
limits = ylims) +
scale_starshape_manual(name = "Compartment corresponding to AED \n or Exposure",
labels = c("maternal liver", "fetal liver",
"maternal thyroid", "fetal thyroid",
"conceptus",
"fetal brain",
"placenta", "exposure"),
values = c('Cfbrain' = 15,
'exposure' = 28,
'Cliver' = 11,
'Cfliver' = 11,
'Cplacenta' = 13,
'Cthyroid' = 1,
'Cfthyroid' = 1,
'Cconceptus' = 23),
drop = F) + # this makes sure all levels of a factor are displayed even if unused
scale_fill_manual(name = "Compartment corresponding to AED \n or Exposure",
labels = c("maternal liver", "fetal liver",
"maternal thyroid", "fetal thyroid",
"conceptus",
"fetal brain",
"placenta", "exposure"),
values = c('Cfbrain' = my_viridis_colors[2],
'exposure' = 'dimgrey',
'Cliver' = my_viridis_colors[1],
'Cfliver' = my_viridis_colors[2],
'Cplacenta' = '#990066',
'Cthyroid' = my_viridis_colors[1],
'Cfthyroid' = my_viridis_colors[2],
'Cconceptus' = '#73D055FF'),
drop = F) +
scale_color_manual(name = "Compartment corresponding to AED \n or Exposure",
labels = c("maternal liver", "fetal liver",
"maternal thyroid", "fetal thyroid",
"conceptus",
"fetal brain",
"placenta", "exposure"),
values = c('Cfbrain' = my_viridis_colors[2],
'exposure' = 'dimgrey',
'Cliver' = my_viridis_colors[1],
'Cfliver' = my_viridis_colors[2],
'Cplacenta' = '#990066',
'Cthyroid' = my_viridis_colors[1],
'Cfthyroid' = my_viridis_colors[2],
'Cconceptus' = '#73D055FF'),
drop = F) +
labs(x = 'Chemical', y = 'Oral AED or Exposure (log10 mg/kg-bw/day)', title = k) +
theme_bw() +
theme(legend.position = "none",
plot.title = element_text(size = 8, face = "bold"),
plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"),
axis.title = element_text(size = 6),
axis.text.y = element_text(size = 5),
axis.text.x = element_text(size = 6),
legend.text = element_text(size = 6),
legend.title = element_text(size = 6, hjust = 0.5))+
coord_flip()
}
# remove extra axis labels
listgg$IYD <- listgg$IYD + labs(x = '', y = '')
listgg$DIO3 <- listgg$DIO3 +
labs(x = '')
listgg$DIO1 <- listgg$DIO1 + labs(y = '')
# arrange DIO1 and DIO2 with equal sizes since they have ~20 chems each (21 vs 24)
pcol1 <- plot_grid(listgg$DIO1, listgg$DIO2,
align = "hv", nrow = 2,
rel_heights = c(1,1.1))
# make DIO3 3x taller than IYD (44 chems vs 4)
pcol2 <- plot_grid(listgg$IYD, listgg$DIO3,
align = "hv", nrow = 2,
rel_heights = c(1,3))
# combine two columns of plots
prow <- plot_grid(pcol1, pcol2,
align = "h", ncol = 2,
rel_widths = c(1,1))
# extract legend from one of the plots
grobs <- ggplotGrob(
listgg$DIO2 +
guides(
starshape = guide_legend(ncol = 2, byrow = F,
keyheight = unit(0.25, "cm")),
fill = guide_legend(ncol = 2, byrow = F,
keyheight = unit(0.25, "cm")),
color = guide_legend(ncol = 2, byrow = F,
keyheight = unit(0.25, "cm")),
) +
theme(legend.position = "bottom")
)$grobs
legend2 <- grobs[[which(sapply(grobs, function(x) x$name) == "guide-box")]]
# add the legend underneath the plots; make it 10% of the height of the combined plots
fig9 <- plot_grid(prow, legend2,
nrow = 2, rel_heights = c(0.88,0.12)) +
theme(plot.background = element_rect(fill = "white", color = "white"))
# ggsave(plot = fig9,
# units = "mm",
# dpi = 300,
# width = 190, height = 150,
# device = "png",
# filename = "./figures/preg_prioritization_by_BER.png")targets <- c('DIO1', 'DIO2', 'DIO3', 'IYD')
# DATA WRANGLING
tissue.bers[min_ber == ber, most_sensitive_tissue := compt]
# table matching order of chemicals in enzyme subplots with all calculated BERs
# every enzyme x chem x min BER constitutes a row
minber.by.chem <- tissue.bers[!is.na(most_sensitive_tissue),
.(chnm, min_ber, most_sensitive_tissue), by = .(enzyme, dtxsid)]
# check if there are multiple most-sensitive tissues per chem in each enzyme group
minber.by.chem[, if(.N > 1) .SD, by = .(enzyme, dtxsid)][, length(unique(dtxsid))] #70 chems
# collapse multiple most-sensitive tissues per chem in each enzyme group
min.ber.reduced <- minber.by.chem[, most_sensitive_compts := paste(most_sensitive_tissue, collapse = ", "), by = .(enzyme, dtxsid)][, .(chnm, min_ber, most_sensitive_compts), by = .(enzyme, dtxsid)]
min.ber.reduced <- unique(min.ber.reduced)
# min BER matrix
min.ber.mat <- data.table::dcast(min.ber.reduced[, 1:4],
dtxsid + chnm ~ enzyme, value.var = c('min_ber'))
min.ber.mat[, min_across_enzyme := pmin(DIO1, DIO2, DIO3, IYD, na.rm = T)]
min.ber.mat <- min.ber.mat[order(min_across_enzyme)]
# find the target corresponding to the min_across_enzyme
min.ber.mat[, enzyme := mapply(function(x) targets[x], apply(min.ber.mat[, targets, with = FALSE], 1, which.min))]
# merge data on most sensitive tissues reflecting min BER
new.ranking.tissues <- merge.data.table(min.ber.mat,
minber.by.chem[, .(dtxsid, chnm, enzyme, min_ber, most_sensitive_tissue)],
by = c('dtxsid', 'chnm', 'enzyme'),
all.x = TRUE)
setnames(new.ranking.tissues, old = 'enzyme', new = 'target_min')
# round all BERs to 2 sigfigs for easy table viewing
ranked.by.plasma <- ivive.moe.tb[, lmoe50 := signif(log.pBER50, digits = 2)][order(lmoe50)]
final.new.ranking <- new.ranking.tissues[, -c("target_min", "min_ber")]
final.new.ranking <- final.new.ranking[order(min_across_enzyme)]
# comparing tissue BERs from full preg model with plasma BERs from 3compss model
plasma.tissue.bers <- merge.data.table(ranked.by.plasma[, c('dtxsid', 'chnm', 'lmoe50')],
final.new.ranking[, -c('chnm')],
by = c('dtxsid'))
setnames(plasma.tissue.bers, old = colnames(plasma.tissue.bers)[3:8],
new = c('plasma_BER50', 'DIO1_BER', 'DIO2_BER', 'DIO3_BER', 'IYD_BER', 'min_BER'))
plasma.tissue.bers <- plasma.tissue.bers[order(min_BER)]# merge min BER with plasma BER for each chem with unique row for each chem
# ignoring multiple sensitive tissues
ber.comp <- merge.data.table(ranked.by.plasma[, c('dtxsid', 'chnm', 'lmoe50')],
min.ber.mat[, -c("chnm")],
by = c("dtxsid"))
setnames(ber.comp, old = colnames(ber.comp)[3:8],
new = c('plasma_BER50', 'DIO1_BER', 'DIO2_BER', 'DIO3_BER', 'IYD_BER', 'min_BER'))
max_ber <- max(ber.comp[, c('plasma_BER50', 'min_BER')], na.rm = T)
min_ber <- min(ber.comp[, c('plasma_BER50', 'min_BER')], na.rm = T)
ber.lims <- c(floor(min_ber), ceiling(max_ber))
# chems with labels
chem.labels <- ber.comp[min_BER < 0][order(min_BER)][1:6, chnm]
# theme for the next 2 figures
my_theme <- theme_bw() +
theme(axis.title = element_text(size = 8, face = "bold"),
axis.text = element_text(size = 7))
preg.vs.nonpreg.pp <- ggplot(data = ber.comp, aes(x = plasma_BER50, y = min_BER)) +
geom_point(alpha = 0.75, size = 1,
show.legend = T) +
geom_abline(slope = 1, linewidth = 0.5, linetype = 'solid', color = "#666666") +
geom_abline(aes(slope = 1, intercept = 1), linetype = 'longdash', color = 'grey', linewidth = 0.5) + # default linesize = 1
geom_abline(aes(slope = 1, intercept = -1), linetype = 'longdash', color = 'grey', linewidth = 0.5) +
geom_label_repel(aes(label = chnm),
data = ber.comp[chnm %in% chem.labels],
max.overlaps = Inf,
force = 50,
size = 2,
label.padding = 0.1,
label.size = 0.1,
label.r = 0.08,
segment.size = 0.2,
min.segment.length = 0) +
my_theme +
scale_x_continuous(breaks = seq(ber.lims[1], ber.lims[2], 1),
limits = ber.lims) +
scale_y_continuous(breaks = seq(ber.lims[1], ber.lims[2], 1),
limits = ber.lims) +
labs(x = 'plasma BER using 3compss model', y = 'min BER across DIO/IYD targets and tissues')
# how many chems to the right of y = x?
ber.comp[plasma_BER50 > min_BER, .N] #89 chems
# plot SEEM3 uncertainty
ivive.moe.tb$uncertainty = log10(ivive.moe.tb$seem3.u95/ivive.moe.tb$seem3)
ber.comp$uncertainty <- ivive.moe.tb$uncertainty[match(ber.comp$dtxsid, ivive.moe.tb$dtxsid)]
min.min.ber <- min(ber.comp$min_BER)
max.min.ber <- max(ber.comp$min_BER)
min.ber.lims <- c(floor(min.min.ber), ceiling(max.min.ber))
max.uncertainty <- max(ber.comp$uncertainty)
seem3pp <- ggplot(data = ber.comp,
mapping = aes(x = min_BER,
y = uncertainty)) +
geom_point(size = 1, alpha = 0.8) +
geom_label_repel(aes(label = chnm),
data = subset(ber.comp, chnm %in% chem.labels),
max.overlaps = Inf,
force = 75,
size = 2,
label.padding = 0.1,
label.size = 0.1,
label.r = 0.08,
segment.size = 0.2,
min.segment.length = 0) +
my_theme +
scale_x_continuous(breaks = seq(min.ber.lims[1], min.ber.lims[2], 1),
limits = min.ber.lims) +
scale_y_continuous(breaks = seq(0, ceiling(max.uncertainty), 1),
limits = c(0, ceiling(max.uncertainty))) +
labs(x = 'min BER by Chemical', y = TeX(r'($log_{10}ExpoCast_{u95} - log_{10}ExpoCast_{med}$)'))
median(ber.comp$min_BER)
#> [1] 2.9
fig10 <- plot_grid(preg.vs.nonpreg.pp, seem3pp,
labels = "AUTO", label_size = 10)
# ggsave(plot = fig10,
# units = "mm",
# dpi = 300,
# width = 190, height = 100,
# device = "png",
# filename = "./figures/BER_uncertainty.png")thyroid.ac50s <- ac50.dt
truong25.seem3 <- SEEM
save(thyroid.ac50s, truong25.seem3,
file="./httk/Truong2025Vignette.RData")