Meade (submitted): High Throughput Dermal Exposure Model

Annabel Meade, Celia Schacht, Alex George

2025-07-23

Please send questions to or

from “Development and analysis of high throughput physiologically based pharmacokinetic/toxicokinetic (PBPK/TK) dermal exposure model”

Annabel E. Meade, Celia M.Schacht, Marina V Evans, Alex George, Rachael Cogbill, Risa R. Sayre, and John F. Wambaugh

Abstract

Dermal absorption of chemicals represents an important route of exposure in pharmaceutical, occupational, and environmental settings. There are thousands of chemicals in use with little or no toxicity or toxicokinetic data, and dermal absorption data are lacking, unsurprisingly, for many of these chemicals. One alternative is to estimate human toxicokinetics using high-throughput methods. Accordingly, the aim of this study was to develop a generalized physiologically based pharmacokinetic/toxicokinetic (PBPK/TK) dermal exposure model for the R package httk. For dermal exposures this model can be used in a high-throughput manner to estimate human blood and tissue concentrations and estimate risk for many chemicals. The structure of the dermal PBPK/TK model was based on Campbell, Clewell, Gentry, Anderson, and Clewell, Computational Toxicology, 929, 439 (2012). Chemical-specific metabolism and protein binding data were obtained from the literature as collected by R package httk. The physiochemical property data needed for the model were obtained from the EPA CompTox Chemicals Dashboard. Here we have compared two different in vitro methods for estimating chemical-specific and vehicle-specific permeability in the skin: 1) Potts-Guy [Potts and Guy, Pharm Res, 9, 663 (1992)] and 2) Surrey [Chen, Han, Saib, and Lian, Pharm Res, 32, 1779 (2015)]. The model was constructed to allow comparison of dermal, oral, and intravenous exposures. Over 26 exposure scenarios across 13 chemicals were modeled and compared to published concentration-time in vivo data from the EPA CvT database [Sayre, Wambaugh, and Grulke, Scientific Data, 7, 122 (2020)]. Of these 13 chemicals, four are pharmaceuticals while the other 9 are occupationally or environmentally relevant, 9 of the 13 are lipophilic (log of the octanol to water partition coefficient, logP > 2), three of the 13 have low water solubility (< 10 mg/L), four of the 13 have high water solubility (> 1,000 mg/L), and eight of the 13 are considered volatile. The Root Mean Squared Log10 Error (RMSLE) between log-transformed simulated and observed concentrations was calculated to determine how well the model captures the behavior of this data when using the methods for dermal permeability. Based on these results, the Potts-Guy methods slightly out-perform the Surrey method of calculating dermal permeability for more volatile chemicals and the vehicle is a large factor in method success, where Potts-Guy/octanol and Surrey/water combinations generated the best fits. We also found that exhalation is an important route of excretion for volatiles, in particular. The views expressed in this abstract are those of the authors and do not necessarily represent the views or the policies of the U.S. Environmental Protection Agency.

HTTK Version

This vignette was created with httk v2.3.3. 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/

Table of Contents (added by AG 01/02/2024)

Prepare for session Result 1: Model Verification Usint In Vivo Data Step One: Load Data Display which chemicals are being looked at Step Two: Clean and Format Data

Prepare for session

eval = FALSE

If you are using the RMarkdown version of this vignette (has the file extension .RMD) you will be able to see that several chunks of code in this vignette have the statement eval = FALSE. The next chunk of code, by default, sets run.simulations = 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 CRAN has limits 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, open the .RMD file in RStudio and either change run.simulations to TRUE or press “play” (the green arrow) on each chunk in RStudio.

We set a variable run.simulations. This will run simulations and regenerate data files. If this is set to TRUE, make sure you have the following data file:

and if run.simulations = FALSE, make sure to have the following files in your working directory:

# Set whether or not the following chunks will be executed (run):
run.simulations <- FALSE # Change FALSE to TRUE here to make the vignette work
make.plots <- FALSE # change to TRUE to make the plots
run.simulations.2 <- FALSE #for second part of vignette
make.plots.2 <- FALSE # change to TRUE to make the plots

Set working directory and install packages

Here, we set our working directory (the location of this file and the necessary data files) and load necessary packages. Please reset your working directory below to loc.wd or the code will produce an error.

We also set the location of the R package httk using the name loc.httk so that we can run the functions modelinfo_dermal.R and modelinfo_dermal_1subcomp.R file later in the code. Note that if you have run.simulations = FALSE, the variable loc.httk is not used, so you do not need to provide or change it.

# Set working directory
loc.wd <- "C:/Users/jwambaug/OneDrive - Environmental Protection Agency (EPA)/Profile/Documents/Research Projects/DermalTK/vingette"
#loc.wd <- "/Users/ameade/Library/CloudStorage/GoogleDrive-aemeade7@gmail.com/My Drive/EPA/HTTK-Dermal/HTTK-Dermal-Vignette"

#loc.wd = "C:/Users/cschacht/OneDrive - Environmental Protection Agency (EPA)/Profile/Documents/CCTE_2023_2024/dermal/thisone"

# Set location of httk package
#loc.httk <- "C:/Users/jwambaug/git/httk-dev/httk"


# Load packages
packages <- c("httk","dplyr",
              "tidyverse","xlsx","Metrics", #for data manipulation
              "ggplot2","ggforce","ggpubr","ggrepel","viridis", #for plotting
              "knitr", "ggpubr","grid","ggh4x","readr","ggforce","knitr","tidyr",
              "stringr") #to make pretty visual tables

sapply(packages, require, character.only=TRUE)
#> Loading required package: dplyr
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
#> Loading required package: tidyverse
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ forcats   1.0.0     ✔ stringr   1.5.1
#> ✔ lubridate 1.9.4     ✔ tibble    3.2.1
#> ✔ purrr     1.0.4     ✔ tidyr     1.3.1
#> ✔ readr     2.1.5     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#> Loading required package: xlsx
#> 
#> Loading required package: Metrics
#> 
#> Loading required package: ggforce
#> 
#> Loading required package: ggpubr
#> 
#> Loading required package: ggrepel
#> 
#> Loading required package: viridis
#> 
#> Loading required package: viridisLite
#> 
#> Loading required package: knitr
#> 
#> Loading required package: grid
#> 
#> Loading required package: ggh4x
#>      httk     dplyr tidyverse      xlsx   Metrics   ggplot2   ggforce    ggpubr 
#>      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE 
#>   ggrepel   viridis     knitr    ggpubr      grid     ggh4x     readr   ggforce 
#>      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE 
#>     knitr     tidyr   stringr 
#>      TRUE      TRUE      TRUE

# Make sure we're using a clean version of the httk chemical library:
#reset_httk()

Result 1: Model Verification using In Vivo Data

First, we want to verify our model by testing it against in vivo data from literature.

Step One: Load Data

We have some data which is taken from the CvT Database (Sayre, Wambaugh, and Grulke, 2020) and some data which has been prepared for the database but is not yet added to the Database, denoted as noSQL since it is not yet in the SQL space. The CvT Database contains structured data in the form of separate datasets linked with foreign keys (“fk”). Below is the entity-relationship for each of the databases.

To learn more about where these datasets come from, please see R_CvT_Database folder, which contains files used to generate the data below.

# Load CvT table database
load(paste0(loc.wd,"/R_CvT_Database/CvT_all_tables_2022-08-17.Rdata")) #Necessary file that is not included yet
  # Load CvT data that is not yet in CvT database (Annabel did this)
  load(paste0(loc.wd,"/R_CvT_Database/CvT_all_tables_2023-07-05_noSQL.Rdata")) #Necessary file that is not included yet
  
  #Add cumulative_amount column in SQL DB_series (we don't use any cumulative amounts from SQL data in final results)
  DB_series <- DB_series %>% mutate(cumulative_amount=NA)

# Bind data sets from database and noSQL collection.
DB_documents_noSQL<-DB_documents_noSQL[colnames(DB_documents_noSQL) %in% colnames(DB_documents)] #Keep only columns from DB_documents_noSQL that are also found in DB_documents
DB_documents <- rbind(DB_documents,DB_documents_noSQL)
DB_studies_noSQL<-DB_studies_noSQL[colnames(DB_studies_noSQL) %in% colnames(DB_studies)]
DB_studies <- rbind(DB_studies,DB_studies_noSQL)
DB_chemicals <- rbind(DB_chemicals,DB_chemicals_noSQL)
DB_subjects <- rbind(DB_subjects,DB_subjects_noSQL)
DB_series <- rbind(DB_series,DB_series_noSQL)
DB_conc_time_values_noSQL<-DB_conc_time_values_noSQL[colnames(DB_conc_time_values_noSQL) %in% colnames(DB_conc_time_values)]
DB_conc_time_values <- rbind(DB_conc_time_values,DB_conc_time_values_noSQL)

Next, we need to clean the data and filter out NA results. First, we join the needed data sets using their foreign keys. * DB_conc_time_values * DB_series * DB_studies * DB_subjects * DB_documents: Note that only the PubMed ID (PMID) of each extraction article is saved from this data set.

Next, we remove columns that are duplicated in each data set and are not used here. These columns in each data set, listed below, tell the user who created and updated the SQL files and additional comments by the data curator or article author. * created_by * updated_by * rec_update_dt * rec_created_dt * curator_comment * author_comment ( not in DB_conc_time_values or DB_subjects) We also remove any data that: * does not have a defined DTXSID for the chemical * does not have a recorded concentration * does not have a defined recorded medium (e.g., plasma, blood, urine, etc.) * where the analyte measured does not match the test substance used to dose the individual

Finally, after cleaning this data set, called df.full, we filter so that we only have data for chemicals that have dermal data. Our final data set is called df.

