Content

1 Uncertainty and Variability

2 Monte Carlo simulation - Prediction

3 Markov Chain Monte Carlo - Calibration

4 Hands-on Exercise

Deterministic vs Probabilistic

Traditional - Deterministic

Deterministic vs Probabilistic

Modeling in Risk Assessment

https://www.epa.gov/risk/about-risk-assessment#whatisrisk

Uncertainty vs. Variability

Variability in Risk Assessment

• Reduce chances using a strain that is a "poor" model of humans
• Obtaining information about "potential range" to inform risk assessment

Variability & Uncertainty

Monte Carlo Simulation

Uncertainty in Risk Analysis

Two major types of uncertainty need to be differentiated:

(1) Uncertainty due to physical variability

(2) Uncertainty due to lack of knowledge in

Modeling uncertainty

Deterministic Simulation

• Define exposure unit & calculate point estimate

1-D Monte Carlo Simulation: Uncertainty

• Identify probability distributions to simulate probabilistic outputs

2-D Monte Carlo Simulation: Uncertainty & Variability

• Bayesian statistics to characterize population uncertainty and variability

Uncertainty in parameter

• The parameter is an element of a system that determine the model output.

• Parameter uncertainty comes from the model parameters that are inputs to the mathematical model but whose exact values are unknown and cannot be controlled in physical experiments.

$$y = f(x_i) $$

Mayer J, Khairy K, Howard J. Drawing an elephant with four complex parameters. Am. J. Phys. 78, 648 (2010) https://doi.org/10.1119/1.3254017

Uncertainty in PBPK model parameter

Simulation in GNU MCSim

Monte Carlo simulations

• Perform repeated (stochastic) simulations across a randomly sampled region of the model parameter space.

Used to: Check possible simulation (under given parameter distributions) results before model calibration

SetPoints simulation

• Solves the model for a series of specified parameter sets. You can create these parameter sets yourself or use the output of a previous Monte Carlo or MCMC simulation.

Used to: Posterior predictive check, Local/global sensitivity analysis

Currently, the Bayesian Markov chain Monte Carlo (MCMC) algorithm is an effective way to do population PBPK model calibration.

It is a powerful tool, Because...

Frequentist vs. Bayesian

Bayes' rule

\(y\): Observed data

\(\theta\): Observed or unobserved parameter

\(p(\theta)\): Prior distribution of model parameter

\(p(y|\theta)\): Likelihood of the experiment data given by a parameter vector

\(p(\theta|y)\): Posterior distribution

\(p(y)\): Likelihood of data

Frequentist vs. Bayesian

Probabilistic Modeling

https://doi.org/10.1016/j.ddmod.2017.08.001

Markov Chain Monte Carlo

The algorithm was named for Nicholas Metropolis (physicist) and Wilfred Keith Hastings (statistician). The algorithm proceeds as follows.

Initialize

1. Pick an initial parameter sets \(\theta_{t=0} = \{\theta_1, \theta_2, ... \theta_n\}\)

Iterate

1. Generate: randomly generate a candidate parameter state \(\theta^\prime\) for the next sample by picking from the conditional distribution \(J(\theta^\prime|\theta_t)\)
2. Compute: compute the acceptance probability \(A\left(\theta^{\prime}, \theta_{t}\right)=\min \left(1, \frac{P\left(\theta^{\prime}\right)}{P\left(\theta_{t}\right)} \frac{J\left(\theta_{t} | \theta^{\prime}\right)}{J\left(\theta^{\prime} | \theta_{t}\right)}\right)\)
3. Accept or Reject:
   1. generate a uniform random number \(u \in[0,1]\)
   2. if \(u \leq A\left(x^{\prime}, x_{t}\right)\) accept the new state and set \(\theta_{t+1}=\theta^{\prime}\), otherwise reject the new state, and copy the old state forward \(\theta_{t+1}=\theta_{t}\)

Simulation in GNU MCSim

Markov-chain Monte Carlo (MCMC) simulation

• Performs a series of simulations along a Markov chain in the model parameter space.
• They can be used to obtain the Bayesian posterior distribution of the model parameters, given a statistical model, prior parameter distributions and data for which a likelihood function can be computed. Used to Model calibration

Calibration & evaluation

Prepare model and input files

• Need at least 4 chains in simulation

Check convergence & graph the output result

• Parameter, log-likelihood of data
• Trace plot, density plot, correlation matrix, auto-correlation, running mean, ...
• Gelman–Rubin convergence diagnostics

Evaluate the model fit

• Global evaluation
• Individual evaluation

Example - Linear model

## linear.model.R ####
Outputs = {y}
# Model Parameters
A = 0; # 
B = 1;
CalcOutputs { y = A + B * t); }
End.

## linear_mcmc.in.R ####
MCMC ("MCMC.default.out","", # name of output 
     "",           # name of data file
     2000,0,       # iterations, print predictions flag,
     1,2000,       # printing frequency, iters to print
     10101010);    # random seed (default )
Level {
  Distrib(A, Normal, 0, 2); # prior of intercept 
  Distrib(B, Normal, 1, 2); # prior of slope 
  Likelihood(y, Normal, Prediction(y), 0.05);
  Simulation {
    PrintStep (y, 0, 10, 1); 
    Data  (y, 0.01, 0.15, 2.32, 4.33, 4.61, 6.68, 
                7.89, 7.13, 7.27, 9.4, 10.0);
  }
}
End.

Example - MCMC simulation

model <- "models/linear.model"
input <- "inputs/linear.mcmc.in"
set.seed(1111) 
out <- mcsim(model, input)

