Kapraun et al. (2022): Generic Human Gestational Model"

In httk: High-Throughput Toxicokinetics

Please send questions to [e-mail]

from "Evaluation of a Rapid, Generic Human Gestational Dose Model"

Dustin F.Kapraun, Mark Sfeir, Robert Pearce, Sarah Davidson, Annie Lumen, Andre Dallmann, Richard Judson, John F.Wambaugh

Reproductive Toxicology, 2022

https://doi.org/10.1016/j.reprotox.2022.09.004

Abstract

Chemical risk assessment considers potentially susceptible populations including pregnant women and developing fetuses. Humans encounter thousands of chemicals in their environments, few of which have been fully characterized in terms of internal human dosimetry and mode-of-action. Toxicokinetic (TK) information is needed to relate chemical exposure to potentially bioactive tissue concentrations. Observational data describing human gestational exposures are unavailable for most chemicals, but physiologically based TK (PBTK) models estimate such exposures. However, development of chemical-specific PBTK models themselves requires considerable time and resources. As an alternative, generic PBTK approaches describe a standardized physiology and characterize chemicals with a set of standard physical and TK descriptors -- primarily plasma protein binding and hepatic clearance. Here we report and evaluate a generic PBTK model of a human mother and developing fetus. We used a previously published set of formulas describing the major anatomical and physiological changes that occur during pregnancy to augment the High-Throughput Toxicokinetics (httk) software package. We performed simulations to estimate the ratio of concentrations in maternal and fetal plasma and compared these estimatesto literature in vivo measurements. We evaluated the model with literature in vivo time-course measurements of maternal plasma concentrations in pregnant and non-pregnant women. Finally, we demonstrated conceptual prioritization of chemicals measured in maternal serum based on predicted fetal brain concentrations. This new generic model can be used for TK simulations of any of the 856 chemicals with existing human-specific in vitro data that were found to be within the domain of the model, as well as any new chemicals for which such data become available. With sufficient evaluation, this gestational model may allow for in vitro to in vivo extrapolation of point of departure doses relevant to reproductive and developmental toxicity.

HTTK Version

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

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.

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.

# 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.

#
#
# Setup
#
#
#library(readxl)
library(ggplot2)
library(httk)
library(scales)
library(gridExtra)
library(cowplot)

Prepare custom functions

scientific_10 <- function(x) {                                  
  out <- gsub("1e", "10^", scientific_format()(x))              
  out <- gsub("\\+","",out)                                     
  out <- gsub("10\\^01","10",out)                               
  out <- parse(text=gsub("10\\^00","1",out))                    
}  

RMSE <- function(x)
{
  mean(x$residuals^2)^(1/2)
}

Number of chemicals

The number of chemicals with in vitro TK data ($Cl_{int}$ and $f_{up}$) that are also non-volatile and non-PFAS can be found using:

length(get_cheminfo(model="fetal_pbtk",suppress.messages=TRUE))

Data sets were curated from the literature to allow evaluation of the gestational PBTK model. In all cases chemical identities from the original publications were mapped onto DTXSID's from the CompTox Chemicals Dashboard (https://comptox.epa.gov/dashboard) [Williams et al., 2017] (https://doi.org/10.1186/s13321-017-0247-6). Statistical testing for correlation between predictions and observations was performed using R function "lm" and p-values were calculated according to an F-distribution.

Aylward cord-blood data

Aylward et al., 2014 compiled measurements of the ratio of maternal to fetal cord blood chemical concentrations at birth for a range of chemicals with environmental routes of exposure, including bromodiphenyl ethers, fluorinated compounds, organochlorine pesticides, polyaromatic hydrocarbons, tobacco smoke components, and vitamins. The PBTK model does not have an explicit cord blood compartment, so the ratio between maternal and fetal venous plasma concentrations was used as a surrogate when making comparisons of model predictions to these data.. For each chemical three daily oral doses (every eight hours) total 1 mg/kg/day were simulated starting from the 13th week of gestation until full term (40 weeks). Simulations were made both with and without the correction to $f_{up}^f$.

Load/format the data

#MFdata <- read_excel("Aylward-MatFet.xlsx")
MFdata <- httk::aylward2014

cat(paste("summarized data from over 100 studies covering ",
  length(unique(MFdata$DTXSID)[!(unique(MFdata$DTXSID)%in%c("","-"))]),
  " unique chemicals structures\n",sep=""))

# We don't want to exclude the volatiles and PFAS at this point:
MFdata.httk <- subset(MFdata,DTXSID %in% get_cheminfo(
  info="DTXSID",
  suppress.messages=TRUE))
MFdata.httk[MFdata.httk$Chemical.Category=="bromodiphenylethers",
  "Chemical.Category"] <- "Bromodiphenylethers"  
MFdata.httk[MFdata.httk$Chemical.Category=="organochlorine Pesticides",
  "Chemical.Category"] <- "Organochlorine Pesticides"  
  MFdata.httk[MFdata.httk$Chemical.Category=="polyaromatic hydrocarbons",
  "Chemical.Category"] <- "Polyaromatic Hydrocarbons"  

colnames(MFdata.httk)[colnames(MFdata.httk) == 
  "infant.maternal.conc...Central.tendency..calculate.j.k..or.report.paired.result."] <-
  "MFratio"
colnames(MFdata.httk)[colnames(MFdata.httk) == 
  "PREFERRED_NAME"] <-
  "Chemical"
colnames(MFdata.httk)[colnames(MFdata.httk) == 
  "details.on.matrix.comparison...e.g...cord.blood.lipid..maternal.serum.lipid..or.cord.blood.wet.weight..maternal.whole.blood.wet.weight"] <-
  "Matrix"

# Format the columns:
MFdata.httk$MFratio <- as.numeric(MFdata.httk$MFratio)
MFdata.httk$Chemical <- as.factor(MFdata.httk$Chemical)  
MFdata.httk$Matrix <- as.factor(MFdata.httk$Chemical)  
MFdata.httk$Chemical.Category <- as.factor(MFdata.httk$Chemical.Category)  

colnames(MFdata.httk)[15] <- "infant"
colnames(MFdata.httk)[16] <- "maternal"
colnames(MFdata.httk)[17] <- "obs.units"

MFdata.httk$infant <- suppressWarnings(as.numeric(MFdata.httk$infant))
MFdata.httk$maternal <- suppressWarnings(as.numeric(MFdata.httk$maternal))
MFdata.httk$AVERAGE_MASS <- as.numeric(MFdata.httk$AVERAGE_MASS)

Convert the units:

convert1.units <- c("ng/ml","ng/mL","ug/L","ug/l","ng/mL serum","ng/g",
  "ng/g wet wt.","ppb")

MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] <- 
  MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] / # ng/ml = ug/L
  MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"]  # ug/L -> uM
MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] <- 
  MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] / # ng/ml = ug/L
  MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"]  # ug/L -> uM
MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"obs.units"] <- "uM" 

convert2.units <- c("mg/L","ppm")

MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] <- 
  MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] * 1000 / # mg/L = ug/L
  MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"]  # ug/L -> uM
MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"] <- 
  MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"]* 1000 / # mg/L = ug/L
  MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"]  # ug/L -> uM
MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"obs.units"] <- "uM" 

Compare HTTK predictions with maternal:fetal observations from Aylward 2014

Note that we can't do an absolute scale comparison (for example, fetal:fetal or maternal:maternal) because we don't know the dose for the Aylward data but we assume that the maternal:fetal ratio is independent of dose and so we use the plasma ratio at full term (40 weeks) resulting from a 1 mg/kg/day exposure rate starting in week 13.

# Simulate starting from the 13th week of gestation until full term (40 weeks):
times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))

# For each chemical with maternal-fetal ratio data and HTTK in vitro data:  
for (this.id in unique(MFdata.httk$DTXSID))
{
  print(this.id)
  p <- parameterize_steadystate(dtxsid=this.id,
                        suppress.messages=TRUE)
  # skip chemicals where Fup is 0:  
  if (p$Funbound.plasma>1e-4)
  {

    # Load the chemical-specifc paramaters:
    p <- parameterize_fetal_pbtk(dtxsid=this.id,
                                                  fetal_fup_adjustment =TRUE,
                                                  suppress.messages=TRUE)
    # Skip chemicals where the 95% credible interval for Fup spans more than 0.1 to
    # 0.9 (that is, Fup is effectively unknown):
    if (!is.na(p$Funbound.plasma.dist))
    {
      if (as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][3])>0.9 & 
          as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][2])<0.11)
      {
        skip <- TRUE
      } else skip <- FALSE
    } else skip <- FALSE

    if (!skip)
    {
      # Solve the PBTK equations for the full simulation time assuming 1 total daily
      # dose of 1 mg/kg/day spread out over three evenly spaces exposures:
      out <- solve_fetal_pbtk(
        parameters=p,
        dose=0,
        times=times,
        daily.dose=1,
        doses.per.day=3,
        output.units = "uM",
        suppress.messages=TRUE)
      # Identify the concentrations from the final (279th) day:
      last.row <- which(out[,"time"]>279)
      last.row <- last.row[!duplicated(out[last.row,"time"])]
      # Average over the final day:
      MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred"] <- mean(out[last.row,"Cplasma"])
      MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred"] <- mean(out[last.row,"Cfplasma"])
      MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred"] <- 
        mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])

      # Repeat the analysis without the adjustment to fetal Fup:
      p <- parameterize_fetal_pbtk(dtxsid=this.id,
                                                    fetal_fup_adjustment =FALSE,
                                                    suppress.messages = TRUE)
      out <- solve_fetal_pbtk(
        parameters=p,
        dose=0,
        times=times,
        daily.dose=1,
        doses.per.day=3,
        output.units = "uM",
        maxsteps=1e7,
        suppress.messages = TRUE)

      last.row <- which(out[,"time"]>279) # The whole final day
      last.row <- last.row[!duplicated(out[last.row,"time"])]
      MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred.nofup"] <- mean(out[last.row,"Cplasma"])
      MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred.nofup"] <- mean(out[last.row,"Cfplasma"])
      MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred.nofup"] <- 
        mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
    }
  }
}

