Efficient sensitivity analysis for biomechanical models with correlated inputs

In most variance‐based sensitivity analysis (SA) approaches applied to biomechanical models, statistical independence of the model input is assumed. However, often the model inputs are correlated. This might alter the interpretation of the SA results, which may severely impact the guidance provided during model development and personalization. Potential reasons for the infrequent usage of SA techniques that account for input correlation are the associated high computational costs, especially for models with many parameters, and the fact that the input correlation structure is often unknown. The aim of this study was to propose an efficient correlated global sensitivity analysis method by applying a surrogate model‐based approach. Furthermore, this article demonstrates how correlated SA should be interpreted and how the applied method can guide the modeler during model development and personalization, even when the correlation structure is not entirely known beforehand. The proposed methodology was applied to a typical example of a pulse wave propagation model and resulted in accurate SA results that could be obtained at a theoretically 27,000× lower computational cost compared to the correlated SA approach without employing a surrogate model. Furthermore, our results demonstrate that input correlations can significantly affect SA results, which emphasizes the need to thoroughly investigate the effect of input correlations during model development. We conclude that our proposed surrogate‐based SA approach allows modelers to efficiently perform correlated SA to complex biomechanical models and allows modelers to focus on input prioritization, input fixing and model reduction, or assessing the dependency structure between parameters.


| INTRODUCTION
7][8] This resulted in, among others, guidelines documents for assessing the adequacy of concluded verification and validation activities, which demonstrates the credibility of computational models of medical devices, by the American Society of Mechanical Engineers. 9][11][12] Sensitivity analysis is an integrated part of this uncertainty quantification exercise, because it identifies the relative contribution of uncertain inputs to the model output, either by the input itself or by its interactions with other model inputs.This information is relevant for model development and personalization.First, for model simplification, by removing redundant modules and/or fixing parameters that have no significant effect on the output. 11,13Additionally, this information is relevant for guiding the measurement protocol by measuring the most significant input values (input prioritization). 14,15ariance-based sensitivity analysis is currently considered best practice. 16The variance-based method is based on quantifying the model output uncertainty by its output variance and subsequently attributing fractions of this variance to the individual model inputs and their uncertainties, or their contribution due to possible interactions with other uncertain inputs.][19][20] Most variance-based sensitivity analyses assume that all model inputs are statistically independent, that is, the value of one input is not associated with the value of another one.However, in many cases model inputs are statistically dependent and this may affect the importance ranking of the sensitivity analysis dramatically. 21Consequently, the modeler might make misguided decisions in different phases of the modeling process if statistical independence is assumed.
Previous research has proposed uncertainty quantification and sensitivity analysis techniques that take correlations into consideration. 22,23However, these methods have only been able to provide an overall importance measure for one input variable.Brault et al. 24 did consider spatial aortic stiffness correlation among neighboring sections of the aorta, but, they still assumed independency between the parameters describing the spatial correlations and other parameters, such as the heart rate.
The sparsity of sensitivity analysis studies in literature that consider correlations between inputs may be attributed to various factors.These include (1) the considerable computational expenses associated with obtaining accurate sensitivity analysis metrics, (2) the challenges involved in finding the correlation structure among parameters, particularly in the clinical context where measurements are scarce, and until recently, and (3) the absence of universally applicable techniques for decomposing variance contributions into their structural and correlative components for variance-based sensitivity indices.
A solution for the absence of universally applicable techniques has been proposed by several studies, [25][26][27] which extended the variance-based sensitivity analysis for models with correlated inputs.Out of all these studies, Li et al. 26 proposed a sensitivity analysis method that generalized the variance-based method for correlated inputs.They demonstrated that the same mathematical definition as the well-known main and total Sobol sensitivity indices can be used. 28owever, these indices should be interpreted in a different way.In addition, the new method allows for apportioning the independent and correlated parts of the variance contributions of the inputs when considering statistical dependencies.Unfortunately, this method is still computationally expensive for biomechanical models, especially for complex models with large sets of inputs since high numbers of model evaluations are required. 26he aim of this study is to adapt the method of Li et al. in order to increase its computational efficiency and improve its practical applicability within the field of (patient-specific) biomedical modeling.We do this by introducing a twostage approach in which first an efficient surrogate model of the original model is created.Subsequently, the surrogate model is used instead of the original model during sensitivity analysis.Additionally, we demonstrate how this new method can guide the modeler during model development and personalization, even when the correlation structure is not completely known beforehand.
The surrogate models are used to optimize the sensitivity analysis methodology with regard to computational costs.A surrogate model can quickly evaluate the input-output relation of a biomechanical model.The surrogate models are obtained through the application of kernel-based methods, that allow for the creation of non-linear surrogate models based on data samples that make fast and accurate predictions when evaluated with new inputs.An advantage is that kernel methods do not need excessively large amounts of data. 29,30ubsequently, the correlated sensitivity analysis method of Li et al. 26 is applied to the surrogate model.We will refer to our methodology as a surrogate-based sensitivity analysis (SSA).After verification of the different components of our methodology, we applied it to a one-dimensional pulse wave propagation model to demonstrate its efficacy, and utility to guide personalization of the model, even when the correlation structure is not known.2][33][34][35][36] In addition, the pulse wave propagation model has a large input space ( > 20).

| Surrogate modeling approach
The surrogate modeling approach used in this article is a kernel-based method, called the Vectorial Kernel Orthogonal Greedy Algorithm (VKOGA), which was introduced by Santin and Haasdonk and allows for the generation of a surrogate model from a finite set of input and output variables. 30The outputs can be either scalar-or vector-valued outputs.
Here, we will briefly explain the VKOGA method.For a detailed mathematical description and the algorithm, we refer to References 30,37.