# Join necessary DBs: conc_time_values, series, studies, subjects
df.full <- DB_conc_time_values %>% 
  left_join(DB_series,by=c("fk_series_id"="id")) %>%
  left_join(DB_studies,by=c("fk_study_id"="id")) %>%
  left_join(DB_subjects,by=c("fk_subject_id"="id")) %>%
  left_join(DB_documents[,c("id","pmid")],by=c("fk_extraction_document_id" = "id")) %>%
  dplyr::rename(extraction_document_pmid=pmid) %>%
  select(!(ends_with(".x") | ends_with(".y"))) %>%
  filter(!is.na(test_substance_dtxsid)) %>% #remove series w/o defined chemicals (Changed from study to series AG)
  filter( !is.na(conc)) %>% #filter out NA data
  filter(!is.na(conc_medium_normalized)) %>% #remove series w/o defined medium (i.e., plasma, blood, urine, etc.)
  filter(test_substance_dtxsid==analyte_dtxsid | is.na(analyte_dtxsid)) #analyte (measured substance) and test substance must be the same

df.full <- distinct(df.full)

# Remove empty cells
df.full <- type_convert(df.full)
df.full <- df.full %>% filter(!is.na(conc))

# Get list of chemicals that have dermal data
df.dermal <- df.full %>% filter(administration_route_normalized=="dermal")
dermal_dtxsid <- unique(df.dermal$test_substance_dtxsid)

# Only keep data for chemicals that have dermal data
df <- df.full %>% filter(test_substance_dtxsid %in% dermal_dtxsid)

Display which chemicals are being looked at