# Something is wrong with cotinine:
MFdata.httk <- subset(MFdata.httk,Chemical!="Cotinine")

max.chem <- MFdata.httk[which(MFdata.httk$MFratio==max(MFdata.httk$MFratio)),]
min.chem <- MFdata.httk[which(MFdata.httk$MFratio==min(MFdata.httk$MFratio)),]
cat(paste("The minimum observed ratio was ",
  signif(min.chem[,"MFratio"],2),
  " for ",
  min.chem[,"Chemical"],
  " and the maximum was ",
  signif(max.chem[,"MFratio"],2),
  " for ",
  max.chem[,"Chemical"],
  ".\n",sep=""))

Aylward Ratio figures:

Clean up repeated observations:

MFdata.main <- NULL
MFdata.outliers <- NULL
for (this.id in unique(MFdata.httk$DTXSID))
{
  this.subset <- subset(MFdata.httk,DTXSID==this.id)
  this.row <- this.subset[1,]
  this.row$N.obs <- dim(this.subset)[1]
  this.row$MFratio <- median(this.subset$MFratio)
  this.row$MFratio.Q25 <- quantile(this.subset$MFratio,0.25)
  this.row$MFratio.Q75 <- quantile(this.subset$MFratio,0.75)
  MFdata.main <- rbind(MFdata.main,this.row)
  this.subset <- subset(this.subset,
    MFratio<this.row$MFratio.Q25 |
    MFratio>this.row$MFratio.Q75)
  MFdata.outliers <- rbind(MFdata.outliers,this.subset) 
}

Cord Blood Ratio Predictions without Fup Adjustement:

FigCa  <- ggplot(data=MFdata.main) +
  geom_segment(color="grey",aes(
    x=MFratio.pred.nofup,
    y=MFratio.Q25,
    xend=MFratio.pred.nofup,
    yend=MFratio.Q75))+    
  geom_point(aes(
    x=MFratio.pred.nofup,
    y=MFratio,
    shape=Chemical.Category,
    color=Chemical.Category),
    size=4)   +
  scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+ 
  geom_point(data=MFdata.outliers,aes(
    x=MFratio.pred.nofup,
    y=MFratio,
    shape=Chemical.Category,
    color=Chemical.Category),
    size=2)   +
  xlim(0,2) +
  ylim(0,3) +
#   geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
#   scale_y_log10(label=scientific_10) +
#,limits=c(10^-7,100)) +
#   scale_x_log10(label=scientific_10) +
#   ,limits=c(10^-7,100)) +
#    annotation_logticks() + 
  geom_abline(slope=1, intercept=0) + 
#    geom_abline(slope=1, intercept=1,linetype="dashed") + 
#    geom_abline(slope=1, intercept=-1,linetype="dashed") + 
  ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) + 
  xlab("Generic PBTK Predicted Ratio") +
#    scale_color_brewer(palette="Set2") + 
  theme_bw()  +
  theme(legend.position="bottom")+
  theme( text  = element_text(size=14))+ 
  theme(legend.text=element_text(size=10))+ 
  guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+ 
  guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))

print(FigCa)

Generate text for results section:

cat(paste("In Figure 4 we compare predictions made with our high-throughput \
human gestational PBTK model with experimental observations on a per chemical \
basis wherever we had both in vitro HTTK data and in vivo observations (",
length(unique(MFdata.main$DTXSID)),
" chemicals).\n",sep=""))


repeats <- subset(MFdata.main,N.obs>1)

cat(paste("Multiple observations were available for ",
  dim(repeats)[1],
  " of the chemicals,\n",sep=""))


max.chem <- as.data.frame(repeats[which(repeats$MFratio==max(repeats$MFratio)),])
min.chem <- as.data.frame(repeats[which(repeats$MFratio==min(repeats$MFratio)),])

cat(paste("However, among the chemicals with repeated observations, the median \
observations only ranged from ",
  signif(min.chem[,"MFratio"],2),
  " for ",
  min.chem[,"Chemical"],
  " to ",
  signif(max.chem[,"MFratio"],2),
  " for ",
  max.chem[,"Chemical"],
  ".\n",sep=""))

max.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==max(MFdata.main$MFratio.pred,na.rm=TRUE)),])
min.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==min(MFdata.main$MFratio.pred,na.rm=TRUE)),])

cat(paste("The predictions for all chemicals ranged from ",
  signif(min.chem[,"MFratio.pred"],2),
  " for ",
  min.chem[,"Chemical"],
  " to ",
  signif(max.chem[,"MFratio.pred"],2),
  " for ",
  max.chem[,"Chemical"],
  ".\n",sep=""))



fit1 <- lm(data=MFdata.main,MFratio~MFratio.pred.nofup)
summary(fit1)
RMSE(fit1)

fit1sub <- lm(data=subset(MFdata.main,
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons"))),
  MFratio~MFratio.pred.nofup)
summary(fit1sub)
RMSE(fit1sub)

Cord Blood Ratio Predictions with Fup Adjustement:

FigCb  <- ggplot(data=MFdata.main) +
  geom_segment(color="grey",aes(
    x=MFratio.pred,
    y=MFratio.Q25,
    xend=MFratio.pred,
    yend=MFratio.Q75))+    
  geom_point(aes(
    x=MFratio.pred,
    y=MFratio,
    shape=Chemical.Category,
    color=Chemical.Category),
    size=3)   +
  scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+ 
  geom_point(data=MFdata.outliers,aes(
    x=MFratio.pred,
    y=MFratio,
    shape=Chemical.Category,
    color=Chemical.Category),
    size=1)   +
  xlim(0,2) +
  ylim(0,3) +
#   geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
#   scale_y_log10(label=scientific_10) +
#,limits=c(10^-7,100)) +
#   scale_x_log10(label=scientific_10) +
#   ,limits=c(10^-7,100)) +
#    annotation_logticks() + 
  geom_abline(slope=1, intercept=0) + 
#    geom_abline(slope=1, intercept=1,linetype="dashed") + 
#    geom_abline(slope=1, intercept=-1,linetype="dashed") + 
  ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) + 
  xlab(expression(paste(italic("In vitro")," Predicted Ratio"))) +
#    scale_color_brewer(palette="Set2") + 
  theme_bw()  +
  theme(legend.position="bottom")+
  theme( text  = element_text(size=14))+ 
  theme(legend.text=element_text(size=10))+ 
  guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+ 
  guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))

print(FigCb)

Examine performance when excluding certain chemical classes:

# Mean logHenry's law constant for PAH's:
mean(subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFdata.main,Chemical.Category=="Polyaromatic Hydrocarbons")$DTXSID)$logHenry)

nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$DTXSID
fluoros <- chem.physical_and_invitro.data$DTXSID[regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))!=-1]

fit2 <- lm(data=MFdata.main,MFratio~MFratio.pred)
summary(fit2)
RMSE(fit2)

