Quantification of Hemodynamics in Abdominal Aortic Aneurysms During Rest and Exercise Using ...

26 Pages·2010·2.33 MB·English

AnnalsofBiomedicalEngineering,Vol.38,No.4,April2010((cid:1)2010)pp.1288–1313
DOI:10.1007/s10439-010-9949-x
Quantification of Hemodynamics in Abdominal Aortic Aneurysms During
Rest and Exercise Using Magnetic Resonance Imaging and Computational
Fluid Dynamics
ANDREA S. LES,1 SHAWN C. SHADDEN,2 C. ALBERTO FIGUEROA,1 JINHA M. PARK,3 MAUREEN M. TEDESCO,4
ROBERT J. HERFKENS,5 RONALD L. DALMAN,4 and CHARLES A. TAYLOR1,6
1Departmentof Bioengineering,Stanford University,Stanford,CA, USA; 2Departmentof Mechanical andAerospace
Engineering,IllinoisInstituteofTechnology,Chicago,IL,USA;3DepartmentofRadiology,UniversityofSouthernCalifornia,
LosAngeles, CA, USA; 4Division of Vascular Surgery, StanfordUniversity, Stanford,CA, USA; 5Department ofRadiology,
StanfordUniversity, Stanford,CA, USA; and6JamesH. ClarkCenter, RoomE350B, 318Campus Drive,Stanford,
CA94305-5431, USA
(Received20February2009;accepted26January2010;publishedonline9February2010)
AssociateEditorKennethR.Lutchenoversawthereviewofthisarticle.
Abstract—Abdominalaorticaneurysms(AAAs)affect5–7% ABBREVIATIONS
of older Americans. We hypothesize that exercise may slow
AAAgrowthbydecreasinginflammatoryburden,peripheral AAA Abdominal aortic aneurysm
resistance,andadversehemodynamicconditionssuchaslow, DBP Diastolic blood pressure
oscillatory shear stress. In this study, we use magnetic
IR Infrarenal
resonance imaging and computational fluid dynamics to
Mid-An Mid-aneurysm
describe hemodynamics in eight AAAs during rest and
exercise using patient-specific geometric models, flow wave- MRI Magnetic resonance imaging
forms,andpressuresaswellasappropriatelyresolvedfinite- MWSS Mean wall shear stress
element meshes. We report mean wall shear stress (MWSS) OSI Oscillatory shear index
and oscillatory shear index (OSI) at four aortic locations
SB Suprabifurcation
(supraceliac, infrarenal, mid-aneurysm, and suprabifurca-
SBP Systolic blood pressure
tion) and turbulent kinetic energy over the entire computa-
tional domain on meshes containing more than an order of SC Supraceliac
magnitude more elements than previously reported results TKE Turbulent kinetic energy
(mean:9.0-millionelements;SD:2.3M;range:5.7–12.0M).
MWSS was lowest in the aneurysm during rest 2.5dyn/cm2
(SD: 2.1; range: 0.9–6.5), and MWSS increased and OSI
INTRODUCTION
decreased at all four locations during exercise. Mild turbu-
lence existed at rest, while moderate aneurysmal turbulence
was present during exercise. During both rest and exercise, Abdominal aortic aneurysms (AAAs), defined as a
aorticturbulencewasvirtuallyzerosuperiortotheAAAfor 1.5-fold or more enlargement of the abdominal aorta,
seven out of eight patients. We postulate that the increased affect 5–7% of Americans over age 60 and kill
MWSS, decreased OSI, and moderate turbulence present approximately 9,000 people each year.17,18,49 Risk
duringexercise mayattenuate AAA growth.
factors for AAA include male gender, advanced age,
cigarette smoking, heritable predisposition, and ath-
Keywords—Turbulence, Mean wall shear stress, Oscillatory
erosclerotic disease.31 Because no effective medical
shear index, Mesh independence, Flow waveforms, Blood
therapy exists for AAA patients, and because on
pressure, Windkessel boundarycondition, Patient-specific.
average AAAs enlarge at a predictable rate,8 most
AAApatientspostponetreatmentuntiltheiraneurysm
diameter reaches the threshold for surgical repair:
4.5 cmforwomenand5.5 cmformen.9Thisthreshold
represents the point at which the risk of rupture is
thought to outweigh the risk of surgical intervention.
Address correspondence to Charles A. Taylor, James H. Clark
Inflammation, increased peripheral resistance, and
Center,RoomE350B,318CampusDrive,Stanford,CA94305-5431,
USA.Electronicmail:[email protected] adverse/pro-inflammatory localized hemodynamic
1288
0090-6964/10/0400-1288/0(cid:1)2010BiomedicalEngineeringSociety
Hemodynamics in AAAsDuringRest andExercise 1289
conditions, including low and oscillatory wall shear randomized to an exercise or usual-activity control
stress, are all thought to play a critical role in AAA cohort and followed for 3 years.13
formation. On the microscale, inflammatory processes While gene and protein expression analysis and
such as recruitment of leukocytes into the aortic wall population-based studies provide insight into the
mediated by elastin degradation products19 and inflammatory and peripheral resistance theories of
chemokines,29 phagocyte infiltration resulting in AAA development, the use of noninvasive medical
pro-inflammatory cytokines40 and reactive oxygen imaging and computational fluid dynamics methods
species,33 and the production of matrix metallopro- enable an in-depth study of the hemodynamic condi-
teinases such as MMP-9 by macrophages59 are all tionsspecifictolargearteries,includingtheabdominal
thought to instigate or promote AAA formation and aorta.55 While important progress has been made in
growth.12 On the macroscale, patients with increased modeling blood flow in arteries over the last decade,
peripheral resistance caused by chronic spinal cord much of this study has utilized zero-pressure, or
injury,67severeperipheralvasculardisease,1andabove prescribed pressure or flow boundary conditions that
thekneeamputation63havebeenshowntohavea2.9-, resultinnon-physiologicpressureand/orflowfields.In
4.2-, and 5.3-fold increased rate of AAA, respectively, addition, some studies use severely under-resolved
compared to control populations. Lastly, adverse/pro- meshes, idealized AAA geometry, or neglect the
inflammatory hemodynamic conditions such as low nearby mesenteric, renal, and iliac arteries that greatly
andoscillatorywallshearstress,aswellasthepresence affect the flow entering and exiting the AAA. Addi-
of recirculation zones are also hypothesized to aug- tionally,comparativelylittleprogresshasbeenmadein
ment AAA development. In the human abdominal usingthisunderstandingofhemodynamicmechanisms
aorta, a correlation between low and/or oscillatory to guide the development of an effective medical
wall shear stress and atherosclerosis, which shares therapy for AAAs.
many of the same disease mechanisms as AAAs, has In this article, we use magnetic resonance imaging
been demonstrated by autopsy,11 as well as by exper- (MRI) and computational fluid dynamics (CFD) to
imental36andcomputationalmethods.53Furthermore, study hemodynamics under rest and exercise condi-
a correlation between regions of low wall shear stress tions in eight AAApatients. The useof MRI provides
and growth of intracranial aneurysms has recently a means to model blood flow in a highly patient-spe-
beenelucidated.7Inaddition,ithasbeendemonstrated cificmannerbecauseitenablesthenearlysimultaneous
thatregionsofslowand/orrecirculatingflowcorrelate acquisition of anatomic and physiologic flow data,
with regions of thrombus deposition in basilar aneu- whereas the use of computational simulations enables
rysms,45 and a similar relationship has been proposed the study of blood flow at resolution many times finer
for abdominal aortic aneurysms.5,48 A more complete than any imaging modality and can capture transient
review of biomechanical factors in aneurysm forma- phenomena unobservable by cardiac-gated imaging
tioncanbefoundinarecentarticlebyHumphreyand sequences. Our hemodynamic analyses include quan-
Taylor.25 tification of pressure and velocity, normal stresses
WehaveproposedanewmedicaltherapyforAAAs acting on the luminal surface of the AAAs, mean wall
whereby lower limb exercise may slow or halt the shear stress (MWSS), oscillatory shear index (OSI),
growth of AAAs by decreasing inflammatory burden, and turbulent kinetic energy (TKE). We aim to
peripheral resistance, and adverse hemodynamic con- quantitatively and qualitatively capture the baseline
ditions in the abdominal aorta.13 Aside from reducing hemodynamic conditions in AAAs using physiologic
all-causemortality,38exercisehasbeenshowntoeffect boundary conditions, appropriately resolved meshes,
potentially salutary increases in shear stress in the andpatient-specificgeometricmodels,flowwaveforms,
abdominal aortas ofboth young53,54and old people.10 and pressures. Additionally, we aim to examine some
Furthermore, Nakahasi et al. and Hoshina et al. of the mechanisms by which exercise might slow the
showed that increasing infrarenal flow through growth of AAAs in humans in order to guide and
experimentallycreatedAAAsviaadownstreamarterio- evaluate the effectiveness of exercise as a medical
venous fistula increased antioxidative enzyme expres- therapy for AAA disease.
sion, restricted reactive oxygen species production,39
and limited AAA enlargement in rats.23,39 Despite
these promising studies, it is not known whether
METHODS
exercise may slow AAA growth in humans. This
questionis thebasis of aNationalInstitutes ofHealth
Magnetic Resonance Imaging
(NIH)-sponsored randomized study examining the
effects of a prescribed exercise regimen on a group Sevenmaleandonefemalepatientswithameanage
of AAA patients. In this study, AAA patients are of 73.3 years (SD: 5.70; range: 66–84) with known
1290 LES etal.
small AAAs (<5 cm) were imaged in the supine po- Image Processing
sitionusinga1.5TGESignaMRscanner(GEMedical
First, gradient nonlinearities arising during acqui-
Systems, Milwaukee WI) with an 8-channel cardiac
sition were corrected in the MRAs to avoid image
coil. Imaging studies were conductedunder a protocol
distortion in the slice direction of the magnetic
approved by the institutional review board, informed
field.14 Next, three-dimensional patient-specific geo-
consentwasobtainedfromallsubjects,andallpatients
metric models were constructed from the MRA data
were screened for contraindications to MR and gado-
usingcustomsoftware.65Themodelsbeganatthelevel
linium usage, including renal insufficiency.60 A 3D
of the diaphragm and included the aorta, the hepatic
gadolinium-enhancedmagneticresonanceangiography
and splenic arteries, the superior mesenteric artery
(MRA) sequence was used to image the lumen of the
(SMA), the right and left renal arteries, and the right
aorta and nearby branch vessels. Imaging parameters
and left internal and external iliac arteries. Approxi-
included: 512 9 192 acquisition matrix (reconstructed
mate centerline paths through the vessels of interest
to 512 9 512), a 40 cm2 square field of view, 3.0–
were defined, and the image data were resampled in
3.3 ms TR (repetition time), 0.7–0.8 ms TE (echo
planes perpendicular to the centerline paths. Vessel
time), and 25(cid:2) flip angle. A 3-mm slice thickness was
boundariesweresegmentedusingalevel-setmethod.65
used, 124 slices were included in each volume, and
Solid models were lofted for each vessel and then
all MRAs were acquired coronally. A cardiac-gated
unionedtocreateasinglesolidmodelrepresentingthe
1-component phase contrast sequence (PC-MRI) was
flowdomain(Fig. 1).Avisualcomparisonbetweenthe
used to measure blood flow velocity perpendicular
maximum intensity projections (MIPs) of the MRAs
to the aorta at the supraceliac (SC) and infrarenal
andthe3DcomputermodelscanbeseeninFig. 2.For
(IR) levels. Imaging parameters included: 256 9 192
each patient, the PC-MRI images acquired at the SC
acquisition matrix (reconstructed to 256 9 256), a
and IR levels were segmented at each of the 24 time
24 9 24 to 34 9 34 cm2 square field of view, 5 mm
points using a level-set technique, baseline corrected
slicethickness,11.9–12.7 msTR(repetitiontime),4.5–
using a linear correction algorithm based on the aver-
5.3 ms TE (echo time), 20(cid:2) flip angle, and two signal
age background noise, and then assembled into volu-
averages. A velocity encoding gradient of 150 cm/s
metric flow waveforms.
through-plane was used, and 24 time points were
reconstructed. The temporal resolution of the velocity
data equaled two times the TR (23.8–25.4 ms).
Computational Fluid Dynamics Analysis
Immediately after the MR exam, a brachial systolic
blood pressure (SBP) and a diastolic blood pressure The Navier–Stokes and continuity equations repre-
(DBP)measurementwereacquiredusinganautomatic senting conservation of momentum and mass of a
pressure cuff (Omron Healthcare Inc., Bannockburn, Newtonianincompressiblefluidwereusedtorepresent
Illinois). blood flow in an arbitrary numerical domain X with
1. MRA 2. Pathlines 3. Segmentation 4. Lofting 5. Solid Model 6. Mesh
FIGURE1. BuildingamodelfromtheMRAdata(1)entailsgeneratingpathlines(2),segmentingalongthepathlines(3),loftingthe
segmentationstocreatevessels(4),unioningthevesselstocreateasingle,solidmodelrepresentingtheflowdomain(5),andthen
discretizingthemodelintoafinite-elementmesh(6).
Hemodynamics in AAAsDuringRest andExercise 1291
FIGURE2. Themaximumintensityprojection(MIP)oftheMRAiscomparedtothe3Dcomputermodelforalleightpatients.
boundaryC.TheboundaryCissuchthatC¼C [C In order to define the variational form correspond-
i w
[C ; C \C \C ¼[; where C represents an inlet ing to the Eqs. (1)–(5), we used the following trial
h i w h i
boundary where a Dirichlet condition is prescribed on solution and weighting function spaces:
thevelocityfield;Cwrepresentsthevesselwallboundary; S ¼n*u(cid:7)(cid:7)*uð(cid:1);tÞ2H1ðXÞnsd; t2½0;T(cid:3); *uð(cid:1);tÞ¼*uin
and C represents all the outflow boundaries where a
h
Neumann condition is prescribed. A well-posed Initial onC; *uð(cid:1);tÞ¼*0onC o
i w
Boundary Value Problem is defined by the aforemen-
tioned partial differential equations subject to suitable W ¼nw*(cid:7)(cid:7)w*ð(cid:1);tÞ2H1ðXÞnsd; t2½0;T(cid:3); w*ð(cid:1);tÞ¼*0
(cid:7)
boundaryandinitialconditions,viz.:
q*u;tþqu*(cid:1)ru*¼(cid:2)rpþdivðs~Þþ*f9>= inX onPC¼i;(cid:8)w*pðj(cid:1)p;tðÞ(cid:1);¼tÞ*02oHn1CðXwÞo; t2½0;T(cid:3)(cid:9) ð6Þ
divð*uÞ¼0 >; ð1Þ where H1 represents the usual Sobolev space of func-
s ¼2lD with D ¼1(cid:1)ru*þr*uT(cid:2) tionswithsquare-integrablevaluesandfirstderivatives
~ ~ ~ 2 in X and nsd represents the number of spatial dimen-
sions.Consideringthis,thevariationalformulationfor
*uðx*;tÞ¼*uinðx*;tÞ x* 2C ð2Þ this problem can be written as:
i
B ðw*;q;*u;pÞ¼0
G
*uðx*;tÞ¼*0 x* 2Cw ð3Þ BGðw*;q;*u;pÞ¼Z nw*(cid:1)(cid:1)qu*;tþq*u(cid:1)ru*(cid:2)*f(cid:2)
X
(cid:3) (cid:4) (cid:1) (cid:2)o Z
*t (cid:1)x*;t(cid:2)¼ (cid:2)pI þs n*¼*h(cid:1)x*;t(cid:2) x* 2C ð4Þ þrw* : (cid:2)pIþs dx*(cid:2) rq(cid:1)u*dx*
*n ~ ~ h Z ~ ~ XZ
(cid:2) w*(cid:1)ð(cid:2)pIþsÞ(cid:1)*nds þ qu*(cid:1)n*ds
*uðx*;0Þ¼*u0ðx*Þ x* 2X ð5Þ Ch ~ ~ C
ð7Þ
Here, the primary variables are the blood velocity
*u¼(cid:5)u ;u ;u (cid:6)andthepressurep.Theblooddensityq We adoptedthe Coupled-Multidomainformulation
x y z described in Vignon-Clementel et al.61 This method is
and viscosity l areassumed constant (q = 1.06 g/cm3,
* basedonavariationalmultiscaledecompositionofthe
l = 0.04 P). The external body force f representing
domainXintoanupstream3DnumericaldomainX^and
gravitywasneglectedandthevesselwallboundarywas
adownstreamreduced-orderanalyticaldomainX0;such
assumed to be rigid, as Eq. (3) indicates. Please note
that*hinEq.(4)representsthetractionvectorimposed that X^ \X0 ¼[ and X^ [X0 ¼X: The numerical
on the outflow boundaries C . andtheanalyticaldomainsareseparatedbyaninterface
h
1292 LES etal.
C . This interface is what we refer to as ‘‘outlet’’ or Note that the solution in the numerical domain X^
B
‘‘outlets’’ (of the 3D computational domain) through- depends on the operators defined by the mathematical
out this article. Furthermore, the method considers a model chosen for the downstream analytical domain,
disjointdecompositionfortheflowvariables:thesolu- expressedasafunctionofthesolutionvariablesofthe
tion vector V* ¼n*u;poTis separated into a component numerical domain n^u*;p^o. For a three-element Wind-
defined within the numerical domain and another kessel (R , C, R ) model of the downstream vascula-
p d
componentdefinedwithintheanalyticaldomain,viz: ture, the operators M and H become62:
V* ¼(cid:10)V^* þV*0(cid:11); with V^*(cid:7)(cid:7)(cid:7) ¼0 and V*0(cid:7)(cid:7)(cid:7) ¼0: Mm(cid:1)^u*;p^(cid:2)¼RpZ ^u*(cid:1)*^nds
(cid:7)X0 (cid:7)X^ ~ CB
ð8Þ Z te(cid:2)ðt(cid:2)t1Þ=sZ ^ ^ ^ ^
þ u*ðt Þ(cid:1)n*dsdt þn*(cid:1)s^(cid:1)n*(cid:2)s^
This decomposition satisfies V*^ ¼V*0 at the interface 0 C CB 1 1 ~ ~
(cid:12) Z (cid:13)
CB,andisalsoappliedtotheweightingfunctionsw* and Hm ¼ P^ð0Þ(cid:2)Rp ^u*ð0Þ(cid:1)*^ndC(cid:2)P^dð0Þ e(cid:2)t=sþP^dðtÞ
q.Undertheseassumptions,Eq.(7)canberewrittenas: ~ CB
Z w^*(cid:1)(cid:1)q*^u þq*^u(cid:1)r^u*(cid:2)*f(cid:2)þrw^* :(cid:12)(cid:2)p^I þs^(cid:13)dx* M*c(cid:1)^u*;p^(cid:2)¼^u*
;t
X^ ~ ~ H* ¼0* ð12Þ
Z (cid:12) (cid:13) Z (cid:12) (cid:13) c
(cid:2) w^*(cid:1) (cid:2)p^I þs^ (cid:1)^n*dsþ w*0(cid:1) (cid:2)p0I þs0 (cid:1)n*0ds
C^h ~ ~ CB ~ ~ where s = RdC. Here, Rp represents the proximal
Z Z Z resistance, C the capacitance, and R the distal resis-
(cid:2) rq^(cid:1)^u*dx* þ q^*^u(cid:1)*^nds (cid:2) q0u*0(cid:1)*n0ds ¼0 ð9Þ tance of the downstream vasculaturde; P is a time-
d
X^ C^ CB varying function representing the terminal capillary
The integral terms defined on C represent the con- pressure. P(0), Q(0), and P (0) are the initial condi-
B d
tinuityoftractionandfluxattheinterfacebetweenthe tions for pressure, flow, and capillary pressure for the
numerical and analytical domains. These terms can be downstream analytical domain, respectively.
approximatedbytheoperatorsM¼nM ;M* o(cid:7)(cid:7) and The variational form defined by Eqs. (11) and (12)
n o(cid:7) ~ m c (cid:7)CB and the boundary and initial conditions (2)–(5) are
H¼ H ;H* (cid:7) : then solved using a custom stabilized finite-element
m c (cid:7)
~ CB solver which utilizes equal-order interpolation for
Z (cid:12) (cid:13) Z velocity and pressure.51,61,64 All simulations were per-
w*0(cid:1) (cid:2)p0I þs0 (cid:1)n*0ds (cid:2) q0u*0(cid:1)*n0ds (cid:4) formed using direct numerical simulations, without
CB ~ ~ CB any assumptions that impose laminar or fully devel-
Z w*0(cid:1)(cid:12)M (cid:1)u*0;p0(cid:2)þH (cid:13)(cid:1)*n0ds ð10Þ oped turbulent regimes on the flow. This is justified
m m
CB ~ ~ because most turbulence models assume fully devel-
(cid:2)Z q0(cid:1)M* (cid:1)u*0;p0(cid:2)þH* (cid:2)(cid:1)*n0ds oped turbulent conditions, which is not the case for
c c pulsatile flow in AAAs: The flow, undoubtedly tur-
CB
bulent during late systole, largely relaminarizes during
TheoperatorsMandHaredefinedonthedomainX0 early systole and diastole.57 Computations were per-
based on the reduced-order model chosen to represent
formed on a 20.6 teraflops supercomputer with a total
flow and pressure in the downstream domain (i.e.,
of 2,208 cores (276 Dell PowerEdge 1950 compute
morphometry-based impedance function, lumped-
nodes, consisting of dual-socket quad-core processors,
parameter model, etc.). The final weak form of the
featuring 16 GB of memory each).
Coupled-Multidomainmethodresultsin61:
Z w^*(cid:1)(cid:1)q*^u þq^u*(cid:1)r*^u(cid:2)*f(cid:2)þrw^* :(cid:12)(cid:2)p^Iþs^(cid:13)dx* Numerical Simulation Details
;t
X^ ~ ~ Below we describe the simulation process for a
Z ^ (cid:12) (cid:13) ^ singlepatient;thesimulationprocesswasrepeatedina
(cid:2) w*(cid:1) (cid:2)p^Iþs^ (cid:1)*nds
patient-specificmannerforallpatientsunderbothrest
C^h ~ ~
Z ^ (cid:12) (cid:1)^ (cid:2) (cid:13) ^ Z ^ and exercise conditions. First, the model was discret-
(cid:2) w*(cid:1) M u*;p^ þH (cid:1)*nds (cid:2) rq^(cid:1)*udx* izedinto aninitialisotropic finite-elementmesh witha
m m
CB ~ ~ X^ maximum edge size of 1.5 mm using a commercial
þZ q^*^u(cid:1)*^nds þZ q^(cid:1)M* (cid:1)^u*;p^(cid:2)þH* (cid:2)(cid:1)^n*ds ¼0 ð11Þ meshing kernel (MeshSimTM, Simmetrix, Troy, NY).
c c
C^ CB Next, steady simulations with mean SC flow as the
Hemodynamics in AAAsDuringRest andExercise 1293
inlet boundary condition and resistance outlet bound- A measured (rest) or extrapolated (exercise) SC
ary conditions were performed on the initial isotropic waveform was decomposed into a steady nonzero
mesh under both rest and exercise conditions for each mean flow which was mapped to a Poiseuille velocity
patient. Resistances for each outlet were equal to the profile,andapulsatilezero-meanflowwavewhichwas
R described later. The results of this steady simu- mapped to a ten-frequency Womersley velocity pro-
total
lation were used to adapt the mesh,37 creating an file.66 The Poiseuille and Womersley velocity profiles
average intermediate anisotropic mesh consisting of were then linearly superimposed to define the inlet
922,525lineartetrahedralelements(SD:238,766;range velocity profile on C. At each outlet, a three-element
i
595.558–1,323,275) for rest, and 967,862 linear Windkessel model consisting of a proximal resistance
tetrahedral elements (SD: 211,180; range 623,630– (R ), capacitance (C), and distal resistance (R ) was
p d
1,219,620)forexercise.Apulsatilesimulationwasthen used to represent the resistance and capacitance of the
performedonthisnewlyadaptedmeshforfourcardiac proximalvessels(thearteries),andtheresistanceofthe
cycles,andthecapacitance,C,andR ateachoutlet distal vessels (the arterioles and capillaries), down-
total
wereadjustedslightlyasneededtoattainmeanflowsat stream of each outlet (Fig. 3).61,62 Below we describe
theoutletsandSBPandDBPattheinletofthemodel the R , C, and R specification for a single patient
p d
to within 5% of target values. Lastly, a final isotropic under rest and exercise conditions; this process was
mesh with an edge size of 0.5 mm was used for each repeated in a patient-specific manner for all subjects.
patient. This resulted in an average mesh size of
9,048,902 linear tetrahedral elements (SD: 2,347,967;
range: 5,731,780–11,951,555). A pulsatile simulation R , C, R Specification for Rest
p d
was then run for eight cardiac cycles using a time step
SpecificationofR ,C,R isrequiredforeachoutlet
sizeof1/1000thofthecardiaccyclewithfivenonlinear p d
of the numerical domain in order to fully define the
iterationspertimestep.Thusthemeantimestepsizewas
Coupled-Multidomain analytical representation of the
0.000974 s (SD: 0.000176; range: 0.000769–0.001333)
downstream vasculature. The values for these param-
forrestand0.000649 s(SD:0.000117;range:0.000513–
eters corresponding to resting conditions were chosen
0.000889)forexercise.Whileall1000timepointsineach
using a combination of patient-specific data and liter-
cardiaccycleweresolvedfor,50timepointswerewritten
ature values as outlined in the following steps:
outpercardiaccycleduetodatastorageconstraints.At
thetimeofsimulation,theresidualforthelastiteration 1. The brachial SBP and DBP acquired immedi-
of each time step was recorded. For each patient, the ately after the exam were assumed to be equal
samemeshwasusedforbothrestandexercise.Thefirst to the intraaortic SBP and DBP at the level of
threecardiaccycleswereignored,asthepressurewaves the diaphragm, which corresponds to the inlet
had not yet achieved periodicity. The fourth through of our model.22
eighth cardiac cycles were used to calculate pressures, 2. Using a three-element Windkessel analog for
velocities, MWSS, OSI, and TKE. For both rest and the system (3D model plus downstream ves-
exercise, outlet velocity profiles were constrained to a sels), we determined the total resistance and
parabolic profile using the augmented lagrangian totalcapacitanceofthesystem,R and
total-system
method developed by Kim et al.28 These constraints C ,byiterativelysolvingforapressure
total-system
were used on an as-needed basis to prevent solution waveform as a function of the input SC flow
divergence, a common numerical problem when an waveform and the parameters of the Wind-
outlet area experiences reverse or complex flow struc- kessel model. R and C were
total-system total-system
tures leaving the three-dimensional computational do- adjusted until the measured SBP and DBP
main. It should be noted that this technique does not were achieved. In this process, we assumed an
affect the amount of flow leaving an outlet, but only R /R ratio of 5.6%.30 The
p total-system total-system
constrainstheshapeofthevelocityprofileattheoutlet. mean blood pressure was calculated as
Use of constraints on branch vessels of the abdominal R times mean SC flow. This is similar
total-system
aorta has been shown to affect the flow, pressure, tothemethodusedbyStergiopolusetal.witha
MWSS,andOSIintheabdominalaortaverylittle.28 two-element Windkessel model.52
Boundary conditions for each patient were derived 3. The target mean flows to each outlet were
from literature and from the measured patient-specific calculated as follows: Target mean upper
SCandIRflowsandthemeasuredbrachialarterySBP branch vessel flows (UBVF) were determined
and DBP. At the wall, a no-slip condition was pre- by subtracting the mean measured IR flow
scribed, as indicated in Eq. (3). The inlet velocity from the mean measured SC flow and then
profile for each patient was defined as follows: distributingthisUBVFbasedontheliterature:
1294 LES etal.
4. Thetotalresistance(R = R + R )ateach
total p d
outlet was calculated as mean blood pressure
dividedbytargetmeanflowateachoutlet.For
the purposes of estimating the initial R va-
total
lue, we assumed that the mean pressure at the
inlet of the model was equal to the mean
pressure at the outlets.
5. Except for the renal arteries, the proximal
resistance, R , for each outlet was arbitrarily
p
chosen to be 5.6% of R .30 For the renal
total
arteries,R wasarbitrarilychosentobe28%of
p
R to account for the unique features of re-
total
nal flow: an absence of backflow throughout
the cardiac cycle and relatively high and con-
stant diastolic flow.3,21 These flow features re-
sult from the relatively low resistance of the
distal vascular beds of the kidneys,41 a feature
captured in the comparatively low R value
d
assigned to the renal arteries.
6. Lastly, capacitance, C, for each outlet was
calculated by distributing the C pro-
total-system
portional to flow at each outlet.
R , C, R and SC Inflow Waveform Specification
p d
for Exercise
Resting waveforms, target SBP and DBP, target
mean flows, and the R , C, R values were altered in
p d
accordance with the literature to simulate moderate
FIGURE3. Apatient-specificSCflowwaveformmeasuredby lower limb exercise.
phase-contrastMRI(restingconditions)ismappedtothein-
flow face using a Womersley velocity profile. A patient-spe- 1. Each patient was categorized as hypertensive or
cific3-elementWindkesselmodelwithaproximalresistance
normotensive,andatargetexerciseSBPandDBP
(R ), capacitance (C), distal resistance (R ) boundary condi-
p d
tionisusedtorepresenttheresistanceandcomplianceofthe werecalculatedfora50%increaseinheartrateby
vasculaturedownstreamofeachoutlet.Theoutletlabelsare linearly interpolating the hypertensive and nor-
as follows: S5splenic artery, LR5left renal artery, LEI5
motensive data from Montain et al.34
left external iliac artery, LII5left internal iliac artery, RII5
right internal iliac artery, REI5right external iliac artery, 2. Target mean flow to each outlet was calculated as
SMA5superior mesenteric artery, ARR5accessory right follows: The target mean resting IR flow was
renal artery (Note: accessory renals are not found in the
increased 5.44-fold, in accordance with in vivo
majorityofpatients),RR5rightrenalartery,H5hepaticar-
tery.Pleasenotethattypically(asinthispatient),thesplenic measurements by Cheng et al.10 and target mean
andhepaticarteriesbranchfromacommonceliactrunk. exercise UBVF was calculated by decreasing the
total flow to each upper branch vessel 21% from
33.0%totheceliac(whichwasthendistributed rest. This is in agreement with Cheng et al.10 and
50% to the hepatic and 50% to the splenic the assumption that flow to the renal, celiac, and
arteries),and22.3%totheSMA,theleftrenal mesenteric system decreases with moderate lower
artery, and the right renal artery.35 In the limb exercise asa resultof the vasoconstriction of
presence of accessory renal arteries, renal flow these vascular beds during exercise. Lastly, the
was divided proportional to outlet area. The mean IR flow was divided equally to the two
mean IR flow was divided equally to the two common iliac arteries; 70% of this flow was di-
common iliac arteries; 70% of this flow was verted to the external iliac arteries and 30% was
diverted to the external iliac arteries and 30% assigned to the internal iliac arteries, thus pro-
was assigned to the internal iliac arteries, thus viding the target mean flow to each iliac artery.
providing the target mean flow to each iliac 3. The target mean exercise IR flow and the target
artery.51 mean exercise UBVF were summed to generate a
Hemodynamics in AAAsDuringRest andExercise 1295
targetmeanexerciseSCflow.Becauseeachpatient
10. Similar to the rest case, the total capacitance,
hadaslightlydifferentmeasuredmeanrestingSC/
C , was distributed proportional to flow
IR flow ratio, the calculated target mean exercise total-system
at each outlet.
SCflowincreasedanaverageof2.57-fold(SD:0.41;
range:1.90–3.07)fromresttoexercise,similartothe
2.61-foldSCincreasefoundinChengetal.10
Dimensionless Variables Governing Flow
4. The diastolic portion of the resting SC flow
waveform was truncated such that heart rate The mean and systolic Reynolds numbers and the
increased 50%,10 and the remaining flow wave- Womersley numbers were calculated from the mea-
form was shifted up to match the target mean sured SC and IR flows as follows:
exercise SC flow. Thus, the resulting exercise SC
q4Q
waveform lacked diastolic backflow and the sys- Reynolds ¼ mean ð13Þ
mean pD l
tolic portion of the waveform represented a mean
greater portion of the cardiac cycle, a result simi-
q4Q
lar to experimental measurements.10 Reynolds ¼ systole ð14Þ
systole pD l
5. This extrapolated patient-specific exercise SC mean
waveform was mapped to the inflow face of the
model using a Womersley velocity profile. (cid:12)xq(cid:13)1=2
WomersleyNumber¼R ð15Þ
6. Similartorest,athree-elementWindkesselanalog l
for the system under exercise conditions (3D
where Q is the mean flow over the cardiac cycle,
model plus downstream vessels) was used to mean
Q isthesystolic(peak)flowoverthecardiaccycle,
determine the total resistance and total capaci- systole
D isthemeandiameteroftheaortaattheSCorIR
tanceofthesystem,R andC ,by mean
total-system total-system
location, l is the dynamic viscosity of 0.04 P, q is the
iteratively solving for a pressure waveform as a
density (1.06 g/cm3), and x is the angular frequency
function of the input extrapolated SC flow wave-
(equal to 1/cardiac cycle length).
formandtheparametersoftheWindkesselmodel.
R andC wereadjusteduntilthe
total-system total-system
measured SBP and DBP were achieved. In this
Normal Stresses
process, we assumed an R /R
p total-system total-system
ratio of 5.6%. The mean blood pressure was cal- The distribution of normal stresses (including both
culated as Rtotal-system times target mean SC exer- the static and dynamic viscous components) acting on
cise flow. the aneurysm lumen at peak flow was plotted for all
7. Using the target exercise pressures and mean eight patients for rest and exercise. We described the
flows, the resting R , C, R parameters were temporalandspatialvariationofthetotalnormalstress
p d
altered to simulate exercise. First, R was cal- acting on the aneurysm luminal wall: First, we calcu-
total
culated for each outlet as the mean pressure latedthespatialrangeofthenormalstressatpeakflow
divided by the target mean exercise flow for each for both rest and exercise. Then, we computed the
branch vessel. temporal range of the normal stresses by subtracting
8. Rp was kept the same as rest; as the arteries just themaximumandminimumaveragenormalstressesin
distal to the outlets of the model were assumed to the aneurysm over the cardiac cycle for both rest and
neither dilate nor constrict significantly during exercise. Last, we calculated the ratio between the
exercise. spatial and temporal ranges of the total normal stress
9. R was calculated as R minus R . Thus, R foralleightpatientsduringrestandexercise.Thetime
d total p d
reflected the physiologic states of the distal vas- points corresponding to peak flow and maximum nor-
cular beds during lower limb exercise: Rd mal stresses in the aneurysms were not necessarily the
increased in the upper branch vessels during same, aspressure generallylags flow.
exercise, reflecting the vasoconstriction of the
arteriolesandcapillaries inceliac, mesenteric,and
Shear Stress
renal systems during exercise, while R dramati-
d
callydecreasedintheiliac arteries,21reflectingthe MWSS and OSI were averaged over the fourth
vasodilation of arterioles and recruitment of the through eighth cardiac cycles and then quantified in
capillaries during exercise. This unchanging R 1-cm strips at four aortic locations: 1 cm above the
p
and dramatically reduced R during exercise is celiac artery (supraceliac, SC), 1 cm below the most
d
supported by experimental data reported by distal renal artery (infrarenal, IR), at the largest
Laskey et al.30 diameter of the aneurysm (Mid-An), and 1 cm above
1296 LES etal.
the aortic bifurcation (suprabifurcation, SB). When This ensemble average is then extended over n
averaging over the fourth through eighth cardiac cardiac cycles to cover the same temporal domain as
cycles, the MWSS, *smean; was defined as: the original velocity field u*(cid:1)x*;t(cid:2). This extension is
*smean ¼(cid:7)(cid:7)(cid:7)(cid:7)51TZ 8T*tsdt(cid:7)(cid:7)(cid:7)(cid:7): ð16Þ accomplished by repeating Du*E(cid:1)x*;s(cid:2) n times to
3T generatetheperiodicaverageU*(cid:1)x*;t(cid:2).Mathematically,
Here, T is the period of one cardiac cycle, *t is the
s
wallshearvectordefinedasthein-planecomponentof
(cid:1) (cid:2) D E(cid:1) (cid:2)
the surface traction vector *t (i.e., *t ¼*t (cid:2)(cid:1)*t (cid:1)*n(cid:2)*n). U* x*;t ¼ *u x*; tmoduloT ; ð19Þ
s
Thesurfacetraction vector*t iscalculatedbysolvinga wheretmoduloTgivestheremainderwhentisdivided
newvariationalproblemwherethesolutionforvelocity byT.Now,thefluctuationvelocityfieldcanbedefined
and pressure is known, and the goal is to evaluate the by subtracting the periodic average from the original
tractionattheboundariesofthedomain.Thisapproach velocity field:
is known as variationally consistent calculation of
boundary fluxes.24 When averaging over the fourth u*~(cid:1)x*;t(cid:2)¼u*(cid:1)x*;t(cid:2)(cid:2)U*(cid:1)x*;t(cid:2): ð20Þ
througheighthcardiaccycles,OSI,whichmeasuresthe
unidirectionalityofshearstress,wascalculatedas: Finally,themeanfluctuatingkineticenergyperunit
0 (cid:7) (cid:7)1 mass, often referred to as the turbulent kinetic energy,
(cid:7)1 R8T*t dt(cid:7)
1 (cid:7)5T 3T s (cid:7) is defined as:
OSI¼2B@1(cid:2)51TR38TT(cid:7)(cid:7)(cid:7)*ts(cid:7)(cid:7)(cid:7)dtCA ð17Þ k(cid:1)x*;s(cid:2)¼1qh(cid:14)u~2(cid:15)(cid:1)x*;s(cid:2)þ(cid:14)u~2(cid:15)(cid:1)x*;s(cid:2)þ(cid:14)u~2(cid:15)(cid:1)x*;s(cid:2)i;
2 1 2 3
where T isagain theperiod ofone cardiac cycle andt
s 8s2½0;TÞ; ð21Þ
is the in-plane component of the surface traction vec-
tor. OSI values range from 0 to 0.5, where 0 indicates wherethedensity,qis1.06 g/cm3andu~ denotestheith
i
unidirectional shear stress, and 0.5 indicates shear component of *u~; and
stress with a time average of zero.20,56
Turbulent Kinetic Energy (cid:14)u~2i(cid:15)ðx*;sÞ¼1nXn(cid:2)1u~2iðx*;sþkTÞ; 8s2½0;TÞ: ð22Þ
k¼0
InordertostudyturbulenceintheAAAmodels,we
Similarly,theensemble-averagedkineticenergy,KE,
investigated the cycle-to-cycle variation of the velocity
isdefinedas:
fieldbycomputingtheturbulentkineticenergy(TKE).
A brief summary of the derivation of TKE is shown
below. k(cid:1)x*;s(cid:2)¼1qh(cid:14)u2(cid:15)(cid:1)x*;s(cid:2)þ(cid:14)u2(cid:15)(cid:1)x*;s(cid:2)þ(cid:14)u2(cid:15)(cid:1)x*;s(cid:2)i;
The TKE computation is centered upon a decom- 2 1 2 3
(cid:1) (cid:2) 8s2½0;TÞ ð23Þ
position of the velocity field u* x*;t into an averaged
part U*(cid:1)x*;t(cid:2) plus a fluctuating part u*~(cid:1)x*;t(cid:2). First, we wherethedensityqis1.06 g/cm3anduidenotestheith
component of u, and
(cid:1) (cid:2)
assume that the velocity field, *u x*;t , is nominally
periodic with respect to the cardiac cycle, where x* (cid:14)u2(cid:15)ðx*;sÞ¼1Xn(cid:2)1u2ðx*;sþkTÞ; 8s2½0;TÞ ð24Þ
i n i
denotes space and t denotes time. Next, an ensemble k¼0
average is computed by averaging out the cycle-to-
ToquantifytheTKEandKEwithintheabdominal
cycle changes in n cardiac cycles worth of data and
aorta and aneurysm, the abdominal aorta was first
mapping them to a single cardiac cycle. Mathemati-
(cid:1) (cid:2) isolatedbytrimmingthemodelsuperiorlyjustdistalto
cally,givenncardiaccyclesofthevelocityfieldu* x*;t the lowest renal artery and inferiorly just proximal to
where t2½0;nTÞ; we define the ensemble average theaorticbifurcation.ThemeanTKEandKEatpeak
D E(cid:1) (cid:2) systole, mid-deceleration, and mid-diastole were
*u x*;s as:
quantified over cardiac cycles 4 though 8 on the final
meshes for all eight patients during both rest and
D E(cid:1) (cid:2) 1Xn(cid:2)1 (cid:1) (cid:2)
*u x*;s ¼ u* x*;sþkT ; 8s2½0;TÞ: ð18Þ exercise.
n
k¼0
Hemodynamics in AAAsDuringRest andExercise 1297
Mesh and Temporal Independence Study RESULTS
In order to ensure the mesh-independence of our Morphology of AAAs
solutions, three isotropic meshes consisting of
The MRAs showed that AAA luminal morphology
2,175,303 (edge size: 0.82 mm), 8,343,663 (edge size of
varied widely from patient to patient, with shapes
0.5 mm), and 31,831,619 linear tetrahedral elements
varying from fusiform (patients 1, 4, and 5), single-
(edge size 0.3 mm) were generated for patient 1
lobed (patients 2 and 7), and bi-lobed (patients 3, 6,
(Fig. 4).Hereafter,thesemesheswillbeabbreviatedas
and 8) (Fig. 2).
the 2.2, 8.3, and 31.8-million element meshes. These
mesheswerechosensuchthatratioofelementfromthe
smalltointermediateandintermediatetolargemeshes Measured Flows and Pressures
would be approximately the same. Simulations were
The mean flows measured from the PC-MRI scans
performed for eight cardiac cycles using the same
were 3.51 L/min (SD: 0.50; range: 2.51–4.32) for SC
patient-specific exercise boundary conditions. The
and 1.31 L/min (SD: 0.34; range: 0.84–1.75) for IR.
resultant instantaneous velocity field and MWSS and
The mean measured SBP was 155.8 mmHg (SD: 12.2;
OSIvalueswerecompared.Additionally,toensurethe
range: 136–170) and the mean measured DBP was
temporal independence of our solutions, the patient 1
90.9 mmHg (SD: 6.8, range: 80–100), indicating the
exercise simulation on the final mesh (8.3-million ele-
presence of hypertension.
ments) was extended for an additional five cardiac
cycles.TheresultantMWSS,OSI,andTKEcalculated
from cardiac cycles four through eight were compared Simulation
to those calculated from cardiac cycles nine through
Forthefinalfiniteelementmeshes,themeanabsolute
thirteen.
residual for cycles four through eight was 2.9e-4 (SD:
1.7e-4, range: 1.0e-5–6.0-e4) under resting conditions
and 1.8e-3 (SD: 1.2e-3; range: 1.5e-4–3.7e-3) under
Comparison of Measured and Simulated Infrarenal
exercise conditions. For the mesh independence study,
Flow Waveforms
the average absolute residual was 4.6e-3 for the 2.2-
The measured and simulated IR flow waveforms millionelementmesh,2.3e-3forthe8.3-millionelement
were compared. The simulated IR waveforms resulted mesh,and7.1e-4forthe31.8millionelementmesh.For
purely from the inlet SC waveforms and the outlet the temporal independence study of patient 1 under
boundary conditions; the shape of the IR waveforms exercise conditions, the mean residual was 2.3e-4 for
was never prescribed or specified in the simulation bothcardiaccyclesfourthrougheightandcardiaccycles
process, and thus the IR waveforms provided a nine through thirteen. For each of the 16 simulations
good measure of how well simulated flow conditions (eight patients, rest and exercise conditions), the
matched the measured flow conditions for each boundaryconditiontuningprocesstookapproximately
patient. 1week,andthefinalsimulationtookapproximately3–7
Mesh Cross Section At The AAA
2.2-million elements 8.3-million elements 31.8-million elements
FIGURE4. Cross sections of the 2.2-, 8.3-, and 31.8-million finite-element meshes are shown at the mid-aneurysm level for
patient1.