Introduction

Chemical risk assessment depends on knowledge of inherent chemical hazard, the dose-response relationship, and the extent of exposure that occurs (NASEM 2017). High throughput screening (HTS) for in vitro bioactivity allows characterization of potential hazard for thousands of chemicals for which no other testing has occurred (Judson et al., 2009).

Toxicokinetics (TK) describes the Absorption, Distribution, Metabolism, and Excretion (ADME) of a chemical by the body. TK relates external exposures to internal tissue concentrations of chemical and therefore informs the chemical dose-response relationship (Coecke et al., 2013). TK requires chemical-specific information that is traditionally obtained in vivo. However, most chemicals do not have TK data, so a new approach methodology (NAM) uses in vitro TK methods adapted from the pharmaceutical industry to provide the necessary chemical-specific information (Wambaugh et al., 2019, Chang et al., 2022). High(er) throughput toxicokinetics (HTTK) involves the combination of chemical-specific in vitro TK data and chemical-independent (“generic”) TK models to predict TK (Breen et al., 2021).

The primary goal of HTTK is to provide a human dose context for bioactive in vitro concentrations from HTS (that is, in vitro-in vivo extrapolation, or IVIVE) (Wetmore, 2015). IVIVE is Utilization of in vitro experimental data to predict phenomena in vivo. In pharmaceutical development, HTTK methods allow IVIVE to estimate therapeutic doses for clinical studies – predicted concentrations are typically on the order of values measured in clinical trials (Wang, 2010). High throughput screening + IVIVE can predict a daily chemical dose rate (mg/kg bw/day) that might be adverse (Paul Friedman et al., 2019).

A secondary goal is to provide open source data and models for evaluation and use by the scientific community (Pearce et al., 2017).

The following material is from a lecture taught as part of the course “Computational Methods in Chemical Risk Assessment” at Rutgers University in September, 2019.

knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
options(rmarkdown.html_vignette.check_title = FALSE)

rm(list=ls()) 

Background

If you are familiar with computer files and R you might want to skip directly to “Step 1: Loading In Vitro Bioactivity Data”.

install.packages("httk")

install.packages("readxl")

install.packages("ggplot2")

library(httk)
library(readxl)
library(ggplot2)

NSAMP <- 25

Step One: Loading In Vitro Screening Data

The ToxCast and Tox21 research programs employ batteries of high-throughput assays to assess chemical bioactivity in vitro. Not every chemical is tested through every assay (Richard et al., 2016). Most assays are conducted in concentration response, and each corresponding assay endpoint is analyzed statistically to determine if there is a concentration-dependent response or “hit” using the ToxCast Pipeline (Filer et al., 2017). Most assay endpoint-chemical combinations are nonresponsive. Here, only the hits are treated as potential indicators of bioactivity.

In vitro bioactivity does not necessarily indicate adversity or hazard. Biological models can be used to make predictions of toxicity based on in vitro data (for example, Sipes et al., 2011, Browne et al., 2015 and Kleinstreuer et al., 2017). However, here we make a more conservative, precautionary assumption that concentrations too low to cause in vitro bioactivity are more likely to be safe (Paul Friedman et al., 2019). Among all of the assay hits for each chemical, we choose to use the lower 10th percentile of the \(\mu M\) potencies (50% active concentration or \(AC_{50}\)). We are assuming that the 10th percentile represents a low concentration for activating multiple assays (assuming 10 or more bioactive assays were observed for a given chemical). However, this bioactivity does not necessarily have a direct toxicological interpretation.

setwd("FOLDERPATH")

load("invitrodb_3_5_mc5.Rdata")

toxcast.httk <- subset(mc5, dsstox_substance_id %in% get_cheminfo(
         info="DTXSID",
         suppress.messages=TRUE))

Critical Step: Select your chemicals of interest

We might have chemicals we are interested in for one reason or another – you could type any chemical ID’s you want into the following my.chems vector, or even load it from a file.

Here we’ll pick 50 chemicals at random from among the ToxCast chemicals:

set.seed(1234)
my.chems <- sample(mc5$dsstox_substance_id,50)
example.toxcast <- as.data.frame(mc5[mc5$dsstox_substance_id %in% my.chems,])

Unfortunately for this vignette there are too many ToxCast data to fit into a 5mb R package. So we will subset to just the selected chemicals and distribute only those data. In addition, out of 78 columns in the data, we will keep only eight. Download the full data following the instructions above.

