INTRODUCTION

Physiologically based pharmacokinetic (PBPK) models are increasingly used in drug development for their ability to predict drug concentrations in clinically relevant tissues (e.g. target site of the drug) and to extrapolate across species, experimental conditions (e.g. co-administration with a metabolic inductor/inhibitor) and special populations (e.g. children). Since PBPK models are mechanistic in nature, parameters can be defined a piori using knowledge from the anatomy/physiology literature (system-specific parameters) and preclinical experiments (drug-specific parameters), and can thus be used for exploratory simulations. Nevertheless, it is often desirable to fit the model to animal or clinical data to improve the estimates of the drug-specific parameters. A whole-body PBPK (WBPBPK) model is typically expressed as a system of ordinary differential equations (ODEs), in which each state represents a tissue or an organ space. During clinical studies, tissue concentrations are usually not measured due to ethical constraints. As a consequence, fitting a WBPBPK model to clinical data commonly yields a numerically unstable analysis as many parameters (e.g. tissue-to-blood partition coefficients) cannot be informed just by plasma data. To stabilise the analysis, prior knowledge on the parameter values can be incorporated in the model using a Bayesian approach. Such an approach also allows biological/physiological constraints to be maintained when estimating the parameters. A Bayesian population analysis produces posterior distributions of both the population and individual parameters, which can be approximated using Markov Chain Monte Carlo (MCMC) simulation (1). Optimising complex systems, such as WBPBPK models, using MCMC techniques can be time consuming, especially when analysing a large amount of sparse and noisy data typically collected during clinical trials (2). In addition, depending on the informativeness of the priors assigned to the model parameters, the statistical model might not be identifiable for a high-dimensional system (3, 4).

The dimension of a dynamical system can be reduced using model order reduction methods (5). Lumping is one of these methods that allows all states of the original system to be transformed to fewer pseudo-states, thus reducing the number of ODEs and parameters accordingly (6). Proper lumping is a specific case of lumping methods where each state of the original model contributes to only one of the states of the reduced model. Therefore, the parameters of the reduced model can be directly related to those of the original model. This is of particular value in a Bayesian framework since prior knowledge on the original model can still be used to construct a prior distribution for the parameters of the reduced model. WBPBPK models are usually reduced in a way that the compartments of clinical interest remain unlumped. The reduced model should be able to describe the main features of the original model (e.g. response in blood and target site) while perhaps neglecting some less important aspects (e.g. response in the splanchnic organs). An automatic algorithm for proper lumping that can handle constraints (to avoid all possible combination of states) has been proposed by Dokoumetzidis and Aarons (7). In general, a model is reduced specifically for a particular set of parameter values. However, parameters always carry uncertainty and should hence be considered as random variables following statistical distributions rather than single values. Dokoumetzidis and Aarons proposed a robust Bayesian method to account for parameter uncertainty during the lumping process (8).

We have recently reported the MCMC analysis of mavoglurant (MVG) Phase-I clinical data using a generic WBPBPK model (9). We have found that in spite of assigning informative prior distributions to all population median parameters, the posterior distribution of the parameters could not converge to equilibrium, which suggested that the model was probably numerically unidentifiable. We therefore proposed to estimate only the population median parameters that have a significant influence on the plasma response (6 out of 15 parameters) and fix the others based on a sensitivity analysis of the model. However, although this approach permitted convergence to be achieved, it can be deemed too subjective as we assumed a priori that the available plasma data would not inform the majority of the parameters. In addition, fixing parameters while optimising some others can distort the covariance structure of the parameters and produce biased estimates (10).

The aim of the current study was thus to propose an alternative approach to stabilise the Bayesian analysis of clinical data with a PBPK model, which is the use of a reduced model obtained by proper lumping of a WBPBPK model in order to retain the mechanistic nature of the system. The WBPBPK model for MVG was used to illustrate this approach. The second goal of this study was to investigate whether the reduced yet mechanistic model could still be used to extrapolate MVG pharmacokinetics from adults to children as it was successfully done in our previous study of the WBPBPK model (9).

METHODS

The work flow for reducing the WBPBPK model for MVG, fitting the reduced model to clinical data and using the optimised reduced model for extrapolation from adults to children, includes the following steps:

  1. 1.

    Development of the original model

  2. 2.

    Model order reduction using proper lumping

  3. 3.

    Sensitivity analysis of the reduced model

  4. 4.

    Bayesian data analysis

  5. 5.

    Extrapolation

Thereby, prior information on the model parameters is propagated all along the modelling work flow and is combined with the information in the data to then perform possibly better predictions.

WBPBPK Model for MVG

The WBPBPK model used to describe the time course of MVG plasma concentrations after intravenous (IV) administration in healthy adult subjects has been described in detail by Wendling et al. (9) and is depicted in Fig. 1. Briefly, the model comprises 14 tissue compartments (lungs (LU), heart (HT), brain (BR), muscle (MU), adipose (AD), skin (SK), spleen (SP), pancreas (PA), liver (LI), stomach (ST), gut (GU), bones (BO), kidneys (KI) and rest-of-body (RB)) and 2 blood compartments (arterial and mixed venous), that is 16 states in total. Each tissue compartment is assumed well-stirred, with the extent of distribution being characterised by the equilibrium tissue-to-blood partition coefficient (K b, T). The rate equation for the tissue compartments can be expressed as follows:

Fig. 1
figure 1

Schematic representation of the WBPBPK model for MVG. See text for definition of symbols. The numbers refer to the state numbers of the system

Full size image
$$ \frac{\mathrm{d}{A}_{\mathrm{T}}}{\mathrm{d}t}=\frac{Q_{\mathrm{T}}}{V_{\mathrm{VEN}/\mathrm{ART}}}\cdot {A}_{\mathrm{VEN}/\mathrm{ART}}-\frac{Q_{\mathrm{T}}}{V_{\mathrm{T}}{K}_{\mathrm{b},\ \mathrm{T}}}\cdot {A}_{\mathrm{T}} $$
(1)

where A T denotes the amount (μg), V T the volume (l) and Q T the blood flow (l/h) for the different tissues, and A VEN/ART and V VEN/ART are the amount (μg) and volume (l), respectively, of either mixed venous blood (for the LU) or arterial blood (all other tissues). Drug elimination is assumed to occur entirely in the liver compartment via oxidative metabolism, as extra-hepatic metabolism and renal excretion are thought to be minor for MVG (11). The rate equation for the liver is thus defined as:

$$ \begin{array}{l}\frac{\mathrm{d}{A}_{\mathrm{LI}}}{\mathrm{d}t}=\frac{Q_{\mathrm{HA}}}{V_{\mathrm{ART}}}\cdot {A}_{\mathrm{ART}}+\frac{Q_{\mathrm{SP}}}{V_{\mathrm{SP}}{K}_{\mathrm{b},\ \mathrm{S}\mathrm{P}}}\cdot {A}_{\mathrm{SP}}+\frac{Q_{PA}}{V_{\mathrm{PA}}{K}_{\mathrm{b},\ \mathrm{P}\mathrm{A}}}\cdot {A}_{\mathrm{PA}}+\frac{Q_{\mathrm{ST}}}{V_{\mathrm{ST}}{K}_{\mathrm{b},\ \mathrm{S}\mathrm{T}}}\cdot {A}_{\mathrm{ST}}\\ {}+\frac{Q_{\mathrm{GU}}}{V_{\mathrm{GU}}{K}_{\mathrm{b},\ \mathrm{G}\mathrm{U}}}\cdot {A}_{\mathrm{GU}}-\frac{Q_{\mathrm{LI}}}{V_{\mathrm{LI}}{K}_{\mathrm{b},\ \mathrm{L}\mathrm{I}}}\cdot {A}_{\mathrm{LI}}-\frac{{\mathrm{CL}}_{\mathrm{int},\ \mathrm{L}\mathrm{I}}{\mathrm{fu}}_{\mathrm{b}}}{V_{\mathrm{LI}}{K}_{\mathrm{b},\ \mathrm{L}\mathrm{I}}}\cdot {A}_{\mathrm{LI}}\end{array} $$
(2)

where Q HA is the blood flow (l/h) from the hepatic artery, CLint, LI is the intrinsic clearance in the liver (l/h) and fub is the fraction unbound in blood.

Since the model is highly mechanistic, all parameters have a biological/physiological interpretation and can be defined a priori based on former experiments. These parameters are of two types: system-specific and drug-specific. The organ/tissue volumes and regional blood flows are the system-specific parameters and are typically gathered from the anatomy/physiology literature. The drug-specific parameters are the CLint, LI and the K b, T s and were extrapolated from an in vitro metabolism experiment and a distribution study in rats, respectively, as described in (9). Both system- and drug-specific parameter values can carry uncertainty due to errors in assumptions, hypotheses, observations, experiments and handling of the system studied (12). However, no uncertainty in the system-specific parameters was considered in our previous Bayesian analysis with the WBPBPK model in order to decrease the computational burden and potential numerical instabilities (9). On the other hand, a multivariate log-normal distribution was assumed for the prior distribution of the drug-specific parameters with uncertainty coefficients of variation (CV) of 26% for the CLint, LI value and of 30% for the K b, T values. The values of the organ/tissue blood flows and volumes as well as the marginal prior distributions of CLint, LI and K b, T s can be found in (9).

Reduction of the WBPBPK Model for MVG

The lumping methodology applied in this study has been described in detail in (7). Nevertheless, to understand how the ODEs for the reduced system were extracted, the basics of this technique are briefly described below.

The WBPBPK model for MVG can be expressed as

$$ \frac{\mathrm{d}y}{\mathrm{d}t}=f(y)=K\cdot y $$
(3)

where y is a vector of n states representing the amount of drug in the organ spaces and tissues, and K is the corresponding matrix of transfer rate constants. y can be reduced to a vector of \( \widehat{n}<n \) pseudo-states ŷ such that:

$$ \widehat{y}=M\cdot y $$
(4)

where M is the lumping matrix of dimension \( \widehat{n}\times n \) that determines the transformation from y to ŷ. The inverse transformation is given by y = M + ⋅ ŷ, where M + is the Moore-Penrose inverse of M. A new set of ODEs that described the reduced system can then be defined as follows:

$$ \frac{\mathrm{d}\widehat{y}}{\mathrm{d}t}=\widehat{f}\left(\widehat{y}\right)=M\cdot f\left({M}^{+}\widehat{y}\right) $$
(5)

or

$$ \frac{\mathrm{d}\widehat{y}}{\mathrm{d}t}=\widehat{K}\cdot \widehat{y} $$
(6)

where \( \widehat{K} \) is the matrix of transfer rate constants for the reduced system and can be calculated for a linear system as in Eq. 7.

$$ \widehat{K}=M\cdot K\cdot {M}^{+} $$
(7)

The lumping matrix M contains only 1 and 0 s, which describe the lumping scheme. The objective of model order reduction is to determine an appropriate matrix M such that the reduced system satisfies a specific property or criterion. We used the algorithm proposed by Dokoumetzidis and Aarons (7) to determine the optimal lumping matrix M such that an objective function accounting for the statistical distributions of the parameters (BOF) is minimised (Eq. 8).

$$ \mathrm{B}\mathrm{O}\mathrm{F}={\displaystyle \underset{S}{\int}\left(P\left(\theta \right)\cdot OF\right)\mathrm{d}\theta } $$
(8)

In Eq. 8, P(θ) is the probability density function of the prior distribution of the parameters θ, S is the parameter space and OF is the local objective function that takes the sum of squares of the differences of the responses between the reduced and original systems for all states (8).

Some advantages of this algorithm are: it is automated and was developed to combat combinatorial explosion; it can handle constraints that have a physiological meaning or force some organs/tissues of clinical interest to remain unlumped; it accounts for the prior distribution of the parameters as in a Bayesian modelling framework. Using this algorithm, the lumping matrix M is calculated for a desired number of pseudo-states \( \widehat{n} \) of the reduced system. Thus, determining the optimal lumping matrix M also involves determining manually the minimal value of \( \widehat{n} \) so that the reduced model still satisfies the property/criterion of interest. This was done in two steps. Firstly the predictions in the lumped states for the original and reduced systems were visually compared. To do so, the states of the original system were lumped in the same way as the reduced system but after the simulation using Eq. 4. Note that at this stage, all simulations and calculations done to determine the optimal lumping scheme were based on amounts as the ODEs for the WBPBPK model were defined to describe the temporal change of drug amount in the organ/tissue compartments. Secondly, the plasma concentration-time profiles predicted with the reduced and original models were visually compared. All calculations/simulations were performed using MATLAB R2014a (The MathWorks, Inc., Natick, Massachusetts, USA) and were speeded up using the Parallel Toolbox on a 4 core Intel® Core™ i5-3570 processor (3.40 GHz, 8 GB of RAM). Once the optimal lumping matrix M was obtained, the ODEs for the reduced model were determined using Eqs. 6 and 7 and the Symbolic Math Toolbox in MATLAB.

The way a PBPK model is reduced is specific to the aim of the subsequent modelling exercise. For MVG, the objective was to have a simplified model that would be able to predict the response in plasma and the brain (target site of the drug) similarly to the WBPBPK model. If no constraints are imposed, any state could be grouped with any other as the lumping method groups states purely based on kinetics (8). Hence, the brain was forced to remain unlumped, and the blood compartments were allowed to be lumped only together as well as with the lungs since these compartments are in series (Fig. 1). Similarly, the liver compartment (elimination site) was forced to remain unlumped and the splanchnic organs (SP, PA, ST and GU) were allowed to be grouped only together to maintain a clear anatomical interpretation. This also allows the incorporation of a mechanistic absorption model to predict pharmacokinetics after oral administration of MVG as it will be done later in this study.