| Kernel methods for surrogate modeling
A surrogate model S : ℝ d !ℝ q can be constructed using a set X n ≔ x j j x j ℝ d È É n j¼1 of inputs and a corresponding set Y n ≔ y j j y j ℝ q n o n j¼1 of outputs of the full order model f : ℝ d !ℝ q .We will start describing the kernel method by only considering the case of scalar-valued functions, that is, q ¼ 1, and at the end of this section, we present the generalization for vector-valued functions.The construction of the surrogate model S : ℝ d !ℝ q makes use of a (strictly) positive definite kernel K : ℝ d Â ℝ d !ℝ. 30 We can then construct (for q ¼ 1) a surrogate model S : ℝ d !ℝ by using the following definition: with α j the jth unknown scalar coefficient.The surrogate model output calculates the output of the full order model f x ð Þ for all x X n , hence: with x i as the ith set of input parameters.This set of equations can be represented as a linear system with in which A K,X n ℝ nÂn is the kernel matrix, for which it holds that A i,j ≔ K x i , x j À Á .If the kernel matrix A K,X n is positive definite, then the linear system of Equation (3) has a unique solution.This is the case when the kernel function K is strictly positive definite.
In many applications, data values are often affected by noise.Therefore, the strict interpolation requirement given in Equation ( 2) is not desired and the solution can be found via the optimization problem 30 : This also helps to avoid possible oscillations within the surrogate model output. 29As can be seen in Equation ( 5), the regularization term λ k Sk 2 H is added.The regularization parameter λ ≥ 0 is used to control the tradeoff between the data-accuracy term and the regularization term.This is done to find the candidate surrogate model (S) with the lowest complexity, that is, penalizing the candidate surrogate models with a large norm in the reproducing Kernel Hilbert space.The representer theorem 38 guarantees that the optimization problem in Equation ( 5) has a unique solution of the form as defined in Equation ( 1), when a positive definite kernel is used.The unknown coefficient α ℝ n is then defined by the solution in which I is the identity matrix and α and y are column vectors.When considering vector-valued functions (or multioutput functions), i:e: q > 1, the same theory described above can be applied to each of the q components by using separable matrix-valued kernels. 37By doing so, the unknown scalar coefficient α j in Equation (1) now becomes an unknown column vector α j ℝ q .The same holds for the output f x i ð Þ.This results in the following matrix notations for solving the linear system of equations given in Equation (6):

| Greedy approximation
In general, the computational cost for the evaluation of Equation ( 1) is related to the number of model input-output pairs, that is, the n elements in the sum.To minimize the number of elements in the sum, only the terms that add significant non-zero coefficients α j to the expansion need to be included.Such a sparse expansion can effectively be found by using greedy algorithms.The greedy algorithm used here is the P-greedy algorithm and is summarized below, while for further information see References 39-41.
The P-greedy algorithm starts with an empty subset X n sub ¼ ; f g.Thereafter, the set is iteratively updated by adding a suitable new input X n sub [ x j È É with x j È É X n X n sub À1 .The P-greedy algorithm selects the input x j that is added in each step by picking the maximum on X n ∖ X n sub À1 of the nth power function: By minimizing the power function, the P-greedy algorithm provides an upper bound for the interpolation error.An exact definition of the power function can be found in Reference 30.Once the new input x j is selected, the quality of the surrogate S is checked by computing the residual with the remaining inputs x X n nX n sub À1 and outputs y Y n nY n sub À1 .More inputs are added until a given tolerance τ has been met, that is, A Newton Basis was used to improve the efficiency of the algorithm since it allows for a nested basis, where the new elements can be incrementally added without the need for recomputing the unknown coefficients every time a new input x j is added.For a detailed explanation of this efficient computational method, we refer to Reference 30.

| Sensitivity analysis
3][44] Moreover, we will give the connotations of the associated sensitivity indices in the case of both independent and correlated input parameters.For complete derivations and validation, we refer to respectively Saltelli et al. 42,44 and Li et al. 26

| Variance-based sensitivity analysis: independent parameters
To derive the global-variance-based sensitivity indices introduced by Sobol, 28 each model was considered as a black box with output Y and k stochastic input parameters X 1 , X 2 , …,X k ð Þthat contain uncertainty and are defined on the k-dimensional input domain Ω.Note that Y can contain either scalar or vector values but for notational convenience, we will further elaborate using the scalar notation.Furthermore, please keep in mind that this Y will turn into the surrogate approximation S x ð Þ, when the surrogate model is used to perform sensitivity analysis.After introducing the black box notation Sobol further expressed the model in terms of functions with increasing dimensionality (High Dimensional Model Representation, HDMR): Þ .After introducing the property that the integrals of all other summands over any of its parameters are zero, it was shown that all summands are orthogonal and that the decomposition is unique. 28The first three summands are for example given by: and Herein Ω Ài and Ω Àij indicate the domain Ω without parameter X i , and the domain without X i and X j , respectively.The vector X Ài represents all parameters but X i .The function p X 1 , …, X k ð Þis the joint probability of all k parameters, whereas p X Ài X Ài ð Þ and p X Àij X Àij À Á are respectively the marginal probability density functions obtained after integration of the joint probability density function over X i , and over X i and X j .For higher order terms similar equations as Equations ( 12) and ( 13) can be derived.
When the input parameters are orthogonal (independent) to each other, Sobol demonstrated by using the HDMR above that the total output variance V Y ð Þ can be decomposed to: in which the partial variances are defined in terms of variances of conditional expected values: Note that in the above derivation, it was assumed that marginal distributions such as After normalization with V Y ð Þ the Sobol decomposition has the following form: The terms S i are called the main Sobol sensitivity indices and can be interpreted as the expected reduction in output variance that can be obtained when the parameter X i would have been known exactly. 42The main Sobol indices are often used to decide (under uncertainty) which model parameters are most rewarding to measure more accurately, that is, parameter prioritization. 11The higher-order terms quantify the contribution to the total variances due to interactions between parameters.Later on, Homma and Saltelli 45 introduced the so-called total sensitivity index to quantify the total contribution of an individual parameter X i .These total sensitivity index is defined as: and includes all terms in which . The total sensitivity index can be interpreted as the expected variance that will be left if all other parameters would have been set to their true values.This index is commonly used to determine which parameters can be fixed on constant values within their uncertainty domain and which parameters are irrelevant to measure, that is, parameter fixing.

| Variance-based sensitivity analysis: correlated parameters
In the case of correlated parameters, the Sobol variance decomposition in Equation ( 14) no longer holds. 26This is due to the fact that the relationship between each part of the variance contribution is now more complex. 46However, the definitions of the main and total sensitivity indices can still be used to quantify the variance contributions of the different parameters but their connotation change.The change in the definition of the connotations of the different sensitivity indices is depicted in Figure 1 and will be described in this section.Li et al. 26 provided a general validation for the accurate meanings of these sensitivity indices.
As can be seen in Figure 1, the main sensitivity index in the case of independent model parameters becomes the total correlated contribution in the case of correlated model parameters, that is, Herein is the conditional expected value of Y given X i expressed by: and its corresponding variance by: The index S TC i represents the fractional contribution to the output variance in which correlations (parameter dependencies) are involved.The total correlated contribution can be subdivided into: in which S C i involves the contributions by interactions associated with correlations of X i and X Ài , and contributions by X i correlated with X Ài , whereas S U i captures the independent contribution of parameter X i .The sensitivity index S U i can be determined by: Given S U i and S TC i , the index S C i can be found by using Equation (24).
As can be seen in Figure 1, the expression of the total sensitivity index in the case of independent parameters should be considered as the total uncorrelated contribution in the case of correlated parameters.The total uncorrelated contribution is therefore given by: and can be further subdivided into two different indices: in which S IU i is the independent fractional contribution to the output variance of interactions between X i and X Ài .This latter index can be calculated as soon as S TU i and S U i are known.A visual representation of how the sensitivity indices are related to one another is given in Figure 2 in the case of two input parameters X i and X j .
The figure shows how all components are related to one another.The variance contribution associated with their correlation overlaps and from this, we can derive the relationships formulated in Equations ( 24) and (28). 46A full proof of all the sensitivity indices for correlated inputs can be found in Li et al. 26

