Global Sensitivity Analysis of a Homogenized Constrained Mixture Model of Arterial Growth and Remodeling

Growth and remodeling in arterial tissue have attracted considerable attention over the last decade. Mathematical models have been proposed, and computational studies with these have helped to understand the role of the different model parameters. So far it remains, however, poorly understood how much of the model output variability can be attributed to the individual input parameters and their interactions. To clarify this, we propose herein a global sensitivity analysis, based on Sobol indices, for a homogenized constrained mixture model of aortic growth and remodeling. In two representative examples, we found that 54–80% of the long term output variability resulted from only three model parameters. In our study, the two most inﬂuential parameters were the one characterizing the ability of the tissue to increase collagen production under increased stress and the one characterizing the collagen half-life time. The third most inﬂuential parameter was the one characterizing the strain-stiffening of collagen under large deformation. Our results suggest that in future computational studies it may - at least in scenarios similar to the ones studied herein - sufﬁce to use population average values for the other parameters. Moreover, our results suggest that developing methods to measure the said three most inﬂuential parameters may be an important step towards reliable patient-speciﬁc predictions of the enlargement of abdominal aortic aneurysms in clinical practice. future computational studies and save resources in parameter studies.


1I n t r o d u c t i o n
Growth and remodeling of arteries has been researched extensively during the past 15 years. In particular their important role in diseases such as aneurysms, which belong to the most important causes of mortality and morbidity in industrialized countries, has attracted significant attention. Aneurysms are local pathological dilatations of blood vessels that often keep growing over years until the blood vessel ruptures. Understanding and predicting this process is important for planning surgical interventions and researching potential future therapies. In the early 2000s, Watton et al. [1] and Baek et al. [2] proposed the first computational models to understand the natural history and evolution of fusiform aneurysms. A few years later Kroon and Holzapfel studied for the first time growth and remodeling of saccular cerebral aneurysms and identified the "continuous turnover of collagen" as the "driving mechanism in aneurysmal growth" [3]. Together with the constrained mixture theory of growth and remodeling introduced by Humphrey and Rajagopal [4], [1][2][3] inspired a host of increasingly detailed studies of vascular growth and remodeling over the last decade, for example, . While computational constrained mixture models of arterial growth and remodeling have substantially contributed to our understanding of this complex phenomenon, their clinical application, for example, for the computer-aided planning of treatments and surgeries, is still pending. A major difficulty in transferring these models from academic studies to clinical practice is the determination of patient-specific mechanobiological model parameters, which are required to make individualized predictions. It is naturally difficult and in particular potentially very expensive -if possible at all -to determine all these parameters with high accuracy. Therefore, it is important to understand which of these parameters have the most impact on the results of computational predictions. Knowing this, research can focus on the development of novel approaches to measure at least these parameters in a way that is on the one hand compatible with standard clinical workflows and acceptably cheap and on the other hand still sufficiently accurate for meaningful computer-aided predictions. To understand, which parameters in computational models of growth and remodeling are most important, [30,31] proposed parametric studies where single parameters were varied in growth and remodeling model based on the constrained mixture theory introduced in [4]. While the studies presented in [30,31] provided important insights, they also had limitations. Most importantly, they were limited to variations of single parameters, skipping thereby completely the in general important interactions between different parameters [32]. One notable exception was the study of Valentín and Humphrey [33] who investigated the combined influence of two parameters. Another more recent branch of research uses Bayesian methods [34][35][36][37], which are well-known from other areas of applied mechanics to be powerful tools for quantifying the effect of parameter uncertainties. However, what remains missing is a mathematically rigorous global sensitivity analysis ranking the importance of all the different parameters of computational models of growth and remodeling.
In this paper, we are presenting such an analysis for the homogenized constrained mixture models introduced in [38]. Our analysis uses the mathematically rigorous variance-based approach of Sobol and Saltelli [39][40][41][42]. This method decomposes the variance in the model's output upon variation of the input parameters and determines the contribution of each input parameter to the output's variance. Thereby, it allows us to understand not only the importance of single parameters of homogenized constrained mixture models but, for the first time, also the importance of combinations of input parameters. The resulting global and also quantitative understanding of parameter sensitivity allows us to make comprehensive, reliable and even quantitative statements about which parameters are most important in order to make reliable computational predictions. At the moment, most parameters of computational models of growth and remodeling are difficult to determine in a patient-specific way that fits into a standard clinical workflow. To overcome this problem systematically over the next years, one will have to develop step by step more and more methods to this end. The key question thereby is, the development of which methods should be prioritized in order to increase the predictive ability of computational models as fast and as much as possible. The global sensitivity analysis presented in this paper provides an important basis to make rational and well-founded decisions with respect to this question because it provides mathematically rigorous and quantitative statements about the importance of the different model parameters. We thus hope that this paper can help to guide the future biomedical research focused on vascular growth and remodeling.
This paper is organized as follows. In Sect. 2, we first introduce the concept of global sensitivity analysis and Sobol' indices. In Sect. 3, we briefly summarize the homogenized constrained mixture model introduced in [38]. In Sect. 4, we discuss the details of a global sensitivity analysis of this model. The results of this analysis are summarized in Sect. 5 and discussed in Sect. 6. Finally, in Sect. 7, we discuss the conclusions drawn from these results with respect to potential future research into methods for the clinical determination of mechanobiological parameters.

Sobol Indices: An Approach for Global Sensitivity Analysis
Many parameter studies use so-called local methods where only single parameters are varied at a time. However, this approach is not suitable for nonlinear models because it neglects the possibly important interactions between different parameters and tends to underestimate the input space due to the "curse of dimensionality". To overcome these limitations, global sensitivity analysis methods try to infer the global influence of model parameters on the model output by quantifying the amount of uncertainty in the model output caused by the individual parameters including their interactions with other parameters [32].