table <- df %>% mutate(Chemical = get_chem_id(dtxsid=test_substance_dtxsid)$chem.name) %>% dplyr::count(Chemical,test_substance_dtxsid,extraction_document_pmid,sort=F) %>% dplyr::rename(Number.of.Data.Points = n)
row.names(table) <- 1:nrow(table)
knitr::kable(table, row.names=TRUE,caption="List of chemicals with dermal data before cleaning data.")
write.table(table,
            file=paste0(loc.wd,"/tables/ChemicalswDermalData.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")

Next, We look at several chemical factors with respect to our data. Thus, we ran a batch search of the 13 chemicals in our study in the CompTox Dashboard and now load information about these chemicals (Boiling point, Volatility, etc.) and save this as ‘supptab1_meade2023’.

# Load chemical info
supptab1_meade2023 <- read.xlsx(file=paste0(loc.wd,"/CCD-Batch-Search_2022-11-29_07_06_35.xlsx"),
                                sheetName="Main Data") #A separate download was needed for this file AG 12/6/23

# Rename chemical properties and change units
supptab1_meade2023 <- supptab1_meade2023 %>% mutate(Water.Solubility.mg.L = WATER_SOLUBILITY_MOL.L_OPERA_PRED*AVERAGE_MASS*1000,
                                                    Vapor.Pressure.mmHg = VAPOR_PRESSURE_MMHG_OPERA_PRED,
                                                    Boiling.Point.C = BOILING_POINT_DEGC_OPERA_PRED,
                                                    logP = OCTANOL_WATER_PARTITION_LOGP_OPERA_PRED,
                                                    MW = AVERAGE_MASS,
) %>%
  select(c(PREFERRED_NAME,DTXSID,Water.Solubility.mg.L,Vapor.Pressure.mmHg,Boiling.Point.C,logP,MW))

# Add column for "Volatility
supptab1_meade2023 <- supptab1_meade2023 %>% mutate(Volatility = ifelse(Boiling.Point.C<=75,"Very Volatile",
                                                                        ifelse(Boiling.Point.C>400,"Not Volatile",
                                                                               ifelse((Boiling.Point.C>=250)&(Boiling.Point.C<=400),"Semi-Volatile","Volatile"))),
                                                    # Rename long chemicals to nick-names
                                                    PREFERRED_NAME = ifelse(PREFERRED_NAME=="Methyl tert-butyl ether","MBTE",PREFERRED_NAME),
                                                    PREFERRED_NAME = ifelse(PREFERRED_NAME=="3,3',5,5'-Tetrabromobisphenol A","Bromdian",PREFERRED_NAME),
                                                    PREFERRED_NAME = ifelse(PREFERRED_NAME=="4,4'-Sulfonyldiphenol","Bisphenol S",PREFERRED_NAME),
                                                    PREFERRED_NAME = ifelse(PREFERRED_NAME=="Dichloromethane","DCM",PREFERRED_NAME),
                                                    PREFERRED_NAME = ifelse(PREFERRED_NAME=="Tetrachloroethylene","PERC/TCE",PREFERRED_NAME),
                                                    PREFERRED_NAME = ifelse(PREFERRED_NAME=="1-Methylbenzene","Toluene",PREFERRED_NAME))

supptab1_meade2023$Volatility <- factor(supptab1_meade2023$Volatility, levels=c("Not Volatile","Semi-Volatile","Volatile","Very Volatile"))
#save chemical dataframe
save(supptab1_meade2023,file=paste0("meade2023_",format(Sys.time(), "%b_%d_%Y"),".Rdata"))

supptab1_meade2023[,c(3:7)]=signif(supptab1_meade2023[,c(3:7)],3)
knitr::kable(supptab1_meade2023, 
             #row.names=TRUE,
             caption="Chemical properties",
             floating.environment="sidewaystable",
             align="c")
write.table(supptab1_meade2023,
            file=paste0(loc.wd,"/tables/ChemicalProperties.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")

Step Two: Clean and Format Data

Next we clean the data by looking only at data that can be simulated using the httk dermal model. First, we load the most recent QSAR in silico data (Dawson et al., 2021) from httk. #Add a statement about why QSAR values are needed - there are many chemicals for which httk does not have experimentally derived Clint and Fup. - AG 01/01/2024

load_dawson2021() #get QSAR data, overwrite for p-values greater than 0.05

Next, we only load data for routes iv, oral, and of course dermal. We also must remove the chemicals: dibromomethane, halothane, and isoflurane, because these chemicals do not have in silico or in vitro data to inform necessary parameters for the httk dermal model.

# Select route (dermal/iv/oral)
route.ls = c("iv","dermal","oral")

# Reduce data based on available in vitro data
# df.reduced <- df %>% filter(administration_route_normalized %in% route.ls) %>% #only iv, oral, dermal routes
#   filter(test_substance_dtxsid %!in% c("DTXSID4021557","DTXSID4025371","DTXSID3020752")) #no metabolism data for: dibromomethane, halothane, isoflurane
df.reduced <- df %>% filter(administration_route_normalized %in% route.ls) %>% #only iv, oral, dermal routes
  filter(!(test_substance_dtxsid %in% c("DTXSID4021557","DTXSID4025371","DTXSID3020752"))) #no metabolism data for: dibromomethane, halothane, isoflurane

df.reduced <- type_convert(df.reduced) #change dose_level_normalized to numeric 

#Join df.reduced with the Volatility/chemical info from supptab1_meade2023
df.reduced = df.reduced %>% left_join(supptab1_meade2023,by=c("test_substance_dtxsid"="DTXSID"))

df.reduced=subset(df.reduced, !(conc_medium_normalized %in%c("feces","urine")))
# Display which chemicals are being looked at
table <- df.reduced %>% 
  mutate(Chemical = get_chem_id(dtxsid=test_substance_dtxsid)$chem.name) %>% 
  dplyr::count(Chemical,test_substance_dtxsid,extraction_document_pmid,conc_medium_normalized,sort=F) %>% 
  dplyr::rename("Number.of.Data.Points" = n)
row.names(table) <- 1:nrow(table)
knitr::kable(table, 
             row.names=TRUE,
             caption="List of chemicals with dermal data after cleaning data.")
write.table(table,
            file=paste0(loc.wd,"/tables/ChemicalswDermalDataCleaned.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")

Next, we clean up the data so we can match experimental conditions to simulation conditions. We start by normalizing concentration and concentration units (where concentration refers to the experimental concentration vs. time data)

# Set dummy column for Concentration units 
df.reduced$conc_normalized = NA
df.reduced$conc_units_normalized = df.reduced$conc_units_original

allset = which(df.reduced$conc_units_original %in% c("ug/mL","ug/g","mg/kg", "mcg/mL","mg/L"))
df.reduced$conc_units_normalized[allset]="mg/L"; df.reduced$conc_normalized[allset]=df.reduced$conc_original[allset]

ngml = which(df.reduced$conc_units_original=="ng/ml");df.reduced$conc_normalized[ngml]=  df.reduced$conc_original[ngml]*0.001
df.reduced$conc_units_normalized[ngml]="mg/L"

ngl = which(df.reduced$conc_units_original=="ng/l");df.reduced$conc_normalized[ngl]=  df.reduced$conc_original[ngl]*1e-06
df.reduced$conc_units_normalized[ngl]="mg/L"

molunts = which(df.reduced$conc_units_original %in% c("uM","nmol/L","umol/L"))
dtx=unique(df.reduced[molunts,c("test_substance_dtxsid","conc_units_original")])
dtx$mw = get_physchem_param(param="MW", dtxsid=dtx$test_substance_dtxsid)

for(z in 1:dim(dtx)[1]){
  these.rows=which(df.reduced$test_substance_dtxsid==dtx$test_substance_dtxsid[z] & df.reduced$conc_units_original==dtx$conc_units_original[z])
  conversion_val = convert_units(input.units = dtx$conc_units_original[z],output.units = "mg/L",dtxsid=dtx$test_substance_dtxsid[z])
  df.reduced[these.rows,"conc_normalized"] = df.reduced[these.rows,"conc_original"] * conversion_val
  df.reduced$conc_units_normalized[these.rows] = "mg/L"
}

perc.dose=which(df.reduced$conc_units_original %in%c("% dose", "% dose/h"))
df.reduced[perc.dose, "conc_normalized"]=df.reduced[perc.dose, "conc_original"]*df.reduced[perc.dose, "dose_level_original"]
#df.reduced[perc.dose,"conc_units_normalized"]=df.reduced[perc.dose,"dose_level_original_units"]
df.reduced[perc.dose,"conc_units_normalized"]="mg/L"                                                                                       
# ppbv
ppbv = which(df.reduced$conc_units_original=="ppbv")
df.reduced[ppbv, "conc_normalized"] = convert_units(input.units = "ppbv",output.units = "mg/L",dtxsid=unique(df.reduced$analyte_dtxsid[ppbv]),state="gas") * df.reduced[ppbv, "conc_original"]
df.reduced$conc_units_normalized[ppbv] = "mg/L"

We also convert the dose units from the database to be in amounts (not concentrations) using the vehicle volumes. After doing this, all unique dosing scenarios are listed in a table.

# Convert ug/mL to mg
df.reduced <- df.reduced %>% 
  separate(dose_volume,
           into=c("dose_volume","dose_volume_units"),
           sep=" ") %>% #separate dose volume into number and units, units either NA or ml/mL
  mutate(dose_volume=as.numeric(dose_volume)) %>%
  mutate(dose_level_new = 
           ifelse(tolower(dose_level_original_units)=="ug/ml",
                  dose_level_original*dose_volume/1e+3, # ug >- mg, vol in mL
                  dose_level_original), #convert ug/ml to mg
         dose_level_new_units = 
           ifelse(tolower(dose_level_original_units)=="ug/ml",
                  "mg",
                  dose_level_original_units))


# Convert ppm units to mg/L (assuming ppm is in air) - add MW to dataframe to do this
df.reduced <- df.reduced %>% 
  left_join(chem.physical_and_invitro.data[,c("DTXSID","MW")], 
            by=c("test_substance_dtxsid"="DTXSID"))
df.reduced$MW.y=NULL
colnames(df.reduced)=gsub("MW.x","MW",colnames(df.reduced)) #fix MW 
df.reduced  = df.reduced %>% #
  mutate(MW = as.numeric(MW)) %>%
  group_by(test_substance_dtxsid) %>%
  mutate(dose_level_new = 
           ifelse(dose_level_original_units=="ppm" ,
                  dose_level_new *
                    convert_units(input.units="ppmw", #?? whatever
                                  output.units = "mg/l",
                                  MW = unique(MW)),
                  dose_level_new),
         # dose_level_new = 
         #   ifelse(dose_level_original_units=="ppm" & administration_route_original!="dermal vapor",
         #          dose_level_new *
         #            convert_units(input.units="ppmw", #?? whatever
         #                          output.units = "mg/l",
         #                          MW = unique(MW), state="liquid"),
         #          dose_level_new),
         dose_level_new_units = 
           ifelse(tolower(dose_level_original_units)=="ppm",
                  "mg/L",
                  dose_level_new_units))
df.reduced$dose_level_new_units[which(df.reduced$dose_level_new_units=="ppm")]="mg/L"

# Set normalized vehicle
df.reduced <- df.reduced %>% 
  mutate(dose_vehicle_normalized = 
           ifelse(grepl("corn oil",dose_vehicle) | 
                    grepl("ethanol",dose_vehicle),"octanol", 
                  ifelse(!is.na(dose_vehicle) & 
                           !grepl("Cookie",dose_vehicle),
                         "water",
                         NA)))

# Remove added chemical MW column for now
df.reduced <- df.reduced %>% select(!MW)

# Display dosing scenarios
table <- df.reduced %>% 
  mutate(dose_level_original=signif(dose_level_original,3),
         dose_level_new = signif(dose_level_new,3),
         Chemical = get_chem_id(dtxsid=test_substance_dtxsid)$chem.name) %>% 
  count(Chemical,dose_level_original,
        dose_level_original_units,
        dose_level_new,
        dose_level_new_units,
        dose_volume,
        dose_volume_units,
        dose_vehicle,
        dose_vehicle_normalized,
        dermal_applied_area,
        dermal_applied_area_units,sort=F) %>% 
  rename("Number.of.Data.Points" = n)
row.names(table) <- 1:nrow(table)
knitr::kable(table, 
             #row.names=TRUE,
             caption="Dosing Scenarios in CvT Data",
             floating.environment="sidewaystable",
             align="c")

Display which chemicals are being looked at

table <- df %>% mutate(Chemical = get_chem_id(dtxsid=test_substance_dtxsid)$chem.name) %>% dplyr::count(Chemical,test_substance_dtxsid,extraction_document_pmid,sort=F) %>% dplyr::rename(Number.of.Data.Points = n)
row.names(table) <- 1:nrow(table)
knitr::kable(table, row.names=TRUE,caption="List of chemicals with dermal data before cleaning data.")

Now, update certain experimental parameters. If body weight is not reported, set it as the population mean, calculate skin area, and specify the model output name that corresponds to each medium. Also, we need to extract dose duration information for explicit values for the length of duration and the length of duration units (in hours, making conversions when necessary).

df.new=df.reduced

df.new$species = tolower(df.new$species)
# get dose duration in numbers and time units
df.new$dose_duration_num = as.numeric(str_extract(df.new$dose_duration, "[0-9]+"))
df.new$dose_duration_units = str_extract(df.new$dose_duration, "[aA-zZ]+")
df.cases <- df.new %>% mutate(BW = ifelse(!is.na(weight_kg),weight_kg,
                                                    case_when(
                                                      species=="rat" ~ 0.4,
                                                      species=="human" ~ 70,
                                                      species=="mouse" ~ 0.025,
                                                      species=="rabbit" ~ 1.6
                                                    )),
                                        # Set human default height to 175 cm if not listed
                                        height_cm = ifelse(species=="human" & is.na(height_cm),175,height_cm),
                                        totalSA = case_when( 
                                          species=="rat" ~ 9.83*(BW*1000)^(2/3), # formula from Gouma et al., 2012, 195-240grams
                                          species=="human" ~ sqrt(height_cm * BW/3600) * 100^2,
                                          species=="mouse" ~ 9.83*(BW*1000)^(2/3),
                                          species=="rabbit" ~ 100*11*(BW)^(2/3)  # formula from Tadashi et al., 2018
                                        ),
                                        #Assume Fskin_exposed = 10% if not stated
                                        Fskin_exposed = ifelse(!is.na(dermal_applied_area),dermal_applied_area/totalSA,0.1),
                                        #Infinite Dose if dosing units are concentration mg/l AJG 4/13/24: Double check that mg/l are listed in the correct places in the previous chunk.
                                        InfiniteDose = ifelse(dose_level_new_units=="mg/L",1,0), # omg
    
                                        Vvehicle = ifelse(tolower(dose_volume_units)=="ml",dose_volume/1e3,dose_volume),
                                        #Standardize time units for dose duration AJG 5/7/24
                                        dose_duration_units = case_when(
                                          dose_duration_units %in% c("hr","h","hour","hrs","hours") ~ "hours",
                                          dose_duration_units %in% c("min","m","minute","minutes","mins") ~ "mins",
                                          dose_duration_units %in% c("day","d","days") ~ "days"
                                        ),
                              dose_duration_num2 = ifelse(dose_duration_units=="mins", dose_duration_num/60, dose_duration_num),
                              dose_duration_units2 = ifelse(dose_duration_units=="mins", "hours", dose_duration_units)
                              
                              
)
# remove experiments where conc_medium_normalized = feces and set observation to collect
df.cases = subset(df.cases, !(conc_medium_normalized %in%c("urine","feces")))
df.cases = df.cases%>%
  mutate(conc_medium_obs = recode(conc_medium_normalized,
                                  "blood"="Cven",
                                  "liver"="Cliver",
                                  "lung"="Clung",
                                  "kidney"="Ckidney",
                                  "plasma"="Cplasma",
                                  "serum"="Cplasma",
                                  "expired air" = "Aexhaled"
  ))


# Display dosing scenarios
table <- df.cases %>% 
  mutate(dose_level_original=signif(dose_level_original,3),
         dose_level_new = signif(dose_level_new,3),
         Chemical = get_chem_id(dtxsid=test_substance_dtxsid)$chem.name) %>% 
  count(Chemical, species,dose_level_original,
        dose_level_original_units,
        dose_level_new,
        dose_level_new_units,
        dose_volume,
        dose_volume_units,
        dose_vehicle,
        dose_vehicle_normalized,
        dose_duration_num2,
        dermal_applied_area,
        dermal_applied_area_units,
        extraction_document_pmid,sort=F)
row.names(table) <- 1:nrow(table)
knitr::kable(table, 
             #row.names=TRUE,
             caption="Dosing Scenarios in CvT Data",
             floating.environment="sidewaystable",
             align="c")
write.table(table,
            file=paste0(loc.wd,"/tables/CvTdbDosingScenarios.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")

Finally, we clarify exhalation (depending on chemical volatility) and set vehicle volume parameters and save the entire cleaned datasheet as concentration_data_meade2023.

df.cases_new = df.cases
df.cases_new$Kvehicle2water = df.cases_new$dose_vehicle_normalized # set this parameter
df.cases_new =  df.cases_new %>% mutate(
  administration_route_normalized = case_when(
    administration_route_normalized == "oral" ~ "Oral",
    administration_route_normalized == "iv" ~ "IV",
    administration_route_normalized == "dermal" ~ "Dermal"
  ),
    #Add exhalation for volatile compounds
    Exhalation = case_when(
      levels(Volatility)[Volatility] == "Very Volatile" ~ TRUE,
      levels(Volatility)[Volatility] == "Volatile" ~ TRUE, 
      levels(Volatility)[Volatility] == "Semi-Volatile" ~ TRUE,
      levels(Volatility)[Volatility] == "Not Volatile" ~ FALSE
    ),
    # Set default vehicle volume when it is unknown. Assume a uniform 1 mm layer of vehicle on skin. Then there is 0.0001 L vehicle per cm^2 of exposed skin. For cases with non-dermal routes of exposure, these values are NA.
    Vvehicle = case_when(
      administration_route_normalized=="Dermal" ~ ifelse(!is.na(Vvehicle), Vvehicle, Fskin_exposed*totalSA*0.0001),
      administration_route_normalized=="Oral" ~ NA,
      administration_route_normalized=="IV" ~ NA)
   
)

# save what we want
concentration_data_meade2023 = df.cases_new
cgwtools::resave(concentration_data_meade2023,file=paste0("meade2023_",format(Sys.time(), "%b_%d_%Y"),".Rdata"))

Step Three: Simulate Data

Now, we will simulate the in vivo concentration vs time data using the httk dermal model. Using concentration_data_meade2023, which is clean and filtered above, we create a new data set called df.simulations which will contain simulation data. We simulate 6 different dermal scenarios: * Potts-Guy, vehicle=octanol * Potts-Guy, vehicle=water * Surrey, vehicle=octanol * Surrey, vehicle=water * 2 Compartment, vehicle=octanol * 2 Compartment, vehicle=water

Simulation Set-Up

Specify the units we want returned and set up simulation data.


# Set output units for dermal model to be same as data
output.units.original <- c( "Agut"="mg", "Agutlumen"="mg", "Aliver"="mg", "Aven"="mg", "Alung"="mg", "Arest"="mg", "Akidney"="mg", "Atubules"="mg", "Ametabolized"="mg", "Avehicle"="mg", "Ain"="mg",  
                            "Cgut"="mg/L", "Cliver"="mg/L", "Cven"="mg/L", "Clung"="mg/L", "Cart"="mg/L", "Crest"="mg/L", "Ckidney"="mg/L", "Cplasma"="mg/L", "Cvehicle"="mg/L", "Cvehicle_infinite"="mg/L")

df.simulation = concentration_data_meade2023



# Find experiments with unique dosing regimes/species/etc and average those that had multiple time points
# for the same type of experiment. Rename each experiment with new_id.
keep.cols=c("PREFERRED_NAME","test_substance_dtxsid", "conc_medium_obs", "conc_medium_normalized", "administration_route_normalized", "species", "BW", "dose_level_new", "dose_level_new_units","totalSA","Vvehicle","InfiniteDose","dose_duration_num2","dose_duration_units2","Fskin_exposed","Exhalation","dose_vehicle_normalized","Water.Solubility.mg.L", "Vapor.Pressure.mmHg","Boiling.Point.C","logP", "Volatility", "Kvehicle2water")
df.simulation_red=data.frame(df.simulation%>%group_by(across(all_of(c("time_hr",keep.cols) )))%>%
  summarize(conc_normalized=mean(conc_normalized),.groups="drop")%>%group_by(across(all_of(keep.cols)))%>% arrange(time_hr)%>% mutate(new_id=cur_group_id())%>%
  ungroup())

df.simulation_red2=df.simulation_red%>%mutate(across(c(dose_duration_num2,dose_duration_units2,Vvehicle),~ifelse(is.na(.x),list(NULL),.x)))

Simulations: Nested For Loops

We have nested for loops to account for each aspect of the data that needs to be simulated. Each ID number (or “new_id”) in df.simulation_red2 corresponds to a single experiment (or experiment with averaged concentration data if the subject species, dose, route, etc. was the same), where we have defined the necessary variables: species, chemical, dose, dose duration (when applicable), whether or not an Infinite dose is applicable, route of exposure, the collected medium/tissue being sampled (i.e., plasma, blood, liver) and subject parameters such as body mass (kg) and surface area of skin. Thus, we run through each unique ID and first parameterize the model using parameterize_dermal_pbtk() and then solve the model using solve_dermal_pbtk(). For each ID, we simulate each permeability method (Surrey and Potts-Guy); for each ID/permeability method, we simulate either a water or octanol vehicle; then for each combination, we simulate the model both with and without exhalation. Note that some experimental data indicate the vehicle used, but we still test both vehicles for simulation to compare. Likewise, we calculated volatility for each chemical (which indicate whether or not exhalation actually occurs), but we still generate for all options. Note that solve_dermal_pbtk() can handle IV or oral exposure by setting the “route” input variable. Each model output is saved to a list called conc_predictions_meade2023LIST. Then this list is converted to a data frame.

df.run.sims = df.simulation_red2
fkID = unique(df.run.sims$new_id)
method.perm = c("Potts-Guy", "UK-Surrey")# for
Kvehicle2water=c("octanol","water") # for
exhalation_for_simulation = c("TRUE","FALSE")


conc_predictions_meade2023LIST=list()


for (include.e in exhalation_for_simulation){
  for(i in 1:length(fkID)){
    for(mp in 1:length(method.perm)){
      for(k in 1:length(Kvehicle2water)){
        sub.df = subset(df.run.sims, new_id == fkID[i])
        specs = unique(sub.df[,c("test_substance_dtxsid","conc_medium_obs","administration_route_normalized",
                                 "species","BW","dose_level_new","dose_level_new_units",
                                 "totalSA","Vvehicle","InfiniteDose","dose_duration_num2","dose_duration_units2",
                                "Fskin_exposed","Exhalation")])
        # check for experiments with only 1 timepoint and reset sub.df
        if(length(sub.df$time_hr)==1){
            sub.df=rbind(sub.df, sub.df[rep(1,1),])
            sub.df[1,c("conc_normalized","time_hr")]=0
        }
      
        parms = parameterize_dermal_pbtk(dtxsid = specs$test_substance_dtxsid,
                                 species = specs$species,
                                 totalSA = specs$totalSA,
                                 default.to.human=T,
                                 Kvehicle2water = Kvehicle2water[k],
                                 InfiniteDose = specs$InfiniteDose,
                                 method.permeability = method.perm[mp],
                                 
                                 model.type="dermal_1subcomp",
                                 clint.pvalue.threshold = Inf,
                                 suppress.messages = T
                                 )
        parms$BW = specs$BW
        parms$Fskin_exposed = specs$Fskin_exposed
      
        # Remove exhalation or keep exhalation (except for collected air which is )
        if(include.e == FALSE) {parms$Qalvc = 0}
        
        if(unique(sub.df$PREFERRED_NAME!="DCM"))
          { # tolerance is bad for DCM, so set tolerance to default values for all other compounds.
           rtol=1e-5; atol=1e-5
        }else{
           rtol=1e-4; atol=1e-4
        }
        
    
        sub.df=sub.df[order(sub.df$time_hr),]
        timey = unique(c(0,sub.df$time_hr)) #time in hours
        gc() # garbage collection (memory)
        ok=data.frame(solve_dermal_pbtk(parameters = parms,
                        times = timey/24, #convert to days
                         route = tolower(specs$administration_route_normalized),
                         model.type="dermal_1subcomp",
                         plot = FALSE,
                         output.units = output.units.original,
                         washoff = TRUE,
                         input.units = specs$dose_level_new_units,
                         initial.dose=specs$dose_level_new,
                         dose.duration = unlist(specs$dose_duration_num2),
                         dose.duration.units = unlist(specs$dose_duration_units2),
                         Vvehicle = unlist(specs$Vvehicle),
                        #dosing.dermal = dosing.dermal,
                         Kvehicle2water = Kvehicle2water[k],
                         suppress.messages = TRUE,
                        atol = atol, rtol= rtol,
                        method="lsoda"
                         ))
         # extract only times we need
          ok.time = subset(ok, time%in%c(ok$time[!(is.na(match(ok$time, sub.df$time_hr/24)))]))
          sub.df$Prediction = ok.time[,c(specs$conc_medium_obs)]
          sub.df$exhalation_for_simulation = include.e
          
          # specify vehicle and perm. method and vehicle if dermal
          if(specs$administration_route_normalized != "Dermal"){
            sub.df$method.perm = NA
            sub.df$Kvehicle2water = NA
            sub.df$method.name = specs$administration_route_normalized
          }else{
            sub.df$method.perm = method.perm[mp]
            sub.df$Kvehicle2water = Kvehicle2water[k]
            sub.df$method.name = paste(method.perm[mp], "- vehicle = ", Kvehicle2water[k])
          }

        conc_predictions_meade2023LIST[[paste0(fkID[i],".",Kvehicle2water[k],".",method.perm[mp],".exh",include.e)]]= sub.df
      }
    }
  }
}
#Bind data together to create a dataframe and organize some columns
 conc_predictions_meade2023DF = do.call(rbind, conc_predictions_meade2023LIST) 
 rownames(conc_predictions_meade2023DF)=NULL



 cgwtools::resave(conc_predictions_meade2023DF, conc_predictions_meade2023LIST, file=paste0("meade2023_",format(Sys.time(), "%b_%d_%Y"),".Rdata"))

Now, if run.simulations = FALSE, we must load data. This load supptab1_meade2023, concentration_data_meade2023, conc_predictions_meade2023DF, conc_predictions_meade2023LIST.

load(paste0(loc.wd, "/meade2023_Jun_13_2025.Rdata"))

Step Four: Calculate Root Mean Squared Log10 Errors (RMSLE)

Now, we want to get the RMSLE for each chemical, comparing the model to the data. Let us group by chemical, method, and whether or not we include exhalation.

conc_predictions_meade2023DF$method.perm = gsub( "UK-Surrey","Surrey", conc_predictions_meade2023DF$method.perm)
conc_predictions_meade2023DF$method.name = gsub( "UK-Surrey","Surrey", conc_predictions_meade2023DF$method.name)
df.sim = conc_predictions_meade2023DF
# Convert aexhaled air from mg to mg/l (assume a room of x liters, for humans - humans 
# the only species in these exhaled experiments)
df.sim[which(df.sim$conc_medium_obs=="Aexhaled"),"Prediction"]=df.sim[which(df.sim$conc_medium_obs=="Aexhaled"),"Prediction"]/10000

df.sim$conc_normalized[which(df.sim$conc_normalized==0)]=1e-10
df.sim$Prediction[which(df.sim$Prediction==0)]=1e-10
df.sim.rmsle=df.sim %>% group_by(PREFERRED_NAME, method.name, exhalation_for_simulation)%>%
  mutate(cmax_obs = max(conc_normalized), cmax_pred=max(Prediction)) %>%
  mutate(RMSLE = rmse(log10(conc_normalized),log10(Prediction)),
    Cmax_RMSLE = rmse(log10(cmax_obs), log10(cmax_pred)))


# Linear fit
ddf.sim.rmsle <- df.sim.rmsle %>% group_by(method.name) %>% do(
  {
    Adjusted.R2 = summary(lm(log10(Prediction) ~ log10(conc_normalized),data=.))$adj.r.squared
    Coefficient = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$coefficients[2]
    Intercept = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$coefficients[1]
    Residuals = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$residuals
    Fitted.Values = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$fitted.values
    data.frame(.,Adjusted.R2,Coefficient,Intercept,Residuals,Fitted.Values)
  }
)

Step Five: Plot Simulations

Let us compare the RMSLE and the volatility (i.e., boiling point) for each compound (assuming exhalation).

# Plot RMSLE vs Boiling Point

data = subset(df.sim.rmsle, exhalation_for_simulation==TRUE)

data=data%>%group_by(PREFERRED_NAME, Boiling.Point.C,method.name)%>%
  summarize(Chemical.cmax.RMSLE = rmse(log10(cmax_obs),log10(cmax_pred)),
    Chemical.RMSLE=rmse(log10(conc_normalized),log10(Prediction)))

data.rect = data.frame(xmin=c(0,75,250,400),
                       xmax=c(75,250,400,Inf),
                       ymin=rep(-Inf,4),
                       ymax=rep(Inf,4),
                       Volatility=c("Very Volatile",
                                    "Volatile",
                                    "Semi-Volatile",
                                    "Not Volatile"))
data.rect$Volatility <- factor(data.rect$Volatility, levels=data.rect$Volatility)
viridis.colors <- viridis(4)
plot <- ggplot(data=data,aes(x=`Boiling.Point.C`,y=Chemical.RMSLE)) +
  geom_point()+
  geom_smooth(method="lm",formula='y~x') +
  geom_text_repel(data=data,aes(x=`Boiling.Point.C`,y=Chemical.RMSLE,label=PREFERRED_NAME),
                  max.overlaps=Inf)+
  facet_wrap(~method.name,scales="free") +
  #scale_x_log10() +
  geom_rect(data=data.rect,
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,fill=Volatility),
            alpha=0.5,
            inherit.aes=FALSE) +
  scale_fill_manual(values=c(`Not Volatile`=viridis.colors[1],
                             `Semi-Volatile`=viridis.colors[2],
                             `Volatile`=viridis.colors[3],
                             `Very Volatile`=viridis.colors[4]))+

  labs(y="RMSLE",title="Boiling Point vs. RMSLE",x="Boiling Point (C)")+
  theme_bw() +
  theme(text=element_text(size=33),
        plot.title=element_text(hjust=0.5),
        legend.position="bottom")
plot
ggsave(paste(getwd(), "/Figures_May1/BPvsRMSLE_suppl.png",sep=""),
       width = 20, height = 10, units = "in")

plot2 <- ggplot(data=data,aes(x=`Boiling.Point.C`,y=Chemical.cmax.RMSLE)) +
  geom_point()+
  geom_smooth(method="lm",formula='y~x') +
  geom_text_repel(data=data,aes(x=`Boiling.Point.C`,y=Chemical.cmax.RMSLE,label=PREFERRED_NAME),
                  max.overlaps=Inf)+
  facet_wrap(~method.name,scales="free") +
  #scale_x_log10() +
  geom_rect(data=data.rect,
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,fill=Volatility),
            alpha=0.5,
            inherit.aes=FALSE) +
  scale_fill_manual(values=c(`Not Volatile`=viridis.colors[1],
                             `Semi-Volatile`=viridis.colors[2],
                             `Volatile`=viridis.colors[3],
                             `Very Volatile`=viridis.colors[4]))+

  labs(y="RMSLE",title="Boiling Point vs. Cmax RMSLE",x="Boiling Point (C)")+
  theme_bw() +
  theme(text=element_text(size=33),
        plot.title=element_text(hjust=0.5),
        legend.position="bottom")
#annotate_figure(plot2,top=text_grob("Boiling Point vs RMSLE",size=40,face="bold"))
ggsave(paste(getwd(), "/Figures_May1/BPvsRMSLE_cmax_suppl.png",sep=""),
       width = 20, height = 10, units = "in")

Now, we compare simulations with exhalation to those without exhalation (as a function of boiling point).

df.exh.compare=df.sim.rmsle

Exh = data.frame(df.exh.compare %>% group_by(PREFERRED_NAME,method.name, Boiling.Point.C)%>%
  arrange(exhalation_for_simulation) %>%
  summarize(Cmax=Cmax_RMSLE[exhalation_for_simulation==FALSE]-Cmax_RMSLE[exhalation_for_simulation==TRUE], Full_time_course = RMSLE[exhalation_for_simulation==FALSE]-RMSLE[exhalation_for_simulation==TRUE]))

Exh$method.name= sub(" - ", " -\n", Exh$method.name)



#stack data
stackExh = data.frame(Exh %>%
  pivot_longer(cols = c("Cmax", "Full_time_course"), 
               names_to = "ind", 
               values_to = "values"))

data.text = unique(stackExh[,c("Boiling.Point.C","values","PREFERRED_NAME","ind","method.name")])
data.rect = data.frame(xmin=c(0,75,250,400),
                       xmax=c(75,250,400,Inf),
                       ymin=rep(-Inf,4),
                       ymax=rep(Inf,4),
                       Volatility=c("Very Volatile",
                                    "Volatile",
                                    "Semi-Volatile",
                                    "Not Volatile"))
data.rect$Volatility <- factor(data.rect$Volatility, 
                               levels=data.rect$Volatility)
viridis.colors <- viridis(4)

# rename labels for facets
# labeller=labeller(ind = c(Cmax="Cmax",
#               Full_time_course = "Full time-course"))

plot <- ggplot(data=stackExh,aes(x=Boiling.Point.C,y=values)) +
  geom_point()+ #geom_line()+
  theme_bw()+
  geom_smooth(method='lm') +
  facet_grid(ind~method.name, scales="free", labeller=labeller(ind = c(Cmax="Cmax",
              Full_time_course = "Full time-course")))+
  geom_rect(data=data.rect,
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,fill=Volatility),
            alpha=0.5,
            inherit.aes=FALSE) +
  scale_fill_manual(values=c(`Not Volatile`=viridis.colors[1],
                             `Semi-Volatile`=viridis.colors[2],
                             `Volatile`=viridis.colors[3],
                             `Very Volatile`=viridis.colors[4]))+
  geom_text_repel(data=data.text,aes(x=Boiling.Point.C,y=values,label=PREFERRED_NAME),
                  max.overlaps=Inf)+
  #facet_grid(cols=vars(Route.or.Method.if.Dermal), rows=vars(RMSLE.Type),scales="free_y")+
  #scale_x_log10() +
  geom_hline(aes(yintercept=0)) +
  labs(y=expression(Delta~"RMSLE"),x="Boiling Point (C)",title="Change in RMSLE when Exhalation is Added")+
  scale_color_brewer(type="qual",palette="Paired") +
  theme(axis.text.x = element_text(angle=45,vjust=1,hjust=1),
        plot.title=element_text(hjust=0.5),
        text=element_text(size=30),
        legend.position="bottom")
plot

ggsave(paste(getwd(), "/Figures_May1/BP_RMSLE_Exhalation_main.png",sep=""),
       width = 20, height = 10, units = "in")

# P-values for regression lines
stackExhStats = stackExh%>% group_by(ind, method.name) %>%
  summarize(model = list(lm(values~Boiling.Point.C, data=cur_data())),
            .groups="drop")%>%
  rowwise()%>%
  mutate(msummary=list(summary(model)),
         rsquared=msummary$r.squared,
         p_val= pf(msummary$fstatistic[1], msummary$fstatistic[2], msummary$fstatistic[3], lower.tail=FALSE)
        )%>%select(ind,method.name,rsquared,p_val)

Now, let’s generate CvT plots for each chemical and color by whether or not exhalation was included.

df.sim2 = df.sim.rmsle

chem = unique(df.sim2$PREFERRED_NAME)

for(chx in 1:length(chem)){
  each.chem = subset(df.sim2, PREFERRED_NAME==chem[chx])
  #each.chem$dose = paste(each.chem$dose_level_new, each.chem$dose_level_new_units)
  each.chem$dose=paste(each.chem$dose_level_new, each.chem$dose_level_new_units)
  # each.chem$dose = factor(each.chem$dose,
  #                         levels=unique(each.chem$dose)[order(unique(each.chem$dose_level_new))])
  
  each.chem$method.name= sub(" - ", " -\n", each.chem$method.name)
  each.chem$subject=paste(each.chem$species, "BW",each.chem$BW, "kg")
  
  
 ggplot(each.chem)+
    geom_point(data=each.chem, aes(x=time_hr, y=conc_normalized,colour=dose))+
    geom_line(data=each.chem, aes(x=time_hr, y=Prediction,colour=dose,linetype=exhalation_for_simulation),linewidth=1)+
    labs(title=paste(chem[chx], "\n Log-transformed Concentrations v time"), x=" Time (hr)", y="Log Concentrations (mg/L)", colour="")+
    facet_nested(~conc_medium_normalized+method.name,scales="free")+
   scale_y_log10() +
     theme(legend.position = "bottom")
    

  ggsave(paste(getwd(), "/Figures_May1/cvt_eachchem_by_exhalation_suppl", unique(each.chem$test_substance_dtxsid),".png",sep=""),
         width = 8, height = 6, units = "in")
}

Now plot the data in Predicted vs Observed (PvO) plots for each of the 13 chemicals. First, we extract only simulated results that match the correct exhalation conditions by finding simulations where exhalation_for_simulations = T if Exhalation (which was calculated from volatility parameters), is TRUE and FALSE if FALSE.

df.new = subset(df.sim.rmsle, exhalation_for_simulation==Exhalation)

chem = unique(df.new$PREFERRED_NAME)

for(chx in 1:length(chem)){
  each.chem = subset(df.new, PREFERRED_NAME==chem[chx])
  #each.chem$dose = paste(each.chem$dose_level_new, each.chem$dose_level_new_units)
  each.chem$dose=paste(each.chem$dose_level_new, each.chem$dose_level_new_units)
  each.chem$dose = factor(each.chem$dose,
                          levels=unique(each.chem$dose)[order(unique(each.chem$dose_level_new))])
  
  each.chem$method.name= sub(" - ", " -\n", each.chem$method.name)
  each.chem$subject=paste(each.chem$species, "BW",each.chem$BW, "kg")
  
  ggplot(each.chem, aes(x=conc_normalized,y=Prediction))+
    geom_point(aes(colour=dose,shape=conc_medium_normalized))+
    labs(title=paste(chem[chx], "\n Log-transformed Concentrations"), x="Observed Concentrations (mg/L)", x="Predicted Concentrations (mg/L)", colour="", shape="")+
    #facet_wrap(~dose,scales="free")+
    facet_wrap(~method.name,scales="free")+
    geom_smooth(method='lm',formula='y~x',colour="black")+
     geom_abline(linetype=2)+
     #geom_smooth(aes(color=method.name),method='lm',formula='y~x')+
          scale_x_log10() +
    scale_y_log10() +
    theme_bw() + theme(text=element_text(size=11),
                       strip.text=element_text(size=13))
    #theme(legend.position = "bottom")

  ggsave(paste(getwd(), "/Figures_May1/pvo_bydose_suppl", unique(each.chem$test_substance_dtxsid),".png",sep=""),
         width = 6, height = 6, units = "in")
  
  
}


# Linear fit for these.
linear.fit.new <- df.new %>% group_by(method.name) %>% do(
  {
    Adjusted.R2 = summary(lm(log10(Prediction) ~ log10(conc_normalized),data=.))$adj.r.squared
    Coefficient = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$coefficients[2]
    Intercept = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$coefficients[1]
    Residuals = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$residuals
    Fitted.Values = lm(log10(Prediction) ~ log10(conc_normalized),data=.)$fitted.values
    data.frame(.,Adjusted.R2,Coefficient,Intercept,Residuals,Fitted.Values)
  }
)
# RMSLE grouped by method type
rmsle.all = df.new%>%group_by(method.name)%>%
  summarize(  Adjusted.R2 = summary(lm(log10(Prediction) ~ log10(conc_normalized)))$adj.r.squared,
              RMSLE_total =rmse(log10(Prediction),log10(conc_normalized)),
                CMAX_RMSLE_total =rmse(log10(cmax_pred) ,log10(cmax_obs)    )
              
              )

now plot log-transformed predicted vs. observed values and group by permeability method, vehicle (or route of exposure), and tested tissue.

df.new = subset(df.sim, exhalation_for_simulation==Exhalation)
df.sim.rmsle2=df.new %>% group_by(PREFERRED_NAME, method.name)%>%
  mutate(cmax_obs = max(conc_normalized), cmax_pred=max(Prediction)) %>%
  mutate(RMSLE = rmse(log10(conc_normalized),log10(Prediction)),
    Cmax_RMSLE = rmse(log10(cmax_obs), log10(cmax_pred)))


col=c("firebrick", "indianred", "blue4","steelblue1","#2ca25f", "#addd8e")
plot.pvo=ggplot(df.sim.rmsle2, aes(x=conc_normalized, y=Prediction))+
  geom_point(aes( colour=method.name,shape=PREFERRED_NAME), alpha=1,size=2)+
    labs(title=paste0("(A) Concentrations"),
       y = "Simulated (mg/L)", x="Observed (mg/L)",
       color="Route or\n Method",shape="Chemical")+
    scale_shape_manual(values=1: length(unique(df.sim.rmsle2$PREFERRED_NAME)))+
    geom_smooth(aes(colour=method.name),method="lm",se=F)+
  #scale_color_brewer(type="qual",palette="Paired",direction=-1) +
  scale_colour_manual(values=col)+
  theme_bw() + 
  scale_x_log10() + scale_y_log10() +
  geom_abline(lty="dashed") +
  #coord_fixed(xlim=lims,ylim=lims) +
  theme(plot.title=element_text(hjust=0.5),
        text=element_text(size=30))

plot.cmax = ggplot(df.sim.rmsle2, aes(x=cmax_obs, y=cmax_pred))+
  geom_point(aes( colour=method.name,shape=PREFERRED_NAME), alpha=1,size=2)+
    labs(title=paste0("(B) Max Concentrations"),
       y = "Simulated (mg/L)", x="Observed (mg/L)",
       color="Route or\n Method", shape="Chemical")+
  #scale_color_brewer(type="qual",palette="Paired",direction=-1) +
   scale_colour_manual(values=col)+
     scale_shape_manual(values=1: length(unique(df.sim.rmsle2$PREFERRED_NAME)))+
  theme_bw() + 
  scale_x_log10() + scale_y_log10() +
  geom_abline(lty="dashed") +
  #coord_fixed(xlim=lims,ylim=lims) +
  theme(plot.title=element_text(hjust=0.5),
        text=element_text(size=30))

ggarrange(plot.pvo, plot.cmax, common.legend=T,legend="right")
ggsave(paste(getwd(), "/Figures_May1/PvO_main.png",sep=""),
       width = 20, height = 10, units = "in",dpi=300)
  
 
method.rmsle.only=df.new %>% group_by(PREFERRED_NAME, method.name)%>%
    mutate(cmax_obs = max(conc_normalized), cmax_pred=max(Prediction)) %>%
    summarize(RMSLE = rmse(log10(conc_normalized),log10(Prediction)),
           Cmax_RMSLE = rmse(log10(cmax_obs), log10(cmax_pred))) %>%group_by(method.name)%>%summarize(mean(Cmax_RMSLE),
                                                                                                      mean(RMSLE))

Now let’s plot a heat map of the RMSLE and Cmax RMSLE. We will subset the dataframe to only plot the simulations appropriate for exhalation (that is, for a volatile chemical where Exhalation = TRUE, we will extract the solutions when exhalation_for_simulation = TRUE).We will also make note of the experiments for which a dosing vehicle was specified (water or octanol). For dermal, we will plot the results using both permeability methods and both vehicle methods, but will mark boxes over the cases where the vehicle matched to the simulated vehicle. We also want to plot the chemicals in order of volatility (based on boiling point).

# Extract correct exhalation

# df = conc_predictions_meade2023DF
# df$conc_normalized[which(df$conc_normalized==0)]=1e-10
# df$Prediction[which(df$Prediction==0)]=1e-10
df=df.sim
dfmatchE =  subset(df, exhalation_for_simulation==Exhalation)
# order chemicals by volatility
vol.chem = unique(dfmatchE[,c("PREFERRED_NAME","Boiling.Point.C")])
vol.chem2=vol.chem[order(vol.chem$Boiling.Point.C,decreasing=F),]

#also want to extract table to identify  experiments that actually did list the vehicle as either water or octanol (or similar)
dfmatchE2 = data.frame(dfmatchE%>%group_by(PREFERRED_NAME,method.name,dose_vehicle_normalized,Kvehicle2water,administration_route_normalized)%>%
      mutate(cmax_obs = max(conc_normalized), cmax_pred=max(Prediction)) %>%
  summarize(Cmax_RMSLE = rmse(log10(cmax_obs),log10(cmax_pred)),
         RMSLE = rmse(log10(conc_normalized),log10(Prediction))))

# Reorganize the columns so RMSLE values can be grouped by Cmax or total RMSLE and then take average RMSLE values of oral or iv values so we have one per chemical
dfmatchE3 = data.frame(dfmatchE2 %>%
  pivot_longer(cols = c("Cmax_RMSLE", "RMSLE"), 
               names_to = "ind", 
               values_to = "value") %>% group_by(PREFERRED_NAME,ind,administration_route_normalized) %>%
    mutate(value = if_else(administration_route_normalized !="Dermal", mean(value[administration_route_normalized !="Dermal"]), value)
           ))

# Make PREFERRED_NAME a factor and set the levels to match those from the volatility order.
dfmatchE3$PREFERRED_NAME = factor(dfmatchE3$PREFERRED_NAME, 
                                              levels=vol.chem2$PREFERRED_NAME)

Generate the plot.

ggplot(dfmatchE3, aes(y=PREFERRED_NAME, x=method.name))+
  geom_tile(aes(fill=value))+
  #geom_tile(aes(fill=Chemical.Cmax.RMSLE.NE))+
  facet_wrap(~ind)+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  scale_fill_viridis()+
  geom_text(aes(label=round(value,2)),colour="white",size=4)+
  labs(title="",
       x="Route or\n Method",
       fill="RMSLE",y="Chemical")+
  theme_bw() +
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1),
        axis.text.y=element_text(angle=45,vjust=1.05,hjust=0.9),
        plot.title=element_text(hjust=0.5),
        text=element_text(size=20))+
  geom_tile(data=subset(dfmatchE3, dose_vehicle_normalized==Kvehicle2water & administration_route_normalized=="Dermal"),aes(y=PREFERRED_NAME, x=method.name),
            size=1.5,fill=NA,colour="red" )