| Calculation of the sensitivity indices
To estimate all sensitivity indices for both cases, that is, with and without correlations between parameters, we have implemented the approach introduced by Li et al. 26 into MATLAB R2022a (MathWorks, Natick, MA, USA) and will only briefly explain the key details in this section.

Algorithm in-and outputs
For Gaussian input distributions, the algorithm has as inputs the mean vector and the covariance matrix, whereas the maximum and minimum values and a correlation structure via Gaussian copulas are the inputs when considering uniform distributions. 47The outputs of the algorithm are the variance V Y ð Þ, and the sensitivity indices S TC i , S TU i , S U i , S IU i and S C i .
F I G U R E 2 Components of the computed correlated sensitivity indices.With S U as the uncorrelated contributions of both X i and X j .They both share a correlated contribution S C ij .Note that in this simple case with two parameters: Their shared independent factorial contribution is given by S IU ij .Figure is based on Hao et al. 46 Calculation of the integrals The integrals that define the variance V Y ð Þ, the expected value E Y ð Þ, and all nested integrals in Equations ( 21), (25), and ( 27) are approximated via sparse grid interpolation (SGI) using the Smolyak algorithm. 48A step-by-step description of these estimations can be found in Li et al. 26 Gaussian-Hermite quadrature nodes and weights were used in the case of Gaussian distributions while the zeroes of Gauss-Legendre polynomials were used in the case of uniform distributions.The lower error bound of the Smolyak algorithm scales with the number of dimensions and the accuracy level of the Smolyak algorithm was set to 3, representing 5th order polynomial exactness.To generate the quadrature nodes and weights using the Smolyak algorithm, we have used the MATLAB function nwSpGr that had been made publicly available by Heiss and Winschel. 49ince the sparse grid interpolation relies on the assumption of model input independency, the quadrature points need to be transformed from independent space to correlated space.For Gaussian distributions, this mapping was done analytically by using an orthogonal matrix with the eigenvectors of the covariance matrix as its columns (see Li et al. 26 ), whereas a Gaussian copula was used for mapping quadrature nodes in the case of uniform distributions.

Verification
To verify the correct implementation of the approach of Li et al., 26 we have benchmarked our numerical implementation with two polynomial functions and three variants of the Ishigami functions, as provided by Li et al.Accurate results for all benchmark problems were obtained (Appendix A).

| 1D pulse wave propagation model
The two-step approach described in the previous sections is applied to a 1D pulse wave propagation model to demonstrate the efficacy and utility of the approach for biomechanical models.Such a model was chosen because it has many model parameters but still allows for relatively fast model evaluations.This makes it practically feasible to benchmark our surrogate-based sensitivity analysis approach with state-of-the-art sensitivity analysis methods, given that the evaluation time of the 1D model evaluations for the benchmark method is being kept low.The pulse wave propagation model is governed by the 1D mass and momentum balance and a constitutive law.It is able to simulate pressure and flow waveforms throughout the cardiovascular system.We used the numerical implementation of the 1D-wave equations by Kroon et al. 50The pulse wave propagation model was solved for the domain discretized into non-overlapping two-noded elements in series with a size of 5 mm.The advection term was solved explicitly and a second-order backward difference scheme with a time step of 2.5 ms was used for time discretization.For more details regarding the numerical scheme, we refer to Reference 50.The geometrical and mechanical properties, the topology, the computational domain, and boundary conditions of Boileau et al. 32 were used as inputs.

| Simulations and analysis
In this section, we will briefly explain how we generated the surrogate models for different outputs of interest of our pulse wave propagation model.Afterwards, the various simulations will be discussed to verify our surrogate models and our two-stage methodology.Finally, we investigate how correlations between inputs affect the sensitivity indices, and demonstrate how these sensitivity indices can guide the modeler.

| Surrogate model construction
The VKOGA method, described in Section 2.1, was used to make surrogate models of the 1D pulse wave propagation model.

Training and test set creation
As for the model input space X n ≔ x j È É n j¼1 & ℝ d , similar to Melis et al., 51 we chose to vary the Young's Modulus E, the length l and the radius r 0 of five groups of arteries, namely: aorta, neck, upper limbs, organs, and lower limbs.For the latter four groups, the peripheral resistance, R, and compliance, C, were varied as well.The resulting input space thus contains 23 parameters.The parameter ranges were based on literature and therefore deemed realistic. 20,51All input parameters and their ranges can be found in Table 1.The input space X n of the test and training set was generated by sampling with a Sobol low-discrepancy sequence, using the MATLAB function sobolset. 52In total, 5000 input sets were generated and used to run 5000 simulations with the 1D pulse wave propagation model.The outputs of interest, depicted in Table 2, were stored and constitute the output space Y n ≔ y k f g n k¼1 ℝ.Note that in the case of waveforms q > 1, and thus a time-series surrogate model is generated.The PWV was computed as the foot-to-foot carotidfemoral pulse wave velocity. 53The input and corresponding output space were divided into a training set X tr ,Y tr f g n tr¼1 and an independent test set X ite , Y ite f g n ite¼1 of equal size (n = 2500).The first set will partially (n = 2250) be used for training with k-fold cross validation and partially for preliminary model testing on a small subset (n = 250).The independent test set will be used for independent evaluation of the resulting surrogate models on a large test set (n = 2500).
We chose to add an extra independent test set for an extensive verification of the surrogate models their performance.An overview of all datasets and their sizes is summarized in Figure 3.

Training of the surrogate models
The surrogate models are trained on the training data set while using a Gaussian radial basis function (RBF) kernel, which reads: where, γ represents the shape parameter, defining the shape of the kernel, and … the ℓ 2 -norm.For the P-greedy selection rule algorithm, the maximum number of selected inputs was set to 2500 (the entire training set).The tolerance on the power function P K,X nsubÀ1 and the residual r n was set to 5 Â 10 À8 [À] and 1 Â 10 À4 [À], respectively.As previously mentioned, the hyperparameters used for the training of the surrogate model are the shape parameter γ > 0 of the RBF kernel and the regularization parameter λ ≥ 0. Selecting the optimal hyperparameters can be quite challenging and the optimal combination can differ extensively for different types of biomechanical models.Therefore, in order to find the optimal hyperparameter combination, the model was validated iteratively for different hyperparameter combinations during the training.The ranges for the shape parameter γ were 1 Â 10 À2 ,50 ½ and for the regularization parameter λ were 1 Â 10 À16 ,1 Â 10 À2 ½ .For both hyperparameters, 10 logarithmically spaced points were chosen between both ranges, resulting in a total set of 100 different hyperparameter combinations.This selection procedure of the optimal T A B L E 1 Input parameters and their corresponding uncertainty domains used for the training of the surrogate model.

