---
title: "Kapraun et al. (2022): Generic Human Gestational Model"
author: "Dustin F. Kapraun, Mark Sfeir, Robert Pearce, Sarah Davidson, Annie Lumen, Andre Dallmann, Richard Judson, and John F. Wambaugh"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{f) Kapraun (2022): Human Gestational Model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
*Please send questions to wambaugh.john@epa.gov*

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.
```{r configure_knitr, eval = TRUE}
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
options(rmarkdown.html_vignette.check_title = FALSE)
```

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

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

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

From the R command prompt:

install.packages("X")

Or, if using RStudio, look for 'Install Packages' under 'Tools' tab.

```{r load_packages, eval = execute.vignette}
#
#
# Setup
#
#
#library(readxl)
library(ggplot2)
library(httk)
library(scales)
library(gridExtra)
library(cowplot)
```

### Prepare custom functions
```{r load_code, eval = execute.vignette}
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:
```{r chem_numbers, eval=FALSE}
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](http://dx.doi.org/10.1080/10937404.2014.884956) 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
```{r load_aylward_data, eval = execute.vignette}
#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:
```{r process_aylward_data, eval = execute.vignette}
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.
 
```{r make_aylward_preds, eval = execute.vignette}
# 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:
```{r cleanup_repeated_aylward_figure, eval = execute.vignette}
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:
```{r aylward_figure_no_fup, eval = execute.vignette}
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:
```{r aylward_analysis_text, eval = execute.vignette}

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:
```{r aylward_figure_with_fup, eval = execute.vignette}
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:
```{r voc_analysis, eval = execute.vignette}
# 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](http://dx.doi.org/10.1007/s40262-017-0594-5) 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. 

```{r dallmann2018_data, eval = execute.vignette}
#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.

```{r dallmann2018_preds, eval = execute.vignette}
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
```{r dallmann2018_figure, eval = execute.vignette}
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](http://dx.doi.org/10.1016/j.tiv.2007.09.010) as described by [Pearce et al., 2017](http://dx.doi.org/10.1007/s10928-017-9548-7). 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](https://doi.org/10.1080/00039896.1969.10666901) 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. 

```{r curley1969_data, eval = execute.vignette}
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.

```{r curley1969_calcpcs, eval = execute.vignette}
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](http://dx.doi.org/10.1007/s00204-002-0339-5) for Bisphenol A and Diadzen. 

```{r csanady2002_pcs, eval = execute.vignette}
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](https://doi.org/10.1080/00039896.1969.10666901) were modeled by [Weijs et al., 2013](http://dx.doi.org/10.1021/es400386a) using the same partition coefficients for mother and fetus. The values used represented "prior knowledge" summarizing the available literature.


```{r weijs2013_loadPCs, eval = execute.vignette}
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
```{r curley1969_makepreds, eval = execute.vignette}
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:
```{r curley1969_figure_compounds, eval = execute.vignette}
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:
```{r curley1969_figure_tissues, eval = execute.vignette}
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:
```{r read_pksim_pcs, eval = execute.vignette}
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"])

```
```{r compare_pksim_pcs, eval = execute.vignette}
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
```{r wang2018_loaddata, eval = execute.vignette}
#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
```{r wang2018_makepreds, eval = execute.vignette}

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:
```{r brain_uncertainty_prop_table, eval = execute.vignette}
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):
```{r Rfetmat_uncertainty, eval = execute.vignette}
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
```{r Rfetbrain, eval = execute.vignette}
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
```{r Rfetbrain_uncertainty, eval = execute.vignette}

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
```{r Wang_Uncertainty_Propagation, eval = execute.vignette}

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:
```{r wang2018_figure, eval = execute.vignette}
#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:
```{r wang2018_details, eval = execute.vignette}

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:
```{r impact_of_model, eval = execute.vignette}
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
```{r wang2018_makepreds_nofup, eval = execute.vignette}

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:
```{r brain_uncertainty_prop_table_nofup, eval = execute.vignette}
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):
```{r Rfetmat_uncertainty_nofup, eval = execute.vignette}
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
```{r Rfetbrain_nofup, eval = execute.vignette}
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
```{r Rfetbrain_uncertainty_nofup, eval = execute.vignette}

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


```{r compare_fup_correction_case_study,eval=FALSE}
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](http://dx.doi.org/10.1208/ps040104),

$$
	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](http://dx.doi.org/10.1208/ps040104)). 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. 

```{r fup_table, eval = execute.vignette}
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:
```{r fup_figure, eval = execute.vignette}
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:
```{r MFratio_allhttk_preds, eval = execute.vignette}
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
```{r MFratio_allhttk_figure, eval = execute.vignette}
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
```{r MFratio_allhttk_stats, eval = execute.vignette}
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:
```{r MFratio_allhttk_preds_nofup, eval = execute.vignette}
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
```{r MFratio_allhttk_figure_nofup, eval = execute.vignette}
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
```{r MFratio_allhttk_stats_nofup, eval = execute.vignette}
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
```{r create_distributable_data, eval = execute.vignette}
aylward2014 <-MFdata
pregnonpregaucs <- TKstats
fetalpcs <- Curley.pcs
wang2018 <- Wangchems

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