Definition and Interpretation of Sobol Indices
This section introduces variance-based global sensitivity measures for general, nonlinear models, which are often referred to as Sobol indices [39][40][41][42][43][44]. To keep the notation simple, we abstain in the following from a notational distinction between random variables and their realizations and trust that the difference is evident from the context. Thus, let x = {x 1 ,x 2 ,...,x n }∈ denote a continuous random vector whose components x i are random variables. By x ∼i , we denote the random vector of all components except x i ,thatis, The expected value of a function y of x is denoted by where p x (x) is the probability density function of x. The corresponding variance follows as Let the model of interest be represented by an integrable function y with n mutually independent input parameters x ={x 1 ,x 2 ,...,x n } and scalar output such that For the sake of readability, we commit a slight abuse of notation by not distinguishing between the function y and its value y(x) in the following. Herein we assume, without loss of generality (cf. [42]o r [ 45,Chap. 15]) that the model parameters are distributed uniformly within the n-dimensional unit hypercube K n ,t h a ti s ,x i ∼ U(0, 1),w h e r eU(a, b) denotes the continuous uniform distribution on the interval [a,b] with −∞ <a<b<∞.D u et o the mutual independence of the x i , the joint probability density function is calculated as One can show that, under the above assumptions, there exists a unique decomposition of y, often called analysis of variance (ANOVA) representation [39,41] or high-dimensional model representation (HDMR) [45], such that where y 0 is constant and the components y ij k... = y ij k... (x i ,x j ,x k ,...) are functions of as many (up to n) arguments {x i ,x j ,x k ,...} as they exhibit subscripts. Using (5), we can decompose the total variance of y as where V x i x j x k ... y ij k... are the variances of the summands in (5). Dividing (6) where the define variance-based sensitivity measures, called Sobol indices. A Sobol index of order s (written with s subscripts) gives the fraction of the total variance of the model output that can be attributed to the interaction between the s different input parameters x i ,x j ,x k ,... alone. In general, the computation of all Sobol indices of a model is very expensive. For this reason, in practice one often computes only the first order sensitivity indices S i and the so-called total sensitivity indices S T i . The latter is defined as the sum of all Sobol indices where parameter x i is involved: One can prove [39][40][41][42] that the first order sensitivity indices and total sensitivity indices can be computed equivalently to the definitions (8)and (9)as where E x ∼i [y|x i ] and E x i [y|x ∼i ] are the expected values given the component x i respectively the vector x ∼i . The latter two equations are often used for the efficient computation of Sobol indices (see Appendix).
The first order sensitivity index S i is often also called the main effect of parameter x i . It describes the fraction of the variance of y that can directly be linked to an uncertainty in x i alone. In other words, it describes by which fraction the variance of y would reduce if the component x i were known exactly. The total sensitivity index S T i is also called the total effect of parameter x i . It describes the expected fraction of the variance of the output y that would remain if all parameters except for x i were known exactly. It is a measure of the combined influence of x i alone (i.e., its first order effect) together with all higher order interaction terms where x i is involved.
While the S i quantify the influence of each parameter alone on the model output variance, the S T i additionally quantify the influence of the interactions into which each parameter is involved. In practice, computing these two types of indices and omitting the other higher order indices defined above has been found to be a good trade-off between computational cost and insight into the characteristic properties of the model of interest [43]. For example, the S i and S T i can be used for the following analyses.
Linearity analysis: if n i=1 S i is close to one, the model is largely linear, whereas if this sum is close to zero, the model is dominated by nonlinear interaction terms. Similarly, on the individual parameter level, the difference S T i − S i gives the amount of variance of y due to all interactions where parameter x i is involved.
Parameter priority analysis: let us assume, we seek to reduce the uncertainty of our model as much as we can by measuring one of the input parameters exactly. Then the above delineated theory tells us that we have to focus on the parameter with the highest first order sensitivity index S i .
Parameter fixation: to simplify the execution of computations, it is often helpful to choose reasonable ad-hoc values for a couple of parameters without spending too much effort on measurements or analyses of their exact value. Such simplifications are acceptable, however, only if the effect of where exactly a parameter is fixed within a certain reasonable range can be trusted to be small. As pointed out by [43] this is the case for parameters with S T i ≈ 0. For practical purposes it is important to choose a smart strategy to compute the sensitivity indices because brute-force approaches quickly become prohibitively expensive as the number of model input parameters increases. The Appendix briefly summarizes such a smart strategy that enables the evaluation of a full set of first order and total sensitivity indices at the cost of just N tot = N(n+ 2) model evaluations, where N is the number of Monte Carlo sample points used for integration and n the number of model input parameters [42,46].

Continuum Mechanical Framework
In this section, we briefly summarize the concept of homogenized constrained mixture models of vascular growth and remodeling as developed in [38,47,48]. Thereby, we focus on the simple (but for this paper sufficient) special case that the blood vessel is modeled as a thin membrane. The deformation of this membrane is modeled on the basis of the theory of nonlinear continuum mechanics [49]. Thereby we model the artery as a continuum with some reference configuration 0 . Mechanical loading as well as growth and remodeling can result in a deformation of the artery over time t into some current configuration t .T h i s deformation translates each material point X in the reference configuration to at time t to a current position x(X,t)= X + u(X,t), where u is the so-called displacement field (Fig. 1). A key quantity to describe this deformation within the theory of nonlinear continuum mechanics is the deformation gradient Constrained mixture models assume that a mechanical body consists in general of m different constituents, distinguished in the following by m different superscripts i ∈ I gathered in an index set I . These different constituents share each differential volume element. They form a compound and thus deform together. However, the single constituents may exhibit different stress-free configurations [4]. In the theory of nonlinear continuum mechanics this concept can be modeled as follows. As the constituents deform together, they all exhibit the same deformation gradient F . However, for the different constituents it is in general multiplicatively decomposed into different elastic parts F i e and inelastic parts F i gr : The inelastic part changes only by growth and remodeling of the respective constituent. For each constituent, F and F i gr define its respective state of elastic deformation and thus its strain energy i (F i e ). The total strain energy of the constrained mixture is assumed to be the sum of the contributions of the individual constituents: Here i denotes a strain energy per unit (referential) volume of the i-th constituent, W i its strain energy per unit mass, and ρ i 0 its (referential) mass density. The mechanical stress in the continuum can be defined by the 1st Piola Kirchhoff stress where the partial derivative with respect to F should be understood in a way that all F i gr are kept constant so that only the F i e and thus also F may vary. Growth and remodeling in vascular tissue occurs on very long time scales so that inertia can be neglected. As typically also body forces such as gravitation are negligible in vascular tissue compared to the mechanical loading from blood pressure, the balance of linear momentum reduces to div P (F ) = 0. (17) Solving this equation renders at each point in time and space the a priori unknown deformation. This is possible with standard methods such as a finite element discretization of the blood vessel geometry as long as the inelastic parts F i gr of the deformation gradient are known. These can be computed by the mathematical model of growth and remodeling as pointed out below in Sect. 3.3.

