Building your first PBPK model:the basics

In previous posts we referred exclusively to modelling using the top-down, population approach. However, in recent years, physiology-based, bottom-up approaches are getting more attention from both industry and regulators.

Population and physiology-based approaches share some common ground: both approaches describe the body as a system of compartments connected by rates to describe drug disposition. If, population approaches use 1, 2 and rarely 3 or 4 compartments to describe the PK of a drug, the PBPK approaches have at least one compartment for every major organ.

The data required as input for building each type of model is different. For example, a population PK model is based on concentration-time data, whereas, PBPK models are based on more diverse data:

  • drug-specific parameters (e.g., molecular weight, pKa, lipophylicity, etc.)
  • system-specific parameters (e.g., tissue volumes and blood perfusion, lipid, water and protein content, etc.)
  • drug biological properties, that are a combination between a drug and a system parameter (e.g. active processes dependent on the number of enzymes in the tissue and the affinity of the substrate to that specific enzyme)

In this case, all these parameters are required to predict the PK profile of the drug of interest.

Each organ is described either as perfusion-rate limited or permeability-rate limited. For perfusion limited tissues, the distribution of the drug into the tissue and the time to reach steady-state are limited only by the blood flow through the tissue. At steady-state the drug concentration in the tissue is in equilibrium with the concentration in the circulation. In this case, the equilibrium is reached instantly and is strictly dependent on tissue volume, tissue perfusion and the tissue to plasma partitioning coefficient (Kp). For permeability limited tissues, time to reach steady-state is limited by the permeability of the drug across membranes and/or by active processes (i.e., transporters). The membrane divides the tissue into extra- and intra-cellular space. Capturing the distribution into these tissues is more advanced and out of the scope of this post.

In addition, PBPK organs can be either non-eliminating (e.g. adipose tissue) or eliminating (e.g., liver where metabolism can occur). The elimination of the drug from the eliminating organs includes a “rate of elimination ” in addition to blood flow only.

The whole PBPK model in the mrgsolve library incorporates all these concepts and more. This was originally published by Jones and Rowland-Yeo and the code in Berkeley-Madonna is available as a supplement. I’ll quickly go through the code and give more explanations alongside it.

You can find the model code under R -> library -> mrgsolve -> models ->pbpk.cpp and open it in a text editor (like notepad).

First, you can see the problem block [PROB] where it is stated that we are looking at a “Human PBPK Model” and the original reference for it. Also, the default settings for the simulation are given under [SET]: it ends after 24h (end = 24) with one data point simulated every 15 minutes (delta = 0.25).

[ PROB ]


1: Jones H, Rowland-Yeo K. Basic concepts in physiologically based pharmacokinetic modeling in drug discovery and development. CPT Pharmacometrics Syst Pharmacol. 2013 Aug 14;2:e63. doi: 0.1038/psp.2013.41. PubMed PMID: 23945604; PubMed Central PMCID: PMC3828005.

[ SET ] end = 24, delta = 0.25

Next, you define the required compartments: one for dosing and one compartment for every major organ/tissue that make up the PBPK model:

[ CMT ]
D //; dose
Aad //; adipose
Abo //; bone
Abr //; brain
Agu //; gut
Ahe //; heart
Aki //; kidney
Ali //; liver
Alu //; lung
Amu //; muscle
Ask //; skin
Asp //; spleen
Ate //; testes
Ave //; venous blood
Aar //; arterial blood
Are //; rest of body

You should notice that the last compartment defines the rest of the body. This can be used as an “additional organ” that we are interested in and is not represented. For example, it could represent a cancerous tumour when we are interested in the distribution of the chemotherapeutic agent inside the tumour.

In the [PARAM] block we start defining the fractions of the whole body weight that will be needed later to obtain the volumes of each tissue. To keep this post short, I will remove the values in between and mark this action with […].

BW = 70 //; BW (kg)
// {Fractional tissue volumes}
FVad = 0.213 //; adipose
FVbo = 0.085629 //; bone
FVpl = 0.0424 //; plasma
FVrb = 0.0347 //; erythrocytes
FVre = 0.099771 //; rest of body

