M a r k V a n D o o rm a a l M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada
Y u -Q in g Z h o u M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada e-m ail: yqzhou@ m ouseim aging.ca
X ia o li Z h a n g M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada e-m ail: xzhang@ m ouseimaging.ca
D a v id A. S te in m a n Mechanic and Industrial Engineering, University of Toronto, 5 King’s College Road, Toronto, ON M 5 S 3G 8, Canada e-m ail: steinman@ mie.utoronto.ca
R . M a r k H e n k e lm a n 1 M em , ASME M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada e-m ail: m henkel@ mouseim aging.ca
Inputs for Subject-Specific Computational Fluid Dynamics Simulation of Blood Flow in the Mouse Aorta Mouse models are an important way for exploring relationships between blood hemody namics and eventual plaque formation. We have developed a mouse model of aortic regurgitation (AR) that produces large changes in plaque burden with charges in hemo dynamics [Zhou et al., 2010, "Aortic Regurgitation Dramatically Alters the Distribution o f Atherosclerotic Lesions and Enhances Atherogenesis in Mice," Arterioscler. Thromb. Vase. Biol., 30(6), pp. 1181-1188], In this paper, we explore the amount of detail needed for realistic computational fluid dynamics (CFD) calculations in this experimental model. The CFD calculations use inputs based on experimental measurements from ultrasound (US), micro computed tomography (CT), and both anatomical magnetic resonance imaging (MRI) and phase contrast MRI (PC-MRI). The adequacy of five different levels of model complexity (a) subject-specific CT data from a single mouse; (b) subject-specific CT centerlines with radii from US; (c) same as (b) but with MRI derived centerlines; (d) average CT centerlines and averaged vessel radius and branching vessels; and (e) same as (d) but with averaged MRI centerlines) is evaluated by demonstrating their impact on relative residence time (RRT) outputs. The paper concludes by demonstrating the neces sity o f subject-specific geometry and recommends for inputs the use of CT or anatomical MRI for establishing the aortic centerlines, M-mode US for scaling the aortic diameters, and a combination of PC-MRI and Doppler US for estimating the spatial and temporal characteristics o f the input wave forms. [DOI: 10.1115/1.4028104] Keywords: CFD, mouse aorta, subject-specific, ultrasound, microCT, MRI, RRT
1 Introduction A novel method for inducing atherosclerotic plaques in other wise lesion free areas of the descending thoracic and abdominal aorta of mice has been developed [1]. This method creates AR by surgically damaging the aortic valve of the mouse. Preliminary studies of this model of atherosclerosis in mice have shown a clear link between the hemodynamic environment, which is modified by the induction of AR, and atherosclerotic plaque patterns. This model is especially attractive because the degree of AR can be controlled and measured and is the only independent variable. Diet and genetics are kept constant. Previous work using Doppler US and oil red O lesion staining has shown that despite oscillatory flow throughout the descending aorta, plaques tend to form in spe cific, consistent areas of the vessel wall [1]. Subsequent work by Hoi et al. showed that the areas in the descending aorta most prone to plaques are also areas subject to low time averaged wall shear stress (TAWSS), high oscillatory wall shear stress index (OSI), and high relative residence time (RRT, a measure propor tional to the amount of time that particles spend close to a particu lar point on the wall) [2]. Since RRT correlated best with pattern of plaque, it is the metric used throughout this paper. The areas prone to plaque formation in mice with AR were also associated with regions of curvature in the aorta [2], 'Corresponding author. Manuscript received November 15, 2013; final manuscript received July 18, 2014; accepted manuscript posted July 30, 2014; published online August 12, 2014. Assoc. Editor: Alison Marsden.
Journal of Biomechanical Engineering
The previous work [2] studying the hemodynamics of mice with AR has not taken into account mouse-specific flow (AR severity) or mouse-specific geometric boundary conditions (curva ture of the aorta) that are likely to impact both local wall shear stress, RRT and atherosclerotic plaque patterns. Additionally, the focus of the previous modeling work by Hoi et al. [2] has been on the descending aorta and not on the aortic arch, which is prone to atherosclerotic lesions on the inner curvature of the arch of the mouse regardless of the presence of AR. In order to determine subject-specific hemodynamic patterns in the aorta of mice, with and without AR, using CFD, two types of initial conditions are required: (1) Subject-specific aortic geometry, specifically from the ascending aorta distal to the aortic sinus and then proceed ing to the descending thoracic aorta down to the abdominal aorta at the level of the celiac artery, including the three major branches of file arch. (2) Subject-specific flow velocity data from at least the ascend ing aorta, three major branches of the aortic arch, and the descending aorta just proximal to the celiac artery. Addi tional flow velocity data at intermediate points in the descending aorta can also be used to quantify flow rate attenuation and flow out of the intercostal arteries. Even though the eventual plaque distributions can be read out using optical projection tomography [3,4] or microCT of fixed samples, to allow development of plaque over time, subjectspecific initial boundary conditions must be obtained in live mice.
Copyright ©2014 by ASME
OCTOBER 2014, Vol. 136 / 101008-1
For example, a combination of ex vivo CT scanning of a corrosion cast which is a terminal procedure and in vivo PC-MRI (e.g., Ref. [5]) would not allow for plaque development and quantifica tion and is hence not a viable solution for obtaining the initial conditions. Although ex vivo microCT has been used extensively to deter mine the geometry of the mouse aorta [5-7], it cannot be consid ered the “gold standard.” Ideally, the gold standard method for acquiring the geometry of the mouse aorta would accurately cap ture the in vivo curvature, vessel radius, branch angles, and vessel cross-sectional shape. With this in mind, in vivo microCT would seem to fit the criteria; however, problems do exist with this method. Previous work by Vandeghinste has shown differences in vessel diameter and branch angles in the mouse aorta between ex vivo microCT (corrosion casting) and in vivo microCT [8]. MR angiography could also be used to establish initial bound ary conditions for detailed CFD calculations [9,10]. However, the spatial resolution obtainable in the mouse is not sufficient to satis factorily delineate the geometries of all vessel walls. In this paper, we describe a method that uses in vivo MRI in combination with M-mode US to determine the geometry of the ascending aorta, aortic arch, and descending aorta. This method is also capable of providing in vivo blood velocity measurements that can be used for subject-specific CFD studies, and does not preclude the subsequent use of ex vivo methods to determine plaque patterns at the end of the experiment. This method was compared with current ex vivo CT based methods of acquiring mouse aortic geometry data in order to understand the method and justify its use. We compared the RRT patterns using the described method to the ex vivo CT based method in both the aortic arch and the descending aorta in a group of mice. The differences in RRT patterns due to imaging modality are compared to the differ ences in RRT patterns due to AR severity in order to help quantify the impact of differing modalities in the context of a study with varying degrees of AR. The relative merits of each method in the aortic arch and the descending aorta of the mouse are discussed. Additionally, we investigated whether a geometry created as an average of several mice can be adequately used to substitute for subject-specific geometries in CFD studies. The paper concludes by demonstrating the necessity of subject-specific geometries for CFD calculations in the aorta and shows the adequacy of a combi nation of PC-MRI and Doppler US for estimating input wave forms. 2
Methods
2.1 Data Acquisition and Image Analysis. The aortas of 20 Ldlr—/ — mice (14 weeks of age) were imaged using several imaging modalities. All mice were imaged with all of the imag ing modalities. Of these, the data from 14 mice were determined subjectively to be of sufficient quality to be included in the study. The excluded data included severe overlapping motion artifacts, two MRI technical failures, one mouse death and one incomplete Microfil® (Flow Tech, Carver, MA) perfusion. First, M-mode US recordings (VisualSonics 2100 with a lateral resolution of 115/un and an axial resolution of 50/un) were made at the beginning of ascending aorta, proximal section of the descending thoracic aorta, and abdominal aorta in order to determine the time varying aortic diameter [1], Second, MR an giography of the aorta with spatial resolution of 120 x 120 yum2 and slice thickness of 1 mm was conducted using a 7 T scanner using a 2D spin echo sequence with electrocardiogram (ECG) and respiratory gating. Localization images were acquired and 22 slices along the aorta were imaged perpendicular to the vessel centerline as shown in Fig. 1: six slices on the aortic arch (top panel) and 16 slices on the descending aorta (bottom panel). Finally, a Microfil® perfusion of the mouse vasculature was cre ated followed by partial dissection of the heart, aorta, and spine. This sample was then scanned using microCT [11] to obtain detailed anatomy. 101008-2 / Vol. 136, OCTOBER 2014
The microCT scanner was a Bruker 1172 (formerly SkyScan) specimen scanner operated at a spatial resolution of 28 /im. Although microCT angiography would have provided more physi ologically realistic geometry, it would not have provided suffi ciently high definition of the luminal wall. The diameters of the ascending aorta, proximal descending thoracic aorta, and proximal abdominal aorta were determined from the M-mode images by using a maximum gradient method to determine the wall locations, for three levels corresponding to MRI slices (Fig. 1): (1) The proximal ascending aorta: The M-mode cursor line was placed at the level immediately distal to the aortic sinus. The aortic sinus is able to be identified in both US imaging and live MRI. (2) The initial segment of the descending thoracic aorta: The M-mode cursor line was placed at the level where the curved aortic arch became straight at the descending tho racic aorta. It is approximately the level of the junction between the atria and ventricles which was used as a land mark in both US and MRI. (3) The proximal abdominal aorta: The M-mode cursor line was placed at the level prior to the celiac artery (the first major branch from the abdominal aorta). The diameter of the abdominal aorta between the diaphragm and the celiac artery is relatively consistent; therefore, slight shift of slice level should not cause significant variation of diameter measurement. The results of the MR angiography were segmented manually and the vessel centerline (in three-dimensional (3D) space) was created by linking the geometric centers of the segmented vessel slices. A smooth spline representation of the centerline was cre ated using the Numpy spline interpolation routines (numpy.splrep, numpy.scipy.org). The microCT data was segmented using the level set method of the vascular modeling toolkit (VMTK, vmtk.org). We then used VMTK to determine vessel centerline and radius based on the perimeter of the vessel cross section (as opposed to maximum ra dius of the inscribed sphere). A perimeter based method was used because in some cases insufficient perfusion pressures resulted in slightly noncircular vessel cross sections. “Average” mouse centerlines were determined for both MR and microCT geometries as follows. First, the subject-specific centerlines were aligned using rigid body rotation and translation of the centerlines by minimizing the root mean square displacement of the minimum distance between centerlines. To accomplish this, a series of 1 mm spaced points along the centerline were chosen and for each point, the largest sphere that did not intersect the other centerline was taken as a measure of the distance between the cen terlines. The root mean square of the sum of these distances was taken as a measure of the displacement between the centerlines. This displacement was minimized over the six degrees of freedom of the rigid body rotation. The alignment achieved was independ ent of which centerline was chosen as the target. Next, the 3D locations of individual points on the average centerline were cal culated. For each abscissa value (distance from the aortic root along the centerline), the average 3D location of this point was calculated from the 14 mice. The average radius along the centerline was also determined for the CT data. Average branch centerlines were determined in a similar man ner to the aorta centerlines as described above. Additionally, using VMTK, branch vessel locations (distance along centerline from aortic root), in-plane vessel angle, and out-of-plane vessel angles were determined for each branch for each mouse. For each branch, the vessel location, in-plane angle, and out-of-plane angle were averaged across all 14 mice. This data subsequently was used to create both idealized subject-specific geometries and a mouse-averaged geometry. In-plane and out-of-plane angles were defined as the angles between the daughter vessel (the branch) and the parent vessel (the aorta before the branch) as seen in Fig. 2. Transactions of the ASME
Fig. 1 Graphic prescriptions from MR angiography determ ine the aortic centerline and for in vivo US M-m ode recordings to m easure aortic diam eters. Two localization images are required to localize the aortic arch and descending aorta images separately, with six slices (AAr_1-AA rJ>) acquired on aortic arch (upper left) and 16 slices (DAo_15-DAo_16) on descend ing aorta (low er left). On the right, B-mode US images along the centerline of the aorta at three cross-sectional levels are shown. US M-m ode recordings were made at the proximal ascending aorta im m ediately distal to the aortic sinus (AAr_1) (upper right), at the initial segm ent of the descending thoracic aorta where the curved aortic arch becam e straight (Dao„1) (m iddle right), and at the proximal abdom inal aorta prior to the celiac artery (Dao_15) (lower right). LV is the left ventricle and Li is the liver.
2.2 Creating Vessel Geometries. Five different geometries were created incorporating different degrees of individual anat omy, and with the aortic centerline obtained from either CT or MRI. The third geometry (c) is the proposed method and the other methods (a, b, d, and e) serve as comparisons. We will refer to these as geometry construction methods throughout the text. The five geometry construction methods are: (a) Subject-specific CT data from a single mouse, with no radius corrections. This method ignores all the MRI and US imaging data. (b) Subject-specific CT derived centerline from a single mouse, with a circular cross section and radius corrected to subjectspecific M-mode US measurements. (c) Subject-specific MRI derived centerlines from a single mouse, with circular cross section and radius corrected to Journal of Biomechanical Engineering
subject-specific M-mode US measurements, and mouseaveraged aortic arch branch centerlines, locations, and branch angles. At the conclusion of the paper, this is the preferred geometry. (d) Mouse-averaged CT derived centerline with mouseaveraged vessel radius, aortic arch branch centerlines, loca tions, and branch angles. This is a single representation of the average behavior of all of the aortas with no individual subject specific data remaining. (e) Mouse-averaged MRI derived centerline with mouseaveraged vessel radius, aortic arch branch centerlines, loca tions, and branch angles. The individual single mouse chosen for geometries a-c was selected based on high CT data quality as well as being represen tative of a mouse with aortic arch branch angles and locations OCTOBER 2014, Vol. 136 / 101008-3
Fig. 2 In-plane (left) and out-of-plane (right) angles shown for the first branch in an example mouse aorta. The right side view is normal to the left side view.
being approximately one standard deviation from the mean i.e., not representative of an average mouse. This single mouse allowed us to determine eventually that reliable CFD results could not be obtained if geometry from an average over all mice was used. The microCT and MRI geometries (geometries b and c) were constructed from the radius and centerline data using VMTK. Here, the mouse-specific CT radius was scaled so the radius at the proximal descending thoracic and abdominal aorta matched the cycle averaged radius measured with US at both locations. Radius data for the subject-specific MRI data was determined by using the average CT radius data along the centerline, also corrected with cycle averaged subject-specific US data at the proximal descending thoracic aorta and abdominal aorta. 2.2.1 Additional Mouse-Specific Geometries. After investi gating the five methods of constructing vessel geometry (a-e), and concluding that only b and c were sufficiently adequate methods, we chose five different mice and generated geometries by methods b and c for each of them. These geometries were constructed using methods b and c as described above but without including the arch and the branches in the arch. The arch was included only as a flow extension to ensure a realistic velocity profile at the proximal portion of the descending thoracic aorta (DTAo, as in Ref. [2]). 2.3 Flow Boundary Conditions. Traction free boundary con ditions were imposed at the descending aorta outlet while Womersley profiles were imposed at the three aortic arch branches according to the temporal flow rate waveform. Rigid, no slip walls, and a constant kinematic blood viscosity of 3.5 x 106 m2/s and a density of 1.06 g/cm3 were assumed. This assumption of a Newtonian fluid is appropriate given the high shear rates present in the mouse aorta as per Feintuch et al. [10]. We chose to impose AR flow conditions [1] in this study because it provides a more complicated velocity field and should result in a higher probability of identifying differences between the methods tested [12]. However, flow splits in the three major branches of the aortic arch in mice with AR are not yet known. In order to estimate the flow rates and waveforms at the three branches, Doppler US was completed in the innominate artery, the left common carotid artery, and the left subclavian artery in five mice with moderate AR. Using intensity weighted mean velocity tracing, the relative flow split to each branch was deter mined. In all five mice, in all three branches, despite retrograde flow in diastole in the aortic arch, there was very little retrograde flow in the branch arteries. Hence, in our simulations an idealized temporal flow waveform was used in each branch with the systolic portion having the same shape as the systolic portion of the wave form in the descending aorta scaled to ensure proper flow splits. The diastolic portion of the waveform was simply set to zero to reflect the lack of flow in diastole. The temporal flow waveform at the aortic root (inlet) was set to be the sum of the three branch
101008-4
/ Vol. 136, O C T O B E R 2014
waveforms and the waveform used by Hoi et al. in the descending aorta [2], thus ensuring that the flow waveform in the descending aorta was equal to that used by Hoi et al. At the aortic root we imposed a flat velocity inlet profile. Although this is not the ideal boundary condition for determining wall shear stress in the arch [5], it is a reasonable simplification and our focus in this study was on comparing methods, not determining specific wall shear stress patterns in the aortic arch of mice with AR. In addition, two cases of reduced AR severity were simulated in the mouse-averaged CT derived geometry. The purpose of these simulations was to determine the relative effects of AR severity compared to varying the geometry construction method. The temporal flow waveform imposed at the aortic root was the same except the diastolic portion of the waveform (i.e., the retro grade portion) was reduced such that the diastolic flow was a frac tion (70% and 85%) of the original value. The flow waveforms to the three branches remained unchanged. 2.4 CFD Simulations. ICEM-CFD was used to generate quadratic tetrahedral finite element meshes of approximately 1.000. 000 elements for the five different geometries. These 1.000. 000 elements correspond to roughly 8 x 106 linear tetrahedra with over 1 x 106 nodes. This is about 10 x the number of nodes reported being used by Trachet at al. [13] who report satis factory convergence. The unsteady Navier-Stokes equations were solved using the AR flow waveforms detailed above via a previ ously validated in-house solver [2,14]. In order to ensure that initial transients were not present in the solution, at least three cardiac cycles were simulated with at least 4800 time steps included per cardiac cycle to ensure stability. The choice of the number of nodes and the number of time steps is probably excessively large and more than is necessary. If a large number of CFD calculations were being contemplated, investiga tions into reducing this computer burden would need to be under taken. We chose to focus on RRT, as this measure correlated best with plaque patterns in a previous study of the AR model [2], Although other measures of wall shear such as SAWSS or TAWSS or OSI could have been used, we have shown previously that OSI or RRT are very similar and all the measures are highly correlated; therefore, we have chosen to show RRT as it is recom mended as a “robust single metric of low and oscillatory shear” [15], Himburg et al. [16] define RRT = T/ | J tw
Y u -Q in g Z h o u M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada e-m ail: yqzhou@ m ouseim aging.ca
X ia o li Z h a n g M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada e-m ail: xzhang@ m ouseimaging.ca
D a v id A. S te in m a n Mechanic and Industrial Engineering, University of Toronto, 5 King’s College Road, Toronto, ON M 5 S 3G 8, Canada e-m ail: steinman@ mie.utoronto.ca
R . M a r k H e n k e lm a n 1 M em , ASME M ouse Imaging Centre, Hospital for Sick Children, 25 Orde Street, Toronto, ON M 5T 3H 7, Canada e-m ail: m henkel@ mouseim aging.ca
Inputs for Subject-Specific Computational Fluid Dynamics Simulation of Blood Flow in the Mouse Aorta Mouse models are an important way for exploring relationships between blood hemody namics and eventual plaque formation. We have developed a mouse model of aortic regurgitation (AR) that produces large changes in plaque burden with charges in hemo dynamics [Zhou et al., 2010, "Aortic Regurgitation Dramatically Alters the Distribution o f Atherosclerotic Lesions and Enhances Atherogenesis in Mice," Arterioscler. Thromb. Vase. Biol., 30(6), pp. 1181-1188], In this paper, we explore the amount of detail needed for realistic computational fluid dynamics (CFD) calculations in this experimental model. The CFD calculations use inputs based on experimental measurements from ultrasound (US), micro computed tomography (CT), and both anatomical magnetic resonance imaging (MRI) and phase contrast MRI (PC-MRI). The adequacy of five different levels of model complexity (a) subject-specific CT data from a single mouse; (b) subject-specific CT centerlines with radii from US; (c) same as (b) but with MRI derived centerlines; (d) average CT centerlines and averaged vessel radius and branching vessels; and (e) same as (d) but with averaged MRI centerlines) is evaluated by demonstrating their impact on relative residence time (RRT) outputs. The paper concludes by demonstrating the neces sity o f subject-specific geometry and recommends for inputs the use of CT or anatomical MRI for establishing the aortic centerlines, M-mode US for scaling the aortic diameters, and a combination of PC-MRI and Doppler US for estimating the spatial and temporal characteristics o f the input wave forms. [DOI: 10.1115/1.4028104] Keywords: CFD, mouse aorta, subject-specific, ultrasound, microCT, MRI, RRT
1 Introduction A novel method for inducing atherosclerotic plaques in other wise lesion free areas of the descending thoracic and abdominal aorta of mice has been developed [1]. This method creates AR by surgically damaging the aortic valve of the mouse. Preliminary studies of this model of atherosclerosis in mice have shown a clear link between the hemodynamic environment, which is modified by the induction of AR, and atherosclerotic plaque patterns. This model is especially attractive because the degree of AR can be controlled and measured and is the only independent variable. Diet and genetics are kept constant. Previous work using Doppler US and oil red O lesion staining has shown that despite oscillatory flow throughout the descending aorta, plaques tend to form in spe cific, consistent areas of the vessel wall [1]. Subsequent work by Hoi et al. showed that the areas in the descending aorta most prone to plaques are also areas subject to low time averaged wall shear stress (TAWSS), high oscillatory wall shear stress index (OSI), and high relative residence time (RRT, a measure propor tional to the amount of time that particles spend close to a particu lar point on the wall) [2]. Since RRT correlated best with pattern of plaque, it is the metric used throughout this paper. The areas prone to plaque formation in mice with AR were also associated with regions of curvature in the aorta [2], 'Corresponding author. Manuscript received November 15, 2013; final manuscript received July 18, 2014; accepted manuscript posted July 30, 2014; published online August 12, 2014. Assoc. Editor: Alison Marsden.
Journal of Biomechanical Engineering
The previous work [2] studying the hemodynamics of mice with AR has not taken into account mouse-specific flow (AR severity) or mouse-specific geometric boundary conditions (curva ture of the aorta) that are likely to impact both local wall shear stress, RRT and atherosclerotic plaque patterns. Additionally, the focus of the previous modeling work by Hoi et al. [2] has been on the descending aorta and not on the aortic arch, which is prone to atherosclerotic lesions on the inner curvature of the arch of the mouse regardless of the presence of AR. In order to determine subject-specific hemodynamic patterns in the aorta of mice, with and without AR, using CFD, two types of initial conditions are required: (1) Subject-specific aortic geometry, specifically from the ascending aorta distal to the aortic sinus and then proceed ing to the descending thoracic aorta down to the abdominal aorta at the level of the celiac artery, including the three major branches of file arch. (2) Subject-specific flow velocity data from at least the ascend ing aorta, three major branches of the aortic arch, and the descending aorta just proximal to the celiac artery. Addi tional flow velocity data at intermediate points in the descending aorta can also be used to quantify flow rate attenuation and flow out of the intercostal arteries. Even though the eventual plaque distributions can be read out using optical projection tomography [3,4] or microCT of fixed samples, to allow development of plaque over time, subjectspecific initial boundary conditions must be obtained in live mice.
Copyright ©2014 by ASME
OCTOBER 2014, Vol. 136 / 101008-1
For example, a combination of ex vivo CT scanning of a corrosion cast which is a terminal procedure and in vivo PC-MRI (e.g., Ref. [5]) would not allow for plaque development and quantifica tion and is hence not a viable solution for obtaining the initial conditions. Although ex vivo microCT has been used extensively to deter mine the geometry of the mouse aorta [5-7], it cannot be consid ered the “gold standard.” Ideally, the gold standard method for acquiring the geometry of the mouse aorta would accurately cap ture the in vivo curvature, vessel radius, branch angles, and vessel cross-sectional shape. With this in mind, in vivo microCT would seem to fit the criteria; however, problems do exist with this method. Previous work by Vandeghinste has shown differences in vessel diameter and branch angles in the mouse aorta between ex vivo microCT (corrosion casting) and in vivo microCT [8]. MR angiography could also be used to establish initial bound ary conditions for detailed CFD calculations [9,10]. However, the spatial resolution obtainable in the mouse is not sufficient to satis factorily delineate the geometries of all vessel walls. In this paper, we describe a method that uses in vivo MRI in combination with M-mode US to determine the geometry of the ascending aorta, aortic arch, and descending aorta. This method is also capable of providing in vivo blood velocity measurements that can be used for subject-specific CFD studies, and does not preclude the subsequent use of ex vivo methods to determine plaque patterns at the end of the experiment. This method was compared with current ex vivo CT based methods of acquiring mouse aortic geometry data in order to understand the method and justify its use. We compared the RRT patterns using the described method to the ex vivo CT based method in both the aortic arch and the descending aorta in a group of mice. The differences in RRT patterns due to imaging modality are compared to the differ ences in RRT patterns due to AR severity in order to help quantify the impact of differing modalities in the context of a study with varying degrees of AR. The relative merits of each method in the aortic arch and the descending aorta of the mouse are discussed. Additionally, we investigated whether a geometry created as an average of several mice can be adequately used to substitute for subject-specific geometries in CFD studies. The paper concludes by demonstrating the necessity of subject-specific geometries for CFD calculations in the aorta and shows the adequacy of a combi nation of PC-MRI and Doppler US for estimating input wave forms. 2
Methods
2.1 Data Acquisition and Image Analysis. The aortas of 20 Ldlr—/ — mice (14 weeks of age) were imaged using several imaging modalities. All mice were imaged with all of the imag ing modalities. Of these, the data from 14 mice were determined subjectively to be of sufficient quality to be included in the study. The excluded data included severe overlapping motion artifacts, two MRI technical failures, one mouse death and one incomplete Microfil® (Flow Tech, Carver, MA) perfusion. First, M-mode US recordings (VisualSonics 2100 with a lateral resolution of 115/un and an axial resolution of 50/un) were made at the beginning of ascending aorta, proximal section of the descending thoracic aorta, and abdominal aorta in order to determine the time varying aortic diameter [1], Second, MR an giography of the aorta with spatial resolution of 120 x 120 yum2 and slice thickness of 1 mm was conducted using a 7 T scanner using a 2D spin echo sequence with electrocardiogram (ECG) and respiratory gating. Localization images were acquired and 22 slices along the aorta were imaged perpendicular to the vessel centerline as shown in Fig. 1: six slices on the aortic arch (top panel) and 16 slices on the descending aorta (bottom panel). Finally, a Microfil® perfusion of the mouse vasculature was cre ated followed by partial dissection of the heart, aorta, and spine. This sample was then scanned using microCT [11] to obtain detailed anatomy. 101008-2 / Vol. 136, OCTOBER 2014
The microCT scanner was a Bruker 1172 (formerly SkyScan) specimen scanner operated at a spatial resolution of 28 /im. Although microCT angiography would have provided more physi ologically realistic geometry, it would not have provided suffi ciently high definition of the luminal wall. The diameters of the ascending aorta, proximal descending thoracic aorta, and proximal abdominal aorta were determined from the M-mode images by using a maximum gradient method to determine the wall locations, for three levels corresponding to MRI slices (Fig. 1): (1) The proximal ascending aorta: The M-mode cursor line was placed at the level immediately distal to the aortic sinus. The aortic sinus is able to be identified in both US imaging and live MRI. (2) The initial segment of the descending thoracic aorta: The M-mode cursor line was placed at the level where the curved aortic arch became straight at the descending tho racic aorta. It is approximately the level of the junction between the atria and ventricles which was used as a land mark in both US and MRI. (3) The proximal abdominal aorta: The M-mode cursor line was placed at the level prior to the celiac artery (the first major branch from the abdominal aorta). The diameter of the abdominal aorta between the diaphragm and the celiac artery is relatively consistent; therefore, slight shift of slice level should not cause significant variation of diameter measurement. The results of the MR angiography were segmented manually and the vessel centerline (in three-dimensional (3D) space) was created by linking the geometric centers of the segmented vessel slices. A smooth spline representation of the centerline was cre ated using the Numpy spline interpolation routines (numpy.splrep, numpy.scipy.org). The microCT data was segmented using the level set method of the vascular modeling toolkit (VMTK, vmtk.org). We then used VMTK to determine vessel centerline and radius based on the perimeter of the vessel cross section (as opposed to maximum ra dius of the inscribed sphere). A perimeter based method was used because in some cases insufficient perfusion pressures resulted in slightly noncircular vessel cross sections. “Average” mouse centerlines were determined for both MR and microCT geometries as follows. First, the subject-specific centerlines were aligned using rigid body rotation and translation of the centerlines by minimizing the root mean square displacement of the minimum distance between centerlines. To accomplish this, a series of 1 mm spaced points along the centerline were chosen and for each point, the largest sphere that did not intersect the other centerline was taken as a measure of the distance between the cen terlines. The root mean square of the sum of these distances was taken as a measure of the displacement between the centerlines. This displacement was minimized over the six degrees of freedom of the rigid body rotation. The alignment achieved was independ ent of which centerline was chosen as the target. Next, the 3D locations of individual points on the average centerline were cal culated. For each abscissa value (distance from the aortic root along the centerline), the average 3D location of this point was calculated from the 14 mice. The average radius along the centerline was also determined for the CT data. Average branch centerlines were determined in a similar man ner to the aorta centerlines as described above. Additionally, using VMTK, branch vessel locations (distance along centerline from aortic root), in-plane vessel angle, and out-of-plane vessel angles were determined for each branch for each mouse. For each branch, the vessel location, in-plane angle, and out-of-plane angle were averaged across all 14 mice. This data subsequently was used to create both idealized subject-specific geometries and a mouse-averaged geometry. In-plane and out-of-plane angles were defined as the angles between the daughter vessel (the branch) and the parent vessel (the aorta before the branch) as seen in Fig. 2. Transactions of the ASME
Fig. 1 Graphic prescriptions from MR angiography determ ine the aortic centerline and for in vivo US M-m ode recordings to m easure aortic diam eters. Two localization images are required to localize the aortic arch and descending aorta images separately, with six slices (AAr_1-AA rJ>) acquired on aortic arch (upper left) and 16 slices (DAo_15-DAo_16) on descend ing aorta (low er left). On the right, B-mode US images along the centerline of the aorta at three cross-sectional levels are shown. US M-m ode recordings were made at the proximal ascending aorta im m ediately distal to the aortic sinus (AAr_1) (upper right), at the initial segm ent of the descending thoracic aorta where the curved aortic arch becam e straight (Dao„1) (m iddle right), and at the proximal abdom inal aorta prior to the celiac artery (Dao_15) (lower right). LV is the left ventricle and Li is the liver.
2.2 Creating Vessel Geometries. Five different geometries were created incorporating different degrees of individual anat omy, and with the aortic centerline obtained from either CT or MRI. The third geometry (c) is the proposed method and the other methods (a, b, d, and e) serve as comparisons. We will refer to these as geometry construction methods throughout the text. The five geometry construction methods are: (a) Subject-specific CT data from a single mouse, with no radius corrections. This method ignores all the MRI and US imaging data. (b) Subject-specific CT derived centerline from a single mouse, with a circular cross section and radius corrected to subjectspecific M-mode US measurements. (c) Subject-specific MRI derived centerlines from a single mouse, with circular cross section and radius corrected to Journal of Biomechanical Engineering
subject-specific M-mode US measurements, and mouseaveraged aortic arch branch centerlines, locations, and branch angles. At the conclusion of the paper, this is the preferred geometry. (d) Mouse-averaged CT derived centerline with mouseaveraged vessel radius, aortic arch branch centerlines, loca tions, and branch angles. This is a single representation of the average behavior of all of the aortas with no individual subject specific data remaining. (e) Mouse-averaged MRI derived centerline with mouseaveraged vessel radius, aortic arch branch centerlines, loca tions, and branch angles. The individual single mouse chosen for geometries a-c was selected based on high CT data quality as well as being represen tative of a mouse with aortic arch branch angles and locations OCTOBER 2014, Vol. 136 / 101008-3
Fig. 2 In-plane (left) and out-of-plane (right) angles shown for the first branch in an example mouse aorta. The right side view is normal to the left side view.
being approximately one standard deviation from the mean i.e., not representative of an average mouse. This single mouse allowed us to determine eventually that reliable CFD results could not be obtained if geometry from an average over all mice was used. The microCT and MRI geometries (geometries b and c) were constructed from the radius and centerline data using VMTK. Here, the mouse-specific CT radius was scaled so the radius at the proximal descending thoracic and abdominal aorta matched the cycle averaged radius measured with US at both locations. Radius data for the subject-specific MRI data was determined by using the average CT radius data along the centerline, also corrected with cycle averaged subject-specific US data at the proximal descending thoracic aorta and abdominal aorta. 2.2.1 Additional Mouse-Specific Geometries. After investi gating the five methods of constructing vessel geometry (a-e), and concluding that only b and c were sufficiently adequate methods, we chose five different mice and generated geometries by methods b and c for each of them. These geometries were constructed using methods b and c as described above but without including the arch and the branches in the arch. The arch was included only as a flow extension to ensure a realistic velocity profile at the proximal portion of the descending thoracic aorta (DTAo, as in Ref. [2]). 2.3 Flow Boundary Conditions. Traction free boundary con ditions were imposed at the descending aorta outlet while Womersley profiles were imposed at the three aortic arch branches according to the temporal flow rate waveform. Rigid, no slip walls, and a constant kinematic blood viscosity of 3.5 x 106 m2/s and a density of 1.06 g/cm3 were assumed. This assumption of a Newtonian fluid is appropriate given the high shear rates present in the mouse aorta as per Feintuch et al. [10]. We chose to impose AR flow conditions [1] in this study because it provides a more complicated velocity field and should result in a higher probability of identifying differences between the methods tested [12]. However, flow splits in the three major branches of the aortic arch in mice with AR are not yet known. In order to estimate the flow rates and waveforms at the three branches, Doppler US was completed in the innominate artery, the left common carotid artery, and the left subclavian artery in five mice with moderate AR. Using intensity weighted mean velocity tracing, the relative flow split to each branch was deter mined. In all five mice, in all three branches, despite retrograde flow in diastole in the aortic arch, there was very little retrograde flow in the branch arteries. Hence, in our simulations an idealized temporal flow waveform was used in each branch with the systolic portion having the same shape as the systolic portion of the wave form in the descending aorta scaled to ensure proper flow splits. The diastolic portion of the waveform was simply set to zero to reflect the lack of flow in diastole. The temporal flow waveform at the aortic root (inlet) was set to be the sum of the three branch
101008-4
/ Vol. 136, O C T O B E R 2014
waveforms and the waveform used by Hoi et al. in the descending aorta [2], thus ensuring that the flow waveform in the descending aorta was equal to that used by Hoi et al. At the aortic root we imposed a flat velocity inlet profile. Although this is not the ideal boundary condition for determining wall shear stress in the arch [5], it is a reasonable simplification and our focus in this study was on comparing methods, not determining specific wall shear stress patterns in the aortic arch of mice with AR. In addition, two cases of reduced AR severity were simulated in the mouse-averaged CT derived geometry. The purpose of these simulations was to determine the relative effects of AR severity compared to varying the geometry construction method. The temporal flow waveform imposed at the aortic root was the same except the diastolic portion of the waveform (i.e., the retro grade portion) was reduced such that the diastolic flow was a frac tion (70% and 85%) of the original value. The flow waveforms to the three branches remained unchanged. 2.4 CFD Simulations. ICEM-CFD was used to generate quadratic tetrahedral finite element meshes of approximately 1.000. 000 elements for the five different geometries. These 1.000. 000 elements correspond to roughly 8 x 106 linear tetrahedra with over 1 x 106 nodes. This is about 10 x the number of nodes reported being used by Trachet at al. [13] who report satis factory convergence. The unsteady Navier-Stokes equations were solved using the AR flow waveforms detailed above via a previ ously validated in-house solver [2,14]. In order to ensure that initial transients were not present in the solution, at least three cardiac cycles were simulated with at least 4800 time steps included per cardiac cycle to ensure stability. The choice of the number of nodes and the number of time steps is probably excessively large and more than is necessary. If a large number of CFD calculations were being contemplated, investiga tions into reducing this computer burden would need to be under taken. We chose to focus on RRT, as this measure correlated best with plaque patterns in a previous study of the AR model [2], Although other measures of wall shear such as SAWSS or TAWSS or OSI could have been used, we have shown previously that OSI or RRT are very similar and all the measures are highly correlated; therefore, we have chosen to show RRT as it is recom mended as a “robust single metric of low and oscillatory shear” [15], Himburg et al. [16] define RRT = T/ | J tw