fit2sub <- lm(data=subset(MFdata.main,
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons"))),
  MFratio~MFratio.pred)
summary(fit2sub)

RMSE(fit2sub)

cat(paste("When volatile and fluorinated chemicals are omitted only ",
  dim(subset(MFdata.main,
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons"))))[1],
  " evaluation chemicals remain, but the R2 is ",
  signif(summary(fit1sub)$adj.r.squared,2),
  " and the RMSE is ",
  signif(RMSE(fit1sub),2),
  " without the correction. When the fetal plasma fraction unbound correction is used, the predictivity decreases, slightly: R2 is ",
  signif(summary(fit2sub)$adj.r.squared,2),
  " and the RMSE is ",
  signif(RMSE(fit2sub),2),
  " for the non-volatile, non-fluorinated chemicals.\n",
  sep=""))



cat(paste("We compare the RMSE for our predictions to the standard deviation \
of the observations ",
  signif(sd(MFdata.main$MFratio)[1],2),
  " (",
  signif(sd(subset(MFdata.main,
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons")))$MFratio),2),
  " for non PAH or fluorinated compounds).\n",sep=""))

cat(paste("The average standard deviation for chemicals with repeated observations was ",
  signif(sd(subset(MFdata.main,N.obs>1)$MFratio)[1],2),
  " (",
  signif(sd(subset(MFdata.main,
  N.obs > 1 &
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons")))$MFratio),2),
  " for non PAH or fluorinated compounds).\n",sep=""))

fit3 <- lm(data=repeats,MFratio~MFratio.pred.nofup)
summary(fit3)
RMSE(fit3)

fit3sub <- lm(data=subset(MFdata.main, N.obs > 1 &
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons"))),
  MFratio~MFratio.pred.nofup)
summary(fit3sub)


fit4 <- lm(data=subset(MFdata.main,N.obs > 1),MFratio~MFratio.pred)
summary(fit4)
RMSE(fit4)

fit4sub <- lm(data=subset(MFdata.main, N.obs > 1 &
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons"))),
  MFratio~MFratio.pred)
summary(fit4sub)

cat(paste("Overall, we observed relatively poor correlation (R2 = ",
  signif(summary(fit1)$adj.r.squared,2),
  ", RMSE = ",
  signif(RMSE(fit1),2),
  ") without our fetal fup correction that was unchanged with that assumption (R2 = ",
  signif(summary(fit2)$adj.r.squared,2),
  ", RMSE = ",
  signif(RMSE(fit2),2),
  ").\n",sep=""))

repeats <-subset(MFdata.main,N.obs > 1)
cat(paste("The RMSE of the predictions for the ",
  dim(subset(repeats,!(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons"))))[1],
  " non-PAH and non-fluorinated compounds with repeated observations is ",
  signif(RMSE(fit4sub),2),
  " with the fup correction and ",
  signif(RMSE(fit3sub),2),
  " without.\n",sep=""))

cat(paste("These values are close to the standard deviation of the mean but the p-value for the chemicals with repeated observations is ",
  signif(pf(
    summary(fit4sub)$fstatistic[1],
    summary(fit4sub)$fstatistic[2],
    summary(fit4sub)$fstatistic[3]),2),    
" indicating some value for the predictive model, albeit for only ",
dim(subset(repeats,!(Chemical.Category %in% c(
  "Fluorinated compounds",
  "Polyaromatic Hydrocarbons"))))[1],
 " chemicals",sep=""))



Fig4.table <- data.frame()
Fig4.table["Number of Chemicals","All Fig 4"] <- length(unique(MFdata.httk$DTXSID))
Fig4.table["Number of Observations","All Fig 4"] <- dim(MFdata.httk)[1]
Fig4.table["Observed Mean (Min - Max)","All Fig 4"] <- paste(
  signif(mean(MFdata.httk$MFratio),2)," (",
  signif(min(MFdata.httk$MFratio),2)," - ",
  signif(max(MFdata.httk$MFratio),2),")",sep="")
Fig4.table["Observed Standard Deviation","All Fig 4"] <- signif(sd(MFdata.httk$MFratio),2) 
Fig4.table["Predicted Mean (Min - Max)","All Fig 4"] <-  paste(
  signif(mean(MFdata.main$MFratio.pred,na.rm=TRUE),2)," (",
  signif(min(MFdata.main$MFratio.pred,na.rm=TRUE),2)," - ",
  signif(max(MFdata.main$MFratio.pred,na.rm=TRUE),2),")",sep="") 
Fig4.table["RMSE","All Fig 4"] <- signif(RMSE(fit2),2)
Fig4.table["R-squared","All Fig 4"] <- signif(summary(fit2)[["adj.r.squared"]],2)
Fig4.table["p-value","All Fig 4"] <- signif(summary(fit2)[["coefficients"]]["MFratio.pred",4],2)

MFdata.sub1 <- subset(MFdata.httk,
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons")))

Fig4.table["Number of Chemicals","No PFAS/PAH"] <- length(unique(MFdata.sub1$DTXSID))
Fig4.table["Number of Observations","No PFAS/PAH"] <- dim(MFdata.sub1)[1]
Fig4.table["Observed Mean (Min - Max)","No PFAS/PAH"] <- paste(
  signif(mean(MFdata.sub1$MFratio,na.rm=TRUE),2)," (",
  signif(min(MFdata.sub1$MFratio,na.rm=TRUE),2)," - ",
  signif(max(MFdata.sub1$MFratio,na.rm=TRUE),2),")",sep="")
Fig4.table["Observed Standard Deviation","No PFAS/PAH"] <- signif(sd(MFdata.sub1$MFratio),2) 

MFdata.sub2 <- subset(MFdata.main,
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons")))

Fig4.table["Predicted Mean (Min - Max)","No PFAS/PAH"] <-  paste(
  signif(mean(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," (",
  signif(min(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," - ",
  signif(max(MFdata.sub2$MFratio.pred,na.rm=TRUE),2),")",sep="") 
Fig4.table["RMSE","No PFAS/PAH"] <- signif(RMSE(fit2sub),2)
Fig4.table["R-squared","No PFAS/PAH"] <- signif(summary(fit2sub)[["adj.r.squared"]],2)
Fig4.table["p-value","No PFAS/PAH"] <- signif(summary(fit2sub)[["coefficients"]]["MFratio.pred",4],2)


MFdata.sub3 <- subset(MFdata.main,N.obs > 1 &
  !(Chemical.Category %in% c(
    "Fluorinated compounds",
    "Polyaromatic Hydrocarbons")))

Fig4.table["Number of Chemicals","Replicates Only"] <- length(unique(MFdata.sub3$DTXSID))
Fig4.table["Number of Observations","Replicates Only"] <- dim(MFdata.sub3)[1]
Fig4.table["Observed Mean (Min - Max)","Replicates Only"] <- paste(
  signif(mean(MFdata.sub3$MFratio),2)," (",
  signif(min(MFdata.sub3$MFratio),2)," - ",
  signif(max(MFdata.sub3$MFratio),2),")",sep="")
Fig4.table["Observed Standard Deviation","Replicates Only"] <- signif(sd(MFdata.sub3$MFratio),2) 

Fig4.table["Predicted Mean (Min - Max)","Replicates Only"] <-  paste(
  signif(mean(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," (",
  signif(min(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," - ",
  signif(max(MFdata.sub3$MFratio.pred,na.rm=TRUE),2),")",sep="") 
Fig4.table["RMSE","Replicates Only"] <- signif(RMSE(fit4sub),2)
Fig4.table["R-squared","Replicates Only"] <- signif(summary(fit4sub)[["adj.r.squared"]],2)
Fig4.table["p-value","Replicates Only"] <- signif(summary(fit4sub)[["coefficients"]]["MFratio.pred",4],2)

write.csv(Fig4.table[1:6,],file="Fig4Table.csv")

Evaluation of Area Under the Curve Predictions

Dallmann et al., 2018 includes descriptions of toxicokinetics summary statistics, including time-integrated plasma concentrations (area under the curve or AUC) for six drugs administered to a sample of subjects including both pregnant and non-pregnant women. Additional data from X and Y for two chemicals among the httk library were collected from the literature.

#TKstats <- as.data.frame(read_excel("MaternalFetalAUCData.xlsx"))
TKstats <- httk::pregnonpregaucs


ids <- unique(TKstats$DTXSID)
cat(paste(sum(ids %in% get_cheminfo(
  model="fetal_pbtk",
  info="dtxsid",
  suppress.messages=TRUE)),
          "chemicals for which the model fetal_pbtk can run that are in the Dallmann data set."))


TKstats.Cmax <- subset(TKstats,DTXSID!="" & Parameter=="Cmax")
TKstats <- subset(TKstats,DTXSID!="" & Parameter %in% c("AUCinf","AUClast"))

ids <- unique(TKstats$DTXSID)
cat(paste(sum(ids %in% get_cheminfo(
  model="fetal_pbtk",
  info="dtxsid",
  suppress.messages=TRUE)),
          "chemicals for which the model fetal_pbtk can run that have AUCs in the Dallmann data set."))


# Assumed body weight (kg) from Kapraun 2019
BW.nonpreg <- 61.103

#TKstats$Dose <- TKstats$Dose/BW
#TKstats$Dose.Units <- "mg/kg"
colnames(TKstats)[colnames(TKstats)=="Observed Pregnant"] <- "Observed.Pregnant"
colnames(TKstats)[colnames(TKstats)=="Observed Pregnant5"] <- "Observed.Pregnant5"
colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant"] <- "Observed.Non.Pregnant"
colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant2"] <- "Observed.Non.Pregnant2"
colnames(TKstats)[colnames(TKstats)=="Dose Route"] <- "Dose.Route"

Predict AUC

The circumstances of the dosing varied slightly between drugs and pregnancy status required variation in simulated dose regimen as summarized in Table 12. The function solve_fetal_pbtk was used to determine the time-integrated plasma concentration (area under the curve, or AUC) for the mothers both when pregnant and non-pregnant.

for (this.id in unique(TKstats$DTXSID))
{
  if (any(regexpr("ng",TKstats[TKstats$DTXSID==this.id,"Dose Units"])!=-1))
  {
  }  
  if (this.id %in% get_cheminfo(
    info="DTXSID",
    model="fetal_pbtk",
    suppress.messages=TRUE))
  {
    this.subset <- subset(TKstats,DTXSID==this.id)
    p <- parameterize_pbtk(
      dtxsid=this.id,
      suppress.messages=TRUE)
    p$hematocrit <- 0.39412 # Kapraun 2019 (unitless)
    p$Rblood2plasma <- calc_rblood2plasma(
      parameters=p,
      suppress.messages=TRUE)
    p$BW <- BW.nonpreg 
    p$Qcardiacc <- 301.78 / p$BW^(3/4) # Kapraun 2019 (L/h/kg^3/4)
    for (this.route in unique(this.subset$Dose.Route))
    {
      this.route.subset <- subset(this.subset, Dose.Route==this.route)
      if (this.route.subset[1,"Gestational.Age.Weeks"] > 12)
      {
        this.dose <- this.route.subset$Dose
        # Non-pregnant PBPK:
        out.nonpreg <- solve_pbtk(
          parameters=p,
          times = seq(0, this.route.subset[1,"NonPreg.Duration.Days"],
                      this.route.subset[1,"NonPreg.Duration.Days"]/100),
          dose=this.dose/BW.nonpreg, # mg/kg
          daily.dose=NULL,
          iv.dose=(this.route=="iv"),
          output.units="uM",
          suppress.messages=TRUE)
        # Pregnant PBPK:
        BW.preg <- suppressWarnings(calc_maternal_bw(
          week=this.route.subset[1,"Gestational.Age.Weeks"]))
        out.preg <- solve_fetal_pbtk(
          dtxsid=this.id,
          times = seq(
            this.route.subset[1,"Gestational.Age.Weeks"]*7, 
            this.route.subset[1,"Gestational.Age.Weeks"]*7 + 
              this.route.subset[1,"Preg.Duration.Days"], 
            this.route.subset[1,"Preg.Duration.Days"]/100),
          dose=this.dose/BW.preg, # mg/kg
          iv.dose=(this.route=="iv"),
          daily.dose=NULL,
          output.units="uM",
          suppress.messages=TRUE)
        if (any(regexpr("AUC",this.subset$Parameter)!=-1))
        {
          TKstats[TKstats$DTXSID==this.id &
                    TKstats$Dose.Route == this.route &
                    regexpr("AUC",TKstats$Parameter)!=-1, 
                  "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"AUC"])*24 #uMdays->uMh                        
          TKstats[TKstats$DTXSID==this.id &
                    TKstats$Dose.Route == this.route &
                    regexpr("AUC",TKstats$Parameter)!=-1, 
                  "Predicted.Pregnant.httk"] <- max(out.preg[,"AUC"])*24 # uM days -> uM h 
        }
        if (any(regexpr("Cmax",this.subset$Parameter)!=-1))
        {
          TKstats[TKstats$DTXSID==this.id &
                    TKstats$Dose.Route == this.route &
                    regexpr("Cmax",TKstats$Parameter)!=-1, 
                  "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"Cplasma"])
          TKstats[TKstats$DTXSID==this.id &
                    TKstats$Dose.Route == this.route &
                    regexpr("Cmax",TKstats$Parameter)!=-1, 
                  "Predicted.Pregnant.httk"] <- max(out.preg[,"Cfplasma"])        
        }
      }
    }
  }
}

TKstats$Ratio.obs <- TKstats$Observed.Pregnant / TKstats$Observed.Non.Pregnant
TKstats$Ratio.httk <- TKstats$Predicted.Pregnant.httk / TKstats$Predicted.Non.Pregnant.httk
TKstats <- subset(TKstats, !is.na(TKstats$Ratio.httk))

write.csv(TKstats,file="Table-Dallmann2018.csv",row.names=FALSE)

Generate AUC Predicted vs. Observed Figure

FigEa  <- ggplot(data=TKstats) +
  geom_point(aes(
    y=Observed.Non.Pregnant2,
    x=Predicted.Non.Pregnant.httk,
    shape=Dose.Route,
    alpha=Drug,
    fill=Drug),
    size=5)   +
  geom_abline(slope=1, intercept=0) +
  geom_abline(slope=1, intercept=1, linetype=3) + 
  geom_abline(slope=1, intercept=-1, linetype=3) + 
  scale_shape_manual(values=c(22,24))+ 
  xlab("httk Predicted (uM*h)") + 
  ylab("Observed") +
  scale_x_log10(#limits=c(10^-3,10^3),
    label=scientific_10) +
  scale_y_log10(#limits=c(10^-3,10^3),
    label=scientific_10) +
  ggtitle("Non-Pregnant") +
  theme_bw()  +
  theme(legend.position="none")+
  theme( text  = element_text(size=12)) + 
  annotate("text", x = 0.1, y = 1000, label = "A", fontface =2)

#print(FigEa)
cat(paste(
  "For the Dallman et al. (2018) data we observe an average-fold error (AFE)\n\ of",
  signif(mean(TKstats$Predicted.Non.Pregnant.httk/TKstats$Observed.Non.Pregnant2),2),
  "and a RMSE (log10-scale) of",
  signif((mean((log10(TKstats$Predicted.Non.Pregnant.httk)-log10(TKstats$Observed.Non.Pregnant2))^2))^(1/2),2),
  "for non-pregnant woman.\n"))

FigEb  <- ggplot(data=TKstats) +
  geom_point(aes(
    y=Observed.Pregnant5,
    x=Predicted.Pregnant.httk,
    shape=Dose.Route,
    alpha=Drug,
    fill=Drug),
    size=5)   +
      geom_abline(slope=1, intercept=0) +
  geom_abline(slope=1, intercept=1, linetype=3) + 
  geom_abline(slope=1, intercept=-1, linetype=3) +
  scale_shape_manual(values=c(22,24))+ 
  xlab("httk Predicted (uM*h)") + 
  ylab("Observed") +
  scale_x_log10(#limits=c(10^-5,10^3),
                label=scientific_10) +
  scale_y_log10(#limits=c(10^-5,10^3),
                label=scientific_10) +
  ggtitle("Pregnant")+
  theme_bw()  +
  theme(legend.position="none")+
  theme( text  = element_text(size=12)) + 
  annotate("text", x = 0.1, y = 300, label = "B", fontface =2)

#print(FigEb)
cat(paste(
  "We observe an AFE of",
  signif(mean(TKstats$Predicted.Pregnant.httk/TKstats$Observed.Pregnant5),2),
  "and a RMSE (log10-scale) of",
  signif((mean((log10(TKstats$Predicted.Pregnant.httk)-log10(TKstats$Observed.Pregnant5))^2))^(1/2),2),
  "for pregnant woman.\n"))



FigEc  <- ggplot(data=TKstats) +
  geom_point(aes(
    x=Predicted.Non.Pregnant.httk,
    y=Predicted.Pregnant.httk,
    shape=Dose.Route,
    alpha=Drug,
    fill=Drug),
    size=5)   +
      geom_abline(slope=1, intercept=0) +
  geom_abline(slope=1, intercept=1, linetype=3) + 
  geom_abline(slope=1, intercept=-1, linetype=3) +
  scale_shape_manual(values=c(22,24),name="Route")+ 
  ylab("httk Pregnant (uM*h)") + 
  xlab("httk Non-Pregnant (uM*h)") + 
  scale_x_log10(#limits=c(10^-1,10^2),
                label=scientific_10) +
  scale_y_log10(#limits=c(10^-1,10^2),
                label=scientific_10) +
  ggtitle("Model Comparison") +
  theme_bw()  +
  theme(legend.position="left")+
  theme( text  = element_text(size=12))+ 
  theme(legend.text=element_text(size=10))+
  guides(fill=guide_legend(ncol=2,override.aes=list(shape=21)),alpha=guide_legend(ncol=2),shape=guide_legend(ncol=2))
 print(FigEc)   
FigEc <- get_legend(FigEc)


FigEd  <- ggplot(data=subset(TKstats,!is.na(Ratio.httk))) +
  geom_point(aes(
    y=Ratio.obs,
    x=Ratio.httk,
    shape=Dose.Route,
    alpha=Drug,
    fill=Drug),
    size=5)   +
  scale_shape_manual(values=c(22,24))+ 
  xlab("httk Predicted") + 
  ylab("Observed") +
  scale_x_continuous(limits=c(0.25,3)) +
  scale_y_continuous(limits=c(0.25,3)) +
  geom_abline(slope=1, intercept=0) +
 geom_abline(slope=1, intercept=1, linetype=3) + 
  geom_abline(slope=1, intercept=-1, linetype=3) + 
  ggtitle("Ratio") +
  theme_bw()  +
  theme(legend.position="none")+
  theme( text  = element_text(size=12)) + 
  annotate("text", x = 0.5, y = 2.75, label = "C", fontface =2)

dev.new()
grid.arrange(FigEa,FigEb,FigEc,FigEd,nrow=2)

write.csv(subset(TKstats,
  !is.na(Ratio.httk)),
  file="DallmanTable.txt")

Evaluation of Partition Coefficient Predictions

For each tissue, the partition coefficient (which describes the ratio of the concentration in the tissue to concentration of unbound chemical in the plasma at equilibrium) is a constant. We calculate each partition coefficient using the method of Schmitt, 2008 as described by Pearce et al., 2017. The partition coefficient for any given type of tissue (for example, liver tissue) depends on fraction unbound in plasma ($f_{up}^m$ or $f_{up}^f$), so in general these differ for mother and fetus.

Curley et al., 1969 compiled data on the concentration of a variety of pesticides in the cord blood of newborns and in the tissues of infants that were stillborn.

Curley <- as.data.frame(read_excel("Curley1969.xlsx"))
dim(Curley)
Curley.compounds <- Curley[1,4:13]
Curley <- Curley[4:47,]
colnames(Curley)[1] <- "Tissue"
colnames(Curley)[2] <- "N"
colnames(Curley)[3] <- "Stat"

The ratio of chemical concentration in tissue ($C_y^f$) vs. blood ($C_{venb}^f$) was related to the tissue-to-unbound-plasma concentration partition coefficients used in the PBTK model as $$ (C_y^f )/(C_{venb}^f )=(C_y^f)/(R_{b:p}^f \times C_{plas}^f )=(C_y^f \times f_{up}^f)/(R_{b:p}^f \times C_{up}^f)=(f_{up}^f)/(R_{b:p}^f) \times K_y^f $$ where $C_{up}^f$ denotes the concentration of substance unbound in the fetal plasma.

Curley.pcs <- NULL
cord.blood <- subset(Curley, Tissue == "Cord Blood") 
suppressWarnings(
for (this.tissue in unique(Curley$Tissue))
  if (this.tissue != "Cord Blood")
  {
    this.subset <- subset(Curley, Tissue == this.tissue)
    for (this.chemical in colnames(Curley)[4:13])
    {
      if (!is.na(as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical])) &
        !is.na(as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])))
      {
        this.row <- data.frame(
          Compound = Curley.compounds[,this.chemical],
          DTXSID = this.chemical,
          Tissue = this.tissue,
          PC = as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical]) /
            as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
          )
        Curley.pcs <- rbind(Curley.pcs,this.row)
      } else if (!is.na((as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]))) &
        !is.na((as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical]))))
      {
        this.row <- data.frame(
          Compound = Curley.compounds[,this.chemical],
          DTXSID = this.chemical,
          Tissue = this.tissue,
          PC = as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]) /
            as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
          )
        Curley.pcs <- rbind(Curley.pcs,this.row)
      }
    }
  }
)
Curley.pcs[Curley.pcs$Tissue=="Lungs","Tissue"] <- "Lung"

pearce2017 <- read_excel("PC_Data.xlsx")
pearce2017 <- subset(pearce2017,DTXSID%in%Curley.pcs$DTXSID)
any(pearce2017$DTXSID%in%Curley.pcs$DTXSID)
print("None of the Curley chems were in the Pearce 2017 PC predictor evaluation.")

Partition coefficients were measured for tissues, including placenta, in vitro by Csanady et al., 2002 for Bisphenol A and Diadzen.

csanadybpa <- read_excel("Csanady2002.xlsx",sheet="Table 2")
csanadybpa$Compound <- "Bisphenol A"
csanadybpa$DTXSID <- "DTXSID7020182"
csanadybpa <- csanadybpa[,c("Compound","DTXSID","Tissue","Mean")]
colnames(csanadybpa) <- colnames(Curley.pcs)

csanadydaid <- read_excel("Csanady2002.xlsx",sheet="Table 3",skip=1)
csanadydaid$Compound <- "Daidzein"
csanadydaid$DTXSID <- "DTXSID9022310"
csanadydaid <- csanadydaid[,c("Compound","DTXSID","Tissue","Mean...6")]
colnames(csanadydaid) <- colnames(Curley.pcs)

Curley.pcs <- rbind(Curley.pcs,csanadybpa,csanadydaid)
Curley.pcs <- subset(Curley.pcs,Tissue!="Mammary gland")

Three of the chemicals studied by Curley et al., 1969 were modeled by Weijs et al., 2013 using the same partition coefficients for mother and fetus. The values used represented "prior knowledge" summarizing the available literature.

weijstab3 <- read_excel("Weijs2013.xlsx",sheet="Table3",skip=1)
weijstab3 <- subset(weijstab3, !is.na(Compound) & !is.na(Tissue))
weijstab4 <- read_excel("Weijs2013.xlsx",sheet="Table4",skip=1)
weijstab4 <- subset(weijstab4, !is.na(Compound) & !is.na(Tissue))

for (this.compound in unique(weijstab3$Compound))
  for (this.tissue in unique(weijstab3$Tissue))
  {
    Curley.pcs[
      Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
      "Weijs2013"] <- weijstab3[weijstab3$Compound==this.compound &
                                  weijstab3$Tissue==this.tissue,"value"]

  }

for (this.compound in unique(weijstab4$Compound))
  for (this.tissue in unique(weijstab4$Tissue))
  {
    Curley.pcs[
      Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
      "Weijs2013"] <- weijstab4[weijstab4$Compound==this.compound &
                                  weijstab4$Tissue==this.tissue,"value"]

  }

write.csv(Curley.pcs,"PCs-table.csv",row.names=FALSE)

Predict partition coefficients

Curley.pcs <- httk::fetalpcs
dtxsid.list <- get_cheminfo(
  info="DTXSID",
  model="fetal_pbtk",
  suppress.messages=TRUE)

suppressWarnings(load_sipes2017())
for (this.chemical in unique(Curley.pcs$DTXSID))
  if (this.chemical %in% dtxsid.list)
  {
    this.subset <- subset(Curley.pcs,DTXSID==this.chemical)
    p <- parameterize_fetal_pbtk(dtxsid=this.chemical,
      fetal_fup_adjustment = FALSE,
      suppress.messages=TRUE,
      minimum.Funbound.plasma = 1e-16)
    fetal.blood.pH <- 7.26   
    Fup <- p$Fraction_unbound_plasma_fetus
    fetal_schmitt_parms <- parameterize_schmitt(
      dtxsid=this.chemical,
      suppress.messages=TRUE,
      minimum.Funbound.plasma = 1e-16)
    fetal_schmitt_parms$plasma.pH <- fetal.blood.pH
    fetal_schmitt_parms$Funbound.plasma <- Fup
    Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Fup"] <- Fup
    # httk gives tissue:fup (unbound chemical in plasma) PC's:
    fetal_pcs <- predict_partitioning_schmitt(
      parameters = fetal_schmitt_parms,
      suppress.messages=TRUE,
      model="fetal_pbtk",
      minimum.Funbound.plasma = 1e-16)
    fetal_pcs.nocal <- predict_partitioning_schmitt(
      parameters = fetal_schmitt_parms,
      regression=FALSE,
      suppress.messages=TRUE,
      model="fetal_pbtk",
      minimum.Funbound.plasma = 1e-16)
    out <- solve_fetal_pbtk(
      dtxsid = this.chemical,
      fetal_fup_adjustment =FALSE,
      suppress.messages=TRUE,
      minimum.Funbound.plasma = 1e-16)
    Rb2p <- out[dim(out)[1],"Rfblood2plasma"]
    Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Rb2p"] <- Rb2p
    # Convert to tissue:blood PC's:
    for (this.tissue in this.subset$Tissue)
      if (tolower(this.tissue) %in% 
        unique(subset(tissue.data,Species=="Human")$Tissue))
      {
        Curley.pcs[Curley.pcs$DTXSID==this.chemical &
          Curley.pcs$Tissue == this.tissue, "HTTK.pred.t2up"] <-
          fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]
        Curley.pcs[Curley.pcs$DTXSID==this.chemical &
          Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal.t2up"] <-
          fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]
        Curley.pcs[Curley.pcs$DTXSID==this.chemical &
          Curley.pcs$Tissue == this.tissue, "HTTK.pred"] <-
          fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
        Curley.pcs[Curley.pcs$DTXSID==this.chemical &
          Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal"] <-
          fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
      } else {
       print(this.tissue)
      }  
  } else print(this.chemical)
reset_httk()

cat(paste(
  "For the two placental observations (",
  signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"PC"],2),
  " for Bisphenol A and ",
  signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"PC"],2),
  " for Diadzen) the predictions were ",
  signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"HTTK.pred"],2),
  " and ",
  signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"HTTK.pred"],2),
  ", respectively.\n",sep=""))

Create Observed vs. Predicted Plot

First indicate compound:

FigFa  <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred.nocal))) +
  geom_point(size=2,aes(
    y=Weijs2013,
    x=HTTK.pred,    
    shape=Compound,
    color=Compound)) +
  geom_point(size=5,aes(
    y=PC,
    x=HTTK.pred,
    shape=Compound,
    color=Compound)) +
  geom_abline(slope=1, intercept=0, size=2) +
  geom_abline(slope=1, intercept=1, linetype=3, size=2) + 
  geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
  scale_shape_manual(values=rep(c(15,16,17,18),4))+
  xlab("Predicted Tissue:Blood Partition Coefificent") + 
  ylab("Observed") +
  scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
  scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
  theme_bw()  +
  theme(legend.position="bottom")+
  theme( text  = element_text(size=28))+ 
  theme(legend.text=element_text(size=12))+ 
   theme(legend.title = element_blank())+
  guides(shape=guide_legend(nrow=3,byrow=TRUE))

print(FigFa)

fitFa <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred))
RMSE(fitFa)
fitFb <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred.nocal))
RMSE(fitFb)