Next, we need the fractions of the cardiac output to obtain the blood flows through each organ. All these parameters have been published and are relatively easy to find in literature for humans, as well as for a few animal models.

Note: There is no fraction of cardiac output to define the blood flow through plasma or erythrocytes.

// {Fractional tissue blood flows}
FQad = 0.05 //; adipose
FQbo = 0.05 //; bone
FQre = 0.103855 //; rest of body

These fractions define the main system-specific parameters required by the PBPK model. This model is based only on perfusion rate limited organs which means that their kinetics are dependent on tissue volumes, blood flows and the tissue-to-plasma partitioning coefficients (Kp) for each organ. Kp values are defined below, under {Compound specific parameters}. These coefficients can be obtained in silico and are highly dependent on fraction unbound, lipophylicity and pKa of the compound as well as the water/lipid/protein composition of the tissue. More information about how to obtain Kp values were published here.

// {Tissue to plasma partition coefficients}
Kpad = 0.191 //; adipose
Kpbo = 0.374 //; bone
Kpre = 0.600 //; rest of body

The fraction of unbound drug in plasma, the blood-to-plasma ratio and the fraction unbound in microsomes (to be discussed together with in vitro-in vivo extrapolation of clearance) are measured in vitro and added to the model. These parameters determine the concentration of drug available in plasma for passive diffusion, metabolism, transport or excretion.

// {In vitro binding data}
fup = 0.681 //; fraction unbound in plasma
BP = 0.98 //; blood to plasma ratio
fumic = 1 //; fraction unbound in microsomes

In this case, drug clearance was derived from in vitro experiments with human liver microsomes (HLM). This can give a good estimate of the CYP450-mediated metabolism in the liver. This drug is neither metabolised nor excreted unchanged by the kidneys (CLrenal = 0).

// {Clearances}
HLM_CLint = 8 //; HLM CLint apparent (ul/min/mg)
CLrenal = 0 //; CLint renal (L/hr)

In this example, the drug was administrated orally, for which absorption rate (Ka) and bioavailability (F) are given as parameters. Alternatively, absorption can be described as a first-order absorption model, which requires effective permeability and the intestinal surface area available for drug absorption. Effective permeability can be measured in vivo with a jejunal perfusion system or in vitro in cell cultures (e.g., Caco-2 cells, MDCK, LLC-PK, etc.) and then scaled accordingly.

// {Absorption}
Ka = 2.18 //; Ka (hr-1)
F = 1.00 //; fraction absorbed
CO = 108.33 //; cardiac output (ml/s)

In the [MAIN] block tissue volumes (L), blood flows (L/h) and venous and arterial plasma volumes are calculated. Notice that the blood flow from the artery to the liver is dependent on the blood flow through the gut and the spleen.

[ MAIN ]

// {Total tissue volumes – L}
double Vad = BW*FVad; // adipose

double Vbo = BW*FVbo; // bone


double Vre = BW*FVre; // rest of body

double Vplas_ven = Vpl*Vve/(Vve + Var) ; // venous plasma

double Vplas_art = Vpl*Var/(Vve + Var) ; // arterial plasma

// {Total tissue blood flows – L/hr}
double QC = CO/1000*60*60 ; // cardiac output (L/hr)

double Qad = QCFQad ; // adipose

double Qbo = QCFQbo ; // bone


double Qh = QC*FQh ; // hepatic (venous side)
double Qha = Qh – Qgu – Qsp; // hepatic artery


double Qre = QC*FQre ; // rest of body

Under [ODE] we calculate the concentrations in each organ from the amount and the organ volume and we derive the unbound concentrations in the liver and kidney tissues, by multiplying the total concentration with the fraction unbound in plasma (as measured in vitro).


double Cadipose = Aad/Vad; // adipose

double Cbone = Abo/Vbo; // bone


// {Calculation of free concentrations – mg/L}
double Cliverfree = Cliverfup; // liver

double Ckidneyfree = Ckidneyfup; // kidney