ggsave(paste(getwd(),"/Figures_May1/RMSLE_compare_main.png",sep=""),width=11,height=8,units="in")

Results 2: IVIVE

Now, we perform in vitro to in vivo extrapolation using the dermal model to obtain some occupational hazardous levels.

Step One: Load Data

First we load ToxCast Data data from the EPA-ORD-CCTE Clowder. Some of the variables in this dataset of importance include:

# Load in vitro data
load("invitrodb_3_5_mc5.Rdata")

# Get chemicals of interest 
load_sipes2017(overwrite=F)

Step Two: Load Chemical Data

# Using list, download chem info from CompTox Dashboard for chem.list chemicals, and load.
df.chem <- read.xlsx(file="CCD-Batch-Search_2023-03-08_12_23_25.xlsx",
                     sheetName="Main Data")# %>%filter(!is.na(BOILING_POINT_DEGC_OPERA_PRED))
df.chem <- type_convert(df.chem,na="N/A")
df.chem <- type_convert(df.chem,na="N/A")

Step Three: Clean Data (Only Non-Volatiles and a few of the Semi-Volatiles from the validation study)

#Get list of chemicals that can run in dermal model:
chem.list <- get_cheminfo(info="DTXSID",species="Human",model="dermal_1subcomp")
write.table(chem.list,
            file=paste0(loc.wd,"/tables/ChemicalList.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")

# Take out volatile chemicals
chems.from.invivo <- c("DTXSID0020232","DTXSID6022923","DTXSID8024151","DTXSID3022409","DTXSID1026081")
df.chem <- df.chem %>% filter(!is.na(BOILING_POINT_DEGC_OPERA_PRED))%>% 
  filter((BOILING_POINT_DEGC_OPERA_PRED >= 400) | (DTXSID %in% chems.from.invivo))
chem.list <- df.chem$DTXSID

#get sample to test code with
set.seed(100)
#chem.list <- sample(chem.list,50)
toxcast.dermal <- mc5 %>% subset(dsstox_substance_id %in% chem.list)

Step Four: Compute Q10.AC50 (and other parameters of interest)

# CONCENTRATIONS OF INTEREST - https://doi.org/10.1093/bioinformatics/btw680 ----------------
toxcast.table <- NULL
old.time <- Sys.time()
num.chem <- length(chem.list)
for(k.chem in 1:num.chem){ #for each chemical
  this.id <- chem.list[k.chem]
  # Modulus operator
  if(k.chem %% 100==0){
    new.time <- difftime(Sys.time(), old.time,units="secs")
    cat(paste0("Run ", k.chem,":",round(new.time)," secs \n"))
  }
  toxcast.dermal.chem <- subset(toxcast.dermal, dsstox_substance_id == this.id) #subset over chemical
  toxcast.dermal.chem.hits <- subset(toxcast.dermal.chem, hitc==1) #only look at hits
  if(dim(toxcast.dermal.chem.hits)[1]>0){ #if there were any hits
    this.row <- data.frame(Chemical = as.character(toxcast.dermal.chem.hits[1,"chnm"]),
                           DTXSID = this.id,
                           Total.Assays = dim(toxcast.dermal.chem)[1],
                           Unique.Assays = length(unique(toxcast.dermal.chem$aeid)),
                           Total.Hits = dim(toxcast.dermal.chem.hits)[1],
                           Unique.Hits = length(unique(toxcast.dermal.chem.hits$aeid)),
                           Low.AC50 = min(toxcast.dermal.chem.hits$modl_ga),
                           Low.AC10 = min(toxcast.dermal.chem.hits$modl_ac10),
                           Low.ACC = min(toxcast.dermal.chem.hits$modl_acc),
                           Q10.AC50 = quantile(toxcast.dermal.chem.hits$modl_ga,probs=0.1),
                           Q10.AC10 = quantile(toxcast.dermal.chem.hits$modl_ac10,probs=0.1),
                           Q10.ACC = quantile(toxcast.dermal.chem.hits$modl_acc,probs=0.1),
                           Med.AC50 = quantile(toxcast.dermal.chem.hits$modl_ga,probs=0.5),
                           Med.AC10 = quantile(toxcast.dermal.chem.hits$modl_ac10,probs=0.5),
                           Med.ACC = quantile(toxcast.dermal.chem.hits$modl_acc,probs=0.5),
                           Q90.AC50 = quantile(toxcast.dermal.chem.hits$modl_ga,probs=0.9),
                           Q90.AC10 = quantile(toxcast.dermal.chem.hits$modl_ac10,probs=0.9),
                           Q90.ACC = quantile(toxcast.dermal.chem.hits$modl_acc,probs=0.9)
                           )
    toxcast.table <- rbind(toxcast.table, this.row)
  }
}
rownames(toxcast.table) <- seq(1,dim(toxcast.table)[1]) # set rownames to be sequential numbers

#View table
knitr::kable(head(toxcast.table[,1:6]), caption = "Summarized ToxCast Data",
             floating.environment="sidewaystable")
write.table(toxcast.table,
            file=paste0(loc.wd,"/tables/ToxCastSummaryTable.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")

Step Five: Conduct In Vitro to In Vivo Extrapolation


# Get Cmax = max(Cplasma) from dermal model for each chemical in toxcast table
old.time <- Sys.time()
include.these=which(toxcast.table$DTXSID %in% get_cheminfo(info="dtxsid", suppress.messages=T))
toxcast.table2 = subset(toxcast.table,DTXSID %in% toxcast.table$DTXSID[include.these])

atol <- rtol <- 1e-5

num.chem <- length(toxcast.table2$DTXSID)
plot.each=F
# Set Dermal Solutions
method.ls <- c("Potts-Guy","UK-Surrey")
toxcast.ls=list()
for(k.chem in 1:num.chem){
  this.id = toxcast.table2$DTXSID[k.chem]
  suppress.messages = TRUE

    for(k.method in 1:2){
      this.method <- method.ls[k.method]; 
    #  if (this.method=="2-Compartment") this.method.input <- "NULL"
      this.model.type <-"dermal_1subcomp"
      end.time=0
      parms=parameterize_dermal_pbtk(dtxsid=this.id, 
                                     model.type = this.model.type,
                                     method.permeability = this.method,
                                     clint.pvalue.threshold=Inf,
                                     suppress.messages=TRUE,
                                     Kvehicle2water = 1,
                                     species="Human",
                                     default.to.human = T)
      parms$InfiniteDose=1
      #if(end.time<2){ #make sure solver finishes
        out=try( solve_dermal_pbtk(parameters=parms,
                                     model.type=this.model.type,
                                     method.permeability=this.method,
                                     #Kvehicle2water=1, #vehicle is water
                                     days=5,
                                     initial.dose = 1,
                                     input.units = "ppm",
                                     dose.duration=8,
                                     dose.duration.units="hr",
                                     washoff=T,
                                    # InfiniteDose=T,
                                     atol = atol, rtol= rtol,
                                     method="lsoda",
                                     suppress.messages=suppress.messages) )#put SA for hands!!!
         if (is(out,"try-error")){
           Cmax=0
           browser()
         }
       
        else{end.time <- out[nrow(out),"time"]
        

      Cmax <- max(out[,"Cplasma"])
   
        }
       toxcast.ls[[paste0(k.chem,".",this.method)]]=cbind(toxcast.table2[toxcast.table2$DTXSID==this.id,],Cmax=Cmax,Permeability = this.method)
    }
  
}

toxcast2.all=do.call(rbind,toxcast.ls)

#     # Calculate the EquivDose in units same as input.units (ppm)
toxcast2.all$EquivDose = signif(10^toxcast2.all$Q10.AC50 /
                                  toxcast2.all$Cmax,
                                4)
# Edit data
toxcast3=toxcast2.all%>%
  rename(`Permeability Type` = Permeability,
         `Equivalent Dose`=EquivDose)%>%
  mutate(Chemical = ifelse(Chemical=="Methyl tert-butyl ether","MTBE",Chemical),
         Chemical = ifelse(Chemical=="3,3',5,5'-Tetrabromobisphenol A","Bromdian",Chemical),
         Chemical = ifelse(Chemical=="4,4'-Sulfonyldiphenol","Bisphenol S",Chemical),
         Chemical = ifelse(Chemical=="Dichloromethane","DCM",Chemical),
         Chemical = ifelse(Chemical=="Tetrachloroethylene","PERC/TCE",Chemical),
         Chemical = ifelse(Chemical=="1-Methylbenzene","Toluene",Chemical))

# Remove unused columns
ivive_meade2023 <- toxcast3 %>% select(Chemical,DTXSID,`Permeability Type`,`Equivalent Dose`,Cmax)
chem.props <- subset(chem.physical_and_invitro.data,
                     DTXSID %in% ivive_meade2023$DTXSID)[
                       ,c("DTXSID","MW","logWSol")]


ivive_meade2023 <- merge(ivive_meade2023,
                         all.x = TRUE,
                         chem.props,
                         by="DTXSID")
ivive_meade2023 <- within(ivive_meade2023, 
                          WS.ppm <- 
signif(MW*10^logWSol*1.001142303,4))

# Save data
 cgwtools::resave(ivive_meade2023,
                  file=paste0(
                    loc.wd,
                    "/meade2023_",
                    format(Sys.time(), "%b_%d_%Y"),
                    ".Rdata"))

Step Six: Plot Results

# Modify for plotting
ivive_meade2023$`Permeability Type` = gsub("UK-Surrey","Surrey",ivive_meade2023$`Permeability Type`)
ivive_meade2023[ivive_meade2023==Inf] = 1e15
write.table(subset(ivive_meade2023,
                   ivive_meade2023[,"Equivalent Dose"] < 
                     ivive_meade2023[,"WS.ppm"]),
            file=paste0(loc.wd,"/tables/IVIVEAchievable.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")       
ivive_meade2023$`Equivalent Dose` <- log10(ivive_meade2023$`Equivalent Dose`)
ivive_meade2023$`Permeability Type` = paste0("Permeability Method is ",ivive_meade2023$`Permeability Type`)
write.table(ivive_meade2023,
            file=paste0(loc.wd,"/tables/IVIVEResults.txt"),
            col.names=TRUE,
            row.names=FALSE,
            quote=FALSE,
            sep="\t")
# Plot Histogram
data = ivive_meade2023 %>% group_by(`Permeability Type`) %>% mutate(Mean = mean(`Equivalent Dose`,na.rm=TRUE)) %>% mutate(Median = median(`Equivalent Dose`,na.rm=TRUE))
AED.labels <- c("Needs Protection (AED caused by <1ppm)",
                "AED Possibly Acheivable",
                "\"Safe\" (Unachievable AED)",
                "No Dermal Penetration")

# Shaded rectangles with cutoff at 10,000 ppm:
data.rect.1e4ppm = data.frame(xmin=c(-Inf,0,4,12),
                       xmax=c(0,4,12,Inf),
                       ymin=rep(-Inf,4),
                       ymax=rep(Inf,4),
                       Bioactivity=AED.labels)
# Shaded rectangles with cutoff at 1,000,000 ppm:
data.rect.1e6ppm = data.frame(xmin=c(-Inf,0,6,12),
                       xmax=c(0,6,12,Inf),
                       ymin=rep(-Inf,4),
                       ymax=rep(Inf,4),
                       Bioactivity=AED.labels)
data.rect.1e4ppm$Bioactivity <- factor(data.rect.1e4ppm$Bioactivity, 
                                       levels=data.rect.1e4ppm$Bioactivity)
data.rect.1e6ppm$Bioactivity <- factor(data.rect.1e6ppm$Bioactivity, 
                                       levels=data.rect.1e6ppm$Bioactivity)

data.med = unique(data[,c("Permeability Type","Mean","Median")]) %>% pivot_longer(cols=c("Mean","Median"),names_to="Statistic",values_to="Statistic.Value")
viridis.colors <- viridisLite::viridis(4)
plot <- ggplot(data,aes(x=`Equivalent Dose`)) +
  geom_histogram(na.rm=TRUE,alpha=0.5,position="identity") +
  geom_freqpoly(na.rm=TRUE) +
  geom_vline(data = data.med, aes(xintercept=Statistic.Value,linetype=Statistic)) +
  geom_rect(data=data.rect.1e4ppm,
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,fill=Bioactivity),
            alpha=0.5,
            inherit.aes=FALSE) +
  #scale_x_log10() + 
  scale_fill_manual(values=setNames(viridis.colors[4:1], AED.labels))+
  labs(title=paste0("Distribution of concentrations causing AED\nfor hands submerged for 8 hours a day\nfor five days"),
       y = "Number of Chemicals", x="Log10 of Administrered Equivalent Doses (AEDs) ppm in water") +
  facet_col(vars(`Permeability Type`),scales="fixed") +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5))