Input parameters
with n the number of elements of X, S the surrogate model that is being evaluated, and … the ℓ 2 -norm.This is repeated iteratively for all 10 disjoint validation subsets.After k-fold cross-validation, the mean error is computed of all 10 surrogates, quantifying the validation error.This in turn is done for all 100 different combinations of hyperparameters and the combination resulting in the lowest validation error is deemed the best hyperparameter choice and is chosen for the training of the final surrogate model.The remaining 10% of the training set X te , Y te f g(n = 250) was then used to evaluate the test error of the resulting surrogate model since it is imperative to test the quality of the surrogate on an independent test set in order to prevent overfitting.The test error was again computed using the ε RMSE , but now defined as ε RMSE S, X te , Y te ð Þ .Eventually, the VKOGA method gives as output: the surrogate model S x ð Þ, the optimal hyperparameters λ and γ, and the test error ε RMSE S, X te ,Y te ð Þ .evaluating them with the large independent test set X ite , Y ite f g n ite¼1 of n = 2500 different simulations.The evaluation of the required training set size was done for the Pao, Pao s , Qcca, and the PWV .This includes both scalar and waveform outputs.Normalization was achieved by dividing the ε RMSE by the mean of all model outputs, including the outputs of the surrogate model and the outputs given in the independent test set.To assess the mean of the model outputs, the mean difference is computed in order to avoid outliers and the absolute mean is taken so negative values will not lead to a mean value around zero, which can cause spikes in the determination of the error.

| Surrogate model performance on a large independent test set
After the surrogate models were generated, as described in Section 2.4.1, using the entire training set (n = 2500), the performance of the surrogate models was assessed by comparing the results of the surrogates again to the generated large independent test set X ite , Y ite f g n ite¼1 .Bland-Altman plots were generated for the scalar outputs in Table 2 (Pao s , Pao d , and PWV ) using a 95% limits of agreement (μ AE 2 Á SD) and confidence interval.The test error was quantified for all outputs of interest by computing the ε RMSE S,X ite ,Y ite ð Þ(Equation 30).Additionally, the normalized root mean square error (ε NRMSE ) was obtained for better comparison.

| Verification of the SSA
To verify the two-stage methodology, the SSA approach was first performed assuming input independence.These statistically independent sensitivity analysis results were compared to the established adaptive generalized polynomial chaos expansion (agPCE) by Quicken et al., 54 which also assumed statistical independence of the input space.The sensitivity indices of the independent input parameters and their effect on two typical scalar features, namely the systolic and diastolic aortic pressure, were computed with both SSA and agPCE considering uniform distributions.Additionally, we have benchmarked our SSA methodology (Appendix B) for correlated cases on one variant of the strongly non-linear and non-monotonous Ishigami function (Equation A6) provided by Li et al. 26

| Effect of correlations on the total correlated and total uncorrelated sensitivity indices
As proof of principle, correlated sensitivity analysis was performed on a pulse wave propagation model to demonstrate the effect of correlations on the total correlated and total uncorrelated indices.To investigate these effects, the sensitivity indices of input parameters and their influence on the pulse wave velocity (PWV ) were computed.Uniform input distributions were assumed and a Gaussian copula was used to describe the dependency structure.The output of interest PWV was chosen as it depends on the radii of the arteries, which is likely to be correlated to, among others, the length of the arteries in reality. 55,56Sun et al. found the correlation coefficient for the aortic length and radius to be ρ = 0.66.However, since the exact correlation between the input parameters would not always be known in real-life applications, we study the full range of the correlation coefficient ρ À1, 1 ð Þ between l aorta and r 0aorta .As such our method is more generally applicable to common scenarios in biomechanical modeling.

| Interpretation of the indices and use for model guidance
To explain the results obtained in the previous subsection, we thoroughly analyze the independent (S U and S IU ) and correlated (S C ) parts of the total correlated and total uncorrelated sensitivity indices.Additionally, the correlated sensitivity analysis of Section 2.4.5 is used to demonstrate how correlated sensitivity analysis can guide the modeler during model development and personalization.In this demonstration, first, the sensitivity indices of the uncorrelated case (ρ = 0) are analyzed.Subsequently, the same analysis is performed for the correlated case.Finally, the importance of taking into account input correlations is investigated by comparing the interpretation of the correlated and uncorrelated sensitivity analysis results.

| Surrogate model construction
To show the variation within the training data set, the pressure and flow waveforms of the aorta and common carotid artery have been depicted.A large variation was observed in the pressure and flow waveforms of the aorta and the common carotid artery (Figure 4).For the aortic waveforms, the maximal pressure difference of the upper and lower boundary of the envelope of the waveform, taken over all the time steps, is 52.7 mmHg, and for the flow 51.8 mL/s.For the common carotid artery waveforms, the maximal pressure difference is 55.7 mmHg, and for the flow 40.4 mL/s.Moreover, for both pressure and flow, a subset of 10 waveforms is depicted in Figure 4. From these, the inter-variability of the waveforms can also be deduced.

| Evaluation of required training set size
Figure 5 shows the computed ε NRMSE for various training set sizes.It can be seen for the aortic pressure waveform, common carotid arterial flow waveform, and the aortic systolic pressure that the normalized root mean square error decreases when the training data set increases.Around a training set size of 500, the most significant surrogate model performance improvement has been achieved when trained on these three outputs of interest.For the pulse wave velocity, the most improvement is obtained when the size of the data set is increased up to only 75.Afterwards, only an improvement of 1:74 Â 10 À3 [À] was achieved, which is no longer significant compared to the other outputs of interest.Moreover, for the aortic pressure waveform, the aortic systolic pressure, and the pulse wave velocity, the error is already small ( < 2 Â 10 À2 ) when only 25 input-output combinations are used.Figure 7 shows three Bland-Altman plots of the scalar-valued outputs of interest obtained when using the sets of inputs of the independent test set on the surrogate models compared to their actual outputs generated with the pulse wave propagation model.As can be seen in Figure 7, the differences between the output of the surrogates and the actual results generated by the pulse wave propagation model are relatively small.For the systolic pressure, the mean difference is 8:5 Â 10 À3 mmHg, with À4:4 Â 10 À1 and 4:2 Â 10 À1 mmHg as the limits of agreement.For diastolic pressure, the mean difference is À1:5 Â 10 À2 mmHg, whereas the limits of agreement are À3:2 Â 10 À1 and 3:1 Â 10 À1 mmHg.A wider spread can be seen for the comparison regarding the PWV .The mean difference for PWV is 2:2 Â 10 À3 m/s and the limits of agreement are À1:3 Â 10 À1 and 1:3 Â 10 À1 m/s.