Constitutive Equations
Following the lines of [50] and subsequently in the context of growth and remodeling, for example [30,31,38], we assume herein that vascular tissue can be modeled as a constrained mixture of elastin and a number of collagen and smooth muscle fiber families.
Forelastin,weuseasuperscripti = el and assume a Neo-Hookean strain energy per unit mass with the trace operator tr(·) and a material parameter μ.
Smooth muscle and collagen are modeled as uniaxial fiber families aligned in the reference direction with some unit vector a i 0 . Unlike in [30,31,38], for simplicity we do not distinguish herein between collagen and smooth muscle and rather assume that the fiber families in our model describe a mixture of both together. Their elasticity is governed by a Fung exponential function, which for the i-th constituent takes on the form Here ⊗ denotes a dyadic tensor product, the colon a double contraction product, C i e = F i e T F i e is the elastic right Cauchy-Green tensor, and a i gr = F i gr a i 0 / F i gr a i 0 .Notethatherein we assume that the two material parameters k 1 and k 2 of the Fung exponential function are identical for all the fiber families in our model. Vascular tissue can often be modeled as an incompressible material. Specifically when blood vessels are modeled as membranes as we are doing it herein, this constitutive assumption can easily be implemented by requiring that the out-of-plane component of F always equals the inverse of the product of the two principle stretches in in-plane direction. Doing so, the related elasticity problem can be solved by a discretization with two-dimensional membrane finite elements whose thickness parameter is simply always adjusted such that it preserves the tissue volume under elastic deformation, see also [30,31,38].

Remark 1
Note that herein we neglect for simplicity active smooth muscle tension because we primarily aim at aortic aneurysms. Unlike cerebral vessels, the aorta is an elastic rather than a muscular artery so that smooth muscle tension can be assumed to play only a minor role there compared to passive elasticity.

Growth and Remodeling
Living tissues are subject to a continuous mass turnover during which extant collagen and smooth muscle tissue is degraded and new such tissue is deposited. In mathematical models of growth and remodeling, it is often assumed that deposition and degradation balance in case the tissue is in a so-called homeostatic state. The exact definition of the homeostatic state remains controversial. For simplicity most mathematical models define it as a preferred (homeostatic) stress or strain of the collagen and smooth muscle tissue. Herein, we follow this approach and assume that no growth and remodeling takes place if the absolute value of the Cauchy stress σ i in collagen and smooth muscle fibers equals a preferred value σ h . Deviations from this preferred stress value are assumed to have two consequences.
First, they result in a net mass production, which is assumed to be governed for each constituent subject to growth and remodeling by the evolution equatioṅ where k σ is a dimensionless gain parameter for mass deposition and T represents the average life time of fibers during the continuous turnover process of deposition and degradation in the tissue. It is worth mentioning that the turnover time T corresponds to a half-life time of ln (2)T . This simple evolution equation assumes that the production of collagen and smooth muscle is directly proportional to the amount of extant tissue, which appears reasonable, because production and degradation are mainly driven by cells, whose number can be assumed to scale in living tissues under typical conditions roughly linearly with the amount of tissue. Moreover, (20) assumes that the net mass production scales linearly with the deviation of the current stress from the homeostatic value, which can always be considered a straightforward first order approximation of reality, following directly from the rational of Taylor expansion around a homeostatic state with zero net mass production. In a nonlinear continuum mechanical membrane model of incompressible vascular tissue as assumed herein, net mass production and degradation can directly be included via an inelastic net change of the membrane thickness by a factor where the denominator is equal to the constant current mass density ρ = i∈I ρ i 0 (t = 0). For reasons discussed in [38], in membrane models of vascular growth and remodeling λ g (t) is typically assumed to govern the inelastic deformation due to changes of mass of all constituents alike.
The second result of a deviation between the current and the homeostatic stress in collagen and smooth muscle tissue is remodeling. Remodeling can be understood as a reorganization of the tissue microstructure, which results in an inelastic deformation of the tissue on the macroscale, in many respects similar (though not exactly identical [51]) to a viscoelastic deformation. As discussed in [38], for uniaxial fiber families as used herein to model collagen and smooth muscle this inelastic deformation can be captured by an inelastic part λ i r of the total fiber stretch λ i along the fiber axis, and the evolution equation for this inelastic fiber stretch can be shown to bė In (22) it is assumed that fibers are aligned in the in-plane direction of the membrane representing the vessel wall and that growth, that is, the above λ g results in a nonlinear deformation only in the referential out-of-plane direction a ⊥ 0 so that the in-plane fiber stretch λ i = λ i e λ i r can always be multiplicatively decomposed into an inelastic remodeling part λ i r and an elastic part λ i e . With these assumptions (20), (21)and (22) together define the temporal evolution of the inelastic part of the deformation gradient, which itself can be computed at each point in time from The first term on the right-hand side captures the thickening in out-of-plane direction by net mass deposition and the potential simultaneous transverse contraction that may result in case of an inelastic fiber remodeling stretch due to the assumption of incompressibility. The second term captures the inelastic fiber remodeling stretch, and the third term the resulting in-plane transverse contraction. The algebraic and differential equations (12) through (23) form a closed system that defines at each point in time and space the deformation of vascular tissue due to growth and remodeling if the material parameters and the parameters characterizing the vascular geometry are known. Both types of parameters form the set of input parameters to our model of growth and remodeling. It is the objective of this paper to analyze the sensitivity of the output of our model to variations of these input parameters. Details of this sensitivity analysis are discussed in the subsequent section.

Global Sensitivity Analysis of Arterial Growth and Remodeling
In this section, we discuss how the global sensitivity analysis framework from Sect. 2 can be applied to the homogenized constrained mixture model of growth and remodeling from Sect. 3. Thereby we focus on a generic, idealized model of the abdominal aorta described in the following subsection.