Then indicate tissues:

FigFb  <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred))) +
  geom_point(size=2,aes(
    y=Weijs2013,
    x=HTTK.pred,
    shape=Tissue,
    color=Tissue)) +
  geom_point(size=5, aes(
    y=PC,
    x=HTTK.pred,
    shape=Tissue,
    color=Tissue))   +
  geom_abline(slope=1, intercept=0, size=2) +
  geom_abline(slope=1, intercept=1, linetype=3, size=2) + 
  geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
  scale_shape_manual(values=rep(c(15,16,17,18),4))+
  xlab("Predicted Tissue:Blood Partition Coefificent") + 
  ylab("Observed") +
  scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
  scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
  theme_bw()  +
  theme(legend.position="bottom")+
  theme( text  = element_text(size=28))+ 
  theme(legend.text=element_text(size=20))+ 
   theme(legend.title = element_blank())+
  guides(shape=guide_legend(nrow=3,byrow=TRUE))

print(FigFb)

Compare httk partition coefficients with PKsim:

pksim.pcs <- as.data.frame(read_excel("PartitionCoefficients.xlsx"))
dim(pksim.pcs)
for (this.id in unique(pksim.pcs$DTXSID))
{
  httk.PCs <- predict_partitioning_schmitt(dtxsid=this.id,
                                           suppress.messages = TRUE)
  p <- 
    parameterize_fetal_pbtk(dtxsid=this.id, 
    suppress.messages = TRUE)
  httk.PCs[["Kplacenta2pu"]] <- p[["Kplacenta2pu"]] 
  httk.PCs[["Fup"]] <- p[["Funbound.plasma"]]
  lapply(httk.PCs,as.numeric)
  this.subset <- subset(pksim.pcs,DTXSID==this.id)
  for (this.tissue in unique(this.subset$Tissue))
  {
    if (any(regexpr(tolower(this.tissue),names(httk.PCs))!=-1))
    {
      pksim.pcs[pksim.pcs$DTXSID==this.id & pksim.pcs$Tissue==this.tissue,
        "HTTK.pred"] <-     
        httk.PCs[regexpr(tolower(this.tissue),names(httk.PCs))!=-1][[1]]*
        httk.PCs[["Fup"]][[1]]         
    }
  }
}
colnames(pksim.pcs)[colnames(pksim.pcs) ==
                      "Tissue-to-plasma partition coefficient"] <- "PKSim.pred"