| Surrogate model performance on a large independent test set
The ε RMSE and ε NRMSE were computed for all output of interest and are depicted in Table 3.It can be seen in Table 3 that the surrogate models outputs resulted in small errors compared to the actual order of magnitude of the outputs of the pulse wave propagation model.The largest ε NRMSE is obtained when training the surrogate model with the common carotid arterial flow waveform as output.However, this error is only 3.1% of the average flow.All other outputs of interest resulted in lower errors.

| Verification of the SSA
The sensitivity indices of both the surrogate model and the polynomial function of the agPCE were estimated using both the systolic and diastolic pressure as an output of interest.The results are shown in Table 4.
It was found that both methods resulted in similar numerical values for sensitivity indices.The maximally observed difference between the two methods was only 6:7 Â 10 À4 [À].Moreover, SSA was applied to the Ishigami benchmark problem, assuming correlated input variables, and SSA produced accurate results with a maximal difference of 1:6 Â 10 À2 [À].See Appendix B for details.

| Effect of correlations on the total correlated and total uncorrelated sensitivity indices
Of all parameters, the total correlated sensitivity indices S TC of the correlated parameters, the aortic length (l aorta ) and aortic radius (r 0aorta ), are most impacted by the change in correlation coefficient ρ (see Figure 8A).An increase of 1:94 Â 10 À1 [À] can be seen when the correlation coefficient is changed from À1 to +1.It was observed that E aorta has the highest total correlated sensitivity index S TC for each value of ρ.The S TC of E aorta decreases only by 1:08 Â 10 À1 [À] when ρ increases from À1 to +1.The parameter with the second highest S TC is l neck .The total correlated sensitivity index of l neck decreases by 3:55 Â 10 À2 [À] when ρ is varied from À1 to +1.For this parameter, it is observed that the total correlated sensitivity index is larger than those of r 0aorta and l aorta until ρ ¼ 0:6 and ρ ¼ 0:8, respectively.In general, for every other non-correlated variable, the total correlated sensitivity index is hardly impacted (variation < 2%) by the correlations between r 0aorta and l aorta .
Figure 8B shows the total uncorrelated sensitivity indices (S TU ) for various coefficients ρ.The parameter E aorta has again the highest total uncorrelated sensitivity index for every value of ρ and now also the largest range in obtained sensitivity indices for various ρ's (1:08 Â 10 À1 [À]).The parameter l neck always has the second highest sensitivity index and decreases by 3:50 Â 10 À2 [À] when the correlation coefficient is changed from À1 to +1.The total uncorrelated sensitivity indices of the two correlated parameters, l aorta and r 0aorta , vary parabolically, where the highest total uncorrelated sensitivity index is found for ρ ¼ 0, the uncorrelated case.Their variation in total uncorrelated sensitivity indices is 5:36 Â 10 À2 [À] and 5:86 Â 10 À2 [À], respectively.Again, all other variables are hardly impacted by the correlation between r 0aorta and l aorta (variation < 2%).

| Interpretation of the indices
From the decomposition of the variance contributions (Figure 9), for both l aorta and r 0aorta , it can be seen that, when their correlation starts to become negative, the observed decrease in S TC can first be attributed to the decrease in the S U and the negative values observed for S C (see for ρ ¼ À0:5).However, when the correlation coefficient ρ goes towards À1, the S C becomes slightly positive for both parameters.Then the reduction in the S TC can mainly be attributed to the decrease in the S U .When the applied correlation between l aorta and r 0aorta becomes positive, an increase can be seen in the correlated part S C .The uncorrelated part S U decreases again, but the net result causes the total correlated contribution S TC to increase.Note that when ρ > 0, S TU is smaller than S TC .It was found for all uncorrelated input parameters, that the correlated part S C was zero.Therefore, the uncorrelated part S U only contributed to S TC for these parameters.Since the correlated part S C was zero for the uncorrelated input parameters, these indices were not further explored during the interpretation of the variance decomposition.Moreover, the largest S IU found for all 23 input parameters was 6:24 Â 10 À3 [À].This tells us that in our analysis since the contributions of interactions between all parameters are negligible, the correlated contributions are therefore purely based on the correlation between l aorta and r 0aorta .This explains the almost exact similarity of S C for both l aorta and r 0aorta .Additionally, since S IU is negligible for all input parameters, the uncorrelated part S U mainly contributed to the S TU .

| Example for model guidance
As anticipated, in the uncorrelated case of ρ ¼ 0 (Figure 8A), S C is zero.As such, the S TC equals S U and can therefore be interpreted the same as the more generally applied main Sobol index.Furthermore, the total uncorrelated index S TU is again similar to the total Sobol index.Under these circumstances, the equivalent main indices guide the modeler to prioritize from most to less important: E aorta , l neck , E lowerlimbs , E neck , r 0aorta , and l aorta .All other 17 input parameters have a main sensitivity index lower than 0:05 and are therefore not rewarding to estimate more accurately.For the equivalent total indices, the following importance ranking applies E aorta , l neck , E lowerlimbs , r 0aorta , E neck , and l aorta , which all have total sensitivity indices > 0:05 (Figure 8B).This indicates that the 17 other parameters can be fixed to any value within their ranges without considerably impacting model output variance.However, when potential correlations are accounted for, the interpretation of the total correlated S TC and total uncorrelated S TU can no longer be interpreted as the main and total Sobol indices, so we focus on the decomposition of the independent (S U and S IU ) and correlated (S C ) parts for model guidance.When there is correlation between l aorta and r 0aorta (Figure 9), the relative contributions of the independent (S U and S IU ) and correlated parts (S C ) change.
The independent variance contributions (S U and S IU ) of all uncorrelated parameters are hardly impacted.In descending order, E aorta , l neck , E lowerlimbs , and E neck their independent variance contributions (S U ) are still relatively high in F I G U R E 8 Total correlated (A) and total uncorrelated (B) sensitivity indices of all input parameters and their effect on the pulse wave velocity with varying correlation coefficients defined between r 0aorta and l aorta .The yellow box highlights the results when no correlations were present.comparison to the other indices.In contrast, when the correlation structure between l aorta and r 0aorta varies to À1 or +1, the independent variance contributions (S U ) of l aorta and r 0aorta become less.So if we, naively, would focus only on the independent contributions to the output variance (S U and S IU ), we might decide to fix l aorta and r 0aorta based on the resulting S TU .However, for a positive correlation structure, their S C increases, which indicates an increase in the contribution of the correlated part to the output variance.On the other hand, when the dependency structure becomes negative, the correlated sensitivity analysis tells the modeler that the contributions of S U , S C , and S IU of l aorta and r 0aorta are close to zero, indicating to ignore or even fix l aorta and r 0aorta .Accordingly, for all other uncorrelated parameters, when S U , S C , and S IU are close to zero, the correlated sensitivity analysis guides the modeler to fix these parameters as well.
Since the importance ranking of the most important parameter (E aorta ) based on both the total correlated and total uncorrelated sensitivity index does not change, the sensitivity analysis still guides the modeler to accurately measure this input parameter.However, the effect of correlation between l aorta and r 0aorta does drastically affect their sensitivity indices.Hence, once the modeler wants to address more than only the effect of E aorta , the results demonstrate that it would thus be imperative for the modeler to first properly assess the correlation structure between l aorta and r 0aorta , before obtaining more accurate measurements for any of these other five relevant (sensitivity index above 0.05) input parameters (l neck , E lowerlimbs , E neck , r 0aorta , and l aorta ).
F I G U R E 9 Decomposition of the variance of l aorta and r 0aorta and their effect on the pulse wave velocity with varying correlation coefficients.

