IDENTIFICATION OF MATERIAL AND GEOMETRIC PARAMETERS OF
ARTERIAL WALL FOR PATIENT-SPECIFIC VASCULAR GROWTH AND
REMODELING MODELS
By
Shahrokh Zeinali-Davarani

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2011

ABSTRACT
IDENTIFICATION OF MATERIAL AND GEOMETRIC PARAMETERS OF
ARTERIAL WALL FOR PATIENT-SPECIFIC VASCULAR GROWTH AND
REMODELING MODELS
By
Shahrokh Zeinali-Davarani

The uncertainty associated with an individual abdominal aortic aneurysm (AAA) rupture
carries a considerable amount of health risks as well as social and economic burden. There is
a need for advanced technologies to mitigate this burden by early detection, patient-specific
risk assessment and clinical management. Computational vascular mechanics has been of
great interest along with the recent advances in medical imaging, experimental technique
and computational simulations. In particular, computational modeling of vascular growth
and remodeling (G&R) has provided further understanding of the G&R process that occurs
in vascular diseases. Despite rapid expansion of our knowledge of vascular G&R, developing
patient-specific models of AAA evolution is still an open problem and subject to numerous
challenges.
In this study, a framework is presented where, first the identification of the intrinsic material and geometric parameters for patient-specific modeling is addressed and then a finite
element model of AAA G&R that accounts for medical image-based geometrical models is
introduced. With regard to the material parameters, estimation of parameters of a multifiber family model for passive arteries in the presence of the measurement noise is considered.
First, the uncertainty propagation due to the errors in variables is carefully characterized
using the constitutive model. Then, the parameter estimation of the artery model is formulated into nonlinear least squares optimization with an appropriately chosen weight from
the uncertainty model. The proposed technique is evaluated using multiple sets of synthesized data with fictitious measurement noises. The results of the estimation are compared
with those of the conventional nonlinear least squares optimization without a proper weight

factor. The proposed method significantly improves the quality of parameter estimation as
the amplitude of the errors in variables becomes larger. We also investigate model selection
criteria to decide the optimal number of fiber families in the multi-fiber family model with
respect to the experimental data. The effect of multiple models for mechanical behavior of
arteries is also investigated.
Distribution of the geometric parameters, being wall thickness and fiber orientations, are
estimated using a 2-D parameterization of the vessel wall surface and a global approximation
scheme integrated within an inverse optimization routine. Two conditions determine the
objective of the optimization. First, the fundamental assumption in adaptation models,
namely the existence of vascular homeostasis in normal vessels, should be maintained in
a vessel model built from medical images. Second, the deviation of vessel wall from the
original/image geometry subject to the normal pressure should be minimized. The same
inverse technique is used to investigate the consequence of different homeostasis assumptions
on the optimized distribution of parameters. The numerical optimization method results in
a considerable improvement in both satisfying homeostatic condition and minimizing the
deviation of geometry from the original shape based on in-vivo images.
Then, a finite element model of stress-mediated G&R of arteries based on medical imagebased geometries is presented for simulation of AAAs. Degradation of elastin initiates a local
dilatation of the aorta while stress-mediated turnover of collagen and smooth muscle compensates the loss of elastin. Stress distributions and expansion rates during the aneurysm growth
are studied for multiple spatial distributions of elastin degradation and kinetic parameters
of the growth. Temporal variations of the elastin degradation and kinetic parameters are
also investigated with either direct time-dependent degradation or stretch-induced degradation as possible biochemical and biomechanical mechanisms for elastin degradation. The
results show that this computational model has the capability to capture the complexities of
aneurysm progression due to variations of geometry, extent of damage, and stress-mediated
turnover as an important step towards patient-specific modeling.

To my parents whose endless patience, unwavering love and selflessness have been
inspirational to me.

iv

ACKNOWLEDGEMENTS

I would like to thank Dr. Seungik Baek for his support, expertise and dedication throughout my research at the Cardiovascular Tissue Mechanics Laboratory. His passion in pursuing
solutions to the challenging problems in Biomechanics has been a motivation for my research
career. I am grateful to Dr. Guy Raguin for his support and advice that have helped me
improve this work. Also, I would like to thank Drs. Alejandro Diaz, Thomas Pence and
Jung-Wuk Hong for serving on my committee and for their insightful suggestions. I am also
grateful to other colleagues from whom I have benefited over the years; Dr. Jungsil Kim,
Chris Hunley, Azadeh Sheidaei, Byron Zambran, Gui DeAraujo and Alex Dupay.
I would like to sincerely acknowledge all financial supports that I have received from
the Department of Mechanical Engineering, the College of Engineering and the Graduate
School at Michigan State University. It could have been a cumbersome task to conduct and
develop this research without the support of the Michigan State University High Performance
Computing Center and the Institute for Cyber Enabled Research.

v

Table of Contents
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv
1 INTRODUCTION . . . . . . . . . . . . . .
1.1 Anatomy, structure and function of the aorta .
1.2 Mechanical behavior of the aortic tissue . . . .
1.3 Abdominal Aortic Aneurysm . . . . . . . . . .
1.4 Computational analysis of AAA . . . . . . . .
1.5 Computational modeling of vascular growth .
1.6 Motivations . . . . . . . . . . . . . . . . . . .
1.7 Objectives . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

1
1
2
4
5
5
6
7

2 ESTIMATION OF MATERIAL PARAMETERS OF AR-

TERIAL WALL TISSUE USING A WEIGHTED NONLINEAR LEAST SQUARES METHOD . . . . . . . . . . . . . . . . . .
2.1
2.2

2.3

2.4

8
INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
METHODS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Constitutive relations for the mechanical behavior of the passive artery 10
2.2.2 WNLS optimization for a biaxial test of the artery . . . . . . . . . . 12
2.2.3 Model selection criteria for optimal number of parameters . . . . . . 14
RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.1 Comparison between the proposed WNLS optimization and the conventional NLS optimization . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.2 Optimal number of fiber families . . . . . . . . . . . . . . . . . . . . 20
DISCUSSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3 COMPARISON OF DIFFERENT CONSTITUTIVE MOD-

ELS OF VASCULAR TISSUE MECHANICAL BEHAVIOR
FROM PARAMETER ESTIMATION PERSPECTIVE . . .
3.1
3.2
3.3

3.4

Constitutive modeling of arterial wall
Holzapfel-type model . . . . . . . . .
Models for elastin contribution . . . .
3.3.1 Zulliger model . . . . . . . . .
3.3.2 Demiray model . . . . . . . .
Models for collagen contribution . . .
vi

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

25
25
26
28
28
28
30

3.5

3.4.1 Zulliger model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Lanir model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30
32
32

4 ESTIMATION OF THE ARTERIAL WALL THICKNESS

AND ANISOTROPY FOR IMAGE-BASED MODELS OF
ARTERIES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1
4.2

4.3
4.4

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . .
METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Estimation of material constitutive parameters . . . . . . . . .
4.2.2 Inverse optimization problem statement . . . . . . . . . . . . .
4.2.3 Parameterizing the aortic wall surface with two variables . . .
4.2.4 Global approximation approach . . . . . . . . . . . . . . . . .
4.2.5 Optimization algorithm . . . . . . . . . . . . . . . . . . . . . .
RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Finding the optimal distributions of thickness and anisotropy .
DISCUSSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

34
34
35
35
40
43
46
46
49
49
53

5 TESTING DIFFERENT HYPOTHESES OF VASCULAR

HOMEOSTASIS USING IMAGE-BASED GEOMETRIC MODELS AND AN INVERSE OPTIMIZATION METHOD . . . 57
5.1
5.2
5.3
5.4
5.5

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . .
Hypotheses on vascular mechanical homeostasis . . . . . .
Constitutive relations and an inverse optimization method
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

57
59
61
62
67

6 A FINITE ELEMENT MODEL OF STRESS-MEDIATED

VASCULAR ADAPTATION: APPLICATION TO ABDOMINAL AORTIC ANEURYSMS . . . . . . . . . . . . . . . . . . . . . . 71
6.1
6.2

6.3

6.4

6.5

Introduction . . . . . . . . . . . . . . . . . . . . . .
Theoretical modeling of arterial G&R . . . . . . . .
6.2.1 Kinematics and important configurations . .
6.2.2 Stress response and constitutive assumptions
6.2.3 Stress-mediated G&R . . . . . . . . . . . . .
Computational considerations . . . . . . . . . . . .
6.3.1 Local Cartesian coordinate system . . . . .
6.3.2 Finite element formulation . . . . . . . . . .
6.3.3 Numerical solutions . . . . . . . . . . . . . .
Simulation of an AAA growth . . . . . . . . . . . .
6.4.1 A geometric model and mesh generation . .
6.4.2 Initialization for G&R simulation . . . . . .
6.4.3 Simulation Cases . . . . . . . . . . . . . . .
6.4.4 Results . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . .
vii

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

71
73
73
76
79
80
80
82
87
88
88
88
91
92
97

7 MEDICAL IMAGE-BASED SIMULATION OF ABDOMINAL AORTIC ANEURYSM GROWTH . . . . . . . . . . . . . . .
7.1 A framework for modeling of AAA growth based on medical images . . . . .
7.2 Simulation cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

107
107
110
111
113

8 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
A Equilibrium equations considered for inflation-extension test
of tubular vessels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
B Nelder-Mead Simplex
B.1 Introduction . . . . . . .
B.2 Algorithm statement . .
B.3 Stopping criteria . . . .

Algorithm . . . . . . . . . . . . . . . . . . . . 130
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

viii

List of Tables
2.1