example.toxcast <- example.toxcast[, c("chnm",
              "dsstox_substance_id",
              "spid",
              "hitc",
              "modl",
              "aeid",
              "modl_ga",
              "modl_ac10",
              "modl_acc")]
# reduce precision to decrease space:
cols <- c("modl_ga", "modl_ac10", "modl_acc")
for (this.col in cols)
  example.toxcast[, this.col] = signif(example.toxcast[, this.col], 3)
save(example.toxcast,file="introtoivive-toxcast.Rdata",version=2)

example.toxcast <- httk::example.toxcast
knitr::kable(head(example.toxcast), caption = "ToxCast In Vitro Bioactivity Data",
             floating.environment="sidewaystable")

toxcast.table <- NULL
for (this.id in unique(example.toxcast$dsstox_substance_id))
{
  this.subset <- subset(example.toxcast, dsstox_substance_id == this.id)
  these.hits <- subset(this.subset, hitc==1)
  if (dim(these.hits)[1]>0){
      this.row <- data.frame(Compound=as.character(this.subset[1,"chnm"]),
                         DTXSID=this.id,
                         Total.Assays = dim(this.subset)[1],
                         Unique.Assays = length(unique(this.subset$aeid)),
                         Total.Hits = dim(these.hits)[1],
                         Unique.Hits = length(unique(these.hits$aeid)),
                         Low.AC50 = signif(min(these.hits$modl_ga),3),
                         Low.AC10 = signif(min(these.hits$modl_ac10),3),
                         Low.ACC = signif(min(these.hits$modl_acc),3),
                         Q10.AC50 = signif(quantile(these.hits$modl_ga,probs=0.1),3),
                         Q10.AC10 = signif(quantile(these.hits$modl_ac10,probs=0.1),3),
                         Q10.ACC = signif(quantile(these.hits$modl_acc,probs=0.1),3),
                         Med.AC50 = signif(median(these.hits$modl_ga),3),
                         Med.AC10 = signif(median(these.hits$modl_ac10),3),
                         Med.ACC = signif(median(these.hits$modl_acc),3),
                         Q90.AC50 = signif(quantile(these.hits$modl_ga,probs=0.9),3),
                         Q90.AC10 = signif(quantile(these.hits$modl_ac10,probs=0.9),3),
                         Q90.ACC = signif(quantile(these.hits$modl_acc,probs=0.9),3)
                         )
    toxcast.table <- rbind(toxcast.table, this.row)
  }
}
rownames(toxcast.table) <- seq(1,dim(toxcast.table)[1])
knitr::kable(head(toxcast.table[,1:6]), caption = "Summarized ToxCast Data",
             floating.environment="sidewaystable")

Step Two: High Throughput Toxicokinetics from In Vitro Data

For each chemical in your subset of ToxCast, add the HTTK plasma steady-state concentration if it is available. calc_mc_css() gives us the predicted steady-state plasma concentration resulting from an ongoing 1 mg/kg/day exposure.

for (this.id in unique(toxcast.table$DTXSID))
{
# get_cheminfo() gives a list of all the CAS numbers for which HTTK will work:
  if (this.id %in% get_cheminfo(info="dtxsid", suppress.messages=TRUE))
  {
# Set a random number generator seed so that the Monte Carlo will always give
# the same random sequence:
    set.seed(12345)
    toxcast.table[toxcast.table$DTXSID==this.id,"Css"] <- 
      calc_mc_css(dtxsid=this.id,
                  output.units="uM",
                  samples=NSAMP,
                  suppress.messages=TRUE)
    toxcast.table[toxcast.table$DTXSID==this.id,"Css.Type"] <- "in vitro"
  }
}

knitr::kable(toxcast.table[1:10,c("Compound","Q10.AC50","Css","Css.Type")], 
                                caption = "Summarized ToxCast Data",
             floating.environment="sidewaystable")

Step Three: High Throughput Toxicokinetics with QSPR Predictions

For chemicals where no in vitro toxicokinetic data are available, we can sometimes load predictions from in silico models. We do not provide in silico models for \(f_{up}\) and \(Cl_{int}\) themselves (in many cases those models are much larger in terms of file size than the whole httk package). However we do have tables of pre-made predictions form a variety of models. For example, you can load the table of quantitative structure-property relationship (QSPR) model predictions from the Sipes et al. (2017) supplemental material:

load_sipes2017()
#> Loading CLint and Fup predictions from Sipes et al. (2017) for 8719 chemicals.
#> Existing data are not being overwritten.
#> Please wait...

Now go through the list of chemicals again but only calculate if we didn’t already have a value (NOTE: load_sipes2017 will not overwrite the in vitro data by default, so you’ll still get the same answer if you use the same random number seed, but you would also be wasting time):

for (this.id in unique(toxcast.table$DTXSID))
{
  if (this.id %in% get_cheminfo(info="dtxsid", suppress.messages=TRUE) &
    is.na(toxcast.table[toxcast.table$DTXSID==this.id,"Css"]))
  {
# Set a random number generator seed so that the Monte Carlo will always give
# the same random sequence:
    set.seed(12345)
    toxcast.table[toxcast.table$DTXSID==this.id,"Css"] <- 
      calc_mc_css(dtxsid=this.id,
                  output.units="uM",
                  samples=NSAMP,
                  suppress.messages=TRUE)
    toxcast.table[toxcast.table$DTXSID==this.id,"Css.Type"] <- "in silico"
  }
}

Take another look at our table – many of the gaps are now filled in. Any Css values of “NA” (not a number) indicate that we don’t currently have in vitro TK data or in silico predictions from Sipes et al. (2017) for those chemicals (note that there are also predictions available from Pradeep et al. (2020) (load_pradeep2020) and Dawson et al. (2021) (load_dawson2021):

knitr::kable(toxcast.table[1:10,c("Compound","Q10.AC50","Css","Css.Type")], 
                                caption = "Summarized ToxCast Data",
             floating.environment="sidewaystable")

Step Four: Reverse Dosimetry In Vitro-In Vivo Extrapolation

Most IVIVE calculations in httk take advantage of the fact that, to date, the toxicokinetic models included in httk are linear in dose. What this means is that your dose rate is \(X~mg/kg/day\) and the \(C_{ss}\) for \(1~mg/kg/day\) is \(Y~uM\), then the \(C_{ss}\) for \(X~mg/kg/day\) is \(X*Y~uM\), that is: \[C_{ss}\left(X~mg/kg/day\right) = X * C_{ss}\left(1~mg/kg/day\right)\]. Using toxicokinetics to predict tissue concentrations resulting from a known dose is known as “forward dosimetry”. Using toxicokinetics to estimate the dose that might have led to a known tissue cocentration is known as “reverse dosimetry” (Clewell, et al., 2008). httk allows the use of reverse dosimetry to calculate the equivalent steady-state dose to produce a plasma concentration equal to a bioactive concentration identified in vitro (Breen et al., 2021). For example, we can start from the tenth percentile ToxCast \(AC_{50}\) and find the administered equivalent dose that would cause that concentration in the plasma. We do this by dividing the in vitro concentration by the \(C_{ss}\) from 1 mg/kg/day, ending up with a dose rate with units of mg/kd/day that will product the in vitro concentration: \[AED = \frac{AC_{50}}{C_{ss}\left(1~mg/kg/day\right)}~mg/kg/day\] where b0oth \(AC_{50}\) and \(C_{ss}\) must have the same units. Do not forget that the ToxCast concentrations are given on the log-10 scale:

toxcast.table$EquivDose <- signif(10^toxcast.table$Q10.AC50 / toxcast.table$Css,
                                  3)

Look at the table again:

knitr::kable(toxcast.table[1:10,c("Compound","Q10.AC50","Css","EquivDose")], 
                                 caption = "Summarized ToxCast Data",
              floating.environment="sidewaystable")          

Step Five: Compare with Estimated Exposure Rates

The U.S. National Academy of Sciences, Engineering, and Medicine advises calculating chemical risk posed to the public health on the basis of hazard, the dose-response relationship, and exposure (NRC, 1983). We can use in vitro bioactivity as a surrogate for hazard and httk as a surrogate for the dose-reponse relationship. However, to prioritize chemicals on the basis of putative risk we need some estimate of human exposure. Because most chemicals lack measured data characterizing population intake rates (mg/kg body weight/ day) (Egeghy, 2012), a consensus prediction from EPA’s Systematic Empirical Evaluation of Models (SEEM) is often the only option available for exposures. The SEEM consensus prediction is composed of multiple exposure predictors aggregated on the basis of likely exposure pathways (Ring, 2018). The individual exposure predictors are weighted based upon their ability to predict intake rates that could be inferred from biomonitoring data (Ring, 2017; Wambaugh, 2014).
The SEEM exposure forecasts are uncertain, and the upper boundary of the 95% credible interval is used here. Although a 95% credible interval is calculated, this reflects uncertainty but not variability – the (Ring, 2018) model predicts population median exposure rates for the U.S. population.

There are three ways to access the SEEM predictions. If you can access the journal article online, we can load the SEEM predictions from Supplemental Table 1 of the Supporting Information for Ring et al. (2018). Don’t forget that you have to “unzip” a file that ends in “.zip” and place the “.txt” file in your working directory with your other files.

SEEM <- read.csv("SupTable-all.chem.preds-2018-11-28.txt",stringsAsFactors=F)

If you want, you can also download the predictions for specific chemicals from the Comptox Chemicals Dashboard. Be sure to use Batch Search and select the option for “NHANES/Predicted Exposure”. You can download the file in a “comma separated values” or “.CSV” format, in which case you use “read.csv”. Or you can download an Excel “.xls” file which requires load the package readXL.

SEEM <- read_excel("CCD-Batch-Search_DATE.xlsx",sheet=2)

Perhaps most easily we can grab SEEM predictions already in RData format from GitHub Download the file Ring2018Preds.RData and load it:

load("Ring2018Preds.RData")
SEEM <- Ring2018.preds

Again we do not have the space to distribute all the SEEM predictions within this R package, but we can give you our example chemicals:

example.seem <- as.data.frame(subset(SEEM, 
                              dsstox_substance_id %in% toxcast.table$DTXSID))
for (this.col in 4:36) 
  example.seem[,this.col] <- signif(example.seem[,this.col],4)
save(example.seem,file="introtoivive-seem.Rdata",version=2)

Merge the consensus model predictions with your table:

example.seem <- httk::example.seem
toxcast.table <- merge(toxcast.table,
                       example.seem[,c(
                         "dsstox_substance_id",
                         "seem3",
                         "seem3.l95",
                         "seem3.u95",
                         "Pathway",
                         "AD")],
                       by.x="DTXSID",
                       by.y="dsstox_substance_id",
                       all.x = TRUE)

Now we can calculate a bioactivity:exposure ratio (BER), which is a risk-like metric since it takes hazard, dose-response, and exposure into account (Wetmore, 2015):

toxcast.table$BER <- signif(toxcast.table$EquivDose/toxcast.table$seem3.u95, 3)

Reorder your table by BER:

toxcast.table <- toxcast.table[order(toxcast.table$BER),]
knitr::kable(toxcast.table[1:10,c("Compound","EquivDose","seem3.u95","Pathway","BER")], 
                                 caption = "Bioactivity:Exposure Ratios",
              floating.environment="sidewaystable")  

toxcast.table$Compound <- factor(toxcast.table$Compound,
                                    levels=toxcast.table$Compound)

BER.plot <- ggplot(data=toxcast.table) +
  geom_segment(aes(x=Compound,
                   xend=Compound,
                   y=(10^Q10.AC50)/Css,
                   yend=(10^Q90.AC50)/Css),
               size=1,
               color="red",
               alpha=0.5) + 
  geom_segment(aes(x=Compound,
                   xend=Compound,
                   y=seem3.l95,
                   yend=seem3.u95),
               size=1,
               color="blue",
               alpha=0.5) +   
  geom_point(aes(x=Compound,y=(10^Med.AC50)/Css),shape=3,color="red") +
  geom_point(aes(x=Compound,y=seem3),shape=3,color="blue") +
  theme_bw() +
  theme(axis.title.x = element_text(size=8)) + 
  theme(axis.title.y = element_text(size=8)) + 
  scale_y_log10() + 
  theme(axis.text.x = element_text(
    face = "bold", 
    hjust = 1, 
    vjust = 1, 
    size = 4, 
    angle = 45)) +
  ylab("Bioactive Dose & Exposure\nmg/kg BW/day") + 
  xlab("Chemicals Ranked By BER")
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
print(BER.plot)
#> Warning: Removed 16 rows containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Warning: Removed 5 rows containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Warning: Removed 7 rows containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Warning: Removed 16 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Warning: Removed 5 rows containing missing values or values outside the scale range
#> (`geom_point()`).