## 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 ‘four-fibre 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 non-parametric bootstrap with a modal clustering method provides improved confidence intervals for estimated best-fit values of the eight associated constitutive parameters. It is suggested that this constitutive relation captures the well-known 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 (Choke *et al*. 2005; Shimizu *et al*. 2006; Vorp 2007). 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. Shah *et al*. 1998; Watton *et al*. 2004; Baek *et al*. 2006; Watton & Hill 2009). Computations also suggest that wall stress is a better predictor of rupture potential than the commonly used clinical metric of maximum diameter (Shah *et al*. 1998; Fillinger *et al*. 2003); 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 ‘four-fibre family’ constitutive relation that is motivated by microscopic data on the organization of arterial collagen (Wicker *et al.* 2008) and that captures differences in the mechanical behaviour of carotid arteries that result from diverse structural abnormalities in genetically modified mice (Gleason *et al.* 2008; Eberth *et al.* 2009), provides new insight into origins of residual stress and axial prestretch (Cardamone *et al.* 2009), and describes growth and remodelling in response to altered haemodynamics (Valentin & Humphrey 2009). Moreover, we show that combining the technique for statistical inference known as non-parametric bootstrap with a modal clustering method yields improved confidence intervals for best-fit values of the model parameters, particularly given the possible non-uniqueness in nonlinear estimations.

## 2. Material and Methods

### 2.1. Constitutive formulation

We re-analysed biaxial data on the mechanical behaviour of human AA and AAA reported by Vande Geest *et al*. (2004, 2006). 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 in-plane 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 elastin-dominated 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. Carlson & Secomb 2005) were absorbed naturally by the circumferential family of fibres. Hence, consider a strain energy function of the form (Baek *et al*. 2007)
2.1
where *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.2
where superscript *d* denotes ‘diagonal’. Diagonal families of collagen are typically regarded as mechanically equivalent, hence . This four-fibre 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. Wicker *et al*. 2008), 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.3
and
2.4
which enables a straightforward analysis of biaxial data.

Motivated by Gasser *et al*. (2006), 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.5
where *κ* is a scalar parameter representing fibre distributions in an integral sense. The associated Green-type strain experienced by the *i*th family of fibres, having mean orientation **M**^{i}, is characterized by *E*_{i} = **H**_{i}∶**C** − 1, namely
2.6
Thus, the second strain energy function considered herein has the form,
2.7
from which principal Cauchy stresses can be determined similar to equations (2.3) and (2.4). Note that mathematically admissible values are *κ* ∈ [0,1/2] (Li & Robertson 2009), where *κ* = 1/3 models an isotropic distribution of fibres and *κ* = 0 recovers the original ideal preferential alignment (Gasser *et al*. 2006). Hence, equation (2.7) is a natural extension of our four-fibre family model (equation 2.1), which is recovered as a special case. Whereas Gasser *et al*. (2006) used a single dispersion parameter *κ* to characterize distributions of two families of fibres, Pandolfi *et al*. (2008) 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, cross-linking, 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 best-fit values are physically realistic (Humphrey 2002). For example, the eight parameters associated with equation (2.1) must satisfy plus 0 ≤ *α*_{0} ≤ *π*/2 owing to symmetry of the diagonal fibres (Holzapfel *et al*. 2000). 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*= 2

*m*(i.e. biaxially) measured values of stretch within the central region, namely 2.8 with

*λ*_{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 (Ogden

*et al*. 2004), 2.9 with 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 built-in 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 pseudo-random 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}; Hurvick & Tsai 1989),
2.10
AIC_{c} converges to the AIC (Akaike 1974) as *n* becomes large, hence it can be employed regardless of sample size (Burnham & Anderson 2002); the preferred model is often the one with the lowest value of the AIC_{c}. This criterion for model selection was used previously by Zeinali-Davarani *et al*. (2009) to compare models with an increasing number of fibre families (from 2 to 10).

### 2.3. Non-parametric bootstrap

