Kapraun et al. (submitted): Generic Human Gestational Model

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

2022-09-08

wambaugh.john@epa.gov

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.

Prepare for session

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

Load the relevant libraries

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.

Prepare custom functions

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:

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:

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)

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

Compare httk partition coefficients with PKsim:

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:

Make prioritization figure:

Look at ionization state:

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

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:

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"]
}

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:

Maternal-Fetal Predictions across HTTK Chemical Library:

Histogram of maternal-fetal ratio

Statistics on maternal-fetal ratio for full HTTK library

Examine Maternal-Fetal ratio Predictions without fetaul Fup correction:

Supplemental Histogram of maternal-fetal ratio without Fup correction

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

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)