colnames(pksim.pcs)[colnames(pksim.pcs) ==
        "Method for calculating tissue-to-plasma partition coefficients"] <- 
  "Method"
pksim.pcs[,"Method"] <- as.factor(pksim.pcs[,"Method"])
pksim.pcs[,"Drug"] <- as.factor(pksim.pcs[,"Drug"])
pksim.pcs[,"Tissue"] <- as.factor(pksim.pcs[,"Tissue"])


pksim.pcs[,"PKSim.pred"] <- as.numeric(
  pksim.pcs[,"PKSim.pred"])

pksim.pcs <- httk::pksim.pcs

pksim.fit1 <- lm(data=pksim.pcs,
                log10(PKSim.pred)~log10(HTTK.pred))
summary(pksim.fit1)

pksim.fit2 <- lm(data=pksim.pcs,
                log10(PKSim.pred)~log10(HTTK.pred)+
                  Drug+Tissue+Method)
summary(pksim.fit2)

Chemical Prioritization on the Basis of Fetal Brain Concentration

#Wangchems <- read_excel("Wang2018.xlsx",sheet=3,skip=2)
Wangchems <- httk::wang2018
maternal.list <- Wangchems$CASRN[Wangchems$CASRN %in% 
  get_cheminfo(model="fetal_pbtk",
  suppress.messages=TRUE)]
nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$CAS
nonfluoros <- chem.physical_and_invitro.data$CAS[
  regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))==-1]
maternal.list <- maternal.list[maternal.list %in% intersect(nonvols,nonfluoros)]

For the plot we need a data frame with a column "Ratio.pred" broken down by the different contributions to uncertainty being plotted (columne "Uncertainty"). We build this plot by making multiple versions of pred.table and concatonating them together at the end.

Make predictions for fetal:maternal plasma ratio using three daily doses from week 13 to full term

pred.table1 <- subset(get_cheminfo(
  info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
  model="fetal_pbtk"),
  CAS %in% maternal.list,
  suppress.messages=TRUE)
pred.table1$Compound <- gsub("\"","",pred.table1$Compound)    

for (this.cas in maternal.list)
{
  Fup <- subset(get_cheminfo(
    info="all",
    suppress.messages=TRUE),
    CAS==this.cas)$Human.Funbound.plasma
  # Check if Fup is provided as a distribution:
  if (regexpr(",",Fup)!=-1)
  {
    if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
      (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 & 
      as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
    {
      skip <- TRUE
    } else skip <- FALSE
    Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
    Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
    # These results are actually from a Bayesian posterior, but we can approximate that
    # they are normally distributed about the median to estimate a standard deviation:
    Fup.sigma <- (Fup.upper - Fup)/1.96
  # If it's not a distribution use defaults (see Wambaugh et al. 2019)
  } else if (as.numeric(Fup)== 0) 
  {
    skip <- TRUE
  } else {
    skip <- FALSE
    Fup <- as.numeric(Fup)
    Fup.sigma <- Fup*0.4
    Fup.upper  <- min(1,(1+1.96*0.4)*Fup)
  }

  # Only run the simulation if we have the necessary parameters:
  if (!skip)
  { 
    print(this.cas)
    out <- solve_fetal_pbtk(
      chem.cas=this.cas,
      dose=0,
      daily.dose=1,
      doses.per.day=3,
      fetal_fup_adjustenment=TRUE,
      suppress.messages=TRUE)
    last.row <- which(out[,"time"]>279)
    last.row <- last.row[!duplicated(out[last.row,"time"])]
    pred.table1[pred.table1$CAS==this.cas,"fup"] <- signif(Fup,3)
    pred.table1[pred.table1$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
    pred.table1[pred.table1$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
    pred.table1[pred.table1$CAS==this.cas,"Ratio.pred"] <- 
      signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
  }
}

Create final table holding all predictions for paper:

PC.table <- pred.table1
colnames(PC.table)[colnames(PC.table)=="Ratio.pred"] <- "R.plasma.FtoM"
pred.table1$Uncertainty <- "Predicted F:M Plasma Ratio"

Error in fetal:maternal ratio (inverse of maternal:fetal ratio):

pred.table3 <- pred.table1
pred.table3$Uncertainty <- "Plasma Error (Fig. 4)"
empirical.error <- RMSE(fit2sub)
for (this.cas in maternal.list)
{

  pred.table3[pred.table3$CAS==this.cas,"Ratio.pred"] <- signif(
    1/((1/pred.table1[pred.table3$CAS==this.cas,"Ratio.pred"])*
         (1-1.96*empirical.error)), 3)
}
# Update final table for paper:
PC.table$RMSE <- signif(RMSE(fit2sub),3)
PC.table$R.plasma.FtoM.upper <- pred.table3$Ratio.pred

brain:plasma partitioning

pred.table4 <- pred.table1
pred.table4$Uncertainty <- "Fetal Brain Partitioning"
for (this.cas in maternal.list)
{
  print(this.cas)
  p <- parameterize_fetal_pbtk(chem.cas=this.cas,
      fetal_fup_adjustment=TRUE,
      suppress.messages=TRUE)
  Kbrain2pu <- p$Kfbrain2pu
  fup <- pred.table4[pred.table4$CAS==this.cas,"fup"]
  pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"] <- signif(
     PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]*
    Kbrain2pu * fup, 3)
  PC.table[PC.table$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
  PC.table[PC.table$CAS==this.cas,"R.brain.FtoM"] <- 
    pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"]
}

brain:plasma uncertainty partitioning

pred.table5 <- pred.table1
pred.table5$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
for (this.cas in maternal.list)
{                                             
  p <- parameterize_fetal_pbtk(chem.cas=this.cas,
      fetal_fup_adjustment=TRUE,
      suppress.messages=TRUE)
  Kbrain2pu <- p$Kfbrain2pu

  fup <- pred.table5[pred.table5$CAS==this.cas,"fup"]
  sigma.fup <- pred.table5[pred.table5$CAS==this.cas,"fup.sigma"]
  Rmatfet <- 1/PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]
  Rbrain2matblood <- Kbrain2pu * fup / Rmatfet

# From Pearce et al. (2017) PC paper:
  sigma.logKbrain <- 0.647
  Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)

  error.Rmatfet <- PC.table[PC.table$CAS==this.cas,"RMSE"]
  sigma.Rbrain2matblood <- Rbrain2matblood * 
    (log(10)^2*sigma.logKbrain^2 +
    sigma.fup^2/fup^2 +
    error.Rmatfet^2/Rmatfet^2)^(1/2)
  Rbrain2matblood.upper <- Rbrain2matblood * 
    (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
    sigma.fup^2/fup^2 +
    error.Rmatfet^2/Rmatfet^2)^(1/2))
  pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"] <- 
    signif(Rbrain2matblood.upper,3)
  PC.table[PC.table$CAS==this.cas,"sigma.logKbrain"] <- 
    signif(sigma.logKbrain, 3)
  PC.table[PC.table$CAS==this.cas,"Kbrain2pu.upper"] <- 
    signif(Kbrain2pu.upper, 3)
  PC.table[PC.table$CAS==this.cas,"sigma.Rbrain2matblood"] <- 
    signif(sigma.Rbrain2matblood, 3)
  PC.table[PC.table$CAS==this.cas,"R.brain.FtoM.upper"] <- 
    pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"]
}

Combine different ratio predicitons together for figure

pred.levels <- pred.table5$Compound[order(pred.table5$Ratio.pred)]

pred.table <- rbind(
  pred.table1,
#  pred.table2,
  pred.table3,
  pred.table4,
  pred.table5)
pred.table$Compound <- factor(pred.table$Compound,
  levels = pred.levels)

pred.table$Uncertainty <- factor(pred.table$Uncertainty, 
  levels = c(pred.table1[1,"Uncertainty"],
#    pred.table2[1,"Uncertainty"],
    pred.table3[1,"Uncertainty"],
    pred.table4[1,"Uncertainty"],
    pred.table5[1,"Uncertainty"]))

Make prioritization figure:

#Wang 2018 confirmed 6 chemical identities:
confirmed.chemicals <- c(
  "2,4-Di-tert-butylphenol",
  "2,4-Dinitrophenol",
  "Pyrocatechol",
  "2'-Hydroxyacetophenone",
  "3,5-Di-tert-butylsalicylic acid",
  "4-Hydroxycoumarin"
  )
confirmed.chemicals <- c(
  "96-76-4",
  "19715-19-6",
  "51-28-5",
  "120-80-9",
  "118-93-4",
  "1076-38-6")


FigG  <- ggplot(data=pred.table) +
  geom_point(aes(
    x=Ratio.pred,
    y=Compound,
    color = Uncertainty,
    shape = Uncertainty),
    size=3)   +
    scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+ 
  scale_x_log10(limits=c(10^-2,10^3),label=scientific_10)+
  ylab(expression(paste(
    "Chemicals Found in Maternal Plasma by Wang   ",italic("et al.")," (2018)"))) + 
  xlab("Predicted Ratio to Maternal Plasma") +
  theme_bw()  +
#  theme(legend.position="bottom")+
  theme( text  = element_text(size=14))+ 
  theme(legend.text=element_text(size=10))#+ 
 # guides(color=guide_legend(nrow=4,byrow=TRUE))+ 
  #guides(shape=guide_legend(nrow=4,byrow=TRUE))
  #+
 # theme(legend.justification = c(0, 0), legend.position = c(0, 0))

print(FigG)

Look at ionization state:

