Population PBPK modeling using parametric and nonparametric methods of the Simcyp Simulator, and Bayesian samplers.

Physiologically-based pharmacokinetic (PBPK) models usually include a large number of parameters whose values are obtained using in vitro to in vivo extrapolation. However, such extrapolations can be uncertain and may benefit from inclusion of evidence from clinical observations via parametric inference. When clinical interindividual variability is high, or the data sparse, it is essential to use a population pharmacokinetics inferential framework to estimate unknown or uncertain parameters. Several approaches are available for that purpose, but their relative advantages for PBPK modeling are unclear. We compare the results obtained using a minimal PBPK model of a canonical theophylline dataset with quasi-random parametric expectation maximization (QRPEM), nonparametric adaptive grid estimation (NPAG), Bayesian Metropolis-Hastings (MH), and Hamiltonian Markov Chain Monte Carlo sampling. QRPEM and NPAG gave consistent population and individual parameter estimates, mostly agreeing with Bayesian estimates. MH simulations ran faster than the others methods, which together had similar performance.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

No single software platform allows practitioners to compare parametric and nonparametric estimates of calibrated physiologically‐based pharmacokinetic (PBPK) model parameters.

WHAT QUESTION DID THIS STUDY ADDRESS?

How do a maximum likelihood parametric (quasi‐random parametric expectation maximization) and nonparametric (nonparametric adaptive grid estimation) algorithms, and two Bayesian numerical methods (Hamiltonian Markov Chain Monte Carlo and MH) compare in results and timing for calibration of a PBPK model?

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

We demonstrate, for the theophylline dataset studied, a strong coherence of the four approaches tested. Computation times differ largely, with MH being fastest, but this is partly implementation dependent.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

Practitioners will be able to test parametric assumptions with a nonparametric method when calibrating PBPK model parameters and will have greater confidence when making important decisions in a drug development process.

INTRODUCTION

Physiologically‐based pharmacokinetic (PBPK) models usually integrate a large number of parameters because of their mechanistic nature. The common practice workflow for PBPK modeling in clinical studies involves in vitro to in vivo extrapolation (IVIVE) of parameter values. However, such extrapolations can be uncertain and may benefit from inclusion of evidence from clinical observations, via parametric inference on some model parameters. When clinical interindividual variability is high, or the data sparse, it is essential to use a population pharmacokinetics framework when conducting such parametric inference. Several approaches and algorithms are available for that purpose, and their relative advantages for PBPK modeling are unclear. A recent article compared the application of two Markov Chain Monte Carlo (MCMC) numerical samplers to a diazepam PBPK model: the Metropolis‐Hastings sampler (MH), and the Hamiltonian MCMC sampler (HMCMC). Both produced very good fits at the individual and population levels. HMCMC outperformed MH in sampling efficiency, due to its almost uncorrelated sampling. Yet, MH was much faster per iteration, and because of that it converged faster and produced more samples per unit time. A previous paper compared NONMEM (using a frequentist approach) and Winbugs (using a Bayesian approach) on the same data set.

Physiologically‐based pharmacokinetic models include species‐specific physiological parameters, drug‐specific parameters, and parameters defining the trial design. It is essential to document distributions or best values for this vast number of parameters (up to several hundreds), and their covariances. The Simcyp Population‐Based PBPK Simulator compound database holds more than 100 drugs for which mean parameter values and, in most cases, interindividual variance (determined via meta‐analysis of the available data , , , ) are available. Similarly, distributional estimates for physiological parameters are available for more than 20 populations (e.g., Chinese, Japanese, North European, pediatric, and the obese). In the Simcyp Simulator, the covariance between parameters is typically modeled via covariate submodels when prior information is available. For several compounds, parameter values are often extrapolated from in vitro measurements. In such cases, simulations may not accurately match the available clinical data, and parametric inference based on the latter may yield better parameter estimates. Population estimation methods are computationally intensive, so unless the number of model parameters is already small, it may be necessary to confine population fitting to a small subset of them. For the rest of the parameters, which have a mechanistic definition, literature values supported by IVIVE have been used. When the number of uncertain parameters is large, sensitivity analysis is recommended to decide which ones are the most influential and should be estimated.