head(out)

##   iter    A.1.      B.1.    LnPrior    LnData LnPosterior
## 1    0 5.17187 -0.421849  -6.820405 -591.1577   -597.9781
## 2    1 5.17187 -0.421849  -6.820405 -591.1577   -597.9781
## 3    2 5.43465 -0.421849  -7.168811 -565.2515   -572.4203
## 4    3 8.62813 -0.421849 -12.782460 -493.2532   -506.0356
## 5    4 7.97506 -0.421849 -11.427070 -471.4773   -482.9044
## 6    5 7.51347 -0.421849 -10.533410 -467.4057   -477.9391

tail(out, 4)

##      iter     A.1.     B.1.   LnPrior    LnData LnPosterior
## 1998 1997 0.706138 0.957902 -3.286722 -19.06480   -22.35153
## 1999 1998 0.706138 0.957902 -3.286722 -19.06480   -22.35153
## 2000 1999 0.706138 0.974017 -3.286585 -19.13584   -22.42242
## 2001 2000 0.462481 0.974017 -3.250992 -18.92306   -22.17405

plot(out$A.1., type = "l", xlab = "Iteration", ylab = "")
lines(out$B.1., col = 2)
legend("topright", legend = c("Intercept", "Slope"), 
       col = c(1,2), lty = 1)

Example - Posterior check

plot(out$A.1., out$B.1., type = "b",
     xlab = "Intercept", ylab = "Slope")

str <- ceiling(nrow(out)/2) + 1
end <- nrow(out)
j <- c(str:end)
plot(out$A.1.[j], out$B.1.[j], type = "b",
     xlab = "Intercept", ylab = "Slope")

Example - Posterior Predictive Checks

out$A.1.[j] %>% density() %>% plot(main = "Intercept")
plot.rng <- seq(par("usr")[1], par("usr")[2], length.out = 1000)
prior.dens <- do.call("dnorm", c(list(x=plot.rng), c(0,2)))
lines(plot.rng, prior.dens, lty="dashed")

out$B.1.[j] %>% density() %>% plot(main = "Slope")
plot.rng <- seq(par("usr")[1], par("usr")[2], length.out = 1000)
prior.dens <- do.call("dnorm", c(list(x=plot.rng), c(1,2)))
lines(plot.rng, prior.dens, lty="dashed")

Example - Evaluation of prediction

# Observed
x <- seq(0,10,1)
y <- c(0.0, 0.15, 2.32, 4.33, 4.61, 6.68, 
       7.89, 7.13, 7.27, 9.4, 10.0)
# Expected
dim.x <- ncol(out)
for(i in 1:11){
  out[,ncol(out)+1] <- out$A.1. + out$B.1.*x[i]
}
# Plot
plot(x, y, pch ="")
for(i in 1901:2000){
  lines(x, out[i,c((dim.x+1):ncol(out))],col="grey")
}
points(x, y)

Example - Evaluation of prediction

# Use prediction from the last iteration
i <- 2000
Expected <- out[i, c((dim.x+1):ncol(out))] %>% 
  as.numeric()
Observed <- c(0.01, 0.15, 2.32, 4.33, 4.61, 6.68, 
              7.89, 7.13, 7.27, 9.4, 10.0)
# Plot
plot(Expected, Observed/Expected, pch='x')
abline(1,0)

General Bayesian-PBPK Workflow

1 Model constructing or translating

2 Verify modeling result

• Compare with published result
• Mass balance

3 Uncertainty (and sensitivity) analysis

4 Model calibration and validation

• Markov chain Monte Carlo
  • Diagnostics (Goodness-of-fit, convergence)

Summary

Hands on Exercise

Task 1. Uncertainty analysis on PK model (code: https://rpubs.com/Nanhung/SRP19_6)

• Before model calibration, we need to learn how to conduct Monte Carlo simulation to set the proper parameter distribution

Task 2. Model calibration (code: https://rpubs.com/Nanhung/SRP19_7)

• After the uncertainty analysis, we can calibrate the model parameters by Markov chain Monte Carlo technique

Task 3. Monte Carlo Simulation for PBPK model (code: https://rpubs.com/Nanhung/SRP19_8)

• Learn how parameter effect on model variable in PBPK model

Task 4. PBPK model in MCSim (code: https://rpubs.com/Nanhung/SRP19_9)

• Sometimes, the simulation process in R is very computational expensive. We need to solve it.

Hands on Exercise

Task 1: Uncertainty analysis on PK model

Hands on Exercise

Task 2: Model calibration

Hands on Exercise

Task 3: Monte Carlo Simulation for PBPK model

• Reproduce the published Monte Carlo analysis result in Bois and Brochot (2016)*
  • Testing parameter include body mass (BDM), pulmonary flow (Flow_pul), partition coefficient of arterial blood (PC_art) and metabolism rate (Kmetwp)
• Construct the relationship between body mass and quantity in fat

Hands on Exercise

Task 4: Monte Carlo Simulation for PBPK model (MCSim)

• The Monte Carlo Simulation take a little bit longer with ode function in deSolve package.
• Therefore we want to improve the computational speed. Now, rewrite the R model code to MCSim and conduct Monte Carlo Simulation with the same parameter setting.
• The goal of this exercise is to compare the computational time and output (MCSim vs. R).

Distrib (BDM, Normal, 73, 7.3);
Distrib (Flow_pul, Normal, 5, 0.5);
Distrib (PC_art, Normal, 2, 0.2);
Distrib (Kmetwp, Normal, 0.25, 0.025);