Abstract
The abdominal aorta (AA) in older individuals can develop an aneurysm, which is of increasing concern in our ageing population. The structural integrity of the ageing aortic wall, and hence aneurysm, depends primarily on effective elastin and multiple families of oriented collagen fibres. In this paper, we show that a structurally motivated phenomenological ‘fourfibre family’ constitutive relation captures the biaxial mechanical behaviour of both the human AA, from ages less than 30 to over 60, and abdominal aortic aneurysms. Moreover, combining the statistical technique known as nonparametric bootstrap with a modal clustering method provides improved confidence intervals for estimated bestfit values of the eight associated constitutive parameters. It is suggested that this constitutive relation captures the wellknown loss of structural integrity of elastic fibres owing to ageing and the development of abdominal aneurysms, and that it provides important insight needed to construct growth and remodelling models for aneurysms, which in turn promise to improve our ability to predict disease progression.
1. Introduction
Abdominal aortic aneurysms (AAAs) are focal dilatations of the infrarenal aorta that typically have a fusiform shape and often progress to rupture, which is increasingly responsible for mortality and morbidity in our ageing population. AAAs are characterized histologically by a loss of medial elastin and smooth muscle cells, mechanically by an increased stiffness and decreased distensibility, and clinically by a complex shape, a maximum diameter over 150 per cent of normal and often the presence of an intraluminal thrombus [1–3]. Biomechanics plays fundamental roles in both the enlargement and rupture of these lesions. In particular, computations suggest that aneurysms enlarge, in part, because remnant intramural cells remodel the wall in an attempt to restore intramural stresses towards normal values despite the marked alterations in geometry and haemodynamic loads (cf. [4–7]). Computations also suggest that wall stress is a better predictor of rupture potential than the commonly used clinical metric of maximum diameter [4,8]; rupture occurs when stress exceeds strength locally, which probably results from an imbalance in the production (synthesis) and removal (degradation) of fibrillar collagen. Nevertheless, all prior computational models of stresses in AAAs are limited by diverse simplifying assumptions, including neglect of spatial variations in the evolving mechanical properties of the wall.
The goal of this paper is to identify a single, structurally motivated phenomenological, nonlinear constitutive relation that describes the biaxial mechanical behaviour of ageing human abdominal aorta (AA) and AAAs. Towards that end, we employ a ‘fourfibre family’ constitutive relation that is motivated by microscopic data on the organization of arterial collagen [9] and that captures differences in the mechanical behaviour of carotid arteries that result from diverse structural abnormalities in genetically modified mice [10,11], provides new insight into origins of residual stress and axial prestretch [12], and describes growth and remodelling in response to altered haemodynamics [13]. Moreover, we show that combining the technique for statistical inference known as nonparametric bootstrap with a modal clustering method yields improved confidence intervals for bestfit values of the model parameters, particularly given the possible nonuniqueness in nonlinear estimations.
2. Material and methods
2.1. Constitutive formulation
We reanalysed biaxial data on the mechanical behaviour of human AA and AAA reported by Vande Geest et al. [14,15]. Briefly, they excised nearly planar specimens from the anterior portion of the infrarenal aorta at open surgical repair (26 AAA samples) or autopsy within 24 h of death (21 AA samples). Specimens were loaded inplane along axial (z) and circumferential (ϑ) directions such that the ratio of biaxial tensions was constant at prescribed values (0.5 : 1.0, 0.75 : 1.0, 1.0 : 1.0, 1.0 : 0.5 or 1.0 : 0.75). Because of the lack of any reported hysteresis, we used the loading portion of the stress–stretch curves following preconditioning.
In contrast to prior approaches, we modelled these specimens as a constrained mixture of an elastindominated amorphous matrix and four families of locally parallel collagen fibres (i = 1 = axial, 2 = circumferential and 3 and 4 = symmetric diagonal). Because of the experimental protocol followed by Vande Geest and colleagues, we neglected possible active stresses resulting from smooth muscle tone; possible passive contributions by smooth muscle (cf. [16]) were absorbed naturally by the circumferential family of fibres. Hence, consider a strain energy function of the form [17] 2.1where c, and are material parameters, I_{C} = tr C is the first invariant of the right Cauchy–Green tensor and represents the square of the stretch of the i^{th}fibre family (i.e. (λ^{i})^{2}). The fibre orientations are defined in a convenient reference configuration by the unit vectors M^{i}, which depend on angles defined between directions of fibre reinforcement and the axial direction. Hence, axial and circumferential fibres are fixed at and , respectively, while the diagonal fibres are accounted for by . Thus, 2.2where superscript d denotes ‘diagonal’. Diagonal families of collagen are typically regarded as mechanically equivalent, hence . This fourfibre family constitutive model thereby requires estimation of eight unknown parameters . Multiple special cases can be recovered, as, for example, by allowing axial and circumferential families of fibres to be mechanically equivalent (cf. [9]), which reduces the number of parameters to be determined via nonlinear regression to six.
Assuming a homogeneous deformation within the central region of biaxially tested specimens of the form F = diag (λ_{z}, λ_{ϑ}, λ_{r}), and assuming both incompressibility λ_{r} = 1/(λ_{z} λ_{ϑ}) and plane stress (t_{rr} = 0), this stored energy function yields Cauchy stresses of the form 2.3and 2.4which enables a straightforward analysis of biaxial data.
Motivated by Gasser et al. [18], we also considered a strain energy function that accounts for possible dispersion of collagen fibres about the four primary fibre directions. Each fibre family is then characterized by a structural tensor H_{i} that represents distributed collagen fibres in a continuum sense, with 2.5where κ is a scalar parameter representing fibre distributions in an integral sense. The associated Greentype strain experienced by the ith family of fibres, having mean orientation M^{i}, is characterized by E_{i} = H_{i} : C − 1, namely 2.6Thus, the second strain energy function considered herein has the form, 2.7from which principal Cauchy stresses can be determined similar to equations (2.3) and (2.4). Note that mathematically admissible values are κ ∈ [0,1/2] [19], where κ = 1/3 models an isotropic distribution of fibres and κ = 0 recovers the original ideal preferential alignment [18]. Hence, equation (2.7) is a natural extension of our fourfibre family model (equation 2.1), which is recovered as a special case. Whereas Gasser et al. [18] used a single dispersion parameter κ to characterize distributions of two families of fibres, Pandolfi & Holzapfel [20] suggested that each fibre family can have a separate dispersion (e.g. κ_{1} and κ_{2} for axial and circumferential families and κ_{3} = κ_{4} = κ_{3,4} for diagonal families). Introducing either a single or multiple dispersion parameters, equation (2.7) contains nine or 11 parameters. Unlike the material parameters and , the structural parameters α_{0} and κ (or κ_{i}) can be identified either by fitting experimental data that describe macroscopic mechanical responses or from an associated histological analysis. We adopted the former approach owing to the lack of suitable histological data for the samples studied and because the model is only motivated by the underlying structure, not based on inherent microstructural complexities such as different types of collagen, crosslinking, physical entanglements or other interactions with the many other extracellular constituents.
2.2. Multivariate nonlinear regression
Specific inequality constraints must be respected when estimating constitutive parameters to ensure that bestfit values are physically realistic [21]. For example, the eight parameters associated with equation (2.1) must satisfy plus 0 ≤ α_{0} ≤ π/2 owing to symmetry of the diagonal fibres [22]. The same restrictions hold for equation (2.7), which has the additional constraint 0 ≤ κ_{i} ≤ 1/2.
Regardless of the particular constraints, one must minimize simultaneously the differences between ‘theoretically predicted’ and ‘experimentally measured’ Cauchy stresses in the two stretching directions at each equilibrium configuration k = 1,2, … , m. Hence, let the vector λ collect the n = 2m (i.e. biaxially) measured values of stretch within the central region, namely 2.8with λ_{z} and λ_{ϑ} denoting all values in z and ϑ, respectively, for all test protocols (i.e. tension ratios). Similarly, let be the corresponding vector containing all Cauchy stresses calculated from the experimentally measured applied forces and geometry. The vector of all theoretically predicted Cauchy stresses, for all measured stretches λ, can be written t^{th} = t^{th}(λ;θ), where represents the set of p constitutive parameters to be estimated. The objective function e is then defined as [23], 2.9with the optimization problem given by min_{θ} e(θ). Because the theoretically predicted stress t^{th}(λ; θ) is a nonlinear function of the parameter vector θ, the regression analysis was performed using the builtin function lsqnonlin in the Optimization Toolbox of Matlab.
Initial guesses of the values of the unknown parameters were assigned by generating a vector of p uniformly distributed pseudorandom numbers every time the regression code was executed; these guesses were restricted by lower and (when appropriate) upper bounds to ensure a physically meaningful search space. The regression algorithm stopped when either TolFun (the tolerance placed on the objective function) was less than 1e − 10 or TolX (the tolerance placed on the estimated parameter values) was less than 1e − 8; these values were selected based on preliminary convergence studies. As a measure of goodness of fit, we used the traditional coefficient of determination R^{2} ∈ [0,1]; a value of R^{2} ≥ 0.9 is usually an indicator of a good fit. Because we compared different versions of the model having different numbers of parameters, we also used the corrected Akaike information criterion (AIC_{c}; [24]), 2.10AIC_{c} converges to the AIC [25] as n becomes large, hence it can be employed regardless of sample size [26]; the preferred model is often the one with the lowest value of the AIC_{c}. This criterion for model selection was used previously by ZeinaliDavarani et al. [27] to compare models with an increasing number of fibre families (from 2 to 10).
2.3. Nonparametric bootstrap
In addition to variations that arise in arterial properties owing to common adaptations, ageing or disease, inherent specimentospecimen variations exist even within gender and agematched samples [21]. It is often challenging, therefore, to identify correlations between biological processes and mechanics; robust parameter estimation becomes fundamental to such analysis of data.
Confidence intervals for estimated parameters can be important measures of the precision of a nonlinear regression [28]. We used the nonparametric bootstrap method and a standard technique test both to estimate confidence intervals and to reveal possible nonuniqueness of the bestfit values of the parameters (figure 1). The bootstrap [29] uses ‘resampling with replacement’ of original data, based on an empirical distribution function, to estimate a parameter's probability distribution, which in turn can be used to estimate a confidence interval. To apply the nonparametric bootstrap to our multivariate regression of biaxial data, we indicate by x = (x_{1}, x_{2}, … , x_{n}) a dataset containing n ‘points’, where . Data collected in x are bootstrapped (i.e. resampled) B times so that each resampled dataset, indicated by , contains the same number of observations n as the original dataset but consists of a randomized version of x (in nonparametric bootstrapping, each pair x_{j} has the same probability of being resampled and thus may be included one or more times, or not at all). The constrained nonlinear regression is then executed for each bootstrap dataset, again starting from a randomly generated initial guess, thus yielding an estimate , called a bootstrap replication of , the vector that collects estimated values of the p parameters. Each of the estimated parameters is thus obtained as a collection of B bootstrap replications that is used to approximate the probability distribution for that parameter. Confidence intervals are determined from the percentiles of these distributions.
Without a loss of generality, henceforth consider a representative parameter θ, which can be any component of the vector , with the estimate of θ based on regression of experimental data and obtained via bootstrapping.
2.4. Modality testing and clustering
Possible clustering of parameter replications around different values, which reflect local minima in the objective function, complicates the construction of confidence intervals. That is, if clusters of values around nonglobal minima are not eliminated from the parameter replications that are used to calculate confidence intervals, the resulting intervals will underestimate the precision of the bestfit parameters and negate advantages of nonparametric inference based on the bootstrap. To address this problem while keeping the bootstrap algorithm fully automatic and consistent with a nonparametric estimation framework, we implemented a modal clustering algorithm.
From the histogram of replicates, it may or may be not clear whether the distribution of has a single mode or multiple modes. Hence, we used a standard (Silverman's) modality test based on the kernel density estimate (KDE) of the probability density function. Applied to the present case, given by B bootstrap replications of the estimated parameter θ, the KDE is given by 2.11where K(·) is the kernel function and h is a smoothing parameter called bandwidth. Larger values of h produce a smoother density estimate with fewer modes. Silverman [30] proposed the use of the KDE and bootstrap for testing multimodality, and showed that if the kernel is a standard normal density function, , then the number of modes is a right continuous decreasing function of h. Briefly, in testing the null hypothesis that the true density f has exactly q modes against the alternative that f has more than q modes, one can use as a test statistic the critical bandwidth h_{crit}, which is defined as the smallest bandwidth yielding a KDE with at most q modes and is obtained from bootstrap replications of . This critical bandwidth places the KDE exactly on the boundary between the null and alternative regions; thus larger values of h_{crit} argue against the null hypothesis. Within this framework, the bootstrap is used to approximate the null distribution of h_{crit} and thereby determine the significance of , the critical bandwidth obtained from the bootstrap replications of the parameter . This procedure is called a smooth bootstrap because, instead of sampling with replacement from the set of parameter replications, the sample is drawn from a smooth and rescaled density estimate (for details, refer to the original paper). To reduce computational cost, it is common to assess statistical significance by applying the density estimator with bandwidth to each bootstrap sample and then to calculate the pvalue as the proportion of the resulting density estimates with more than q modes. The null hypothesis that f has q modes is rejected at the confidence level α (herein set to 0.05) if p ≥ 1 − α [31].
Because we did not know a priori the number of modes of the parameter distributions, we tested sequentially the existence of up to nine modes in the distribution of each parameter and, among the number of modes that were not rejected, chose the one corresponding to the highest statistical significance. Once the number of modes was known, we then performed a clustering analysis to identify distributions corresponding to local minima and obtained the distribution of parameters centred about the global minimum to allow an easier calculation of the confidence intervals. The clustering algorithm was implemented via the function cluster in the Statistic Toolbox of Matlab; it builds clusters using the command pdist on the ‘seuclidean’ option, which stands for standardized Euclidean distance (each parameter's sum of squares is divided by its sample variance). In this way, clustering does not depend on measurement units for the parameters (kPa for c and , radians for α_{0}, dimensionless for and κ). This grouping was realized using the command linkage. In a minority of cases, when the two (or more) modes detected by the modality test were close enough to overlap partially, the average linkage method performed better. In these cases, setting the linkage function on the ‘average’ option defined the distance between two clusters as the average of distances between all pairs of objects.
Finally, the group that contained the mode with the largest area in the probability distribution was considered to be most representative of the global minimum of the objective function with respect to that particular parameter. At this point, the ‘filtered’ set of parameter replications (i.e. those falling within the selected group) was used to estimate the confidence intervals. Hence, the actual number of bootstrap replications used to approximate confidence intervals was reduced from B to B_{red}, the number of points within the selected cluster. To distinguish these new parameter replications from the original bootstrapped ones, we indicate them with where r denotes values within a selected subset of {1,2, … , B}.
2.5. BC_{a} confidence intervals
Following DiCiccio & Efron [32], the BC_{a} (biascorrected and accelerated) method approximates confidence intervals for a parameter θ from empirical percentiles of the bootstrap empirical distribution . Compared with the standard percentile method [29], the BC_{a} method corrects the interval for any bias and skewness of the empirical distribution, thus yielding smaller approximation errors. The upper endpoint of a onesided levelα BC_{a} interval is defined by 2.12where represents the cumulative distribution function (CDF) of B_{red} bootstrap replications , and Φ is the standard normal CDF, with z^{(α)} = Φ^{−1} (α). The percentiles that determine endpoints of the BC_{a} interval depend on two numbers, a and z_{0}, called acceleration (or sometimes skewness correction) and bias correction, respectively. The easiest implementation of the BC_{a} algorithm estimates these two constants as 2.13and 2.14where is the jackknife replication of , that is, an estimate of θ based on the reduced dataset x_{(}_{j}_{)} = (x_{1}, x_{2}, … , x_{j}_{−1}, x_{j}_{+1}, … , x_{n}), and . Note that measures the median bias of and measures the rate of change of the standard error on a normalized scale [29]. Using equation (2.12), with nonparametric estimates of and , the central 95 per cent BC_{a} confidence interval is given by 2.15
The advantage of using BC_{a} confidence intervals over conventional intervals is that result (2.15) has two fundamental properties: it is secondorder accurate and it is rangepreserving, which is particularly important because we constrained the parameter search space, and in some cases, the closeness of the estimated parameter to zero would have otherwise lead to a negative lower endpoint (physically unrealistic) when using conventional methods.
2.6. Bootstrap sample size
The number of bootstrap replications required to estimate the 95th percentile of the parameters' distributions is high because the percentile depends on the tail of the distribution where fewer samples occur. Usually B = 1000 provides accurate estimates, but B should be twice as large if the biascorrection constant z_{0} is estimated from the bootstrap distribution [33]. Hence, we estimated the BC_{a} confidence intervals with B = 2000. To implement the smooth bootstrap to estimate the null distribution of h_{crit}, we used B_{smooth} = 500 ([29]; §16.5). Preliminary simulations of 2000 bootstrap replications of biaxial data with 500 smooth bootstrap samples drawn from the KDE of each parameter revealed computational times of the order of 20 min per sample using a standard desktop PC (Intel Pentium E2180, 2.0 GHz).
2.7. Subjective metrics
The patient population from which samples were obtained by Vande Geest et al. [14,15] was very diverse, including differences in age and gender, two of the main factors in determining the pathophysiology of AA. It was difficult, therefore, to compare models based solely on average AIC_{c} and R^{2} values. Below, we discuss the observed results and motivations behind the final selection of the constitutive model.
3. Results
The progressive increase in complexity going from standard fourfibre family models, having six or eight parameters, to those allowing dispersed collagen fibres, having nine or 11 parameters, lead to a slightly better description of the AA data, as revealed by slowly decreasing values of AIC_{c} and comparable mean values of R^{2} (results not shown). This observation suggested that four key versions of equations (2.1) and (2.7) were able to describe data from AA specimens, with little improvement achieved by increasing the number of parameters. The situation was very different for AAA specimens, however. There was a marked reduction in goodness of fit from the eightparameter model (mean R^{2} = 0.90 ± 0.10, AIC_{c} = 1010.6 ± 501.9) to the associated nineparameter model obtained by adding a single dispersion parameter κ for all fourfibre families (mean R^{2} = 0.76 ± 0.16, AIC_{c} = 1180.6 ± 502). This diminished descriptive capability was reflected by poor fits (not shown). Indeed, little improvement was obtained by allowing different dispersion parameters for the different families of collagen fibres, that is, using the 11parameter model (mean R^{2} = 0.86 ± 0.12, AIC_{c} = 1076.6 ± 508.5). Regardless of the diminished goodness of fit observed after introducing dispersion parameter(s), the decision to not consider equation (2.7) further was made primarily because in the few cases of acceptable results for the AAA data (and similarly in several AA datasets), bestfit values of κ (or κ_{i}) were less than 1e − 10, suggesting that better fits resulted when considering only ideally aligned fibres. This result is consistent with experimental observations of Wicker et al. [9] that collagen fibres are organized in the plane defined by the axial and circumferential directions with few to no fibres in the radial direction (perhaps because the artery does not need to carry tensile load in that direction). This type of inplane fibre dispersion is not described by a distribution having rotational symmetry (cf. [18]), which is consistent with the small values of κ_{i} that we found.
Allowing axial and circumferential families of fibres to be mechanically equivalent (cf. [9]) lead to a simpler sixparameter model that fit AAA data adequately (mean R^{2} = 0.88 ± 0.10, AIC_{c} = 1061.8 ± 490.7). Because this goodness of fit was similar to that for the standard eightparameter model (mean R^{2} = 0.90 ± 0.10, AIC_{c} = 1010.6 ± 501.9), we compared these two models further by using the bootstrap method to identify possible crosscorrelations between the parameters that were forced to be equal (i.e. and ). No strong correlations were found, thus it was concluded that the standard eight parameter version of the fourfibre family model gave an overall better description of data from both AA and AAA.
AA data were grouped based on age of the donors [14], with five subjects less than 30 years of age (group 1), eight between 30 and 60 years of age (group 2) and eight over 60 years of age (group 3). The eightparameter fourfibre family model fit data from group 1 very well (figure 2 and table 1) while providing a good description of data from groups 2 and 3 (tables 2 and 3). Notwithstanding high standard deviations for individual parameters within the three groups, owing to aforementioned diversity in the study groups and hence data, note that the mean value of the neoHookean parameter c decreased with increasing age of subjects.
AAA data displayed even higher variability in both mechanical response and estimates of the bestfit values of the parameters (table 4), as might be expected because of the complexity of the disease progression and the inability to know the natural history in individual patients. In addition, these data were subject to high measurement noise (see original plots in [15]) because the marked increase in stiffness reduced displacements of the tracking markers during the experiment to be of the order of the resolution of the video system. Nevertheless, important observations can be inferred from table 4 and the associated stress–stretch plots. For example, in all but one case, the bestfit neoHookean parameter c was practically zero. Moreover, as reported in the original paper, the AAA specimens exhibited greater anisotropy with preferential stiffening in the circumferential direction, which was captured by the fourfibre family model (figure 3). Hence, it appears that the microstructural motivation of this form of the strain energy function enables it to describe both the slightly sigmoidal behaviour of a nearly isotropic and highly extensible specimen from a young donor (figure 2) and the stiffer anisotropic response of an aneurysmal tissue excised from a significantly older subject (figure 3). That the model captured the progressive changes in biaxial mechanical behaviour of the AA and subsequent aneurysm is revealed well by representative fits to data for four different samples (figure 4). Finally, median values of the bestfit parameters represented well the average mechanical behaviour of each group (results not shown). Conversely, mean values of the parameters provided conservative estimates of behaviour, favouring the stiffer specimens.
Differences in estimated values of model parameters for groups of AA and AAA were confirmed by the bootstrap analysis. Figure 5 shows calculated distributions of parameter values and associated BCa confidence intervals. The primary and most striking difference surfaced for the nonnegative neoHookean parameter c, which decreased with increasing age and went toward zero for AAAs. The observed distribution of bootstrap replications of c confirmed that negative values were excluded from the allowed parameter space. Another general tendency was seen in the distributions of collagenrelated parameters: tended to be larger than for the young AA specimens, whereas tended to be smaller than for the older AA specimens and AAAs (tables 1–4). This tendency was confirmed by visual comparison of bootstrap distributions, which further revealed a higher variability in AAA parameters and the presence of multiple local minima for the AAAs compared with normal AAs. Albeit not shown, there were no apparent correlations between parameter sets and AAA patient characteristics such as gender or diameter of the lesion. Hence, the only overall trend in parameter values was that observed with age and the presence of an aneurysm.
4. Discussion
Vande Geest et al. [14] reported biaxial mechanical data for normal human AA as a function of age: less than 30, between 30 and 60 and over 60 years of age. They suggested that specimens in the youngest group were well described (R^{2} = 0.97) by a threedimensional (purely phenomenological) isotropic stored energy function of the form 4.1where I_{C} = tr C and D_{i} are material parameters. Specimens from both groups older than 30 years of age were well described (R^{2} = 0.82 and 0.90 for the middle and oldest groups, respectively) by a twodimensional (purely phenomenological) anisotropic Fung relation, 4.2where E_{ΘΘ} and E_{ZZ} are components of the Green strain tensor E = (C − I)/2 and c, A_{i} are material parameters. Although they noted that a ‘limitation of the current study is the use of two separate strain energy functions to model the same tissue’, Vande Geest et al. [15] subsequently suggested that specimens from both the oldest group of AA and AAAs (mean ages being 70.6 and 72.3 years old, respectively) were best modelled using yet another twodimensional (purely phenomenological) anisotropic stored energy function, namely the Choi–Vito form 4.3where b_{i} are material parameters. Fits to data were very good for normal aorta (R^{2} = 0.95) and similar for the 26 aneurysmal specimens (R^{2} = 0.90). In this paper, we showed that a single (structurally motivated, phenomenological) stored energy function provided a better description of the same data for normal AA at all three age groups (R^{2} of 0.996, 0.949 and 0.958 for youngest to oldest). Our results (tables 1–3) revealed a progressive decrease in the material parameter c (associated with the isotropic contribution attributed to an elastindominated amorphous matrix) and increase in the parameters and (associated with the predominant families of collagen fibres), each consistent with the wellknown progressive stiffening of the aorta with age that associates with increased collagen to elastin ratios [34]. Moreover, this same stored energy function fit data from all 26 AAA specimens well (R^{2} of 0.899), and suggested a further loss of elastin effectiveness and increased stiffening owing to collagen.
Many of these biaxial data have also been studied by other investigators. Rodriguez et al. [35] compared the ability of an isotropic (purely phenomenological) Demiraytype model and an anisotropic (structurally motivated, phenomenological) Gassertype model to fit data from a representative AAA reported by Vande Geest et al. [15]. The anisotropic model fit the data better than did the isotropic model, as expected, with an overall fit (R^{2} = 0.86) similar to that achieved by the Choi–Vito model despite a visually less appealing fit (see their fig. 3). Rodriguez et al. [36] contrasted these findings with the ability of an isotropic polynomial stored energy function and a twofibre family model to fit data, though fits were reported for different (uniaxial or biaxial) data. The best fit to data for the same aneurysm (R^{2} = 0.89) was achieved by allowing the fibre orientation (in the twofibre family model) to be found via regression, although similar results were obtained by fixing the fibre orientation based on histological data from another study.
Rissland et al. [37] fit the mean AAA data reported by Vande Geest et al. [15] using a modified twofibre family model (combining Demiray and polynomial functions for the isotropic contribution). They emphasized that the fit to data was much better with their eightparameter anisotropic model than with a previously used isotropic polynomial model [38], but again the fit to data was not visually appealing (cf. their fig. 2) despite a reported R^{2} = 0.999. Basciano & Kleinstreuer [39] proposed a threedimensional (phenomenological) anisotropic stored energy function of the form 4.4where I_{C} = tr C, IV_{C} = M^{1} · CM^{1}, VI_{C} = M^{2} · CM^{2}, with M^{1}, M^{2} denoting unit vectors that define twofibre directions similar to the original twofibre family model [22], and α and β are material parameters. Mean values of the bestfit parameters for a single sample were α = 74.3 kPa β = 3.30 MPa for AA and α = 74.0 kPa β = 65.5 MPa for AAAs, with the diagonal fibres oriented approximately 45° from circumferential in both cases. Although the reported fit was very good (e.g. ; cf. their fig. 3), note that the parameter for the leading ‘isotropic’ term (cf. [40]) was unchanged despite elastin being significantly reduced in AAAs. Moreover, most data show preferential circumferential stiffening of AAAs, suggesting that collagen becomes progressively oriented more towards the circumferential direction (cf. [6] and figure 4). That is, although this phenomenological model described one dataset well, it does not offer interpretive value with regard to histological changes. Most recently, Haskett et al. [41] employed a variation of the twofibre family model proposed by Gasser et al. [18], including fibre dispersion, to fit biaxial data from multiple regions of the human aorta at multiple ages (from 0–30, 31–60 and 61 years old and above). They did not report R^{2} values or show graphical fits, however.
In this paper, we suggested yet another form of the stored energy function (equation (2.1)), one motivated directly by microstructural information on arterial collagen obtained using multiphoton microscopy [9], by excellent prior fits to data from diverse arteries [10,11,17], by new interpretive value related to residual and prestresses [12] and by utility in growth and remodelling theories [13,42]. We submit that this functional form fits available human biaxial data as well as or better than prior models (e.g. mean R^{2} ranging from 0.899 to 0.996 for all four available groups from AAs and AAAs), but with added advantages of a structural motivation without unnecessary complexity. That is, although actual orientations of collagen fibres probably range smoothly from axial to circumferential, predominantly within constant radial planes (cf. [9]), addition of rotational fibre dispersion about the primary orientations did not improve the fit to data enough to justify increasing the number of structural parameters. Rather, it appears that the primary advantage of the fourfibre family model (cf. [17]) over the original twofibre family model [22] is explicit inclusion of axially and, especially, circumferentially oriented families of fibres. For example, loss of elastin coupled with either reorientation or remodelling of collagen towards the circumferential direction would render the twofibre family model less capable of capturing the observed marked axial stiffness of AAA specimens (cf. figure 4). The fourfibre family model is similarly expected to be free of limitations of the twofibre family model identified by Gasser et al. [18] via simulations of uniaxial stretching tests, but this was not tested directly.
As should be expected from the consistent histopathological finding of reduced elastin in ageing [34,43] and AAAs (e.g. [44,45]), the model parameter associated with the isotropic term (motivated by the contribution of elastin; [22]) in the present model appropriately decreased with increasing age for AA specimens and decreased markedly for AAA specimens (tables 1–4)^{1}. We reported a similar finding when fitting carotid artery data from 9weekold male mgR/mgR mice compared with wildtype controls [11]. The mgR/mgR mice have 70 to 85 per cent less fibrillin1, a microfibrillar extracellular matrix glycoprotein thought to augment the structural stability of elastin over long periods of cyclic loading. It appears, therefore, that the neoHookean assumption for the elastindominated amorphous matrix is sufficient when modelling conduit arteries [40]. Finally, in parallel to the marked decrease in the bestfit value of c for AAA specimens, values of the parameters for the collagen families both increased and changed character relative to those for AA specimens: for specimens from young AAs, but for specimens from AAAs, thus capturing the observed biaxial reduction in extensibility/distensibility (figure 4). Relative changes in the values of these parameters may suggest a progressive change in the undulation of the collagen fibre families at physiological pressure, perhaps owing to diminished elastin which is originally highly prestretched owing to its predominant deposition during the perinatal period, its extreme elasticity and its long halflife. That is, our model suggests that elastin not only contributes to load bearing directly, it also influences the loadcarrying capability of the collagen (cf. [12]). Finally, possible active stresses were not accounted for because of the minimal smooth muscle in AAAs and the lack of information on activation in the autopsy samples of AA. Because of its potential importance in arterial mechanics [21], there is a need for information on possible changes in smooth muscle activity in human AA.
Because of the paucity of human data and the common concern of nonuniqueness in nonlinear parameter estimations, we combined the nonparametric bootstrap with a modal clustering method to identify, for the first time, secondorder accurate confidence intervals for bestfit values of the estimated model parameters. That is, albeit seldom done, determination of the precision of bestfit estimates should be integral to the process of parameter estimation when analysing arterial data from humans. We showed that local minima can exist in the nonlinear parameter estimations, hence use of randomly generated initial guesses and the nonparametric bootstrap can be very helpful in ensuring that a global minimum has been reached.
Whereas the present nonlinear regression and bootstrapping focused on estimating global bestfit values of model parameters suitable for performing stress analyses (tables 1–4), some of these values may be interpreted further given the underlying structural motivation. For example, consistent with [12], one could consider and , where ϕ^{α} ∈ [0,1] are, in principle, histologically measurable constituent mass fractions, ω^{i} ∈ [0,1] represent the percentage of collagen associated with each collagen family i (=1, 2, 3, 4) and are constituent specific material parameters. For example, He & Roach [45] reported that mass fractions for elastin, smooth muscle and collagen/GAGs were ϕ^{e} = 0.227, ϕ^{m} = 0.226 and ϕ^{c/g} = 0.548 in normal human AA, but ϕ^{e} = 0.024, ϕ^{m} = 0.022 and ϕ^{c/g} = 0.954 in AAAs (with the ratio of collagen to GAGs remaining approximately 2 : 1); Rizzo et al. [44] also reported a marked decrease in elastin (92% less) and increase in collagen (54% increase) in human AAAs, and LopezCandales et al. [46] reported an approximately 80 per cent depletion of smooth muscle cells in human AAAs that they speculated was due to increased apoptosis associated with the degradation of elastin. Although such histopathological findings can aid in a reinterpretation of bestfit values of the material parameters reported herein, there is also a need for further information. For example, we need to determine ‘effective mass fractions’ for the structurally sound (i.e. load bearing) constituents, not just the presence of that constituent. For example, it may be more appropriate to consider , where ξ^{e} ∈ [0,1] is a damage parameter that captures the loadcarrying capability of elastin, which may be present but fragmented; similar relations could be determined for the collagen. In other words, augmenting parameter estimation with histological information may allow even better structurally motivated modelling of macroscopic behaviours, thus emphasizing the need for better data. Of particular concern is better histological information on the oriented collagen fibres. For example, collagen I fibres may be stiffer than collagen III fibres, or fibres having larger diameters may be stiffer than those of the same type having smaller diameters. Until such information is available, it is probably better to assess fibre orientation via regressions of macroscopic mechanical data rather than simple histological analyses (e.g. smallangle light scattering); this caveat may address some of the concerns raised in Rodriguez et al. [36] and Haskett et al. [41]. Finally, having wellvalidated functional forms for the constitutive relations and realistic bounds on the values of the model parameters will facilitate future parameter estimations based on in vivo data alone (cf. [47]), for which patientspecific histological information will not be available.
Perhaps most importantly, however, we submit that the present approach to modelling is motivated by and consistent with the need for a new paradigm for stress analyses of AAAs. Notwithstanding many important findings from prior finiteelement stress analyses—including pivotal work by Fillinger et al. [8] showing that stresses can be better indicators of rupture potential than the traditionally used maximum diameter—all finiteelement models of AAAs are limited by the assumption of material homogeneity (e.g. [8,35–37,48–51]). Whereas this assumption stems from the lack of sufficient biaxial data for modelling regional variations in properties, limited mechanical data reveal regional variations [52–54]. Indeed, it is intuitive that if intramural cells seek to remodel the arterial wall in an attempt either to maintain stresses at or to restore stresses towards normal values (cf. [55]), the complex geometry of and applied haemodynamic loads acting on AAAs suggest that stresses (and thus remodelling) must vary from region to region (cf. [5,6]), consistent with measured variations in collagen content [56]. Indeed, computational work on the biomechanics of intracranial saccular aneurysms suggested years ago that regional variations are both likely and fundamental [4], and recent experimental findings on ascending thoracic aortic aneurysms similarly reveal regional variations [57]. There is, therefore, a pressing need to include regional variations in material properties within finiteelement analyses that seek to identify regions within AAAs that have the highest ratio of stress to strength and to correlate these findings with experimental data on rupture [3]. With regard to this general need, we submit that the present model can be useful in two regards. First, our model may facilitate correlations between measured histology and mechanical properties, which in turn could allow properties to be predicted regionally based primarily on histological data, which are arguably easier to obtain regionally given excised aorta. Second, our model can motivate appropriate growth and remodelling simulations (cf. [58]), which in turn can yield information on pointwise changes in mass fractions, orientations of collagen fibres and net multiaxial mechanical properties. That is, growth and remodelling codes, based on constitutive models validated using limited data from particular regions, may be able to predict regional variations in stresses and strength as a function of the evolution of a lesion [59]. Regardless, just as we have seen our understanding of AAAs deepen as we use increasingly better models (e.g. including nonlinear rather than linear material behaviour, anisotropic rather than isotropic properties, regional variations in thickness rather than a uniform thickness, tractionfree rather than a diastolic pressure loaded ‘stressfree’ reference condition and the presence of locally calcified regions or associated intraluminal thrombus), our understanding of the natural history of AAAs will increase significantly when we account properly for regional variations in material properties (e.g. resulting from differences in mass fractions, orientations and crosslinking of extracellular matrix), and the present model may provide an important step towards this important goal.
Acknowledgements
This work was supported, in part, by NIH grant HL86418 through the programme, Collaborations with National Centers for Biomedical Computing (SimBios at Stanford University). We are grateful in this regard for the support of Professor C. A. Taylor at Stanford. We also acknowledge Professors J. Vande Geest and M. Sacks who collected and reported the original biaxial data, Mr J. Muthu who manages this database at University of Pittsburgh and Ms M. Sand Enevoldsen for plotting some results useful for interpretation.
Footnotes

↵1 Because of the high potential for experimental errors (e.g. not orienting samples perfectly along circumferential and axial directions and possible bending of samples to perform inplane stretching tests) as well as specimentospecimen variations in testing human samples, general conclusions are best inferred from median or mean results, not single samples.
 Received June 2, 2010.
 Accepted June 30, 2010.
 This journal is © 2010 The Royal Society