In this work, we compared results obtained using a minimal PBPK model of a classical theophylline dataset with two frequentist asymptotic methods—quasi‐random parametric expectation maximization (QRPEM), , nonparametric adaptive grid estimation (NPAG), both implemented in the Simcyp Simulator—and two Bayesian numerical sampler, MH, and Hamiltonian HMCMC.

INFERENCE METHODS AND ALGORITHMS

Statistical parameter estimation uses frequentist point‐estimation or Bayesian distribution‐based methods. Frequentist methods produce confidence intervals around point‐estimates, which should lead to approximately the same decisions as Bayesian methods, given enough data. An independent classification differentiates parametric (assuming fixed distributions) from nonparametric methods. The methods investigated here cover the parametric and nonparametric frequentist cases, as well as the parametric Bayesian case.

Quasi‐random parametric expectation maximization

The QRPEM algorithm (Figure 1) is a population pharmacokinetic likelihood maximizer first implemented in the Phoenix NLME software and available in Simcyp version 19. In QRPEM, importance sampling of the posteriors is based on quasi‐random (“low discrepancy” or Sobol) sequences with increased sampling efficiency over pseudo‐random sequences. QRPEM also uses sampling importance resampling to draw individual random effects from their conditional empirical Bayes distribution (EBD). The overall marginal log‐likelihood function evaluated in QRPEM is: where represents the observations for individual , the structural parameter vector for individual (). The individual parameters are assumed to be normally or log‐normally distributed. The is a vector of fixed effects (population means), is the population parameters’ variance–covariance matrix, and is the residual error variance. The integrand of Equation (1) is usually written as an expectation with respect to proposal distribution g. The EBD is defined as: where is a normalization factor defined by the integral in (1), is the individual likelihood of the observation , and is the conditional distribution of the individual parameters with population mean and covariance. is approximated by a discrete distribution on support points where the probability on each support point is: with In Simcyp version 20, is a mixture of two multivariate normal distributions. The first distribution is where and are the individual ’s Empirical Bayes estimates of the mean and the covariance of current iteration of the EM algorithm and is a dilation factor. Distribution is wider, and being respectively the current population mean and covariance matrix estimates. Equal numbers of quasi‐random samples for individual parameters are drawn from these two distributions to calculate a current (discrete) EBD via: The QRPEM maximization step is done via optimization of the residual error model parameters with a Newton method or L‐BFGS‐B. A fast bootstrap method is used to calculate the confidence regions for estimated parameters and predictions. Simcyp QRPEM calculates individual estimates as means (eventually in log‐space): Individual covariance matrices using: Population mean: The population variance–covariance matrix is calculated as:

FIGURE 1

Key components of the quasi‐random parametric expectation maximization (QRPEM) algorithm in Simcyp version 19/20

Nonparametric adaptive grid estimation

The NPAG algorithm is a likelihood maximizer that does not rely on distributional assumptions of individual values. The key problem in the nonparametric maximum likelihood (NPML) method is to estimate the optimal distribution function, , that maximizes the likelihood of the observations among the set of discrete distributions with no more support points than the number of subjects in the population. The fact that the number of support points must be limited to is proven by the Carathéodory Theorem.

In the approach proposed by Leary and Burke, the original NPEM algorithm is applied to a modestly sized grid to produce an solution. A primal‐dual interior point algorithm with near‐quadratic convergence is used to solve the NPML estimation problem. Support points with very small probability are deleted from ; new support points are added around each of these remaining support points in leading to a new grid and an improved likelihood (Figure 2). The process is repeated until the convergence of likelihood is achieved.

FIGURE 2