plot  
ggsave(paste(getwd(), "/Figures_May1/histogram2.png",sep=""),
       width = 20, height = 12, units = "cm")

# WATER SOLUBILITY
ivive_meade2023 <- 
  ivive_meade2023 %>% 
  mutate(WSol = case_when( ((MW*10^logWSol)>=10) ~ "Water Soluble (WSol>10g/L)",
                           ((MW*10^logWSol)<10) ~ "Not Soluble"
                           ))
#ivive_meade2023ALL <- ivive_meade2023WS; ivive_meade2023ALL$WSol = "Non-Water Soluble"
ivive_meade2023WS <- ivive_meade2023 %>% filter(WSol=="Water Soluble (WSol>10g/L)")
ivive_meade2023NonWS <- ivive_meade2023 %>% filter(WSol=="Not Soluble")
# Different bands for soluble/non-soluble:
data.rect.1e4ppm$WSol <- "Not Soluble"
data.rect.1e6ppm$WSol <- "Water Soluble (WSol>10g/L)"
data.rect <- rbind(data.rect.1e4ppm,data.rect.1e6ppm)
data.rect[,"Permeability Type"] <- "Permeability Method is Potts-Guy"
data.rect2 <- data.rect
data.rect2[,"Permeability Type"] <- "Permeability Method is Surrey"
data.rect <- rbind(data.rect,data.rect2)