Bayesian Analysis of MVG Clinical Data Using the Reduced Model

Data

The simplified PBPK model for MVG was optimised based on data from a Phase-I clinical study in healthy adult subjects. The design of the study and characteristics of the subjects enrolled have been described in (13). Briefly, plasma concentration-time data were collected from 120 healthy volunteers after a single 10-min IV infusion of 25, 37.5 and 50 mg of MVG. The subjects were young (median age of 31 years) Caucasian (62%) males (87%) with median body weight (BW) of 83 kg and median body mass index of 27 kg/m2. The study was conducted according to the ethical principles of the Declaration of Helsinki, and the protocol was approved by the Independent Ethics Committee for the study centre. Participants were males and nonpregnant females, and all provided full written informed consent prior to inclusion in the study.

Structural Identifiability of the Reduced Model

It has been shown that a generic WBPBPK model, with well-stirred compartments and one elimination site, is globally structurally identifiable even when only blood/plasma concentrations are measured (14). This is because for each tissue compartment, the influx rate k in, T is known a priori as it is defined as k in, T = Q T/V VEN/ART (see Eq. 1) and Q T and V VEN/ART are assumed known. This is true whether each tissue sub-system is parameterised in terms of efflux rate k out, T (Eq. 9) or in terms of blood flow Q T, volume V T and partition coefficient K b, T as it is assumed that K b, T is the only unknown variable.

$$ {k}_{\mathrm{out},\ \mathrm{T}}=\frac{Q_{\mathrm{T}}}{V_{\mathrm{T}}{K}_{\mathrm{b},\ \mathrm{T}}} $$
(9)

In general, any PBPK model based on the well-perfused mammillary model should be globally identifiable if the influx rates are known and only the efflux rates are estimated. Using proper lumping, the transfer rate constants for the reduced model, obtained by extracting the matrix \( \widehat{K} \) (Eq. 7), can be directly related to the blood flows, volumes and partition coefficients for the organ/tissue compartments of the original WBPBPK model. Thus, for simulation, the model can still be parameterised in terms of the original model parameters. However, the efflux rates for the lumped compartments are functions of the blood flows, volumes and partition coefficients of all the tissue compartments that were merged together (shown in the ‘RESULTS’) and are therefore functions of more than one unknown variables (partition coefficients). Hence, for the reduced model to be globally structurally identifiable during its optimisation, the lumped compartments have to be parameterised in terms of efflux rates. The unlumped compartments however, can still be parameterised in terms of the original model parameters.

Statistical Model

MVG plasma data were analysed using a three-stage hierarchical model to account for both population variability and uncertainty in the reduced PBPK model parameters. The same statistical model was used and described in detail in our previous analysis with the WBPBPK model (9). Briefly, we assumed at the first stage a log-normally distributed residual error (mean zero, variance σ 2) to account for model misspecification, unaccounted intra-individual variability in the parameters and measurement error. At the second stage, the model parameters were assumed to follow a multivariate log-normal distribution in the population. At the third stage, assumptions on the prior distribution of both the population and individual parameters were made, i.e. a multivariate log-normal distribution for the population median parameters and an inverse-Wishart distribution for the variability variance-covariance matrix (diagonal elements denoted ω 2).

Definition of Priors

The prior distributions of the WBPBPK model parameters (9) were used to construct priors for the reduced model parameters. When fitting a complex mechanistic model to data, providing prior information on the parameters is important to reduce numerical instabilities and impose biological/physiological constraints on the parameter estimates. For the compartments that remained unlumped, the marginal prior distributions of the K b, T s (and of CLint, LI for the liver) could be readily used. However, additional calculations were required to construct the marginal prior distributions of the efflux rates for the lumped compartments. While the means of the distributions could be directly calculated using the relationship between the efflux rates and the original model parameters, calculation of the variances required applying the delta method, which is a technique for approximating the moments of functions of random variables (15), as it will be shown in the RESULTS.

In our previous study, although we had no prior preclinical information on the variability in the WBPBPK drug-specific parameters, we could accurately estimate from the data the variability in CLint, LI using an uninformative prior (9). However, we did not attempt to quantify the variability in the other 14 parameters (partition coefficients) in order to decrease the computational burden. In the present study, we performed a first analysis run under similar conditions, i.e. with a random effect only on CLint, LI. Depending on how much the model and thus the number of parameters could be reduced, we then ran the analysis again with random effects on all parameters. This was to check the impact of the statistical model on the stability of the Bayesian analysis. In each case, an uninformative prior distribution was assigned to the variance-covariance matrix by setting the degrees of freedom for the inverse-Wishart distribution equal to the size of the matrix (16).

Sensitivity Analysis of the Model

A sensitivity analysis of the reduced PBPK model was performed to identify the parameters that have a significant influence on the plasma response and are thus expected to be updated by the clinical data. To account for uncertainty in the parameter values, 1000 sets of parameters were randomly drawn from the multivariate log-normal prior distribution. The sensitivity of the model to each parameter was then assessed by calculating and plotting over time a relative sensitivity coefficient as explained in (9).

MCMC Simulation

The posterior distribution of the parameters was summarised using random draws by Gibbs sampling as implemented in NONMEM 7.3.0 (ICON Development Solutions, Hanover, Maryland, USA). Three independent Markov chains were run in parallel for one million iterations, starting at different diffuse parameter values. Computations were performed on a medium-sized cluster running Red Hat Enterprise Linux 6.5 on nodes with Intel Xeon E5-2670v2 CPUs with some older nodes running with dual Xeon X5670 CPUs. Each chain was parallelized on 12 different nodes. The nodes are equipped with between 24 and 96 GB of RAM and are interconnected via dedicated Bonded 1 Gbit network cards. Convergence of the posterior distribution to approximate equilibrium was monitored using the potential scale reduction statistic proposed by Gelman and Rubin (\( \widehat{R} \)) (17), as well as by visual inspection of the Markov chains. \( \widehat{R} \) values were calculated using the software package CODA (18). The reduced ODE system was solved using the LSODA solver as implemented in the ADVAN13 NONMEM subroutine.

Scaling the Reduced PBPK Model from Adults to Children

In our previous study, the WBPBPK model could extrapolate reasonably well MVG oral pharmacokinetics from adults to children (9). This was done in three steps. Firstly, the WBPBPK model was fitted to IV adult data by estimating only the most influential parameters, as explained in the introduction. Secondly, to be able to predict pharmacokinetics after oral administration of the drug in adults, a mechanistic absorption model incorporating system-specific (e.g. gastric and small intestine transit times), drug-specific (e.g. solubility) and formulation-specific (e.g. dissolution constant) parameters was implemented. The reader is referred to (9) for details on the structure of the model as well as on the distribution of the absorption parameters. Thirdly, the full model was scaled from adult to children and used for Monte Carlo simulation. The predictive performance of the model for children was visually evaluated using data from a Phase-I clinical study of a single 15-mg dose of MVG administered to 21 subjects aged from 3 to 11 years (9).