Idealized Model of Abdominal Aorta
Geometry: we study an idealized abdominal aorta represented by a thin-walled cylinder of diameter d = 2 cm, length L = 18 cm and wall thickness H . Dirichlet boundary conditions are imposed at both ends of the cylinder mimicking the support of the aorta by surrounding tissue and branching vessels such as the renal arteries. Our model aorta is subject to an internal mean blood pressure p = 100 mmHg. The vessel wall is modeled as a constrained mixture of m = 5 constituents, which are an elastin matrix and four fiber families modeling the collagen and smooth muscle fibers. Hence, our index set to distinguish between the different constituents is I = {el, co 1 ,co 2 ,co 3 ,co 4 }, where index el refers to elastin and the other four indices to the four collagen and smooth muscle fiber families. The referential unit direction vector a i 0 of these fiber families can be uniquely defined by the angle α i between them and the circumferential direction of the cylinder. We assume that one fiber family is oriented in circumferential and axial direction, respectively. Moreover, we assume that there exist two diagonal fiber families forming an angle of ±45°with the circumferential direction. This setting is illustrated in Fig. 2. Note that by variations of the mass densities assigned to the different fiber families one can resemble at least in the sense of a good approximation the effect of a great variety of different fiber orientation distributions, which endows our study with a sufficient generality.
Constituent mass in constrained mixture: the mass of the different constituents in the constrained mixture at any point in time is defined by their referential mass density ρ i 0 (t).I n our discussion below, it is convenient to express it in a normalized form in terms of a mass fraction The mass fractions of the different constituents satisfy a partition of unity property. For simplicity, we assume herein that both the mass fractions of the circumferential and axial fiber families and the mass fractions of the two diagonal fiber families are equal, respectively. The partition of the total fiber content between these two sets of fiber families remains, however, a free parameter. Overall, this reduces the number of free input parameters in the model and thus the computational cost of a global sensitivity analysis without compromising too much the generality. We formalize the above assumptions by first introducing the initial mass fraction of elastin ϕ el t 0 = ϕ el (t = 0). In the following, the index t 0 shall always refer to (initial) quantities at t = 0. Next, the initial mass fraction of all four fiber families together is ϕ co t 0 = 1 − ϕ el t 0 . Finally, we introduce the initial fraction of total fiber content attributed to one diagonal fiber family β t 0 ∈[0, 0.5]. In the extreme case of β t 0 = 0.5a l lfi b e rm a s s is initially oriented in diagonal direction. For soft tissue one often assumes as a reasonable approximation a constant spatial mass density ρ = 1050 kg/m 3 of the tissue as a whole. Under this assumption, the initial referential mass density of elastin is ϕ el t 0 ρ, the one of the diagonal fiber families take on the value β t 0 (1−ϕ el t 0 )ρ and the ones of the circumferential and axial fiber families the value In other words, with respect to referential mass densities our model has two independent input parameters, which are ϕ el t 0 and β t 0 .
Initial configuration: we assume that our model aorta is in a homeostatic configuration at t<0, that is, no growth and remodeling of the fiber families takes place until t = 0 because the Cauchy stress of all fibers equals the homeostatic value. In our simulations, we establish such an initial configuration as follows. First, the homeostatic stress σ h of the fiber families is prescribed as an independent input parameter. Second, the balance of linear momentum in the circumferential direction is used to eliminate the model parameter H , which is uniquely determined by this balance equation and the other input parameters of our model.
Elastin degradation and prestretch: during adulthood, no deposition of load-bearing elastin takes place [47]. Rather it is degraded with half-life time of a few decades. Therefore, (20) does not apply to elastin, but the elastin referential mass density is rather assumed to be governed by an evolution equation of the type where D(t) describes a time-dependent damage parameter between zero and one that can be used to model damage and loss of elastin and that is specified in more detail below. As elastin is not subject to growth and remodeling according to (20), its elastic prestretch in the inital (homeostatic) configuration at time t = 0 is not defined by σ h so that we have to define it independently. While there is some evidence that axial and circumferential prestretches of elastin are different in general [18,52], their values are typically found to be very similar in healthy blood vessels (see also [30,53]). In this study, we thus assume both of them to be defined by some (in principle independent) model input parameter λ pre .
The above paragraphs define an initial configuration resembling a healthy blood vessel. We will use this configuration as a starting point for two case studies of arterial growth and remodeling.
Case 1: Hypertension. In this example, we study the growth and remodeling response of our idealized aorta to hypertension, that is, a persistent increase of mean blood pressure p inside the vessel. For this purpose, we track the radial expansion of the model aorta over a period of 15 years. As already in previous work [12,13], we assume hypertension to arise in the following way: during the first year, the mean blood pressure increases linearly from the reference value p to an elevated level ofp = 120 mmHg. The pressure remains at this elevated level for the rest of the simulated time. For simplicity, we do not consider any damage to the elastin matrix in the context of hypertension, i.e., D(t) = 0.
Case 2: Idealized fusiform abdominal aortic aneurysm (AAA). Spontaneous damage to the elastin matrix has previously been hypothesized to trigger the development of aneurysms [13,26,30,38,53]. We adhere to this approach and apply the following damage to the elastin matrix of the model aorta: We use a coordinate system whose center coincides with the center of the cylinder representing our aorta and whose X 3 -axis is aligned with its rotational symmetry axis. The term on the right-hand side describes an abrupt, spatially distributed damage to the elastin matrix with spread parameter L dam . H(t) denotes the Heaviside function. Typically AAAs are accompanied with significant reduction in the elastin content of the tissue over time often until complete depletion [54][55][56][57]. Here, we model this degradation process with a simplified approach by instantly removing the complete elastin content at the center. The geometrical shape of the elastin damage has a large influence on the evolving aneurysm shape [17,19,25]. However, we limit our analysis to L dam and fix the shape as defined by (26)in order to keep the number of overall parameters for the global sensitivity analysis feasible. Again, we track in our simulations the radial expansion over a period of 15 years. The mean blood pressure inside the vessel remains constant at the initial value p.

Output
Naturally, the growth and remodeling response of the two cases specified in the previous section depends on the choice of model parameters. It is our main goal to quantify the sensitivity of the model output to the model input parameters for both cases. Sobol's method for global sensitivity analysis as introduced in Sect. 2 is defined for models with scalar outputs only. By contrast, the solution of the homogenized constrained mixture model of Sect. 3 results in a vectorial displacement field. To resolve this mismatch, we define herein as our model output the maximum current diameter d max of our aorta over time. We define Here the limit value of d max = 8 cm is assumed to represent the diameter where the aneurysm ruptures. This choice is motivated by statistical data showing that unruptured AAAs with a diameter larger than 8 cm are extremely rare [58]. Mapping all diameters greater or equal to 8 cm to the same diameter value of 8 cm is thus reasonable because all aneurysms with such a large diameter indeed represent the same model outcome, that is, a ruptured aneurysm. Of course, the rupture criterion applied here is very simplistic. However, it appears sufficient for the purpose of this paper which does not aim at quantifying the exact moment of rupture for a specific aneurysm but rather the impact of the different model parameters on growth and remodeling in arteries in general.

Known Input Parameters
In our sensitivity analysis, we do not include all the model parameters because some can be assumed to be known in clinical practice and others can be assumed a priori to play such a minor role that they need not be included in our study. Indeed, d and L can be measured in clinical practice by medical imaging, and p andp by standard blood pressure measurements.
On the other hand, the mass density ρ obviously has nearly no impact as we model arteries as thin membranes and define their strain energy per unit fiber mass. Finally, also our choice of α i does not significantly determine the model as long as the mass fractions of the different fiber families are allowed to vary freely. Therefore, we fix all these parameters according to Table 1 and do not include them in our sensitivity analysis.

