---
title: "Characterizing Accuracy of Model Predictions for Chemical
Concentration in High Throughput Screening Assays"
author: "Meredith Scherer, Madison Feshuk, James M. Armitage, Jon A. Arnot, 
David Crizer, Michael Devito, Stephen S. Ferguson, Kimberly Freeman, 
Gregory S. Honda, Josh Harrill, Alessandro Sangion, Caroline Tebes-Stevens, 
Katie Paul Friedman, John F. Wambaugh"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{j) Scherer (Submitted): In Vitro Distribution}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
*Please send questions to wambaugh.john@epa.gov*

from "Characterizing Accuracy of Model Predictions for Chemical 
Concentration in High Throughput Screening Assays"

Meredith Scherer, Madison Feshuk, James M. Armitage, Jon A. Arnot, 
David Crizer, Michael Devito, Stephen S. Ferguson, Kimberly Freeman, 
Gregory S. Honda, Josh Harrill, Alessandro Sangion, Caroline Tebes-Stevens, 
Katie Paul Friedman, John F. Wambaugh

## Abstract

The aim of in vitro - in vivo extrapolation (IVIVE) of dose is to predict 
potential whole animal effects from perturbations observed in vitro. Reliance 
on nominal (administered) concentration (assumed dose) does not account for 
chemical distribution to cellular and non-cellular elements of an in vitro 
assay system. Chemical distribution can reduce the free concentration 
available to cause effects and can result in an inaccurate estimate of the 
intracellular concentration causing any observed perturbations. There are 
mathematical, chemical property-based partitioning models for predicting 
cellular concentrations when not measured experimentally. This work evaluated 
two of these in vitro disposition models, Armitage et al., 2021 and Kramer et 
al., 2010, using a total of 153 experimental intracellular concentration 
measurements of 43 chemicals, along with parameters describing measurement 
conditions, from 12 peer reviewed studies in addition to data generated for 
this study. Intracellular concentrations were more accurately predicted by 
both the Armitage model (root mean squared log10 error (RMSLE) = 1.12) and the 
Kramer model (RMSLE = 1.30) than by the nominal concentration (RMSLE = 1.45). 
Although limited by the amount of available measurement data, these results 
indicate that mathematical modeling of in vitro distribution can improve the 
accuracy of IVIVE for toxicology.

## HTTK Version 

This vignette was created with **httk** v2.X.0. Although we attempt to maintain 
backward compatibility, if you encounter issues with the latest release of 
**httk**
and cannot easily address the changes, historical versions of httk are 
available from: https://cran.r-project.org/src/contrib/Archive/httk/

## Prepare for session
R package **knitr** generates html and PDF documents from this RMarkdown file,
Each bit of code that follows is known as a "chunk". We start by telling 
**knitr** how we want our chunks to look.
```{r configure_knitr, eval = TRUE}
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
options(rmarkdown.html_vignette.check_title = FALSE)
```

### Clear the memory
It is a bad idea to let variables and other information from previous R
sessions float around, so we first remove everything in the R memory.
```{r clear_memory, eval = TRUE}
rm(list=ls()) 
```

### eval = execute.vignette
If you are using the RMarkdown version of this vignette (extension, .RMD) you
will be able to see that several chunks of code in this vignette
have the statement "eval = execute.vignette". The next chunk of code, by 
default, sets execute.vignette = FALSE.  This means that the code is included 
(and necessary) but was not run when the vignette was built. We do this because 
some steps require extensive computing time and the checks on CRAN limit how 
long we can spend building the package. If you want this vignette to work, 
you must run all code, either by cutting and pasting it into R. Or, if viewing 
the .RMD file, you can either change execute.vignette to TRUE or press "play" 
(the green arrow) on each chunk in RStudio.
```{r runchunks, eval = TRUE}
# Set whether or not the following chunks will be executed (run):
execute.vignette <- FALSE
```

### Load the relevant libraries
We use the command 'library()' to load various R packages for our analysis.
If you get the message "Error in library(X) : there is no package called 'X'"
then you will need to install that package: 

From the R command prompt:

install.packages("X")

Or, if using RStudio, look for 'Install Packages' under 'Tools' tab.
```{r load_packages, eval = execute.vignette}
#
#
# Setup
#
#
library(httk)        # run Armitage.R and Kramer.R in vitro distribution models
library(tidyverse)   # data wrangling
library(data.table)  # data table functionality
library(ggplot2)     # plotting
library(ggrepel)     # jitter labels in plots
library(ggpubr)      # arrange plots for publication
library(viridis)     # colorblind palette
library(scales)      # required for label_log

#save figures and tables to working directory: yes/no
save_output = FALSE

#for saving outputs
data.date <- "16JUL2025"
path_out <- "C:/Users/mscherer/OneDrive - Environmental Protection Agency (EPA)
/Profile/Desktop/IVDManuscriptData/manuscript figures 16JUL25/"

```

```{r rmsle_func, eval = execute.vignette}
#function to calculate RMSLE
sle <- function (actual, predicted) 
{
    log.act <- log10(actual)
    log.pred <- log10(predicted)
    good <- !is.na(log.act) & is.finite(log.act) &
            !is.na(log.pred) & is.finite(log.pred)
    log.act <- log.act[good]
    log.pred <- log.pred[good]
    return((log.act - log.pred)^2)
}
msle <- function (actual, predicted) 
{
    mean(sle(actual, predicted))
}
rmsle <-function (actual, predicted) 
{
    sqrt(msle(actual, predicted))
}
```

# Data Collection

Data sets were curated from the literature and combined with two in house 
experiments to evaluate the in vitro distribution models. 

# Load and format the data
```{r load_IVD_Manuscript_data, eval = execute.vignette}

IVD_Manuscript_data <- httk::Scherer2025.IVD

cat(paste("Summarized data from ", 
          length(unique(IVD_Manuscript_data$Citation)),
          " studies covering ", 
          length(unique(IVD_Manuscript_data$ChemicalName)),
          " unique chemicals ", 
          "with ",
          length(IVD_Manuscript_data$CellConcentration_uM),
          " observations total.",
          sep="")
    )

#14 studies, 43 chemicals, 153 observations

#number each row
IVD_Manuscript_data[,X:=1:length(IVD_Manuscript_data$ChemicalName)]

```