ivive_meade2023 <- rbind(ivive_meade2023NonWS,ivive_meade2023WS)
data = ivive_meade2023 %>% group_by(`Permeability Type`,`WSol`) %>% mutate(Mean = mean(`Equivalent Dose`,na.rm=TRUE)) %>% mutate(Median = median(`Equivalent Dose`,na.rm=TRUE))
data.med = unique(data[,c("Permeability Type","WSol","Mean","Median")]) %>% pivot_longer(cols=c("Mean","Median"),names_to="Statistic",values_to="Statistic.Value")
chem.to.label <-  c("DTXSID0020232",
                  # "DTXSID6022923",
                  # "DTXSID8024151",
                   "DTXSID3022409",
                   "DTXSID1026081", "DTXSID9048699",
                   "DTXSID0022858","DTXSID1037","DTXSID0023163",
                   "DTXSID5024845" )
data.text = data %>% filter(DTXSID %in% chem.to.label)
#PLOT WS
plot <- ggplot(data) +
  geom_rect(data=data.rect,
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,fill=Bioactivity),
            alpha=0.5#,
            #inherit.aes=FALSE
            ) +
  geom_histogram(na.rm=TRUE,alpha=0.5,position="identity",aes(x=`Equivalent Dose`)) +
  geom_freqpoly(na.rm=TRUE,aes(x=`Equivalent Dose`)) +
  geom_vline(data = data.med, aes(xintercept=Statistic.Value,linetype=Statistic)) +