True and estimated parameters for the fictitious data with moderate and high
slopes (corresponding to Figures 2.2 and 2.3) and the experimental data from
mouse carotid and rabbit basilar arteries (corresponding to Figures 2.4(a) and
2.4(b)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

Averaged normalized error e at different noise levels using the NLS and WNLS
methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

True and estimated linearized stiffness at Pi = 75mmHg and λz = 1.3 corresponding to Figures 2.2 and 2.3 using the NLS and WNLS methods. . . . . .

21

4.1

Summary of material parameters used in the optimization. . . . . . . . . . .

39

6.1

Constitutive and kinetic parameters for each constituents used in initialization
and G&R simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

2.2

2.3

ix

List of Figures
1.1

1.2

2.1

2.2

2.3

2.4

Cross-section of the aorta, stained for smooth muscle cells using H&E (a)
stained for elastin fibers using VVG(b), stained for collagen fibers using picosirius red in light microscopy (c) or observed under polarized microscope
(d). (For interpretation of the references to color in this and all other figures,
the reader is referred to the electronic version of this dissertation). . . . . . .

2

Typical stress-strain curves of a circumferential strip of an artery during the
uniaxial test (Holzapfel et al., 2000). . . . . . . . . . . . . . . . . . . . . . .

3

Pressure versus radius plots for the 4-fiber-family model for a set of true data
(solid lines) and fictitious noisy data (dotted lines) at noise levels of 0.005 (a)
and 0.01 (b) and different axial stretches λz = 1.2, 1.3, 1.4. . . . . . . . . . .

17

Comparison between the NLS (a and c) and WNLS (b and d) optimizations
for the 4-fiber-family model. Pressure and axial force versus radius are plotted
using true parameters (dotted lines) and estimated parameters (solid lines) for
the data with moderate slopes at different axial stretches λz = 1.2, 1.3, 1.4 .

18

Comparison between the NLS (a and c) and WNLS (b and d) optimizations
for the 4-fiber-family model. Pressure and axial force versus radius are plotted
using true parameters (dotted lines) and estimated parameters (solid lines) for
the data with stiff slopes in the high stretch region at different axial stretches
λz = 1.2, 1.3, 1.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

Pressure-radius-axial stretch curves are plotted from fitting of the experimental data for the mouse carotid artery (a) and rabbit basilar artery (b) using
the WNLS optimization. Circles show data points and solid lines show the
estimated values. AICc and the residual are plotted against the number of
parameters; 4, 6, 8, 11, 14, 17 (or the number of fiber families; 2, 3, 4, 6, 8,
10) for the mouse carotid (c) and rabbit basilar (d) arteries using the WNLS
optimization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

x

3.1

3.2

3.3

3.4

3.5

4.1

4.2

4.3

4.4

4.5

Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries using Holzapfel-type model. Dots show the data points and solid
lines show the estimated values. . . . . . . . . . . . . . . . . . . . . . . . . .

27

Pressure-radius data of mouse carotid artery and the estimated values (left)
with the corresponding normalized residuals (right) using Zulliger model. Dots
show the data points and solid lines show the estimated values. . . . . . . . .

28

Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries using Demiray model. Dots show the data points and solid lines
show the estimated values. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries using Zulliger model. Dots show the data points and solid lines
show the estimated values. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries using Lanir model. Dots show the data points and solid lines
show the estimated values. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

Kinematics of the deformation associated with biaxial mechanical test and
the corresponding deformation gradients. λ1 and λ2 are stretches in circumferential and axial directions during the biaxial test . . . . . . . . . . . . . .

36

Stress versus stretch plots in both directions. Biaxial test data (empty markers) and fitted values (filled markers) using the estimated parameters. Different markers correspond to a different ratios of tension applied simultaneously
in both directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

Geometry of an arbitrary model of the arterial wall with its centerline; Approximation/nodal points with their associated length s = Lj (j = 1, ..., J).
a is an arbitrary vector used in order to find the orientation θ associated with
a point (Xc ) on the wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

Geometry of the vessel wall parameterized with longitudinal and azimuthal
(s and θ) variables. Dots represent center points of all elements on the wall .

48

The effect of variation of the parameter ξ on both GD and SD using (M =
3, N = 3;top) and (M = 6, N = 5;bottom) . . . . . . . . . . . . . . . . . . .

50

xi

4.6

Changes in the objective function and its associated compartments versus optimization iterations using 36 variables (18 variables for approximating thickness and 18 variables for approximating fiber orientation) considering ξ = 0.01 51

4.7

Deviation of the geometry from the in vivo geometry without (a) and with
(b) optimized distributions of thickness and anisotropy (||x − Ximage ||) . . .

52

k
k
Deviation of the stress ((σ k − σh )/σh ) from the target homeostatic stress in
a helical fiber (k = 3) without (a) and with (b) optimized distributions of
thickness and anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

k
k
Deviation of the stress ((σ k − σh )/σh ) from the target homeostatic stress in
a helical fiber (k = 4) without (a) and with (b) optimized distributions of
thickness and anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

4.10 Distributions of thickness (a) and anisotropy (b) obtained from the optimization results using ξ = 0.01, M = 6, and N = 3 . . . . . . . . . . . . . . . . .

54

4.8

4.9

5.1

Change in the objective function and its associated compartments versus optimization iterations in Cases 1 and 2 using 60 variables (M = 6;N = 5;ξ = 0.01). 62

5.2

Optimized distributions of thickness (a, b) and anisotropy (c, d) for the aorta
model obtained in Cases 1 and 2 (M = 6;N = 5;ξ = 0.01). Black arrows
identify the direction of one set of corresponding collagen fibers. . . . . . . .

64

Optimized distributions of thickness (a, b) and anisotropy (c, d) for the internal iliac artery model obtained in Cases 1 and 2 (M = 5;N = 3;ξ = 0.1). .

65

Deviation of the geometry from the in vivo configuration (||x − Ximage ||)
and deviation of stress in a helical fiber direction from the homeostatic value
k
k
((σ k − σh )/σh ) when (M = 5;N = 3) (a, c) and when (M = 6, N = 5) (b, d)
for the iliac artery model. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Deviation of the geometry (a), cyclic stretch in circumferential direction (λc )(b)
and corresponding distributions of thickness (c) and anisotropy (d) for the
aorta model obtained in Case 3 (M = 6;N = 3;ξ = 0.01). Black arrows
identify the direction of one set of corresponding collagen fibers. . . . . . . .

68

5.3

5.4

5.5

6.1

Schematic view of configurations involved in the G&R simulation. A fixed reference configuration is considered for the computation of deformation associated with G&R time τ = [0, s]. It is chosen to coincide with the configuration
of a healthy in vivo artery at time τ = 0 (i.e., κ0 ). We imagine the existence
of a stress-free configuration, κsf , associated with the current configuration κs . 75

xii

6.2

Orthonormal base vectors for a global Cartesian coordinate system are {E1 , E2 , E3 }.
Those for a local Cartesian coordinate system are Ee , Ee (associated with ξ1 ,
1
2
ξ2 ), and Ee = Ne at a point set to be the center point of a linear triangular
3
element tangent to the membrane at Xc . Ee and Ee are orthonormal bases
1
2
on the tangent plane and Ne is the outward unit normal vector. . . . . . . . 81

6.3

c
c
Deviation of the intramural stress from the homeostatic value ((σ k − σh )/σh )
in axial and two helical fiber directions before (a, b, and c) and after (d, e,
and f) 300 days of the initial simulation. . . . . . . . . . . . . . . . . . . . .

90

Magnitude of the displacement (mm) calculated from the reference configuration after 300 days of the initial simulation (a) and the corresponding geometry
before and after 300 days (b). . . . . . . . . . . . . . . . . . . . . . . . . . .

91

Areal mass density of elastin (kg/m2 ) for different simulation Cases (1)-(4).
Cases (1) and (2) correspond to different damages shapes distributed to relatively large area which are applied at different locations (concave and convex
sides) and Cases (3) and (4) correspond to more localized damage shapes, but
at multiple locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Maximum principal stress distribution (kP a) for the simulation Case (3) at
50, 1200, and 2700 days after the initial damage (a, b, and c) and the corresponding areal mass density of collagen (d, e, and f) (kg/m2 ). . . . . . . . .

95

Maximum principal stress (kP a) for simulation Cases (1), (2), and (4) after
50 (a, b, and c) and 2700 (d, e, and f) days. . . . . . . . . . . . . . . . . . .

96

Expansion rates of the lesion in simulation Cases (1)-(4) as the maximum
diameter of the lesion versus time and their best linear fit (Kg = 0.05). . . .

97

Expansion rates of the lesion in simulation Case (2) with different values of
the kinetic parameter (Kg = 0.02, 0.05 and 0.1). . . . . . . . . . . . . . . . .

98

6.10 Areal mass density of elastin (kg/m2 ) for the case of ‘time-dependent’ degradation after 50 and 2700 days (a, b) and the corresponding maximum principal
stress distribution (c, d) (kP a). . . . . . . . . . . . . . . . . . . . . . . . . .

99

6.4

6.5

6.6

6.7

6.8

6.9

6.11 Areal mass density of elastin (kg/m2 ) for the case of ‘stretch-induced’ degradation after 50 and 2700 days (a, b) and the corresponding maximum principal
stress distribution (c, d) (kP a). . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.12 Expansion rates of the lesion in the ‘time dependent’ simulation case associated with different time constants (T = 3, 6, 10, and 15 years). . . . . . . . . 102

xiii

7.1

Schematic diagram of the steps taken for simulating AAAs using a medical
image-based geometry and material parameters of a healthy aorta as the initial
state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

7.2

Variation of damage ratio with respect to time (Rt vs. time) and variation of
the sensitivity parameter with respect to the local damage (Kg vs. R(l, θ, t)).
Solid line represents the time dependent damage ratio and dashed line represents the damage dependent variation of Kg . . . . . . . . . . . . . . . . . . . 110

7.3

left: Areal mass distributions (kg/m2 ) for elastin after introducing damage.
right: the resultant distributions of von Mises stress (kP a) due to G&R. The
arrows indicate the regions of elevated stress during AAA expansion. . . . . 112

7.4

Variations of peak von Mises stress (a) and maximum diameter (b) with G&R
time for different cases. Equations for best fit lines are shown (lines are not
included in the graph). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

7.5

Spatial relations between von Mises stress and both the rate of collagen mass
production as well as the rate of wall displacement for Cases 1 at different
stages of growth. Wall thickness distribution and collagen concentration are
illustrated for each case corresponding to 2700 days of G&R. . . . . . . . . . 115

7.6

Spatial relations between von Mises stress and both the rate of collagen mass
production as well as the rate of wall displacement for Cases 5 at different
stages of growth. Wall thickness distribution and collagen concentration are
illustrated for each case corresponding to 2700 days of G&R. . . . . . . . . . 116

7.7

Spatial relation between the rate of collagen mass production and von Mises
stress for Case 7 using a homogeneous distribution of kinetic parameters (a;
Kg = 0.05) as well as spatially variable kinetic parameter (b; Kg = 0.02 at
the AAA sac and Kg = 0.05 otherwise) at different times. . . . . . . . . . . . 117

7.8

Elastin (a) and collagen mass (b) areal density along with von Mises stress
(c) distributions for Case 8 after 7 years (2700 days) of G&R. Time history
of maximum diameter and peak von Mises stress is also shown (d). . . . . . 118

7.9

Spatial relation between the nodal displacement and the damage ratio for
Cases 1 and 3 at different times. A number of points on the convex and
concave regions (with reference to the healthy geometry) are distinguished
with different markers at 2700 days. . . . . . . . . . . . . . . . . . . . . . . . 120

7.10 Distributions of axial stretch for different aneurysms. . . . . . . . . . . . . . 122
8.1

A comprehensive simulation framework aimed at patient-specific modeling of
AAA G&R. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

xiv

Chapter 1
INTRODUCTION

1.1

Anatomy, structure and function of the aorta

The aorta is the largest artery that supplies oxygenated blood from the heart to other organs
and tissues. It first extends from the heart in upward direction (called ascending aorta)
towards the aortic bend (called aortic arch) and then turns downward into the body towards
the diaphragm (called descending thoracic aorta). The portion of aorta below the diaphragm
(descending abdominal aorta), located in front of spine, extends downward until it bifurcates
into common iliac arteries supplying blood to legs. The aortic wall is comprises three distinct
layers of tunica intima, tunica media, and tunica adventitia. The intima is the innermost
layer of an aorta and consists of a layer of endothelial cells and a subendothelial layer of
connective tissue. The internal elastic lamina separates the intima and the media. The
media is the middle layer consisting if concentrically arranged smooth muscle cells, collagen
fibers and elastin. The external elastic lamina separates the media and the adventitia.
Finally, the outermost layer of an aorta, the adventitia, is primarily composed of fibroblasts,
elastin and longitudinally arranged collagen fibers. The primary structural proteins that
determine the mechanical properties of aortic wall are collagen fibers and elastin, contents
of which vary with age and physiological and pathological conditions and the anatomical
location throughout the body.

1

(a)
(b)
(c)
(d)
Figure 1.1: Cross-section of the aorta, stained for smooth muscle cells using H&E (a) stained
for elastin fibers using VVG(b), stained for collagen fibers using picosirius red in light microscopy (c) or observed under polarized microscope (d). (For interpretation of the references
to color in this and all other figures, the reader is referred to the electronic version of this
dissertation).
A histological image of the cross-section of an aorta with different stains revealing different
constituents is shown in Fig. 1.1. The cell nuclei are blue in (a), elastin is shown in black (b)
while the collagen is red (c). In order to examine the orientation of collagen fibers, polarized
light microscopy can be used (d).

1.2

Mechanical behavior of the aortic tissue

Like many other biological soft tissues, an aorta undergoes relatively large deformations and
its behavior is nonlinear, inelastic, and anisotropic over finite strains. Figure 1.2 shows a
typical stress-strain response of arterial tissue during the uniaxial test. During the first few
loading cycles, the arterial tissue displays stress softening (Fig. 1.2). It then exhibits a nearly
repeatable behavior (Fung et al., 1979), due to which the mechanical behavior of arterial
tissues is commonly described by pseudo-elasticity. Also, as mentioned earlier, an artery is
considered a layered mixture-composite of elastin, collagen fibers and smooth muscle cells
as well as significant amount of water (Humphrey, 2002). Incompressibility has been found
to be reasonable assumption in arterial mechanics (Chuong and Fung, 1984).

2

Elastic Inelastic

Physiological
range

STRESS

Material stress-strain response
Engineering response

STRAIN

Figure 1.2: Typical stress-strain curves of a circumferential strip of an artery during the
uniaxial test (Holzapfel et al., 2000).

3

1.3

Abdominal Aortic Aneurysm

Abdominal aortic aneurysm (AAA) is the focal dilation of the abdominal aorta that enlarges over the years and eventually rupture which is considered as one the leading causes
of mortality. AAAs affect 2 million people in the US alone and are associated with some
environmental risk factors such as smoking, hypertension and genetic predisposition (Powell,
2002). Despite the abundant clinical studies, the etiology and the natural history of AAAs
remains poorly understood. Aortic aneurysm may reach to 3.5-13 cm in diameter with an
expansion rate of 0.1-1.5 cm/year (Humphrey, 2002). Aneurysmal wall is commonly devoid
of elastin and smooth muscle with a higher ratio of mass for collagen (He and Roach, 1994),
thus giving rise to a fibrotic wall with a stiffer structure. Disruption of elastin, either as a
cause or as a consequence of the disease process, is associated with aneurysm dilation, while
collagen failure was found necessary for rupture (Dobrin et al., 1984).
Current treatments of AAA include either open repair which involves invasive surgery
with extensive tissue dissections or minimally invasive endovascular repair (Vorp and Vande
Geest, 2005). Both treatments are costly and associated with significant amount of risks.
The current clinical decision making for surgical treatments is based on “maximum diameter
criterion”. That is, the likelihood of AAA rupture is assumed to be higher than the risk
of surgical treatments when the aneurysm reaches 5-5.5 cm in diameter (Vorp, 2007). Yet
it is not clear when an individual AAA may rupture since a considerable amount of AAAs
may rupture at smaller sizes. As the population of elderly people grows, the social and
economic burden, imposed by the uncertainty of AAA rupture, will increase. In order to
reduce this public health burden, there are crucial needs for advanced technologies that can
provide AAA patients with early detection, patient-specific risk assessment, and safe clinical
interventions.

4

1.4

Computational analysis of AAA

Mechanically rupture of the tissue occurs when the local stress exceeds the tissue strength.
In that sense, computational stress analysis based on 3D Computer Tomography (CT) or
Magnetic Resonance Imaging (MRI) using nonlinear constitutive models of the vessel enabled
researchers to estimate wall stress more accurately (Dorfmann et al., 2010; Fillinger et al.,
2002; Raghavan et al., 2000; Rissland et al., 2009; Speelman et al., 2007) and, hence, led to
a better prediction of AAA rupture risk than the maximum diameter criterion.
However estimation of the stress alone may not provide a reliable estimation of rupture
potential since the rupture potential depends also on the strength (Vorp and Vande Geest,
2005). Strength of AAA wall tissue depends on its integral microstructure which is constantly
evolving in an abnormal fashion. Furthermore, finite element (FE) analysis based on a
patient-specific geometry yields the stress distribution only for a fixed geometry and does
not account for the time evolution of AAA shape and the wall mechanics.

1.5

Computational modeling of vascular growth

During normal or disease conditions (e.g. aortic aneurysms), arterial wall undergoes continuous changes in terms of morphology, histology and biomechanics in an adaptive manner.
Vascular tissue exhibits a remarkable ability to adapt in various physiological and pathological conditions, often thought to be governed by mechanical factors (Driss et al., 1997;
Jackson et al., 2005; Mulvany, 1992). Understanding adaptive processes involving aortic
tissue growth (i.e. change in mass) and remodeling (i.e. change in structure) is instrumental
in developing models capable of simulating AAA evolution.
For the past decade, computational modeling of vascular growth and remodeling (G&R)
have been developed and used to test multiple hypotheses based on experimental and clinical
studies (Watton et al., 2004; Menzel, 2005; Baek et al., 2006; Hariton et al., 2007b; Kroon
and Holzapfel, 2007; Kuhl et al., 2007; Figueroa et al., 2009; Kroon and Holzapfel, 2009a;
5

Watton and Hill, 2009). Mechanical stress/strain is the parameter driving adaptation process
in all models. Mechanical homeostasis is the fundamental assumption of most G&R models,
meaning there is a tendency toward a preferred level of mechanical stress/strain throughout
the arterial tissue. Despite being focused on symmetric or idealized geometries, those models
revealed a great potential in understanding basic features of vascular G&R. However, in-vivo
aortic aneurysms normally evolve into irregular 3-D shapes and idealized models of G&R
are not sophisticated enough to capture patient-specific features of an in-vivo AAA.

1.6

Motivations

Ultimately, for clinical applications, the current stage of modeling needs to be improved
to utilize patient-specific geometries of AAAs. Along with the recent advances in medical imaging, this can revolutionize the state-of-the-art analyses and clinical interventions.
Patient-specific modeling raises new challenges with regard to characterization of patientspecific parameters of the vessel wall or aneurysmal tissue under certain conditions.
With regard to material parameters intrinsic to the material, mechanical experiments
have been widely conducted to develop appropriate constitutive forms for the mechanical
behavior of arteries and to estimate the associated constitutive parameters. Meanwhile, less
attention has been paid to the fact that measurements are not error-free and precise and
unbiased parameter estimation in vascular mechanics requires one to consider the uncertainty in all measured variables when estimating parameters. Furthermore, identification of
geometric parameters, e.g. wall thickness and material anisotropy, becomes challenging for
patient-specific models when certain conditions are to be satisfied. First condition to be met
is the mechanical homeostatic state as the main assumption for G&R models. Second condition stems from the fact that the stress-free configuration is not available when a medical
image-based geometry is used. That is, during the computation, the geometry no longer
coincides with the original or image geometry under the physiological loads.

6

1.7

Objectives

The ultimate goal of this study is to develop patient-specific models of AAA G&R. Toward
that goal, multiple specific aims are set in this study based on the existing challenges as:
• Develop an enhanced least squares method to account for the measurement noise in
estimation of material constitutive parameters.
• Develop an optimization scheme needed to estimate arterial wall thickness and anisotropy
such that the patient-specific conditions are satisfied.
• Develop a 3-D finite element framework capable of modeling AAAs growth using
patient-specific geometries.
In chapter 2 a weighted nonlinear least squares method that incorporates a noise model is
proposed and its advantage over the conventional nonlinear least squares is demonstrated. In
chapter 3, different constitutive relations are tested for parameter estimation of mechanical
behavior of mouse carotid and rabbit basilar arteries. In chapter 4 an inverse optimization
technique is developed to estimate the distributions of wall thickness and material anisotropy
such that the vascular homeostatic as well as the geometric conditions are satisfied simultaneously. Using the proposed technique in chapter 4, different hypotheses for mechanical
homeostasis are examined in chapter 5. In chapter 6 an image-based finite element model
for vascular adaptation is presented and its application to AAA is demonstrated in sample
simulation cases. Chapter 7 presents a framework for simulation of AAA, demonstrates
the potential utility of the framework with additional simulation cases, and discusses its
shortcomings and future developments.

7

Chapter 2
ESTIMATION OF MATERIAL PARAMETERS OF ARTERIAL
WALL TISSUE USING A WEIGHTED NONLINEAR LEAST
SQUARES METHOD

2.1

INTRODUCTION

Mechanical properties of the arterial wall vary with anatomical variations and for different
species. It has also been suggested that the structure and material properties of the arterial
wall alter under various physiological or pathological conditions (Langille, 1993; Humphrey,
2008). In such studies, precise parameter estimation is essential to quantify the difference
in mechanical behavior. In vascular mechanics, mechanical behavior of the arterial wall is
typically described by a nonlinear equation,

y = f (x, Θ),

(2.1)

where x and y are vectors of independent and dependent variables and Θ is a set of unknown
parameters. Experimental studies in vascular mechanics normally utilize a nonlinear least
squares (NLS) technique in which the sum of the squares of the difference between the experimental measurements and the calculated responses of dependent variables y (e.g., pressure
and axial force) is minimized, while independent variables x (e.g., outer diameter and axial

8

stretch) are considered free of error (Pandit et al., 2005; Schulze-Bauer and Holzapfel, 2003;
Wang et al., 2006). All variables, however, are measured with errors and it is known that an
NLS method results in biased parameter estimation when uncertainty exists in all variables
in a constitutive model (Emery et al., 2000; Fadale et al., 1995). In such a case, parameter
estimation should be correctly formulated as the nonlinear least squares optimization with
an appropriately chosen objective function (Fadale et al., 1995; Schwetlick and Tiller, 1985).
It appears that the arterial wall owes its main mechanical characteristics, such as the
progressive stiffening and anisotropy, to collagen fibers and their orientations (Holzapfel
et al., 2000). Many constitutive models have been proposed to account for the distribution
of collagen fibers (e.g. Lanir et al., 1996; Gasser et al., 2006). These models use exponential functions for collagen fibers in their constitutive models and in general fitted well the
progressive stiffening with the increasing stretch. When there exist measurement noises in
experimental data, however, one often experiences difficulty in obtaining a good fit in the
high stretch region, which also causes inaccurate estimation of vascular stiffness. In this
paper, a weighted nonlinear least squares (WNLS) optimization is presented to estimate the
parameters of a multi-fiber model of arteries (Baek et al., 2007a; Hu et al., 2007; Masson
et al., 2008) considering the uncertainty due to the measurement errors in all variables. An
uncertainty model is first derived from the constitutive equation assuming that experimental
measurements are corrupted by the independent and identically distributed Gaussian noise.
Then, the WNLS optimization is formulated using the inverse of the covariance matrix of
the uncertainty as a correct weight factor. The proposed technique is evaluated using multiple sets of fictitious data containing the measurement errors (with multiple noise levels)
in all variables and the estimation results are compared with those of the conventional NLS
optimization (without a proper weight factor).
In parameter estimation, the larger number of parameters for a model provides more
flexibility and generally gives better fitting, i.e., decreases the residual error. However, too
many parameters increase the bias error as well as the complexity of the model. When

9

the degrees of freedom in parameter estimation exceed the information in the data, the
resulting estimated model cannot be generalized beyond the fitting data. In statistics, this
phenomenon is usually referred to as overfitting (Naik et al., 2007). Therefore, there is a
tradeoff between descriptive accuracy and parsimony, which can be addressed by a model
selection criterion (Wagenmakers and Farrell, 2004). Multiple model selection criteria are
utilized to investigate an optimal number of parameters (or fiber families) for two different
arteries.

2.2
2.2.1

METHODS
Constitutive relations for the mechanical behavior of the
passive artery

Cyclic inflation of an arterial segment at multiple fixed lengths is a typical biaxial test in vascular mechanics. The measurements include the axial force and internal pressure (recorded
by the load cell and pressure transducer) versus changes in the diameter and axial stretch
(recorded by the optical system). Experimental data from the test is used to determine
appropriate constitutive relations based on the assumption of an ideal cylindrical geometry of the artery (Humphrey, 2002). In the current study, a microstructurally motivated,
multi-fiber family model is used. That is, a constitutive strain energy including an isotropic
neo-hookean strain energy function (associated with elastin) and strain energy functions due
to multiple collagen fiber families is assumed (Holzapfel et al., 2000; Baek et al., 2007a):
)
}
(
(k) (k) 2
2 −1 ,
− 1)
exp c2 (λ

(k) {

∑ c
c
1
ˆ
W = (I1 − 3) +
(k)
2
4c
k

(k)

where, c, c1

(k)

and c2

(2.2)

2

(k)

are material parameters, such that c, c1

(k)

and c2

0, (Holzapfel,

2006; Holzapfel et al., 2004). I1 = tr C, where C is the right Cauchy-Green deformation

10

tensor. λ(k) is the stretch of the k th fiber family, given by
√
λ(k) =

(

λθ sin α(k)

)2

)2
(
+ λz cos α(k) ,

(2.3)

where α(k) is the orientation of the fiber family, and λz and λθ are the axial and circumferential stretches. The intramural stress can be obtained as
2 ∂W T
ˆ
T= F
F ,
J ∂C

ˆ
T = −pI + T,

(2.4)

where p is a Lagrange multiplier, F is the deformation gradient, and J is its determinant.
For an ideal straight tube, the Eq. (2.4) expands to the following three equations
(
)
}
( r )2 ( r )2 ∑ { (k) (k) 2
2
(k)
Tθθ = −p + c
+
c1 (λ
− 1)exp c2 (λ(k) − 1)2 sin2 α(k) , (2.5)
R
R
k
(
)
}
∑ { (k)
2
2
(k)
c1 (λ(k) − 1)exp c2 (λ(k) − 1)2 cos2 α(k) ,
Tzz = −p + cλz 2 + λz 2
(2.6)
( R )2
Trr = −p + c
.
rλz

k

(2.7)

For a thin membrane, the transmural pressure Pi and the axial force Fz can be approximated
by (See Appendix A)

Pi =

(
)
ˆ
ˆ
h Tθθ − Trr
rm
(

)

(

)

2
2
ˆ
ˆ
Fz = 2πrm h Tzz − Trr − π rm − rc Pi ,

(2.8)
(2.9)

ˆ ˆ
ˆ
where, rm = (ri +ro )/2, h = ro −ri , and rc is the radius of cannula. Trr , Tθθ and Tzz are the
normal components of stress in the radial, circumferential, and axial directions respectively.
Thus, by substituting (3.1)-(2.4) into (2.8) and (2.9), the theoretical relation between a force
vector y = [Fz , Pi ]T and a displacement vector x = [λz , do ]T can be derived as in Eq. (2.1).

11

2.2.2

WNLS optimization for a biaxial test of the artery

In experiments, all variables are measured with errors and we denote the measured variables
ˆ
ˆ
as xn = [λz (tn ), do (tn )]T and yn = [Fz (tn ), Pi (tn )]T at time tn for n = 1, . . . , m. The true
˜
˜
values of the variables are denoted by xn and yn , which are corrupted by the measurement
errors, ϵn and en during the experiment, i.e. the measured variables can be written as

ˆ
˜
xn = xn + ϵn

(2.10)

ˆ
˜
yn = yn + en ,

(2.11)

where ϵn ∼ W N (0, Σϵ ) and en ∼ W N (0, Σe ) are assumed to be independent and identically
distributed Gaussian noises with zero means and corresponding covariance matrices.
In order to consider the measurement errors in all variables, the total least squares estimation problem can be formulated by the following objective function (Beck and Arnold,
1977; Schwetlick and Tiller, 1985):

S=

m
∑

T





ˆ
ˆ
y − f (xn , Θ) 
 yn − f (xn , Θ) 
−1  n

 Σn 
,
ˆ
ˆ
xn − xn
xn − xn
n=1

(2.12)

where Σn is a collective covariance matrix of the measurement errors. Then, we need to
solve for 2m + N (Θ) unknowns to minimize Eq. (2.12), where m is the number of data
points and N (Θ) is the number of parameters in Θ. Considering a nonlinear function, such
as Eq. (2.1), an iterative scheme (e.g., Gauss-Newton method) has to be employed. Solving
the nonlinear regression problem with such a large number of unknown variables is very
difficult. Schwetlick and Tiller (1985) stated that for the NLS optimization, solving the total
least squares estimation problem using a standard software “cannot be recommended unless
the problem is small.” Consequently, in this study, we use the classical parameter estimation

12

method with only N (Θ) unknowns in the objective function:

S=

m
∑(

ˆ
yn − f (ˆ n , Θ)
x

)T

(
)
ˆ
Wn yn − f (ˆ n , Θ)
x

(2.13)

n=1

where Wn are appropriately chosen weight matrices. The problem is, then, to solve for Θ
by minimizing Eq. (2.13) with correct weights obtained by the uncertainty model which
will be referred to as the WNLS optimization. Let the uncertainty model vn represent
ˆ
the uncertainty in yn − f (ˆ n , Θ). The uncertainty in f (ˆ n , Θ) is propagated from the
x
x
measurement noise of ϵn , and it can be approximated by using the Taylor series of f with
respect to ϵn and ignoring the higher order terms of ϵn in the series. Then, the uncertainty
ˆ
model vn for yn − f (ˆ n , Θ) can be obtained as:
x
vn = −

∂f
(˜ n , Θ) ϵn + en .
x
∂x

(2.14)

Note that the uncertainty model for the force measurement now includes the effect of measurement noise in the displacement. Eq. (2.14) shows that the uncertainty increases with
an increase in ∂f /∂x, which is a stiffness term in vascular mechanics. To incorporate the
uncertainty model into Eq. (2.13), the inverse of the covariance matrix of the uncertainty has
to be used as the weight factor (Beck and Arnold, 1977; Emery et al., 2000). The covariance
matrix for vn is derived as:
(
T
Σvn = E(vn vn ) =

∂f
∂x

)

(
Σϵ

)
∂f T
+ Σe ,
∂x

(2.15)

where E(A) is the expectation of A. Finally, using the uncertainty model Eq. (2.14), the
objective function Eq. (2.13) can be written as:

S=

m
∑(

ˆ
yn − f (ˆ n , Θ)
x

)T

n=1

13

(
)
ˆ
x
Σ−1 yn − f (ˆ n , Θ) .
vn

(2.16)

The correct modeling of the measurement error and the formulation of the WNLS optimization will provide the minimum estimation error variance. For the computation, however, the
initial estimates Θo are obtained by the NLS optimization without the uncertainty model.
Then, the covariance matrix Σvn in Eq. (2.16) is approximated using the estimates Θo
ˆ
and xn . With Σvn , the new estimates of Θ1 are obtained by minimizing Eq. (2.16). The
covariance matrix Σvn is then updated iteratively using the previous estimates for the next
optimization. The estimates in each iteration are obtained by using constrained optimization
in Matlab (The Mathworks Inc.) with multiple choices of initial points, and one to three
iterations are used to obtain the final estimates in this paper.
Vascular stiffness changes according to alterations in physiological and pathophysiological
conditions, such as aging (O’Rourke and Hashimoto, 2007), hypertension (Hu et al., 2007),
during pregnancy (Hu et al., 1998), and diabetes mellitus (Oxlund et al., 1989). Accurate
assessment of the arterial wall stiffness in physiological range can play an important role
in understanding the pathophysiology and progression of vascular diseases. The linearized
circumferential stiffness from both the WNLS and NLS methods is calculated and compared
for its accuracy (see Baek et al., 2007a for linearization). The results demonstrate that our
proposed scheme provides a more accurate assessment of the arterial stiffness.

2.2.3

Model selection criteria for optimal number of parameters

In order to find the optimal number of parameters (or number of fiber families) for the multifiber family model (Eq. (3.1)), three different criteria for model selection are utilized; Akaike
information criterion (AIC; Glatting et al., 2007), a modified form of AIC (AICc; Glatting
et al., 2007), and the root mean square error measure (RMS; Holzapfel et al., 2005), given

14

by
(

)
S
+ 2 (N + 1)
AIC = m ln
m
2 (N + 1) (N + 2)
AICc = AIC +
m−N −2
√
S
RMS =
m−N

(2.17)
(2.18)
(2.19)

where S, N , and m are, respectively, the residual, the number of parameters, and the sample
size. AIC was first introduced by Akaike (1974, 1981) based on the concept of entropy to
describe the tradeoff between bias and variance in model construction. AIC has been used as
a model selection criterion that selects an optimal model by considering both the precision
of fitting and the complexity of the model (Anderson et al., 1994). While the first term
in the right hand side of Eq. (2.17) decreases with a decrease in the residual, the second
term penalizes it for increasing the size (number of parameters) of the model and prevents
overfitting. RMS in Eq. (2.19) has also been used as a heuristic criterion in selecting among
models although no statistical justification exists (Myung, 2000). The number of parameters
that minimizes a given criterion is considered as being optimal for the model.
In the multi-fiber family model, the number of parameters increases with an increase
in the number of fiber families. The fibers are assumed to be symmetric with respect to
the axial axis on the vessel wall. In order to evaluate the model selection criteria, the
residual Eq. (2.16) is obtained by increasing the number of fiber families from 2 to 10 for
each data set. Briefly, the 2-fiber-family model has 4 independent parameters (c, c1 = c2 ,
1
1
c1 = c2 and α1 = −α2 ). For the 3-fiber-family model, one more fiber family is added in
2
2
either circumferential direction or axial direction to minimize the residual, which results in
6 (independent) unknown parameters. The 4-fiber-family model has two symmetric fibers in
helical directions, one in the axial direction, and one in the circumferential direction resulting
in 8 independent parameters. From 6 fiber families, two additional symmetric fibers are
added at each step, yielding three additional parameters (an angle and two parameters for
15

the exponential function). Hence, 11, 14, and 17 independent parameters are assumed for
6-, 8-, and 10-fiber-family models respectively.
The model selection criteria are specific to the experimental data. Hence, experimental
data from two rabbit basilar arteries and three mouse carotid arteries are tested in order to
investigate optimal number of fiber families.

2.3
2.3.1

RESULTS
Comparison between the proposed WNLS optimization and
the conventional NLS optimization

In order to demonstrate the effectiveness of the proposed WNLS method, noise-corrupted,
fictitious experimental data were generated and the parameter estimation was performed
using the WNLS method as well as the standard approach (NLS) without an uncertainty
model. Assuming a set of “true” parameters, the synthesized true data was generated for
the inflation test at fixed axial stretches (λz = 1.2, 1.3, 1.4), i.e., the pressure and axial
force data were obtained from Eq. (2.1) for a given radius and axial stretch using the true
parameters. Gaussian noises with given noise levels were numerically generated and added
to the synthesized data of {do , λz , Fz , Pi } to produce fictitious experimental data. The noise
level was defined as the ratio of the standard deviation of the noise to the maximum value
of the data. Figure 2.1 shows the pressure vs. radius plots for λz = 1.2, 1.3, 1.4 for the true
as well as the fictitious data at noise levels of 0.005 and 0.01.
Then, the NLS and proposed WNLS methods were used to estimate parameters for each
fictitious data set. Figure 2.2 shows pressure-radius and axial force-radius plots using the true
material parameters and the data calculated with estimated parameters of the 4-fiber-family
model for both the NLS and WNLS methods. Evidently, the graphs show better fitting
curves for pressure and axial force when the uncertainty model is incorporated (Figures
16

Pressure (kPa)

Pressure (kPa)

14
12 (a)
10
8
6
4
2
0
−2
0.028 0.032 0.036 0.04 0.044
Radius (cm)

14
12 (b)
10
8
6
4
2
0
−2
0.028 0.032 0.036 0.04 0.044
Radius (cm)

Figure 2.1: Pressure versus radius plots for the 4-fiber-family model for a set of true data
(solid lines) and fictitious noisy data (dotted lines) at noise levels of 0.005 (a) and 0.01 (b)
and different axial stretches λz = 1.2, 1.3, 1.4.
2.2(b) and (d)).
Figure 2.3 shows another set of fictitious data with higher slopes in the high stretch region
and the corresponding fitting curves using both the NLS and WNLS methods. The WNLS
significantly improved the slope of the fit especially in the high stretch region (Figures 2.2
and 2.3).
The advantage of the WNLS method over the standard NLS method is more obvious for
the data with higher slopes (Figure 2.3). Table 2.1 summarizes all true and estimated values
of parameters corresponding to Figures 2.2 and 2.3.
For a quantitative measure of estimation errors, the following normalized error was defined as (Baek et al., 2007a):
√∑
(√ ∑
)
2
1
(Pest − Ptrue )
(Fest − Ftrue )2
∑
∑
e=
+
,
2
(Ptrue )2
(Ftrue )2

(2.20)

where Pest and Ptrue are the estimated and true values of pressure. Fest and Ftrue are the
estimated and true values of axial force. Normalized errors, e, were obtained and averaged
over a set of five different random noise sequences at each noise level. Table 2.2 shows
the averaged, normalized errors obtained from both the NLS and WNLS methods at noise
17

Pressure (kPa)

14
Estimated
12
True
10
8
6
4
2
0
(b)
−2
0.03 0.034 0.038 0.042 0.046

14
Estimated
12
True
10
8
6
4
2
0
(a)
−2
0.03 0.034 0.038 0.042 0.046

4

2

2

0

0

−2

−2

−4

−4

−6

Force ( mN )

4

−6

Estimated
(c)
True
−8
0.03 0.034 0.038 0.042 0.046
Radius ( cm )

(d)
Estimated
True
−8
0.03 0.034 0.038 0.042 0.046
Radius ( cm )

Figure 2.2: Comparison between the NLS (a and c) and WNLS (b and d) optimizations
for the 4-fiber-family model. Pressure and axial force versus radius are plotted using true
parameters (dotted lines) and estimated parameters (solid lines) for the data with moderate
slopes at different axial stretches λz = 1.2, 1.3, 1.4

18

Pressure ( kPa )

4
Force ( mN )

14
Estimated
12
True
10
8
6
4
2
0
(b)
−2
0.03 0.034 0.038 0.042 0.046

14
Estimated
12
True
10
8
6
4
2
0
(a)
−2
0.03 0.034 0.038 0.042 0.046

2

4
2

(c)

0

0

−2

−2

−4

−4

−6

(d)

−6

Estimated
True
−8
0.03 0.034 0.038 0.042 0.046
Radius (cm)

Estimated
True
−8
0.03 0.034 0.038 0.042 0.046
Radius ( cm )

Figure 2.3: Comparison between the NLS (a and c) and WNLS (b and d) optimizations
for the 4-fiber-family model. Pressure and axial force versus radius are plotted using true
parameters (dotted lines) and estimated parameters (solid lines) for the data with stiff slopes
in the high stretch region at different axial stretches λz = 1.2, 1.3, 1.4.

19

Table 2.1: True and estimated parameters for the fictitious data with moderate and high
slopes (corresponding to Figures 2.2 and 2.3) and the experimental data from mouse carotid
and rabbit basilar arteries (corresponding to Figures 2.4(a) and 2.4(b)).

c(kPa) c1 (kPa) c2 (kPa) c3 (kPa)
1
1
1

c1
2

c2
2

c3
2

α1 (◦ )

Fig. 2.2

True
NLS
WNLS

0
0
0

0.36
0.71
0.39

7.9
14.2
8.7

1.2
1.1
1.1

4.3
3.72
4.23

1.64 3.07
1.12 3.1
1.56 3.15

40.3
42
41.4

Fig. 2.3

True
NLS
WNLS

0
0
0

0.28
1.3
0.28

1.73
9.34
3.4

0.2
0.22
0.27

4.73
3.05
4.73

2.96 4.92
1.52 4.63
2.36 4.53

45.6
40.5
44.5

Fig. 2.4(a)

WNLS

42

7.29

14

0.0003

0.152

0

3.11

55.9

Fig. 2.4(b) WNLS

0

18

14

6.2

1.37

1.48

5.64

47.9

levels of 0.005, 0.008, 0.01, 0.02, and 0.03. The WNLS method resulted in normalized errors
smaller than those from the NLS method at all noise levels.
True and estimated values of the circumferential stiffness at Pi = 75mmHg and λz = 1.3
corresponding to Figures 2.2 and 2.3 are shown in Table 2.3. The WNLS method resulted in
much better estimation than the NLS method for both cases. For example, corresponding
to the data of Figure 2.2, the error of estimated stiffness using the NLS method was about
15 percent whereas it was less than 3 percent when using the WNLS method.

2.3.2

Optimal number of fiber families

Experimental data from two rabbit basilar arteries (Baek et al., 2007a) and three mouse
carotid arteries (Dye et al., 2007) were used to investigate the optimal number of fiber families
for the proposed model. Figures 2.4(a) and 2.4(b) depict experimental data for these arteries
as well as the best-fit curves using the WNLS optimization with the 4-fiber-family model.
For all sets of data tested, residuals for the WNLS optimization reduced significantly up to

20

Table 2.2: Averaged normalized error e at different noise levels using the NLS and WNLS
methods.

Noise level

0.005

0.008

eN LS
eW N LS

0.0093 0.018
0.009 0.01

eN LS
eW N LS

1.033

1.8

0.01

0.02

0.03

0.03
0.0123

0.07
0.025

0.12
0.078

2.44

2.8

1.54

Table 2.3: True and estimated linearized stiffness at Pi = 75mmHg and λz = 1.3 corresponding to Figures 2.2 and 2.3 using the NLS and WNLS methods.

True stiffness(kPa) NLS stiffness(kPa) WNLS stiffness(kPa)

Fig. 2.2
Fig. 2.3

1114
1580

944
1105

1085
1381

3-fiber-family model (6 parameters), and then slightly decreased for further increase in the
number of fiber families (see e.g. Figures 2.4(c) and 2.4(d)). Interestingly, three different
criteria led to an identical optimal number of parameters for each mouse and rabbit data,
which were found to be 11 parameters (6 fiber families). Estimated parameters corresponding
to Figures 2.4(a) and 2.4(b) are listed in Table 2.1.

2.4

DISCUSSION

Parameter estimation using the NLS optimization has been widely used in characterizing
mechanical behavior of soft tissues from the experiments. In this study, we proposed an
improved parameter estimation technique considering the measurement errors in variables.
If the measurement errors in independent variables are negligible (i.e., Σϵ ∼ 0), then the
proposed method becomes identical to the conventional NLS method. However, in many
21

2500
2300

AICc
Residual

Pressure (kPa)

4
3
2

4000
3500

1
2100
0
4 6 8 10 12 14 16 18
Number of Parameters

(d)

AICc
Residual

7
Residual (x 10 )

AICc

2700

(c)

4500

5

AICc

2900

14
(b)
10
6
2
−2
m)
Ax1.35
s (c 0.045
diu 0.04
ial 1.25
Ra0.035
St
re 1.15 0.03
tch

cm)
us ( 0.04
i
Rad 0.03
0.02

7
Residual (x 10 )

Pressure (kPa)

25
(a)
20
15
10
5
0
2
Ax 1.9
ial 1.8
St 1.7
re 1.6
tch

2.5
2
1.5
1
0.5

0
3000
4 6 8 10 12 14 16 18
Number of Parameters

Figure 2.4: Pressure-radius-axial stretch curves are plotted from fitting of the experimental
data for the mouse carotid artery (a) and rabbit basilar artery (b) using the WNLS optimization. Circles show data points and solid lines show the estimated values. AICc and the
residual are plotted against the number of parameters; 4, 6, 8, 11, 14, 17 (or the number of
fiber families; 2, 3, 4, 6, 8, 10) for the mouse carotid (c) and rabbit basilar (d) arteries using
the WNLS optimization.

22

studies (i.e., arterial inflation and extension tests) such measurement errors are not negligible.
For example, Saravanan et al. (2006) reported that when the deformation was approximated
by a linear polynomial using three markers, the error in the first invariant of the right CauchyGreen deformation tensor was ± 0.06, which is comparable with noise levels we considered
in this study. Although these problems can be treated as total least squares problems or
errors-in-variables models (Beck and Arnold, 1977; Huffel and Vandewalle, 1991), solving
a nonlinear problem with a large number of unknown variables is still very challenging.
Instead, we developed a WNLS technique based on accurate modeling of uncertainty given
by Eq. (2.14). It was shown that the WNLS optimization with a proper uncertainty model
improves the quality of parameter estimation significantly compared to the conventional NLS
optimization at all noise levels. Especially, for a collagenous tissue, ∂f /∂x is larger in the
high stretch region, so it produces a higher level of uncertainty propagation and, hence, the
estimation results are biased when the displacement measurement error is not considered
within a proper uncertainty model (Figure 2.3). The advantage of using the WNLS was
more pronounced at higher noise levels. The WNLS method provided a better fit in the high
stretch region and proved to be advantageous when estimating the linearized stiffness within
the physiological range.
In the current work, synthesized data was used along with Gaussian noises to evaluate the
WNLS optimization. In experiments, however, the measurement noise should be carefully
characterized. Furthermore, although the proposed WNLS optimization helps eliminate the
biased error due to the measurement errors, it is noteworthy that parameter estimation can
be limited by the model error induced by the chosen constitutive relation.
The presented constitutive model is similar to the one developed by Holzapfel et al. (2000)
and Holzapfel et al. (2004), but has been used in thin wall models (Baek et al., 2007a; Hu
et al., 2007; Masson et al., 2008). It has also been utilized in modeling of vascular adaptation
during progression of vascular diseases (Baek et al., 2006, 2007b). However, the vessel wall
presents more dispersed fiber orientation in the adventitia and intimal layers and the use of

23

an orientation density function has been proposed (Lanir et al., 1996; Gasser et al., 2006).
The choice of the functional form of the constitutive relation may involve multiple factors
such as the anatomic location of the artery, available microstructural information, and the
specific application and, hence, it is beyond the scope of the present work. Instead, our
focus has been on the situation that one has a functional form of a constitutive relation and
needs to find optimal number of parameters and best parameters from experimental data.
The presented parameter estimation method is, however, general enough and can be utilized
with various models for vascular mechanics.
Three criteria were used to evaluate optimal number of parameters (fiber families) for the
chosen constitutive model of arteries. The best model between several competing models
is one that provides an adequate account of the data while using a minimum number of
parameters. Based on the available data and using three different criteria, we found that the
model with 11 parameters (6 fiber families) minimized our criteria.
In closing, the need for the optimal design of experiments and optimal sampling protocols
in vascular mechanics is emphasized (for example, see Lanir et al., 1996). In optimally
designed experiments, the effect of parameter changes is maximized with respect to the
noise and, hence, the estimation error variance can be reduced.

24

Chapter 3
COMPARISON OF DIFFERENT CONSTITUTIVE MODELS
OF VASCULAR TISSUE MECHANICAL BEHAVIOR FROM
PARAMETER ESTIMATION PERSPECTIVE

3.1

Constitutive modeling of arterial wall

Due to large deformation of vascular tissue, constitutive modeling based on finite elasticity
has been extensively used to understand the mechanical behavior of healthy or diseased arteries. Constitutive modeling of arteries requires estimation of constitutive parameters from
ex-vivo biaxial or inflation-extension experiments. In chapter 2, an aspect of the parameter estimation in vascular mechanics was discussed with regard to the measurement error.
The form of the constitutive relation/model, assumed to describe the mechanical behavior
of the tissue, was considered to be known and fixed. Nonetheless mechanical behavior of
arteries may vary between species and even may be different within a subject depending on
the anatomical location (See Fig. 2.4a and b, for example). Thus the question remained
is which constitutive relation adequately describes the mechanical behavior arteries or even
whether one constitutive form is adequate to model all arteries. Variety of constitutive models, either phenomenological or microstructurally-motivated, has been suggested for arteries
in literature, some of which have been reviewed by Holzapfel et al. (2000). In this chapter,
given two sets of data mouse carotid and rabbit basilar arteries, we test different constitutive
25

relations by comparing their parameter estimation results.
Briefly, fibers and ground substance (mainly elastin) are main constituents that determine the mechanical behavior of the vascular tissues (Holzapfel and Ogden, 2010; Humphrey,
2002). Under low loading condition, crimped collagen fibers contribute little to the mechanical response while elastin is mainly responsible for the mechanical behavior. Arterial elastin
(and ground substance in general) has been generally assumed to be isotropic and modeled
as the neo-Hookian material (Gundiah et al., 2007). At higher loads the undulated fibers
begin to stretch out and contribute to the progressive stiffening of the tissue in the fiber
direction. That has motivated exponential forms of the strain energy function for collagen
with anisotropy effect (Fung, 1967; Holzapfel et al., 2000). An additive split for the strain
energy function is suggested as W = We + Wc , where We represents the strain energy due to
isotropic elastin (independent of fiber alignment) and Wc represents the anisotropic strain
energy due to collagen fibers (Holzapfel et al., 2000; Holzapfel and Ogden, 2010). Here we
consider different functional forms of the strain energies for isotropic and anisotropic parts in
contrast with the popular Holzapfel-type model. First we assume the anisotropy part as an
exponential form and use different models for isotropic part. Then we consider the isotropic
part as neo-Hookian and test different models for anisotropic part.

3.2

Holzapfel-type model

A multi-fiber family model, as introduced in chapter 2, has been suggested by Hu et al.
(2007) and Baek et al. (2007a) which is mainly inspired by Lanir (1983) and Holzapfel et al.
(2000). The strain energy function is rewritten as
(k) {

∑ c
c
1
ˆ
W = (I1 − 3) +
(k)
2
4c
k

(
)
}
(k) (k) 2
2 −1 ,
exp c2 (λ
− 1)

2

26

(3.1)

Normalized Residual

20
15
10
5
0

400
200
300
Radius (micron)

Pressure (kPa)

12
10
8
6
4
2
0
−2
0.03

0.2
0.1
0
−0.1
−0.2
−0.3

Normalized Residual

Pressure (kPa)

25

λ= 1.65
λ= 1.80
λ= 1.95
200
300
400
Radius (micron)

0.1
0

−0.1
λ=1.2
λ=1.3
λ=1.4
−0.3
0.03 0.035 0.04 0.045
Radius (cm)
−0.2

0.035 0.04 0.045
Radius (cm)

Figure 3.1: Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries
using Holzapfel-type model. Dots show the data points and solid lines show the estimated
values.
Accordingly, Fig. 3.1 shows the pressure-radius data and the estimated lines along with the
normalized residuals between data and estimated values (i.e., (Pest − Pdata )/Pnominal ) for
mouse carotid and rabbit basilar arteries. Existing trends of variation in residuals around
the zero, underscores the signature behavior associated with the model used to fit both sets
of data. However, the model shows a better fit to basilar artery data in general.

27

20
15
10
5
0

Normalized Residual

Pressure (kPa)

25

0.2
0.1
0

−0.1
−0.2

400
200
300
Radius (micron)

λ= 1.65
λ= 1.80
λ= 1.95
400
200
300
Radius (micron)

Figure 3.2: Pressure-radius data of mouse carotid artery and the estimated values (left)
with the corresponding normalized residuals (right) using Zulliger model. Dots show the
data points and solid lines show the estimated values.

3.3
3.3.1

Models for elastin contribution
Zulliger model

Zulliger et al. (2004a) proposed a modification to the strain energy function of elastin which
we combine with exponential part of the Holzapfel-type model as
(k) {

(
)
}
(k) (k) 2
2 −1 .
exp c2 (λ
− 1)

∑ c
c
1
ˆ
W = (I1 − 3)3/2 +
(k)
2
4c
k

(3.2)

2

Figure 3.2 show the data and fitted lines as well as the normalized residuals for carotid
artery. There is no improvement in terms of fitting and the trend in residuals with respect
to neo-Hookian form of strain energy for elastin.

3.3.2

Demiray model

Ogden and Saccomandi (2007) suggested the exponential model of isotropic part, proposed
by Demiray (1972) and Fung (1967), to be combined with anisotropic part of Holzapfel-type

28

model as
(
)
}
)
} ∑ c(k) {
c { (
(k) (k) 2
2 −1
1
ˆ
exp b(I1 − 3) − 1 +
W =
exp c2 (λ
− 1)
(k)
2b
4c

Normalized Resdiual

25
20
15
10

0.2
0.1
0

−0.1

5
0

Pressure (kPa)

2

−0.2

200
300
400
Radius (micron)

12
10
8
6
4
2
0
−2
0.03

0.035 0.04 0.045
Radius (cm)

λ= 1.65
λ= 1.80
λ= 1.95
200
300
400
Radius (micron)

Normalized Residual

Pressure (kPa)

k

(3.3)

0.05
0
−0.05
−0.1
−0.15 λ = 1.2
λ = 1.3
−0.2
−0.25 λ = 1.4
0.03 0.035 0.04 0.045
Radius (cm)

Figure 3.3: Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries
using Demiray model. Dots show the data points and solid lines show the estimated values.

29

3.4
3.4.1

Models for collagen contribution
Zulliger model

As described before, collagen fibers are originally undulated and under loading condition
more collagen fibers become straightened and engaged in carrying the increasing load (since
each collagen fiber is assumed to become straightened at a different stretch ratio). Based
on that Zulliger et al. (2004a) assumed a log-logistic probability distribution function for
engagement strain. Inspired by them, we consider
∑
c
ˆ
νk Ψk (ϵk )
W = (I1 − 3) +
2

(3.4)

k

∫∞

where ϵk = λk − 1 and Ψk (ϵk ) = ψf ∗ ρf = −∞ ψf (ϵk − x)ρf (x)dx. The distribution
function ρf and the strain energy of a single fiber ψf are given as

 0




ρf (ϵk ) =

ψf (ϵk ) =







(

)
ϵk ak −1

ak
bk
bk [ ( ϵ )ak ]2
1+ k
bk



 0

for ϵ ≤ 0
(3.5)
for ϵ>0

for ϵ ≤ 0

 c (ϵ − log(ϵ + 1)) for ϵ>0
 f k
k

(3.6)

where ak and bk are the parameters associated with kth fiber family (k = 1, ..., 4) and cf
is an intrinsic parameter for a single fiber. bk is a scaling parameter while ak defines the
shape of the distribution (Rezakhaniha et al., 2011). Figure 3.4 shows the fitting results and
the normalized residuals for both mouse carotid and rabbit basilar arteries. Apparently the
parameter estimation has not been improved using this model.

30

15
10

200
300
400
Radius (micron)

12
10
8
6
4
2
0
−2
0.03 0.035 0.04 0.045
Radius (cm)

Pressure (kPa)

0.2
0.1
0

−0.1 λ= 1.65
−0.2 λ= 1.80
λ= 1.95

5
0

Normalized Residual

20

Normalized Residual

Pressure (kPa)

25

200
300
400
Radius (micron)

0.05
0
−0.05
−0.1
−0.15
−0.2 λ= 1.2
λ= 1.3
−0.25 λ= 1.4
0.03 0.035 0.04 0.045
Radius (cm)

Figure 3.4: Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries
using Zulliger model. Dots show the data points and solid lines show the estimated values.

31

3.4.2

Lanir model

Again based on the idea of gradual engagement of collagen fibers, Lanir (1979) considered a
normal distribution for the fiber engagement and its load bearing capacity. He also considered
the engaged individual fibers to be linearly elastic. Similarly, we consider
∑
c
ˆ
W = (I1 − 3) +
νk Ψk (ϵk ),
2

Ψk (ϵk ) =

k

∫ λ
k
1.0

ψf (x)ρf (x)dx

(3.7)

where ψf (x) = 1/2cf ((λk − x)/x)2 is the strain energy of the a single fiber and ρf is the
normal distribution defined as
( (x − c )2 )
k
,
ρf (x) = bk exp −
2a2
k

where

bk =

1
√

ak 2π

.

(3.8)

Figure 3.5 depicts the fitting results and the normalized residuals for both mouse carotid
and rabbit basilar arteries.

3.5

Discussion

Concerning the difference between species, Holzapfel-type model showed relatively a better
fit to the basilar artery data (as well as lower residuals). This can be due to the simple
exponential trend observed for rabbit basilar artery whereas mouse carotid artery showed a
sigmoid-like trend emphasizing more contribution for elastin under low loads. Comparison
between different models did not show a clear advantage for other constitutive models over
the Holzapfel-type model in terms of parameter estimation. Holzapfel-type together with
Demiray-type showed the best results for the parameter estimation and fitting residuals. For
an ideal model, residuals should be oscillating around zero without any obvious signature
trend. Whereas the signature behavior can be obviously seen in residuals of all models used
in the study.
However, as mentioned earlier, Holzapfel-type model represents multiple fiber families of
32

15
10
5

Pressure (kPa)

0

200
300
400
Radius (micron)

12
10
8
6
4
2
0
−2
0.03 0.035 0.04 0.045
Radius (cm)

Normalized Residual

20

Normalized Residual

Pressure (kPa)

25

0.2
0.1
0
−0.1
−0.2

λ= 1.65
λ= 1.80
λ= 1.95
200
300
400
Radius (micron)

0.05
0
−0.05
−0.1
−0.15
−0.2
−0.25

λ= 1.2
λ= 1.3
λ= 1.4
0.03 0.035 0.04 0.045
Radius (cm)

Figure 3.5: Pressure-radius data and the estimated values (left) with the corresponding
normalized residuals (right) for mouse carotid (top) and rabbit basilar (bottom) arteries
using Lanir model. Dots show the data points and solid lines show the estimated values.
collagen fibers along with elastin. This is favorable to be used in G&R models where arteries
are considered as a mixture of constituents which are produced and removed at their own
turnover rates.

33

Chapter 4
ESTIMATION OF THE ARTERIAL WALL THICKNESS AND
ANISOTROPY FOR IMAGE-BASED MODELS OF ARTERIES

4.1

INTRODUCTION

Many of the computational models of vascular adaptation (Baek et al., 2006, 2007b; Figueroa
et al., 2009; Kroon and Holzapfel, 2009b; Watton and Hill, 2009) have been built upon
the theoretical framework of modeling tissue G&R presented by Humphrey and Rajagopal
(2002). They introduced a constrained mixture approach focusing on stress-mediated mass
production and removal in evolving stressed configurations. They also offered key remarks
that are central to guiding the later development of theories of soft tissue G&R. One of
their key remarks is “Normal growth and remodeling tends to be a stable dynamical process,
one that seeks to optimize structure and function with respect to yet unidentified parameters. In comparison to processes during development, there appear to be genetic and perhaps
epigenetic constraints on this optimization process during maturity”. Furthermore, they
emphasized a need to identify both a set of optimization parameters and the associated constraints. Most of the previous computational simulations of vascular adaptation, however,
have been developed using idealized geometries for which the identification of homogeneous
parameters does not pose a problem.
This chapter addresses two technical challenges associated with patient-specific modeling

34

of AAA evolution and proposes possible solutions. First, as stated earlier, theory of G&R
is based on a key assumption, the existence of mechanical homeostasis (Humphrey, 2008;
Kassab, 2008), whereas it is difficult to prescribe the in vivo parameters such that the
assumption of a homeostatic state is satisfied at every point in the vessel wall model. For an
idealized model, where the blood vessel is assumed to be an ideal thin hollow cylinder, the
in vivo material properties are typically assumed to be uniform over the domain. When a
medical image-based geometric model is used, however, it is not a trivial task to prescribe the
distribution of material and structural parameters such as thickness and fiber orientations.
Second, another difficulty associated with using an image-based model stems from the
fact that the in vivo image is obtained under the pressure and the stress-free configuration
is not available. Hence, it is difficult to maintain the original patient-specific model in
a computational simulation under the in vivo pressure. Inverse elastostatic methods have
been pursued to estimate the stress-free state from a pre-deformed in vivo geometry obtained
medical images (Lu et al., 2008; Zhou et al., 2010). Others have used a Lagrangian-Eulerian
formulation to obtain the meaningful prestressed state (Gee et al., 2010, 2009).
In this chapter, first, the constitutive parameters intrinsic to the material are estimated
by fitting the ex vivo biaxial mechanical test data of a healthy human aorta using the same
parameter estimation approach presented in chapter 2. Second, an optimization problem is
solved to estimate the distributions of the wall thickness and anisotropy such that homeostasis is maintained while the geometry deviates minimally from the in vivo configuration.

4.2
4.2.1

METHOD
Estimation of material constitutive parameters

As the first step, the constitutive parameters are estimated by fitting biaxial mechanical test
data of a healthy human aorta (Vande Geest et al., 2004, 2006). Here, the kinematics of the
biaxial test and the constitutive relations are briefly explained.
35

κ1 ( 0 )
n

1

κ2 ( 0)
n

κ in ( 0)

G 2 (0)

G (0)

G i (0)
κR

P
F = F I FR

FR
κI

λ2

FI

λ1

Figure 4.1: Kinematics of the deformation associated with biaxial mechanical test and the
corresponding deformation gradients. λ1 and λ2 are stretches in circumferential and axial
directions during the biaxial test

36

Figure 4.1 shows a schematic drawing for the kinematics of deformation related to a
biaxial test of a healthy aorta. The in vivo configuration of a healthy aorta is assumed
to be the prestressed reference configuration κR , whereas κI represents the intermediate
configuration of the square-cut sample under the traction-free condition. The deformation
gradient FR corresponds to the mapping from κR to κI . It is assumed that there is no
active tone presented during the biaxial test. The deformation gradient FI corresponds
to the mapping from κI to the deformed configuration during the biaxial test, resulting in
F = FI FR . Assuming incompressibility in an ideal geometry,
{

1
R
R
FR = diag F1 , F2 , R R
F1 F2
{
}
1
I = diag λ , λ ,
F
,
1
2
λ1 λ2

}

(4.1)

R
R
where F1 , F2 < 1.0 and λ1 , λ2 > 1.0.

The arterial wall is assumed to be a mixture of constituents ‘i’ such as elastin (i = e),
multiple collagen families (i = 1, ..., k, ..., 4), and smooth muscle (i = m). The strain energy
∑
∑
m
of the mixture per unit reference area is w = i wi = we + k wk + wm + wact and the
membrane stress is given as (Baek et al., 2006; Humphrey, 2002)

T=

2 ∂w T
F
F ,
J ∂C

(4.2)

where J is a determinant of the 2-D deformation gradient F and C = FT F. The stretches of
the smooth muscle (SM) and collagen fiber ‘k’ from their natural (stress-free) configuration
to the current configuration are given as
λk = Gc λk
n
h

(4.3)

λm = G m λ1 ,
n
h

(4.4)

37

where Gm and Gc are homeostatic stretches of SM and collagen. We define a new tensor
h
h
˜
Ge = diag

{
}
1
e , Ge ,
G1
2 Ge Ge ,
1 2

(4.5)

which represents a mapping from the natural configuration of elastin to the reference configuration such that,
˜
Fe = FGe ,
n

˜
˜
Ce = Fe T Fe = [Ge ]T CGe .
n
n
n

(4.6)

Strain energies of the constituents i per unit reference area, wi , are given as
)
c ( e
1
e
we (Ce (t)) = M e 1 Cn[11] + Cn[22] + e
−3
n
e
e
2
Cn[11] Cn[22] − Cn[12] 2
}
)2 ]
c { [ (
wk (λk ) = M k 2 exp c3 (λk )2 − 1
−1
n
n
4c3
{ [ (
}
)2 ]
c
wm (λm ) = M m 4 exp c5 (λm )2 − 1
−1
n
n
4c5
1 (λM − λ1 )3 }
S{
m
λ1 +
wact = M m
,
ρ
3 (λM − λo )2

(4.7)
(4.8)
(4.9)
(4.10)

e
e
e
where M i is the mass per unit reference area for the constituent i. Cn[11] , Cn[22] and Cn[12]

are components of Ce . λM and λo are stretches at which the SM contraction is maximum
n
and at which active force generation ceases, S is the stress at the maximum contraction of
SM.
Components of FR are obtained by considering stress as a function of deformation graˆ
dient, i.e. T = T(F), and assuming that membrane stresses vanish at
F = FR such that
ˆ
T(FR ) = 0.

(4.11)

Material parameters (summarized in Table 4.1) are determined in three different ways. The
first set of parameters, such as density (ρ), mass fraction of constituents (ν i ) and (λ0 , λM ,
S), are prescribed from the literature (He and Roach, 1994; Holzapfel et al., 2002; Menashi
38

Table 4.1: Summary of material parameters used in the optimization.

Elastin:

c1 = 50.6 N m/kg, Ge = 1.22, Ge = 1.23,
1
2
e = 0.2.
ν

Collagen:

c2 = 3195 N m/kg, c2
= 0.1c2 , c3 = 25.0,
c = 1.034, σ c = 143 kP a,
Gh
h

(comp)

ν k = [0.06, 0.06, 0.24, 0.24],
αk = [0◦ , 90◦ , 45◦ , 135◦ ].
c4 = 16.45 N m/kg, c5 = 14.14, Gm = 1.165,
h

Smooth muscle:

m
ν m = 0.2, σh = 81 kP a, λM = 1.4,

λ0 = 0.8, S = 54 kP a.
cyc

Homeostatic cyclic stretch: λh = 1.02.
Density:
ρ = 1050 kg/m3 .
et al., 1987). The second set of parameters (c1 , c2 , c3 , c4 , c5 , Ge , Ge , Gc , Gm ) are estimated
1
2
h
h
via the weighted nonlinear least squares optimization, proposed in chapter 2 using the biaxial
mechanical test data of healthy human aorta (Vande Geest et al., 2004, 2006). Figure 4.2
shows the biaxial data as well as the fitted values using the estimated parameters.
Although the existence of mechanical homeostasis in vasculature is generally accepted,
the theoretical formulation that describes vascular adaptations in response to diverse stimuli
is not completely established yet. Nevertheless, we utilize scalar measures of stress as the
intramural stress of constituents (Baek et al., 2006; Figueroa et al., 2009)
∑
σ k = ||(
ν k σ k )nk ||,

σ m = ||σ m nm ||,

(4.12)

k

where σ k and σ m are the stresses of the k th collagen fiber and SM, respectively, and nk
and nm are unit vectors in the directions of the k th collagen fiber and SM. Then, using the
estimated parameters and assuming an idealized geometry, the target homeostatic values
m c
(σh , σh ) are calculated.

39

Circumferential stress, t11(kPa)

70
60
50
40
30
20
10
0

1

1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08

Cirumferential stretch,λ1

Axial stress, t22 (kPa)

70
60
50
40
30
20
10
0

1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08

Axial stretch, λ2
Figure 4.2: Stress versus stretch plots in both directions. Biaxial test data (empty markers) and fitted values (filled markers) using the estimated parameters. Different markers
correspond to a different ratios of tension applied simultaneously in both directions.

4.2.2

Inverse optimization problem statement

As the next step, the distributions of wall thickness and material anisotropy are estimated
using an inverse optimization method where both the deviation of geometry from the in vivo
configuration and the deviation of stress from the homeostatic value are minimized. Then,

40

the objective function to minimize is,
∫

∫
W =

k
2
∑ ν i (σ i (h, αk ) − σ i )2 dA
Ω ∫
Ω ||x(h, α ) − Ximage || dA + ξ
h
∫
i )2 dA
¯ 2 dA
Ω ||Ximage − X||
Ω (σh
i

(4.13)

where i = m, 1, ..., k and x is the FE solution for the nodal position vector (see chapter 6
for details of the image-based FE model of the pressurized artery). Ximage constitutes the
target geometry and is the nodal position vector of the FE mesh that is generated based on
the 3-D models of arteries. The geometric model of the artery is reconstructed from MRI
images of a normal aorta and the computational mesh is generated using triangular elements
¯
(see Sheidaei et al. (2011)). X is the geometric center of the artery and is used as a fixed
reference point for normalizing the first term. σ i is a scalar measure of stress in the direction
i
of the constituent i obtained from the FE analysis. σh and ν i are the homeostatic stress

and mass fraction assumed for the constituent i. (h, αk ) are the unknown wall thickness
and anisotropy, i.e. orientation of the collagen fiber k. The deviation of stress for each
constituent is scaled by its mass fraction and integrated over the computational domain.
Therefore, constituents such as fiber families in helical directions will be given more weight
in minimizing the objective function (see Table 4.1). The objective function is composed
of two additive terms and a weight parameter ξ; first term is related to the deviation of
geometry (named “GD” hereafter) and the second term is related to the deviation of stress
(named “SD” hereafter). Since we are interested in minimizing both terms, each term
is normalized and a weight parameter ξ (or a Lagrange multiplier) is used to adjust the
minimization weight for each normalized term.
However, solving this optimization problem for the thickness and anisotropy at all nodal
points of the FE model is not practical, even if possible. Thickness and anisotropy distributions can be approximated with a smaller (I) number of variables with associated base

41

functions, independent from the FE mesh as

h(x, y, z) =

I
∑

h
{βj ϕj (x, y, z)}

(x, y, z) ∈ Ω

j=1

αk (x, y, z) =

I
∑

k
{βj ψj (x, y, z)}

(x, y, z) ∈ Ω,

(4.14)

j=1
h k
where (βj , βj ) are variables for thickness and anisotropy associated with the approximation

point j. ϕj (x, y, z) and ψj (x, y, z) are basis/approximation functions defined on the computational domain Ω. The objective function then can be rewritten with respect to the new
design variables as
∫

∫
W =

h k
2
h k
i
∑ ν i Ω (σ i (βj , βj ) − σh )2 dA
Ω ||x(βj , βj ) − Ximage || dA
∫
∫
+ξ
.
i 2
¯ 2
Ω ||Ximage − X|| dA
Ω (σh ) dA
i

(4.15)

To facilitate the approximation in (4.14), the computational domain (the mid-surface of
the vessel wall) can be parameterized by two spatial variables (s, θ) where s and θ represent,
respectively, the longitudinal distance and azimuthal position on the arterial wall (explained
in next section). Then, equation (4.14) can be rewritten as

h(s, θ) =

I
∑

h
{βj ϕj (s, θ)}

j=1

αk (s, θ) =

I
∑

k
{βj ψj (s, θ)}.

(4.16)

j=1
h
Towards solving the optimization problem (equation 4.15), we use initial values of (βj ,
k
k
βj ) that approximate a homogenous field of thickness and anisotropy (h0 , α0 ). That is, the

42

initial values are obtained by solving the following sets of least-squares optimizations

Sh =

Ne
∑
e=1

Sk =

Ne
∑
e=1






I
∑

2
h
βj ϕj (se , θe ) − h0 

j=1
I
∑

(4.17)

2
k
k
βj ψj (se , θe ) − α0  ,

(4.18)

j=1

where Ne and I are the number of elements and approximation points, respectively.

4.2.3

Parameterizing the aortic wall surface with two variables

A point on the vessel wall can be parameterized by two variables, one that characterizes its
longitudinal position (s) and the other which characterizes its orientation (θ) with respect
to a reference direction. To do so, we need to approximate the centerline of the vessel
considering some of the points on the centerline as nodal points (Figure 4.3) and

X(s) =

∑

Φi (s)Xi ,

(4.19)

i

where Xi and Φi are the position vector and interpolation function corresponding to the
nodal point i on the centerline. X is the position vector of any point on the centerline as a
function of s. A fourth order interpolation function is assumed with the general form of
Φ(s) = c(s − a)2 (s − b)2 .

43

(4.20)

The interpolation functions associated with nodal points j = 1, ..., J can be defined as
(s − L3 )2 (s + L3 )2
L1 ≤ s < L3
(L1 − L3 )2 (L1 + L3 )2
(s − L1 )2 (s − L4 )2
Φ2 (s) =
L 1 ≤ s < L4
(L2 − L1 )2 (L2 − L4 )2
(s − Lk−2 )2 (s − Lk+2 )2
k (s) =
Φ
L
≤ s < Lk+2
(Lk − Lk−2 )2 (Lk − Lk+2 )2 k−2
Φ1 (s) =

ΦJ−1 (s) =
ΦJ (s) =

(4.21)
(4.22)
(4.23)

(s − LJ−3 )2 (s − LJ )2
L
≤ s < LJ
(LJ−1 − LJ−3 )2 (LJ−1 − LJ )2 J−3

(4.24)

(s − LJ−2 )2 (s + LJ−2 )2
(LJ − LJ−2 )2 (LJ + LJ−2 )2

(4.25)

LJ−2 ≤ s < LJ

where k = 3, ..., J − 2 and Lj is the value of s at the nodal point j (Figure 4.3). These
∑
interpolation functions, however, do not satisfy the condition J Φj (s) = 1. In order to
j=1
provide this condition, we need to normalize interpolation functions as
Φi (s)
ˆ
Φi (s) = ∑J
.
Φj (s)

(4.26)

j=1

Now, using the interpolation in (4.19) we can find the parameter s associated with any point
on the artery, e.g. center point of a triangular element on the surface (Xc ). That is, for a
given point on the aortic wall, the variable s is calculated by minimizing the distance from
the point on the wall to the centerline (||X(s) − Xc ||). The function to be minimized is given
as

d(s) =

(∑

ˆ
Φi (s)xi − xc

i

+

(∑

)2

+

(∑

ˆ
Φi (s)y i − yc

)2

i

ˆ
Φi (s)z i − zc

i

44

)2

.

(4.27)

Minimizing d(s) with respect to s results in
(∑
)∑
∂d(s)
ˆ
ˆ ,s
= 2
Φi (s)xi − xc
Φi x i
∂s
i
i
(∑
)∑
ˆ
ˆ
+2
Φi (s)y i − yc
Φi y i
,s
+2

i
(∑

ˆ
Φi (s)z i − zc

i
)∑

i

ˆ
Φi z i = 0.
,s

(4.28)

i

Numerical solution of the nonlinear equation (4.28) is obtained using Newton-Raphson
method which also requires the second derivative of the function. The iterative scheme
for the Newton-Raphson is formulated as
∂d(s)

|s=sn
sn+1 = sn − 2∂s
.
∂ d(s)
2 |s=sn

(4.29)

∂s

This is repeated for any other point of interest on the wall in order to find the corresponding
value of s. If s0 is the solution associated with a center point of an element (Figure 4.3), the
vector v connecting the point on the centerline at s = s0 (X(s = s0 )) and the center point
of the element is given as
v = Xc −

∑

ˆ
Φi (s0 )Xi .

(4.30)

i

The normalized vector n tangent to the centerline at s = s0 is then given by
∂X(s)
∂s |s=s0
n=
∂X(s)
|| ∂s |s=s0 ||

where

∂X(s) ∑ ˆ i i
=
Φ,s X .
∂s

(4.31)

i

The vector n is also a normal vector to the plan of cross section at s = s0 . Projection of an
arbitrary vector a on the plan of cross section (Figure 4.3) can be assumed as the reference
direction
ap = a − (a · n)n.

45

(4.32)

The angle θ between ap and v characterizes the orientation associated with the current point
on the wall (i.e., Xc ). Figure 4.4 illustrates the 3-D geometry of the model of aorta mapped
in 2-D plane of longitudinal (s) and azimuthal (θ) variables.

4.2.4

Global approximation approach

For approximations in equations (4.14), products of Legendre polynomials and periodic functions, respectively, for longitudinal and azimuthal directions are used

h(s, θ) =

m=M −1,n=N −1
∑

h
βmn Pm (s)Fn (θ)

(4.33)

k
βmn Pm (s)Fn (θ),

(4.34)

m,n=0

αk (s, θ) =

m=M −1,n=N −1
∑
m,n=0

where M and N are, respectively, the total number of Legendre polynomials and periodic
functions (i.e. I = M × N ). Pm (s) is a univariate Legendre polynomials of order m such
that P0 (s) = 1, P1 (s) = s and
(
Pm+1 (s) = s

2m + 1
m+1

)

(
Pm (s) −

m
m+1

)
Pm−1 (s).

(4.35)

Also, we consider F0 (θ) = 1 and

F2n−1 = sin(nθ)
F2n = cos(nθ).

4.2.5

(4.36)

Optimization algorithm

The Nelder-Mead Simplex method is employed for the optimization (Lagarias et al., 1998;
Nelder and Mead, 1965). As a direct search method, it does not require gradients of the
function, which is desirable in applications where the calculation of gradients of the function
46

s = LJ

s = LJ−2
n

a
θ
v

X (s = s0)
a
Xc

s = L4

s = L3

s = L2
s = L1
Figure 4.3: Geometry of an arbitrary model of the arterial wall with its centerline; Approximation/nodal points with their associated length s = Lj (j = 1, ..., J). a is an arbitrary
vector used in order to find the orientation θ associated with a point (Xc ) on the wall

47

0.07

Longitudinal Parameter (m)

0.06
0.05
0.04
0.03
0.02
0.01
0
−200

−150

−100

−50

0

50

100

150

200

Azimuthal Parameter (deg)
Figure 4.4: Geometry of the vessel wall parameterized with longitudinal and azimuthal (s
and θ) variables. Dots represent center points of all elements on the wall

48

is computationally expensive. Another feature of the Nelder-Mead Simplex method is the
fast reduction in the objective function after the first few iterations (Wright, 1996). See
Appendix B for details of the algorithm used. A stopping criterion is chosen based on both
the relative size of the simplex and function values at vertices of the simplex as (Torczon,
1989):
1
k
max ||v k − v0 || < δ
∆ 1≤j≤I j

(4.37)

k
k
W (vI ) − W (v0 ) < ϵ ,

(4.38)

k
where vj is the j th vertex of the simplex and a vector comprised of all optimization variables
k
k
at k th iteration. v0 and vI are the “best” and “worst” vertices of the simplex at k th iteration
k
and ∆ = max(1, ||v0 ||).

4.3

RESULTS

As a parametric study, the effect of variation of the weight parameter ξ is first investigated.
Figure 4.5 shows the GD and SD corresponding to minimum values of the objective function obtained with different values of ξ and using two different combinations of Legendre
polynomials and periodic functions (M = 3, N = 3) and (M = 6, N = 5).
In both cases small values of ξ puts more weight on GD to minimize the objective function
and increasing ξ shifts the weight towards SD. The tradeoff choice according to both cases
appears to be ξ = 0.01 such that both parts can be minimized at the same time (Figure 4.5).

4.3.1

Finding the optimal distributions of thickness and anisotropy

In general higher number of approximating polynomials increases the flexibility for better
approximation but at the cost higher computational cost. Then we choose to use six Legendre
polynomials (M = 6) and 3 periodic functions (N = 3) for the approximation assuming

49

10

8

M=3
6 N=3

M=3
N=3

8
6

)
−3

2
0

-3

-2

1

-1

2

Error of Stress ( x 10

−5

Error of Geometry ( x 10 )

4

3

10 10 10 1 10 10 10
8

M=6
6 N=5
4
2
0

4
2
0

-3

-2

1

-1

2

3

10 10 10 1 10 10 10
10

M=6
N=5

8
6
4
2

-3

-2

1

-1

2

0

3

-3

-2

1

-1

2

3

10 10 10 1 10 10 10

10 10 10 1 10 10 10

ξ

ξ

Figure 4.5: The effect of variation of the parameter ξ on both GD and SD using (M =
3, N = 3;top) and (M = 6, N = 5;bottom)

50

10
8
6
0

50

100

150

4
3
2
1
0
0

200

−3

4

Error of Stres (x 10 )

Error of Geometry
−5
(x 10 )

Obj. function
( x 10 −5 )

12

100

200

300

250

300

350

10
8
6
4
0

100

200

300

Iteration

Iteration

Figure 4.6: Changes in the objective function and its associated compartments versus optimization iterations using 36 variables (18 variables for approximating thickness and 18
variables for approximating fiber orientation) considering ξ = 0.01

51

mm
0.5
0.4
0.3
0.2
0.1
0
a

b

Figure 4.7: Deviation of the geometry from the in vivo geometry without (a) and with (b)
optimized distributions of thickness and anisotropy (||x − Ximage ||)
ξ = 0.01. This constitutes 18 variables (I = 18) for thickness and anisotropy, including a
total of 36 variables into the optimization process. Note that fibers oriented in circumferential
and axial directions are considered fixed and only helical fibers orientations are assumed to
be changing (α3 = α4 ). Least-squares estimation of variables associated with a homogenous
field of thickness and anisotropy (e.g., 0.8 mm for thickness and 50.0◦ for anisotropy) yielded
h
k
estimates such as β00 = 0.8, β00 = 50.0 and 0 for all other parameters. Figure 4.6 illustrates

the convergence history of the objective function as well as its compartments, GD and SD,
until the stopping criterion is met. A fast decrease in the objective function during the first
100 iterations is noticeable, which is accompanied by sharp decreases in GD and SD. The
appearance of the plateau regions is associated with the iterations during which searching
the space has not led to a new minimum.
For the sake of comparison, the distributions of thickness are prescribed using the same
method employed by Zeinali-Davarani et al. (2011) and the results were compared with
the current method. Figure 6.4 contrasts the deviation from the in vivo/image geometry
(||x − Ximage ||) using both methods. A significant decrease in the maximum deviation
(about 70%) is achieved using the optimization approach.
k
k
The normalized deviation of stress from the homeostatic value ((σ k − σh )/σh ) in the

52

0.3
0.2
0.1
0
-0.1
b
a
k
k
Figure 4.8: Deviation of the stress ((σ k − σh )/σh ) from the target homeostatic stress in
a helical fiber (k = 3) without (a) and with (b) optimized distributions of thickness and
anisotropy

direction of helical fiber families (k = 3, 4) using both methods are shown in Figures 4.8 and
4.9. For fiber families of both helical directions, the maximum deviations of stress from the
homeostatic value are significantly decreased by 70% using the optimization method.
Figure 4.10 depicts the distributions of wall thickness and anisotropy obtained by the
optimization with ξ = 0.01, M = 6, and N = 3. The resulting spatial variation of anisotropy
is not large although thickness considerably varied especially on the convex and concave
regions with higher values on the concave side and lower values on the convex side.

4.4

DISCUSSION

The existence of the vascular mechanical homeostasis and the subsequent adaptation in
response to mechanical stimuli have been fundamental assumptions in mathematical models
of vascular G&R (Baek et al., 2006, 2007b; Figueroa et al., 2009; Kroon and Holzapfel, 2009b;
Watton and Hill, 2009). There has been a growing interest in using such models on a patientspecific basis (Humphrey and Taylor, 2008; Taylor and Humphrey, 2009). Towards that
goal, image-based arterial geometries have been incorporated into stress-mediated vascular

53

0.4
0.3
0.2
0.1
0
-0.1
b

a

k
k
Figure 4.9: Deviation of the stress ((σ k − σh )/σh ) from the target homeostatic stress in
a helical fiber (k = 4) without (a) and with (b) optimized distributions of thickness and
anisotropy

deg.
57
55
53
51
49
47
45

mm
0.85
0.81
0.77
0.73
0.69
0.65
a

b

Figure 4.10: Distributions of thickness (a) and anisotropy (b) obtained from the optimization
results using ξ = 0.01, M = 6, and N = 3

54

adaptation models (Sheidaei et al., 2011; Zeinali-Davarani et al., 2011). Zeinali-Davarani
et al. (2011) utilized the G&R model itself as an optimization tool to drive the mechanical
state towards the target homeostatic value before the main G&R simulations begun. This
approach, however, alters the in vivo configuration even though it provides a desirable stress
distribution. Rather, the present study provides an optimization technique to minimize both
deviations from the homeostatic stress and the in vivo configuration simultaneously.
Numerous methods have been presented in order to compensate for the lack of information
about stress-free or load-free configurations in patient-specific modeling. Raghavan et al.
(2006) used an optimization technique as an approximate method to find the zero-pressure
geometry assuming consistency of displacement field patterns. Using an inverse elastostatic
method, Lu et al. (2007) were able to determine load-free configuration of an AAA as well as
accurate wall tension in a cerebral aneurysm (Lu et al., 2008). Recently, Zhou and Lu (2009)
used the same inverse technique to estimate the open configuration of vessels. In a different
approach, Gee et al. (2009, 2010) showed the utility of the “modified updated Lagrangian”
method in finding meaningful stress analysis results for complex shapes of aneurysms.
However, all of those studies assumed homogenous distributions of the wall thickness and
anisotropy whereas variation of these parameters can have a great impact on the stress/strain
distribution. Instead of finding the load-free configuration, the presented approach focused
on the in vivo configuration and its associated material and geometric parameters of arteries using an inverse optimization method such that the homeostatic condition was restored
while the deviation of geometry from the original in vivo configuration was minimized. In
a somewhat similar approach, Kroon and Holzapfel (2008a) estimated the distribution of
elastic properties of an inhomogeneous and anisotropic membrane using an inverse optimization method and applied the technique to find material properties of the cerebral aneurysm
(Kroon and Holzapfel, 2008b). They used an element partition method for the robust estimation of properties over the domain. That is, they divided the domain into large sub-domains
and performed the optimization for each sub-domain with homogeneous properties. In the

55

next levels of partitioning, they refined each sub-domain while repeating the estimation
process with updated initial values. Alternatively, in the current method, a global approximation scheme is used in order to reduce the number of unknown variables of optimization
and to facilitate estimation of the inhomogeneous properties in a global fashion. Increasing the number of approximation variable theoretically improves the objective function even
more, but at the cost of more computation time. Deviation of stress from the homeostatic
value in both helical directions was dropped by more than 70%, whereas there was no significant reduction in stress of axial and circumferential fibers (not shown), mainly because
of much lower mass fractions assumed in those directions (See equation 4.15). Results of
the AAA simulations using the optimal material parameters, wall thickness and anisotropy
were generally comparable with Zeinali-Davarani et al. (2011), but more advantageous as
the current method reduced the deviation of geometry from the in vivo configuration before
the G&R process initiated.
Direct validation of the optimal distributions of the wall thickness and fiber orientations
requires more experimental data using animal or human arteries. Nevertheless, the proposed
optimization technique provides a useful initialization step, indispensable to patient-specific
G&R simulations.
In closing, in this chapter a scalar measure of stress was used as a mechanical state
governing the mechanosensitive vascular adaptation (Baek et al., 2006). However, it is
still controversial what quantity is responsible for the mechanical homeostatic state (stress,
strain, material stiffness, or their combination?). Next chapter considers the proposed inverse
method to discriminate among different hypotheses of homeostasis.

56

Chapter 5
TESTING DIFFERENT HYPOTHESES OF VASCULAR
HOMEOSTASIS USING IMAGE-BASED GEOMETRIC
MODELS AND AN INVERSE OPTIMIZATION METHOD

5.1

Introduction

In their seminal paper, Humphrey and Rajagopal (2002) introduced a new theoretical framework, called a constrained mixture model for modeling growth and remodeling (G&R) of soft
tissues. They presented a modeling framework that utilizes ideas from classical mixture and
homogenization theories while avoiding the technical difficulties associated with mixture theory. This allows the model to capture the complexity that occurs during soft tissue G&R such
as deposition of multiple structural components (e.g., fibrous collagen, elastin, and smooth
muscle cells) with different natural configurations and turnover rates, whereas the governing
equation can be solved as if the soft tissue were made of a single constituent. Often not mentioned, but another important contribution of their paper is the philosophy and guidance in
developing a model of G&R of biological tissues, summarized in six remarks in the paper.
Based on their pioneering work, during the past decade, many constrained mixture models
have been developed in the studies of vascular mechanics (Hansen et al., 2009; Gleason and
Humphrey, 2004a; Valent´ et al., 2009a; Wan et al., 2010), progression of cardiovascular
ın
diseases (Baek et al., 2006, 2007b; Figueroa et al., 2009; Gleason and Humphrey, 2004b;
57

Sheidaei et al., 2011) and mechanosensitive cellular behavior (Humphrey et al., 2008; Hsu
et al., 2009).
While most studies after Humphrey and Rajagopal (2002) have focused on modeling the
evolution of biological or engineered tissues after the alteration of physiological and pathological conditions, in this paper we study one of the most fundamental, but often overshadowed,
assumption that soft tissue has an optimal structure during the maintenance or normal G&R
period. Humphrey and Rajagopal (2002) stated in Remark 2.3 that “Normal growth and
remodeling tends to be a stable dynamical process, one that seeks to optimize structure and
function with respect to yet unidentified parameters.” They further stated that “Although
one ultimately seeks parameters that govern the underlying mechanisms of mechanotransduction, it will be sufficient for certain modeling purposes to identify parameters that simply
correlate well with the overall process.”
The question that has long been asked and is yet to be answered in biomechanics research is what are the parameters that correlate well with the G&R process? Are they
stress, strain, or strain rate that cells respond to? Humphrey (2001) claimed that the question may be ill-posed. Stress and strain are merely convenient mathematical concepts and
are not unique observable or physical quantities. He suggested, however, that “the concepts
of stress and strain will continue to be convenient metrics in both empirical correlations
and phenomenological constitutive relations that seek to relate certain cellular responses to
particular stimuli.” A practical question still remains as to what are the constitutive relations for such mechanical quantities, that best describe the maintenance of vascular tissues.
During the maintenance, the turnover of cells and extracellular matrix is balanced and unchanged; hence, there is no net change in mass, structure, or the properties. Therefore, it is
plausible that we may find the right constitutive relation or at least discriminate between different constitutive relations by numerically investigating the consequence of the constitutive
relation in tissue structures and comparing the results with experimental studies. Stress and
strain are interrelated (Humphrey, 2001; Kassab, 2008), and using an idealized symmetric

58

model, such as a straight tube, it may be difficult to distinguish the difference in optimal
structures obtained from different hypotheses.
In this chapter, it is proposed that in medical image-based models different hypotheses
for a mechanical homeostatic state may lead to distinct consequences with regard to the
distributions of wall thickness and anisotropy, giving us the opportunity to better understand
the governing rules of tissue adaptation. Toward this end, an inverse optimization method,
developed in the previous chapter is employed to estimate material and structural parameters
of image-based models of arteries corresponding to a given constitutive relation, the results
were compared with regard to multiple hypotheses of a homeostatic condition.

5.2

Hypotheses on vascular mechanical homeostasis

Vascular tissue tends to adapt in response to changes in its mechanical environment and a
variety of evidence offers diverse hypotheses on this homeostatic tendency. In response to a
sustained pressure increase, the thickness of the blood vessel increases, implying a tendency
toward uniformity of circumferential stress (Matsumoto and Hayashi, 1994; Wolinsky, 1971;
Xu et al., 2000). Blood vessel diameter increases in response to an increase in blood flow,
thereby normalizing the wall shear stress (Kamiya and Togawa, 1980; Langille et al., 1989).
Some studies underscore the importance of strain mechanosensitivity along the coronary
arterial tree (Guo and Kassab, 2004; Lu et al., 2001), whereas others indicate strain rate as
an important component of vascular homeostasis (McKnight and Frangos, 2003). Meanwhile,
the vascular system is under pulsatile forces and several observations emphasize the role of
pulsatility in vascular homeostasis (Cummins et al., 2007; Leung et al., 1976). In an in vivo
study, Eberth et al. (2009) observed a strong correlation between morphology and pulsatility
of pressure and flow rather than mean values, which is consistent with a study by Boutouyrie
et al. (1999).
Stress and strain are both tensorial quantities and previous studies have sought simpler

59

scalar measures of mechanical stimuli for modeling vascular adaptation. For example, Baek
et al. (2006) utilized the magnitude of traction on the plane normal to the fiber direction as
a scalar measure of the fiber stress. In another stress-driven model, Hariton et al. (2007a)
and Driessen et al. (2008) used the normal component of traction on the plane normal to
the fiber direction as the remodeling stimulus. These models assumed that the mechanical
stimuli that vascular cells sense correlate with tension or shear on the fibrous tissue that
they reside on. The stretch of collagen fibers was also used by Watton and Hill (2009) as
the mechanical quantity governing the G&R. Nevertheless, these studies assumed that fibers
have the same mechanical properties and, hence, the uniform stretch of fibers means the
uniform tension.
In a more phenomenological study, Guo and Kassab (2004) examined the distributions
of circumferential stress and strain along the porcine aorta and the coronary arterial tree.
They found that the circumferential stretch ratio (from the zero-stress state to the loaded
state) is relatively uniform compared to the stress and suggested that the vascular system
closely regulates the degree of deformation. It is difficult, however, to explain how the
cells in vascular tissues can regulate the deformation given that cells do not experience
the zero-stress state of the blood vessel in vivo. Lillie and Gosline (2007) suggested that
the strain of elastin during the cardiac cycle is nearly constant along the porcine thoracic
aorta. Our recent study with the porcine thoracic aorta also showed that the cyclic strain
during the cardiac cycle is relatively uniform in the circumferential direction compared to
the circumferential component of stress (Kim and Baek, 2011).
In this chapter, the following three scalar measures are chosen to be associated with the
homeostatic state:
• Case 1: σ k = ||σ k n||
• Case 2: σ k = n · σ k n
sys

• Case 3: λcyc = λ1 /λdias ,
1
60

where σ k and n are the partial stress and the unit vector representing the alignment of
sys
the k th constituent. λ1 and λdias are circumferential stretches at systolic and diastolic
1

pressures. The distributions of thickness and fiber alignment are obtained by an inverse
optimization using the above three cases.

5.3

Constitutive relations and an inverse optimization
method

First, constitutive relations and material parameters of a healthy human artery are chosen based on chapter 4. Then, the distributions of wall thickness and fiber alignment are
estimated using the inverse optimization integrated with the finite element model of two
different inflated blood vessels.
For Cases 1 and 2, the objective function is similar to chapter 4 as
∫
∑ ν i (σ i (h, αk ) − σ i )2 dA
||x(h, αk ) − Ximage ||2 dA
Ω ∫
h
W = Ω∫
,
+ξ
i
¯
||Ximage − X||2 dA
(σh )2 dA
Ω
Ω
i
∫

(5.1)

and the corresponding objective function to minimize in Case 3 is defined as
∫
cyc 2
cyc
k
||x(h, αk ) − Ximage ||2 dA
Ω (λ ∫(h, α ) − λh ) dA
Ω∫
+ξ
W =
,
cyc 2
¯ 2
Ω ||Ximage − X|| dA
Ω (λh ) dA
∫

(5.2)

where λcyc is the cyclic stretch in the circumferential direction. Thickness and anisotropy
distributions are then approximated, using the same technique in chapter 4
I
∑

h
h(s, θ) =
{βj ϕj (s, θ)},
j=1

αk (s, θ) =

I
∑

k
{βj ψj (s, θ)},

j=1

61

(5.3)

Case 1
Case 2

1
0.5

200

400

6
4
2
0
0

500
Iteration

1000

600
Deviation of Stress (x 10 −3 )

0
0

Deviation of Geometry
(x 10−5 )

Objective Function
−4
( x 10 )

1.5

800

1000

10
8
6
4
2
0

0

500
Iteration

1000

Figure 5.1: Change in the objective function and its associated compartments versus optimization iterations in Cases 1 and 2 using 60 variables (M = 6;N = 5;ξ = 0.01).

5.4

Results

Two 3-D geometric models, one from a healthy aorta and the other from a healthy internal
iliac artery, are used as computational domains. For simplicity, fibers oriented in circumferential and axial directions are considered fixed and only helical fiber orientations (α3 , α4 )
are considered as variables of anisotropy such that α3 = −α4 .
Figure 5.1 illustrates the convergence history of the objective function and its compartments, GD and SD, corresponding to Cases 1 and 2 for the aorta model. A sharp decrease
in GD and a similar reduction in SD are noticed for both stress hypotheses when the convergence is achieved. Note that the small plateau regions are associated with the search periods
when a new minimum has not been reached yet. Figure 5.2 contrasts the corresponding op62

timized distributions of wall thickness and anisotropy for the aorta model in Cases 1 and 2.
In both cases, the concave side is found to be relatively thicker than the convex side. Helical
fibers on the convex side tend to orient themselves more in the circumferential direction as
opposed to the concave side, even though the overall anisotropy variation is not as large
as the thickness variation. The similar convergence level for both hypotheses (Figure 5.1)
follows from the similar inhomogeneous distributions of thickness and anisotropy for both
cases. In spite of a similar trend of distribution, a larger circumferential anisotropy variation
is noticeable in Case 2 (Figure 5.2d).
An analogous comparison between Cases 1 and 2 has been made for the internal iliac
artery model (Figure 5.3). In both cases, relatively similar distributions of wall thickness
are observed. However, fiber orientations are distributed with more variability between
the two cases than for the aorta model. For the iliac artery model in Case 1, Figure 5.4
illustrates the geometric deviation from the in vivo configuration (||x − Ximage ||) as well as
k
k
the normalized deviation of stress in a helical fiber from the homeostatic value ((σ k −σh )/σh )

when (M = 5;N = 3) (a, c) and (M = 6;N = 5) (b, d). As expected, the geometric deviation
is minimized (< 0.3 mm) in both conditions (Figure 5.4a,b).
The maximum deviation of fiber stress from the homeostatic value still seems to be
large (< 80%) although on average the deviation has been minimized (Figure 5.4c,d). Some
adaptive provisions in the optimization may be useful to reduce the localized high values
of deviations in future improvements of the proposed technique. Apparently, increasing the
number of approximation points (from I = 15 to I = 30) in the optimization results in
only a minor improvement in reducing the average deviations, while the maximum deviations remains the same. That is, a further increase in approximation points may not be
computationally justifiable.
With regard to Case 3, Figure 5.5 shows the geometric deviation from the original in
vivo configuration (a), the normalized deviation of the circumferential cyclic stretch (b) as
well as the optimized distributions of wall thickness (c) and anisotropy (d) for the aorta

63

(b)

mm
0.88
0.84
0.8
0.76
0.72
0.68

mm
0.88
0.84
0.8
0.76
0.72
0.68

(a)

(c)
56
54
52
50
48
46
(d)

Figure 5.2: Optimized distributions of thickness (a, b) and anisotropy (c, d) for the aorta
model obtained in Cases 1 and 2 (M = 6;N = 5;ξ = 0.01). Black arrows identify the
direction of one set of corresponding collagen fibers.

64

mm
0.46
0.42
0.38
0.34
0.3
0.26

mm
0.46
0.42
0.38
0.34
0.3
0.26

(a)

(b)

deg.

deg.

60
57
54
51
48
45

60
57
54
51
48
45

(c)

(d)

Figure 5.3: Optimized distributions of thickness (a, b) and anisotropy (c, d) for the internal
iliac artery model obtained in Cases 1 and 2 (M = 5;N = 3;ξ = 0.1).

65

mm
0.28
0.24
0.2
0.16
0.12
0.08
0.04

mm
0.28
0.24
0.2
0.16
0.12
0.08
0.04

(a)

(b)

0.8
0.6
0.4
0.2
0
-0.2
-0.4

0.8
0.6
0.4
0.2
0
-0.2
-0.4

(c)

(d)

Figure 5.4: Deviation of the geometry from the in vivo configuration (||x − Ximage ||) and
k
k
deviation of stress in a helical fiber direction from the homeostatic value ((σ k − σh )/σh )
when (M = 5;N = 3) (a, c) and when (M = 6, N = 5) (b, d) for the iliac artery model.

66

model. As required by the objective function, the geometric deviation from the image is
minimized (to a maximum value of 0.05 mm). Cyclic stretch is also homogenized toward
the target homeostatic value (i.e., 1.02) on most parts with the largest deviation on areas
close to the fixed boundaries (maximum normalized deviation of 1.5%). When compared
to Cases 1 and 2 (Figure 5.2), lower levels of thickness are predicted by the cyclic stretch
hypothesis (Figure 5.5c). More circumferential and less longitudinal variations of anisotropy
are noticeable in Case 3 (Figure 5.5), whereas anisotropy is more locally distributed in Cases
1 and 2. Convergence of the optimization is an issue when the cyclic stretch hypothesis is
applied to a more complex geometry such as the iliac artery model (not shown).

5.5

Discussion

In this chapter, a numerical optimization method was utilized to find the optimal distributions of structural parameters of healthy arteries based on different homeostatic assumptions
(i.e., uniform stress and cyclic stretch). It was found that for medical image-based models different homeostatic assumptions lead to different inhomogeneous distributions of wall
thickness and anisotropy. The potential application of the proposed optimization method to
arteries with higher levels of geometric complexity was also illustrated by using the internal
iliac artery model.
Different assumptions for stress homeostasis (Cases 1 and 2) led to homologous overall
distributions of thickness, but differing distributions of anisotropy for both arteries (Figures
5.2 and 5.3). Interestingly, the arterial wall on the concave regions was found to be thicker
than on the convex regions and the variation of anisotropy was small. Possible correlations
between the wall thickness and anisotropy variation and other geometric parameters such as
the surface curvature and tortuosity of the vessel wall can be examined in an independent
study by systemic application of the technique using multiple image-based arterial geometries. Case 3 (stretch-based homeostasis) produced noticeably contrasting results with Cases

67

(a)

(c)

(b)

mm
0.05
0.042
0.034
0.026
0.018
0.01
0.002

mm
0.8
0.76
0.72
0.68
0.64
0.6
0.56
0.52

0.003
-0.001
-0.005
-0.009
-0.013

61
59
57
55
53
51
(d)

Figure 5.5: Deviation of the geometry (a), cyclic stretch in circumferential direction (λc )(b)
and corresponding distributions of thickness (c) and anisotropy (d) for the aorta model
obtained in Case 3 (M = 6;N = 3;ξ = 0.01). Black arrows identify the direction of one set
of corresponding collagen fibers.

68

1 and 2. Indeed, in Cases 3, the wall thickness was predicted to be in general smaller than in
Cases 1 and 2 (stress-based homeostasis), while a less longitudinal anisotropy variation was
observed (Figure 5.5). The assumption of a uniform cyclic stretch resulted in a remarkable
improvement for both terms of the objective function (Figure 5.5a,b), thereby revealing a
stronger potential for the cyclic stretch homeostatic assumption, which remains to be further
validated with experimental data.
Meanwhile, here, a constant proportion of the constituents as well as a uniform distribution of the material parameters were assumed and the more attention was given to the
variation of structural parameters such as the wall thickness and fiber orientation. Some
studies, however, indicated non-uniform deformations of the arterial wall during the cardiac
cycle (Draney et al., 2004) and observed a circumferential variation of material properties
of the aortic wall. In a recent ex vivo inflation test of aorta Kim and Baek (2011) found
that the posterior region of the arterial wall is thinner, but stiffer than the anterior region.
This may be due to the influence of surrounding tissue as well as the local variation in the
content of each constituent. Our study excluded the effect of anatomical variation of boundary condition due to the surrounding tissue which may have influenced our results. This
effect can be more pronounced when a complex in vivo geometry such as the internal iliac
artery (Figure 5.3) is involved. Also, the fixed boundary conditions considered at both ends
are inconsistent with the cyclic stretch homeostatic assumption, thereby resulting in large
deviations at both ends (Figure 5.5b).
Based on the preliminary results, the error or deviation of geometry relative to the image
geometry is less sensitive to the fiber orientations than the wall thickness. Removing the
condition α3 = −α4 (i.e., adding one more set of optimization variables) did not result in
more variability in optimal fiber orientations, whereas it doubled the computational cost.
Considering variable orientations for fibers in the axial and circumferential directions greatly
adds up to the numerical complexity of the optimization.
It is worth mentioning that the mechanical homeostasis of the normal tissue maintenance

69

may also change under various additional conditions, such as aging and vascular diseases.
The application of the current technique to diseased arteries will require more understanding
of tissue growth and remodeling under pathological conditions. Nevertheless, the current
method showed a distinguishing capability in terms of the optimal structure based on various
forms of constitutive relations and homeostatic assumptions. There is a clear need for more
experimental data in order to validate and exploit the current technique in understanding the
multifactorial process of vascular homeostasis and possible implementation in patient-specific
models of vascular diseases.

70

Chapter 6
A FINITE ELEMENT MODEL OF STRESS-MEDIATED
VASCULAR ADAPTATION: APPLICATION TO ABDOMINAL
AORTIC ANEURYSMS

6.1

Introduction

FEM-based vascular G&R simulations typically utilize microstructural information of structural components (collagen, elastin, and smooth muscle) to define the mechanical state at
a given time (either stress or strain), which is iteratively fed back for calculating the evolution of microstructural properties of the components due to mechano-sensitive cellular
activities. Humphrey and coworkers presented FEM models of stress-mediated G&R using a
constrained mixture approach and modeled intracranial aneurysms (Baek et al., 2005, 2006;
Figueroa et al., 2009). Hariton et al. (2007b) studied the stress-driven collagen fibers remodeling in a human carotid bifurcation. On the other hand, strain-mediated models of
vascular adaptation have also been developed and implemented to model cerebral and aortic
aneurysms (Watton et al., 2004; Kroon and Holzapfel, 2007, 2009a; Watton and Hill, 2009).
Driessen et al. (2004, 2008) investigated stretch-mediated remodeling of collagen fibers in
arterial tissues.
Elastin degradation and collagen turnover are main biomechanical processes in enlargement of aneurysms. Abdominal aortic aneurysms (AAAs) have been associated with a
71

decrease in elastin content (Rizzo et al., 1989; He and Roach, 1994; Ghorpade and Baxter, 1996; Powell, 2002) resulting from proteolytic activity. It has also been suggested that
stress(strain)-mediated collagen turnover plays an important role in aneurysm enlargement
and it may differ among individuals depending on physiological and pathological conditions
(Choke et al., 2005). Hence, both elastin degeneration and collagen turnover should be
considered in the computational models of AAAs development. Watton et al. (2004) and
Watton and Hill (2009) were the first to develop a mathematical model of the evolution of
an AAA by accounting for elastin degradation and collagen turnover. In their study, spatial
and temporal alterations of elastin concentration were prescribed and collagen remodeling
acted to maintain the strain in fibers to an equilibrium value.
Previous modeling of evolving vascular diseases have mainly focused on hypothesis-testing
with simple geometries, however, and there is a need to account for patient-specific anatomical and physiological information in FEM simulations of vascular diseases (Taylor and
Humphrey, 2009). In this work, we extend the previous FEM model of stress-mediated
G&R developed by Baek et al. (2005, 2006) to incorporate a non-axisymmetric geometry
obtained from medical images. The computational method is, then, applied to model an
enlarging AAA assuming that an AAA grows by elastin degradation/damage and stressmediated collagen turnover governed by the local state of intramural stress. We contrast
stress distributions and expansion rates of the computationally grown AAAs for multiple
kinetic parameters and different spatial distribution of elastin damage.
Toward this end, we employ a theoretical model of arterial G&R based on a constrained
mixture approach and use a nonlinear FEM with linear triangular elements.

72

6.2
6.2.1

Theoretical modeling of arterial G&R
Kinematics and important configurations

There are two very different timescales in arterial mechanics: a short-term timescale associated with pulsatile motion during the cardiac cycle (seconds) and a long-term timescale for
G&R (days to weeks or months). Following Figueroa et al. (2009), t denotes time associated
with the cardiac cycle, during which we assume no G&R, and s denotes time during G&R.
The configuration κt is the current configuration and x(t) denotes the position of a particle
in κt . Here we make a tacit assumption that there exist a stress-free configuration κsf of
the vessel wall which is fixed at a given s. The vector Xsf (s) denotes the position vector
of the particle in the stress-free configuration and Fsf (t) denotes the deformation gradient
corresponding to the mapping from Xsf (s) to x(t). In arterial mechanics, the effect of the
inertial forces in arterial wall motion are negligible (Humphrey and Na, 2002) and we define
the current (in vivo) configuration κs in the G&R timescale as the configuration κt under
the mean pressure during the cardiac cycle.
For modeling arterial G&R, it appears to be advantageous to introduce additional configurations shown in Figure 6.1 (Baek et al., 2006; Figueroa et al., 2009). The configurations
κτ trace the in vivo configurations of body through time τ ∈ [0, s]. Particularly, the configuration κ0 represents a loaded in vivo configuration of a healthy artery. We assume that each
constituent is pre-stretched when deposited into the tissue at time τ and the tensor Gi (τ )
represents the deposition stretch of constituent i. For computational purposes, we introduce
a configuration κR , called the (computational) reference configuration, that is a fixed configuration. Although the configuration κR is fixed in space, we assume that particles can be
created or removed so that there is one-to-one mapping between κR and κs at time s. This
assumption also implies that the total mass in κR changes with time and becomes always
the same as the total mass in κs . The position of a given particle in κR is denoted by X and
the deformation gradient F(s) is given corresponding to the mapping from κR to κs . The

73

stress-free configuration κsf evolves in the G&R time scale and P is the tensor representing
its evolution. Because the total mass is preserved in the mapping from κR to κs at time s,
the density with respect to the reference configuration, ρR (s), can be calculated by

ρR (s) = J(s)ρ(s),

(6.1)

where J(s) = det [F(s)] and ρ(s) is the mass density in the current configuration at time s.
∑
Let φi (s) denote the mass fraction of constituent i, i.e., i φi (s) = 1. The mass density of
constituent i with respect to the reference configuration, ρi (s), is defined as
R
ρi (s) = φi (s)ρR (s).
R

(6.2)

The deformation gradient for each constituent i at time t, relative to its natural configuration, Fi ) (t), is given as (Baek et al., 2006)
n(τ
Fi ) (t) = F(t)F−1 (τ )Gi (τ ),
n(τ

(6.3)

where F−1 (τ )Gi (τ ) is the tensor representing the pre-stretch of the constituent i that has
been produced at τ with respect to the reference configuration (see Figure 6.1). The deformation gradient Fi ) (t) can also be written as
n(τ
Fi ) (t) = Fsf (t)P(s)F−1 (τ )Gi (τ ),
n(τ

(6.4)

where, now, P(s)F−1 (τ )Gi (τ ) represents the pre-stretch of the constituent i that has been
produced at τ with respect to κsf .

74

2
κ1(0) κ n(0)
n

2
κ 1(τ ) κ n(τ )
n

k
κ n(0)

G 1(τ )

G 2(0)

G 2(τ )
G k (τ )

G k (0)

G 1(0)

k
κ n(τ)

X(s )

κτ

κ0

G&R
κs
F (τ )

F ( 0) = I

F (s )

Fsf (s )

P (s )

X

κR

κ sf

Figure 6.1: Schematic view of configurations involved in the G&R simulation. A fixed
reference configuration is considered for the computation of deformation associated with
G&R time τ = [0, s]. It is chosen to coincide with the configuration of a healthy in vivo
artery at time τ = 0 (i.e., κ0 ). We imagine the existence of a stress-free configuration, κsf ,
associated with the current configuration κs .

75

6.2.2

Stress response and constitutive assumptions

Consider an arterial wall tissue consisting of multiple structural components, e.g., elastin,
multiple collagen families, and smooth muscle (SM). Following previous convention, we use
the subscript ‘i’ to refer each constituent (e.g., elastin, collagen families, muscle) and the
superscript ‘k’ to the k th family of collagen fibers. Thus, i = e, 1, ..., k, ..., m where e denotes
elastin and m denotes SM. Although we adopt the G&R formulation based on the constrained
mixture approach (Humphrey and Rajagopal, 2002), following its later modification (Baek
et al., 2006; Figueroa et al., 2009) we take the point of view of a biological tissue comprising
multiple constituents as a single continuum, where the mass fraction of each constituent and
its pre-stretch are considered as attributes for a particle of the biological tissue. We assume
further that for a given time s the mechanical behavior of the artery can be characterized
by a hyperelastic model and we solve inflation problems with the principle of virtual work.
Hence, let the Cauchy stress be given by (Truesdell and Noll, 1965),

t(t) = 2ρ(t)Fsf (t)

∂Ψ(Csf (t)) T
Fsf (t),
∂Csf (t)

(6.5)

where Csf (t) = FT (t)Fsf (t) and Ψ(Csf (t)) is the stored-energy function (per unit mass).
sf
Since the stress-free configurations κsf change with G&R, we utilize a fixed computational
reference configuration κR . Thus, (6.5) can be rewritten with respect to κR , by using
F(t) = Fsf (t)P(s), (6.1), and the chain rule, as

t(t) =

ˆ
2
∂{ρR (s)Ψ(C(t))} T
F(t)
F (t),
J(t)
∂C(t)

(6.6)

(
)
(
)
ˆ
ˆ
where Ψ C(t) = Ψ PT (s)Csf (t)P(s) = Ψ(Csf (t)). For a membrane model, using J =
J2D h/hR for the thickness h, the membrane stress T can be given as
(
)
˜
∂MR (s)Ψ C2D (t) T
T(t) =
F (t)
F2D (t),
J2D (t) 2D
∂C2D (t)
2

76

(6.7)

(
)
(
)
˜
ˆ
where the areal density MR (s) = hR ρR (s) and Ψ C2D (t) = Ψ C(t) (e.g., see Holzapfel
et al. (2000) for two-dimensional formulation). For simplicity, we omit the subscript ‘2D’
(
)
˜
and define wR (s) = MR (s)Ψ C(t) . Then (6.6) becomes

T(t) =

2
∂w (C(t)) T
F(t) R
F (t).
J(t)
∂C(t)

(6.8)

To simulate arterial G&R, we employ the stored energy equation wR (t; s), the stored energy
(per unit area) due to the deformation of x(t) at a given G&R time s, as (Figueroa et al.,
2009)
}
∫ s
∑{
)
)
i (0)Qi (s)Ψi (Ci
i (τ )q i (s, τ )Ψi (Ci
wR (t; s) =
MR
mR
n(0) (t) +
n(τ ) (t) dτ ,

(6.9)

0

i

(
)
where Ψi Ci ) (t) is the stored energy of constituent i that has been produced Ci ) (t) =
n(τ
n(τ
[ i
]T i
Fn(τ ) (t) Fn(τ ) (t), Qi (s) is the fraction of the constituent i that was present at time 0 and
still remains at time s (i.e., has not yet been removed), mi (τ ) is the mass production rate
R
of the constituent i at time τ per unit reference area, and q i (s, τ ) is its survival function,
that is, the fraction produced at time τ that remains at time s.
We employ constitutive relations for deposition tensors and stored energy functions for
constituents used by Baek et al. (2006, 2007b) and Figueroa et al. (2009). The mechanical
property and the “deposition stretch”, named Gc , of the newly synthesized collagen fibers
h
are assumed to be always the same. Let mk (τ ) be the unit vector in the direction of the
k th collagen fiber produced at time τ . The angle between mk (τ ) and the first principal
direction at time τ is denoted by αk (τ ). We can find a unit vector Mk (τ ) in the reference
configuration that corresponds to mk (τ ), i.e.
Mk (τ ) =

F−1 (τ )mk (τ )
.
|F−1 (τ )mk (τ )|

(6.10)

k
If the angle between Mk (τ ) and the first axis of the coordinate system is denoted by αR (τ ),

77

the stretch (of the fiber produced at time τ ) in the fiber direction from the reference to the
current configuration is given by
λk (t) =

√
F(t)Mk (τ ) · F(t)Mk (τ ) .

(6.11)

The stretch of the k th fiber family from its natural to the current configuration can be, then,
calculated as (Baek et al., 2006)
λk (t)
λk ) (t) = Gc k ,
h
n(τ
λ (τ )

(6.12)

where λk (τ ) is stretch of the unit vector in the k th fiber direction from the reference configuration to the configuration at time τ . A similar approach can be adopted for SM, which
are oriented primarily in the circumferential direction. Thus, the stretch of SM is given
by λm ) (t) = Gm λ2 (t)/λ2 (τ ), where Gm is the homeostatic stretch for the SM. The main
h
h
n(τ
structure of cross-linked elastin is formed at early stages of development, thus it is difficult
to trace its production time and to specify F−1 (τ ) in (6.3). So we define a new tensor
˜
Ge = F−1 (τ )Ge , which represents a mapping from the natural configuration of elastin to
h
˜
˜
the computational reference configuration. Then, Fe = F(t)Ge and Ge is postulated as
n
}
{
˜
Ge = diag Ge , Ge , Ge1 e .
1
2
G
1 2

The stored energy functions for the elastin-dominated amorphous, collagen fiber families
(Ψk = Ψc ) and passive SM are given as




c1  e
1
e
Cn[11] + Cn[22] + e
− 3 ,
e
e
2
Cn[11] Cn[22] − Cn[12] 2
}
)2 ]
c2 { [ ( k
−1 ,
Ψc (λk ) (t)) =
exp c3 λn(τ ) (t)2 − 1
n(τ
4c3
}
)2 ]
c4 { [ ( m
−1 ,
exp c5 λn(τ ) (t)2 − 1
Ψm (λm ) (t)) =
n(τ
4c5
Ψe (Ce (t)) =
n

(6.13)
(6.14)
(6.15)

[ ]T
e
e
e
where Cn[11] , Cn[22] and Cn[12] are components of Ce = Fe Fe . Although the constitutive
n
n
n
78

form is the same for collagen fiber families and SM, SM is much more compliant than collagen
fiber families and has less contribution to the passive mechanical behavior of the wall (Burton,
1954). We use a potential function for the active tone of vascular SM as (Baek et al., 2007b;
Zulliger et al., 2004b)
(
)3
S{
1 λM − λ2 (t) }
Ψm (t) =
λ (t) +
,
act
ρ 2
3 (λM − λo )2

(6.16)

where λM and λo are stretches at which the active force generation is maximum and zero, λ2
is the stretch in circumferential direction at time t, and S is the stress at the maximum conm
traction. Then, the total membrane strain energy becomes wR = wR(passive) + MR (t)Ψm .
act

6.2.3

Stress-mediated G&R

In arteries, constituents can be continuously produced and removed and normal tissue maintains its mass and configuration by balanced turnover of constituents. The rates of production and removal change from their balanced normal (basal) values in response to changes in
mechanical environment. Here, we assume that me (s) = 0 and the rates of mass production
R
for collagen fibers and SM are functions of a scalar measure of intermural stress, given by
(Baek et al., 2006)
c
)
MR (s) ( k k
c
Kg (σ (s) − σh ) + mk
basal
c
MR (0)
m
(
)
m (s) = MR (s) K m (σ m (s) − σ m ) + mm
mR
g
h
basal
m
MR (0)

mk (s) =
R

(6.17)
(6.18)

c
m
where MR (0) and MR (0) are the mass of collagen and SM per reference area of a healthy
i
artery at time 0, respectively. Kg (i = 1, 2, ..., k, ..., m) is a scalar parameter that controls

the stress-mediated growth, mi
basal is a basal rate of mass production for the constituent i,
and
σ k (s) =

||Tc (s)mk (s)||
hc (s)

σ m (s) =

79

||Tm (s)mm (s)||
.
hm (s)

(6.19)

where Tc (s) and Tm (s) are the Cauchy membrane stress contributed by collagen and SM
∑
at time s (i.e. Tc (s) = k Tk (s)). hc (s) and hm (s) are contributions of collagen and SM
to the total thickness at time s. Also, let


( s
)



∫ i

 exp − k (˜)d˜


q τ τ

τ
q i (s, τ ) =






 0



s − τ ≤ ai
max
,

(6.20)

s − τ > ai
max

i τ
where kq (˜) can be a function of circumferential stress, wall shear stress, or other state

variables. ai
max is the maximum life span of the constituent i.
The new collagen is deposited with a preferred alignment. Following Baek et al. (2006),
we assume that the alignment of the newly produced collagen is influenced by the orientation
of the existing collagen and it consequently aligns along the direction of the existing collagen
family.

6.3

Computational considerations

6.3.1

Local Cartesian coordinate system

We assume that the wall is a thin membrane, with X = {X1 , X2 , X3 } and x = {x1 , x2 , x3 }
being the reference and current positions in the global Cartesian coordinate system with
base vectors {E1 , E2 , E3 } (Figure 6.2). We use linear triangular elements for developing
a nonlinear FEM model of a non-axisymmetric cylindrical membrane and define a local
Cartesian coordinate system for each element in order to facilitate calculation of the local
deformation gradient and prescribe material anisotropy. The centroid of an element Xc =
c
c
c
{X1 , X2 , X3 } becomes the origin of the local Cartesian coordinate system and ‘element-

wise’ orthogonal surface coordinates ξ1 and ξ2 are allocated to any point on the element
with respect to the origin (Figure 6.2).
80

Current Configuration

Reference Configuration

X( 2)

( 2)

X

Ne

e
E1

ξ1

Xc
(3)

X

Ee
2

ξ1

X(1)

X(1)

X(3)

ξ2

ξ2
E3
E1
E2
Figure 6.2: Orthonormal base vectors for a global Cartesian coordinate system are
{E1 , E2 , E3 }. Those for a local Cartesian coordinate system are Ee , Ee (associated with ξ1 ,
1
2
ξ2 ), and Ee = Ne at a point set to be the center point of a linear triangular element tangent
3
to the membrane at Xc . Ee and Ee are orthonormal bases on the tangent plane and Ne is
1
2
the outward unit normal vector.

81

A linear triangular element has three nodal points {X(1) , X(2) , X(3) } in the reference and
(k)
1∑
c
{x(1) , x(2) , x(3) } in the current configurations and Xi = 3 k Xi
(i, k = 1, 2, 3). Now,
we introduce local bases Ee and Ee in the plane of a linear triangular element (i.e. tangent
1
2
plane) and Ee ≡ Ne the outward normal to that plane (Figure 6.2). The base vector Ee
3
1
is set to a unit vector aligned toward the axial direction of a model of the artery (a unit
tangent vector to ξ1 at the element center point). The outward normal vector Ne is readily
obtained by cross product of vectors connecting nodal points. The base vector Ee is then
2
e
e
e
calculated as Ee = Ne × Ee . If Xe = {X1 , X2 , X3 } and xe = {xe , xe , xe } are the reference
2
1
1 2 3

and current position vectors in the local Cartesian coordinate system, the two-dimensional
right Cauchy-Green deformation tensor is calculated by (Kyriacou et al., 1996; Park and
Youn, 1998)
C=

( ∂xe ∂xe )
·
Ee ⊗ Ee
α
β
∂ξα ∂ξβ

(6.21)

where α, β = 1, 2.

6.3.2

Finite element formulation

The weak form for the membrane can be obtained from the principle of virtual work (Kyriacou et al., 1996),

∫
δI =
S

∫
δwdA −

s

P n · δxda = 0.

(6.22)

In a FEM model, we seek an approximate solution to (6.22). Let a finite approximation of
the current position be
x = Φxp ,

p

xi = ΦiA xA ,

(6.23)

where xp and Φ are the nodal vector for the current position and shape function matrix,
respectively. Then, the governing equation for an element can be obtained from (6.22) as
{ }e
F P =

∫

)
( ∂w ∂C
αβ
˜i ΦiP dA = 0,
p −P
S e ∂Cαβ ∂xP

82

(6.24)

where for a linear triangular element,
( (2)
(1) )( (3)
(1) )
ϵijk xj − xj
xk − xk
˜
.
Pi = P
( (2)
(1) )( (3)
(1) )
||ϵlmn Xm − Xm Xn − Xn ||

(6.25)

We use the Newton-Raphson method to solve (6.24) and the tangent matrix can be given by
[ ]
K

PQ

=

[ ∂F ]e

∫

(

∂Cαβ ∂Cγω
∂ 2w
p
p
p PQ
∂x
S e ∂Cαβ ∂Cγω ∂xP ∂xQ
)
∂w ∂ 2 Cαβ
˜
+
− ΦiP Pi,Q dA
∂Cαβ ∂xp ∂xp
P
Q
=

(6.26)

Note that (i, j, k, l, m, n = 1, 2, 3), (α, β, γ, ω = 1, 2), and (A, B, M, P, Q = 1, 2, 3, ..., np ),
where np is the number of nodes in an element multiplied by 3, i.e., np = 9.
Any point X = {X1 , X2 , X3 } in the global Cartesian coordinate system can be transe
e
e
formed to Xe = {X1 , X2 , X3 } in the local coordinate system using

e
c
Xi = Ee · (X − Xc ) = Qij (Xj − Xj )
i

(6.27)

where Qij = Ee · Ej . We use linear triangular element which is flat and has 9 translational
i
degrees of freedom to discretize the domain (Flores and Onate, 2001). Now, we define linear
interpolation function for the triangular element with respect to ith nodal point as
ϕi (ξ1 , ξ2 ) = ai ξ1 + bi ξ2 + ci ,

(i = 1, 2, 3)

(6.28)

∑(all node) i
where ai , bi , and ci are constants and i
ϕ (ξ1 , ξ2 ) = 1. The interpolation functions
have the following properties:
e(j)

ϕi (X1

e(j)

, X2

) = δij ,

e(j)

where Xk

83

= Ee · (X(j) − Xc ).
k

(6.29)

Thus, ai , bi , and ci can be obtained by solving (6.29), for example,
 

−1  
e(1)
e(1)
 1 
 
 a 
X2
1   1 
  X1

 


 
 

  
=  X e(2) X e(2) 1 
b1
  0 .
2
  1


 
  

 
 
 1 
e(3)
e(3)
 c 

 0 
X1
X2
1

(6.30)

It is convenient to define 9×1 vectors for the nodal points as
e(1)

Xep = {X1

e(1)

xep = {x1

(1)

e(1)

, X2

e(1)

, x2

(1)

e(1)

, X3

e(1)

, x3

e(2)

, X1
e(2)

, x1

(1)

e(2)

, X2

e(2)

, x2

(2)

e(2)

, x3

(2)

e(2)

, X3

e(3)

, X1

e(3)

, x1

(2)

(3)

(3)

(3)

e(3)

, X2

e(3)

, x2

(3)

e(3) T
}

, X3

e(3) T
}

, x3

(3)

Xp = {X1 , X2 , X3 , X1 , X2 , X3 , X1 , X2 , X3 }T
(1)

(1)

(1)

(2)

(2)

(2)

(3)

(6.31)
(6.32)
(6.33)

xp = {x1 , x2 , x3 , x1 , x2 , x3 , x1 , x2 , x3 }T

(6.34)

c
c
c
c
c
c
c
c
c
Xc = {X1 , X2 , X3 , X1 , X2 , X3 , X1 , X2 , X3 }T

(6.35)

The coordinate transformation between global and local coordinates can be expressed by:
˜
Xep = Q(Xp − Xc ),

˜
xep = Q(xp − Xc )

(6.36)

˜
where a 9×9 matrix Q is given by:




 Q 0 0 


˜
Q =  0 Q 0 .




0 0 Q

(6.37)

Note that the local position vectors of element nodes in the current configuration, xep , are
expressed with respect to the local coordinate system associated with that element in the
reference configuration. The position vector in the local coordinate system is
xe (ξ1 , ξ2 ) = Φxep
84

(6.38)

where




ϕ1



Φ= 0


0

0

0

ϕ2

0

0

ϕ3

ϕ1

0

0

ϕ2

0

0

0

ϕ1

0

0

ϕ2

0

0

0 

ϕ3 0  .


3
0 ϕ

and ϕ1 , ϕ2 , and ϕ3 are linear shape functions of ξ1 and ξ2 for the element. Using (6.36),
(6.38) can be rewritten as
[ e
p
c ]
˜
xi = ΦiA QAB (xB − XB )

˜
xe = ΦQ(xp − Xc ),

(6.39)

For convenience, let us define Φe as
˜
Φe = ΦiA QAB .
iB

(6.40)

The components of the derivative of the local position vector with respect to the local
coordinate are
p

˜
xe = ΦiA,α QAB xB ,
i,α

p

or xe = Φe xA
i,α
iA,α

(6.41)

Thus, components of the right Cauchy-Green tensor (6.21) can be obtained as
p

p

Cαβ = xe xe = Φe xA Φe
i,α i,β
iA,α
iM,β xM ,

85

(6.42)

where the basis is Ee ⊗ Ee . Also each component in (6.26) can be obtained by
α
β
∂Cαβ
p
∂xP

∂ 2 Cαβ
p

p

∂xP ∂xQ
∂Cαβ ∂Cγω
∂ 2w
∂Cαβ ∂Cγω ∂xp ∂xp
P
Q

= xe Φe + xe Φe
i,α iP,β
i,β iP,α

(6.43)

e
e
= Φe Φ e
iP,α iQ,β + ΦiP,β ΦiQ,α

(6.44)

=

∂ 2w
xe xe Φe Φe
∂Cαβ ∂Cγω i,β j,ω iP,α jQ,γ
+

∂ 2w
xe xe Φe Φe
∂Cαβ ∂Cγω i,β j,γ iP,α jQ,ω

+

∂ 2w
xe xe Φe Φe
∂Cαβ ∂Cγω i,α j,ω iP,β jQ,γ

+

∂ 2w
xe xe Φe Φe
∂Cαβ ∂Cγω i,α j,γ iP,β jQ,ω

(6.45)

)
∂w ∂ 2 Cαβ
∂w ( e
e
e
ΦiP,α Φe
p
p =
iQ,β + ΦiP,β ΦiQ,α .
∂Cαβ ∂xP ∂xQ
∂Cαβ

(6.46)

For the linear triangle element, the following relation can be obtained:

˜
P=

P
|(X(2) − X(1) ) × (X(3) − X(1) )|

(x(2) − x(1) ) × (x(3) − x(1) ).

(6.47)

{
}
˜
∂ (x(2) − x(1) ) × (x(3) − x(1) )
∂P
P
=
∂xp
∂xp
|(X(2) − X(1) ) × (X(3) − X(1) )|

(3)
(2)
(3)
(2)
0
x3 − x3
x2 − x2

P  (3)
(2)
(3)
 x − x(2)
=
0
x1 − x1
3
|(·)|  3

(2)
(3)
(3)
(2)
x2 − x2
x1 − x1
0
(3)
(1)
x3 − x3

0
(1)

(3)

(3)

(1)

x3 − x3
x2 − x2

0
(1)

(3)

x1 − x1

(1)
(3)
x2 − x2
(3)
(1)
x1 − x1

(1)
(2)
x3 − x3

x3 − x3

0

x2 − x2

0
(1)

(1)

86

(2)

(2)

0
(2)

(1)

x1 − x1

(2)
(1)
x2 − x2
(1)
(2)
x1 − x1

0

(6.48)

(6.49)



.



6.3.3

Numerical solutions

The spatial integration in (6.24) is approximated by using Gauss integration. At every Gauss
point, the temporal integrations in (6.22) are calculated by using the trapezoidal rule (Press
et al., 1992) over past times. Although the survival function (6.20) is given by an exponential
i
function, we define a maximum lifetime ai
max and truncate the value if s − τ ≤ amax . Thus,

the numerical integration can be done over the time interval [s − ai , s] using fixed number
max
of discretization points. Equations (6.17), (6.18), and (6.24) are solved iteratively for the
nodal positions and rates of mass production at a given time s. Briefly, an initial guess of
mi (s) is made at each Gauss point based on the previous time step, and (6.24) is solved for
R
the current positions using the nonlinear FEM technique prescribed in the previous section.
Then, mi (s) is updated using the FEM solution along with (6.17) and (6.18). In this way,
R
rates of mass production and displacements at the current time s are updated iteratively
until solutions converge to the prescribed tolerance.
For testing the utility of the present work, we apply the computational model for simulation of an AAA. One Gauss point is used for the integration. It has been found that
the nonlinear FEM analysis does not always converge well with (6.26). To remedy this, we
modify the tangent matrix by
[ ]
K

PQ

∫

∂Cαβ ∂Cγω
∂w ∂ 2 Cαβ )
∂ 2w
dA
+
p
p
∂Cαβ ∂xp ∂xp
S e ∂Cαβ ∂Cγω ∂xP ∂xQ
P
Q
∫
˜
−
ΦiP Pi,Q dA,
(

= ζ

Se

where ζ = 2.0 shows a good convergence for all simulations in this work.

87

(6.50)

6.4
6.4.1

Simulation of an AAA growth
A geometric model and mesh generation

A 3D computational geometry of a healthy aorta is reconstructed from magnetic resonance
(MR) images of a healthy subject using Simvascular (Cardiovascular lab, Stanford University). The computational domain is extended to the upper part of abdominal aorta (proximal
side) and iliac branches (distal side) for future use in hemodynamic simulations (Sheidaei
et al., 2011). For simulation of the aneurysm, we use only the central region of the geometric
model which is separately meshed with triangular elements (1584 elements) using Gambit
(Lebanon, NH, USA).

6.4.2

Initialization for G&R simulation

In normal physiological conditions, production and removal of each constituent are balanced
such that the vessel maintains its shape under a preferred homeostatic state. For an idealized model of the blood vessel, the in vivo material properties are typically assumed to
be constant over the domain with parameter values estimated from experimental data while
the homeostatic assumption is imposed as a constraint (Baek et al., 2007b; Figueroa et al.,
2009; Valent´ et al., 2009b). When medical image-based geometric models are used, howın
ever, it is not a trivial task to prescribe the distribution of material and geometric parameters
while maintaining the homeostatic mechanical state under physiological pressure. For an approximation, we first prescribe material constitutive parameters shown in Table 6.1. Other
parameters, c1 , c2 , and c4 , are calculated inversely assuming that all constituents have the
same homeostatic stress value (i.e., σ i = σh ) in a healthy state. For that, four discrete
fiber families are assumed to be initially aligned in 0, 90, 45, and -45 (axial, circumferential,
and helical directions, respectively) and constant volume fractions for all constituents are

88

Table 6.1: Constitutive and kinetic parameters for each constituents used in initialization
and G&R simulations.

Elastin: c1 = 112 N m/kg, Ge = 1.25, Ge = 1.25
1
2
c = 1.07, k c = 0.02
Collagen: c2 = 917 N m/kg, c3 = 25, Gh
q
m
SM: c4 = 27 N m/kg, c5 = 8.5, Gm = 1.2, kq = 0.02,
h
S = 50 kP a, λM = 1.2, λ0 = 0.7

Others:

c
m
ρ = 1050 kg/m3 , σh = σh = 135 kP a

prescribed. Then, using strain energy functions (6.13)-(6.16), c1 , c2 , and c4 are obtained by
{ (
)}−1
e 2 − (Ge Ge )−2
c1 = σh ρ G1
,
1 2
{
}−1
4
{
} ∑ Mk
R sin2 αk
c2 = σh ρGc 2 (Gc 2 − 1) exp c3 (Gc 2 − 1)2
,
h
h
h
c
MR

(6.51)
(6.52)

k=1

c4 =

λ −1
σh − S{1 − ( λ M−λ )2 }
0
M
{
}.
ρGm 2 (Gm 2 − 1) exp c5 (Gm 2 − 1)2
h
h
h

(6.53)

Regional wall thickness is initially approximated by estimating the cross sectional mean
radius and using the Laplace equation, σ = P r/h. Then, using the G&R simulation, we let
the vessel wall adapt to an equilibrium state with a high value of stress-mediated parameters
k
m
(Kg = Kg = 1). Deviation of stress in three families of collagen fibers (i.e., axial and helical
c
c
directions) from the desired homeostatic stress ((σ k − σh )/σh ) is plotted before and after

300 “days” of this initial simulation in Figure 6.3. Stress deviation of the fiber family along
the circumferential direction was less than the other fibers (not shown). Figures 6.4(a) and
(b), respectively, depict the displacement and the change in the geometry from the reference
configuration after 300 days of this initial simulation. The stress deviation of each fiber
family is reduced substantially after the initialization process while causing little change in
the geometry from its original/reference configuration (maximum displacement is about 1.5
mm).

89

(a)

(d)

(b)

(e)

(c)

0.25
0.15
0.05
-0.05
-0.15
-0.25

(f)

c
c
Figure 6.3: Deviation of the intramural stress from the homeostatic value ((σ k − σh )/σh ) in
axial and two helical fiber directions before (a, b, and c) and after (d, e, and f) 300 days of
the initial simulation.

90

After
Before

(a)

mm

(b)

1.6
1.4
1.2
1
0.8
0.6
0.4
0.2

Figure 6.4: Magnitude of the displacement (mm) calculated from the reference configuration
after 300 days of the initial simulation (a) and the corresponding geometry before and after
300 days (b).

6.4.3

Simulation Cases

After initialization, an AAA is initiated by introducing damage to the aorta, where elastin
is degraded/removed from the blood vessel. The spatial and temporal function for elastin
damage is given as

(
)
D(X, s) = 1 − 0.01s/T f (X),

(6.54)

where T is a time constant and f (X) defines the spatial distribution of damage. D(X, s) ∈
[0, 1] is the ratio of the degenerated elastin to the initial amount, where D(X, s) = 1 for
complete degradation. First, we assume a limiting case when T → 0+ such that D(X, s)
becomes only a spatial function, f (X). That is, the damage (elastin removal) is applied immediately after the initial simulation for a homeostatic condition and the amount of damage
is kept constant during G&R. Four spatial functions are used for representing different cases
of elastin damage. The areal density of elastin after applying different damages is plotted
in Figure 6.5. Cases (1) and (2) correspond to damage shapes distributed to a relatively
large area on the concave and convex sides, respectively. Cases (3) and (4) correspond to
more localized damage shapes applied at multiple locations. Second, the ‘time dependent’
elastin degradation is investigated using (6.54) with a spatial function of Case (3) considering

91

T = 10 × 365 (days). For final simulation, we assume ‘stretch-induced’ degradation, where
additional removal of elastin occurs during aneurysm growth due to stretch-induced damage.
Specifically, the function for elastin damage can be written as,
]
( e
)[
D(X, s) = f (X) + g I1 (X, s) 1 − f (X) ,

(6.55)

e
where I1 (X, s) is the first invariant of the Ce and
n



 1



))
( (


e
π 7.48−I1
e
g(I1 ) =
1 − sin
2×1.48






 0

e
(I1 − 3) ≥ 4.48
e
3 ≤ (I1 − 3) < 4.48 .

(6.56)

e
(I1 − 3) < 3

It is believed that elastin plays an important role in controlling SM cells migration/proliferation
(Karnik et al., 2003; Li et al., 1998) as well as their phenotype modulation (Ailawadi et al.,
2009) by stabilizing extracellular matrix. The amount of SM in our model is changed proportionally with the initial elastin degradation which reduces both passive and active contributions of SM to the wall mechanical properties.

6.4.4

Results

Distribution of the maximum principal stress and the areal mass density of collagen at 50,
1200, and 2700 days of G&R for Case (3) are shown in Figure 6.6. The damage to elastin
introduced at s = 0 causes a sharp increase in stress at the location of the damage. Stressmediated collagen production increases the areal density of collagen as the lesion enlarged,
compensating the loss of elastin, and the value of the peak stress initially starts to decrease.
As the aneurysm enlarges further, however, the stress level is shown to increase again (see
Figure 6.6(c)). The shape of aneurysm and stress distributions resulted from Cases (1), (2),
and (4) are plotted at 50 and 2700 days of the G&R simulation in Figure 6.7. The simulation

92

Case1
0.3
0.26
0.22
0.18
0.14
0.1
0.06
0.02

Case 3

Case 4

Case2

Figure 6.5: Areal mass density of elastin (kg/m2 ) for different simulation Cases (1)-(4).
Cases (1) and (2) correspond to different damages shapes distributed to relatively large area
which are applied at different locations (concave and convex sides) and Cases (3) and (4)
correspond to more localized damage shapes, but at multiple locations.

93

suggests that the aneurysmal wall tends to enlarge more in the convex side of the vessel than
the concave side.
For example, Figure 6.7(d) shows a large amount of dilatation in the convex side even
with the initial damage introduced in the concave side. Different spacial functions for elastin
damage result in variety of shapes with different expansion rates although in all cases the
same kinetic parameter (Kg = 0.05) was used. Figure 6.8 plots expansion rates (maximum
diameter increase per unit time) for simulations with four different cases. Cases (1) and (2)
result in higher expansion rates compared to Cases (3) and (4), which might be due to the
wider area of the damage prescribed for Cases (1) and (2). The effect of kinetic parameter
(Kg =0.02, 0.05, and 0.1) on the expansion rate of the aneurysm is shown for Case (2) in
Figure 6.9. Apparently, simulation with Kg =0.02 results in an increase in the aneurysm
expansion rate, simulation with Kg = 0.05 causes a linear enlargement over time (constant
expansion rate), and a higher value of the kinetic parameter (Kg =0.1) even stabilizes the
aneurysm growth.
Figure 6.10 shows the areal density of elastin and the maximum principal stress for the
‘time dependent’ case at 50 and 2700 days of the G&R simulation. At 50 days, changes in
stress are small because of the small amount of elastin degradation (Figure 6.10(a)). The
stress increases gradually as the amount of elastin degradation increases and at 2700 days the
maximum principal stress reaches a similar value as in the Case (3) at 2700 days (see Figures
6.6(c) and 6.10(d)). For the last simulation, the damage shape is initially the same as in Case
(3) (Figure 6.11(a)) but elastin is degraded further with the ‘stretch-induced’ elastin damage.
The simulation shows that the affected area gradually increases with stretch-induced damage
(Figure 6.11(b)) and, hence, the principal stress increases to a level higher than Case (3) at
2700 days without the stretch-induced damage (see Figures 6.6(c) and 6.11(d)).

94

500
450
400
350
300
250
200
150
100

(a)

1.65

(d)

1.35
1.05
0.75
0.45
0.15

(b)

(e)

(c)

(f)

Figure 6.6: Maximum principal stress distribution (kP a) for the simulation Case (3) at 50,
1200, and 2700 days after the initial damage (a, b, and c) and the corresponding areal mass
density of collagen (d, e, and f) (kg/m2 ).

95

500
450
400
350
300
250
200
150
100

(d )

(a)

(b)

(e)
(c)

(f)

Figure 6.7: Maximum principal stress (kP a) for simulation Cases (1), (2), and (4) after 50
(a, b, and c) and 2700 (d, e, and f) days.

96

Max. Diameter (mm)

35

Case 1: y = 2.5x + 11
Case 2: y = 2x + 12
Case 3: y = 1.3x + 13
Case 4: y = 0.9x + 14

30
25
20
15
10

2

4
6
Time (Years)

8

Figure 6.8: Expansion rates of the lesion in simulation Cases (1)-(4) as the maximum diameter of the lesion versus time and their best linear fit (Kg = 0.05).

6.5

Discussion

In this work, a FEM model of vascular adaptation is presented and applied to model the
enlargement of an AAA without considering thrombus formation. We use a membrane model
of arterial wall. Although a membrane model has some limitations, it is still preferable in
many G&R simulations (Gleason et al., 2004; Baek et al., 2006; Watton and Hill, 2009).
Watton and Hill (2009) suggested that the ratio of thickness of the wall to the diameter of the
aneurysm decreases as an AAA enlarges and also the remodeling process tends to naturally
maintain a uniform strain (or stress) field through the thickness. Hence, a membrane model
suits to model the deformation of the abdominal aorta and the development of an aneurysm
at physiological pressure. Moreover, a full 3D model of vascular adaptation needs more
information about the variation of constituent properties through the thickness and their
evolution during the vascular adaptation.
Following Watton et al. (2004) and Watton and Hill (2009), elastin degradation was assumed to initiate our G&R simulations and a similar form was used for spatial and temporal
97

28

Max. Diameter (mm)

26

K g = 0.02
K g = 0.05
K g = 0.1

24
22
20
18
16
14
1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

Time (Years)
Figure 6.9: Expansion rates of the lesion in simulation Case (2) with different values of the
kinetic parameter (Kg = 0.02, 0.05 and 0.1).

98

500
450
400
350
300
250
200
150
100
(c)

0.3
0.26
0.22
0.18
0.14
0.1
0.06
0.02
(a)
(b)

(d)

Figure 6.10: Areal mass density of elastin (kg/m2 ) for the case of ‘time-dependent’ degradation after 50 and 2700 days (a, b) and the corresponding maximum principal stress distribution (c, d) (kP a).

99

500
450
400
350
300
250
200
150
100

0.3
0.26
0.22
0.18
0.14
0.1
0.06
0.02
(a)

(c)

(d)

(b)

Figure 6.11: Areal mass density of elastin (kg/m2 ) for the case of ‘stretch-induced’ degradation after 50 and 2700 days (a, b) and the corresponding maximum principal stress distribution (c, d) (kP a).

100

distribution of degradation. We considered, however, multiple spatial functions for elastin
degradation which resulted in variety of aneurysm shapes with different stress distributions
and growth rates indicating the potential clinical application of computational simulation
of G&R with realistic geometries. Consistent with Watton and Hill (2009), it appears that
collagen production tends to compensate for the loss of elastin. The computations suggest,
however, that as a lesion evolves into a more complex shape, stress may be elevated at a
location which was not necessarily a damage site (see Figure 6.6). In Cases 1 and 2, stress
initially increased at locations where damage was introduced. Later, in the course of enlargement, the convex side of both aneurysms became the region of maximum stress leading
to more dilation on this side. The results suggest that the location and the geometry of
the lesion can influence AAA enlargement and rupture. For an AAA in vivo, however, the
effect of geometry may be through its influence on both hemodynamics (e.g., wall shear
stress; Hoi et al. (2004)) and wall stress distribution (Doyle et al., 2009; Vorp et al., 1998).
Asymmetric flow patterns have been associated with formation of aneurysm in lower limb
amputees via asymmetric distribution of wall shear stress (Vollmar et al., 1989; Paes et al.,
1990; Naschitz and Lenger, 2008). As expected, in our simulations, damage shapes which
were more dispersed resulted in more dilation than localized damage shapes (see Figure 6.7).
The ‘time dependent’ simulation for Case (3) did not result in a substantial difference in the
aneurysm shape, but different time constants showed significant effects on the expansion rate
(Figure 6.12). Although we used simple spatial and temporal functions for elastin degradation, elastin degradation during AAA growth involves multiple biological and mechanical
parameters. It has been suggested that in AAAs, elastin degradation is due to the proteolytic activity which may have several causes including abnormal distribution of wall shear
stress (Miller, 2002; Hoshina et al., 2003; Sho et al., 2004), circumferential stress (Humphrey,
2002), influx of inflammatory cells (Choke et al., 2005; Pearce and Shively, 2006; Shimizu
et al., 2006; Middleton et al., 2007), and the formation of the intraluminal thrombus layer
(Vorp et al., 2001; Fontaine et al., 2002; Vorp and Vande Geest, 2005). In our following

101

Max. Diameter (mm)

20
19
18

T = 3 years
T = 6 years
T =10 years
T =15 years

17
16
15
14
1

1.5 2

2.5 3 3.5 4
Time (Years)

4.5 5

5.5

Figure 6.12: Expansion rates of the lesion in the ‘time dependent’ simulation case associated
with different time constants (T = 3, 6, 10, and 15 years).
study (Sheidaei et al., 2011) the current model is coupled with hemodynamic simulation in
an iterative manner and the effects of wall shear stress on the elastin degradation and, hence,
on the aneurysm growth are investigated.
In our final simulation, we assumed that further degradation of elastin is induced by the
stretch of elastin. Although it is not clear how much mechanical factors influence elastin
degradation during the aneurysmal growth, previous ex vivo studies show that arterial elastin
fails under uniaxial stretch (Gosline et al., 2002; Lillie and Gosline, 2007). Then, it might not
be unreasonable to hypothesize that the over-stretch of elastin accelerates elastin degradation. Previously, Wulandana and Robertson (2005) presented a model of a cerebral arterial
tissue that accounts for a stretch-dependent failure mechanism of elastin. Recently, Li and
Robertson (2009) have utilized a strain-induced damage model in a more structurally complicated model for balloon angioplasty.
We assumed that the amount of SM reduced proportionally with the initial elastin dam-

102

age in our model. The results do not indicate significant role of SM on AAA progression
possibly because we did not incorporate direct effects of SM loss on the extracellular matrix
turnover. However, there are increasing evidence on the critical role of vascular SM cells
in the aneurysm pathobiology through their activity and quantity (Curci, 2009). Vascular
SM cells are capable of producing high levels of matrix degrading enzymes in AAAs (Patel
et al., 1996; Crowther et al., 2000) as well as their inhibitors preventing the degeneration
(Allaire et al., 1998). Stretch/stress induced synthesis of extracellular matrix by SM cells
has also been documented (O’Callaghan and William, 2000; Sumpio et al., 1988). Apparently, the imbalance between proteolytic and synthetic activities of SM cells contributes to
the structural deterioration of arterial wall. In addition, advanced stages of AAA have been
associated with marked apoptosis of SM cells (L´pez-Candales et al., 1997; Zhang et al.,
o
2003; Thompson et al., 1997) which can exacerbate the weakening process. We need, however, more data to build a better model and predict the progressive weakening of the wall
structure and its failure.
The half-life of collagen in the arterial wall is reported to be 60-70 days in normal conditions (Humphrey, 2002) which can be reduced to 17 days in pathological conditions (Nissen
et al., 1978). Watton and Hill (2009) showed that reducing the half-life time nonlinearly
increased the expansion rate. Our preliminary results also showed that expansion rate is
inversely related to half-life time. We chose the half-life time to be about 35 days in our
simulations. Our results demonstrated that the kinetic parameter, Kg , has a key effect
on the expansion rate (Figure 6.9) consistent with other studies (Baek et al., 2005, 2006).
Nonetheless, the collagen half-life as well as the kinetic parameters were assumed to be fixed
during the evolution of aneurysm in our model, whereas in an actual AAA, these parameters may change during aneurysm growth due to pathological changes. Dynamic changes
of turnover parameters are multifaceted and complex and, hence, more studies are required
to quantify these pathological changes in aneurysm growth. Also, individuals may have
different degrees of mechano-sensitivity and stress-mediated turnover of collagen depending

103

of their physiological and pathological conditions. In fact, findings on the collagen content
variation during AAA growth are highly variable, e.g., decreased (Sumner et al., 1970), increased (Ghorpade and Baxter, 1996), and no change (Rizzo et al., 1989). The simulation
showed that Kg = 0.05 provides almost linear growth for all damage shapes (Figure 6.8)
and, especially, for Cases (1) and (3), growth rates were within the range observed for small
aneurysms (Nevitt et al., 1989; Baxter et al., 2008). With Kg ≤ 0.05 the lesion enlarged
continuously implying that the stress-mediated collagen turnover was not enough to return
the stress back to the homeostatic level although it reduced the local damage-induced stress
(Figure 6.9).
Results from the current model represent the early stage of aneurysms. As AAAs grow,
other factors such as intraluminal thrombus and perivascular boundary conditions should
also be taken into account for the model. Intraluminal thrombus not only has a direct
mechanical effect (Wang et al., 2002) but also changes the chemomechanical environment
and influences the strength of the lesion and its progression (Vorp et al., 2001; Taylor and
Humphrey, 2009).
Remodeling of the constituents and the mechanism involved in their deposition can have
a great impact on the mechanical properties of the evolved tissue. Computational models
of collagen remodeling assumed that either the principal strains (Driessen et al., 2004, 2008;
Boerboom et al., 2003) or stresses (Hariton et al., 2007b; Baek et al., 2006) govern the
orientation of collagen fibers. With regard to AAA, marked increase of anisotropy was found
in the diseased tissue (Vande Geest et al., 2006) suggesting structural changes in aneurysm
formation. Watton et al. (2004) assumed fixed fiber directions for collagen fibers in their
AAA growth and remodeling simulations. In our model, collagen aligned towards existing
fiber directions. With newly deposited collagen during enlargement, anisotropy increased
in the circumferential direction, which is consistent with Vande Geest et al. (2006). In
addition, we chose the mechanical property of newly synthesized collagen fibers to be the
same (Gleason et al., 2004; Baek et al., 2006, 2007b; Humphrey, 1999). Elastin degradation

104

alone does not conduce to rupture in normal vessels (Dobrin et al., 1984). Impaired collagen
networking and its microstructural defects have been associated with advanced AAAs despite
an increase in the collagen concentration (Lindeman et al., 2010). Therefore, in addition to
alteration of fiber distribution, it might also be desirable to account for changes of collagen
type, cross-linking and fiber thickness in the model (Driessen et al., 2008). A future extension
of the current model can be towards accounting for alteration of collagen structure during
the AAA progression and different mechanical properties of the newly synthesized collagen
fibers.
Prescribing the in vivo material and geometric parameters that satisfy the balance of
momentum and homeostatic conditions is a major challenge for patient-specific G&R simulations. In the current work we approximated the thickness by the Laplace relation, σh = P r/h
where r is the first principal radius of curvature. For a more complex geometry, however,
there is a need for developing numerical techniques to identify the spatial distribution of
parameters and integrating them with the developed vascular adaptation model.
Previous computational studies of vascular adaptation have considered mostly axisymmetric and simple geometries, while clinical application of the G&R simulation demands
patient-specific anatomical information. Watton et al. (2004) simulated an asymmetric
aneurysm by assuming an axisymmetric degradation and considering an effective pressure
(instead of constant pressure) which resembles the contact with spine. Although varying the
effective pressure resulted in an asymmetric aneurysm growth, their initial configuration was
still axisymmetric. Recently, Kuhl et al. (2007) simulated stress-induced growth of a medical
image-based model of human aorta based on the concept of incompatible growth combined
with open system thermodynamics. They simulated the effect of balloon or stent exposure
in atherosclerotic patients by a high internal pressure locally applied to the wall. Their
work, however, was based on a single constituent approach and did not consider growth and
remodeling of multiple constituents. On the other hand, many studies have used patientspecific geometries for estimating accurate stress levels and the rupture risk of aneurysms

105

without a G&R mechanism (Raghavan and Vorp, 2000; Fillinger et al., 2002, 2003; Speelman et al., 2007). However, Vorp and Vande Geest (2005) suggested that “despite recent
reports, it should be noted that evaluation of rupture potential based on only one of these
parameters –stress or strength – is not sufficient because a region of the AAA wall that
is under elevated wall stress may also have a higher wall strength”. We submit that, by
coupling a G&R model with a patient-specific geometric model, a computational simulation
can provide information about both stress and structural strength during aneurysm growth
and, hence, better prediction of the rupture can be acquired. Therefore, the computational
model presented here will provide a useful foundation towards a patient-specific modeling of
AAAs.

106

Chapter 7
MEDICAL IMAGE-BASED SIMULATION OF ABDOMINAL
AORTIC ANEURYSM GROWTH

7.1

A framework for modeling of AAA growth based
on medical images

Here we summarize our efforts in previous chapters and present a framework for simulation of
AAA enlargement with more simulation examples. Figure 7.1 illustrates a schematic diagram
of the steps taken. First, a 3-D computational geometry of a healthy aorta is constructed by
segmentation of medical images of a healthy aorta. For the aneurysm simulation, we use a
central region of the geometric model. Considering the vessel wall as a thin membrane, it is
convenient to parameterize the 3-D wall surface geometry with only two spatial parameters.
That is, any point on the wall surface can be characterized by its longitudinal position (l)
along the centerline of the artery as well as its azimuthal orientation (θ) with respect to a
reference direction (See the 2-D plot in Fig. 7.1 and chapter 4 for technical details).
Here, we employ the proposed inverse optimization (chapter 4) to estimate spatial distribution of geometric properties, e.g. wall thickness and anisotropy, such that a target
mechanical homeostasis is satisfied and the in-vivo geometry is maintained under in-vivo

107

Segmentation of
Medical Images & 3-D
Model Reconstruction

Computational Mesh
Generation

2-D Parameterization
of a 3-D Surface

Longitud. param (m).

0.07
0.06
0.05
0.04
0.03
0.02
0.01
Azimuthal Param. (deg.)
0
−200 −100
0
100 200
Identification of
in vivo Parameters

Anisotropy Distribution

Wall Thickness
Distribution

Induction of a Pathologic
Degeneration to Initiate G&R

Altered Mass
Distribution

Wall Stress
Distribution

G&R Simulations &
Stress Analysis

Figure 7.1: Schematic diagram of the steps taken for simulating AAAs using a medical
image-based geometry and material parameters of a healthy aorta as the initial state.

108

condition. The general form of the optimization objective function is given as
∫
J

=
Ω

||u(h, αk ) − uimage ||2 dA + ξ

∫ ∑( k
σ (h, αk )
Ω k

σh

)2
−1

dA

(7.1)

where u is the FE solution for displacement field and uimage is the displacement field from
the biomedical image. That is uimage = 0 when the in vivo configuration obtained from the
medical image is assumed to be the reference configuration. σ k is a scalar measure of stress
of the constituent k obtained from the FE analysis and σh is its homeostatic value. ξ adjusts
the minimization weight on each additive term in the objective function. (h, αk ) are the
unknown wall thickness and the orientation of the major diagonal fiber family controlling
the tissue anisotropy, approximated by

h(l, θ) =

I
∑

h
{βj ϕj (l, θ)}

(7.2)

j=1

αk (l, θ) =

I
∑

k
{βj ϕj (l, θ)},

(7.3)

j=1
h k
where (βj , βj ) are unknown variables associated with their base functions, ϕj (l, θ), which are

defined on the 2-D computational domain. A direct search method (Nelder-Mead Simplex)
is used to estimate spatial distributions of thickness and collagen fiber orientation that
minimize the objective function (See chapter 4 for details).
Elastin has been widely observed to be deficient in ruptured and unruptured AAAs.
Even though the factors that trigger AAA G&R have not been completely identified yet,
elastin degeneration is considered the pathologic factor contributing to the growth of AAAs
in this study. Then, as in chapter 6, G&R simulations are triggered by inducing damage/degradation to elastin with different spatial and temporal distributions as

R(X, t) = R(l, θ, t) = Rs (l, θ)Rt (t),

109

(7.4)

where the damage ratio, R ∈ [0, 1], is the ratio of the degenerated elastin mass to the initial
mass. Note that the spatial distribution, Rs , is defined on the mapped domain of (l, θ).

7.2

Simulation cases

In the first set of simulations, AAAs are initiated by inducing instantaneous damage to
elastin with different spatial shapes considering homogeneous distribution of kinetic parameters (Cases 1-6). Case 7 considers an instantaneous damage but with either homogeneous
i
or inhomogeneous distribution of Kg (= Kg ). In Case 8, after a portion of elastin is instan-

taneously removed, further time dependent degradation (Rt ) occurs while Kg is assumed to
depend on the local damage at the given time, i.e. R(l, θ, t) (See Fig. 7.2).

0
1

Damage Ratio, R(l, θ, t)
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0.05
0.045

0.9

0.04
Rt (t)

0.8
Kg

0.035
0.7
0.03
0.6

0.5
0

0.025

500

1000

1500
2000
Time (days)

2500

0.02
3000

Figure 7.2: Variation of damage ratio with respect to time (Rt vs. time) and variation
of the sensitivity parameter with respect to the local damage (Kg vs. R(l, θ, t)). Solid
line represents the time dependent damage ratio and dashed line represents the damage
dependent variation of Kg .

110

240
220
200
180
160
140
120
100
80

Case 3

240
220
200
180
160
140
120
100
80

0.26
0.22
0.18
0.14
0.1
0.06
0.02

7.3

1700 Days

Case 2
0.26
0.22
0.18
0.14
0.1
0.06
0.02

2700 Days

240
220
200
180
160
140
120
100

0.26
0.22
0.18
0.14
0.1
0.06
0.02

2700 Days

Case 1

Results

Figure 7.3 shows a set of distinct shapes of spatial damage (distributions of elastin areal mass
density) and the resulting AAAs with the corresponding von Mises stress distributions.
The time histories of the variations of the peak von Mises stress during enlargement of
5 aneurysms are plotted in Fig. 7.4a. Evidently, inducing damage is followed by a sharp
increase in the peak stress (on regions that match the degenerated regions) which, after
a time delay, begins to gradually diminish via stress-mediated collagen production. The

111

Case 5

240
220
200
180
160
140
120
100
80

0.26
0.22
0.18
0.14
0.1
0.06
0.02

Case 6

240
220
200
180
160
140
120
100
80

0.26
0.22
0.18
0.14
0.1
0.06
0.02

2700 Days
2700 Days

240
220
200
180
160
140
120
100
80

2700 Days

Case 4
0.26
0.22
0.18
0.14
0.1
0.06
0.02

Figure 7.3: left: Areal mass distributions (kg/m2 ) for elastin after introducing damage. right:
the resultant distributions of von Mises stress (kP a) due to G&R. The arrows indicate the
regions of elevated stress during AAA expansion.

112

decline of peak stress is initially sharper but becomes slower during the AAA expansion.
In some cases the stress slightly increases again after some time (e.g. around 4 years) on
the regions where the degradation is not necessarily maximal (See regions specified by black
arrows in Fig. 7.3). The maximum diameter of aneurysms in all cases changes linearly over
time but at different rates depending on the damage shape (Fig. 7.4b).
Figures 7.5 and 7.6 compare spatial correlations between von Mises stress and both the
rate of mass production and the rate of nodal displacement for Cases 1 and 5. The maximum
rates of mass production and displacement gradually increase with aneurysm expansion while
the peak stress reduces. In general, the rate of mass production increases with the stress,
but regions of high stress correspond to a wider range of mass production and displacement
rates especially in advanced stages. Compared to the rate of mass production, the rate of
displacement (i.e., radial expansion) has a weaker correlation with wall stress.
In a different simulation case (Case 7), an aneurysm is generated with a homogeneous
(Fig. 7.7a) as well as inhomogeneous (Fig. 7.7b) distribution of the kinetic/sensitivity
parameter. Assuming lower value for the kinetic parameter (i.e. Kg = 0.02) on the region
where aneurysm sac forms leads to lower rates of the local mass production and further
enlargement. Insufficient local mass production leads to an increasing trend in the local wall
stress during enlargement (Fig. 7.7b).
Figure 7.8 shows the distributions of elastin and collagen areal mass density and the von
Mises stress for Case 8 where the elastin damage is time-dependent and Kg is assumed to be
damage-dependent (relations shown in Fig 7.2). Increasing degradation along with gradually
decreasing Kg accelerate the AAA enlargement and heightens wall stress (up to 450 kPa).

7.4

Discussion

Consistent with previous studies (Watton and Hill, 2009), collagen turnover showed to compensate for the loss of elastin and tended to normalize wall stress via the stress-mediated

113

Peak Wall Stress (kPa)

Max. Diameter (mm)

450

Case 1
Case 3
Case 4
Case 5
Case 6

400
350
300
250
200
150
100
0

34
32
30
28
26
24
22
20
18
16
14
0

1

2

3 4 5 6
Time (years)

(a)
7 8

Case 1: y=2.5x + 14
Case 3: y=1.6x + 15
Case 4: y=1.3x + 15
Case 5: y=2.3x + 15
Case 6: y=2x + 14

(b)
1

2

3 4 5 6
Time (years)

7

8

Figure 7.4: Variations of peak von Mises stress (a) and maximum diameter (b) with G&R
time for different cases. Equations for best fit lines are shown (lines are not included in the
graph).

114

Rate of mass prod. (kg/m2 .day)

-3
Displacement rate (10 mm/day)

0.07
0.06
0.05
0.04
0.03

2
kg/m

mm

3
2.6
2.2
1.8
1.4
1
0.6
0.2

1.2
1.1
1
0.9
0.8
0.7
0.6
0.5

Wall Thickness
0.02 @ 2700 Days

Collagen

0.01
0

50

(a)
100 150 200 250 300 350
von Mises stress (kPa)

10
8

300 days
1250 days
1700 days
2700 days

6
4
2
0
−2
50

100

150
200
250
von Mises stress (kPa)

(b)
300

Figure 7.5: Spatial relations between von Mises stress and both the rate of collagen mass production as well as the rate of wall displacement for Cases 1 at different stages of growth. Wall
thickness distribution and collagen concentration are illustrated for each case corresponding
to 2700 days of G&R.

115

Rate of mass prod. (kg/m 2.day)

0.06

Displacement rate (10-3 mm/day)

10

0.05
0.04
0.03
0.02
0.01
0

8
6

50

kg/m2
3
2.6
2.2
1.8
1.4
1
0.6
0.2
Collagen

mm
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5
Wall thickness
@2700 Days

(c)
100 150 200 250 300 350
von Mises stress (kPa)

300 days
1250 days
1700 days
2700 days

4
2
0
−2
(d)
−4
80 100 120 140 160 180 200 220 240 260
von Mises stress (kPa)

Figure 7.6: Spatial relations between von Mises stress and both the rate of collagen mass production as well as the rate of wall displacement for Cases 5 at different stages of growth. Wall
thickness distribution and collagen concentration are illustrated for each case corresponding
to 2700 days of G&R.

116

Rate of mass prod.(kg/m 2.day)
0.
0.
0. 0.0
0.
01
00
02
25
01
5
5

2
kg/m

Rate of mass prod. (kg/m2.day)
0.
0
0
0
00 0.0 .01 0.0 .02 0.0 .03
2
3
5
5
5
5
1

Collagen
2
kg/m
1.15
0.95
0.75
0.55
0.35
0.15

2
kg/m

0.26
0.22
0.18
0.14
0.1
0.06
0.02

1.15
0.95
0.75
0.55
0.35
1
0. 5

Elastin

Collagen
Case 7@1700 Days

0
100

150

(a)
200 250 300 350 400
von Mises stress (kPa)

10 day
295 day
700 day
1245 day
1700 day

@1700 Days

(b)

0
100 150 200 250 300 350 400 450
von Mises stress (kPa)

Figure 7.7: Spatial relation between the rate of collagen mass production and von Mises
stress for Case 7 using a homogeneous distribution of kinetic parameters (a; Kg = 0.05) as
well as spatially variable kinetic parameter (b; Kg = 0.02 at the AAA sac and Kg = 0.05
otherwise) at different times.

117

(c)

Peak Wall Stress (kPa)

2700 Days

kPa
440
400
360
320
280
240
200
160
120
80

kg/m2
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2

(b)

32
400 (d)

28

300

24

200

20

16
100
0 1 2 3 4 5 6 7 8
Time (Years)

Max. Diameter (mm)

kg /m2
0.26
0.22
0.18
0.14
0.1
0.06
0.02
(a)

Figure 7.8: Elastin (a) and collagen mass (b) areal density along with von Mises stress
(c) distributions for Case 8 after 7 years (2700 days) of G&R. Time history of maximum
diameter and peak von Mises stress is also shown (d).

118

kinetics of G&R process. Nevertheless, assuming Kg = 0.05, the peak wall stress never returned to the homeostatic value leading to continuous enlargement with a linear trend which
is typical of small AAAs (Fig. 7.4). The high stress regions in Cases 1 and 5 (Fig. 7.3) were
found to be consistent with other FE studies (Vorp et al., 1998; Scotti et al., 2005; Doyle
et al., 2010; Hua and Mower, 2001; Giannoglou et al., 2006), suggesting regions of inflections
as the high stress regions. We noticed a strong spatial correlation between the rate of collagen production and wall stress in earlier stage of growth. Later, even though the peak wall
stress was reduced, collagen was produced at an increasing rate leading to high rates of wall
displacement (Figs. 7.5 and 7.6). This seems to be mainly because the rate of mass production was assumed to be linearly related to the number of cells (depending on the collagen
mass) at a given time. Although there was a strong spatial correlation between damage
and the rate of mass production at all times (not shown), a weaker correlation was found
between the local damage and wall displacement in advanced AAAs (Fig. 7.9). In Case 1,
wall displacement was more global and not restricted to the damaged regions while in Case
3 the damaged regions were subject to high wall displacement (i.e., large radial expansion)
(Fig. 7.9). Interestingly, in Case 1, although the damage was uniformly distributed around
the circumference (Fig. 7.3), points on the convex side were subject to more displacement
than points on the concave side (See Figs. 7.3 and 7.9).
The time-dependent case was inspired by a possibility of increasing elastin degradation
due to exacerbated hemodynamic conditions in AAAs as the geometry evolves (Mohan et al.,
1999; Walpola et al., 1993). Sheidaei et al. (2011) demonstrated that low wall shear stress
may induce further degradation during aneurysm progression and accelerate the lesion expansion. In addition, the damage-induced weakening of the G&R sensitivity was introduced
because some studies have suggested an increase in apoptosis (Lopez-Candales et al., 1997;
Thompson et al., 1997) or phenotype modulation (Ailawadi et al., 2009) of SM cells during the AAA enlargement. The imbalance between extracellular matrix synthesis (Sumpio
et al., 1988) and degradation (Crowther et al., 2000), mediated by altered SM cells can be

119

Nodal Displacement (mm)

14

Case 1

12
10
8
6
4
2
0
0

0.2

0.4
0.6
0.8
Damage Ratio, R(l, θ)

1

Nodal Displacement (mm)

14
12
10

Pts. on convex side
Pts. on concave side
Case 3
300 Days
1250 Days
2700 Days

8
6
4
2
0

0

0.2

0.4
0.6
Damage Ratio R(l,θ)

0.8

1

Figure 7.9: Spatial relation between the nodal displacement and the damage ratio for Cases
1 and 3 at different times. A number of points on the convex and concave regions (with
reference to the healthy geometry) are distinguished with different markers at 2700 days.

120

the reason for progressive weakening and rupture of the wall. Recently, Maier et al. (2011)
have quantitatively found a spatial correlation between wall stress and the level of in vivo
tissue reactions which can indicate either synthetic or proteolytic activities in the arterial
wall tissue.
In this study, we were able to capture relatively complex 3-D shapes of AAAs by prescribing distinct shapes for spatial damages. Even with a fixed spatial distribution of elastin
damage, we noticed that regions of elevated wall stress can change in the course of expansion
(Fig. 7.3). Interestingly, during the evolution of some AAA shapes, the peak wall stress locally increased again (Fig. 7.4). In Cases 1 and 5, the lesions were larger in diameter while in
Cases 3 and 6 the lesions developed higher wall stress after 2700 days (7 years), implying the
influence of aneurysm shape on the peak wall stress, beyond what the maximum diameter
criterion may solely suggest.
Together, these results emphasize the impact of geometrical complexities on complex
stress distributions during AAA expansion as suggested by Sacks et al. (1999). Other studies have sought correlations between wall stress and quantified AAA shapes, asymmetry,
tortuosity, as alternative indicators of rupture (Shum et al., 2011; Vorp et al., 1998; Georgakarakos et al., 2010; Doyle et al., 2009; Giannoglou et al., 2006). Moreover, we found that
the regions of low wall thickness initially coincide with the damaged regions but they too
change as the lesions evolve. In some cases they may correspond to high stress regions (e.g.
Cases 1 and 5; Figs. 7.3 and 7.5 and 7.6) while in other cases there may be no correspondence
(other cases not shown). Coincidence of the low wall thickness and elevated wall stress could
be indicative of an impending rupture in advanced stages. As an advantage, current study is
capable of capturing morphological complexities during enlargement while at the same time
the evolution of wall properties during G&R is being reflected.
Aneurysms are often accompanied by bending of the proximal neck area which influences
the hemodynamics and possibly intraluminal thrombus formation (Biasetti et al., 2010).
Figure 7.10 depicts the axial stretch distributions for different simulation cases. Aneurysms

121

1700 Days

1700 Days
2700 Days

Case 4
1.45
1.35
1.25
1.15
1.05
0.95
0.85
0.75

Case 2
1.45
1.35
1.25
1.15
1.05
0.95
0.85
0.75
Case 7
1.45
1.35
1.25
1.15
1.05
0.95
0.85
0.75

1700 Days

Case 1
1.45
1.35
1.25
1.15
1.05
0.95
0.85
0.75

Figure 7.10: Distributions of axial stretch for different aneurysms.
enlarge in both circumferential and axial directions and because of the fixed boundary condition, axial stretch on the neck region reduces to values less than 1. Even if not the main
reason, this can be used as an indicator of where bending may occur in later stages (Jackson
et al., 2005). It is noteworthy that we have modeled the vessel wall as a membrane while
models that account for bending of the wall are necessary to study buckling of the vessel
walls.

122

Chapter 8
CONCLUSION

From the practical point of view, we are ultimately interested to model AAA expansion from
an existing AAA, which is an intermediate stage of the disease. While the current model uses
the assumption of homeostatic condition for the healthy aorta to prescribe initial conditions
(e.g. thickness, fiber alignment, kinetic parameters) of the simulation, less information is
available for the spatial and temporal distribution of the parameters for a diseased artery.
Despite this challenge, we believe that the current model can be used for estimating the
material and kinetic parameters of an existing AAA and, hence, predicting patient-specific
AAA expansion if used in a broader framework. For example, using the same model but in
an inverse scheme, it seems feasible to regenerate a patient-specific AAA shape by estimating
spatio-temporal distributions of elastin degradation as well as the G&R sensitivity/kinetic
parameter. Cyclic strain quantified from medical images can also be incorporated into the
inverse scheme such that not only the AAA shape but also the local strain is predicted
more accurately. That requires further development of the model regarding the boundary
condition, e.g. the spine support, which as a growth barrier significantly changes the resulting
AAA shape and stress/strain distributions in the G&R model.
Lack of verification for the estimated distributions of the wall thickness and fiber orientation in healthy condition (prior to G&R) and during the AAA development is a limitation of
the current study. Recent advances in medical imaging and segmentation techniques make

123

Fluid-Solid Growth
(FSG) Simulation
Hemodynamic
G&R
Simulation
Simulation

Identification of
in vivo conditions
In vivo
Material
Properties

Boundary
Conditions

Initialization

ILT Model

(Pathogenesis)
Perivascular
Tissue Model

Vessel
Wall
Mixture
Model

FSI

or
Rigid Wall
Assumption

Figure 8.1: A comprehensive simulation framework aimed at patient-specific modeling of
AAA G&R.
it possible to quantify in-vivo wall thickness (Shum et al., 2010; Martufi et al., 2009) and
verify the patient-specific simulations of AAAs. Another limitation of the current model
is that the role of intraluminal thrombus layer (ILT) is overlooked. Thrombus has been
found to either protect the wall by reducing the wall stress (Inzoli et al., 1993; Thubrikar
et al., 2003; Wang et al., 2002) or compromise its integrity via biochemical environment
(Vorp et al., 2001; Fontaine et al., 2002). Figure 8.1 shows all the required elements of a
comprehensive patient-specific model of AAA G&R. This study basically covers the topics
in black and red boxes, but the focus has been mainly on identification of in-vivo properties.
The model of ILT and perivascular tissue are the steps needed to enhance the G&R simulation which, on the other hand, has a strong coupling with hemodynamic simulation (blue
or green boxes). Hemodynamic simulations in evolving aneurysms have been performed in
collaborative studies either with rigid wall assumption (Sheidaei et al., 2011) or as fluid-solid
interaction framework (Figueroa et al., 2009).
In conclusion the current study is an important step toward patient-specific modeling of
AAAs. Apparently accurate estimation of rupture potential requires further developments
of the model as well as more detailed information about pathophysiology and biochemical
reaction involved.

124

Appendices

125

Appendix A
Equilibrium equations considered for inflation-extension test of
tubular vessels

We use a cylindrical coordinate system to model axial extension and inflation of a thin
artery. The positions of a particle can be expressed by cylindrical coordinates {R, Θ, Z}
and {r, θ, z} in the reference (unloaded) and current (loaded) configuration, respectively
(Humphrey, 2002). The motion due to the inflation and axial stretch is assumed as

r = r(R),

θ = Θ,

z = λz Z.

(A.1)

The deformation gradient for the large deformation is




∂r
 ∂R


F= 0


0

0

0 

r
0 

R

0 λz

(A.2)

We assume that the arterial wall material is incompressible, so that
∂r
R
=
.
∂R
rλz

126

(A.3)

Using the incompressibility, the inner diameter can be obtain as
√
ri =

2
ro −

2
2
(Ro − Ri )
λz

(A.4)

where Ri , Ro and ro are the inner and outer diameter in the reference configuration and
outer diameter in the current configuration. For axisymmetric condition, the equations of
motion in the absence of body force and the force due to acceleration are reduced to
1 ∂(rTrr ) ∂Tzr Tθθ
+
−
= 0
r ∂r
∂z
r
∂Tzθ 1 ∂(rTrθ ) Tθr
+
+
= 0
∂z
r ∂r
r
∂Tzz 1 ∂(rTrz )
+
= 0.
∂z
r ∂r

(A.5)
(A.6)
(A.7)

Furthermore, Trθ = Trz = 0 for the thin tube, and Tzθ = 0 for the absence of torsion. Then,
the remaining equation is
∂Trr Trr − Tθθ
+
= 0
∂r
r
∂Tzz
= 0,
∂z

(A.8)
(A.9)

where assuming incompressible material

ˆ
Trr = −p + Trr

(A.10)

ˆ
Tθθ = −p + Tθθ

(A.11)

ˆ
Tzz = −p + Tzz .

(A.12)

127

The internal pressure Pi is given for the inflated tube, and Trr (r = ri ) = −Pi and Trr (r =
ro ) = 0 where ri and ro are in and out radius for the tube. The integration of (A.8) gives
∫ ro ˆ
ˆ
Trr − Tθθ
Trr (r) =
dξ.
ξ
r

(A.13)

Internal pressure Pi and the axial force Fz are given by integral forms as
∫ ro ˆ
ˆ
Trr − Tθθ
Pi = −
dr
r
ri
∫ ro
2
2
ˆ
(−p + Tzz )(2πr)dr − π(ri − rc )Pi .
Fz =

(A.14)
(A.15)

ri

where rc is inner diameter of cannula. The Lagrange multiplier p can be solved using (A.10)
and (A.13), then (A.15) becomes
∫ ro ( ∫ ro ˆ
)
ˆ
Trr − Tθθ
ˆrr + Tzz rdr − π(r2 − r2 )Pi .
ˆ
Fz = 2π
dξ − T
c
i
ξ
ri
r
Since,

(A.16)

∫ ro ∫ ro ˆ
∫ ro
2
ˆ
ri Pi
Trr − Tθθ
r ˆ
ˆ
r
dξdr =
+
(Trr − Tθθ )dr
ξ
2
ri
r
ri 2

(A.17)

The equation (A.16) results in
∫ ro (
Fz = 2π
ri

ˆ
ˆ
T + Trr )
2
ˆ
Tzz − θθ
rdr + πrc Pi .
2

(A.18)

∫
For a thin membrane, the integration rro f (r)dr can be approximated by f (rm )h where
i
rm = (ri + ro )/2 and h = ro − ri . Hence
ˆ
ˆ
h(Tθθ − Trr )
rm
2
2
ˆ
ˆ
Fz = 2πrm h(Tzz − Trr ) − π(rm − rc )Pi .
Pi =

128

(A.19)
(A.20)

2
ˆ
ˆ
Considering Tθθ = Trr + πrm Pi we also have

2
ˆ
ˆ
ˆ
Fz = πrm h(2Tzz − Trr − Tθθ ) + πrc Pi .

(A.21)

The inner radius (A.4) can be rewritten as
√
ri =

2
ro −

where Rm = (Ro + Ri )/2 and H = Ro − Ri .

129

2Rm H
λz

(A.22)

Appendix B
Nelder-Mead Simplex Algorithm

B.1

Introduction

Since its publication in 1965, the Nelder-Mead simplex algorithm (Nelder and Mead, 1965)
has been widely used for nonlinear unconstrained optimization. The Nelder-Mead, as direct search method, attempts to minimize a scalar-valued nonlinear function of n variables
(f (x); x ∈ Rn ) using only function calculation, without any derivative information. Each
iteration of the algorithm starts with a simplex, specified by its n + 1 vertices and the associated function values. Four scalar parameters must be specified for the algorithm: coefficients
of reflection (ρ), expansion (µ), contraction (γ), and shrinkage (σ). The universal choices
for those parameters are (Lagarias et al., 1998):

ρ = 1, µ = 2, γ = 0.5 and σ = 0.5

B.2

(B.1)

Algorithm statement

At the beginning of the k th a simplex is given consisting of n + 1 vertices (points in Rn ).
Then the following steps are taken (Lagarias et al., 1998):

130

1. ORDER. Order the n + 1 vertices, xk , ..., xk , such that
1
n+1
k
k
k
f1 ≤ f2 ≤ ... ≤ fn+1

where

fik = f (xk ).
i

(B.2)

ˆ
xk = Σn xk /n.
i=1 i

(B.3)

2. REFLECT. Compute the reflection point as
ˆ
xk = xk + ρ(ˆ k − xk )
x
r
n+1

where

k
k
k
k
Evaluate fr = f (xk ). If f1 ≤ fr ≤ fn , accept the reflected point and terminate iteration.
r
k
k
3. EXPAND. If fr < f1 , calculate the expansion point as

ˆ
ˆ
xk = xk + µ(xk − xk ).
e
r

(B.4)

k
k
k
Evaluate fe = f (xk ). If fe < fr , accept xk and terminate the iteration. Otherwise, accept
e
e

xk and terminate the iteration.
r
k
k
4. CONTRACT. If fr ≥ fn , perform a contraction as
k
k
k
a. OUTSIDE. If fn ≤ fr ≤ fn+1 , perform outside contraction as

ˆ
ˆ
xk = xk + γ(xk − xk ).
c
r

(B.5)

k
k
k
Evaluate fc = f (xk ). If fc ≤ fr , accept xk and terminate the iteration. Otherwise, go to
c
c

step 5.
k
k
b. INSIDE. If fr ≥ fn+1 , perform inside contraction as

ˆ
xk = xk − γ(ˆ k − xk ).
x
cc
n+1

(B.6)

k
k
k
Evaluate fcc = f (xk ). If fcc ≤ fn+1 , accept xcc and terminate the iteration. Otherwise, go
cc

to step 5.
k
3. SHRINK. Evaluate f at n points, vi = xk + σ(xk − xk ), i = 2, ..., n + 1. The unordered
1
1
i

131

k
k
vertices for the simplex at the next iteration consist of xk , v2 , ..., vn+1 .
1

B.3

Stopping criteria

A stopping criterion is chosen based on both the relative size of the simplex and function
values at vertices of the simplex as (Torczon, 1989):
1
max ||xk − xk || < δ
1
∆ 2≤j≤n+1 j

(B.7)

f (xk ) − f (xk ) < ϵ ,
n+1
1

(B.8)

where xk is the j th vertex of the simplex and a vector comprised of all optimization variables
j
th
at k th iteration. xk and xk
1
n+1 are the “best” and “worst” vertices of the simplex at k

iteration and ∆ = max(1, ||xk ||).
1

132

BIBLIOGRAPHY

133

Bibliography
Ailawadi, G., Moehle, C. W., Pei, H., Walton, S. P., Yang, Z., Kron, I. L., Lau, C. L., Owens,
G. K., 2009. Smooth muscle phenotypic modulation is an early event in aortic aneurysms.
J Thorac Cardiovasc Surg 138, 1392–1399.
Akaike, H., 1974. New look at statistical-model identification. IEEE Transactions on Automatic Control 19, 716–723.
Akaike, H., 1981. Likelihood of a model and information criteria. Journal of Econometrics
16, 3–14.
Allaire, E., Forough, R., Clowes, M., 1998. Local overexpression of timp-1 prevents aortic
aneurysm degeneration and rupture in a rat model. J Clin Invest 102 (7), 1413–1420.
Anderson, D. R., Burnham, K. P., White, G. C., 1994. AIC model selection in overdispersed
capture-recapture data. Ecology 75, 1780–1793.
Baek, S., Gleason, R. L., Rajagopal, K. R., Humphrey, J. D., 2007a. Theory of small on
large: Possible utility for computations of fluid-solid interactions in arteries. Computer
Methods in Applied Mechanics and Engineering 196, 3070–3078.
Baek, S., Rajagopal, K. R., Humphrey, J. D., 2005. Competition between radial expansion and thickening in the enlargement of an intracranial saccular aneurysm. Journal of
Elasticity 80, 13–31.
Baek, S., Rajagopal, K. R., Humphrey, J. D., 2006. A theoretical model of enlarging intracranial fusiform aneurysms. ASME Journal of Biomechanical Engineering 128, 142–149.
Baek, S., Valent´ A., Humphrey, J. D., 2007b. Biochemomechanics of cerebral vasospasm
ın,
and its resolution: II. constitutive relations and model simulations. Annals of Biomedical
Engineering 35, 1498–1509.
Baxter, B. T., Terrin, M. C., Dalman, R. L., 2008. Medical management of small abdominal
aortic aneurysms. Circulation 117 (14), 1883–1889.
Beck, J. V., Arnold, K. J., 1977. Parameter Estimation in Engineering and Science. Wiley,
New York.
Biasetti, J., Gasser, T., Auer, M., Hedin, U., Labruto, F., 2010. Hemodynamics of the normal
aorta compared to fusiform and sacuular abdominal aoritic aneurysms with emphasis on a
potential thrombus formation mechanism. Annals of Biomedical Engineering 38, 380–390.
134

Boerboom, R. A., Driessen, N. J. B., Bouten, C. V. C., J. M. Huyghe, F. P. T. B., 2003.
Finite element model of mechanically induced collagen fiber synthesis and degradation in
the aortic valve. Ann Biomed Eng 31 (9), 1040–1053.
Boutouyrie, P., Bussy, C., Lacolley, P., Girerd, X., Laloux, B., Laurent, S., 1999. Association between local pulse pressure, mean blood pressure, and large-artery remodeling.
Circulation 100 (13), 1387–1393.
Burton, A. C., 1954. Relation of structure to function of the tissues of the wall of blood
vessels. Physiol Rev 34 (4), 619–642.
Choke, E., Cockerill, G., Wilson, W. R. W., Sayed, S., Dawson, J., Loftus, I., Thompson,
M. M., 2005. A review of biological factors implicated in abdominal aortic aneurysm
rupture. Eur J Vasc Endovasc Surg 30 (3), 227–244.
Chuong, C. J., Fung, Y. C., 1984. Compressibility and constitutive relation of arterial wall
in radial compression experiments. J Biomech 17, 35–40.
Crowther, M., Goodall, S., Jones, J. L., Bell, P. R., Thompson, M. M., 2000. Increased matrix
metalloproteinase 2 expression in vascular smooth muscle cells cultured from abdominal
aortic aneurysms. J Vasc Surg 32, 575–583.
Cummins, P. M., Sweeney, N. V., Killeen, M. T., Birney, Y. A., Redmond, E. M., Cahill,
P. A., 2007. Cyclic strain-mediated matrix metalloproteinase regulation within the vascular
endothelium: a force to be reckoned with. Am J Physiol Heart Circ Physiol 292 (1), H28–
H42.
Curci, J. A., 2009. Digging in the ”soil” of the aorta to understand the growth of abdominal
aortic aneurysms. Vascular 17, S21–S29.
Demiray, H., 1972. A note on the elasticity of soft biological tissues. J Biomech 5, 309–311.
Dobrin, P. B., Baker, W. H., Gley, W. C., 1984. Elastolytic and collagenolytic studies of
arteries. Arch Surg 119 (4), 405–409.
Dorfmann, A., Wilson, C., Edgar, E. S., Peatte, R. A., 2010. Evaluating patient-specific
abdominal aortic aneurysm wall stress based on flow-induced loading. Biomech Model
Mechanobiol 9 (2), 127–139.
Doyle, B. J., Callanan, A., Burke, P. E., Grace, P. A., Walsh, M. T., Vorp, D. A., McGloughlin, T. M., 2009. Vessel asymmetry as an additional diagnostic tool in the assessment of
abdominal aortic aneurysms. J Vasc Surg 49, 443–454.
Doyle, B. J., Cloonan, A. J., Walsh, M. T., Vorp, D. A., McGloughlin, T. M., 2010. Identification of rupture locations in patient-specific abdominal aortic aneurysms using experimental
and computational techniques. J Biomech 43, 1408–1416.
Draney, M., Arko, F., Alley, M., R., H., Pelc, N., Zarins, C., Taylor, C., 2004. Quantification
of vessel wall motion and cyclic strain 337 using cine phase contrast mri: in vivo validation
in the porcine aorta. Magn Reson Med 52, 286–295.
135

Driessen, N., Wilson, W., Bouten, C., Baaijens, F., 2004. A computational model for collagen
fibre remodelling in the arterial wall. Journal of Theoretical Biology 226, 53–64.
Driessen, N. J. B., Cox, M. A. J., Bouten, C. V. C., Baaijens, F. P. T., 2008. Remodelling of the angular collagen fiber distribution in cardiovascular tissues. Biomech Model
Mechanobiol 7 (2), 93103.
Driss, A. G., Benessiano, J., Poitevin, P., Levy, B. I., Michael, J.-B., 1997. Arterial expansive
remodeling induced by high flow rates. Am J Physiol 272 (2), H851–H858.
Dye, W. W., Gleason, R. L., Wilson, E., Humphrey, J. D., 2007. Altered biomechanical
properties of carotid arteries in two mouse models of muscular dystrophy. Journal of
Applied Physiolology 103, 664–672.
Eberth, J. F., Gresham, V. C., Reddy, A. K., Popovic, N., Wilson, E., Humphrey, J. D.,
2009. Importance of pulsatility in hypertensive carotid artery growth and remodeling. J
Hypertens 27 (10), 2010–2021.
Emery, A. F., Nenarokomov, A. V., Fadale, T. D., 2000. Uncertainties in parameter estimation: the optimal experiment design. International Journal of Heat and Mass Transfer 43,
3331–3339.
Fadale, T. D., Nenarokomov, A. V., Emery, A. F., 1995. Uncertainties in parameter estimation:the inverse problem. International Journal of Heat and Mass Transfer 38, 511–518.
Figueroa, C. A., Baek, S., Taylor, C. A., Humphrey, J. D., 2009. A computational framework
for fluid-solid-growth modeling in cardiovascular simulations. Comput Methods Appl Mech
Eng 198, 3583–3602.
Fillinger, M. F., Marra, S. P., Raghavan, M. L., Kennedy, F. E., 2003. Prediction of rupture
risk in abdominal aortic aneurysm during observation: Wall stress versus diameter. J Vasc
Surg 37 (4), 724–732.
Fillinger, M. F., Raghavan, M. L., Marra, S. P., Cronenwell, J. L., Kennedy, F. E., 2002.
In vivo analysis of mechanical wall stress and abdominal aortic aneurysm rupture risk. J
Vasc Surg 36 (3), 589–597.
Flores, F. G., Onate, E., 2001. A basic thin shell triangle with only translational dofs for
large strain plasticity. Int J Numer Meth Engng 51, 57–83.
Fontaine, V., Jacob, M. P., Houard, X., Rossignol, P., Plissonnier, D., Angles-Cano, E.,
Michel, J. B., 2002. Involvement of the mural thrombus as a site of protease release and
activation in human aortic aneurysms. Am J Pathol 161 (5), 1701–1710.
Fung, Y. C., 1967. Elasticity of soft tissues in simple elongation. Am J Physiol 213, 1532–
1544.
Fung, Y. C., Fronek, K., Patitucci, P., 1979. Pseudoelasticity of arteries and the choice
of its mathematical expression. American Journal of Physiology - Heart and Circulatory
Physiology 237, H620–H631.
136

Gasser, T. C., Ogden, R. W., Holzapfel, G. A., 2006. Hyperelastic modelling of arterial
layers with distributed collagen fibre orientations. Journal of The Royal Society Interface
3, 15–35.
Gee, M. W., Forster, C., Wall, W. A., 2010. A computational strategy for prestressing
patient-specific biomechanical problems under finite deformation. Int J Numer Meth
Biomed Engng 26, 52–72.
Gee, M. W., Reeps, C., Eckstein, H. H., Wall, W. A., 2009. Prestressing in finite deformation
abdominal aortic aneurysm simulation. J Biomech 42, 1732–1739.
Georgakarakos, E., Ioannou, C. V., Kamarianakis, Y., Papaharilaou, Y., Kostas, T.,
Manousaki, E., Katsamouris, A. N., 2010. The role of geometric parameters in the prediction of abdominal aortic aneurysm wall stress. Eur J Vasc Endovasc Surg 39, 42–48.
Ghorpade, A., Baxter, B. T., 1996. Biochemistry and molecular regulation of matrix macromolecules in abdominal aortic aneurysms. Ann New York Acad Sci 800, 138–150.
Giannoglou, G., Giannakoulas, G., Soulis, J., Chatzizisis, Y., Perdikides, T., Melas, N.,
Parcharidis, G., Louridas, G., 2006. Predicting the risk of rupture of abdominal aortic
aneurysms by utilizing various geometrical parameters: Revisiting the diameter criterion.
Angiology 57, 487–494.
Glatting, G., Kletting, P., Reske, S. N., Hohl, K., Ring, C., 2007. Choosing the optimal fit
function: Comparison of the akaike information criterion and the F-test. Medical Physics
34, 4285–4292.
Gleason, R., Humphrey, J., 2004a. A mixture model of arterial growth and remodeling in
hypertension: altered muscle tone and tissue turnover. Journal of Vascular Research 41,
352–363.
Gleason, R., Taber, L., Humphrey, J., 2004. A 2-d model of flow-induced alterations in the
geometry, structure, and properties of carotid arteries. Journal of Biomechanical Engineering 126, 371–381.
Gleason, R. L., Humphrey, J. D., 2004b. A mixture model of arterial growth and remodeling
in hypertension: Altered muscle tone and tissue turnover. J Vasc Surg 41, 352–363.
Gosline, J., Lillie, M., Carrington, E., Guerette, P., Ortlepp, C., Savage, K., 2002. Elastic proteins: Biological roles and mechanical properties. Phil Trans R Soc London B
357 (1418), 121–132.
Gundiah, N., Ratcliffe, M. B., Pruitt, L. A., 2007. Determination of strain energy function for
arterial elastin: experiments using histology and mechanical tests. J Biomech 40, 586–594.
Guo, X. M., Kassab, G. S., 2004. Distribution of stress and strain along the porcine aorta
and coronary arterial tree. Am J Physiol Heart Circ Physiol 286, H2361–H2368.

137

Hansen, L., Wan, W., Gleason, R. L., 2009. Microstructurally motivated constitutive modeling of mouse arteries cultured under altered axial stretch. J Biomech Eng 131, 101015.
Hariton, I., deBotton, G., Gasser, T. C., Holzapfel, G. A., 2007a. Stress-driven collagen fiber
remodeling in arterial walls. Biomech Model Mechanobiol 6, 163–175.
Hariton, I., deBotton, G., Gasser, T. C., Holzapfel, G. A., 2007b. Stress-modulated collagen
fiber remodeling in a human carotid bifurcation. J Theor Biol 248 (3), 460–470.
He, C. M., Roach, M. R., 1994. The composition and mechanical properties of abdominal
aortic aneurysms. J Vasc Surg 20 (1), 6–13.
Hoi, Y., Meng, H., Woodward, S. H., Bendok, B. R., Hanel, R. A., Guterman, L. R.,
Hopkins, L. N., 2004. Effects of arterial geometry on aneurysm growth: threedimensional
computational fluid dynamics study. J Neurosurg 101 (4), 676–681.
Holzapfel, G. A., 2006. Determination of material models for arterial walls from uniaxial
extension tests and histological structure. Journal of Theoretical Biology 238, 290–302.
Holzapfel, G. A., Gasser, T. C., Ogden, R. W., 2000. A new constitutive framework for
arterial wall mechanics and a comparative study of material models. J Elasticity 61, 1–48.
Holzapfel, G. A., Gasser, T. C., Ogden, R. W., 2004. Comparison of a multi-layer structural
model for arterial walls with a fung-type model, and issues of material stability. ASME
Journal of Biomechanical Engineering 126, 264–275.
Holzapfel, G. A., Gasser, T. C., Stadler, M., 2002. A structural model for the viscoelastic
behavior of arterial walls: Continuum formulation and finite element analysis. European
Journal of Mechanics A-Solids 21, 441–463.
Holzapfel, G. A., Ogden, R. W., 2010. Constitutive modeling of arteries. Proc R Soc 466,
1551–1597.
Holzapfel, G. A., Sommer, G., Gasser, T. C., Regitnig, P., 2005. Determination of layerspecific mechanical properties of human coronary arteries with nonatherosclerotic intimal
thickening and related constitutive modeling. American Journal of Physiology-Heart and
Circulatory Physiology 289, 2048–2058.
Hoshina, K., Sho, E., Sho, M., Nakahashi, T. K., Dalman, R. L., 2003. Wall shear stress and
strain modulate experimental aneurysm cellularity. J Vasc Surg 37 (5), 1067–1074.
Hsu, H.-J., Lee, C.-F., Kaunas, R., 2009. A dynamic stochastic model of frequency-dependent
stress fiber alignment induced by cyclic stretch. PLoS ONE 4, e4853.
Hu, J., Bjorklund, A., Nyman, M., Gennser, G., 1998. Mechanical properties of large arteries
in mother and fetus during normal and diabetic pregnancy. Journal of Maternal-Fetal
Investigation 8, 185–193.
Hu, J.-J., Baek, S., Humphrey, J. D., 2007. Stress-strain behavior of the passive basilar
artery in normotension and hypertension. Journal of Biomechanics 40, 2559–2563.
138

Hua, J., Mower, W. R., 2001. Simple geometric characteristics fail to reliably predict abdominal aortic aneurysm wall stresses. J Vasc Surg 34, 308–315.
Huffel, S. V., Vandewalle, J., 1991. The Total Least Squares Problem: Computational Aspects and Analysis. Society for Industrial and Applied Mathematics, Philadelphia.
Humphrey, J., Rajagopal, K., 2002. A constrained mixture model for growth and remodeling
of soft tissues. Mathematical Models and Methods in Applied Science 12, 407–430.
Humphrey, J. D., 1999. Remodeling of collagenous tissue at fixed lengths. J Biomech Eng
121, 591–597.
Humphrey, J. D., 2001. Stress, strain, and mechanotransduction in cells. J Biomech Eng 123,
638–641.
Humphrey, J. D., 2002. Cardiovascular Solid Mechanics: Cells, Tissues, and Organs.
Springer-Verlag, New York.
Humphrey, J. D., 2008. Vascular adaptation and mechanical homeostatsis at tissue, cellular,
and sub-cellular levels. Cell Biochemistry and Biophysics 50, 53–78.
Humphrey, J. D., Na, S., 2002. Elastodynamics and arterial wall stress. Ann Biomed Eng
30 (4), 509–523.
Humphrey, J. D., Taylor, C. A., 2008. Intracranial and abdominal aortic aneurysms: Similarities, differences, and need for a new class of computational models. Ann Rev Biomed
Eng 10, 1–26.
Humphrey, J. D., Wells, P., Baek, S., Hu, J.-J., McLeroy, K., Yeh, A., 2008. A theoreticallymotivated biaxial tissue culture system with intravital microscopy. Biomech Model
Mechanobiol 7, 323–334.
Inzoli, F., Boschetti, F., , M. Z., Longo, T., Fumero, R., 1993. Biomechanical factors in
abdominal aortic aneurysm rupture. Eur J Vasc Surg 7, 667–674.
Jackson, Z. S., Dajnoweiec, D., Gotlieb, A. I., Langille, B. L., 2005. Partial off-loading
of longitudinal tension induces arterial tortuosity. Arterioscler Thromb Vasc Biol 25 (5),
957–962.
Kamiya, A., Togawa, T., 1980. Adaptive regulation of wall shear-stress to flow change in the
canine carotid-artery. Am J Physiol 239 (1), H14–H21.
Karnik, S. K., Brooke, B. S., Bayes-Genis, A., Sorensen, L., Wythe, J. D., Schwartz, R. S.,
Keating, M. T., Li, D. Y., 2003. A critical role for elastin signaling in vascular morphogenesis and disease. Development 130 (2), 411–423.
Kassab, G. S., 2008. Mechanical homeostasis of cardiovascular tissue. In: Artmann, G. M.,
Chien, S. (Eds.), Bioengineering in Cell and Tissue Research. Springer, pp. 371–391.

139

Kim, J., Baek, S., 2011. Circumferential variations of mechanical behavior of the porcine
thoracic aorta during the inflation test. J Biomech, (submitted).
Kroon, M., Holzapfel, G. A., 2007. A model for saccular cerebral aneurysm growth by collagen fibre remodelling. J Theor Biol 247 (4), 775–787.
Kroon, M., Holzapfel, G. A., 2008a. Elastic properties of anisotropic vascular membranes
examined by inverse analysis. Comput Methods Appl Mech Engrg 198, 3622–3632.
Kroon, M., Holzapfel, G. A., 2008b. Estimation of the distributions of anisotropic, elastic
properties and wall stresses of saccular cerebral aneurysms by inverse analysis. Proc R Soc
464, 807–825.
Kroon, M., Holzapfel, G. A., 2009a. A theoretical model for fibroblast-controlled growth of
saccular cerebral aneurysms. J Theor Biol 257 (1), 73–83.
Kroon, M., Holzapfel, G. A., 2009b. A theoretical model for fibroblast-controlled growth of
saccular cerebral aneurysms. J Theor Biol 257 (1), 73–83.
Kuhl, E., Maas, R., Himpel, G., Menzel, A., 2007. Computational modeling of arterial wall
growth. attempts towards patient-specific simulations based on computer tomography.
Biomech Model Mechanobiol 6 (5), 321–331.
Kyriacou, S. K., Schwab, C., Humphrey, J. D., 1996. Finite element analysis of nonlinear
orthotropic hyperelastic membranes. Comput Mech 18 (4), 269–278.
Lagarias, J. C., Reeds, J. A., Wright, M. H., Wright, P. E., 1998. Convergence properties of
the nelder-mead simplex method in low dimensions. SIAM J Optim 9 (1), 112–147.
Langille, B. L., 1993. Remodeling of developing and mature arteries: endotheium, smooth
muscle, and matrix. Journal of Cardiovascular Pharmacology 21, 11–17.
Langille, B. L., Bendeck, M. P., Keeley, F. W., 1989. Adaptation of carotid arteries of young
and mature rabbits to reduced carotid blood-flow. Am J Physiol 256 (4), H931–H939.
Lanir, Y., 1979. A structural theory for the homogeneous biaxial stress-strain relationships
in flat collagenous tissues. J Biomech 12, 423–436.
Lanir, Y., 1983. Constitutive equations for fibrous connective tissues. J Biomech 16, 1–12.
Lanir, Y., Lichtenstein, O., Imanuel, O., 1996. Optimal design of biaxial tests for structural
material characterization of flat tissues. ASME Journal of Biomechanical Engineering 118,
41–47.
Leung, D. Y. M., Glagov, S., Mathews, M. B., 1976. Cyclic stretching stimulates synthesis of
matrix components by arterial smooth-muscle cells invitro. Science 191 (4226), 475–477.
Li, D., Robertson, A. N., 2009. A structural multi-mechanism damage model for cerebral
arterial tissue. Journal of Biomechanical Engineering in print.
140

Li, D. Y., Brooke, B., Davis, E. C., Mecham, R. P., Boak, L. K. S. B. B., Eichwald, E.,
Keating, M. T., 1998. Elastin is an essential determinant of arterial morphogenesis. Nature
393 (6682), 276–280.
Lillie, M. A., Gosline, J. M., 2007. Mechanical properties of elastin along the thoracic aorta
in the pig. J Biomech 40 (10), 2214–2221.
Lindeman, J. H. N., Ashcroft, B. A., Beenakker, J. M., van Es, M., Koekkoek, N. B. R.,
Prins, F. A., Tielemans, J. F., Abdul-Hussien, H., Bank, R. A., Oosterkamp, T. H.,
2010. Distinct defects in collagen microarchitecture underlie vessel-wall failure in advanced
abdominal aneurysms and aneurysms in marfan syndrome. Proc Natl Acad Sci U S A
107 (2), 862–865.
L´pez-Candales, A., Holmes, D. R., Liao, S., Scott, M. J., Wickline, S. A., Thompson,
o
R. W., 1997. Decreased vascular smooth muscle cell density in medial degeneration of
human abdominal aortic aneurysms. Am J Pathol 150 (3), 993–1007.
Lopez-Candales, A., Holmes, D. R., Liao, S., Scott, M. J., Wickline, S. A., Thompson,
R. W., 1997. Decreased vascular smooth muscle cell density in medial degeneration of
human abdominal aortic aneurysms. Am J Pathol 150, 993–1007.
Lu, J., Zhou, X., Raghavan, M. L., 2007. Inverse elastostatic stress analysis in pre-deformed
biological structures: Demonstration using abdominal aortic aneurysms. J Biomech 40,
693–696.
Lu, J., Zhou, X., Raghavan, M. L., 2008. Inverse method of stress analysis for cerebral
aneurysms. Biomech Model Mechanobiol 7, 477–486.
Lu, X., Zhao, J. B., Wang, G. R., Gregersen, H., Kassab, G. S., 2001. Remodeling of the
zero-stress state of femoral arteries in response to flow overload. Am J Physiol Heart Circ
Physiol 280 (4), H1547–H1559.
Maier, A., Essler, M., Gee, M. W., Eckstein, H. H., Wall, W. A., Reeps, C., 2011. Correlation
of biomechanics to tissue reaction in aortic aneurysm assessed by finite elements and 18Ffluordeoxyglucose-PET/CT. Int J Num Meth Biomed Engng (submitted).
Martufi, G., Di Martino, E. S., Amon, C. H., Muluk, S. C., Finol, E. A., 2009. Threedimensional geometrical characterization of abdominal aortic aneurysms: image-based
wall thickness distribution. J Biomech Eng 131, No. 061015.
Masson, I., Boutouyrie, P., Laurent, S., Humphrey, J. D., Zidi, M., 2008. Characterization
of arterial wall mechanical behavior and stresses from human clinical data. Journal of
Biomechanics 41, 2618–2627.
Matsumoto, T., Hayashi, K., 1994. Mechanical and dimensional adaptation of rat aorta to
hypertension. J Biomech Eng 116 (3), 278–283.
McKnight, N. L., Frangos, J. A., 2003. Strain rate mechanotransduction in aligned human
vascular smooth muscle cells. Ann Biomed Eng 31 (3), 239–249.
141

Menashi, S., Campa, J. S., Greenhalgh, R. M., Powell, J. T., 1987. Collagen in abdominal
aortic aneurysm: Typing, content, and degradation. J Vasc Surg 6, 578–582.
Menzel, A., 2005. Modelling of anisotropic growth in biological tissues. a new approach and
computational aspects. Biomech Model Mechanobiol 3 (3), 147–171.
Middleton, R. K., Lloyd, G. M., Bown, M. J., Cooper, N. J., London, N. J., Sayers, R. D.,
2007. The pro-inflammatory and chemotactic cytokine microenvironment of the abdominal
aortic aneurysm wall: A protein array study. J Vasc Surg 45 (3), 574–580.
Miller, F. J., 2002. Aortic aneurysms: It’s all about the stress. Arterioscler Thromb Vasc
Biol 22 (12), 1948–1949.
Mohan, S., Mohan, N., Valente, A. J., Sprague, E. A., 1999. Regulation of low shear flowinduced HAEC VCAM-1 expression and monocyte adhesion. Am J Physiol Cell Physiol
276, c1100–7.
Mulvany, M. J., 1992. Vascular growth in hypertension. J Cardiovasc Pharmacol 20, S17–11.
Myung, I. J., 2000. The importance of complexity in model selection. Journal of Mathematical Psychology 44, 190–204.
Naik, P. A., Shi, P., Tsai, C., 2007. Extending the akaike information criterion to mixture
regression models. Journal of the American Statistical Association 102, 244–254.
Naschitz, J. E., Lenger, R., 2008. Why traumatic leg amputees are at increased risk for
cardiovascular diseases. QJM 101 (4), 251–259.
Nelder, J. A., Mead, R., 1965. A simplex-method for function minimization. Computer
Journal 7 (4), 308–313.
Nevitt, M. P., Ballard, D. J., Hallett, J. W., 1989. Prognosis of abdominal aortic aneurysms.
a population-based study. N Engl J Med 321 (15), 1009–1014.
Nissen, R., Cardinale, G. J., Udenfriend, S., 1978. Increased turnover of arterial collagen in
hypertensive rats. Proc Natl Acad Sci 75 (1), 451–453.
O’Callaghan, C. J., William, B., 2000. Mechanical strain induced extracellular matrix production by human vascular smooth muscle cells. Hypertension 36 (3), 319–324.
Ogden, R. W., Saccomandi, G., 2007. Introducing mesoscopic information into constitutive
equations for arterial walls. Biomech Model Mechanobiol 6, 333–344.
O’Rourke, M. F., Hashimoto, J., 2007. Mechanical factors in arterial aging: A clinical perspective. Journal of the American College of Cardiology 50, 1–13.
Oxlund, H., Rasmussen, L. M., Andreassen, T. T., Heickendorff, L., 1989. Increased aortic
stiffness in patients with type 1 (insulin-dependent) diabetes mellitus. Diabetologia 32,
748–752.
142

Paes, E. H., Vollmar, J. F., Pauschinger, P., Mutschler, W., Henze, E., Friesch, A., 1990. Late
vascular damage after unilateral leg amputation. Z Unfallchir Versicherungsmed 83 (4),
227–236.
Pandit, A., Lu, X., Wang, C., Kassab, G. S., 2005. Biaxial elastic material properties of
porcine coronary media and adventitia. American Journal of Physiology-Heart and Circulatory Physiology 288, 2581–2587.
Park, H. C., Youn, S., 1998. Finite element analysis and constitutive modelling of anisotropic
nonlinear hyperelastic bodies with convected frames. Comput Methods Appl Mech Eng
151, 605–618.
Patel, M. I., Melrose, J., Ghosh, P., Appleberg, M., 1996. Increased synthesis of matrix
metalloproteinases by aortic smooth muscle cells is implicated in the etiopathogenesis of
abdominal aortic aneurysms. J Vasc Surg 24 (1), 82–92.
Pearce, W. H., Shively, V. P., 2006. Abdominal aortic aneurysm as a complex multifactorial
disease: interaction of polymorphisms of inflammatory genes, features of autoimmunity,
and current status of MMPs. Ann New York Acad Sci 1085, 117–132.
Powell, J. T., 2002. An introduction to vascular biology, 2nd Edition. Cambridge University
Press, Cambridge, Ch. Abdominal aortic aneurysm.
Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling, W. T., et al., 1992. Numerical
Recipes. Cambridge University Press.
Raghavan, M., Ma, M., Fillinger, M., 2006. Non-invasis determination of zero-pressure geometry of arterial aneurysms. Annals of Biomedical Engineering 34, 1414–1419.
Raghavan, M. L., Vorp, D. A., 2000. Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model
and evaluation of its applicability. J Biomech 33 (4), 475–482.
Raghavan, M. L., Vorp, D. A., Federle, M., Makaroun, M. S., Webster, M. W., 2000. Wall
stress distribution on three dimensionally reconstructed models of human abdominal aortic
aneurysm. J Vasc Surg 31, 760–769.
Rezakhaniha, R., Fonck, E., Genoud, C., Stergiopulos, N., 2011. Role of elastin anisotropy
in structural strain energy functions of arterial tissue. Biomech Model Mechanobiol 10,
599–611.
Rissland, P., Alemu, Y., Einav, S., Ricotta, J., Bluestein, D., 2009. Abdominal aortic
aneurysm risk of rupture: Patient-specific FSI simulations using anisotropic model. J
Biomech Eng 131, 031001.
Rizzo, R. J., McCarthy, W. J., Dixit, S. N., Lilly, M. P., Shively, V. P., Flinn, W. R., Yao,
J. S. T., 1989. Collagen types and matrix protein content in human abdominal aortic
aneurysms. J Vasc Surg 10 (4), 365–373.

143

Sacks, M. S., Vorp, D. A., Raghavan, M. L., Federle, M. P., Webster, M. W., 1999. In vivo
three-dimensional surface geometry of abdominal aortic aneurysms. Ann Biomed Eng 27,
469–479.
Saravanan, U., Baek, S., Rajagopal, K. R., Humphrey, J. D., 2006. On the deformation
of the circumflex coronary artery during inflation tests at constant length. Experimental
Mechanics 46, 647–656.
Schulze-Bauer, C. A. J., Holzapfel, G. A., 2003. Determination of constitutive equations for
human arteries from clinical data. Journal of Biomechanics 36, 165–169.
Schwetlick, H., Tiller, V., 1985. Numerical methods for estimating parameters in nonlinear
models with errors in the variables. Technometrics 27, 17–24.
Scotti, C. M., Shkolnik, A. D., Muluk, S. C., Finol, E. A., 2005. Fluid-structure interaction
in abdominal aortic aneurysms: effects of asymmetry and wall thickness. Biomed Eng
online 4, 64:1–22.
Sheidaei, A., Hunley, S., Zeinali-Davarani, S., Raguin, L., Baek, S., 2011. Simulation of
abdominal aortic aneurysm growth with updating hemodynamic loads using a realistic
geometry. Medical Engineering & Physics 33, 80–88.
Shimizu, K., Mitchell, R. N., Libby, P., 2006. Inflammation and cellular immune responses
in abdominal aortic aneurysms. Arteriscler Thromb Vasc Biol 26 (5), 987–994.
Sho, E., Sho, M., Hoshina, K., Kimura, H., Nakahashi, T. K., Dalman, R. L., 2004. Hemodynamic forces regulate mural macrophage infiltration in experimental aortic aneurysms.
Exp Mol Pathol 76 (2), 108–116.
Shum, J., Di Martino, E. S., Goldhammer, A., Goldman, D. H., Acker, L. C., Patel, G.,
Ng, J. H., Martufi, G., Finol, E. A., 2010. Semiautomatic vessel wall detection and quantification of wall thickness in computed tomography images of human abdominal aortic
aneurysms. Med Phys 37, 638–648.
Shum, J., Martufi, G., Martino, E. D., Washington, C. B., Grisafi, J., Muluk, S. C., Finol,
E. A., 2011. Quantitative assessment of abdominal aortic aneurysm geometry. Ann Biomed
Eng 39 (1), 277–286.
Speelman, L., Bohra, A., Boosman, E. M. H., Schurink, G. H. W., van de Vosse, F. N.,
Makaroun, M. S., Vorp, D. A., 2007. Effects of wall calcifications in patient-specific wall
stress analyses of abdominal aortic aneurysms. J Biomech Eng 129 (1), 105–109.
Sumner, D. S., Hokanson, D. E., Strandness, D. E., 1970. Stress-strain characteristics and
collagen-elastin content of abdominal aortic aneurysms. Surg Gynecol Obstet 130 (3),
459–466.
Sumpio, B. E., Banes, A. J., Link, W. G., Johnson, G. J., 1988. Enhanced collagen production
by smooth muscle cells during repetitive mechanical stretching. Arch Surg 123, 1233–1236.

144

Taylor, C. A., Humphrey, J. D., 2009. Open problems in computational vascular biomechanics: Hemodynamics and arterial wall mechanics. Comput Methods Appl Mech Eng 198,
3514–3523.
Thompson, R. W., Liao, S. X., Curci, J. A., 1997. Vascular smooth muscle cell apoptosis in
abdominal aortic aneurysms. Coron Artery Dis 8, 623–631.
Thubrikar, M. J., Robicsek, F., Labrosse, M., Chervenkoff, V., Fowler, B. L., 2003. Effect of
thrombus on abdominal aortic aneurysm wall dilation and stress. J Cardiovasc Surg 44,
67–77.
Torczon, V. J., 1989. Multi-directional search: A direct search algorithm for parallel machines. Ph.D. thesis, Rice University.
Truesdell, C., Noll, W., 1965. The nonlinear field theories of mechanics. Springer-Verlag.
Valent´ A., Cardamone, L., Baek, S., Humphrey, J. D., 2009a. Complementary vasoactivın,
ity and matrix remodeling in arterial adaptations to altered flow and pressure. J R Soc
Interface 6, 293–306.
Valent´ A., Cardamone, L., Baek, S., Humphrey, J. D., 2009b. Complementary vasoactivın,
ity and matrix remodeling in arterial adaptations to altered flow and pressure. J R Soc
Interface 6 (32), 293–306.
Vande Geest, J. P., Sacks, M. S., Vorp, D. A., 2004. Age dependency of the biaxial biomechanical behavior of human abdominal aorta. J Biomech Eng 126, 815–822.
Vande Geest, J. P., Sacks, M. S., Vorp, D. A., 2006. The effects of aneurysm on the biaxial
mechanical behavior of human abdominal aorta. J Biomech 39, 1324–1334.
Vollmar, J. F., Pauschinger, P., Paes, E., Henze, E., Friesch, A., 1989. Aortic aneurysms as
late sequelae of above-knee amputation. Lancet 334 (8667), 834–835.
Vorp, D. A., 2007. Biomechanics of abdominal aortic aneurysm. J Biomech 40.
Vorp, D. A., Lee, P. C., Wang, D. H. J., Makaroun, M. S., Nemoto, E. M., Ogawa, S., ,
Webster, M. W., 2001. Association of intraluminal thrombus in abdominal aortic aneurysm
with local hypoxia and wall weakening. J Vasc Surg 34 (2), 291–299.
Vorp, D. A., Raghavan, M. L., Webster, M. W., 1998. Stress distribution in abdominal aortic
aneurysm: influence of diameter and asymmetry. J Vasc Surg 27, 632–639.
Vorp, D. A., Vande Geest, J. P., 2005. Biomechanical determinants of abdominal aortic
aneurysm rupture. Arterioscler Thromb Vasc Biol 25 (8), 1558–1566.
Wagenmakers, E. J., Farrell, S., 2004. AIC model selection using akaike weights. Psychonomic
Bulletin and Review 11, 192–196.

145

Walpola, P. L., Gotlieb, A. I., Langille, B. L., 1993. Monocyte adhesion and changes in
endothelial-cell number, morphology, and F-actin distribution elicited by low shear-stress
in vivo. Am J Pathol 142, 1392–1400.
Wan, W., Hansen, L., Gleason, R. L., 2010. A 3-d constrained mixture model for mechanically
mediated vascular growth and remodeling. Biomech Model Mechanobiol 9, 403–419.
Wang, C., Garcia, M., Lu, X., Lanir, Y., Kassab, G. S., 2006. Three-dimensional mechanical
properties of porcine coronary arteries: a validated two-layer model. American Journal of
Physiology-Heart and Circulatory Physiology 291, 1200–1209.
Wang, D. H. J., Makaroun, M. S., Webster, M. W., Vorp, D. A., 2002. Effect of intraluminal
thrombus on wall stress in patient-specific models of abdominal aortic aneurysm. J Vasc
Surg 36, 598–604.
Watton, P., Hill, N., Heil, M., 2004. A mathematical model for the growth of the abdominal
aortic aneurysm. Biomechanics and Modeling in Mechanobiology 3, 98–113.
Watton, P. N., Hill, N. A., 2009. Evolving mechanical properties of a model of abdominal
aortic aneurysm. Biomech Model Mechanobiol 8 (1), 25–42.
Wolinsky, H., 1971. Effects of hypertension and its reversal on thoracic aorta of male and
female rats - morphological and chemical studies. Circ Res 28, 622–637.
Wright, M. H., 1996. Direct search methods: once scorned, now respectable. In: Griffiths,
D. F., Watson, G. A. (Eds.), Numerical Analysis 1995. Addison-Wesley Longman, Reading,
MA, pp. 191–208.
Wulandana, R., Robertson, A. M., 2005. An inelastic multi-mechanism constitutive equation
for cerebral aretrial tissue. Biomech Model Mechanobiol 4 (4), 235–248.
Xu, C. P., Zarins, C. K., Bassiouny, H. S., Briggs, W. H., Reardon, C., Glagov, S., 2000.
Differential transmural distribution of gene expression for collagen types i and iii proximal
to aortic coarctation in the rabbit. J Vasc Res 37 (3), 170–182.
Zeinali-Davarani, S., Sheidaei, A., Baek, S., 2011. A finite element model of stressmediated vascular adaptation: application to abdominal aortic aneurysms. Comput Methods Biomech Biomed Eng 14 (9), 803–817.
Zhang, J., Schmidt, J., Ryschich, E., H. Schumacher, J. R. A., 2003. Increased apoptosis and
decreased density of medial smooth muscle cells in human abdominal aortic aneurysms.
Chin Med J 116 (10), 1549–1552.
Zhou, X., Lu, J., 2009. Estimation of vascular open configuration using finite element inverse
elastostatic method. Engineering with Computers 25, 49–59.
Zhou, X., Raghavan, M. L., Harbaugh, R. E., Lu, J., 2010. Specific wall stress analysis in
cerebral aneurysms using inverse shell model. Ann Biomed Eng 38 (2), 478–489.

146

Zulliger, M. A., Frideza, P., Hayashib, K., Stergiopulos, N., 2004a. A strain energy function
for arteries accounting for wall composition and structure. J Biomech 37, 989–1000.
Zulliger, M. A., Rachev, A., Stergiopulos, N., 2004b. A constitutive formulation of arterial
mechanics including vascular smooth muscle tone. Am J Physiol Heart Circ Physiol 287 (3),
H1335–H1343.

147