| DISCUSSION
The aim of this study was to develop an efficient two-step sensitivity analysis methodology for correlated inputs and to demonstrate how this method can help modelers during model development and personalization of biomechanical models.This was done using kernel-based surrogate models to efficiently evaluate the input-output relation of a biomechanical model, in this case, a 1D pulse wave propagation model.Afterwards, this newly devised methodology, known as SSA, was used to study the effect of correlations between inputs on the sensitivity indices of the PWV output of the pulse wave propagation model and the accompanying obtained guidance.

| Performance of the surrogate model
It was shown that the VKOGA method introduced by Santin and Haasdonk 30 allowed for the creation of accurate surrogate models that could correctly mimic the input-output relation of the pulse wave propagation model.This was observed for both scalar and time-series outputs.In the Bland-Altman plots of the scalar outputs, the widest spread was found for the PWV .However, the mean difference is still very small (2:2 Â 10 À3 m/s).For the time-series outputs, the worst resembling pressure and flow curves still agree very well with each other.A maximal difference of only 1.64 mmHg and 7.27 mL=s was found, outperforming the clinical accuracy of blood pressure and flow measurement devices. 57,58Moreover, when using the VKOGA method, our results showed that an extensive training data set was not needed since a training set of about 500 samples already resulted in accurate surrogate models for most outputs of interest, which was in agreement with literature. 29,30Similar conclusions were drawn by Koeppl et al. who also demonstrated that the VKOGA surrogate model could serve as an excellent replacement for a pulse wave propagation model.In addition, we observed a decrease in error for pressure and flow that was similar to the errors reported by Koeppl et al. 29 In our study, we observed that the surrogate model had the most difficulties with accurately predicting flow, especially in the common carotid artery (NRMSE of 0.031).This is as well in accordance with the findings of Koeppl et al., 29 who also showed that it was more difficult to train the surrogate model with the flow as output.As such, we considered the surrogate model as an accurate representation of the pulse wave propagation model to be used for performing correlated sensitivity analysis using our two-step approach.
Increasing the number of training data samples has been proven to be effective in reducing the error and improving the accuracy of the surrogate model for most outputs of interest.It was noted, however, for the PWV there was not much improvement after using more than 75 input-output combinations.An explanation for the lack of improvement when training the PWV on more samples could be due to the fact that the resolution of the PWV is limited by the chosen time step size Δt used for the discretization of the model.We chose Δt ¼ 2:5 ms, thus this resolution accordingly serves as the lower boundary of the computable error.Most surrogates depicted in Figure 5 already performed well (average NRMSE of 1:44 Â 10 À2 [À]), even when the minimal data set size of 25 input-output samples was used for training, which makes this methodology relevant for cases where limited training data can be acquired/generated, which is in alignment with the statements of Koeppl et al. 29 and Santin et al. 30

| Verification of SSA approach
To verify the SSA approach, sensitivity indices of all model inputs on the systolic and diastolic pressure were computed and compared with those obtained with the agPCE method, 54 whilst assuming independent input parameters.The comparison showed almost perfect agreement between the two methods.Additionally, in Appendix B it was shown that SSA was also able to correctly compute the correlated sensitivity indices of the adapted version of the Ishigami function.Therefore, these results give a strong indication that there is no longer a need for computationally-demanding model evaluations and that SSA could serve as an effective and veracious manner for efficiently estimating the sensitivity indices of a biomechanical model.
The surrogate model of the 1D pulse wave propagation method could be evaluated within 0.05 ms, allowing for the sensitivity analysis to be performed within an hour.To train a surrogate model of the 1D pulse wave propagation model, we used a training set that contained a total of 2500 simulations.One evaluation of the pulse wave propagation model takes about 30 min.Therefore the generation of the training data set took about 1250 h.The training of the surrogate took on average 34 h to complete, resulting in a total of 1284 h of preparatory work.However, to perform the sensitivity analysis of the pulse wave propagation model more than 70 million model evaluations were necessary.The SSA approach proposed, therefore, reduced the time required to perform the correlated sensitivity analysis by approximately a factor of 27,000 [À], as compared to directly applying the method of Li et al. 26 on the full pulse wave propagation model.Please note, that this speed-up is specific for the pulse wave propagation model, different models might lead to different reductions in computational time.
Due to the broader applicability of our method, SSA would potentially be perfectly suited to also perform uncertainty quantification and sensitivity analysis on even more computationally demanding models, such as ones relying on 3D computational fluid dynamics.Overall, SSA has shown to be a promising method to perform quick and accurate sensitivity analyses on biomechanical models.