In the current study, the same three steps were followed to evaluate the ability of the reduced model to predict MVG oral pharmacokinetics in children. Although at the third step, the reduced model predictions were checked not only against the paediatric data but also against the WBPBPK model predictions. For both the reduced and whole-body model, the Monte Carlo simulation were performed exactly as in (9) using MATLAB. Since the liver compartment was forced to remain unlumped during the reduction process, the absorption model could be integrated in the reduced PBPK model as for the WBPBPK model.

A generic WBPBPK model can be easily scaled from adults to children by incorporating in the model the age-related changes in the physiological parameters (i.e. organ/tissue volumes and blood flows) as well as in some drug-specific parameters such as the intrinsic clearance and the fraction unbound in plasma, as it was done for MVG in our previous study (9). The same can be done for a formally reduced model as the parameters retain mechanistic meanings. Yet, the efflux rates for the lumped compartments are a mix of physiological and drug-specific parameters and are therefore also expected to vary with age. It is thus not sufficient to scale the model from adult to children by only accounting for the age-related changes in the physiological parameters, but the efflux rates should also be scaled, especially if most compartments are pseudo-states resulting from lumping. One approach to scale efflux rates across age groups is to use the following allometric scaling relationship (19):

$$ {k}_{\mathrm{i}}={k}_{\mathrm{std}}{\left(\frac{{\mathrm{BW}}_i}{70}\right)}^{-0.25} $$
(10)

where k i is the efflux rate for the ith individual with body weight BWi, and k std is the efflux rate for a standard individual in the population. Hence, the approach that we propose to scale a formally reduced PBPK model across age groups involves both what we call ‘physiological scaling’ and allometric scaling. The method used for physiological scaling is described in detail in (9).

RESULTS

Reduction of the WBPBPK Model for MVG

To account for parameter uncertainty during the order reduction of the WBPBPK model for MVG, the Bayesian objective function BOF (Eq. 8) implemented in the lumping algorithm that we used was minimised based on 1000 parameter sets drawn by Latin hypercube sampling from the prior distribution of the parameters. Calculation of the objective function value for each of the 1000 samples is computationally intensive as it requires solving the original ODE system 1000 times. Nevertheless, parallelising the optimisation on 4 cores helped to reduce the run time to less than an hour.

The lumping schemes determined for reduction of the 16 original states to 7 and 6 states are shown in Table I, columns 2 and 3, respectively. For both schemes, the blood compartments were lumped together with the lungs to form a central compartment from which plasma concentrations are derived. Figure 2 shows plots of the reduced model predictions versus the original model predictions for each lumped state. Using a seven-state system, the predictions in most lumped states are reasonably consistent with the original system predictions except in the 2nd pseudo-state, which is a lump of the HT, SK, KI, BO and RB tissues, and in the 6th pseudo-state, which is the splanchnic organs (SP, PA, ST and GU) lumped together (Fig. 2a). However, predictions of drug amount in these two lumped compartments are not of any particular clinical interest for MVG. In addition, Fig. 3 shows that the deviation from the original system in these lumped compartments barely affects the time course of drug concentration in plasma. Reducing the system one step further resulted in lumping the MU and AD tissues together. The predictions in this additional lumped state are on average deviating from the original model predictions (Fig. 2b) which seems to significantly affect MVG plasma concentration-time profile as seen on Fig. 3, especially in the terminal phase of the profile, i.e. from 36 h post-dose. As a consequence, the system reduced to seven states was selected as a simplified PBPK model for MVG. The structure of the model is depicted in Fig. 4. The matrix of transfer rates for this model was calculated using Eq. 7 and can be expressed as follows:

Table I Lumping Scheme of Sizes 7 and 6 States
Full size table
Fig. 2
figure 2

Drug amount in the lumped states predicted with a reduced system of size 7 (a) and size 6 (b) plotted against the original system predictions. The simulation was performed every 0.1 h for 48 h using a 50-mg IV dose (used in the clinical study) for 1000 parameter sets sampled from the prior distribution. The subtitles refer to the state numbers of the original system (see Fig. 1; Table I)

Full size image
Fig. 3
figure 3

Plots of the confidence intervals (90%) around MVG median plasma concentration-time profiles after an IV dose of 50 mg. Confidence intervals produced with the original WBPBPK model and with the seven-state reduced model are represented by the solid black lines and solid red lines (superimposed), respectively, whereas the intervals are shown by the dashed blue lines for the six-state reduced model. The simulation was performed every 0.1 h for 48 h, for 1000 parameter sets sampled from the prior distribution. Plasma concentrations were derived from the amount of drug in the venous blood compartment for the original system and in the central compartment for the reduced systems

Full size image
Fig. 4
figure 4

Schematic representation of the seven-state reduced PBPK model selected for the Bayesian analysis of MVG clinical data. The numbers refer to the state numbers of the reduced system. CEN central compartment, LUM dummy lumped compartment, SPL splanchnic compartment. See text for definition of the other symbols