The NPAG algorithm starts from a set of initial grid points obtained with only a few iterations of QRPEM. The grid is updated (removing and adding points) according to probabilities re‐calculated in every iteration of the algorithm, as explained in the text. NPAG, nonparametric adaptive grid estimation; QRPEM, quasi‐random parametric expectation maximization

NPAG results in a conditional probability matrix with one row per subject and one column per final support point , where is the likelihood of subject ’s data given the parameter values of the support point. Given , the primal dual algorithm estimates an optimal probability vector with one element per support point. The population estimate of a parameter is calculated as the expectation over the final support points, , where is the value of the parameter vector at the grid point. Parameter estimates for subject are computed using the discrete empirical Bayes distribution for that subject. The probabilities in that distribution form a vector with elements equal to: The parameter vector estimate for subject is then given by the expectation .

In Simcyp version 19, NPAG was added as a population estimation algorithm. It uses the QRPEM to generate the NPAG initial grid. The population covariance matrix, calculated using individual estimates, and the probability distribution of the final grid are reported with other numerical and graphical performance measures. The coefficients and of the combination error model‐ estimated within the initial QRPEM run (see Equation (15) below) are input to the NPAG algorithm and remain fixed.

Metropolis‐Hastings sampling

The MH algorithm , is an MCMC sampling method. It generates a sequence of random draws from any distribution for which the density can be computed up to a constant. Let be a set of parameters of interest. According to Bayes’ theorem, the joint posterior distribution of given data is proportional to the product of its prior distribution by the data likelihood: To simplify the notation, let be the product . At iteration zero of the sampler, the MH algorithm starts from an arbitrary point (the exponent is used here for indexing). It then samples a “candidate” point using a so‐called “instrumental” (i.e., arbitrary but judiciously chosen) conditional distribution. The conditioning is, in general, on the previously sampled value of (in this case ), hence the appellation “Markov chain.” is often a multivariate normal distribution centered on the previous draw of . The candidate is selected or not, according to the following procedure:

Compute the ratio

If r ≥1, accept ; otherwise, accept it with probability r by sampling a uniform random number u between zero and one and accepting if u ≤ r ;

If is accepted, keep and call it 1; otherwise discard it and make 1 =  0.

The algorithm continues sampling proposed values and accepting or rejecting them, according to the value of , as long as needed. It can be shown that after an infinite number of draws, the retained values of form a random sample from the desired posterior distribution. In practice, it is enough to run several chains from different starting points 0 and check that they approximately converge to the same distribution. Note that if J( ) is symmetrical and centered on the previous value, the ratio J( | ′)/J( ′| ) in Equation 12 is always equal to 1 and does not need to be evaluated.

Instead of sampling the whole parameter vector at once, it is common practice to split it into multiple components using conditional independence arguments and update these components one by one. In the simplest case, each component is a scalar, and its proposal distribution J is univariate. We selected that sampling scheme for our MH computations.

Hamiltonian Markov Chain Monte Carlo sampling

This MCMC sampler avoids random walk. It uses Hamiltonian dynamics to generate proposal values. Imagine a particle with some position and momentum, moving over time in a frictionless space. The model parameters correspond to the position vector and an auxiliary vector represents the momentum vector. The potential and kinetic energy of the particle constitute the total energy (equal to the Hamiltonian in the HMCMC context), which remains constant over time due to lack of friction. The position and momentum of the particle can be calculated by solving Hamilton's equations (a system of partial differential equations). ,

An advantage of HMC is that it informs the proposal about the target distribution through numerical approximation of its gradient. Therefore, HMCMC performs better in high dimensional problems where the posterior has high curvature compared to other MCMC methods. Yet, per step, HMCMC sampling is much heavier computationally.

MATERIALS AND METHODS

Experimental data