In addition to variations that arise in arterial properties owing to common adaptations, ageing or disease, inherent specimen-to-specimen variations exist even within gender- and age-matched samples (Humphrey 2002). 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 (Roy 1994). We used the non-parametric bootstrap method and a standard technique test both to estimate confidence intervals and to reveal possible non-uniqueness of the best-fit values of the parameters (figure 1). The bootstrap (Efron & Tibshirani 1993) 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 non-parametric 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 non-parametric 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 non-global minima are not eliminated from the parameter replications that are used to calculate confidence intervals, the resulting intervals will underestimate the precision of the best-fit parameters and negate advantages of non-parametric inference based on the bootstrap. To address this problem while keeping the bootstrap algorithm fully automatic and consistent with a non-parametric 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.11
where *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 (1981) 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 *p*-value 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 − *α* (Henderson *et al*. 2008).

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 (1996), the BC_{a} (bias-corrected and accelerated) method approximates confidence intervals for a parameter *θ* from empirical percentiles of the bootstrap empirical distribution . Compared with the standard percentile method (Efron & Tibshirani 1993), 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 one-sided level-*α* BC_{a} interval is defined by
2.12
where 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.13 and 2.14 where 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 (Efron & Tibshirani 1993). Using equation (2.12), with non-parametric 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 *second-order accurate* and it is *range-preserving*, 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 bias-correction constant *z*_{0} is estimated from the bootstrap distribution (Efron 1987). 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 (Efron & Tibshirani 1993; §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*. (2004, 2006) 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 four-fibre 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 eight-parameter model (mean *R*^{2} = 0.90 ± 0.10, AIC_{c} = 1010.6 ± 501.9) to the associated nine-parameter model obtained by adding a single dispersion parameter *κ* for all four-fibre 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 11-parameter 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), best-fit 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*. (2008) 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 in-plane fibre dispersion is not described by a distribution having rotational symmetry (cf. Gasser *et al*. 2006), which is consistent with the small values of *κ*_{i} that we found.

Allowing axial and circumferential families of fibres to be mechanically equivalent (cf. Wicker *et al*. 2008) lead to a simpler six-parameter 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 eight-parameter 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 cross-correlations 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 four-fibre family model gave an overall better description of data from both AA and AAA.

AA data were grouped based on age of the donors (Vande Geest *et al*. 2004), 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 eight-parameter four-fibre 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 neo-Hookean parameter *c* decreased with increasing age of subjects.

