Considerations and Caveats when Applying Global Sensitivity Analysis Methods to Physiologically Based Pharmacokinetic Models

Three global sensitivity analysis (GSA) methods (Morris, Sobol and extended Sobol) are applied to a minimal physiologically based PK (mPBPK) model using three model drugs given orally, namely quinidine, alprazolam, and midazolam. We investigated how correlations among input parameters affect the determination of the key parameters influencing pharmacokinetic (PK) properties of general interest, i.e., the maximal plasma concentration (Cmax) time at which Cmax is reached (Tmax), and area under plasma concentration (AUC). The influential parameters determined by the Morris and Sobol methods (suitable for independent model parameters) were compared to those determined by the extended Sobol method (which considers model parameter correlations). For the three drugs investigated, the Morris method was as informative as the Sobol method. The extended Sobol method identified different sets of influential parameters to Morris and Sobol. These methods overestimated the influence of volume of distribution at steady state (Vss) on AUC24h for quinidine and alprazolam. They also underestimated the effect of volume of liver (Vliver) for all three drugs, the impact of enzyme intrinsic clearance of CYP2C9 and CYP2E1 for quinidine, and that of UGT1A4 abundance for midazolam. Our investigation showed that the interpretation of GSA results is not straightforward. Dismissing existing model parameter correlations, GSA methods such as Morris and Sobol can lead to biased determination of the key parameters for the selected outputs of interest. Decisions regarding parameters’ influence (or otherwise) should be made in light of available knowledge including the model assumptions, GSA method limitations, and inter-correlations between model parameters, particularly in complex models. Graphical abstract Electronic supplementary material The online version of this article (10.1208/s12248-020-00480-x) contains supplementary material, which is available to authorized users.


3) Regression and correlation based method
Regression and Correlation based GSA models, such as standardised regression correlation (SRC), partial correlation coefficient (PCC), and partial ranked correlation coefficient (PRCC), are also popular approaches [6,7]. SRC and PCC were not suitable for non-linear situation or non-monic problems as they were developed based on linear and monotonic assumption [1,3]. PRCC is robust for nonlinear problem, but still depends on a monotonic assumption between inputs and outputs. Hence it will be inaccuracy when non-monotonicities in the model are present [4].

4) Method based on Monte Carlo filtering
Multi-parametric sensitivity analysis (MPSA) is one type of GSA method based on Monte-Carlo filtering, which evaluates the parameter sensitivity based on Kolmogorov-Smirnov statistics by classifying the parameter sets [1]. The overall impact of the input variables on the model output can be investigated by MPSA method. However, the estimation can be subjective as an acceptable threshold and a proper distribution of variables are required, which can be subject to high inter-observer variability.

5) Meta model-based methods
Meta-modelling method is to analyse the impact of input parameter on model output by replacing the original models with statistical or experimental design methods [7,8].

GSA method for models with correlated input variables
Although variance-based method is robust to nonlinear and non-monotonic system, it does assume the input space are not correlated. Recently a few studies have been proposed to extend the variance-based method to cope with correlated input variables [9][10][11][12][13].

1) Extension of regression method
Xu and Gertner developed a variable decomposition method based on regression by considering the correlation among input variables, which works well when there is an approximately linear relationship between response and inputs [13]. Three indices corresponding to uncorrelated, correlated, and total (=uncorrelated + correlated) effect will be derived to assess parameter importance.

Extension of FAST
An extension of FAST was developed by Xu and Gertner [12] to models with correlated variables. The characteristic frequency of a parameter is exploited to capture both the uncertainties of the parameter itself and the dependent variations of other variables. Similar to FAST, only main effect index can be quantified.

Extension of HDMR
Li et al [10] established a new framework of global sensitivity analysis, named as structural and correlative sensitivity analysis (SCSA), for models with independent and/or correlated inputs. SCSA relies on covariance decomposition of the unconditional variance of the output, and split the contribution of an input to structural and correlative parts. Three indices will be derived corresponding to structural, correlative, and total (= structural + correlative) effect respectively. When inputs are independent this proposed framework will reduce to a random sampling-high dimensional model representation (RS-HDMR) method.
Also, Zuniga et al. [14] extended the random sampling (RS) or Quasi Monte-Carlo sampling (QMC)-HDMR methods to coupe with model with dependent variables. This method was claimed to be more efficient than the extended Sobol method by Kucherenko et al. [9] in evaluations of main effect.

