A two-stages global sensitivity analysis by using the δ sensitivity index in presence of correlated inputs: application on a tumor growth inhibition model based on the dynamic energy budget theory

Global sensitivity analysis (GSA) evaluates the impact of variability and/or uncertainty of the model parameters on given model outputs. GSA is useful for assessing the quality of Pharmacometric model inference. Indeed, model parameters can be affected by high (estimation) uncertainty due to the sparsity of data. Independence between model parameters is a common assumption of GSA methods. However, ignoring (known) correlations between parameters may alter model predictions and, then, GSA results. To address this issue, a novel two-stages GSA technique based on the δ index, which is well-defined also in presence of correlated parameters, is here proposed. In the first step, statistical dependencies are neglected to identify parameters exerting causal effects. Correlations are introduced in the second step to consider the real distribution of the model output and investigate also the ‘indirect’ effects due to the correlation structure. The proposed two-stages GSA strategy was applied, as case study, to a preclinical tumor-in-host-growth inhibition model based on the Dynamic Energy Budget theory. The aim is to evaluate the impact of the model parameter estimate uncertainty (including correlations) on key model-derived metrics: the drug threshold concentration for tumor eradication, the tumor volume doubling time and a new index evaluating the drug efficacy-toxicity trade-off. This approach allowed to rank parameters according to their impact on the output, discerning whether a parameter mainly exerts a causal or ‘indirect’ effect. Thus, it was possible to identify uncertainties that should be necessarily reduced to obtain robust predictions for the outputs of interest. Supplementary Information The online version contains supplementary material available at 10.1007/s10928-023-09872-w.