| SSA to guide model development and personalization with correlated inputs
Lastly, we have applied SSA to investigate the effects of correlations on the sensitivity indices and the resulting guidance for model personalization.The results showed that taking into account and varying the correlation between input parameters can have a significant effect on the computed sensitivity indices.The interpretation of these new sensitivity indices is not as straightforward as the original main and total Sobol indices.Information drawn from the total correlated and total uncorrelated sensitivity indices should be interpreted carefully by accurately considering their individual components (S U , S C , and S IU ).It can be noted that the changes in all these sensitivity indices found when adapting the correlation coefficient show similar behavior as the ones found in Xiao and Duan, 21 who applied correlated sensitivity analysis to two analytical models, an additive model, and a non-additive model.
From the correlated sensitivity analysis results of this article can be concluded that while the decomposition of the total correlated contribution of the input variables can inform the engineers whether the correlated variations among specific variables or the variable itself are important to the output variance, the decomposition of the total uncorrelated contribution can provide useful information about the structural relationship between model inputs and outputs (e.g., contributions of the parameters themselves and interactions between them).This is in line with the definition provided by Li et al. 59 Moreover, it was found that the contributions to the variance due to the correlated part S C were the same for l aorta and r 0aorta .This is well in accordance with the findings of Li et al., 26 who stated that when only two parameters are correlated, their correlated contributions would be equal only if they had a similar interaction structure.This was the case for both parameters since their S IU were both nearly zero.If more correlations or interactions would be present, these indices would not be the same, because S C is defined as the correlated variance contribution of parameter X i , which includes the contribution by X i correlated with other parameters and the contribution by interactions associated with correlations of X i and all other parameters. 26Additionally, the behavior of S C and S IU when correlations are regarded allows for the total uncorrelated sensitivity index to become smaller than the total correlated sensitivity index (see Figure 9 for ρ > 0), which is not possible for the main and total Sobol indices.This occurs when S IU < S C .
If we apply SSA, which can deal with possible correlations, we see that the sensitivity indices of correlated inputs can vary extensively (about 0.2 [À] for the S TC of r 0aorta and l aorta ).When the indices are significantly impacted by the possible correlation structures, the obtained guidance for model personalization can vary extensively.For the input parameters l aorta and r 0aorta of the 1D pulse wave propagation model, depending on the chosen value of ρ, the distinction in obtained guidance can vary from advising: more accurate measurements, or fixing the parameters.This can considerably affect the reliability of the model.Therefore, to reduce the variance contribution of l aorta and r 0aorta , the SSA method advocates to first properly assess their correlation structure.Assuming input independence, as has been done in many studies before, 20,51,54 might, therefore, be unrealistic and could lead to incorrect sensitivity indices.These findings tell us that it is imperative to find out which input parameters are correlated, and then estimate the correct correlation coefficient.In case it is not possible to discover which parameters are correlated based on literature or available data, varying the correlation coefficients (which can be efficiently done using SSA) could provide more insight into the relative contribution of possible correlations to the output variance (uncertainty), and help the modeler decide whether considering correlations is necessary for certain inputs.To properly assess for which inputs of the 1D pulse wave propagation model the correlation structure needs to be assessed, the correlation structure between all input parameters, where any correlation is expected, should be varied and analyzed.However, varying multiple correlations between inputs simultaneously should be done with care to ensure that the well-defined mathematical constraints of the correlation matrix are not compromised.If the constraints are not met, the sensitivity analysis cannot be executed.Therefore, as a practical approach, it is advised to vary correlation coefficients only per input pair, resulting in a very sparse correlation matrix per sensitivity analysis.Although such an analysis is incomplete, it at least provides some idea regarding the relative importance of correlations on the output uncertainty.The real correlation structure between the input parameter pairs needs to be identified for those input parameter pairs whose different correlation structures affect the obtained guidance of the sensitivity analysis.Once the real correlation structure is known for impacted parameters, a more complete and accurate correlated sensitivity analysis can be performed and further guide the modeler with respect to model personalization and development.Performing correlated sensitivity analysis thus gives the modeler a better insight into what is most relevant for the next modeling phase, that is, should they invest time in parameter fixing, parameter prioritization, and/or whether they should invest time in unraveling the dependency structure of their input parameters.

| Performance of the surrogate model
The surrogates of the pulse wave propagation model performed really well.However, when applying SSA to more complex and computationally expensive models with a more complicated input-output relationship, the training of the surrogate might still become too computationally demanding for this method to be feasible.Applying SSA to more complex models could provide us insight into the feasibility of our methodology for such models.It should also be noted that by using surrogates instead of the actual model, a small error is made with regard to the computed sensitivity indices.Applying SSA onto multi-compartment systems, where an individual surrogate is made of each compartment, could lead to further propagation of this error.In such a situation, caution is advised.

| Verification of SSA approach
The parameter ranges during this study were based on literature and therefore, deemed realistic. 20,51Compared to Melis et al. 51 we have opted for a smaller range for the peripheral compliances C (À30%, 30% as opposed to À50%, 50%).Nonetheless, as the sensitivity indices of C are negligible, we found similar sensitivity analysis results for ρ ¼ 0.Moreover, for the purpose of this article, exact assessment of the input parameter ranges was not needed, since the obtained sensitivity indices just serve as a means of discerning the effect of correlations.It should be noted, however, that the results of sensitivity analysis are problem-dependent and thus depend on the output of interest, the model under consideration, and the ranges of the chosen input space. 60Therefore, for new sensitivity analyses in practice, the modeler should properly characterize the uncertainty ranges of their input parameters. 61n this article, the Sobol indices were only computed for scalar cases, even though surrogates were also generated for the time-series cases.If one wants to perform SSA with the time-series data as output, which is often the case for biomechanical outputs, several approaches can be used.An often applied approach is the computation of the Sobol indices pointwise in time.However, this approach disregards the temporal correlation between the data points. 62This might be especially problematic when considering time data over a short interval.3][64] The excellent ability of the VKOGA method to create timedependent surrogate models of the pulse wave propagation model makes the inclusion of temporal correlations into the proposed correlated SSA method likely practically feasible.However, whether taking into account both temporal correlation and correlations between inputs can be done simultaneously should be demonstrated in future work.
Additionally, it was assumed that the variance is an adequate measure for the uncertainty of the pulse wave propagation model.However, other studies have indicated that other methods might be preferable if the variance is not a good measure.For instance, entropy-based methods and methods that consider the entire probability density function are favorable when there is a strong skewness within the output distribution or when the output is multi-modal. 65,66here was, however, no strong skewness found in our outputs of interest, so there is a strong indication that the variance is a good proxy for the uncertainty of our model.Nonetheless, one should consider this before applying SSA to other applications.

| SSA to guide model development and personalization with correlated inputs
It was noted when going towards a negative correlation, the correlated contributions S C of l aorta and r 0aorta became negative at first, but eventually positive again.Similar behavior was observed by Hao et al. 46 Further research could provide insight into the underlying mechanics of this change.
It is known that in the cardiovascular system, correlations between various parameters are present. 67,68However, studies that quantify these correlations are sparse, and in clinical settings often only single measurements are done to minimize patient burden.A significant limitation that threatens the application of performing correlated SSA, therefore, is finding enough data to assess the input correlation in a robust manner.This problem aggravates when the number of input parameters increases.However, in practice, as long as data availability might be lacking, we believe that applying our methodology would be advised for the input parameters that are expected to be correlated, that is, varying their correlation coefficient from À1 to +1.