Next, the in vitro measurements from the HLM system are extrapolated to the corresponding in vivo value by in vitro-in vivo extrapolation (IVIVE). The metabolic clearance in the liver for this compound is extrapolated from the measured in vitro intrinsic clearance to the whole organ intrinsic clearance that serves as input in the PBPK model. Intrinsic clearance (CLint) is defined as the intrinsic capacity of a compound to be metabolised by relevant enzymes without the influence of fraction unbound and blood flow through the organ. CLint is a drug-biological parameter as it is influenced by
both drug properties (i.e., substrate specificity to CYP enzymes as derived from the in vitro HLM experiment) and system properties (i.e., mg of microsomal proteins per gram of liver (MPPGL) and the volume of the liver (LV). Fraction unbound in microsomes, as measured in the in vitro system (fumic), MPPGL and LV are used as scalars for IVIVE. The methods used to perform IVIVE are dependent on the in vitro system, as different scalars are required for different set-ups. Also, something one should always pay attention to when performing IVIVE are the units. The units for the scaled parameter should be consistent with the units of the rest of the parameters making up the PBPK model.

/ {Clearance calculations}
double MPPGL = 45; // mg microsomal protein per g liver
double CLmet = (HLM_CLint/fumic)*MPPGL*Vli*60/1000; // CLint scaled (L/hr)

Once all parameters are derived it’s time to write up the differential equations that make up the PBPK model. The blood flow through the organs goes from the arterial compartment to the venous compartment, except for the lung. The drug concentration that ends-up in the venous compartment is defined separately to make the rest of the code easier to read, as it sums up all the concentrations going out of all the compartments. The only rate going out of the venous compartment (in this case) is the blood flow to the lung.

dxdt_Ave = Venous – Qlu*Cvenous; // venous blood

Oral absorption follows the classic Ka*Dose*F and goes into the gut compartment. Even if F was set to 1 in the {Absorption} section of the [PARAM] block, the authors decided to remove it when defining the Absorption parameter. To change the administration from oral to IV bolus, you only need to add the Absorption into the Venous compartment and remove Ka. Don’t forget to remove the Absorption from the gut compartment.

The most common collected sample is plasma. Therefore we are interested in capturing that as an output in a table.

[ CAPTURE ] Cvenous = Ave/Vve
capture Cp = Cvenous/BP ; // venous plasma

The venous concentration represents the drug concentration in venous blood, therefore, the blood-to-plasma ratio (BP) is included in all differential equations to ensure that the tissue concentrations are related to blood concentrations. To obtain the plasma concentrations we divide the blood venous concentration by BP. In short, BP relates drug concentrations in plasma and blood at equilibrium and it is dependent on the Kp for blood cells (drug-specific), the hematocrit (system-specific) that is used to calculate the volume of blood cells and fraction unbound in plasma.

The advantage of having this PBPK model in the library is that you can easily change it and use it for simulations in R. You just need to load the mrgsolve library and get started.

# — read in pbpk model from library:
mod <- mread(“pbpk”, modlib())
# — simulate from basic pbpk model:
mod %>%
ev(amt = 100) %>%
mrgsim %>%

Basic simulation with the PBPK model showing plasma concentration vs. time for a 70 kg individual.

I quickly went through the tutorial presented on mrgsolve intro page and I saw that it is quite easy to do a local sensitivity analysis – just by changing one parameter at a time and assessing the influence of that change on your output, in this case Cp.

# — How to simulate for a range of parameter values (local sensitivity analysis)
idata <- expand.idata(Kpli = seq(4,20,2))
mod %>%
ev(amt = 100) %>%
idata_set(idata) %>%
mrgsim %>%

Drug plasma concentrations over time obtained for a 70 kg individual by varying kp,liver between 4 and 20.

In previous blog posts, it was shown how to quickly change the dosing event, parameters values on the fly. All of that is applicable for the PBPK model too.


THIS POST WAS WRITTEN BY Sinziana Cristea FROM Leiden University.
She works on predictive models to study maturation of renal elimination processes, ultimately to guide dosing of renally excreted drugs in children.

Leave a Reply

Your email address will not be published. Required fields are marked *