#  geom_point(data=data.text,aes(x=`Equivalent Dose`,y=0),size=2) + 
#  geom_text_repel(data=data.text,aes(x=`Equivalent Dose`,y=0,label=Chemical),
#                  max.overlaps=14,
#                  hjust = "right",nudge_x=-1,nudge_y=0.2)+
  scale_fill_manual(values=setNames(viridis.colors[4:1], AED.labels))+
  #ylim(0,135)+
  #geom_text(x=) +
  labs(title=paste0("Distribution of concentrations causing AED\nfor hands submerged for 8 hours a day\nfor five days"),
       y = "Number of Chemicals", x="Log10 of Needed Concentration (ppm) in water") +
  facet_wrap(vars(`Permeability Type`,`WSol`),scales="free_y",ncol=2) +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5),
        text=element_text(size=20))
plot  
ggsave(paste(getwd(), "/Figures_May1/histogram_all_main.png",sep=""),
       width = 15, height = 10, units = "in")

Include a table with the equivalent doses, grouped by permeability method

data2=data[order(data$`Equivalent Dose`),]
data2$`Equivalent Dose`= exp(data2$`Equivalent Dose`)
data3=data.frame(data2%>%group_by(Chemical,DTXSID, `Permeability Type`,WSol,logWSol)%>%
  mutate(row=row_number())%>%
  ungroup() %>%
pivot_wider(id_cols=c(Chemical, DTXSID, row, WSol,logWSol),
  names_from = `Permeability Type`,
                    values_from = c(`Equivalent Dose`),
                    names_sep=""))