#In vitro assay descriptors (optional step)
In vitro assay descriptors for the measurements in Supplemental Materials 
Table S2 were added as additional assays into 
httk:invitro.assay.params using the AENM “IVD” + Author + Year 
(that is, Smith et al. (2000) would be AENM “IVDSmith2000”.
```{r invitrodescriptors, eval = execute.vignette}

### Create assay_component_endpoint_name columns for assay parameters ###

#create assay_name column
IVD_Manuscript_data[, assay_name := 
                      gsub(", ", "", IVD_Manuscript_data$Citation)]

#create the assay_component_endpoint_name column
IVD_Manuscript_data[, assay_component_endpoint_name := 
                      paste("IVD", gsub(", ", "", assay_name), sep="")]

#tack on the "interesting" part for each assay to create the endpoint_name 
IVD_Manuscript_data[assay_component_endpoint_name=="IVDBellwon2015b" |
                      assay_component_endpoint_name == "IVDBroeders2013" |
                      assay_component_endpoint_name == "IVDWei2010",
                    assay_component_endpoint_name := 
                      paste(assay_component_endpoint_name, CellType, sep="_")]
#these two have multiple concentrations or time points
IVD_Manuscript_data[assay_component_endpoint_name=="IVDMundy2004" |
                      assay_component_endpoint_name=="IVDLundquist2014",
                assay_component_endpoint_name := 
                  paste(paste(assay_component_endpoint_name, 
                              paste(ExposureDuration_hrs, "hr", sep = ""), 
                              sep="_"), paste(FBSf, "FBSf", sep = ""), 
                        sep = "_")]

#create a copy with the assay parameters for the options table (Table S7)
options.IVD_Manuscript_data<-data.table::copy(IVD_Manuscript_data)

### Pull in httk assay parameter table ###
invitro.assay.params <- data.table::copy(invitro.assay.params)

#copy the experimental data frame to pull out the parameter lines
parameters_only<-data.table::copy(IVD_Manuscript_data)

### Set up assay parameters to match with invitro.assay.params.RData ###

#rename columns to align with the names in invitro.assay.params.RData
parameters_only[, timepoint_hr := ExposureDuration_hrs] %>%
  .[, media_serum := SerumType] %>%
  .[, cell_type := CellType] 

# get the list of columns we need
ivpnames<-names(invitro.assay.params)

#add on the extra columns to match invitro.assay.params and parameters_only
if(!all(ivpnames%in%names(parameters_only))){
    parameters_only[,ivpnames[!(ivpnames %in% names(parameters_only))]] <- 
      as.double(NA)}

#filter parameters_only to match up with ivpnames and only the unique rows
parameters_only<-parameters_only %>% select(as.character(ivpnames)) %>% unique()

# reassign own parameters table to pass to armitage and kramer functions
user_assay_parameters<-parameters_only
# should have 29 + 1649 = 1678 rows

#remove parameter columns since we will be pulling from "user_assay_parameters"
ivpnames[2]<-NA #remove endpoint_name bc we need to be able to reference it
IVD_Manuscript_data[,ivpnames[(ivpnames %in% names(IVD_Manuscript_data))]]<- 
  NULL

```


# Run the in vitro distribution models
```{r run_IVD_models, eval = execute.vignette}

# run Armitage model, output concentrations in umol/L (e.g. uM)
armitageOutput_data.dt <- armitage_eval(tcdata=IVD_Manuscript_data,
                                       surface.area.switch = TRUE,
                                       user_assay_parameters = 
                                         user_assay_parameters,
                                       restrict.ion.partitioning = TRUE)



# run Kramer model, output concentrations in umol/L (e.g. uM)
kramerOutput_data.dt <- kramer_eval(tcdata=IVD_Manuscript_data, 
                                    surface.area.switch = TRUE,
                                    user_assay_parameters = 
                                      user_assay_parameters,
                                    restrict.ion.partitioning = TRUE)

```

# Manipulate data and prep for plotting
```{r manipulate_for_plotting, eval = execute.vignette}

#delete nomconc column
Armitage.dt<-subset(armitageOutput_data.dt, select=-c(nomconc))

#pull ioc types from armitage since kramer doesn't generate them
IOC_Type_output<-Armitage.dt[,.(ChemicalName, IOC_Type)] %>% distinct()

#delete nomconc column and add in the IOC type column for Kramer
Kramer.dt_subset<-subset(kramerOutput_data.dt, select=-c(nomconc))
Kramer.dt<-Kramer.dt_subset[IOC_Type_output, on = "ChemicalName"]

#delete cellconc column and add in the IOC type column for Nominal
Nominal.dt_subset<-subset(kramerOutput_data.dt, 
                          select = -c(concentration_cells))
Nominal.dt<-Nominal.dt_subset[IOC_Type_output, on = "ChemicalName"]

#add tags
Armitage.dt[,label:= "Armitage"]
Kramer.dt[,label:= "Kramer"]
Nominal.dt[,label:= "Nominal"]

#combine data for plotting (datasets are side by side)
combodata<-merge(Armitage.dt, kramerOutput_data.dt, by = "X", all.x = TRUE) 
  #kramer output still has nomconc values 

if(save_output){
  write.csv(combodata, file = paste0(path_out,"ManuscriptOutput_",
                                     data.date,"_",Sys.Date(),".csv"))
}

```

```{r ggplot_set_up, eval = execute.vignette}

#set default theme elements for visually consistent manuscript plots
personal_theme <- function(){ 
  theme_bw() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(colour = "black", linewidth=1),
        axis.title = element_text(size = rel(1.25)),
        axis.text = element_text(size = rel(1.15)),
        plot.title = element_text(size = rel(1.5),
                                  hjust = 0.5),
        panel.grid = element_blank(),
        legend.position = "none") 
}

#set plot colors to correspond with each chemica
plotcolors<-scales::viridis_pal()(length(
  unique(armitageOutput_data.dt$ChemicalName)))
names(plotcolors)<-armitageOutput_data.dt$ChemicalName %>% unique() %>% sort()

```

# Main Paper Plots:

## Figure 1: 
In Vitro Distribution Diagram (no code)

## Figure 2: 
Predicted and experimental intracellular concentrations from the 
Cnominal (left panel), Armitage (middle panel), and Kramer (right panel) models.
Dashed line shows unity.
```{r joint_intracel_plot_fig2, eval = execute.vignette}

#ARMITAGE:

#Armitage RMSLE:
arm_rmsle_long<-rmsle((Armitage.dt$CellConcentration_uM),(Armitage.dt$ccells))
arm_rmsle<-round(arm_rmsle_long, 2)
#1.12

#Armitage r^2:
ml_arm_overall = lm(CellConcentration_uM~ccells, data = Armitage.dt)
summary(ml_arm_overall)$r.squared 
#0.33

#plot set up
yrng <- range(Armitage.dt$CellConcentration_uM)
xrng <- range(Armitage.dt$ccells)
n_chems<-n_distinct(Nominal.dt$ChemicalName)
n_data<-nrow(Nominal.dt)

# now plot
Armitage_ccell<-ggplot(Armitage.dt) +
  geom_point(aes(ccells, CellConcentration_uM, color = ChemicalName)) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  xlab(expression("Armitage Predicted "*italic("C"["cell"])~"(\u03BCM)")) +
  ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) +
  scale_x_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  scale_y_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste(
      "N[data] ", "~'='~", n_data)), hjust = 1, vjust = 0, size = 4) +
  annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste(
      "N[chemicals] ", "~'='~", n_chems)), hjust = 1, vjust = 1.2, size = 4)+
  annotate(geom = "text", x = xrng[2], y = yrng[1], 
           label = paste("RMSLE =", arm_rmsle), 
           hjust = 1, vjust = 3, size = 4) +
  personal_theme() +
  scale_color_manual(values=plotcolors)

#KRAMER:

#Kramer RMSLE:
kram_rmsle_long<-
  rmsle((Kramer.dt$CellConcentration_uM),(Kramer.dt$concentration_cells))
kram_rmsle<-format(round(kram_rmsle_long, digits=2), nsmall = 2) 
#1.30

#Kramer r^2:
ml_kram_overall = lm(CellConcentration_uM~concentration_cells, data = Kramer.dt)
summary(ml_kram_overall)$r.squared 
#0.03

#plot
Kramer_ccell<-ggplot(Kramer.dt)+
  geom_point(aes(concentration_cells, CellConcentration_uM, 
                 colour = ChemicalName)) + 
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  scale_x_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  scale_y_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  xlab(expression("Kramer Predicted "*italic("C"["cell"])~"(\u03BCM)")) +
  ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) +
  annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste(
      "N[data] ", "~'='~", n_data)), hjust = 1, vjust = 0, size = 4) +
  annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste(
      "N[chemicals] ", "~'='~", n_chems)), hjust = 1, vjust = 1.2, size = 4)+
  annotate(geom = "text", x = xrng[2], y = yrng[1], 
  label = paste("RMSLE =", kram_rmsle), hjust = 1, vjust = 3, size = 4) +
  personal_theme() +
  scale_color_manual(values=plotcolors)

#NOMINAL:

#Nominal RMSLE:
nom_rmsle_long<-rmsle((Nominal.dt$CellConcentration_uM),(Nominal.dt$nomconc))
nom_rmsle<-round(nom_rmsle_long, 2)
#1.45

Nominal_ccell<-ggplot(Nominal.dt) +
  geom_point(aes(nomconc, CellConcentration_uM, colour = ChemicalName))+
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  scale_x_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  scale_y_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  xlab(expression(italic("C"["nominal"])~"(\u03BCM)")) +
  ylab(expression("Experimental "~italic("C"["cell"])~"(\u03BCM)")) +
  labs(colour = "Chemical")+
  annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste(
      "N[data] ", "~'='~", n_data)), hjust = 1, vjust = 0, size = 4) +
  annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste(
      "N[chemicals] ", "~'='~", n_chems)), hjust = 1, vjust = 1.2, size = 4) +
  annotate(geom = "text", x = xrng[2], y = yrng[1], 
           label = paste("RMSLE =", nom_rmsle), 
           hjust = 1, vjust = 3, size = 4) +
  personal_theme() +
  theme(legend.position = "bottom",
        legend.key.size = unit(2, 'mm'),
        legend.spacing = unit(10, "pt"))+
        # legend.text = element_text(size = 8))+
  guides(col = guide_legend(ncol = 5)) +
  scale_color_manual(values=plotcolors)

#Combine the three plots w RMSLE
Ccells_comboplot<-ggarrange(Nominal_ccell, Armitage_ccell, Kramer_ccell, 
                            ncol=3 , legend = "bottom", common.legend = TRUE)

#save
if(save_output){
  ggsave(file = paste0("Figure2_",data.date,"_",Sys.Date(),".png"),
         path = path_out,
         Ccells_comboplot, width = 10, height = 5.5, dpi = 300)
}

print(Ccells_comboplot)

```

## Figure 3: 
Predicted and experimental intracellular concentrations from the Armitage 
(left panel) and Kramer (right panel) models for the ten chemicals with 
multiple observations. Solid lines connect multiple measurements for a 
single chemical. The dashed line shows unity. 
```{r multiobs_intracel_plot_fig3, eval = execute.vignette}

#Armitage.multobs.dt %>% count(ChemicalName, NominalDose_uM, Citation)
#some of these have super close nomconcs because they are experimentally 
#measured - not really relevant for this comparison because they are replicates
#chemicals to remove:
#Hexachlorobenzene, Malathion, Pentachlorophenol, Propiconazole 
#(all of the chemicals from Stadnicka, 2014)

multiObsChemnameList<-Armitage.dt %>% 
  dplyr::filter(Citation != "Stadnicka, 2014") %>% 
  count(ChemicalName, NominalDose_uM) %>% 
  count(ChemicalName) %>% dplyr::filter(n>=2) %>% 
  select(ChemicalName) %>% distinct()
#10 chemicals with multiple nominal doses measured

#filter measured data to just these chemicals
Armitage.multobs.dt_og <- Armitage.dt %>% 
  dplyr::filter(ChemicalName %in% multiObsChemnameList$ChemicalName)

Kramer.multobs.dt_og <- Kramer.dt %>% 
  dplyr::filter(ChemicalName %in% multiObsChemnameList$ChemicalName)

Nominal.multobs.dt_og <- Nominal.dt %>% 
  dplyr::filter(ChemicalName %in% multiObsChemnameList$ChemicalName)
#each have 89 observations

#merge the lines with multiple observations for the same nomconc (ie tox21)
Armitage.multobs.dt<-Armitage.multobs.dt_og %>% 
  group_by(ChemicalName, NominalDose_uM) %>% 
  mutate(mean_measuredccell = mean(CellConcentration_uM),
         mean_predictedccell = mean(ccells)) %>% 
  select(ChemicalName, NominalDose_uM, mean_measuredccell, 
         mean_predictedccell) %>% distinct()

Kramer.multobs.dt<-Kramer.multobs.dt_og %>% 
  group_by(ChemicalName, NominalDose_uM) %>% 
  mutate(mean_measuredccell = mean(CellConcentration_uM),
         mean_predictedccell = mean(concentration_cells)) %>% 
  select(ChemicalName, NominalDose_uM, mean_measuredccell, 
         mean_predictedccell) %>% distinct()

Nominal.multobs.dt<-Nominal.multobs.dt_og %>% 
  group_by(ChemicalName, NominalDose_uM) %>% 
  mutate(mean_measuredccell = mean(CellConcentration_uM),
         mean_predictedccell = mean(NominalDose_uM)) %>% 
  select(ChemicalName, NominalDose_uM, mean_measuredccell, 
         mean_predictedccell) %>% distinct()
#down to 50 observations because the duplicates have been averaged

### Plotting ###

#plot set up
yrng_mult <- range(Armitage.multobs.dt$mean_measuredccell)
xrng_mult <- range(Armitage.multobs.dt$mean_predictedccell)
n_chems_mult<-n_distinct(Armitage.multobs.dt$ChemicalName)
n_data_mult<-nrow(Armitage.multobs.dt)

#ARMITAGE

#Armitage RMSLE:
arm_mult_rmsle_long<-rmsle((Armitage.multobs.dt$mean_measuredccell),
                           (Armitage.multobs.dt$mean_predictedccell))
arm_mult_rmsle<-round(arm_mult_rmsle_long, 2)
#1.27

#r squared 
ml_arm <- lm(mean_measuredccell~mean_predictedccell, data = Armitage.multobs.dt)
rsq_armitage_mult <- summary(ml_arm)$r.squared
#0.87

# now plot
Armitage_ccell_mult<-ggplot(Armitage.multobs.dt) +
  geom_point(aes(mean_predictedccell, mean_measuredccell, 
                 colour = ChemicalName)) +
  geom_smooth(aes(mean_predictedccell, mean_measuredccell, 
                  colour = ChemicalName), method=lm, se = FALSE) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  scale_x_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  scale_y_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  xlab(expression("Armitage Predicted "*italic("C"["cell"])~"(\u03BCM)")) +
  ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) +
  labs(colour = "Chemical")+
  annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], 
           label = parse(text = paste("N[data] ", "~'='~", n_data_mult)), 
           hjust = 1, vjust = 0, size = 4) +
  annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], 
           label = parse(text = paste("N[chemicals] ", "~'='~", n_chems_mult)), 
           hjust = 1, vjust = 1.2, size = 4)+
  annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], 
  label = paste("RMSLE =", arm_mult_rmsle), hjust = 1, vjust = 3, size = 4) +
  personal_theme()  +
  scale_color_manual(values=plotcolors)

#KRAMER

#Kramer RMSLE:
kram_mult_rmsle_long<-rmsle((Kramer.multobs.dt$mean_measuredccell),
                            (Kramer.multobs.dt$mean_predictedccell))
kram_mult_rmsle<-round(kram_mult_rmsle_long, 2)
# 1.58

#r squared 
ml_kram <- lm(mean_measuredccell~mean_predictedccell, data = Kramer.multobs.dt)
rsq_kramer_mult <- summary(ml_kram)$r.squared
#0.87

# now plot
Kramer_ccell_mult<-ggplot(Kramer.multobs.dt) +
  geom_point(aes(mean_predictedccell, mean_measuredccell, 
                 colour = ChemicalName)) +
  geom_smooth(aes(mean_predictedccell, mean_measuredccell, 
                  colour = ChemicalName), method=lm, se = FALSE) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  scale_x_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  scale_y_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), 
                labels = label_log(digits = 2)) +
  xlab(expression("Kramer Predicted "*italic("C"["cell"])~"(\u03BCM)")) +
  ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) +
  labs(colour = "Chemical")+
  annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], 
           label = parse(text = paste("N[data] ", "~'='~", n_data_mult)), 
           hjust = 1, vjust = 0, size = 4) +
  annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], 
           label = parse(text = paste("N[chemicals] ", "~'='~", n_chems_mult)), 
           hjust = 1, vjust = 1.2, size = 4) +
  annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], 
  label = paste("RMSLE =", kram_mult_rmsle), hjust = 1, vjust = 3, size = 4) +
  personal_theme()  +
  scale_color_manual(values=plotcolors) +
  theme(legend.key.size = unit(2, 'mm'),
        legend.text = element_text(size = 6))


#arrange
Ccells_comboplot_mult<-ggarrange(Armitage_ccell_mult, Kramer_ccell_mult, ncol=2,
                                 legend = "bottom", common.legend = TRUE)

if(save_output){
  ggsave(file = paste0("Figure3_",data.date,"_",Sys.Date(),".png"),
         path = path_out,
         Ccells_comboplot_mult, width = 8, height = 4, dpi = 300)
}

print(Ccells_comboplot_mult)

```

## Figure 4: 
Chemical partitioning into medium, air, cell, and plastic compartments for 
the Armitage and Kramer models. Dashed line shows unity.
```{r chemicalPart_intracel_plot_fig4, eval = execute.vignette}

#Cwat_uM -------- Cwat_R
Cwater_comparison<-ggplot(combodata) +
  geom_point(aes(concentration_medium, cwat_s, colour = ChemicalName.x)) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  ggtitle("Medium Concentration") +
  xlab(expression("Kramer Predicted "*italic("C"["medium"])~"(\u03BCM)")) +
  ylab(expression("Armitage Predicted "*italic(" C"["medium"])~"(\u03BCM)")) +
  scale_x_log10(labels = label_log(digits = 2)) +
  scale_y_log10(labels = label_log(digits = 2)) +
  personal_theme() +
  scale_color_manual(values=plotcolors)

#Cair_uM -------- Cair_R
Cair_comparison<-ggplot(combodata) +
  geom_point(aes(concentration_air, cair, colour = ChemicalName.x)) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  ggtitle("Air Concentration") +
  xlab(expression("Kramer Predicted "*italic("C"["air"])~"(\u03BCM)")) +
  ylab(expression("Armitage Predicted "*italic(" C"["air"])~"(\u03BCM)")) +
  scale_x_log10(labels = label_log(digits = 2)) +
  scale_y_log10(labels = label_log(digits = 2)) +
  personal_theme() +
  scale_color_manual(values=plotcolors)

#C_cells_uM ----- Ccells
Ccells_comparison<-ggplot(combodata) +
  geom_point(aes(concentration_cells, ccells, colour = ChemicalName.x)) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  ggtitle("Cell Concentration") +
  xlab(expression("Kramer Predicted "*italic("C"["cell"])~"(\u03BCM)")) +
  ylab(expression("Armitage Predicted "*italic(" C"["cell"])~"(\u03BCM)")) +
  scale_x_log10(lim = c(10^-4, 10^5), labels = label_log(digits = 2), 
                breaks = c(10^-2, 10^0, 10^2, 10^4)) +
  scale_y_log10(lim = c(10^-4, 10^5), labels = label_log(digits = 2), 
                breaks = c(10^-2, 10^0, 10^2, 10^4)) +
  personal_theme() +
  scale_color_manual(values=plotcolors)

#Aplastic_uM_m2 - (Cplastic_R / Sarea_R)
Aplastic_comparison<-ggplot(combodata) +
  geom_point(aes(concentration_plastic, cplastic, colour = ChemicalName.x)) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  ggtitle("Plastic Concentration") +
  xlab(expression("Kramer Predicted "*italic("C"["plastic"])~"(\u03BCM)")) +
  ylab(expression("Armitage Predicted "*italic(" C"["plastic"])~"(\u03BCM)")) +
  labs(colour = "Chemical")+
  scale_x_log10(labels = label_log(digits = 2), 
                breaks = c(10^-6, 10^-4, 10^-2, 10^0)) +
  scale_y_log10(labels = label_log(digits = 2)) +
  personal_theme() +
  scale_color_manual(values=plotcolors) +
  theme(legend.key.size = unit(1, 'mm'),
        legend.spacing = unit(2, "pt"),
        legend.text = element_text(size = 6))


compartment_plots<-ggarrange(Ccells_comparison, Aplastic_comparison, 
                             Cwater_comparison, Cair_comparison, 
                             common.legend = TRUE, legend = "none")

if(save_output){
  ggsave(file = paste0("Figure4_",data.date,"_",Sys.Date(),".png"),
         path = path_out,
         compartment_plots, width = 8, height = 6, dpi= 300)
}

print(compartment_plots)

```

## Figure 5:
Change in per-chemical RMSLE for predicted and experimental intracellular 
concentrations from the Armitage model with predictions made using batch 
downloaded physico-chemical property predictions from OPERA and ChemAxon or 
individually inspected physico-chemical property predictions. Boxed numbers 
show difference in RMSLE pertaining to corresponding chemicals 
and the dashed line shows unity. 

```{r curated_values_plot_fig5, eval = execute.vignette}

#load curated data info
original_curated.data.dt <- httk::Dimitrijevic.IVD

#first - run using the default/dashboard values
default.dt<-data.table::copy(original_curated.data.dt %>% 
                   select(-c( "Arnot_pka", "Arnot_pkb", "log KOW,N", "log KAW,N", 
                              "Predicted_Ccell_µM", "FoA")))

# run the model, output concentrations in umol/L (e.g. uM)
armitageOutput_default.dt <- armitage_eval(tcdata=default.dt,
                                        restrict.ion.partitioning = TRUE,
                                        surface.area.switch = FALSE)

#next- run using curated values (EAS-E Suite)
curated.data.dt<-data.table::copy(original_curated.data.dt)

#overwrite using curated physchem properties
curated.data.dt[, pKa_Donor := as.character(Arnot_pka)] %>% #acidic
  .[, pKa_Accept := as.character(Arnot_pkb)] %>% #basic
  .[, gkow := `log KOW,N`] %>% 
  .[, gkaw_n := `log KAW,N`] 

# run the model, output concentrations in umol/L (e.g. uM)
armitageOutput_curated.data.dt <- armitage_eval(tcdata=curated.data.dt,
                                        restrict.ion.partitioning = TRUE,
                                        surface.area.switch = FALSE)

#save output with these values (Supplemental Materials T9)
defaultvals<- armitageOutput_default.dt %>% 
  select(Name, casrn, gkow_n, gkaw_n, pKa_Donor, pKa_Accept)
curatedvals<-armitageOutput_curated.data.dt %>% 
  select(Name, casrn, gkow_n, gkaw_n, pKa_Donor, pKa_Accept)
#tie the two tables together for the supplement
s9_table<-rbind(curatedvals, defaultvals)

#calculate rmsle

#default values

#individual
armitageOutput_default.dt[,rmsle_ccell_default:=
                            rmsle(Reported_Ccell_µM, ccells), by = Name]

#total
RMSLE_default.dt<-rmsle(armitageOutput_default.dt$Reported_Ccell_µM, 
                        armitageOutput_default.dt$ccells)
RMSLE_default.dt<-round(RMSLE_default.dt, 3)
#0.586

#curated values

#individual
armitageOutput_curated.data.dt[,rmsle_ccell_curated:=
                                 rmsle(Reported_Ccell_µM, ccells), by = Name]

#total
RMSLE_curated.data.dt<-rmsle(armitageOutput_curated.data.dt$Reported_Ccell_µM, 
                             armitageOutput_curated.data.dt$ccells)
RMSLE_curated.data.dt<-round(RMSLE_curated.data.dt, 3)
#0.570

#improvement when using curated data:
RMSLE_default.dt-RMSLE_curated.data.dt
#0.016

## difference between the two on a per-chemical basis ##
curated_data.dt<-merge(armitageOutput_default.dt, 
                       armitageOutput_curated.data.dt, by = "Name")

#calculate rmsle difference for each chemical
curated_data.dt[,rmsle_difference:=
                  rmsle_ccell_curated-rmsle_ccell_default, by = Name]

#which chemicals had the greatest magnitude of difference
#largest improvement when using curated data
head(curated_data.dt %>% 
  select(rmsle_difference, Name) %>% 
  arrange(rmsle_difference), 2)

#largest improvement when using default data
tail(curated_data.dt %>% 
  select(rmsle_difference, Name) %>% 
  arrange(rmsle_difference), 2)

#set custom colors and shapes for plot below
#tricky to tell apart with just colors
custom_shapes <- c(22, 23, 24, 25, 22, 23, 
                   24, 25, 22, 23, 24, 25) # 12 distinct shapes
custom_colors <- scales::viridis_pal()(length(
  unique(curated_data.dt$Name))) # 12 distinct colors

#plot
curated_difference_httk_plot<-ggplot(curated_data.dt) +
  geom_point(aes(ccells.y, Reported_Ccell_µM.y, fill = Name, shape = Name), size = 4) +
  geom_label_repel(aes(ccells.y, Reported_Ccell_µM.y, 
                       label = round(rmsle_difference, 2))) +
  geom_abline(intercept=0,slope=1, linetype = "dashed") +
  xlab(expression("Armitage Predicted "*italic("C"["cell"])~"(\u03BCM)")) +
  ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) +
  labs(fill = "Chemical", shape = "Chemical")+
  scale_x_log10(lim = c(10, 10000), labels = label_log(digits = 2)) +
  scale_y_log10(lim = c(10, 10000), labels = label_log(digits = 2)) +
  personal_theme() +
  scale_shape_manual(values = custom_shapes) +
  scale_fill_manual(values = custom_colors) +
  theme(legend.position = "bottom",
        legend.key.size = unit(2, 'mm'),
        legend.spacing = unit(10, "pt"),
        legend.text = element_text(size = 10))

if(save_output){
  ggsave(file = paste0("Figure5_",data.date,"_",Sys.Date(),".png"),
         path = path_out,
         curated_difference_httk_plot, width = 8, height = 5, dpi = 300)
  
  write.csv(s9_table, file = paste0(path_out,"table_s9_",data.date,"_",
                                    Sys.Date(),".csv"))
}

```

# Main Paper Calculations:

## Over/underestimation
```{r fold_change_calculation, eval = execute.vignette}

# The Armitage model resulted in a slight over prediction of Ccell whereas 
# the Kramer model had the greatest spread but most evenly over and 
# underpredicted Ccell (Figure 2). 
range(Armitage.dt$ccells)
#1.555754e-02 4.845393e+04
range(Kramer.dt$concentration_cells)
#5.838149e-04 2.663947e+04

#Kramer has a larger spread of predictions

##create datatable with both outputs stacked 
stacked.data<-rbind(Armitage.dt, Kramer.dt, fill=TRUE)

#overpredict or underpredict for each model compared to experimental
stacked.data_HL<-data.table::copy(stacked.data)
stacked.data_HL[ccells>CellConcentration_uM, 
                ARMITAGE_ccellcomp:="ARMITAGE higher"] %>%
  .[ccells<CellConcentration_uM, 
    ARMITAGE_ccellcomp:="ARMITAGE lower"] %>%
  .[concentration_cells>CellConcentration_uM, 
    KRAMER_ccellcomp:="KRAMER higher"] %>% 
  .[concentration_cells<CellConcentration_uM, 
    KRAMER_ccellcomp:="KRAMER lower"] 

stacked.data_HL %>% count(ARMITAGE_ccellcomp, KRAMER_ccellcomp)
#Armitage model more often over predicts
#Kramer is pretty even split of over and under predictions

#pick out where armitage prediction or kramer prediction is larger
sidebyside.data_HL<-data.table::copy(combodata)
sidebyside.data_HL[ccells>concentration_cells, ARMITAGE:="ARMITAGE higher"] %>%
  .[ccells<concentration_cells, KRAMER:="KRAMER higher"]

#count which is higher or lower by chemical
sidebyside.data_HL %>% select(ChemicalName.x, ARMITAGE, KRAMER) %>% 
  distinct() %>% count(ARMITAGE, KRAMER)
#armitage: higher for 38 chemicals
#kramer: higher for 5 chemicals
#38+5 = 43 chemicals overall

```

## Effect of water solubility value for Cyclosporin A
```{r water_solublity, eval = execute.vignette}

#Calculating water solubility using OPERA gswat
#pull cyc a values from input data
cycA.dt<-data.table::copy(IVD_Manuscript_data %>% 
                dplyr::filter(ChemicalName == "Cyclosporin A"))

#using regular gswat value:
#run armitage 
armitageOutput_cycA.dt_orig <- armitage_eval(tcdata=cycA.dt,
                                       surface.area.switch = TRUE,
                                       user_assay_parameters =
                                         user_assay_parameters,
                                       restrict.ion.partitioning = TRUE)

#regular gswat value is 2.341266 mg/L
#armitageOutput_cycA.dt_orig$gswat_n

#check gswat value: 
# 1. pull from data.table
chem.physical_and_invitro.data %>% 
  dplyr::filter(CAS == "59865-13-3") %>% 
  select(MW, logWSol, logWSol.Reference)

# 2. convert mol/L to mg/L which is what armitage uses
#from chem.physical_and_invitro.data table:
-3.739 #(mol/L) log10 water solubility
10^-3.739 # un-log = 0.0001823896 (mol/L)
#convert to mg/L
0.0001823896*(1203*convert_units("g", "mg")) #1203 is MW in g/mol
#output: 219.4147 mg/L = swat
log10(219.4147) 
#output: 2.341266 mg/L = gswat

#calculate RMSLE
cycA.dt_rmsle_orig<-round(rmsle((
  armitageOutput_cycA.dt_orig$CellConcentration_uM),
  (armitageOutput_cycA.dt_orig$ccells)),2)
#2.89 - overall rmsle

#Calculating water solubility using experimental gswat
#pull cyc a values from input data
cycA.dt<-data.table::copy(IVD_Manuscript_data %>% dplyr::filter(ChemicalName == "Cyclosporin A"))

#using measured gswat value: mg/L
cycA.dt[, gswat_n := log10(0.0073)]

#run armitage 
armitageOutput_cycA.dt_meas <- armitage_eval(tcdata=cycA.dt,
                                       surface.area.switch = TRUE,
                                       user_assay_parameters =
                                         user_assay_parameters,
                                       restrict.ion.partitioning = TRUE)

#measured gswat value is -2.136677 mg/L
#armitageOutput_cycA.dt_meas$gswat_n

#calculate RMSLE
cycA.dt_rmsle_meas<-round(rmsle((
  armitageOutput_cycA.dt_meas$CellConcentration_uM),
  (armitageOutput_cycA.dt_meas$ccells)),2)
#1.47 - overall rmsle

```

## Effect of experimentally measured nominal dose
Change in accuracy of intracellular concentration predicted by studies 
reporting experimentally measured nominal dose
```{r experimental_nomconc, eval = execute.vignette}

#pull observations that did experimentally measure the nominal dose
#exp_nom.dt<-data.table::copy(IVD_Manuscript_data[!(is.na(ExperimentalDose_uM)),]) 
#6 citations did measure experimental dose
#42 observations

#BUT - Stadnicka 2014 did not report a nominal dose to begin with, 
#only the experimentally measured values
exp_nom.dt<-data.table::copy(IVD_Manuscript_data[!(is.na(ExperimentalDose_uM)) &
                                       Citation != "Stadnicka, 2014",]) 
#unique(exp_nom.dt$Citation)
#5 citations measured experimental dose AND reported nominal dose
#27 observations

#run armitage with nomconc as NominalDose_uM.x (not experimentally measured)
arm_nomconc<-data.table::copy(exp_nom.dt[, nomconc := NominalDose_uM])

# run Armitage model, output concentrations in umol/L (e.g. uM)
armitageOutput_arm_nomconc <- armitage_eval(tcdata=arm_nomconc,
                                            surface.area.switch = TRUE,
                                            user_assay_parameters = 
                                              user_assay_parameters,
                                            restrict.ion.partitioning = TRUE)

#run armitage with nomconc as ExperimentalDose_uM.x (experimentally measured)
arm_expconc<-data.table::copy(exp_nom.dt[, nomconc := ExperimentalDose_uM])

armitageOutput_arm_expconc <- armitage_eval(tcdata=arm_expconc, 
                                            surface.area.switch = TRUE,
                                            user_assay_parameters = 
                                              user_assay_parameters,
                                            restrict.ion.partitioning = TRUE)

#Armitage NOMCONC RMSLE:
arm_nomconc_rmsle_long<-rmsle((armitageOutput_arm_nomconc$CellConcentration_uM),
                              (armitageOutput_arm_nomconc$ccells))
arm_nomconc_rmsle<-round(arm_nomconc_rmsle_long, 2)
#1.87

#Armitage EXPCONC RMSLE:
arm_expconc_rmsle_long<-rmsle((armitageOutput_arm_expconc$CellConcentration_uM),
                              (armitageOutput_arm_expconc$ccells))
arm_expconc_rmsle<-round(arm_expconc_rmsle_long, 2)
#1.75

#improvement between the experimental and nomconc rmsles
arm_nomconc_rmsle- arm_expconc_rmsle
#0.12

```

#Supplement:

## Table S1. 
Physicochemical properties of test chemicals and Ccell/Cnominal ratios 
from each reference.
```{r physchem_data_tableS1, eval = execute.vignette}

#create datatable to hold the info
table1.dt<-data.table::copy(Nominal.dt)

# pull out the columns we need
data.dt<-table1.dt[, c("Citation", "casrn", "ChemicalName", "IOC_Type", 
                       "nomconc", "CellConcentration_uM"), with = FALSE]


# pull dtxsids, log P, henry's law, and pka:

#dtxsids
cheminfo.dt<-as.data.table(get_chem_id(chem.cas=unique(data.dt$casrn)))

cheminfo.dt[, casrn:=chem.cas] %>% #rename
  .[,chem.cas:=NULL] %>% #delete
  .[,chem.name:=NULL] #delete

#add dtxsid info to the storage file
data.dt2<-data.dt[cheminfo.dt,on="casrn"]

#collect log P, henry's law, and pka values
data.dt2[, c("logP","logHenry", "pKa_Donor", "pKa_Accept") := 
           as.data.frame(get_physchem_param(param = c("logP","logHenry", 
                                                      "pKa_Donor", 
                                                      "pKa_Accept"), 
                                            chem.cas = casrn))] 


# calculate median CellConcentration_uM:cnom ratio 
# (for each citation if multiple) and number of observations per chemical

#get ccell : cnom ratio
data.dt2[,cellconc_cnom_ratio:=CellConcentration_uM/nomconc]

#for each chemical/citation combination, calculate the median ccell:cnom ratio
data.dt2[, median_ratio := median(cellconc_cnom_ratio), 
         by=list(Citation, ChemicalName)]

#count the number of observations by citation and chemname
data.dt2[, number_of_obs:=.N, by=list(Citation, ChemicalName)]

#filter down to one observation per citation/casrn combo
final_output<-data.dt2 %>% 
  mutate(logP=round(logP, 2),
         logHenry=round(logHenry, 2),
         median_ratio=round(median_ratio, 3)) %>% 
  select(-c("casrn", "nomconc", "CellConcentration_uM", 
            "cellconc_cnom_ratio")) %>% #get rid of extra rows
  distinct(Citation, ChemicalName, .keep_all = TRUE) 
#1 obs per citation/chem pair

#save the output
if(save_output){
  write.csv(final_output, file = paste0(path_out,"table_s1_",data.date,"_",
                                        Sys.Date(),".csv"))
}

```

## Table S2. In vitro assay descriptors for the measurements

```{r assay_descriptors_tableS2, eval = execute.vignette}

#copy table to manipulate
parameter_table<-data.table::copy(Nominal.dt)

table_s2<-parameter_table[, c("Citation", "ChemicalName", "nomconc", 
                              "NominalDose_uM", "v_total", "v_working", 
                              "sarea", "well_number", "cell_yield", "FBSf", 
                              "FigureNumber", "CellType")] #subset columns

#save the output
if(save_output){
  write.csv(table_s2, 
            file = paste0(path_out,"table_s2_raw_",data.date,"_",
                          Sys.Date(),".csv"))
}

```

## Table S3: 
Overall root mean squared log error (RMSLE) and coefficient of 
determination (r2) for the Armitage and Kramer models. 
See Figure 2 code above. 

## Table S4: 
RMSLE and r^2 for acids and bases.
```{r ionization_RMSLE_tableS4, eval = execute.vignette}

#copy for editing
ionization_rmsle<-data.table::copy(stacked.data)

#calculate RMSLSE and r^2 for each model and ionization type
ionization_rmsle_output<- data.table::copy(ionization_rmsle[label=="Armitage", 
                                                RMSLE_armitage_ionization:= 
                         rmsle(CellConcentration_uM, ccells), 
                       by = "IOC_Type"] %>% #RMSLE
  .[label=="Armitage", rsq_armitage_ionization:= 
      summary(lm(CellConcentration_uM~ccells))$r.squared, 
    by = "IOC_Type"] %>% #r^2
  .[label=="Kramer", RMSLE_kramer_ionization:= 
                         rmsle(CellConcentration_uM, concentration_cells), 
                       by = "IOC_Type"] %>% #RMSLE
  .[label=="Kramer", rsq_kramer_ionization:= 
      summary(lm(CellConcentration_uM~concentration_cells))$r.squared, 
    by = "IOC_Type"] %>% #r^2
  .[,numberofpoints:= .N/2, by = "IOC_Type"] %>% 
  .[, c("IOC_Type", "RMSLE_armitage_ionization",
        "rsq_armitage_ionization", "RMSLE_kramer_ionization", 
        "rsq_kramer_ionization", "numberofpoints")] %>% #select columns we want
  distinct())

print(ionization_rmsle_output)
#number of points is overall, not for each calculation

#save the output
if(save_output){
  write.csv(ionization_rmsle_output, 
            file = paste0(path_out,"table_s4_",data.date,"_",Sys.Date(),".csv"))
}

```

## Table S5:
RMSLE across lipophilicity.
```{r logP_RMSLE_tableS5, eval = execute.vignette}

#cut up into chunks
gkowbins<-stacked.data %>% mutate(LogP_bins = 
                                    cut(gkow, breaks = c(-Inf, 1.30, 3.36, Inf),
  labels = c("Hydrophilic", "Moderate Lipophilic", "Lipophilic")))

#calculate rmsle for each model and logP bin combo
gkowbins_rmsle_output<-suppressWarnings(data.table::copy(
  gkowbins[label == "Armitage", Armitage_RMSLE := rmsle(
  CellConcentration_uM, ccells), by = LogP_bins] %>% 
  .[label == "Kramer", Kramer_RMSLE := rmsle(
  CellConcentration_uM, concentration_cells), by = LogP_bins] %>% 
  .[, numberofpoints := .N, by = list(LogP_bins, label)] %>% 
  .[, c("LogP_bins", "Armitage_RMSLE", "Kramer_RMSLE", 
        "numberofpoints", "label")] %>% 
  distinct()))

print(gkowbins_rmsle_output)
#number of points is overall, not for each calculation

#save the output
if(save_output){
  write.csv(gkowbins_rmsle_output, 
            file = paste0(path_out,"table_s5_",data.date,"_",Sys.Date(),".csv"))
}

```

## Table S6. 
RMSLE across volatility.
```{r kaw_RMSLE_tableS6, eval = execute.vignette}

#source for chunks is table 2 Volatility class in water: 
#https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/
#guidance-reporting-environmental-fate-and-transport#:~:text=is%20very%20low.-,
#Volatility%20from%20Water,the%20reciprocal%20of%20KAW.&text=GMW%20=%20gram%20
#molecular%20weight%20in%20g/mole.&text=This%20value%20is%20the%20inverse,as%20
#Cwater/Cair.&text=Note%20that%20all%20chemicals%20may,volatility%20potential%20
#is%20very%20low.

#cut up kaw into relevant bins
kawbins<-stacked.data %>% 
  mutate(Kaw_bins = cut(DR_kaw, breaks = c(Inf, 10^-2, 10^-3, 10^-5, -Inf),
  labels = c("Rapidly lost from a water surface", 
             "Volatile from a water surface", 
             "Slightly volatile from a water surface", "Non-volatile")))

#calculate rmsle for each model and logP bin combo
kaw_rmsle_output<-suppressWarnings(kawbins[label == "Armitage", 
                                    Armitage_RMSLE := rmsle(
  CellConcentration_uM, ccells), by = Kaw_bins] %>% 
  .[label == "Kramer", Kramer_RMSLE := rmsle(
  CellConcentration_uM, concentration_cells), by = Kaw_bins] %>% 
  .[, numberofpoints := .N, by = list(Kaw_bins, label)] %>% 
  .[, c("Kaw_bins", "Armitage_RMSLE", "Kramer_RMSLE", 
        "numberofpoints", "label")] %>% 
  distinct())


#save the output
if(save_output){
  write.csv(kaw_rmsle_output, 
            file = paste0(path_out,"table_s6_",data.date,"_",Sys.Date(),".csv"))
}


```

## Table S7:
Root mean squared log error (RMSLE) and coefficient of determination (r2) 
for the chemicals with multiple observations. Note: Atrazine, Flusilazole, 
N-Phenyl-1,4-benzenediamine, Rosiglitazone, Thiacloprid, and Triphenyl phosphate
each only reflect two data points.
```{r multobs_RMSLE_tableS7, eval = execute.vignette}

#Armitage

#copy for editing
multobs_rmsle_armitage<-data.table::copy(as.data.table(Armitage.multobs.dt))

#group by chemname, calc RMSLE, r^2, and number of obs
multobs_rmsle_armitage[, RMSLE_armitage_mult:= 
                         rmsle(mean_measuredccell, mean_predictedccell), 
                       by = "ChemicalName"] %>% #rmsle
  .[, rsq_armitage_mult:= 
      summary(lm(mean_measuredccell~mean_predictedccell))$r.squared, 
    by = "ChemicalName"] %>% #r^2
  .[, numberofpoints:= .N, by = "ChemicalName"] 

#numberofpoints = number of datapoints considered (includes averaged points)

#repeat for kramer

#copy for editing
multobs_rmsle_kramer<-data.table::copy(as.data.table(Kramer.multobs.dt))

#group by chemname, calc RMSLE, r^2, and number of obs
multobs_rmsle_kramer[, RMSLE_kramer_mult:= 
                         rmsle(mean_measuredccell, mean_predictedccell), 
                       by = "ChemicalName"] %>% #rmsle
  .[, rsq_kramer_mult:= 
      summary(lm(mean_measuredccell~mean_predictedccell))$r.squared, 
    by = "ChemicalName"] %>% #r^2
  .[, numberofpoints:= .N, by = "ChemicalName"] 


#cut down both to the columns we need and only 1 obs per chemical
mult.dt<-merge(multobs_rmsle_armitage[, -c("NominalDose_uM", 
                                           "mean_measuredccell", 
                                           "mean_predictedccell")] %>% 
                 distinct() , 
               multobs_rmsle_kramer[, -c("NominalDose_uM", 
                                         "mean_measuredccell", 
                                         "mean_predictedccell")] %>% 
                 distinct())

#for higher/lower rmsle comparison
mult_HL<-data.table::copy(mult.dt[RMSLE_kramer_mult>RMSLE_armitage_mult, 
                      higherorlower:="Kramer higher"] %>% 
  .[RMSLE_kramer_mult<RMSLE_armitage_mult, 
    higherorlower:="kramer lower"])

mult_HL %>% count(higherorlower)
#armitage predictions have lower error for 6 chemicals
#kramer predictions are lower error for the other 4

#save the output
if(save_output){
  write.csv(mult.dt, 
            file = paste0(path_out,"table_s3_",data.date,"_",Sys.Date(),".csv"))
}

```

## Table S8:
Effect of "options" on rmsle in the Armitage model.
```{r options_tableS8, eval = execute.vignette}

#list of options that can be changed in the armitage model:

# 		- this.option.kbsa2 Use alternative bovine-serum-albumin partitioning
# 		- this.option.swat2 Use alternative water solubility correction
# 		- this.option.kpl2 Use alternative plastic-water partitioning model
# 		- option.plastic == TRUE (default) gives nonzero surface area
# 		- option.bottom == TRUE (default) includes surface area of the bottom of 
#       the well in determining sarea
# 		- restrict.ion.partitioning 
# 			§ False (default):  Treat entire chemical as neutral and allow it all 
#       to partition normally as such
# 			§ True: Restrict so only the neutral fraction of the chemical is 
#       available to partition normally but the ionized portion is scaled


##filter to just the ones we can run with the surface area calculation since 
#option.plastic and option.bottom require that function:
options.IVD_Manuscript_data_2<-
  options.IVD_Manuscript_data[!is.na(well_number), ] %>% 
  .[,sarea:=NULL]  #delete sarea so the surface area function will run 
#126 out of the original 153 observations
# length(options.IVD_Manuscript_data$Citation %>% unique()) #14 citations
# length(options.IVD_Manuscript_data_2$Citation %>% unique()) #12 citations
# 153-126 data points; 27 data points total removed

# create table to append the results to 
options.table <- NULL

# run basic armitage, only setting the surface area switch
option.basic <- armitage_eval(tcdata=options.IVD_Manuscript_data_2,
                              surface.area.switch = TRUE,
                              user_assay_parameters = user_assay_parameters)

option.basic_RMSLE<-rmsle((option.basic$CellConcentration_uM),
                          (option.basic$ccells))
#1.148

#append to table
options.table<-rbind(options.table, c(round(option.basic_RMSLE ,3), 
                                      "All Default"))

# run with this.option.kbsa2 == TRUE (default is false)
option.kbsa <- armitage_eval(tcdata=options.IVD_Manuscript_data_2,
                             this.option.kbsa2 = TRUE,
                             surface.area.switch = TRUE,
                             user_assay_parameters = user_assay_parameters)

option.kbsa_RMSLE<-rmsle((option.kbsa$CellConcentration_uM),
                         (option.kbsa$ccells))
#1.150

#append to table
options.table<-rbind(options.table, c(round(option.kbsa_RMSLE, 3), 
                                      "option.kbsa=TRUE"))


# run with this.option.swat2 == TRUE (default is false)
option.swat <- armitage_eval(tcdata=options.IVD_Manuscript_data_2,
                             this.option.swat2 = TRUE,
                             surface.area.switch = TRUE,
                             user_assay_parameters = user_assay_parameters)

option.swat_RMSLE<-rmsle((option.swat$CellConcentration_uM),
                         (option.swat$ccells))
#3.585

#append to table
options.table<-rbind(options.table, c(round(option.swat_RMSLE, 3), 
                                      "option.swat2=TRUE"))

# run with this.option.kpl2 == TRUE (Fischer) (default is false/kramer)
option.kpl <- armitage_eval(tcdata=options.IVD_Manuscript_data_2,
                             this.option.kpl2 = TRUE,
                            surface.area.switch = TRUE,
                            user_assay_parameters = user_assay_parameters)

option.kpl_RMSLE<-rmsle((option.kpl$CellConcentration_uM),(option.kpl$ccells))
#1.101

#append to table
options.table<-rbind(options.table, c(round(option.kpl_RMSLE, 3), 
                                      "option.kpl2=TRUE"))

# run with option.plastic == FALSE (default is true)
opt.plast_input<- data.table::copy(options.IVD_Manuscript_data_2)
opt.plast_input[,option.plastic := FALSE]

option.plast <- armitage_eval(tcdata=opt.plast_input,
                              surface.area.switch = TRUE,
                              user_assay_parameters = user_assay_parameters)

option.plast_RMSLE<-rmsle((option.plast$CellConcentration_uM),
                          (option.plast$ccells))
#1.178

#append to table
options.table<-rbind(options.table, c(round(option.plast_RMSLE, 3), 
                                      "option.plastic=FALSE"))

# run with option.bottom == FALSE (default is true)
option.bottom_input<- data.table::copy(options.IVD_Manuscript_data_2)
option.bottom_input[,option.bottom := FALSE]

option.bottom <- armitage_eval(tcdata=option.bottom_input,
                               surface.area.switch = TRUE,
                              user_assay_parameters = user_assay_parameters)

option.bottom_RMSLE<-rmsle((option.bottom$CellConcentration_uM),
                           (option.bottom$ccells))
#1.157

#append to table
options.table<-rbind(options.table, c(round(option.bottom_RMSLE,3), 
                                      "option.bottom=FALSE"))

# run with restrict.ion.partitioning == TRUE (default is FALSE)
option.restion <- armitage_eval(tcdata=options.IVD_Manuscript_data_2,
                             restrict.ion.partitioning = TRUE,
                             surface.area.switch = TRUE,
                              user_assay_parameters = user_assay_parameters)

option.restion_RMSLE<-rmsle((option.restion$CellConcentration_uM),
                            (option.restion$ccells))
#1.152

#append to table
options.table<-rbind(options.table, c(round(option.restion_RMSLE,3), 
                                      "restrict.ion.partitioning = TRUE"))

#save the output
if(save_output){
  write.csv(options.table, 
            file = paste0(path_out,"table_s7_",data.date,"_",Sys.Date(),".csv"))
}

```

## Table S9:
No code

## Table S10:
See Figure 5 code above

## Table S11:
Difference in RMSLE between liver and non-liver cells for 
the Armitage and Kramer models.
```{r cell_type_tableS11, eval = execute.vignette}

#copy for editing
celltype_rmsle<-data.table::copy(stacked.data)

#one set doesnt have cell type annotated, add here:
celltype_rmsle[assay_name == "Tox21 Experiment", cell_type := "MCF-7"]

#check all assays have cell type added (no nas)
unique(celltype_rmsle$cell_type)

#create bins for the cell types
celltype_rmsle[cell_type == "PRH" |
                 cell_type == "PHH" |
                 cell_type == "HepaRG" |
                 cell_type == "HepG2"
                 , cell_type_class:="Liver"] %>% 
  .[cell_type == "Neocortical" |
     cell_type == "PFC" |
     cell_type == "Caco-2" |
     cell_type == "LS180" |
     cell_type == "MCF-7" |
     cell_type == "H295R" |
     cell_type == "PC12" |
     cell_type == "RTgill W1" |
     cell_type == "HEK293" |
     cell_type == "Balb/c 3T3", cell_type_class:="Non-Liver"]

#check all are labeled
celltype_rmsle %>% select(cell_type, cell_type_class) %>% distinct()

#calculate RMSLSE and r^2 for each model and ionization type
celltype_rmsle_output<- data.table::copy(celltype_rmsle[label=="Armitage", 
                                                RMSLE_armitage_cell_type:= 
                         rmsle(CellConcentration_uM, ccells), 
                       by = "cell_type_class"] %>% #RMSLE
  .[label=="Armitage", rsq_armitage_cell_type:= 
      summary(lm(CellConcentration_uM~ccells))$r.squared, 
    by = "cell_type_class"] %>% #r^2
  .[label=="Kramer", RMSLE_kramer_cell_type:= 
                         rmsle(CellConcentration_uM, concentration_cells), 
                       by = "cell_type_class"] %>% #RMSLE
  .[label=="Kramer", rsq_kramer_cell_type:= 
      summary(lm(CellConcentration_uM~concentration_cells))$r.squared, 
    by = "cell_type_class"] %>% #r^2
  .[,numberofpoints:= .N/2, by = "cell_type_class"] %>% 
  .[, c("cell_type_class", "RMSLE_armitage_cell_type",
        "rsq_armitage_cell_type", "RMSLE_kramer_cell_type", 
        "rsq_kramer_cell_type", "numberofpoints")] %>% #select columns we want
  distinct())

print(celltype_rmsle_output)
#number of points is overall, not for each calculation

```

## Figure S1:
Histogram showing chemical-specific RMSLE for Armitage Ccell vs Measured Ccell 
and Kramer Ccell vs Measured Ccell:
```{r intracellular_histogram_plotS1, eval = execute.vignette}

#calculate RMSLE for each chemical
combodata[, RMSLE_armitage_by_chem:= rmsle(CellConcentration_uM.x, ccells), 
          by = "ChemicalName.x"]
combodata[, RMSLE_kramer_by_chem:= 
            rmsle(CellConcentration_uM.x, concentration_cells), 
          by = "ChemicalName.x"]
combodata[, RMSLE_nominal_by_chem:= rmsle(CellConcentration_uM.x, nomconc), 
          by = "ChemicalName.x"]

#cut down to one observation per chemical
histo.data<-data.table::copy(combodata %>% 
                   select(ChemicalName.x, RMSLE_armitage_by_chem, 
                          RMSLE_kramer_by_chem, RMSLE_nominal_by_chem) %>% 
                   unique())

#histo.data %>% arrange(desc(RMSLE_kramer_by_chem))

#plot histograms
nominal_histogram_by_chem <-ggplot(histo.data, aes(x=RMSLE_nominal_by_chem)) + 
  geom_histogram(binwidth = 0.5, boundary = 0) +
  labs(x = "Nominal RMSLE",
       y = "Number of Observations") + 
  xlim(0, 5) + ylim(0, 17) +
  theme_classic() +  
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14))

armitage_histogram_by_chem <-ggplot(histo.data, aes(x=RMSLE_armitage_by_chem)) + 
  geom_histogram(binwidth = 0.5, boundary = 0) +
  labs(x = "Armitage RMSLE",
       y = "Number of Observations") +
  xlim(0, 5) + ylim(0, 17) +
  theme_classic() + 
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14))

kramer_histogram_by_chem <- ggplot(histo.data, aes(x=RMSLE_kramer_by_chem)) + 
  geom_histogram(binwidth = 0.5, boundary = 0) +
  labs(x = "Kramer RMSLE",
       y = "Number of Observations") +
  xlim(0, 5) + ylim(0, 17) +
  theme_classic()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14))

histogram_comboplot_by_chem<-ggarrange(nominal_histogram_by_chem +
                                         theme(axis.ticks.y = element_blank(),
                                               plot.margin = margin(r = 1)), 
                                       armitage_histogram_by_chem +
                                         theme(axis.text.y = element_blank(),
                                               axis.ticks.y = element_blank(),
                                               axis.title.y = element_blank(),
                                               plot.margin = margin(r = 1, 
                                                                    l = 1)), 
                                       kramer_histogram_by_chem +
                                         theme(axis.text.y = element_blank(),
                                               axis.ticks.y = element_blank(),
                                               axis.title.y = element_blank(),
                                               plot.margin = margin(r = 1, 
                                                                    l = 1)), 
                                       ncol=3)

if(save_output){
  ggsave(file = paste0("Figure_S1_by_chem_",data.date,"_",Sys.Date(),".png"),
         path = path_out,
         histogram_comboplot_by_chem, width = 10, height = 5, dpi = 300)
}

print(histogram_comboplot_by_chem)

```