for (this.col in 7:14)
  PC.table[,this.col] <- signif(PC.table[,this.col],3)

PC.table <- PC.table[order(PC.table$R.brain.FtoM.upper,decreasing=TRUE),]

for (this.row in 1:dim(PC.table)[1])
{
  out <- calc_ionization(
    pH=7.26,
    pKa_Donor=PC.table[this.row,"pKa_Donor"],
    pKa_Accept=PC.table[this.row,"pKa_Accept"])
  if (out$fraction_neutral>0.9) PC.table[this.row,"Charge_726"] <- "Neutral"
  else if (out$fraction_positive>0.1) PC.table[this.row,"Charge_726"] <- 
    paste(signif(out$fraction_positive*100,2),"% Positive",sep="")
  else if (out$fraction_negative>0.1) PC.table[this.row,"Charge_726"] <- 
    paste(signif(out$fraction_negative*100,2),"% Negative",sep="")
  else if (out$fraction_zwitter>0.1) PC.table[this.row,"Charge_726"] <- 
    paste(signif(out$fraction_zwitter*100,2),"% Zwitterion",sep="")
}

PC.table <- PC.table[,c(
  "Compound",
  "CAS",
  "DTXSID",
  "logP",
  "Charge_726",
  "R.plasma.FtoM",
  "RMSE",
  "R.plasma.FtoM.upper",
  "Kbrain2pu",
  "fup",
  "R.brain.FtoM",
  "Kbrain2pu.upper",
  "R.brain.FtoM.upper")]

write.csv(PC.table,
  file="WangTable.txt",
  row.names=FALSE)   

Evaluate any value added by PBPK model over just using partition coefficients:

MFbrainfit <- lm(R.brain.FtoM.upper~Kbrain2pu.upper,data=PC.table)
summary(MFbrainfit)

cat(paste("As expected, the predicted fetal brain to maternal blood ratio was strongly correlated with the predicted brain partition coefficient (R2 = ",
  signif(summary(MFbrainfit)$adj.r.square,2),
  ", p-value ",
  signif(summary(MFbrainfit)$coefficients[2,4],2),
  "). However, the PBTK simulation impacted the plasma and tissue concentrations such that ",
  dim(subset(PC.table, rank(R.brain.FtoM.upper) <  rank(Kbrain2pu.upper)))[1],
  " chemicals were ranked higher than they would have been using partitioning alone.",
  sep=""))

Repeat analysis without fetal protein binding change:

Make predictions for fetal:maternal plasma ratio using three daily doses from week 13 to full term

pred.table1.nofup <- subset(get_cheminfo(
  info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
  model="fetal_pbtk",
  suppress.messages=TRUE),
  CAS %in% maternal.list)
pred.table1.nofup$Compound <- gsub("\"","",pred.table1.nofup$Compound)    

for (this.cas in maternal.list)
{
  Fup <- subset(get_cheminfo(
    info="all",
    suppress.messages=TRUE),
    CAS==this.cas)$Human.Funbound.plasma
  # Check if Fup is provided as a distribution:
  if (regexpr(",",Fup)!=-1)
  {
    if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
      (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 & 
      as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
    {
      skip <- TRUE
    } else skip <- FALSE
    Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
    Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
    # These results are actually from a Bayesian posterior, but we can approximate that
    # they are normally distributed about the median to estimate a standard deviation:
    Fup.sigma <- (Fup.upper - Fup)/1.96
  # If it's not a distribution use defaults (see Wambaugh et al. 2019)
  } else if (as.numeric(Fup)== 0) 
  {
    skip <- TRUE
  } else {
    skip <- FALSE
    Fup <- as.numeric(Fup)
    Fup.sigma <- Fup*0.4
    Fup.upper  <- min(1,(1+1.96*0.4)*Fup)
  }

  # Only run the simulation if we have the necessary parameters:
  if (!skip)
  {  
    out <- solve_fetal_pbtk(
      chem.cas=this.cas,
      dose=0,
      daily.dose=1,
      doses.per.day=3,
      fetal_fup_adjustenment=FALSE,
      suppress.messages=TRUE)
    last.row <- which(out[,"time"]>279)
    last.row <- last.row[!duplicated(out[last.row,"time"])]
    pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup"] <- signif(Fup,3)
    pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
    pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
    pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"Ratio.pred"] <- 
      signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
  }
}

Create final table holding all predictions for paper:

PC.table.nofup <- pred.table1.nofup
colnames(PC.table.nofup)[colnames(PC.table.nofup)=="Ratio.pred"] <- "R.plasma.FtoM"
pred.table1.nofup$Uncertainty <- "Predicted F:M Plasma Ratio"

Error in fetal:maternal ratio (inverse of maternal:fetal ratio):

pred.table3.nofup <- pred.table1.nofup
pred.table3.nofup$Uncertainty <- "Plasma Error (Fig. 4)"
empirical.error <- RMSE(fit2sub)
for (this.cas in maternal.list)
{

  pred.table3.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
    1/((1/pred.table1.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"])*
         (1-1.96*empirical.error)), 3)
}
# Update final table for paper:
PC.table.nofup$RMSE <- signif(RMSE(fit2sub),3)
PC.table.nofup$R.plasma.FtoM.upper <- pred.table3.nofup$Ratio.pred

brain:plasma partitioning

pred.table4.nofup <- pred.table1.nofup
pred.table4.nofup$Uncertainty <- "Fetal Brain Partitioning"
for (this.cas in maternal.list)
{
  p <- parameterize_fetal_pbtk(chem.cas=this.cas,
      fetal_fup_adjustment=FALSE,
      suppress.messages=TRUE)
  Kbrain2pu <- p$Kfbrain2pu
  fup <- pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"fup"]
  pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
     PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]*
    Kbrain2pu * fup, 3)
  PC.table.nofup[PC.table.nofup$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
  PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.brain.FtoM"] <- 
    pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"]
}

brain:plasma uncertainty partitioning

pred.table5.nofup <- pred.table1.nofup
pred.table5.nofup$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
for (this.cas in maternal.list)
{                                             
  p <- parameterize_fetal_pbtk(chem.cas=this.cas,
      fetal_fup_adjustment=FALSE,
      suppress.messages=TRUE)
  Kbrain2pu <- p$Kfbrain2pu

  fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup"]
  sigma.fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup.sigma"]
  Rmatfet <- 1/PC.table[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]
  Rbrain2matblood <- Kbrain2pu * fup / Rmatfet

# From Pearce et al. (2017) PC paper:
  sigma.logKbrain <- 0.647
  Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)

  error.Rmatfet <- PC.table.nofup [PC.table.nofup $CAS==this.cas,"RMSE"]
  sigma.Rbrain2matblood <- Rbrain2matblood * 
    (log(10)^2*sigma.logKbrain^2 +
    sigma.fup^2/fup^2 +
    error.Rmatfet^2/Rmatfet^2)^(1/2)
  Rbrain2matblood.upper <- Rbrain2matblood * 
    (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
    sigma.fup^2/fup^2 +
    error.Rmatfet^2/Rmatfet^2)^(1/2))
  pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"] <- 
    signif(Rbrain2matblood.upper,3)
  PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.logKbrain"] <- 
    signif(sigma.logKbrain, 3)
  PC.table.nofup [PC.table.nofup $CAS==this.cas,"Kbrain2pu.upper"] <- 
    signif(Kbrain2pu.upper, 3)
  PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.Rbrain2matblood"] <- 
    signif(sigma.Rbrain2matblood, 3)
  PC.table.nofup [PC.table.nofup $CAS==this.cas,"R.brain.FtoM.upper"] <- 
    pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"]
}

case.study.fup.correct.table <- merge(PC.table,PC.table.nofup,by=c("Compound","CAS","DTXSID"))

MFbrainfit.fup <- lm(R.brain.FtoM.upper.x~R.brain.FtoM.upper.y,
                     data=case.study.fup.correct.table)
summary(MFbrainfit.fup)

cat(paste("The predictions for fetal brain to maternal blood ratio with or without the fetal plasma binding correction (Eqn. 1) were strongly correlated (R2 = ",
  signif(summary(MFbrainfit.fup)$adj.r.square,2),
  "). There were ",
  dim(subset(case.study.fup.correct.table, 
             rank(R.brain.FtoM.upper.x) > rank(R.brain.FtoM.upper.y)))[1],
  " chemicals that were ranked higher with the correction than without, with an average increase of ",
  signif(mean(
    case.study.fup.correct.table$R.brain.FtoM.upper.y /
      case.study.fup.correct.table$R.brain.FtoM.upper.x),2),
  " times when the plasma binding change was included.\n",
  sep=""))

The next two figures take a long time to generate and so are disabled by default.

Protein Binding Figures

The fetal fraction unbound (f_{up}^f i) is calculated from the maternal fraction unbound and the serum protein concentration ratio in infants vs. mothers based on Equation 6 of McNamara et al., 2019,

$$ f_{up}^f =(f_{up}^m )/(f_{up}^m +(P^f /P^m )*(1-f_{up}^m ) ) $$