A canonical theophylline data set , was used (see Table S1). Eighteen healthy male subjects were administered an immediate release formulation of theophylline. For subjects A–K, 320 mg, for subject L, 268 mg, and for subjects M–R, 190 mg oral doses of theophylline were given, respectively. The concentration (in mg/L) of theophylline in patient’s blood plasma samples was measured. The concentrations at time zero were not used (the model estimates would always be zero). Body weight was recorded for only 12 of the subjects (codes A to L) and was missing for the others.

Physiologically‐based pharmacokinetic model

Physiologically‐based pharmacokinetic model equations

The Simcyp minimal PBPK model consists of three well‐stirred compartments, predicting the systemic, portal vein, and liver concentrations (Figure 3). Using a first order absorption model, input from the gut lumen can be obtained analytically. Because the volume of the portal vein blood is very small compared to the others, we assumed that the portal vein drug concentration is always at equilibrium. For theophylline, we also assumed negligible gut wall metabolism. The model equations are provided in Section S2.

FIGURE 3

The minimal physiologically‐based pharmacokinetic (PBPK) model of the Simcyp Simulator

Individuals’ physiological covariates and parameter values

The subject‐specific values of most physiological parameters were scaled using subjects’ age. Because Trembath and Boobis only give the age range (21‐ to 31‐year‐old) of the subjects, we computed their expected age by randomly sampling the age of 1000 White men in that age range from a Weibull distribution, using Simcyp version 20. The average age of the simulated subjects (27.6 years) was then imputed to the 18 subjects studied here. When body mass was given as data, it was used, rather than computed as a function of age.

The absorption rate constant (; 1/h), the volume of distribution of the drug at steady state (; L/kg) and CYP1A2 liver enzyme mediated maximum rate of metabolism () of theophylline (pmol/min/pmol) were estimated for each subject, because they are the most uncertain or variable and sensitive parameters of the model.

Statistical model and estimation procedures settings

The statistical model we consider assume that the theophylline plasma concentration measurement for the subject was normally distributed with a mean a time t given by the structural PBPK model, , and a residual error variance : where was a vector of three estimated subject‐specific parameters ( V ss and V ), was a vector of fixed parameters (given in Table S2), b 0 and b 1 being residual error model parameters (also estimated).

At the population level, the subject‐specific parameters were assumed to be distributed according to a multivariate log‐normal distribution, with population geometric means (estimated) and population covariance matrix Ω in log‐space:

The residual error variance, , describes modeling and measurement errors and was chosen to be proportional to the prediction with an added constant error:

Quasi‐random parametric expectation maximization estimation

Table S3 gives the initial values of the population means and residual error parameters, the bounds applied to and components, and the coefficients of variation specifying the initial population variances. The covariances (off‐diagonal elements of Ω) were initially set to zero. The initial value of b 0 corresponds to 10% of the average data values.

The QRPEM‐specific control parameters’ settings are given in Table S5. The number of QRPEM iterations was set at 100. For the expectation step, 500 values per subject were randomly drawn from the population distribution. An acceptance ratio (α) was used to set the dilation factor (λ; see QRPEM description above) as: where is the number of random effects in the model. L‐BFGS‐B optimization was used to estimate the residual error model parameters in the maximization step. The SUNDIALS CVODES solver, with the settings given in Table S7, was used to integrate the model’s differential equations.

Nonparametric adaptive grid estimation estimation

For NPAG, parameter initial values and bounds were the same as for QRPEM (Table S4). The NPAG‐specific control parameters’ settings are given in Table S6. The number of iterations was set to 100. The initial NPAG grid needs only to hold approximate solutions, so 10 iterations of initial QRPEM were run, with 200 Monte Carlo samples. The SUNDIALS CVODES solver was used for integration, as above (Table S7).

Metropolis‐Hastings estimation