Full size image
$$ \widehat{K}=\left(\begin{array}{ccccccc}\hfill -{Q}_{\mathrm{LU}}/{V}_{\mathrm{CE}\mathrm{N}}\hfill & \hfill {k}_{\mathrm{out},\ \mathrm{L}\mathrm{U}\mathrm{M}}\hfill & \hfill {Q}_{\mathrm{BR}}/\left({K}_{\mathrm{b},\ \mathrm{B}\mathrm{R}}\cdot {V}_{\mathrm{BR}}\right)\hfill & \hfill {Q}_{\mathrm{MU}}/\left({K}_{\mathrm{b},\ \mathrm{M}\mathrm{U}}\cdot {V}_{\mathrm{MU}}\right)\hfill & \hfill {Q}_{\mathrm{AD}}/\left({K}_{\mathrm{b},\ \mathrm{AD}}\cdot {V}_{\mathrm{AD}}\right)\hfill & \hfill 0\hfill & \hfill {Q}_{\mathrm{LI}}/\left({K}_{\mathrm{b},\ \mathrm{L}\mathrm{I}}\cdot {V}_{\mathrm{LI}}\right)\hfill \\ {}\hfill {Q}_{\mathrm{LU}\mathrm{M}}/{V}_{\mathrm{CE}N}\hfill & \hfill -{k}_{\mathrm{out},\ \mathrm{L}\mathrm{U}\mathrm{M}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill {Q}_{\mathrm{BR}}/{V}_{\mathrm{CE}\mathrm{N}}\hfill & \hfill 0\hfill & \hfill -{Q}_{\mathrm{BR}}/\left({K}_{\mathrm{b},\ \mathrm{B}\mathrm{R}}\cdot {V}_{\mathrm{BR}}\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill {Q}_{\mathrm{MU}}/{V}_{\mathrm{CE}\mathrm{N}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{Q}_{\mathrm{MU}}/\left({K}_{\mathrm{b},\ \mathrm{M}\mathrm{U}}\cdot {V}_{\mathrm{MU}}\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill {Q}_{\mathrm{AD}}/{V}_{\mathrm{CE}\mathrm{N}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{Q}_{\mathrm{AD}}/\left({K}_{\mathrm{b},\ \mathrm{AD}}\cdot {V}_{\mathrm{AD}}\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill {Q}_{\mathrm{SPL}}/{V}_{\mathrm{CE}\mathrm{N}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{k}_{\mathrm{out},\ \mathrm{S}\mathrm{P}\mathrm{L}}\hfill & \hfill 0\hfill \\ {}\hfill {Q}_{\mathrm{HA}}/{V}_{\mathrm{CE}\mathrm{N}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{\mathrm{out},\ \mathrm{S}\mathrm{P}\mathrm{L}}\hfill & \hfill -{Q}_{\mathrm{LI}}/\left({K}_{\mathrm{b},\ \mathrm{L}\mathrm{I}}\cdot {V}_{\mathrm{LI}}\right)-\left({\mathrm{CL}}_{\mathrm{int},\ \mathrm{L}\mathrm{I}}\cdot {\mathrm{fu}}_{\mathrm{b}}\right)/\left({K}_{\mathrm{b},\ \mathrm{L}\mathrm{I}}\cdot {V}_{\mathrm{LI}}\right)\hfill \end{array}\right) $$
(11)

with

$$ \begin{array}{l}{V}_{\mathrm{CEN}}=3\cdot {V}_{\mathrm{ART}}\\ {}{Q}_{\mathrm{LUM}}={Q}_{\mathrm{HT}}+{Q}_{\mathrm{SK}}+{Q}_{\mathrm{KI}}+{Q}_{\mathrm{BO}}+{Q}_{\mathrm{RB}}\\ {}{Q}_{\mathrm{SP}\mathrm{L}}={Q}_{\mathrm{SP}}+{Q}_{\mathrm{PA}}+{Q}_{\mathrm{ST}}+{Q}_{\mathrm{GU}}\\ {}{Q}_{\mathrm{HA}}={Q}_{\mathrm{LI}}-{Q}_{\mathrm{SP}\mathrm{L}}\\ {}{k}_{\mathrm{out},\ \mathrm{L}\mathrm{U}\mathrm{M}}={Q}_{\mathrm{HT}}/\left(5\cdot {K}_{\mathrm{b},\ \mathrm{H}\mathrm{T}}\cdot {V}_{\mathrm{HT}}\right)+{Q}_{\mathrm{SK}}/\left(5\cdot {K}_{\mathrm{b},\ \mathrm{S}\mathrm{K}}\cdot {V}_{\mathrm{SK}}\right)+{Q}_{\mathrm{KI}}/\left(5\cdot {K}_{\mathrm{b},\ \mathrm{K}\mathrm{I}}\cdot {V}_{\mathrm{KI}}\right)+{Q}_{\mathrm{BO}}/\left(5\cdot {K}_{\mathrm{b},\ \mathrm{B}\mathrm{O}}\cdot {V}_{\mathrm{BO}}\right)+{Q}_{\mathrm{RB}}/\left(5\cdot {K}_{\mathrm{b},\ \mathrm{R}\mathrm{B}}\cdot {V}_{\mathrm{RB}}\right)\\ {}{k}_{\mathrm{out},\ \mathrm{S}\mathrm{P}\mathrm{L}}={Q}_{\mathrm{SP}}/\left(4\cdot {K}_{\mathrm{b},\ \mathrm{S}\mathrm{P}}\cdot {V}_{\mathrm{SP}}\right)+{Q}_{\mathrm{PA}}/\left(4\cdot {K}_{\mathrm{b},\ \mathrm{P}\mathrm{A}}\cdot {V}_{\mathrm{PA}}\right)+{Q}_{\mathrm{ST}}/\left(4\cdot {K}_{\mathrm{b},\ \mathrm{S}\mathrm{T}}\cdot {V}_{\mathrm{ST}}\right)+{Q}_{\mathrm{GU}}/\left(4\cdot {K}_{\mathrm{b},\ \mathrm{G}\mathrm{U}}\cdot {V}_{\mathrm{GU}}\right)\end{array} $$
(12)

where fub is the fraction unbound in blood calculated as the ratio of the fraction unbound in plasma to the blood-to-plasma ratio, V CEN is the volume (l) for the CEN compartment (from which plasma concentrations are derived), k out, LUM and Q LUM are the efflux rate (h−1) and blood flow (l/h), respectively, for the LUM compartment, Q HA is the blood flow (l/h) from the hepatic artery, and k out, SPL and Q SPL are the efflux rate (h−1) and blood flow (l/h), respectively, for the SPL compartment (see Fig. 4 and Table I for definition of the lumped compartments). It should be stressed that the relationships in Eqs. 11 and 12 are direct results from the matrix operation in Eq. 7.

Bayesian Analysis of MVG Clinical Data Using the Reduced Model

Structural Identifiability of the Reduced Model

In the WBPBPK model for MVG, the unknown variables are the CLint, LI and K b, T s. In the reduced model, the efflux rates for the LUM and SPL compartments are thus functions of more than one unknown variable (Eq. 12) meaning that for the model to be structurally identifiable, it cannot be fully parameterised in terms of the WBPBPK model parameters but has to be expressed as in Eq. 11. Hence, by reducing the system from 16 to 7 states, the number of parameters is reduced from 15 to 7 and the parameters are CLint, LI, k out, LUM, K b, BR, K b, MU, K b, AD, k out, SPL and K b, LI. BW was included as a covariate for the efflux rates k out, LUM and k out, SPL as in Eq. 10, in order to explain part of the population variability in these parameters and to further scale them from adults to children using allometry as explained in the ‘METHODS’.

Priors

The prior distribution of the reduced PBPK model parameters is presented in the first column of Table II. For k out, LUM and k out, SPL, the means of the marginal distributions were simply calculated by using the relationships in Eq. 12, and the blood flow and volume values as well as the prior means of the partition coefficients reported in (9). To calculate the variances, the following relationships were determined by applying the delta method:

Table II Prior and Posterior Distributions of the Drug-Specific Parameter Values for MVG-Reduced PBPK Model (STATMOD 2)
Full size table
$$ \begin{array}{l}\operatorname{var}\left({k}_{\mathrm{out},\ \mathrm{L}\mathrm{U}\mathrm{M}}\right)={\mathrm{CV}}^2\left[{\left(\frac{Q_{\mathrm{HT}}}{5\cdot {K}_{\mathrm{b},\ \mathrm{H}\mathrm{T}}\cdot {V}_{\mathrm{HT}}}\right)}^2+{\left(\frac{Q_{\mathrm{SK}}}{5\cdot {K}_{\mathrm{b},\ \mathrm{S}\mathrm{K}}\cdot {V}_{\mathrm{SK}}}\right)}^2+{\left(\frac{Q_{\mathrm{KI}}}{5\cdot {K}_{\mathrm{b},\ \mathrm{K}\mathrm{I}}\cdot {V}_{\mathrm{KI}}}\right)}^2+{\left(\frac{Q_{\mathrm{BO}}}{5\cdot {K}_{\mathrm{b},\ \mathrm{B}\mathrm{O}}\cdot {V}_{\mathrm{BO}}}\right)}^2+{\left(\frac{Q_{\mathrm{RB}}}{5\cdot {K}_{\mathrm{b},\ \mathrm{R}\mathrm{B}}\cdot {V}_{\mathrm{RB}}}\right)}^2\right]\\ {}\operatorname{var}\left({k}_{\mathrm{out},\ \mathrm{S}\mathrm{P}\mathrm{L}}\right)={\mathrm{CV}}^2\left[{\left(\frac{Q_{\mathrm{SP}}}{4\cdot {K}_{\mathrm{b},\ \mathrm{S}\mathrm{P}}\cdot {V}_{\mathrm{SP}}}\right)}^2+{\left(\frac{Q_{\mathrm{PA}}}{4\cdot {K}_{\mathrm{b},\ \mathrm{P}\mathrm{A}}\cdot {V}_{\mathrm{PA}}}\right)}^2+{\left(\frac{Q_{\mathrm{ST}}}{4\cdot {K}_{\mathrm{b},\ \mathrm{S}\mathrm{T}}\cdot {V}_{\mathrm{ST}}}\right)}^2+{\left(\frac{Q_{\mathrm{GU}}}{4\cdot {K}_{\mathrm{b},\ \mathrm{G}\mathrm{U}}\cdot {V}_{\mathrm{GU}}}\right)}^2\right]\end{array} $$
(13)

It should be noted that the relationships in Eq. 13 could be simplified because we assumed in our previous Bayesian analysis the same coefficient of variation CV for all K b, T s.

Sensitivity Analysis of the Model

The results of the sensitivity analysis for the reduced model shown in Fig. 5 suggest that the model is sensitive to all parameters except to k out, SPL and K b, LI that seem to have a negligible influence on the response in the central compartment and thus on the plasma response. Uncertainties in these parameters are therefore not expected to be reduced by the clinical data. On the other hand, as for the WBPBPK model, CLint, LI, K b, MU and K b, AD are the parameters that have the main influence on the plasma response and that might have their marginal distributions updated by the data.

Fig. 5
figure 5

Absolute value of the relative sensitivity coefficient for the central compartment of the reduced model |RSCEN| (solid grey lines) plotted against time for each parameter of each of the 1000 parameter vectors drawn from the prior distribution. The horizontal solid black line represents the threshold value of 0.1 used to assess the influence of a parameter on the response: a coefficient higher than this value suggests a significant influence on the response. See text for definition of symbols

Full size image

MCMC Simulation

The Bayesian analysis of MVG clinical data with the reduced PBPK model was run with two different statistical models: in the first, referred to as ‘STATMOD 1’, the effect of inter-individual variability was estimated only on CLint, LI and was fixed to a small value (corresponding to a variability CV of 1% for a log-normally distributed variable) for the other parameters; in the second, referred to as ‘STATMOD 2’, we estimated random effects on all seven parameters. For both statistical models, it took on average less than 3 days to run one million iterations on the cluster described in the ‘METHODS’ and by parallelising each Markov chain on 12 different nodes. The total run time was only few hours longer for STATMOD 2 than for STATMOD 1. The posterior distribution of the parameters is summarised in the supplemental online (Table I) for STATMOD 1 and in Table II for STATMOD 2 and was generated by pooling the one million samples from the three Markov chains. Using both models, the posterior distribution of the population and individual parameters converged to approximate equilibrium. However, as indicated by the \( \widehat{R} \) statistic, convergence was achieved after 360,000 iterations for STATMOD 1 versus 8000 iterations for STATMOD 2. In addition, the trace plots of the Markov chains suggest a much slower mixing of the chains for STATMOD 1 (Fig. 1 in the Supplemental online material) than for STATMOD 2 (Fig. 6). Slow mixing of Markov chains suggests that the MCMC samples are highly autocorrelated within each chain, which increases the uncertainty of the estimation of posterior quantities of interest (e.g. means, variances or quantiles). The degree of autocorrelation can be measured by the effective sample size of the sequence of interest. Using the CODA software package, we estimated for each model parameter the effective sample size of the posterior sequences obtained with STATMOD 1 and STATMOD 2 (Table II in the Supplemental online material). The smaller the effective sample size estimate, the higher autocorrelation between the MCMC samples and the higher uncertainty in the posterior quantities. The estimates suggest that the MCMC samples for K b, BR, k out, SPL and K b, LI are highly autocorrelated for STATMOD 1 whereas there are not so much for STATMOD 2 which was thus considered as numerically more stable. Finally, the estimates of the residual variance suggest that STAMOD 2 provides a better description of individual profiles than STATMOD 1. We therefore focus on STATMOD 2 in the following paragraph.

Fig. 6
figure 6

Trace plots of the three Markov chains run for the Bayesian analysis of MVG clinical data (STATMOD 2). Parameter values (log-transformed for the fixed effects) are plotted against the number of iterations (×105) of the MCMC simulation. See text for definition of symbols

Full size image

Monitoring the \( \widehat{R} \) statistic during the MCMC simulation (Table II) suggested that most marginal distributions seemed to have reached equilibrium after 2000 iterations except for K b, BR, k out, SPL, K b, LI and \( {\omega^2}_{K_{\mathrm{b},\ \mathrm{B}\mathrm{R}}} \) for which additional sampling was required. The marginal posterior and prior distributions of the population median parameters are compared in Table II (see also Fig. 2 of the Supplemental online material). All distributions were updated by the plasma data. The smallest deviation from the prior mean was for CLint, LI (18%) whereas the most notable change was for K b, LI (68%). However, uncertainty in the parameter value was only slightly decreased for K b, LI (from 30% to 25% CV) and was even increased for k out, SPL (from 17% to 36% CV). The estimates of the inter-individual variances suggest that K b, MU is the most variable parameter in the population (CV of 72%) and that the data had no valuable information for the variability in K b, BR, k out, SPL and K b, LI (uncertainty CV of 130%, 197% and 156% for \( {\omega^2}_{K_{\mathrm{b},\ \mathrm{B}\mathrm{R}}} \), \( {\omega^2}_{k_{\mathrm{out},\ \mathrm{S}\mathrm{P}\mathrm{L}}} \) and \( {\omega^2}_{K_{\mathrm{b},\ \mathrm{L}\mathrm{I}}} \), respectively). As seen in Fig. 7, the population model seems to provide a good description of both the median trend and the variability in the analysed plasma concentration-time data.

Fig. 7
figure 7

A visual predictive check of the ability of the population reduced PBPK model (STATMOD 2) to describe MVG plasma data after IV administration in adult subjects. Open circles are observed concentrations plotted across time, the solid red line is the median of the simulated concentrations and the grey area represents a 90% prediction interval. Both observed and predicted concentrations were dose - normalised. The insert expands the first 2 h of the concentration-time data plotted on linear scales

Full size image

Scaling the Reduced PBPK Model from Adults to Children

Prior to scaling the reduced PBPK model from adults to children, we checked that the reduced model together with the mechanistic absorption model could predict adequately MVG oral pharmacokinetics in adults (results not presented) as it was done for the WBPBPK model in (9). We then applied physiological scaling to the organ/tissue volumes and blood flows, the CLint, LI, the fraction unbound in plasma and the volume of fluid in the small intestine (absorption parameter), and allometric scaling to the efflux rates k out, LUM and k out, SPL (Eq. 9). For the purpose of comparing the performance of the reduced model with that of the whole-body model in extrapolating MVG pharmacokinetics from adults to children, STATMOD 1 was used for the Monte Carlo simulation with the reduced model as it was also used in our previous simulation with the WBPBPK model (9). The results of the simulation are shown in Fig. 8. The scaled reduced model appears to be consistent with the scaled WBPBPK model. Both models slightly over-predict the median trend in the data in the first five hours post-dose. As explained in our previous study (9), this is likely due to unaccounted age-related changes in the absorption parameters which are the same between the reduced and the original model. Note that the predictions made with the reduced and the whole-body model are not expected to be identical simply because the reduced model is an approximation of the whole-body model. Hence, it is not surprising that analysing the same data yields slightly different posterior estimates of the population median and variance parameters and thus provides different prediction intervals. It should also be noted that the simulation was performed by simply using the posterior means of the population and individual parameters and of the residual error. Although simulating with uncertainty, by sampling the parameters from the posterior distribution, can be of value in drug development, e.g. to select a safe dosing regimen prior to a clinical study, it was deemed not critical in evaluating the predictive performance of the reduced model.

Fig. 8
figure 8

A visual evaluation of the ability of the scaled PBPK model to predict MVG pharmacokinetics after oral administration of a 15-mg dose in children aged from 3 to 11 years. Open circles are observed concentrations plotted across time, the dashed red lines are the 5th, 50th and 95th percentiles of the concentrations simulated with the WBPBPK model from Wendling et al. (9), the solid blue lines are the 5th, 50th and 95th percentiles of the concentrations simulated with the reduced PBPK model (STATMOD 1), and the horizontal dotted line represents the lower limit of quantification of the assay

Full size image

DISCUSSION

Bayesian population physiological modelling represents a useful translational approach to bridge preclinical and clinical pharmacokinetic data during drug development. Among PBPK models, whole-body models are especially suited for between- or within-species extrapolation as all the model parameters (system- and drug-specific) have clear mechanistic meanings. However, fitting a WBPBPK model to sparse and noisy clinical data is challenging due to the high number of parameters and due to the difficulty of defining priors that are informative enough given the complexity of the model. For instance, the WBPBPK model previously developed for MVG appeared to be numerically unidentifiable when analysing dense clinical data, possibly because of severe correlations between parameters and of too diffuse priors given the dimension of the system (9). One approach to reduce the number of parameters and make the model numerically identifiable is to fix the parameters that have a negligible influence on the plasma response to values a priori known from preclinical experiments (9). This approach has the problem of possibly producing biased estimates and underestimating the uncertainty in the system. Using MVG as an example, we therefore evaluated in the present study an alternative approach that is to formally reduce a WPBPK model and its number of parameters, and use the reduced model for the Bayesian data analysis. Finally, we proposed an approach to scale reduced PBPK models across age groups so that the model can still be used for extrapolation.

The WBPBPK model for MVG could be reduced from 16 to 7 states while maintaining the main features of the model, i.e. kinetics in the plasma (Fig. 3) and the brain (Fig. 2a). The fact that lumping the adipose tissue with the muscle tissue (six-state model) significantly affects the plasma kinetics (Fig. 3) is likely due to the extensive distribution of MVG into the adipose tissue which equilibrates much slower than the muscle tissue (see Fig. 3 in the Supplemental online material). Using the proper lumping algorithm proposed by Dokoumetzidis and Aarons (8), we could easily give emphasis to specific tissues (e.g. brain and liver) and organs spaces (e.g. mixed blood) and avoid combinations that have no anatomical meaning (e.g. splanchnic organs lumped with the brain). Thereby, we retained the mammillary structure of the model (Fig. 4) which offered the advantage of facilitating the determination of the structural identifiability of the model as explained in the ‘METHODS’. If the mammillary structure of the model is perturbed, a formal structural identifiability analysis should be performed. Another advantage was that a mechanistic absorption model could be directly integrated and can be used to investigate the bioavailability of the drug (first-pass effect) or even the impact of drug-drug interactions (e.g. inhibition of the hepatic metabolic clearance). Nevertheless, incorporating a drug-drug interaction mechanism in the model often means incorporating a nonlinear mechanism. Depending on the experimental conditions and on the characteristics of the drug, this can significantly affect kinetics in the tissues and hence the validity of the lumping scheme. If the model is meant to be used to predict the impact of drug-drug interactions as was done in our previous analysis with the WBPBPK model (9), we strongly recommend that the inhibition/induction mechanism is firstly incorporated in the model before it is reduced. A difficulty may arise in extracting the ODEs as the system would be nonlinear and Eq. 7 would not be valid anymore, but for simulation purpose, the ODEs of the reduced system are not needed as the system can be derived using Eq. 5.

The lumping algorithm appeared to be suitable for a Bayesian PBPK modelling framework as it takes into account the prior distribution of the parameters. Thereby, the lumping scheme is robust to variation in parameter values due to uncertainty that arises from the preclinical experiments and extrapolation methods. In addition, the parameters of the reduced model retain a clear physiological interpretation (Eqs. 11 and 12). Hence, we could easily propagate the prior information on the system- and drug-specific parameters of the WBPBPK model to all the parameters of the reduced model including the more empirical efflux rates for the lumped compartments. This was deemed important to stabilise the MCMC analysis and produce biologically/physiologically plausible posterior values.

As opposed to the MCMC analysis with the WBPBPK model (when estimating all drug-specific parameters), convergence of the posterior parameter distribution could be achieved with the reduced PBPK model (with both STATMOD 1 and STATMOD 2) without having to fix population median parameters. However, due to highly autocorrelated MCMC samples, intensive Gibbs sampling (360,000 samples) was required to reach approximate equilibrium when the statistical model was limited to a random effect on CLint, LI (STATMOD 1) as in our previous WBPBPK analysis (9). In population modelling, it is desirable to assign random effects to all model parameters to learn from the data about variability in an objective manner. Since the number of parameters was significantly reduced from 15 to 7, we could also assess a ‘full’ statistical model (STATMOD 2) which appeared to considerably decrease autocorrelations between the MCMC samples (see Table II of the supplemental online) and hence accomplish convergence to the target distribution much faster (8000 samples), even though uninformative priors were used for all population variances. Whether the impact of the statistical model on the stability of the analysis is specific to the MCMC algorithm implemented in NONMEM, still needs to be investigated by testing our model using other algorithms/software frequently used for Bayesian pharmacokinetic/pharmacodynamic analyses, e.g. Stan, WinBUGS or MCSim (20). Note that we have also tried a ‘full’ statistical model for the WBPBPK model to see if it would also stabilise the MCMC analysis and allow convergence, but the NONMEM runs for the three chains terminated after a few thousand iterations due to numerical issues (error message ‘OBJECTIVE FUNCTION IS INFINITE. PROBLEM ENDED’), which can probably be explained by the high dimension of the system, high correlations between parameters and the lack of prior information on the variability terms.

The Bayesian analysis of MVG clinical data provided better and more accurate estimates for the population median parameters CLint, LI, k out, LUM, K b, BR, K b, MU, K b, AD and K b, LI than our prior preclinical estimates. The uncertainty in the parameter values was especially reduced for CLint, LI, k out, LUM, K b, MU and K b, AD for which convergence of the marginal distributions was achieved after only 2000 iterations (Table II). In addition, population variability in these four parameters could be accurately estimated in contrast to the variability in K b, BR, k out, SPL and K b, LI. This suggests that the data had information mainly for CLint, LI, k out, LUM, K b, MU and K b, AD which is consistent with the results of the sensitivity analysis (see Fig. 5). No clear explanation could be given for the increase in uncertainty around the population median of k out, LUM although it might be due to numerical issues specific to the MCMC algorithm implemented in NONMEM. The variability estimate for k out, LUM (CV of 46%) suggests that the allometric relationship in Eq. 10 helps to explain only partially the overall variability in this parameter which is probably due to variation in the volumes and blood flows of the HT, SK, KI, BO and RB compartments. Also, the high population variance estimate for K b, MU is likely explained by unaccounted variability in the organ/tissue volumes and blood flows (21) perhaps because BW is not the most appropriate covariate for the volumes and blood flows. We used BW for simplicity, based on (22, 23). However, there is evidence that body height is a better explanatory variable for the population variability in the blood flows and volumes than BW (24). Particular attention to physiologic variability, as well as to parameter uncertainty, should be paid if the model were to be used for exploratory simulations to determine for instance a dosing regimen that best fits efficacy or safety requirements. In that regard, it is informative to quantify the variability in all drug-specific parameters (using a ‘full’ statistical model) in order to potentially refine the covariate model if the variance estimates are thought to be physiologically unrealistic.

Using proper lumping, the biological/physiological interpretation of the parameters is retained in the reduced model. In other words, the model can be expressed using the parameters of the original model (see Eqs. 11 and 12). Hence, scaling across age groups can easily be done by incorporation of the age-related anatomical and physiological changes in the subjects. A slight difficulty arises when the model is fitted to data prior to being used for extrapolation, as the rate equations for the lumped states have to be expressed using the more empirical efflux rate parameters (Eq. 11). Nevertheless, it is common to scale a parameter describing a transfer rate using the allometric relationship expressed in Eq. 10. Hence, a formally reduced PBPK model can be scaled across age groups as a WBPBPK model (changes in blood flows, volumes, intrinsic clearance, fraction unbound in plasma etc.) with additional allometric scaling of the efflux rates for the lumped states. For MVG, the extrapolation workflows with the WBPBPK model (9) and with the reduced PBPK model performed similarly (Fig. 8). Yet, the more the model is reduced, the less mechanistic it is and the more allometric scaling predominates. A major limit of allometric scaling is that it fails to account for age-related maturity changes in active processes such as enzymatic elimination or transporter-mediated efflux. Incorporating maturity differences in a PBPK model is critical to predict pharmacokinetics in young children (25). If appropriate constraints are imposed during the lumping process, emphasis can be given to active processes so that maturity and physiological differences can still be incorporated when scaling the model. The same method could be applied to between-species scaling as anatomical and physiological changes over organisms can be found in the literature. Although, a more robust method would be to account for the distribution of the parameter over different animal species during the lumping process (8). This also applies to within-species scaling like between age groups as it was done in the present study. However, defining the prior distribution of the parameters over different species or age groups can be difficult as the parameter values follow multimodal distributions rather than log-normal distributions.

The Bayesian workflow proposed in this study (see ‘METHODS’) could be applied to any complex mechanistic pharmacokinetic and pharmacodynamic model. The choice of reducing the system depends on its complexity and on the purpose of the subsequent simulation as the ability of a reduced model to extrapolate under different scenarios is limited compared to a fully mechanistic model (e.g. prediction of drug-drug interactions with a PBPK model). The difficulty with pharmacodynamic models is that the systems are often nonlinear due to feedback mechanisms. Hence, extraction of the ODEs to be able to fit the reduced model to data is not trivial. Also, determining the structural identifiability of the simplified model requires a formal identifiability analysis. Yet, such an approach has been recently successfully applied to the coagulation network in human, although a frequentist analysis was performed to further model fibrinogen concentration-time data following brown snake envenomation (26). The choice of reducing a model also depends on the amount of prior information available on the structural parameters as well as on the variability terms. Jonsson et al. have shown an example of the impact of the informativeness of the priors for population variances on the stability of the Bayesian analysis of a mechanistic pharmacokinetic/pharmacodynamic model (3). Recently, Garcia et al. have also shown that increasing the informativeness of the priors for variability terms can help a statistical PBPK model to be identifiable (4). However, in both analyses, the informativeness was increased by naively fixing the variability terms for which no priors were available to zero. If prior information is restricted to the population median parameters, it might be desirable to reduce the model to be able to quantify the variability in all parameters and thus learn more objectively about the overall variability in the system.

CONCLUSIONS

We conclude that the present study represents a first example of the practicality of formally reduced PBPK models for Bayesian inference in clinical pharmacokinetics. The Bayesian proper lumping method allowed propagation of MVG preclinical pharmacokinetic information from the original to the reduced system. Using MCMC techniques, the information could then be combined with clinical data to improve the population model without having to fix neither variability terms nor structural drug-specific parameters. Finally, the reduced yet mechanistic model could still be used for PBPK extrapolation across age groups. Nevertheless, more examples are needed to prove the advantages of such Bayesian mechanistic modelling approaches in drug development.