Exercise 2.4: Monte Carlo Simulation for PBPK model (MCSim)

Task

• 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).

Solution

Package

library(simuloR)

## The package is powered by GNU MCSim v6.1.0.

1,3-butadiene PBPK Model

States  = {Q_fat,             # Butadiene quantity in the fat        (mg)
  Q_wp,              # ~         in well-perfused tissues   (mg)
  Q_pp,
  Q_met};             # ~         in poorly-perfused tissues (mg)


Outputs = {C_alv,             # Butadiene concentration in alveolar air   (mg/l)
  C_exh,             # ~                          exhaled        (ppm)
  C_art,             # ~                          arterial blood (mg/l)
  C_ven,             # ~                          venous blood   (mg/l)
  dQ_ven,            # 
  PC_art_obs,        # Measured blood/air partition coefficient
  Flow_pul_obs};     # ~        pulmonary flow

Inputs  = {C_inh};            # Butadiene concentration inhaled (ppm)


# Constants
# =========
# Molecular weight (in grams)
MW_bu = 54.0914;

# Conversions between ppm and mM (under normal conditions):
ppm_per_mM = 24450;

# Conversions from/to ppm: 24450 ppm = 54.0914 mg/l
ppm_per_mg_per_l = 24450 / 54.0914;
mg_per_l_per_ppm = 1 / ppm_per_mg_per_l;


# Default parameter values
# ========================
# Units:
# Volumes: liter
# Time:    minute
# Flows:   liter / minute

# Physiological and pharmacokinetic parameters
# --------------------------------------------
# Bodyweight (kg)
BDM    = 73;  # standard value
height = 1.6; # meters
Age    = 40;
Sex    = 1;
N_BR   = 12; # nbr of breathes per minute

# Tissue volumes as % of BDM (density 1 assumed)
Pct_BDM_wp  = 0.2;  # well perfused tissue, % of lean mass

# Pulmonary ventilation rate (minute volume) (L/min)
Flow_pul = 5;

Pct_Deadspace = 0.7; # Percent pulmonary deadspace

Pct_Flow_fat = 0.1;
Pct_Flow_pp  = 0.35;

# Blood/air partition coefficient
PC_art = 2;

# Tissue/blood partition coefficients
PC_fat = 22;
PC_wp  = 0.8;
PC_pp  = 0.8;

# Ventilation over perfusion ratio
Vent_Perf = 1.14;

# Hepatic rate constant for metabolism (1/min)
Kmetwp = 0.25;

# Scaled parameters
# -----------------
# The following parameters are calculated from the above values in the
# Scale section before the start of each simulation.
# They are left uninitialized here.

# Total blood flow
Flow_tot = 0;
Flow_alv = 0;

Eff_V_fat = 0;
Eff_V_wp  = 0;
Eff_V_pp  = 0;

Flow_fat = 0;
Flow_wp  = 0;
Flow_pp  = 0;


Initialize {
  
  # Calculate Flow_alv from total pulmonary flow
  Flow_alv = Flow_pul * (1 - Pct_Deadspace);
  
  # Calculate total blood flow from the alveolar ventilation rate and the
  # V/P ratio.
  Flow_tot = Flow_alv / Vent_Perf;
  
  Pct_BDM_fat = (1.2 * BDM / (height * height) - 10.8 *(2 - Sex) + 
                   0.23 * Age - 5.4) * 0.01;
  # 1 male, 2 female
  Eff_V_fat = Pct_BDM_fat * BDM;
  Eff_V_wp  = Pct_BDM_wp  * BDM * (1 - Pct_BDM_fat);
  Eff_V_pp  = 0.9 * BDM - Eff_V_fat - Eff_V_wp;
  
  # Calculate actual blood flows from total flow and percent flows
  Flow_fat = Pct_Flow_fat * Flow_tot;
  Flow_pp  = Pct_Flow_pp  * Flow_tot;
  Flow_wp  = Flow_tot * (1 - Pct_Flow_pp - Pct_Flow_fat);
  
} # end of Initialize


Dynamics {
  
  # Venous blood concentrations at the organ exit
  Cout_fat = Q_fat / (Eff_V_fat * PC_fat);
  Cout_wp  = Q_wp  / (Eff_V_wp  * PC_wp );
  Cout_pp  = Q_pp  / (Eff_V_pp  * PC_pp );
  
  # Sum of Flow * Concentration for all compartments
  dQ_ven = Flow_fat * Cout_fat + Flow_wp * Cout_wp + Flow_pp * Cout_pp;
  
  # Arterial blood concentration
  # Convert input given in ppm to mg/l to match other units
  C_art = (Flow_alv * C_inh / ppm_per_mg_per_l + dQ_ven) /
    (Flow_tot + Flow_alv / PC_art);
  
  # Quantity metabolized in liver
  dQmet_wp = Kmetwp * Q_wp;
  
  # Differentials for quantities
  dt (Q_fat) = Flow_fat * (C_art - Cout_fat);
  dt (Q_wp)  = Flow_wp  * (C_art - Cout_wp) - dQmet_wp;
  dt (Q_pp)  = Flow_pp  * (C_art - Cout_pp);
  dt (Q_met) = dQmet_wp;
  
} # End of Dynamics


CalcOutputs {
  
  # Venous blood concentration
  C_ven = dQ_ven / Flow_tot; 
  C_alv = C_art / PC_art;
  C_exh = (C_alv <= 0 ?
             10E-30     :
             (1 - Pct_Deadspace) * C_alv * ppm_per_mg_per_l + 
             Pct_Deadspace * C_inh);
  
  Flow_pul_obs = Flow_pul / pow(BDM, 0.7);
  PC_art_obs   = PC_art;
  
} # End of CalcOutputs

End.

Inputs

# bd.mtc.in -------------------------

Integrate (Lsodes, 1e-6, 1e-6, 1);

MonteCarlo ("simmc.out", 
            1000, 
            10101010);

# Distribution
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);

Simulation {
  C_inh = NDoses (2, 10, 0, 0, 120);
  Print(Q_fat, 1440);  
}
End.

Modeling

model <- "models/bd.model"
makemcsim(model)

## * Creating executable program, pleas wait...

## * Created executable program 'models/mcsim.bd.model'.

start_time <- Sys.time()
set.seed(1)
bd.mtc.out <- mcsim(model = model, input = "inputs/bd.mtc.in", compile = F)

## Executing...

end_time <- Sys.time()

end_time - start_time

## Time difference of 0.1008842 secs

Output analysis

head(bd.mtc.out)

##   Iter     BDM Flow_pul  PC_art   Kmetwp Q_fat_1.1
## 1    0 82.1029  4.62478 2.27013 0.313347  0.270567
## 2    1 65.9625  5.34661 1.92688 0.251279  0.236769
## 3    2 75.7812  5.47253 2.31939 0.248192  0.302859
## 4    3 66.0072  4.78422 1.95573 0.248159  0.223215
## 5    4 75.1593  4.15566 2.15029 0.238767  0.235095
## 6    5 74.8441  4.46096 1.51331 0.294807  0.199929

BDM <- bd.mtc.out$BDM
hist(BDM)

Q_fat <- bd.mtc.out$Q_fat_1.1
hist(Q_fat)

plot(BDM, Q_fat, xlab="Body mass (kg)", ylab="Quantity in Fat (mg)")