Because of software limitations, a diagonal population covariance matrix was used instead of a full matrix with MH. Bayesian methods require the specification of prior distributions. At the population level, the priors for the population means of (K , V ss, and V max) were uniform within the same bounds as QRPEM and NPAG (Table S3). The priors for the population geometric SDs of K , V ss, and V max were set to normal distributions with mean 1.1, SD 0.3 and truncation bounds 1 and 4.0 (corresponding to a half‐normal for the population SDs in log‐space). The prior for b 0 was half‐normal with a SD of 5, and the prior for b 1 half‐normal with a SD of 0.1.

Four Markov chains of 10,000 iterations were simulated. The initial 5000 samples were discarded and only one in two of the last 5000 samples were used, forming a final sample of 10,000 posterior parameter vectors. The convergence of the chains was assessed visually and confirmed using Gelman and Rubin potential scale reduction factor ; lower than 1.05 indicates well converged chains. The Livermore Lsode implicit solver was used to integrate the structural model, with relative and absolute tolerances of 10−7 and 10−8, respectively.

Hamiltonian Markov Chain Monte Carlo estimation

The population means and error model parameters b 0 and b 1 were assigned the same priors as for MH (above). The full population covariance matrix prior was specified using a scale vector distributed like the populations SDs in MH, and a Lewandowski–Kurowicka–Joe prior on the correlation matrix, with parameter 3.

The standard HMC no‐U‐turn sampler was used. Four Markov chains were simulated with 500 warmup iterations and 1000 sampling iterations per chain, and no thinning. Convergence was assessed visually, with the criterion, and with the effective number of samples. Backward differentiation formula method (an implicit integrator, like for the other methods) was used for solving differential equations, with relative and absolute tolerances both set at 10−6.

Software used

We used the QRPEM and NPAG algorithms implemented in Simcyp version 20. MH simulations were performed with GNU MCSim , available at www.gnu.org/software/mcsim. HMCMC simulations were performed with the R package rstan under R version 4.0.5. Graphs and additional statistical treatment were performed in R. Software was run on a Dell Precision 5540 Notebook computer, with 16 GB RAM and 6 Intel Core i7‐9850H processors clocked at 2.6 GHz.

RESULTS

Model fits to the data

Figure 4 shows the best (maximum likelihood or maximum posterior) predictions versus observations for each of the four parameter estimation methods used. The predicted plasma concentration versus time curves, overlaid with individual observations, are shown in Figures S1–S4 of the Supplementary Material from QRPEM, NPAG, MH, and HMCMC, respectively. All four methods gave reasonable fits to the data and the residual error is probably close to the data measurement error. There is a significant, even if not very high, intersubject viability, as can be seen by the deviations from the population average time curve.

FIGURE 4

Individual observations vs. predictions by all four parameter estimation methods. The solid line is the unity line; the blue and red dashed lines mark, respectively, 1 and 2 standard deviations above and below line of unity. HMCMC, Hamiltonian Markov chain Monte Carlo; MCMC, Markov Chain Monte Carlo; NPAG, nonparametric adaptive grid estimation; QRPEM, quasi‐random parametric expectation maximization

Comparison of population parameter estimates

Figure 5 shows the population means and SDs for V max, k , and V ss, and the error model parameter estimates for each method. The populations SDs are given on the log‐scale, because the three parametric methods assume a log‐normal distribution of individual values in the population. This is not relevant for NPAG, but we computed the SD estimates after log‐transformation of the parameter grid points coordinates for easy comparison with the results of the other methods. MH and HMCMC yield posterior distributions, QRPEM computes confidence intervals, and NPAG does not provide uncertainty estimates. Yet, overall, the four methods gave similar estimates for population means and SDs, and for the proportional error component (b 1) of the measurement error model. The additive error component (b 0) has a long left‐tail as it is poorly identifiable from the data, and the best QRPEM estimate is lower than the others, whereas still falling in a feasible region of the Bayesian posteriors. Note that the prior bounds imposed on the population means were not reached by any of the methods used and the estimates are practically only constrained by the data in all cases. Histograms of the posterior population parameter samples are shown in Figures S5 and S6 for HMCMC and MH, respectively.