| CONCLUSION
Our newly proposed methodology, SSA, allows researchers to efficiently compute the sensitivity indices of a biomechanical model with many input parameters whilst considering correlations.An improvement of a factor of 27,000 was achieved using SSA in comparison to applying sensitivity analysis in which the pulse wave propagation model is evaluated instead of the surrogate model.Additionally, this article demonstrates that taking into account input correlations can have a considerable impact on the sensitivity analysis of biomechanical models.Correlations could significantly alter the interpretation of the sensitivity analysis results and thus seriously affect the obtained guidance concerning input prioritization and fixing during model personalization.

A P P END I X A : ANALYTICAL FUNCTIONS
The same analytical models were used in this study as previously by Li et al. 26

A.1. | Non-linear function
The first analytical model is the non-linear function defined as The analytical variance contributions were computed as shown in Li et al. 26 The sensitivity indices of the non-linear analytical model were computed by applying SGI as defined in Section 2.
Figure A1 shows the variance contributions for each parameter of the non-linear model obtained by applying SGI and solving the corresponding analytical equations.As can be seen in Figure A1, the SGI estimates are analogous to the analytical results of the non-linear model.Furthermore, the results obtained with our model are identical to the results obtained by Li et al. 26 F I G U R E A 1 SGI-based (SGI) and analytical (ANA) estimates of the variance contributions for the non-linear model.

A.2. | Function containing a square term of an input
The second analytical model is the response model containing a square term of an input.The analytical model is of the form Again next to applying SGI of Section 2.2.3 to the analytical model, the analytical variance contributions were computed as shown in Li et al. 26 considering Gaussian input distributions and with a level of accuracy of 3 for the Smolyak algorithm.In accordance with Li et al., we chose Figure A2 shows the variance contributions for each parameter of the squared term model obtained by applying SGI and solving the corresponding analytical equations.As can be seen in Figure A2, the analytical and SGI estimates are identical to one another.Again the results are indistinguishable to the results obtained by Li et al. 26

A.3. | Ishigami function
The third analytical model is the smooth, non-linear and monotonous Ishigami function defined as with X i is uniformly distributed over Àπ, π ½ .In analogy to Li et al., 26 we have chosen a ¼ 7 and b ¼ 0:1, and assumed that X 1 and X 3 were correlated with ρ 13 ¼ 0:4.For a more accurate approach, the accuracy level of the Smolyak algorithm has been increased to 15.This example contained uniformly distributed input variables, hence a Gaussian copula was used to transform the original uncorrelated problem to the case of the correlated normal problem.Figure A3 shows the computed sensitivity indices by the SGI based method.The results are again identical to the results obtained by Li et al. 26 Li et al. rewrote the Ishigami response model into the following two forms F I G U R E A 2 SGI-based (SGI) and analytical (ANA) estimates of the variance contributions for the squared term model.
When applying the SGI method described in Section 2.2.3 to these two forms, again with a ¼ 7, b ¼ 0:1, and ρ 13 ¼ 0:4, the following estimates of the variance contributions were obtained as shown in Table A1.Again the results obtained are quite similar to those obtained by Li et al. 26 The small differences in the estimates can be acknowledged to numerical interpolation and integration errors.

A P P END I X B: EVALUATION OF THE TWO-STEP SSA METHODOLOGY
A quasi-random input set was generated of 2500 by three inputs samples between [Àπ, π] using a Sobol low-discrepancy sequence.Using the rewritten Ishigami response model of Equation (A6) in Appendix A with the a ¼ 7 and b ¼ 0:1, 2500 corresponding outputs were generated.According to the SSA methodology, a surrogate model of the rewritten  Ishigami model was created using the quasi-random input set.Afterwards, correlated sensitivity analysis was performed with ρ 13 ¼ 0:4 and the accuracy level of the Smolyak algorithm again set to 15.The obtained results are depicted in Table B1.Again the results are very similar to the ones obtained by Li et al. 26 and the ones previously computed and depicted in Table A1, thus verifying our implementation.
hyperparameters was done by training the models using k-fold cross-validation with equal subsets (n = 2025) of the training set X trÃ , Y trÃ f gwith size n = 2250.After each training, the resulting surrogate model was evaluated on the remaining validation set X val , Y val f g(n = 225).The candidate surrogate models their performances were quantified by computing the root mean square error (ε RMSE ):

2. 4 . 2 |
Evaluation of required training set size It was assessed if a smaller training set could have been used for training of the surrogate models by incrementally increasing the training set and computing the normalized root mean square error for the trained surrogate models by F I G U R E 3 An overview of the different datasets mentioned within this study used for both the training and evaluation of the surrogate models.

Figure 6
Figure6shows a comparison of the obtained waveforms from the surrogate model and the 1D pulse wave propagation model using the same test set sample.The waveforms depicted are the ones obtained from the test set that showed the worst resemblance with the corresponding surrogate output.As can be seen in Figure6, both in the aorta and in the common carotid artery the pressure and flow waveforms are still very similar using both methods.They almost overlap.The maximal difference is 0.53 and 1.65 mmHg for the Pao and Pcca respectively, and 7.27 and 1.20 mL/s for the Qao and Qcca respectively.

F I G U R E 5
The ε NRMSE of the surrogate model output when increasing the data set size evaluated on the independent test set (X ite ).F I G U R E 6 Comparison of the pressure (red) and flow (blue) waveforms generated by the surrogate models and pulse wave propagation model (PWPM).F I G U R E 7 Bland-Altman plots of the scalar outputs (Pao s , Pao d , and PWV ) generated by the surrogate models and pulse wave propagation model.
2.3 considering Gaussian input distributions and with an accuracy level of 3 for the Smolyak algorithm.In accordance with Li et al., we chose

F
I G U R E A 3 SGI based and analytical estimates of the variance contributions for the Ishigami model.T A B L E A 1 SGI-based estimates of the various variance contributions for the response models of Equations (A6) and (A7).
Outputs of interest obtained from the simulations with the 1D pulse wave propagation model.
d Aortic pressure waveform Pao Aortic flow waveform Qao Common carotid arterial pressure waveform Pcca Common carotid arterial flow waveform Qcca Pulse wave velocity PWV T A B L E 3 Computed ε RMSE and ε NRMSE for all outputs of interest.
T A B L E 4 Computed main sensitivity indices for all input variables using both the agPCE and SSA method.
T A B L E B 1 SSA-based estimates of the various variance contributions for the response models of Equation (A6).