Introduction
The development process of anticancer agents is characterized by a high attrition rate [1,2]. Although the promising results obtained during preclinical assessment, the majority of investigated compounds are withdrawn due to inadequate safety and/or efficacy in late clinical phases [3,4]. For these reasons, taking full advantage of the information coming from preclinical studies in order to reduce the gap between the preclinical and clinical settings and to anticipate the outcomes of the clinical trials is essential [1].
Pharmacometric models can be powerful tools to summarize and quantitatively integrate information collected during the preclinical anticancer drug development, generally based on xenograft experiments, and to support their translation to the clinical setting [5][6][7][8][9]. The success of this model-oriented approach is witnessed by the large number of pharmacokinetic-pharmacodynamic (PK-PD) models describing the Tumor Growth Inhibition (TGI) after anticancer treatment in xenografts [5,[10][11][12][13]. PK-PD TGI models quantitatively relate plasma concentration-time curves to the tumor growth-time curves providing highly valuable metrics of tumor growth and anticancer drug activity [14,15]. These model-derived metrics should be independent from the experimental settings and, therefore, can be successfully used to compare drug candidates, to perform preclinical-to-clinical extrapolations [16] and to anticipate the efficacious exposure to be targeted in clinics, supporting the choice of clinical efficacious doses [17,18]. A novel mathematical model describing the tumor and host body weight dynamics during an anticancer drug treatment in xenograft mice has been recently proposed [15,19,20]. This model is based on the Dynamic Energy Budget (DEB) theory [21][22][23]. This tumor-in-host DEB-TGI model is able to integrate several aspects characterizing the in vivo TGI studies: anticancer drug activity on tumor, onset of drug-related and tumor-related cachexia and anorexia and influence of host conditions on tumor growth [24]. The incorporation of such aspects in a single model allows to investigate body weight loss due to tumor progression and treatment and, at the same time, to obtain unbiased estimates of drug anticancer efficacy. The DEBbased modeling approach was successfully applied on several in vivo preclinical studies that had been performed to assess the efficacy of the investigated compounds and that involved different host species (mice and rats), tumor cell lines, type of anticancer agents and experimental settings, including combination regimens [25]. In addition, metrics with a relevant biological meaning were derived through a mathematical analysis of the tumor-in-host DEB-TGI model in the specific formulation defined for cytotoxic agents [15]. In particular, tumor volume doubling time (TVDT), i.e., the time required by a tumor to reach a twofold volume [26] and minimum threshold concentration necessary for tumor eradication, C T , can be computed as functions of the model parameters. Interestingly, the values of TVDT and C T estimated on xenograft mice resulted predictive of the clinical settings [15] and, therefore, can be used to support the preclinical-to-clinical translation.
Uncertainty and variability are two aspects that significantly impact parameters of pharmacometric models. In the context of the DEB-TGI model, uncertainty typically refers to the estimation error introduced by the model identification procedure. This error can be particularly relevant when the experimental design and the sampling schedule are not expressly designed and optimized for model identification, thus, causing some parameters be estimated with low precision [15]. Variability is in general due to different causes. For example, in xenograft experiments, model parameters can vary among different tumor cell lines (i.e., inter-line variability), drugs (i.e., inter-drug variability) or animal species. Whatever the cause, the variation of the parameters of the model is eventually propagated to the model outputs [27]. Therefore, assessing the impact of the parameter uncertainty/variability on the final predictions of the model is a crucial step to evaluate the quality of the model-based inference. However, model parameters are generally not independent each other but they are linked by statistical dependencies (i.e., a correlation structure). Consequently, the variation of a parameter is linked to the variation of the others, according to their correlations. Thus, in this evaluation process, it is essential to take into account also the correlations between the model parameters in order to obtain a reliable analysis [28,29].
Sensitivity Analysis (SA) and Global sensitivity analysis (GSA) are techniques aiming to investigate how the variation of the model parameters influences the model output predictions [27,30,31]. Both US Food and Drug Administration (FDA) and European Medicine Agency (EMA) stressed the importance of SA and GSA to support the Pharmacometrics modelling [8,9,32,33]. GSA is a branch of SA which relies on the multivariate variation of all the considered input parameters [27]. GSA can be defined as ''the study of how uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model input'' [27]. GSA allows to rank inputs (referred to as 'factors') according to their impact on the output variation and to identify the model parameters whose uncertainty and variability should be reduced in order to obtain more reliable model predictions.
Several GSA methods have been proposed in the literature [34,35]. Among them, the variance-based approach [36], which is based on the decomposition of the variance of the model output ; is one of the most well-established and widely used [37]. This approach is able to rank model inputs (i.e., parameters) based on their impact on model output independently on the shape of the input-output relationship, accounting also for not monotonous and nonlinear relationships. Variance-based GSA relies on the assumption of independence between model inputs [27,36]. In such case, the variance decomposition is unique and reflects the structure of the model itself. When this hypothesis is violated (i.e., in presence of correlated inputs), the variance decomposition of model output is not unique [34]. Different extensions of the variance-based GSA method have been proposed to deal with dependent inputs [38][39][40], but the meaning of their results is not easy to interpret [38,[41][42][43].
Despite the challenges, correlations between model inputs cannot be ignored because they may dramatically alter the model predictions and, consequently, the results of SA and GSA [31,40,44]. In presence of statistical dependencies between inputs, the moment-independent GSA methods [31,34] can provide a good compromise between the interpretability of GSA results and the consistency with theoretical assumptions. These techniques consider the entire distribution of the output rather than a single statistical moment (e.g., variance). In particular, the d sensitivity index [45] describes the impact of each model parameter on the output probability density function. This metric is well-defined also in presence of statistical dependencies between the model inputs. Moreover, being a moment-independent technique, d sensitivity index provides robust results, independently from the shape of the output distribution. Conversely, it was reported that variance-based methods can be misleading when the output distribution is multi-modal [46] or highly-skewed [47], since variance is a sensible measure of the output variation [48].
In this paper, a two-stages approach for GSA exploiting the d sensitivity index to deal with correlated inputs is presented. It was then applied to DEB-TGI model to understand the impact of the parameter estimate uncertainty on the three model-derived metrics: C T , TVDT, and an additional one accounting for the trade-off between drug toxicity and efficacy. Relevant considerations with practical implications were derived from GSA results.

Description of the model structure
The tumor-in-host DEB-based model [15,20,24] describes the growth dynamics of the tumor and host body weight according to the DEB theory [21][22][23]. Briefly, it is assumed that the host body is made of two components: the energy reserves, e(t), and the structural biomass, V(t). The energy and structural biomass dynamics follow from energetic balances between the main physiological processes, such as the assimilation, governed by the food-supply coefficientq b , the energetic consumption and the maintenance and growth of structural biomass which costs are given by the m and g parameters, respectively. The tumor is conceived as an additional component able to appropriate a fraction of the host energy and to use this energy to sustain its maintenance,m u , and growth,g u , costs. The amount of energy subtracted by the tumor depends on the tumor gluttony,l u , a measure of the tumor aggressiveness. As tumor exploits host energy resource, the organism has to reduce its growth rate and, eventually, to degrade its structural biomass (tumor-related cachexia). The degradation rate increases until a certain time, t dVmax , when a maximum value,d Vmax , is reached. Finally, the tumor-related anorexia is modelled assuming that tumor progression inhibits host assimilation (IV u50 parameter).
For cytotoxic agents, the model describes two effects: drug cytotoxicity and drug-related anorexia. Drug cytotoxicity is modelled by assuming that a fraction of tumor cells, proportional to the drug concentration through the anticancer drug potency, k 2 , is hit by the drug and becomes not-proliferating. These cells enter into a mortality chain that is governed by the rate k 1 . In addition, an inhibitory effect (Imax model) on host assimilation is included to account for the drug-related anorexia. The half maximal inhibitory concentration, IC 50 , represents a measurement of drug toxicity on host body.
The complete mathematical formulation of the tumorin-host DEB-TGI model can be found in [15,20], while the comprehensive list of model parameter is reported in Table 1.

Analytical model outputs: TVDT and C T
A mathematical analysis of the tumor-in-host DEB-TGI model allowed to define some relevant metrics quantifying tumor-in-host growth and drug anticancer activity [15]. These model-derived metrics were expressed as functions of the model parameters. Consequently, they can be easily computed once the model has been identified on experimental data without the need of performing model simulations.
First, the exponential tumor growth rate, k, that governs the initial tumor dynamics is given by a combination of host-related and tumor-related parameters, i.e., the maintenance and growth costs of both normal and tumor cells, m; g and m u ; g u respectively, along with the tumor gluttony,l u : From the exponential tumor growth rate, a model-based definition of the TVDT can be immediately derived: Regarding the anticancer drug activity, the model postulates the existence of a concentration threshold for tumor eradication, C T (Eq. 3).
This parameter represents the minimal (constant) concentration level beyond which tumor eradication can be asymptotically reached in xenograft mice. Consequently, it can be considered as reference concentration to be targeted in vivo for achieving a significant anticancer activity. However, keeping the drug concentration level beyond C t is a necessary but not sufficient condition for the eradication of tumor. Indeed, the animal could die before reaching the complete tumor eradication.

Simulation-dependent model output: DT dVmax
Let us introduce an additional secondary parameter for the DEB-TGI model: the time delay, DT dVmax , between control and treated individuals in reaching the maximum degradation rate of structural biomass, d Vmax : DT dVmax , which is graphically represented in Fig. 1, is defined as where T dVmax;C and T dVmax;T represent the time at which maximum degradation rate is reached for the first time in control (untreated) and treated individuals, respectively. DT dVmax allows to evaluate the trade-off between toxicity and efficacy of the administered drug. The tumor exploits the host resources for its growth, inducing the degradation of the structural biomass and, consequently, the body weight loss. This tumor-related cachexia is contrasted by the drug treatment that inhibits tumor growth. However, due to its toxic effect, the drug concurrently reduces the host assimilation (drug-related anorexia), thus, worsening the degradation of host structural biomass (drug-related cachexia). Therefore, a DT dVmax [ 0 underlines a toxic effect by the pharmacological treatment on the host organism as the maximum degradation rate of structural biomass is reached earlier in treated individuals. Conversely, a negative value of this index suggests that the drug inhibits tumor growth without inducing a strong reduction of the host assimilation. Differently from C T and TVDT, DT dVmax has not a close analytical expression and must be calculated through model simulations. Given a drug administration schedule and model parameters values, a model simulation has to be performed both in presence and absence of drug treatment to compute T dVmax;C and T dVmax;T . In this work, NONMEM software tool was exploited for performing simulations. The set-up is summarized in the Table S1.1 of the Supplementary Material S1 which also contains the NONMEM code.
Let us consider a generic model where Y is the scalar output of the model, g represents the inputs-output relationship and X ¼ X 1 ; . . .; X N f gis the R N vector of the input parameters. Within the GSA framework, X is considered as a random variable [27] which is characterized by a joint probability density function (pdf); f X X ð Þ: Therefore, Y is a random variable with a pdf, f Y ðYÞ, which can be calculated through Eq. (5) using input samples extracted from f X X ð Þ: The definition of the d index relies on the following considerations [45]. Suppose that one input X i can be fixed to a certain value x Ã i , then, the conditional pdf of Y given As illustrated in Fig. 2, s X i ð Þ represents the difference of the area underlying f Y ðYÞ and f Y ðYjX i Þ, which corresponds to the impact of fixing X i to x Ã i on f Y ðYÞ. X i is a random variable typically assuming more values than just x Ã i . The d sensitivity index for X i can be computed through the expected value of s X i ð Þ over the entire domain of X i as in Eq. (7), where f X i ðX i Þ is the marginal pdf of X i .
It has been demonstrated in [45] that 0 d i 1, in  The difference s(X i ) between these two distributions is represented by the green shaded area and X i is uncorrelated from the other X j ; with i 6 ¼ j A full description of the properties of this sensitivity index can be found in [45].

Two-stages GSA approach
Although the d sensitivity index is well-defined also in presence of correlations among the model inputs [45], the interpretation of the GSA results remains challenging.
As an example, let us consider a simple linear model, 1Þ and let X 4 be correlated to X 1 with q 1;4 ¼ 0:9: From the equation defining Y, it is evident that (1) X 1 and X 2 are the most impactful terms on Y due to their high coefficients; (2) the spurious term X 4 does not exert a causal effect because it does not directly appear in the output equation and (3) X 4 can exert an 'indirect' effect on Y due to its correlation with X 1 . Now, let us perform a GSA and compute the d indices for the four model inputs first neglecting and then, accounting for the correlation between X 1 and X 4 . As shown in Figure S2. 1 and S2.2 in the Supplementary Material S2, the d indices for X 1 ; X 2 and X 3 assume the same positive value in both the cases with d X 1 = d X 2 [d X 3 : Differently, the value of the d index for X 4 varies in the two situations. When q 1;4 is neglected, d X 4 ¼ 0 and, therefore, it is impossible to quantify the 'indirect' impact of X 4 on Y. Instead, when q 1;4 is considered, d X 4 assumes a positive value higher than d X 3 and comparable to d X 1 and d X 2 , so that, in this case, it is difficult to distinguish which model inputs have a causal impact on the output.
From this simple example, it is clear that, although considering the correlations between the model parameters is essential for a correct evaluation of their impact on the model output [28,40], it could make hard the identification of the model parameters exerting a causal effect. Therefore, it is necessary to perform also the GSA without considering the statistical dependencies between model inputs as it allows to easily detect parameters with a direct impact [49].
Based on these considerations, a two-stages approach to perform a GSA with d index in presence of correlated inputs is here proposed. The analysis is composed by two subsequent steps. In Step 1, a first set of d indices, d 1;i ; is computed under the hypothesis of independence between model parameters so that the parameters exerting a causal effect on f Y ðYÞ can be easily identified. In Step 2, a new set of indices, d 2;i , is computed considering the statistical dependencies between model inputs. In this way, the 'correct' f Y ðYÞ is considered and the 'indirect' contributions of each X i on the output distribution due to the correlations with the other parameters is accounted for. At the end, each parameter X i of the model is characterized by two indices, d 1;i and d 2;i . From their values it is possible to establish that: GSA of the DEB-TGI model

Distribution of model parameters
Within the GSA framework, the model parameters are considered as random variables characterized by a joint pdf accounting for their variability/uncertainty. The tumor-in-host DEB-TGI model is generally identified in subsequent steps [15]. First the tumor-related parameters are estimated on healthy (tumor-free) mice growth data; then, the tumor-related, cachexia-related and drug-related parameters are simultaneously identified on tumor and host body weight data from xenograft animals keeping fixed the parameters related to the host and to the PK model.
The aim of this work is to evaluate the impact of the parameter uncertainty coming from the model identification on experimental data from xenograft studies. Thus, it was necessary to define a joint pdf only for the model parameters that are typically estimated on xenograft data (Table 2). Differently, the parameters related to the host and to the PK model, as well as d Vu ; e 0 and x, that were typically fixed to the values reported in the Table S3. 1 in the Supplementary Material S3, were not considered. In particular, a two-compartment model with linear elimination was adopted to describe drug PK.
Due to their biological meaning (see Table 1), all the parameters can assume only positive values. Thus, a multivariate lognormal distribution, LogNðh; XÞ with mean h and variance-covariance matrix X, was used to characterize the joint pdf of the considered model parameters.
The typical values, h; were taken from the parameter estimates obtained by fitting the model on preclinical data [15] and reported in Table 2. Because in the model estimation, the food-supply coefficient q b ; is parametrized as R b was considered in the GSA.
In the first stage of the GSA, parameter variances were selected based on the coefficients of variation (CVs) of the estimates uncertainty and no correlation was considered [15,19]. In particular, two scenarios were analysed to account for situations of accurate (CV1, lower values) and of less accurate (CV2, higher values) estimates.
In the second stage, correlations between model parameters were introduced based on typical correlation matrices obtained during model identification (internal data, not shown): parameter couples were classified in significantly correlated (with corr = ± 0.95, 0.75 or 0.6), uncorrelated (corr = 0) or with an unknown correlation, as illustrated in Figure S3. 1 in the Supplementary Material S3. In particular, the last class included parameters couples for which correlations spanned over a wide range of values depending on the considerd experimental data.
The R package mvLognCorrEst [50] was used to estimate unknown correlations. The correlation structures obtained for the CV1 and CV2 scenarios are reported in Fig. 3 and Figure S3. 2 of the Supplementary Material S3, respectively. In particular, the Supplementary Material S3 provides a detailed description of how the correlation structure was built in the CV2 case.

Preliminary uncertainty analysis
A preliminary uncertainty analysis [27] was performed for each model output, Y, using a Monte Carlo simulation approach [34]. The uncertainty analysis was executed both in presence and absence of correlations between the model parameters. This was done to assess the changes in f Y ðYÞ when statistical dependencies are neglected. 100,000 samples were extracted from the joint lognormal distribution of X for both CV1 and CV2 scenarios. Then, for each extracted parameters set, the output of interest was computed. Eqs. (2) and (3) were used for TVDT and C T respectively, while DT dVmax , was calculated as explained in the section ''Simulation-dependent model output: DT dVmax ''.

Implementation of the two-stages GSA
After the preliminary uncertainty analysis, the Two-stages GSA approach, here proposed, was implemented. A numerical ''given data'' approach to estimate the d sensitivity index has been proposed by Plischke et al. [51]. This strategy leverages Eq. (7) to compute the d index and the kernel density estimation to characterize both the pdf of the output, f Y ðYÞ, and the conditional pdf of the output given an input X i ; f Y YjX i ð Þ. In this work, the MATLAB implementation of the numerical approach presented in [51] and shared by the authors was used.
For each model output, variability scenario and step of the GSA, 100,000 samples of model parameters were extracted from the joint lognormal distribution of inputs (with or without correlations depending on the step of the GSA) using a pseudorandom strategy (rand MATLAB function [52]). For each parameter set the corresponding value of the output was computed as previously detailed. All the input-output sets were used to calculate sensitivity indices. Uncertainty of the sensitivity indices was calculated with 1000 bootstrap samples [53].
Moreover, in each step of the GSA, a spurious input (i.e., not belonging to the DEB-TGI parameters and uncorrelated from them) extracted from a N 0; 1 ð Þ and named 'Noise', was included in the analysis. d Noise was introduced as negative control to provide a quantitative evaluation of the noise/bias introduced in the sensitivity indices estimation by their numerical computation [51]. The 97.5th percentile of d Noise was used as threshold value for evaluating whether d i can be considered different from zero: in case the median d i for X i was lower than the 97.5th percentile of d Noise , d i was considered equal to zero.

Results for C T and TVDT
In Fig. 4 and Figure S4. 1 of the Supplementary Material S4 the uncertainty distributions of C T and TVDT are reported. From these figures, it is possible to appreciate the dramatic impact that ignoring the parameter correlations The results of the proposed two-stages GSA are summarized in Fig. 5.
Panels A and B show d indices obtained for C T in the CV1 and CV2 scenarios, respectively. In the CV1 scenario, in absence of correlations (Step 1), the drug potency, k 2 , resulted the parameter with the highest causal impact on the C T distribution, followed by moderate effect of the tumor-related parameters l u and g u . As expected, d Vmax , R b , V u10 , W 0 , IV u50 , IC 50 and k 1 , that do not appear in the C T definition (Eq. (3)), did not have a direct impact (d 1 =0) on the model output. However, when the correlation is considered (Step 2) IC 50 gained an 'indirect' impact due to its correlation with k 2 (Fig. 3). Interestingly, although m u appears in the definition of C T and is affected by a high uncertainty (Table 2), its d was equal to 0 in both steps. In the CV2 scenario (panel B of Fig. 5), the tumor growth related parameters l u and g u obtained the highest d 1 indices. However, their importance drastically decreased in Step 2 (d 2;l u = 0) when the correlations smoothed their effect on the model output. Conversely, k 2 , with d 1;k 2 [ 0 and the highest d 2 ; remained the most important parameter. Finally, the 'indirect' impact of the IC 50 was still present in the CV2 scenario as well as the moderate effect of m u .
Results for the TVDT are illustrated in panels C and D of Fig. 5. l u and g u had the highest direct impact in both uncertainty scenarios. Again, although m u appears in the TVDT definition (Eq. (2)), its direct contribution was not relevant, in particular in the CV2 scenario where d 1;m u ¼ 0. Again, d Vmax , R b , V u10 , W 0 , IV u50 , k 2 , IC 50 and k 1 had a d 1 =0 because they did not appear in the equation of this metric. Among them, in both uncertainty case, d Vmax was  the parameter with the highest 'indirect' effect due to its correlations.

Results for DT d Vmax
The results of the uncertainty analysis are summarized in Fig. 6 where the distribution of DT d Vmax is represented through box plots and histograms, respectively. As already observed for C T and TVDT, also for DT d Vmax neglecting correlations between the parameters led to an overestimation of the output uncertainty, especially in the CV2 scenario. In the CV1 scenario the 95% CI of DT d Vmax was [-0.10, 1.90] days when excluding and [-0.10, 0.90] when including correlations, highlighting only a moderate overestimation of the variability in the first case. Conversely, in CV2 scenario (high uncertainty) the distribution of the DT d Vmax became more skewed and its 95% CI moved from [-0.10, 3.10] to [0.00, 28.70] days when correlations were ignored. Differently to C T and TVDT, DT d Vmax range approximately remained within two folds of its mean only in CV1 scenario. This aspect made even more relevant performing a GSA to assess which uncertainty sources need to be reduced to obtain a more reliable estimation of the DT d Vmax : Thus, the two-stages GSA was applied to the DT d Vmax in both CV1 and CV2 case. The obtained results are reported in Fig. 7.
As the DT dVmax has no analytical definition, it is impossible to predict a priori which model parameters will not have or not a direct impact on the distribution of this metric. In the CV1 scenario (panel A of Fig. 7), l u and g u had the highest direct effect on the pdf ofDT dVmax , followed by the initial tumor volumeV u10 . Differently, IV u50 ; IC 50 and k 1 had d 1 indices equal to 0 as their median d were lower than the 97.5th percentile of negative controld 1;Noise . Therefore, it was assessed that they had no causal impact on this output. Among them, only k 1 had a d 2 =0, indicating that also its 'indirect' contribution was not relevant. The parameters with the highest d 2 indices resultedV u10 ,g u and l u : As they also obtained a d 1 [0, it was possible to conclude that they had both a direct and 'indirect' impact on the DT dVmax distribution : In the scenario with higher CVs (panel B of Fig. 7), tumor growth related parameters, l u and g u ; had the highest causal effects. All the other parameters, except for the V u10 ; had d 1 indices 0. From the Step 2, it emerged that, even if the effects of l u and g u were confirmed, the parameters with the highest impact due to correlations wered Vmax ,R b ,W 0 , m u : Differently, IC 50 and k 1 were not impactful even when considering the statistical dependencies. To understand the lack of effect exerted by IC 50 on the DT dVmax metric, it is important to observe that  (Table 2) is three orders of magnitude smaller than the concentration level observed for the drug at the administered dose schedule (Table S1.1 of Supplementary Material S1). This implies that the inhibition of pharmacological treatment on food intake is always maximum for all the sampled values of IC 50 as shown by Fig. 8.

Discussions
Pharmacometric models play a key role within the process of anticancer drug development as they can help to inform future clinical trials (e.g., anticipation of the first in-human dose) starting from preclinical data [5,6,16,18]. Therefore, it is essential to investigate the reliability of the model-based inference. Pharmacometric model parameters can be affected by different sources of uncertainty and/or variability. For example, high uncertainty can arise from the model identification on sparse experimental data. This uncertainty can be propagated to the model outputs and predictions. In this context, uncertainty analysis helps in understanding whether model predictions are too uncertain to be useful [30]. In addition, GSA is strongly recommended [30,31] as it can provide useful insights on what are the parameters mostly responsible for the model output uncertainty.
The majority of GSA techniques relies on the assumption of independent model inputs [27,36]. However, this hypothesis can be a strong limitation in many situations. For example, when the GSA scope is to investigate the impact of the estimation uncertainty, it is necessary to consider the correlation related to the estimation process. Ignoring known correlations between model inputs may dramatically alter the model simulations and, consequently, the GSA results [28,31,40,44]. However, the field of GSA in presence of correlated inputs is still embryonic [31,[39][40][41].  In this paper, a novel framework for performing GSA in presence of correlations between model input parameters was proposed based on the d sensitivity index. The choice of d index was done as this index is well defined also in presence of statistical dependencies between model inputs and provides an easier interpretability of the results [45]. This new method was applied on the tumor-in-host DEB-TGI model [15,20,24] with the aim of evaluating the impact of estimate uncertainty on some relevant modelderived metrics. In particular, the focus was on the concentration threshold for tumor eradication, C T ; and the tumor volume doubling time, TVDT, (Eqs. (2) and (3)) derived by the stability analysis of the model [15]. In addition, the DT dVmax , introduced in this paper, was considered. The analysis was performed in two different scenarios of uncertainty, CV1 (''low uncertainty'') and CV2 (''high uncertainty'').
First Indeed, neglecting the statistical dependencies between model parameters led to an overestimation of the model outputs uncertainty, especially when the inputs were sampled from a very wide space (i.e., case of high parameter uncertainty, CV2 scenario). Differently, the presence of correlations imposed some constraints on the sampled parameters, preventing the extraction of implausible parameter sets and, consequently, reducing the presence of outliers in the distribution of the outputs. The uncertainty analysis allowed to assess that, in the realistic scenario of correlated parameters, the estimation uncertainty did not critically affect C T and TVDT that approximatively remained within the two folds of their respective means, even in the CV2 scenario. The DT dVmax metric resulted more sensitive to the parameter uncertainty, however its variation remained limited at least for moderate CVs of parameters (CV1 scenario).
After this preliminary step, the GSA workflow introduced in this paper was applied to the three outputs of interests. The proposed approach is composed of two steps. In the first one, d indices are computed neglecting correlations between parameters to quantify the causal effects of the model parameters on the outputs. In the second step, a new set of d indices is calculated to taking into account also the 'indirect' effects due to the correlations and the real distribution of the output. By comparing the d values obtained in Step 1 and Step 2, it was possible to evaluate how much the presence of correlations can alter the impact of the parameters on the model outputs and, thus, change the results of the GSA.
A parameter with a low causal effect on the model output (Step 1) can gain relevance due to its correlations with the other model inputs (Step 2). The extreme example of this situation is represented by parameters with d 1 = 0 and d 2 [0, i.e., parameters exerting only an indirect effect due to their correlation with other model inputs: this is the case of IC 50 with respect to the model output C T (Fig. 5,  panels A and B). This was trivial as IC 50 does not appear in the C T definition (Eq. (3)). However, IC 50 is strongly correlated with the drug potency k 2 which, in turn, exerts a strong direct effect on C T in both the uncertainty scenarios. For this reason, in the Step 2 of the two-stage GSA, the IC 50 gained relevance emerging as the second most impactful parameter. Such observation has to be carefully considered in the strategies aiming to reduce the estimation uncertainty ofC T . An accurate estimate of k 2 , the most directly impacting parameter on C T ; may not be enough to decrease the output uncertainty, but it is also necessary to reduce the statistical dependencies of those parameters highly correlated withk 2 , such as IC 50 : For example in this specific case, an experimental design that includes a group of treated tumor-free animals could help to statistically disentangle the two drug-related parameters.
Conversely, the presence of statistical dependencies can mitigate the causal effect of some model parameters. This is the case of inputs with high d 1 index and significant lower d 2 . An example is provided by the parameters l u and g u . They had a high direct effect (i.e., high d 1;l u ; d 1;g u ) on C T , TVDT and DT dVmax in the two uncertainty scenarios (Figs. 5 and 6), however, when correlations were considered (Step 2), their impact was remarkably reduced. This example further emphasizes that correlations must be considered in the GSA framework to avoid a wrong ranking of the parameters.
The proposed two-stages approach represents and appropriate methodology for ranking correlated input parameters according to their impact on the model outputs. Furthermore, it helps to discern and quantify the sources of their effects (i.e., direct or indirect). In particular, the obtained ranking identifies those parameters that mostly affect the output uncertainty and consequently require a better estimation (e.g., by using optimal design strategies [54]) to improve the precision of the model predictions. At the same time, parameters with an irrelevant influence on the outputs of interest can be recognized. Assessing the robustness of the model predictions with respect to highly uncertain parameters is extremely relevant when these parameters are affected by identification problems. For example, parameter m u , representing the tumor maintenance cost, is generally estimated with low precision [15,19] and, consequently, high CVs were considered in both the uncertainty scenarios. The GSA performed on the C T highlighted that, even if m u appears the formula of C T , its contribution to the output uncertainty is always close to 0 (panels A and B of Fig. 5). Thus, it was assessed that the prediction of the C T is quite robust with respect to an inaccurate identification of m u .
This work remarks the usefulness of the GSA which provides a diagnostic tool to evaluate the robustness of the inference process based on pharmacometric models [9]. Here, we proposed and applied a two-stages GSA based on the d sensitivity indices to deal with the estimation uncertainty and its associated correlations. In particular, in this work, the statistical dependencies between model input parameters were fixed according to typical correlation matrices obtained during model identification. However, an additional layer of complexity could be added by considering the uncertainty on the correlations, which becomes random variables with a certain distribution. This would imply that the probability density function of model inputs is uncertain and, consequently, the GSA results (i.e., sensitivity indices computed for each input parameter) too. In such case, a nested two-stages GSA could be performed to quantify the impact of each correlation on every sensitivity index [55]. Unfortunately, this strategy is computationally demanding and difficult to implement with lots of uncertain parameters in the input distributions and in presence of complex model simulations [55].
Although the developed GSA method was applied to analyse the effects of estimation uncertainty and its associated correlations, it can be used to characterize different sources of variations. As an example, in the Supplementary Material S5 it was exploited to assess the impact of the inter-tumor cell line variability of the DEB-TGI parameters on theC T . Further, the investigation of the impact of the inter-drug variability on a single tumor cell line could be another interesting analysis to perform [56]. Note that, the same approach can be applied to PBPK models in which parameters are numerous, strongly correlated and generally affected by a significant uncertainty/variability. Therefore, the proposed GSA method is recommended in presence of strong correlations between the model parameters as long as a sampling strategy, both neglecting and considering statistical dependencies, can be implemented, for example by using the analytical definition of the joint pdf, by sampling from the Bayesian posterior distribution [57] or by using copulas [58].

Conclusions
GSA provides a valid diagnostic tool for supporting modelling activities especially in the oncology field where mathematical models are often used to make inference from preclinical to clinical setting. However, the development of GSA methods well-defined and easy-tointerpret also in presence of statistical dependencies between model parameters is still at its early stage. The aim of this work was to introduce a novel Two-Stages approach based on the d sensitivity index to perform GSA in presence of correlated model inputs. The proposed workflow was used to characterize the impact of the estimation uncertainty on output metrics of a tumor-in-host DEB-TGI model. This approach helped to rank the parameters according to their impact on the output distribution, discerning whether a parameter mainly exerts a direct or 'indirect' effect. In this way, it was possible to identify uncertainties that should be necessarily reduced to obtain robust predictions for the outputs of interest.
Author contributions ADC planned and implemented the methodology and wrote the manuscript; EMT planned the methodology, reviewed results and the manuscript; NM planned the methodology and reviewed the manuscript; PM supervised the work and reviewed the manuscript.
Funding Open access funding provided by Università degli Studi di Pavia within the CRUI-CARE Agreement.

Declarations
Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.