FIGURE 5

Population (mean, SD) and error model (b 0, b 1) parameters posterior distributions (a) HMC; (b) MH or maximum likelihood estimates (c) QRPEM, with 95% confidence intervals when available; (d) NPAG. HMCMC, Hamiltonian Markov chain Monte Carlo; , absorption rate constant; MCMC, Markov Chain Monte Carlo; NPAG, nonparametric adaptive grid estimation; QRPEM, quasi‐random parametric expectation maximization; V max, maximum rate of metabolism; V ss, volume of distribution of the drug at steady state

Figure 6 shows the (best) population distributions densities according to QRPEM, MH, and HMCMC for V max, k , and V ss, overlaid with the NPAG individual estimates. The NPAG estimates do not point to any significant deviation from lognormality of the population distribution, but there are only 18 subjects to check that. Note that the final NPAG grid of support points can have at most 18 points, but in fact only eight have a probability greater than 10−5. The population distributions are consistent across the four methods, even though QRPEM and NPAG population mean estimates for V ss are a bit lower than those of MH and HMCMC. Detailed population parameter estimates are given in Tables S8–S11 for QRPEM, NPAG, HMCMC, and MH, respectively. Detailed individual parameter estimates are given in Tables S12–S15. The last two tables also give the convergence values obtained with HMCMC and MH, respectively. Approximate convergence of the simulated Markov chains was always achieved.

FIGURE 6

Best population distributions density estimates according to QRPEM, MH, and HMCMC for V max, k , and V ss, overlaid with the NPAG individual estimates (points). HMCMC, Hamiltonian Markov chain Monte Carlo; , absorption rate constant; MCMC, Markov Chain Monte Carlo; NPAG, nonparametric adaptive grid estimation; QRPEM, quasi‐random parametric expectation maximization; V max, maximum rate of metabolism; V ss, volume of distribution of the drug at steady state

Computation time

Run times of QRPEM and NPAG were 4.5 h and 1 h 20 min respectively, and 1 h 15 min and 1.5 min for HMCMC and MH, respectively. To check whether the much shorter run time of MH was partly due to a simpler (diagonal) matrix estimation, we ran separately an HMCMC estimation with a diagonal matrix. The HMCMC run time was the same as with a full covariance matrix and the results were the same as those of MH (data not shown), so the better performance of MH is not due to simpler matrix estimation.

DISCUSSION

In this work, we used the parametric QRPEM and the nonparametric NPAG methods, both implemented in Simcyp version 20, together with two Bayesian numerical methods (HMCMC and MH) to calibrate a PBPK model using population data on theophylline plasma concentration following oral administration of an immediate release form. We used a minimal PBPK model consisting of only two ordinary differential equations and a first‐order absorption rate but complemented with a rich covariate model.

We could have used a more complex PBPK model, but such a model was not needed to describe theophylline pharmacokinetics. Theophylline undergoes negligible intestinal metabolism and exhibits almost 100% bioavailability, with a time to peak plasma level of about 45 min in humans given theophylline orally as a solution. , Those justify using a first‐order model for theophylline absorption and estimating an absorption rate. Theophylline distributes mostly into body water. Its apparent volume of distribution is variable across subjects (which justifies estimating it, as we did here) with a mean of about 0.48 L/kg. This low volume of distribution justifies using a minimal PBPK model. The baseline values of our theophylline model parameters were extracted from the Simcyp compound library file for that drug. Metabolism of theophylline by hepatic CYP1A2 accounts for 70–80% of its total metabolism. Out of the three theophylline metabolic pathways mediated by CYP1A2, we selected the main N3‐demethylation pathway and V max for that pathway was fitted, given its relatively high expected interindividual variability. Its value is only apparent and includes contributions from all pathways and enzymes.