in which the maternal fraction unbound, $f_{up}^m$, is assumed to be equal to the in vitro measured value for fraction unbound in plasma, and the ratio of protein concentrations $P^f /P^m$ depends on the identity of the dominant binding protein for the chemical (assumed to be either albumin or AAG). Lacking data to model the gestational kinetics of albumin and AAG concentrations, we used the concentrations at birth ([McNamara et al., 2019] (http://dx.doi.org/10.1208/ps040104)) to calculate a constant fetal $f_{up}$, using $P^f/P^m =0.777$ for albumin and $P^f/P^m = 0.456$ for AAG (McNamara et al., 2019). We determine the charge state of a compound separately for maternal and fetal plasma as a function of plasma pH (7.38 for maternal and 7.28 for fetal ([K.H. Lee, 1972] (http://dx.doi.org/10.1136/pgmj.48.561.405)) and chemical-specific predictions for ionization affinity (that is, pKa ([Strope et al., 2018] (http://dx.doi.org/10.1016/j.scitotenv.2017.09.033)) using the "httk" function "calc_ionization" ([Pearce et al., 2017] (http://dx.doi.org/10.1007/s10928-017-9548-7)). If the fraction of a chemical that is predicted to be in positive ionic form is greater than 50%, we treat the chemical as a base (which is in its conjugate acid form) and use only the maternal-to-infant ratio of AAG concentrations. Otherwise, we assume the chemical is neutral or an acid and use the ratio of albumin concentrations.

fup.table <- NULL
all.chems <- get_cheminfo(
  model="fetal_pbtk",
  info="all",
  suppress.messages=TRUE) 
# Get rid of median fup 0:
all.chems <- subset(all.chems,
  as.numeric(unlist(lapply(strsplit(
    all.chems$Human.Funbound.plasma,","),function(x) x[[1]])))!=0)
for (this.chem in all.chems[,"CAS"])
{
  temp <- parameterize_fetal_pbtk(chem.cas=this.chem,suppress.messages = TRUE)
  state <- calc_ionization(
      pH=7.26,
      pKa_Donor=temp$pKa_Donor,
      pKa_Accept=temp$pKa_Accept)
  if (state$fraction_positive > 0.5) this.charge <- "Positive"
  else if (state$fraction_negative > 0.5) this.charge <- "Negative"
  else this.charge <- "Neutral"
  this.row <- data.frame(DTXSID=all.chems[all.chems[,"CAS"]==this.chem,"DTXSID"],
    Compound=all.chems[all.chems[,"CAS"]==this.chem,"Compound"],
    CAS=this.chem,
    Fup.Mat.Pred = temp$Funbound.plasma,
    Fup.Neo.Pred = temp$Fraction_unbound_plasma_fetus,
    Charge = this.charge
    )
  fup.table <- rbind(fup.table,this.row)
}

fup.table[fup.table$Charge=="Positive","Charge"] <- paste("Positive (n=",
  sum(fup.table$Charge=="Positive"),
  ")",sep="")
fup.table[fup.table$Charge=="Negative","Charge"] <- paste("Negative (n=",
  sum(fup.table$Charge=="Negative"),
  ")",sep="")
fup.table[fup.table$Charge=="Neutral","Charge"] <- paste("Neutral (n=",
  sum(fup.table$Charge=="Neutral"),
  ")",sep="")

Plot the fetal protein binding changes predicted:

FigA  <- ggplot(data=fup.table) +
  geom_point(alpha=0.25, aes(
    x=Fup.Mat.Pred,
    y=Fup.Neo.Pred,
    shape=Charge,
    color=Charge),
    size=3)   +
  geom_abline(slope=1, intercept=0) + 
  ylab(expression(paste("Adjusted Neonate ",f[up]))) + 
  xlab(expression(paste(italic("In vitro")," Measured Adult ",f[up]))) +
   scale_x_log10(label=scientific_10) +
   scale_y_log10(label=scientific_10) +
  theme_bw()  +
  theme( text  = element_text(size=14)) 

print(FigA) 

Maternal-Fetal Predictions across HTTK Chemical Library:

times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))

MFratio.pred <- NULL
all.chems <- get_cheminfo(
  model="fetal_pbtk", 
  info=c("Compound","DTXSID","CAS","Funbound.plasma"),
  suppress.messages=TRUE)
for (this.cas in all.chems$CAS)
  if ((this.cas %in% nonvols) & 
    !(this.cas %in% fluoros))
{
  this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
  Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma 
  if (regexpr(",",Fup)!=-1)
  {
    if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
      (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 & 
      as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
    {
      skip <- TRUE
    } else skip <- FALSE
  } else if (Fup== 0) 
  {
    skip <- TRUE
  } else skip <- FALSE

  if (!skip)
  {  
    p <- parameterize_fetal_pbtk(dtxsid=this.id,
    fetal_fup_adjustment =TRUE,
    suppress.messages=TRUE)
    out <- solve_fetal_pbtk(
      parameters=p,
      dose=0,
      times=times,
      daily.dose=1,
      doses.per.day=3,
      output.units = "uM",
      suppress.messages=TRUE)
    last.row <- which(out[,"time"]>279)
    last.row <- last.row[!duplicated(out[last.row,"time"])]
    new.row <- data.frame(
      Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
      DTXSID = this.id,
      Mat.pred = mean(out[last.row,"Cplasma"]),
      Fet.pred = mean(out[last.row,"Cfplasma"]),
      MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
      )
    MFratio.pred <- rbind(MFratio.pred,new.row)
  }
  }

Histogram of maternal-fetal ratio

FigD <- ggplot(data=MFratio.pred)+
   geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+ 
  xlab("Maternal:Fetal Plasma Concentration Ratio") +
  ylab("Number of chemicals")+
    theme_bw()+
    theme( text  = element_text(size=14))


print(FigD)

Statistics on maternal-fetal ratio for full HTTK library

max.chem <- MFratio.pred[which(
  MFratio.pred$MFratio.pred==max(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
min.chem <- MFratio.pred[which(
  MFratio.pred$MFratio.pred==min(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
cat(paste("In Figure X we examine the ratios predicted for the ",
  dim(MFratio.pred)[1],
  " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
  sep=""))


cat(paste("We observe a median ratio of ",
  signif(median(MFratio.pred$MFratio.pred,na.rm=TRUE),3),
  ", with values ranging from ",
  signif(min.chem[,"MFratio.pred"],3),
  " for ",
  min.chem[,"DTXSID"],
  " to ",
  signif(max.chem[,"MFratio.pred"],3),
  " for ",
  max.chem[,"DTXSID"],
  ".\n",sep=""))

# Check out phys-chem > 1.6, < 1:
highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred>1.6)$DTXSID)

cat(paste("There are",
  dim(highratio)[1],
  "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))

# all highly bound
highratio$Compound
suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=TRUE)))


lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred<0.9)$DTXSID)
# No obvious pattern

Examine Maternal-Fetal ratio Predictions without fetaul Fup correction:

times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))

MFratio.pred.nofup <- NULL
all.chems <- get_cheminfo(
  model="fetal_pbtk", 
  info=c("Compound","DTXSID","CAS","Funbound.plasma"),
  suppress.messages=TRUE)
for (this.cas in all.chems$CAS)
  if ((this.cas %in% nonvols) & 
    !(this.cas %in% fluoros))
{
  this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
  Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma 
  if (regexpr(",",Fup)!=-1)
  {
    if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
      (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 & 
      as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
    {
      skip <- TRUE
    } else skip <- FALSE
  } else if (Fup== 0) 
  {
    skip <- TRUE
  } else skip <- FALSE

  if (!skip)
  {  
    p <- parameterize_fetal_pbtk(dtxsid=this.id,
    fetal_fup_adjustment =FALSE,
    suppress.messages=TRUE))
    out <- suppressWarnings(solve_fetal_pbtk(
      parameters=p,
      dose=0,
      times=times,
      daily.dose=1,
      doses.per.day=3,
      output.units = "uM",
      suppress.messages=TRUE)
    last.row <- which(out[,"time"]>279)
    last.row <- last.row[!duplicated(out[last.row,"time"])]
    new.row <- data.frame(
      Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
      DTXSID = this.id,
      Mat.pred = mean(out[last.row,"Cplasma"]),
      Fet.pred = mean(out[last.row,"Cfplasma"]),
      MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
      )
    MFratio.pred.nofup <- rbind(MFratio.pred.nofup,new.row)
  }
  }

Supplemental Histogram of maternal-fetal ratio without Fup correction

FigSupD <- ggplot(data=MFratio.pred.nofup)+
   geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+ 
  xlab("Maternal:Fetal Plasma Concentration Ratio (No Fup Corr.)") +
  ylab("Number of chemicals")+
    theme_bw()+
    theme( text  = element_text(size=14))


print(FigSupD)

Statistics on maternal-fetal ratio for full HTTK library without Fup correction

max.chem <- MFratio.pred.nofup[which(
  MFratio.pred.nofup$MFratio.pred==max(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
min.chem <- MFratio.pred.nofup[which(
  MFratio.pred.nofup$MFratio.pred==min(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
cat(paste("In Figure X we examine the ratios predicted for the ",
  dim(MFratio.pred)[1],
  " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
  sep=""))


cat(paste("We observe a median ratio of ",
  signif(median(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE),3),
  ", with values ranging from ",
  signif(min.chem[,"MFratio.pred"],3),
  " for ",
  min.chem[,"DTXSID"],
  " to ",
  signif(max.chem[,"MFratio.pred"],3),
  " for ",
  max.chem[,"DTXSID"],
  ".\n",sep=""))

# Check out phys-chem > 1.6, < 1:
highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred>1.6)$DTXSID)

cat(paste("There are",
  dim(highratio)[1],
  "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))

# all highly bound
highratio$Compound
suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=2)))


lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred<0.9)$DTXSID)
# No obvious pattern

Code used to create data distributed with vignette

aylward2014 <-MFdata
pregnonpregaucs <- TKstats
fetalpcs <- Curley.pcs
wang2018 <- Wangchems

save(aylward2014,pregnonpregaucs,fetalpcs,wang2018,pksim.pcs,
     file="Kapraun2022Vignette.RData",version=2)

httk documentation built on March 7, 2023, 7:26 p.m.