AAA data displayed even higher variability in both mechanical response and estimates of the best-fit 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 Vande Geest *et al*. 2006) 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 best-fit neo-Hookean 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 four-fibre 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 best-fit 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 non-negative neo-Hookean 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 collagen-related 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*. (2004) 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 three-dimensional (purely phenomenological) isotropic stored energy function of the form
4.1
where *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 two-dimensional (purely phenomenological) anisotropic Fung relation,
4.2
where *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*. (2006) 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 two-dimensional (purely phenomenological) anisotropic stored energy function, namely the Choi–Vito form
4.3
where *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 elastin-dominated amorphous matrix) and increase in the parameters and (associated with the predominant families of collagen fibres), each consistent with the well-known progressive stiffening of the aorta with age that associates with increased collagen to elastin ratios (Greenwald 2007). 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*. (2008) compared the ability of an isotropic (purely phenomenological) Demiray-type model and an anisotropic (structurally motivated, phenomenological) Gasser-type model to fit data from a representative AAA reported by Vande Geest *et al*. (2006). 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*. (2009) contrasted these findings with the ability of an isotropic polynomial stored energy function and a two-fibre 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 two-fibre 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*. (2009) fit the mean AAA data reported by Vande Geest *et al*. (2006) using a modified two-fibre family model (combining Demiray and polynomial functions for the isotropic contribution). They emphasized that the fit to data was much better with their eight-parameter anisotropic model than with a previously used isotropic polynomial model (Raghavan & Vorp 2000), but again the fit to data was not visually appealing (cf. their fig. 2) despite a reported *R*^{2} = 0.999. Basciano & Kleinstreuer (2009) proposed a three-dimensional (phenomenological) anisotropic stored energy function of the form
4.4
where *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 two-fibre directions similar to the original two-fibre family model (Holzapfel *et al*. 2000), and *α* and *β* are material parameters. Mean values of the best-fit 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. Watton *et al*. 2009) 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. Baek *et al*. 2006 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*. (in press) employed a variation of the two-fibre family model proposed by Gasser *et al*. (2006), 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 multi-photon microscopy (Wicker *et al*. 2008), by excellent prior fits to data from diverse arteries (Baek *et al*. 2007; Gleason *et al*. 2008; Eberth *et al*. 2009), by new interpretive value related to residual and pre-stresses (Cardamone *et al*. 2009) and by utility in growth and remodelling theories (Figueroa *et al*. 2009; Valentin & Humphrey 2009). 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. Wicker *et al*. 2008), 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 four-fibre family model (cf. Baek *et al*. 2007) over the original two-fibre family model (Holzapfel *et al*. 2000) 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 two-fibre family model less capable of capturing the observed marked axial stiffness of AAA specimens (cf. figure 4). The four-fibre family model is similarly expected to be free of limitations of the two-fibre family model identified by Gasser *et al*. (2006) 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 (Greenwald 2007; O'Rourke & Hashimoto 2007) and AAAs (e.g. Rizzo *et al*. 1989; He & Roach 1994), the model parameter associated with the isotropic term (motivated by the contribution of elastin; Holzapfel *et al*. 2000) 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 9-week-old male mgR/mgR mice compared with wild-type controls (Eberth *et al*. 2009). The mgR/mgR mice have 70 to 85 per cent less fibrillin-1, a microfibrillar extracellular matrix glycoprotein thought to augment the structural stability of elastin over long periods of cyclic loading. It appears, therefore, that the neo-Hookean assumption for the elastin-dominated amorphous matrix is sufficient when modelling conduit arteries (Watton *et al*. 2009). Finally, in parallel to the marked decrease in the best-fit 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 half-life. That is, our model suggests that elastin not only contributes to load bearing directly, it also influences the load-carrying capability of the collagen (cf. Cardamone *et al*. 2009). 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 (Humphrey 2002), 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 non-uniqueness in nonlinear parameter estimations, we combined the non-parametric bootstrap with a modal clustering method to identify, for the first time, second-order accurate confidence intervals for best-fit values of the estimated model parameters. That is, albeit seldom done, determination of the precision of best-fit 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 non-parametric 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 best-fit 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 Cardamone *et al*. (2009), 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 (1994) 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*. (1989) also reported a marked decrease in elastin (92% less) and increase in collagen (54% increase) in human AAAs, and Lopez-Candales *et al*. (1997) 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 best-fit 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 load-carrying 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. small-angle light scattering); this caveat may address some of the concerns raised in Rodriguez *et al*. (2009) and Haskett *et al*. (in press). Finally, having well-validated 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. Stalhand & Klarbring 2005), for which patient-specific 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 finite-element stress analyses—including pivotal work by Fillinger *et al*. (2003) showing that stresses can be better indicators of rupture potential than the traditionally used maximum diameter—all finite-element models of AAAs are limited by the assumption of material homogeneity (e.g. Raghavan *et al*. 2000; Wang *et al*. 2002; Fillinger *et al*. 2003; Lu *et al*. 2007; Speelman *et al*. 2007; Rodriguez *et al*. 2008, 2009; Rissland *et al*. 2009). Whereas this assumption stems from the lack of sufficient biaxial data for modelling regional variations in properties, limited mechanical data reveal regional variations (Thubrikar *et al*. 2001; Vallabhaneni *et al*. 2004; Raghavan *et al*. 2006). 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. Humphrey 2008), the complex geometry of and applied haemodynamic loads acting on AAAs suggest that stresses (and thus remodelling) must vary from region to region (cf. Watton *et al*. 2004; Baek *et al*. 2006), consistent with measured variations in collagen content (Menashi *et al*. 1987). Indeed, computational work on the biomechanics of intracranial saccular aneurysms suggested years ago that regional variations are both likely and fundamental (Shah *et al*. 1998), and recent experimental findings on ascending thoracic aortic aneurysms similarly reveal regional variations (Iliopoulos *et al*. 2009). There is, therefore, a pressing need to include regional variations in material properties within finite-element 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 (Vorp 2007). 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. Humphrey & Taylor 2008), which in turn can yield information on point-wise 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 (Brady *et al*. 2004). 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, traction-free rather than a diastolic pressure loaded ‘stress-free’ 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 cross-linking 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 HL-86418 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 in-plane stretching tests) as well as specimen-to-specimen 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.

- © 2010 The Royal Society