Given the unavoidable uncertainty affecting the results, the four methods gave comparable estimates of all model parameters at the population and subject levels. The posterior fit to the data were uniformly very good. Because we used uniform prior distributions for population means in Bayesian estimation, it is no surprise that the results are similar to those of maximum likelihood methods. It is reassuring that the choice of method does not strongly condition the data analysis results. The population estimates of k , are close to 2 h−1, which correspond to an absorption half‐time of about 30 min. The population mean V ss estimate is very close to the published value 0.48 L/kg, and the population mean for V max is close to 10 pmol/min/pmol of CYP1A2, which is compatible with the value used in the Simcyp compound library file (6 pmol/min/pmol) because, as mentioned above, our estimate includes contributions from other metabolic pathways which are modeled separately in Simcyp. Population variability estimates correspond to about 80% coefficient of variation (CV) for k , 20% CV for V ss, and 50% CV for V max. A correlation at the population level between V ss and V max seems to exist, but it is quite uncertain (according to HMCMC, which gives a measure of its uncertainty), and MH (which neglect those correlations) gave estimates consistent with those of the other methods. We suggest that the amount of information and covariance modeling that enters the PBPK model leaves their parameters quite uncorrelated. This has already been noticed by others. , This feature allows the use of diagonal population covariance matrices, simpler to estimate.

In terms of computation time, MH was the fastest, followed by NPAG and HMCMC, and then by QRPEM. However, QRPEM and NPAG are implemented in the Simcyp Simulator and we used only a subset of its minimal PBPK model functionalities. It is likely that some overhead in QRPEM is due to unused capabilities of the model.

Population estimation of PBPK model parameters has been previously reported using various Bayesian or maximum likelihood methods implemented in Gnu MCSim, Matlab, Stan, PK‐Sim, and NONMEM. , , , , The current study is the first to compare the performance of parametric Bayesian to those of the parametric QRPEM and non‐parametric NPAG frequentist methods in the context of PBPK models.

The population parameter estimation module in Simcyp offers both a parametric method (QRPEM) and a nonparametric one (NPAG). This allows the usual (log)normality assumptions made by QRPEM to be immediately checked by NPAG. For theophylline, with 18 subjects, NPAG did not point to any flaw in that assumption. This could indeed be different with other drugs and populations.

Population inference can be used with complex models if they are supported by sufficient individual data. Obviously, computing power should be scaled up accordingly. Complex models often include mechanistic parameters for which informative priors can be defined. In a Bayesian context, that helps the estimation by increasing identifiability. Still, only the parameters to which the data predictions are sensitive should be included to limit identifiability problems. That is partly why global sensitivity analysis was recently implemented in Simcyp.

The possibility to use population inference for complex PBPK model opens interesting perspectives for improved simulations of clinical trials for drug–drug interactions, special populations, or bioequivalence assessments. The ability to estimate population variability in drug‐specific PBPK parameters in Simcyp should help users in capturing the variability seen in available data when running clinical trial simulations. Population means and CV estimates given by QRPEM or NPAG can be input for the parameters of the compound of interest. However, one also needs to consider the variability affecting all the other PBPK model parameters that were kept fixed during the estimation process. For those, informative prior distributions determined by meta‐analysis or IVIVE should be used when running predictive checks of the inference and further simulations. On a more fundamental level, better estimates of population variability for some parameters and processes can point to the need to better understand the biological determinants of that variability.

CONFLICT OF INTEREST

All authors are employees and potential shareholders of Certara.

AUTHOR CONTRIBUTIONS

J.W., A.A., S.C., H.M., F.B., R.L., M.D., R.M., K.A., and M.J. designed the research. F.B., R.L., and M.D. contributed new analytical tools. J.W., A.A., S.C., H.M., F.B., R.L., and M.D. performed the research. J.W., A.A., S.C., H.M., and F.B. analyzed the data. J.W., F.B., A.A., S.C., and H.M. wrote the manuscript.

Appendix S1

Click here for additional data file.