Unknown Input Parameters
Fixing some parameters according to Table 1, there remain nine parameters in case 1 (hypertension) and 10 parameters in case 2 (idealized AAA). These are in both cases the homeostatic Cauchy fiber stress σ h , the turnover time T , the gain parameter k σ , the stiffness parameters μ, k 1 and k 2 , the initial mass fraction of elastin ϕ el t 0 , the initial fraction of the fiber mass attributed to the diagonal fiber families β t 0 and the prestretch of elastin λ pre (equal in axial and circumferential direction). In case 2 (idealized AAA), we additionally study the spatial damage spread L dam .

Remark 2
Due to the assumption of an initial homeostatic configuration (see Sect. 4.1)and the definition of stress (16), the choice of the collagen material parameters k 1 , k 2 and the homeostatic Cauchy fiber stress σ h as independent input parameters implies that the deposition stretch of collagen becomes a dependent parameter that is varied implicitly with these three parameters.
These parameters can typically not be measured in clinical practice at the moment and yet it cannot be excluded that they have a considerable impact on the model output. However, while patient-specific measurements of these parameters in vivo are not possible, we have information from various experimental and clinical studies, for example, from mechanical testing of tissue samples ex vivo. Therefore, we can at least define a reasonable range for the above parameters within which they can be assumed to vary. In the following, we will do so for all the above parameters and assume for simplicity (and typically lacking more specific information) a uniform distribution of the parameters within the defined bounds in our sensitivity analysis. The derived distributions are summarized in Table 2. gain parameter (case 2) k σ U(0.05, 0.150) - [ 12] Elasticity and mass fraction of elastin: The elastin content in healthy aortic tissue has been investigated in several studies [53,59,60]. In these studies, consistent values between 0.227 ± 0.057 [59]a n d0 .224 ± 0.031 [60] have been reported. Based on these findings, we set the bounds for the initial mass fraction of elastin ϕ el t 0 to [0.2, 0.3]. Elastin is predominantly deposited during early life [61] and has a very long mean life time of approximately 101 y, cf. Table 1, [62]. During normal biological growth, elastin therefore undergoes significant mechanical deformation which has been hypothesized to result in a considerable level of prestretch in the healthy aorta [63,64]. Mean values reported in the literature are between 1.18 and 1.37 [18,53,[65][66][67][68]. We follow these studies and assume that the prestretch of elastin λ pre varies in the range [1.2, 1.4]. Physiological ranges of mechanical wall properties of arterial tissue in health and disease have been studied extensively in the literature [57,[69][70][71][72][73][74][75][76]. However, the quantitative comparison of material parameters derived from different experimental studies is a challenging task. One reason for this is the fact that often the experimental data is fitted to different constitutive models. A general method to solve this problem is beyond the scope of this paper. As a simplistic solution, we tried to relate at least reported values for closely related constitutive models in a reasonable way to each other. For example, we made the parameter values from homogeneous models of the arterial wall comparable to the values reported for constrained mixture models by correcting them by a factor accounting for the mass fractions of the different constituents that can typically be assumed (see also [77]). Pooling in such a way the data reported in [18,53,68,73,75,77] -under the assumption of ρ as in Table 1 and ϕ el = 0.25 -we were able to define for the material parameter μ of elastin a range 40-80 J kg . For the spatial spread parameter L dam ,we choose 0.5-2 cm to mimic both localized as well as considerably spread elastin damages.
Elasticity and mass fraction of collagen: A quantitative comparison of the material parameters k 1 and k 2 for collagen from different sources is even more challenging than for the parameter μ of elastin. For modeling the constitutive behavior of collagen generally includes also several other structural parameters [78] whose choice may then also affect the values of k 1 and k 2 reported. These other parameters are, for example, the number of fiber families and their orientation or the fiber dispersion. Naturally, structural parameters and material parameters co-depend nonlinearly rendering approximate conversions, as suggested for the elastin case above, almost impossible [76]. Therefore, we had to limit our focus on a choice of papers using very similar constitutive models [30,53,68]. From these, we derive the parameter range 450-600 J kg for k 1 and 7-30 for k 2 . Research concerning the structural parameters of collagen in arterial tissue, like fiber orientation and dispersion, is a vibrant field [57,71,[79][80][81][82][83][84][85]. In particular is known that fiber orientation and dispersion may vary considerably in health and disease. To ensure a sufficient scope of our analysis, we thus allowed the initial fraction of collagen and smooth muscle fibers in the diagonal direction to vary in the theoretically maximal range, that is, β t 0 ∈[0-0.5].
Growth and remodeling: the half-life time of collagen is in the range of 60-70 d for healthy aortic tissue [86,87]. It can however change drastically due to a change of mechanical loading or during disease [86,88,89]. Therefore, we consider in our study an extended range of 25-140 d for the turnover time T . Note, that there is a linear dependence between half-life and turnover time by a factor of ln 2. The exact nature of the homeostatic state of soft tissue remains controversial to date [47,[89][90][91]. Thus only little information is available about a reasonable range for σ h . The studies of [92]a n d [ 93] suggest for the homeostatic stress of arterial tissue a range of 150-300 kPa. [47] derived a homeostatic stress range of around 200-300 kPa from theoretical considerations. Motivated by these studies we assume herein for the collagen and smooth muscle fibers alone a range σ h ∈[125-250 kPa].
Experimental data on the gain parameter k σ is very limited. Based on the concept of mechanobiological stability, [11,12] estimated typical values in health and disease. Following these considerations, we assume herein for case 1 (hypertension) k σ ∈[0.12-0.42] and for case 2 (idealized AAA) k σ ∈[0.05-0.15].

Implementation and Discretization
The homogenized constrained mixture model for thin-walled (membranous) anisotropic volumetric growth described in Sect. 3 was implemented in our in-house research code BACI (written in C++) [94]. An explicit time integration scheme is used to solve the evolution equations at each time step (see, for example, Appendix 3 of [13]).
We note that the two cases introduced in Sect. 4.1 exhibit both a reflection symmetry with respect to the cross-sectional plane in the center of the vessel and a rotational symmetry around the cylinder axis. To reduce the computational cost, we exploited these symmetries. That is, we simulated only half of the cylinder in axial direction and only a wedge with an opening angle of 11.25°in circumferential direction. The application of suitable Dirichlet boundary conditions enforcing the respective symmetries enables this reduction. Figure 3 illustrates the reduced computational domain in comparison with the full domain for one exemplary simulation of case 2 (idealized AAA). We discretized the computational domain with 50 standard quadrilateral nonlinear membrane finite elements in axial direction and one in circumferential direction. Within each finite element, the direction vectors of the fiber families were assumed to be constant. In all simulations, we used a timestep size of 10 d.
The implementation of the algorithm to compute the Sobol indices according to Sect. 2.1 has been adapted from the open-source project SAlib [95]. The adapted code was included in the QUEENS code project (written in Python). QUEENS is a general purpose framework for large scale uncertainty quantification and simulation analytics of complex computational models [96]. For each sensitivity analysis, we use N = 6000 Monte Carlo samples which results in a total of 66000 model evaluations for case 1 (hypertension) and 72000 for case 2 (idealized AAA). These can be split into 12000 independent -drawn from the distributions defined in Table 2 -plus 54000 or 60000 cross-sampled samples, respectively (cf., Appendix). Figure 4 shows the probability density functions (PDFs) of the maximum current diameter d max for both case 1 (hypertension) and case 2 (idealized AAA) for three points in time illustrating their evolution in time. The densities are approximated by kernel density estimation (KDE) with Epanechnikov kernels based on the 12000 independent samples of each case study.

Probability Distributions of Model Output
In case 1 (hypertension), the increase in mean blood pressure of 20 mmHg generally leads to minor dilatation of the vessel that largely stabilizes after around 10 years. By contrast, the elastin damage in case 2 (idealized AAA) typically entails a substantial dilatation, which surpasses the dilatation threshold of 3 cm -the clinical criterion for an aneurysm -in more than 18% of cases. A considerable number of simulated aneurysms does not stabilize even after a decade but rather keeps enlarging, which in reality typically results in rupture at some point, if not treated clinically. In total 3.83% of the 72000 simulations performed for case 2 (idealized AAA) numerically failed within the simulated 15 years. A detailed analysis revealed that this phenomenon was exclusively linked to buckling close to the clamped ends of the cylinder. In accordance with the remarks in Sect. 4.1, the missing results of the numerically failed simulations were mapped to d max = 8 cm because in these simulations we also observed an excessive volume of the aneurysm, which is in practice typically associated with excessive stresses and thus rupture.

Sensitivity Analysis
In this section, we present the evolution of the first order and total sensitivity indices of the maximum current diameter d max for each free input parameter over a period of 15 years.  Table 3.Only four parameters have noticeable total indices. These are the turnover time of collagen T , the gain parameter k σ , the initial fraction of total fiber content attributed to each diagonal fiber family β t 0 and the collagen material parameter k 2 . These four parameters can be further separated where k 2 , β t 0 ,andk σ sustain considerably larger, long-term total indices compared to T . The total indices of the remaining five parameters are all below 0.01 and many are practically zero. Therefore, their influence on the variability of d max appears to be negligible compared to the other four parameters. As explained in Sect. 2.1, these five parameters are prime candidates for parameter fixation.
The four indices of considerable magnitude change drastically over time. The turnover time of collagen T influences d max only during the first years. Its total index quickly decreases from 0.325 to 0.045 within the first five years. After 10 years, T has become negligible. With a total index of 0.359, β t 0 is the most influential parameter after the first year. However, its total index decreases almost exponentially and seems to stabilize at approximately 0.125 after 15 years. With this value β t 0 remains the third most influential parameter. The difference between first order and total indices of β t 0 is very small and increases only minimally to a maximum of 0.026 at 15 years indicating that interactions of β t 0 with other parameters are minor. The second most influential parameter is the material parameter k 2 . For k 2 the total index drops from 0.189 after the first years to a value of 0.175 before it rises to 0.197 after 15 years again. In the later years, the first order index of k 2 remains constant such that the influence of the interaction of k 2 with the other parameters increases over time. Figure 5 clearly shows that the gain parameter k σ is by far the most influential parameter except for the first year. From an initial value of 0.201, its total index rises to 0.716 within 5 years. It peaks at 0.795 after 15 years. The evolution of the total index of k σ over time suggests that it converges to a value close to this. With a value of 0.688 for the first order index of the gain parameter k σ at 15 years, an overwhelming amount of the total variance of d max after 15 years can be explained by the uncertainty in k σ alone. From five through 15 years, the interaction terms of k σ account for 10% of the output variance which is approximately equal to the overall amount of interactions (compare to 1 − n i=1 S i ). The majority of interactions occur between k σ and k 2 alone. Ultimately, this reveals, however, that k σ is the cornerstone of the relevant interactions between all parameters.
The sum of all first order indices is shown in the last row of Table 3. At all times, it is n i=1 S i > 0.902, which is close to the theoretical limit of 1.0 indicating that the variability in the model output is dominated by linear terms while interaction between the parameters seem to play a minor role.
Case 2 (idealized fusiform abdominal aortic aneurysm): in the sensitivity analysis of the radial expansion of the idealized AAA, we investigated the influence of 10 parameters on the variability of the model output d max .Table4 summarizes the values of Sobol indices for four selected points in time, namely after one, five, 10, and 15 years. Figure 6 shows the evolution of the sensitivity indices in time over a period of 15 years evaluated annually. The collagen material parameter k 2 , turnover time T and gain parameter k σ all have considerably higher total indices compared to the rest of the parameters. While, some of the less-influential parameters have non-zero total indices of up to 0.149 (β t 0 ) in the first year, their total indices decrease over time; in many cases until they are almost zero. Interestingly, most of the less  In the first years, the variability of the model output is dominated by T as indicated by the comparably large total index of 0.484. The total index of T evolves markedly over time. Overall, it decreases slightly. A minimal value of 0.420 makes T the second most influential parameter at 15 years. By contrast, the first order index of T decreases from 0.360 to 0.117 indicating that T is increasingly involved in higher order interactions. The total index of the material parameter k 2 rises over time from 0.137 to a final value of 0.334 at 15 years. Within 2 years, k 2 becomes and stays the third most influential parameter. However, the first order index of k 2 remains almost constant between 0.074-0.102 indicating that a rise of higher order interactions is responsible for the increase of the total order index. The gain parameter k σ quickly becomes the most influential parameter. Initially, its total index is very small (0.062) but it rapidly grows to a maximum value of 0.693 at 15 years. While its first order index follows this trend, the difference between the two increases considerably over time.
With values between 0.415 at 10 years and 0.370 at 15 years, k σ has the highest amount of interactions in the last five year period. In fact, these values are close to the total interaction values indicating that k σ is involved in nearly all interactions.
Interestingly, interactions seem to occur almost exclusively between the three most influential parameters (k 2 , T and k σ ). Generally, higher order interactions between parameters play a significant role in explaining output variability, in particular in later years, as indicated by the large total interaction value 1 − n i=1 S i in this period.

Input-Output Relations
The results of the sensitivity analysis in Sect. 5.2 reveal that in particular in case 2 (idealized AAA) there remains a substantial fraction of output variability that can only be explained by interactions between the input parameters. Sensitivity analysis can identify which parameters influence the model output but it is beyond its scope to identify how the parameters influence the model output. In general, such an analysis is a very complex task in particular when studying the nature of nonlinear interactions between the parameters. In this section, we use the model evaluations carried out during the sensitivity analysis to partially answer this question. We do so by studying so-called parallel coordinate plots for case 2 (idealized AAA) as depicted in Fig. 7.

Remark 3
For an interactive experience of this section, the reader may refer to the electronic supplementary material. It provides an interactive version of the parallel coordinate plot discussed in this section in a standalone HTML file. The user can constrain arbitrary combinations of parameter and output ranges. Additionally the order of coordinates can be adjusted. The figure was created with the visualization software plotly [98]. Image created by Sebastian Brandstaeter and licensed under the Creative Commons Attribution 4.0 International License, https://creativecommons.org/licenses/by/4.0/.
In Fig. 7a, all 12000 independent samples evaluated are shown together with the corresponding model outputs: each line connects the respective parameter values with their output. The color additionally illustrates the value of the model output. The parameters are arranged by ascending total index after 15 years from left to right. The objective is to study whether specific parameter combinations can be identified as mainly responsible for certain output ranges. From the orange to dark red lines in Fig. 7a, we see that the unifying property of samples that lead to very large model outputs d max ≥ 6 cm is a combination of relatively small turnover time T<100 d and most prominently small gain parameter k σ < 0.1. At the same time, very large values of k 2 > 26 inhibit extreme radial expansions. This is also highlighted by Fig. 7b which shows all samples that lead to aneurysmatic dilatations with d max ≥ 3 cm. Figure 7b shows for example that also large T can lead to aneurysms but only in combination with very small k σ and small k 2 . Albeit not shown here, the same is true for large values of k 2 which may lead to aneurysms but only in combination with very small T and k σ . Indeed, these three parameters have the highest total sensitivity indices with the largest fraction of interactions S T i − S i (see Table 4). The other parameters do not show a clear trend as indicated by the random arrangement of lines on the left side of Fig. 7a. Although less visible from Fig. 7a, it is truly the interaction between T and k σ that leads to elevated d max . Choosing a small value for only one of these two parameters generally does not yet produce a large, aneurysmatic dilatation. Figure 7c illustrates this by restricting the  Table 2. d max is given in cm. The colorbar relates to d max and applies to all figures. Coordinates are arranged according to increasing total index from left to right. (a) shows all 12000 independent samples. The smaller figures extract a subset of samples as indicated by the magenta colored areas; (b) shows all samples that lead to d max ≥ 3 cm, i.e., all samples that can be classified as aneurysms [97]; (c) shows only samples with k σ ≤ 0.075; (d) depicts d max < 3cm;(e) limits the samples to k σ ≥ 0.12  Fig. 4b). The pale lines show the distributions resulting if all parameters except the three most influential ones (collagen material parameter k 2 , turnover time T and gain parameter k σ ) were fixed at their mean values. The good agreement between the distributions confirms the results of the sensitivity analysis: a restriction of the uncertainty problem to only the three most influential input parameters suffices for a good approximation of the output uncertainty range of k σ to k σ ≤ 0.075. We find that, on the one hand, small k σ can result in almost the full range of model output, which means that small k σ is a necessary condition for large d max but not sufficient. On the other hand, large values of k σ ≥ 0.12 alone guarantee smaller radial expansions. Figure 7eshowsthatd max < 3cmfork σ ≥ 0.12 irrespective of the other parameters. The reverse statement is, however, not necessarily true. Figure 7ds h o w sa l l samples with d max < 3 cm indicating again that also small k σ can lead to non-aneurysmatic dilatations if the other parameter values are favorable. The missing lines between small k σ and small T illustrate again that the combination of these always result in aneurysmatic radial expansion.

Validation of Parameter Fixation
As described in Sect. 2.1, the parameter fixation paradigm suggests that parameters with S T i ≈ 0 can be fixed anywhere within the range studied without significantly affecting the model output. Thereby, the complexity of the uncertainty quantification problem can be reduced substantially. In this section, we compare the PDFs of the maximum current diameter for case 2 (idealized AAA) resulting from a parameter fixation, with the fully resolved results as presented in Sect. 5.1. Specifically, we fixed all parameters that were identified with negligible total indices to the mean value of their respective probability distributions defined in Table 2. We then evaluated the model for 12000 independent samples from the distributions of the remaining 3 parameters, that is, collagen material parameter k 2 , collagen turnover time T and gain parameter k σ . Figure 8 shows the KDE of the resulting distributions of the maximum current diameter d max after 5, 10 and 15 years compared to the fully resolved ones with 10 uncertain parameters. As predicted by the sensitivity analysis, the distributions of the reduced uncertainty model are indeed very good approximations of the full model in particular after 10 and 15 years. The densities depicted in Fig. 8 agree very well. Both the mean values and the variances are almost identical. Due to the quasi-random nature of the Sobol sequences used to draw the samples (cf., Appendix), we were able to directly relate the samples of the full model with the ones from the reduced model. This enabled the computation of the mean relative error of d max (t = 15 y) caused by the fixation of the non-influential parameters at their mean values. This mean relative error was computed as 6.54%. Overall, these results justify the fixation of the 7 non-influential parameters with small total indices as predicted by the parameter fixation paradigm.

Discussion
The choice of our examples allows a systematic interpretation of our results. Both in case 1 (hypertension) and case 2 (idealized AAA), an inelastic deformation by growth and remodeling is triggered by some perturbation. However, while in case 1 our system quickly converges to a new steady state after an initial period of growth and remodeling, in case 2,a continued, unbounded enlargement of the vessel is observed in many cases. That is, in case 1, we study a mechanobiologically stable system in the sense of [11,12], whereas in case 2, we study a potentially mechanobiologically unstable system.
In both cases, the turnover time T determines the time scale of growth and remodeling. Therefore, it substantially affects the model output in both cases directly after the perturbation when the system is subject to highly dynamic growth and remodeling. In case 2,t h i s remains so for the whole simulated period at least for mechanobiologically unstable vessels, which are subject to a continued process of growth and remodeling and which, as a subclass of the vessels studied in case 2, have a large impact on the overall variance of the model due to their unbounded enlargement. By contrast, in case 1 the influence of T quickly declines over time. Because after a while all the vessel attain a new stable mechanobiological equilibrium configuration, which depends on k σ but only to a minor extent on T . The latter can be seen from combining Eq. (79) and Eq. (82) in [11]. This reveals that in mechanobiologically stable blood vessels the residual deformation, that remains in the long term after a perturbation, depends only to a relatively small extent on the turnover time T (for parameter values typical for the aorta). In short, T is typically a parameter of high impact on the model output as long as a blood vessel is subject to an ongoing growth and remodeling dynamics. Therefore, it can be expected to be of particular importance for clinical predictions of aneurysmal enlargement.
The two only elastic parameters of major importance for the model output in our study are β t 0 and k 2 . The former describes the orientation of the collagen fibers in the vessel wall.
In both examples, it was found to be of major importance directly after the perturbation of the systems. In this very early stage, elastic deformation is still much larger than inelastic deformation due to growth and remodeling. Therefore, it is not surprising that the orientation of the collagen fibers, which substantially affects the elastic properties of the vessel wall, has a large impact on the model output. However, the larger the ratio between inelastic and elastic deformation of the vessel, the less important becomes β t 0 .I ncase 1, it can always retain some importance because the overall deformation due to growth and remodeling remains small. However, in mechanobiologically unstable systems such as the one studied in case 2, the substantial inelastic deformation due to growth and remodeling soon makes the impact of β t 0 negligible. This indicates that β t 0 is an important parameter for the elastic deformation of the system but not for its growth and remodeling dynamics. The latter is affected only by one elastic parameter, k 2 , which determines the strain-stiffening of the tissue. Large k 2 are associated with substantial strain-stiffening, which apparently can effectively reduce the inelastic deformation due to growth and remodeling.
Both in the early and late stage of case 1 and case 2, k σ plays a dominant role. Again, this is coherent with the theory of mechanobiological stability [11,12], which reveals that this parameter has not only a major impact on the rates of inelastic deformation during growth and remodeling but also on the steady state itself that is reached in mechanobiologically stable systems in the long term.

7C o n c l u s i o n s
We presented a variance-based global sensitivity analysis for the homogenized constrained mixture model of arterial growth and remodeling in two case studies. In case 1 (hypertension), we investigated the adaptation of an idealized vessel to an elevated mean blood pressure level. In case 2 (idealized AAA), we considered radial expansion of an idealized vessel induced by spontaneous damage to the elastin component. For each case, we studied how the uncertainty in the model input parameters contributes to the uncertainty in the chosen model output, that is, the maximum current diameter of the vessel.
A striking difference between case 1 and case 2 is the fact that in the latter, interactions between the parameters were observed to play a much larger role. This underlines the necessity for global sensitivity analysis and suggests that previous simple parameter studies where only one parameter at a time was varied [30,33,53] could reveal in principle only a limited part of the parameter sensitivities.
As discussed in Sect. 6, the results of our simulations can be interpreted in a systematic way. In the cases studied herein, only three model parameters were found to have a major impact on the inelastic deformation of the blood vessel due to growth and remodeling. These parameters were k σ , T and k 2 . While k σ can be identified with the ability of the tissue to increase collagen production as stress in the vessel wall increases, T is the average life time of collagen fibers, directly linked to their half-life time. The elastic parameter k 2 describes the strain stiffening of the collagen tissue.
Our results may have important implications both for future computational studies and for the directions that appear promising in clinical research.
For future computational studies of growth and remodeling, they mean that it can be acceptable to fix many of the parameters to values close to the population mean reported by the literature without bothering about case-or patient-specific values. This can substantially simplify the design of future computational studies and save resources in parameter studies.
It remains an important goal of clinical research to predict which enlargement can be expected for specific aneurysms. Our global sensitivity analysis clearly suggests that significant progress could be made if ways are found to measure or at least estimate the ability of the vascular tissue to produce collagen and ideally also the collagen half-life time. At the moment there are no clinical imaging protocols available to measure both. Our study suggests that developing such protocols, for example on the basis of functional magnetic resonance tomography, might be one of the most effective steps that could be taken to improve our ability to predict the future enlargement of aneurysms. Additionally, model predictive capabilities might be enhanced by improving the measurement of collagen material properties under large deformation, which, however, may be non-trivial in a non-invasive manner.
While this seems to indicate a promising direction of future research, also some limitations of the presented study should be mentioned. First, only scalar output quantities can be investigated with the variance-based method chosen herein. Hence, we focused on the maximal diameter as the single output quantity of interest. However, input parameters that have little effect on one output quantity might be very important for another one. To overcome this limitation of our study one may combine in the future multiple scalar outputs as suggested in [44] or use more sophisticated methods as suggested in [99][100][101]. Next, it is important to underline that our study is based on simplified physical model assumptions neglecting the influence of blood flow and associated wall shear stresses. For larger, elastic arteries such as the aorta considered herein, also several other groups have apparently considered this simplification acceptable [6,19,30,53,67]. However, at least for smaller arteries such as cerebral arteries blood flow and wall shear stress are frequently discussed as key quantities [5,24,102,103]. In the future, more detailed work should therefore aim at including also such effects. It should also be noted that in clinical practice many larger AAAs contain intraluminal thrombi, which were not considered herein. Their effect on aneurysmal expansion is, however, often not thought to be primarily a mechanical one [104][105][106]b u t rather a biochemical one [21,29,56]. That is, it should translate into changes of the parameter T and k σ in our model and is thus at least implicitly covered by this study to some extent. Nevertheless, an extended sensitivity analysis on the basis of a model including the intraluminal thrombus explicitly would be required to make strong statements about its role.
Finally, we note that, as is the nature of global sensitivity analysis, our study focused on but one computational model of vascular growth and remodeling, the homogenized constrained mixture model. Therefore, additional research is needed to corroborate the results of our study and to ensure that the conclusions drawn here are not artifacts of a specific model but are in fact revealing important aspects of the real physiology.

Appendix: Algorithm for the Computation of First and Total Order Indices
The computation of an estimate for the first order and total order Sobol indices can be carried out by the following algorithm which is similar to the one presented in [44] In other words, the a j i and b j i are independently drawn samples from U(0, 1). In practice, instead of using pseudo-random numbers one uses quasi-random samples from lowdiscrepancy sequences such as the Sobol sequence for improved performance [107][108][109]. In the following, we denote the j -th row of A and B by a j and b j , respectively.

Based on
3. Evaluate the model y(x) at every row of the (n + 2) matrices {A, B, A (1) B , A (2) B ,...,A (n) B }. This results in a total of N tot = N(n+ 2) model evaluations. In the following, we denote by y(A) a vector whose j -th component is y(a j ),t h a ti sy evaluated at the j -th row of A. Likewise, we define the vectors y(B) and y(A (i) B ). Note that keeping the ordering between the sampling matrices and the result vectors is essential. For ease of notation, we introduce the vector y(C) ∈ R 2N denoting the concatenated vectors of y(A), y(B), that is 4. Based on the above model evaluations the following estimators for the first-order and total-order sensitivity indices can be computed: and where the total variance of y is estimated with Supplementary Information The online version contains supplementary material available at https://doi. org/10.1007/s10659-021-09833-9.