Extension of Sobol method (exSobol)
An extended Sobol method was proposed by Most et al. [11] assuming a linear correlation between model input variables. This method divided the contribution of an input to variance into uncorrelated (independent) and correlative parts, hence 4 sensitivity indices would be available, e.g. uncorrelated and correlated firstorder and total indices. Similarly, Mara et al. [15] developed a new variance-based sensitivity analysis method based on the Gram-Schmidt decorrelation procedure and the polynomial chaos expansion (PCE). A set of new sensitivity indices, called full sensitivity indices, can be derived using the new set of independent variables after decorrelating the input variables [16]. Besides, Kucherenko et al. [9] developed a copula-based method to calculate the main (first-order) and total effect analogous to standard Sobol indices. This method would need a priori knowledge of parameter probability distribution, but does not use surrogate models, data-fitting procedures or orthogonalization of the input factor space.

Test function 1 -Ishigami-Homma function
Ishigami-Homma function is a non-linear and non-monic function commonly used in the Global Sensitivity Analysis [6].
The analytical solutions for first and total order indices are  occur, due to the random sampling for Monte-Carlo estimation, which was also observed in other studies [18].

Test function 2 -Sobol g-function
Sobol g-function is another function widely used in the Global Sensitivity Analysis, particularly due to its strong nonlinearity and non-monotonicity [6]: Sample standard deviations ( ) where ai are non-negative parameters, and Xi ~ [0, 1], for all i = 1, …, k.

GSA analysis for model with correlated parameter
Simple linear function was adopted from Kucherenko et al. [9] to evaluate the impact of correlation on estimation of sensitivity indices. Input variables xi have normal distribution with mean μ=0 and covariance matrix as where ρ is the correlation coefficient, and σ is the standard deviation. The analytical solution for first and total order sensitivity indices were given as Figure S4. Comparison of estimated sensitivity indices to analytical solutions for test function 3: σ=2 (a and b); Correlation coefficient ρ23=-0.5, σ=2 (c and d).
Overall, extended Sobol recovered both main and total effects well for test function 3 (Fig. S4). However, the accuracy of Sobol, developed for models with non-correlated inputs, will be compromised by the presence of correlation among input parameters.

Test function 4 -Ishigami-Homma function
The non-linear and non-monotonic Ishigami-Homma function in section 1.1 was adopted for validation of , and a=7 and b=0.1 [9]. The estimated first order and total order sensitivity indices agree well with results presented by Kucherenko et al. [9] with random sample size N = 2 13 (Fig. S5). Minor differences would be expected due to the random sampling of input parameters.   For a total of r sample points, the Morris indices for the parameter Xi are defined as:

Results of Morris method for Quinidine, Alprazolam, and Midazolam
A high μ or μ* indicates a parameter with an important overall influence on the model outputs; a high σ indicates either a parameter interacting with other factors or that its effect is non-linear. The magnitudes of μ and σ for each model parameter are relative to the others and have no absolute meaning [21].

Sobol method
Sobol method is a variance-based type GSA method, which decompose the variance of the model outputs into sums of variances for combinations of input parameters of increasing dimensionality [2]. Although there is no assumption about the relationship between the model inputs and outputs in variance-based GSA methods, they do assume that the input parameters are independent. Generally, when using Sobol, three sensitivity indices are calculated to determine the importance of input parameters: -A first-order "main effect" sensitivity index evaluating only the main influence of each parameter without considering the interaction with others; -A "total-effect" sensitivity index to assess the impact of each parameter including all possible interactions with others; -An "interaction" index, which is the difference between total effect and main effect, representing only the contribution of parameters interactions.

Extended Sobol method
The GSA method proposed by Kucherenko et al. [9] can consider models where the input parameters are correlated. The main (first-order) and total effect sensitivity indices, analogous to standard Sobol indices, are calculated using a copula-based method. Briefly, for a function f(x1,..., xn) where the random vector z is a sample from the conditional distribution p(z|y) and E is the expectation function over f with respect with its subscript. The corresponding first-order effect index of the subset y is defined as: While the total-effect index of the subset y is: Collectively, and in the case of independent variables are called Sobol's Indices.
The first order effect Sy can be estimated as [9]: Here p(y) is the marginal distribution function of y, the random vector (y, ′ ) is a second random vector drawn from the conditional probability distribution p(y, | ), and (y, z) is a sample from the joint distribution p(y, z).
The total order index is given as: Correlation analysis can be adopted to provide a priori knowledge of probability distribution function p(y, z) required by the extended Sobol method [16]. If a linear correlation between two input parameters is assumed, say, Xi and X j, the correlation coefficient, ρi j, is calculated as follows: The covariance Cov(Xi, Xj) between Xi, X j is: where ̅ , ̅ and σXi, σXj are means and standard deviations of Xi and Xj, respectively.
If the correlation matrix is the identity matrix (hence if there are no correlations), the extended Sobol will be reduced to a simple Sobol method using the Jensen-Sobol' formula [9]. I.e., the conditional distributions in   Most's method • Linear correlation is assumed among model input parameters.
• The contribution of an input to variance will be split into uncorrelated (independent) and correlative parts.