data4=data3[,-3]
data4=data4[order(data4$Permeability.Method.is.Surrey),]

Also, calculate permeability for each method for humans

pmat=list()
pm=c("UK-Surrey","Potts-Guy")
kw=c("water","octanol")
for(i in 1:length(supptab1_meade2023$DTXSID)){
  for(m in 1:2){ #each method
    for(k in 1:2){ # each vehicle
  p=parameterize_dermal_pbtk(dtxsid=supptab1_meade2023$DTXSID[i],
                             method.permeability =pm[m] ,species="Human",
                             Kvehicle2water = kw[k])
  pmat[[paste(i,".",m,".",k)]]=cbind(p$P, pm[m], kw[k],p$MW, supptab1_meade2023$DTXSID[i])
  
  }
  }
}

pmat2=do.call(rbind,pmat)
colnames(pmat2)=c("Perm", "Method","Vehicle","MW","DTXSID")
pmat2=data.frame(pmat2)
pmat2$Method=gsub("UK-Surrey","Surrey",pmat2$Method)
pmat2$Vehicle=paste("Vehicle is", pmat2$Vehicle)

pmat3=merge(pmat2,supptab1_meade2023, by="DTXSID")


ggplot(pmat3,aes(x=as.numeric(MW.x),y=as.numeric(Perm)))+
  geom_point()+
  facet_grid(Method~Vehicle,scales="free")+
  labs(x="MW",y="Kp - Permeability")
ggsave(paste(getwd(), "/Figures_May1/MWvPerm.png",sep=""),
       width = 4, height = 3.5, units = "in")

ggplot(pmat3,aes(x=as.numeric(logP),y=as.numeric(Perm)))+
  geom_point()+
  facet_grid(Method~Vehicle,scales="free")+
  labs(x="log(Kow)",y="Kp- Permeability")
ggsave(paste(getwd(), "/Figures_May1/logKowvPerm.png",sep=""),
       width = 4, height = 3.5, units = "in")
save.image(paste0(loc.wd,
                      "MeadeWorkspace",
                      Sys.Date(),
                      ".RData"))