Nonlinear Anisotropic Stress Analysis of Anatomically Realistic Cerebral Aneurysms Baoshun Ma1*, Jia Lu2, Robert E Harbaugh3, Madhavan L. Raghavan1 1 Department of Biomedical Engineering, 1402 SC, University of Iowa, Iowa City, IA 52242 Tel: (319) 335-5704 E-mail: [email protected] 2 Department of Mechanical and Industrial Engineering, University of Iowa, Iowa City, IA 52242 Tel: (319) 335-6405 Email: [email protected] 3 Departments of Neurosurgery, and Engineering Science and Mechanics, Penn State University, Hershey, PA 17033 Tel: (717) 531-8807 E-mail: [email protected] * Current address: Department of Biomedical Engineering, Boston University, Boston, MA 02215 E-mail: [email protected] Submitted to: Journal of Biomechanical Engineering, Transactions of ASME Corresponding author: Madhavan L. Raghavan Department of Biomedical Engineering 1422 SC, University of Iowa Iowa City, IA 52242 Tel: (319) 335-5704 E-mail: [email protected] Fax: (319) 335-5631 1 Abstract Background. Static deformation analysis and estimation of wall stress distribution of patient- specific cerebral aneurysms can provide useful insights into the disease process and rupture. Method of Approach. The three-dimensional geometry of saccular cerebral aneurysms from twenty-seven patients (eighteen unruptured and nine ruptured) was reconstructed based on CTA images. The aneurysm wall tissue was modeled using a nonlinear, anisotropic, hyperelastic material model (Fung-type) which was incorporated in a user subroutine in ABAQUS. Effective material fiber orientations were assumed to align with principal surface curvatures. Static deformation of the aneurysm models were simulated assuming uniform wall thickness and internal pressure load of 100 mmHg. Results. The numerical analysis technique was validated by quantitative comparisons to results in the literature. For the patient-specific models, in-plane stresses in the aneurysm wall along both the stiff and weak fiber directions showed significant regional variations with the former being higher. The spatial maximum of stress ranged from as low as 0.30 MPa in a small aneurysm to as high as 1.06 MPa in a giant aneurysm. The patterns of distribution of stress, strain and surface curvature were found to be similar. Sensitivity analyses showed that the computed stress is mesh independent and not very sensitive to reasonable perturbations in model parameters, and the curvature-based criteria for fiber orientations tend to minimize the total elastic strain energy in the aneurysms wall. Within this small study population, there were no statistically significant differences in the spatial means and maximums of stress and strain values 2 between the ruptured and unruptured groups. However, the ratios between the stress components in the stiff and weak fiber directions were significantly higher in the ruptured group than those in the unruptured group. Conclusions. A methodology for nonlinear, anisotropic static deformation analysis of geometrically realistic aneurysms was developed, which can be used for a more accurate estimation of the stresses and strains than previous methods and to facilitate prospective studies on the role of stress in aneurysm rupture. Keywords: Intracranial aneurysm, biomechanics, hyperelasticity, stress, strain, curvature I. Introduction Cerebral aneurysms are localized dilatations of the intracranial arterial wall that usually occur on or near the circle of Willis. Rupture of intracranial aneurysms is the leading cause of subarachnoid hemorrhage (SAH) and often causes significant mortality and morbidity [1]. Solid biomechanical analysis of patient-specific cerebral aneurysms and estimations of wall stress distribution can provide useful insights into the disease process and help quantify geometric features within a physics-based framework. Investigations on solid biomechanics of Intra-Cranial Aneurysms (ICA) have been scarce unlike for abdominal aortic aneurysms [2-7]. On experimental investigations of tissue behavior, Scott et al. [8] reported in vitro pressure-volume relationships for ICA, assuming the lesions to be spherical. Steiger et al. [9] and Toth et al. [10] both performed uniaxial extension tests on excised ICA tissue specimen. These studies demonstrated the nonlinear behavior over finite 3 strains that are typical for soft biological tissues. Hsu et al. [11, 12] developed a multiaxial experimental system and studied the pressure-deformation behavior of two ICAs harvested as a whole from cadavers. Using their sophisticated test system, they were able to capture multi-axial deformations at multiple locations under varying pressures. Despite the rather small study population, these studies remain the most comprehensive characterization of elastic behavior of these lesions. They proposed the nonlinear theory of elastic membranes as an appropriate theoretical framework for saccular aneurysms. Later work by Kyriacou et al. [13] and Seshaiyer et al. [14, 15] found that the Fung-type anisotropic pseudostrain hyperelasticity model described the experimental results well, for which they reported estimated material parameter values. Reports on computational modeling of ICA solid mechanics has also been scarce. Simplified numerical static deformation analyses using a rubber-like [16] or linear material model [17, 18] have been reported. Analytical stress analysis [18, 19] using the law of Laplace have also been reported, but these were limited to spherical lesions. Kyriacou et al. [20] and Shah et al. [21] performed nonlinear, isotropic and anisotropic stress analyses of a class of axisymmetric saccular aneurysms assuming uniform internal pressure, uniform wall thickness and a Fung-type constitutive model for the wall material. They found that mechanical wall stress is affected by the shape of these lesions. Further, the axisymmetric models were used to study finite strain elastodynamics in ICA [22, 23]. Based on their results, the authors suggested that aneurysm expansion and rupture may not be explained by the limit point instability phenomenon that is known to occur during inflation of rubber balloons [24]. 4 It is noteworthy that previous solid mechanical analyses have stopped short of patient- specificity. The ability to perform patient-specific analyses depends on our ability to collect patient-specific information on aneurysm surface geometry, wall thickness and properties. The 3D aneurysm surface geometry is obtainable using current diagnostic tools such as computed tomography angiography (CTA) as reported recently by our group [25]. But patient-specific information on wall thickness and intrinsic material property (both of which can vary regionally) is virtually unattainable by current diagnostic tools. Indeed, even explant tissue-based population estimates for these properties are scarce. Within this state of limited knowledge, we believe that it is worthwhile to incorporate at least the obtainable information (namely, patient-specific geometry) into aneurysm-wall stress computations while making reasonable assumptions on the sparsely available pieces of information (namely, patient-specific thickness and intrinsic properties). Stress thus computed is unlikely to be highly accurate, but may still have two valuable applications: 1) It may help us develop experimentally testable hypotheses on how thickness and intrinsic properties may vary regionally or between patients and thereby facilitate a better understanding of disease pathogenesis; and 2) it may serve as a physics-sensitive index of patient-specific aneurysm surface geometry and thereby allow us to test whether geometry is a clinically measurable risk factor for rupture. In this study, we propose a methodology to compute the distribution of mechanical wall stresses in geometrically realistic models of ICAs using a nonlinear, anisotropic material model and apply it on a study population of aneurysm patients. We further demonstrate how the process can be effective toward the above two applications. 5 II. Material model 2.1 Constitutive model of cerebral aneurysms In the current study, the ICA wall tissue was modeled by the Fung-type strain energy density function as follows: ( ) w=c eQ −1 (Eq.1) Q = c E2 +c E2 +2c E E +c E2 1 11 2 22 3 11 22 4 12 where w is strain energy density function; E is component of the Green strain tensor; c , c , c , ij 1 2 c , and c are material model parameters; c and c are measures of material stiffness along the 3 4 1 2 stiff and weak material fiber directions and hence representative of the level of anisotropy; and c 4 is a measure of shear stiffness. Based on experimental pressure-inflation testing of two excised lesions, Seshaiyer et al. [14] proposed a functional form similar to (Eq.1), except for the shear term c E 2 (because they ignored in-plane shear during data analysis and subsequent modeling), 4 12 and reported best fit values for c , c , c and c based on their experimental data. A total of 11 1 2 3 sets of parameters were reported for five regions from two human aneurysms for pressures from 3 to 80, 120 and 150 mmHg, respectively. Of these, four sets of parameter values corresponded to non-atherosclerotic regions for pressures between 3 and 120 mmHg. In this current study, the averages of these four sets were assumed to represent the mean homogeneous material properties for all ICAs. Values of c and c were unequal, implying anisotropic properties for realistic 1 2 aneurysms. 6 We have included the shear term in our material model to accommodate in-plane shear that will inevitably exist – likely in diminished levels – owing to our use of arbitrary lesion geometry. Since there are no reported values for the parameter associated with in-plane shear, we assumed c = (c +c )/6 borrowing from the analogous relationship in linearized incompressible 4 1 2 elasticity (shear modulus = 1/3 of Young’s modulus). The following parameter values were assumed to represent population-averaged aneurysm properties and used for all anatomically realistic ICAs in this study: c = 0.28 N/mm2, c = 17.58, c = 12.19, c = 7.57, c = 4.96. 1 2 3 4 There are two non-essential differences between our choices and the report by Seshaiyer et al. [14] that are worth noting. They reported values for c in units of tension, not stress. We divided the c values by the mean wall thickness of each aneurysm from that report to convert it to unit of stress. Additionally, the values of c and c as reported by Seshaiyer et al. were 1 2 interchanged in the current study so that the stiff fiber direction corresponds to the local 1 direction in the local coordinate system. 2.2 Incorporating the Fung-type material model into ABAQUS The commercial structural analysis software ABAQUS (version 6.3, HKS Inc, Pawtucket, RI) was used for the stress analysis. Anisotropic finite elastic models have not been incorporated into the standard interface, so the user defined function “umat” was used to implement the Fung- type model (Eq. 1). The stress and strain were expressed in vector form as: E = (E , E , E , E , 11 22 33 12 E , E )T , σ = (σ , σ , σ , σ , σ , σ )T , where E is strain vector, E is engineering strain, σ is 13 23 11 22 33 12 13 23 ij stress vector, and σ is the component of Cauchy stress. Rewriting Q from (Eq. 1) in the form: ij 7 1 Q= ET ⋅A⋅E (Eq. 2) 2 where A is a 6x6 matrix with A =2c , A = A =2c , A =2c , A =0.5c and all other components 11 1 12 21 3 22 2 44 4 equal to zero. The Cauchy stress can be expressed as: 1 ∂w ceQ σ = F⋅ ⋅FT = F⋅(A⋅E)⋅FT (Eq. 3) J ∂E J where F is deformation gradient, J is the determinant of F. From (Eq. 1), the material Jacobean matrix (stiffness) in the current coordinate system can be expressed as: 1 ∂2w d = F⋅F: :FT ⋅FT (Eq. 4) J ∂E∂E Incorporating (Eq. 3) and (Eq. 4) into the “umat” subroutine will allow ABAQUS to use the Fung-type material model in computations. The subroutine was programmed in FORTRAN to work with general and axisymmetric membrane element, 3-D shell element, and 3-D solid element. The aneurysm wall was modeled using a shell element whose transverse shear stiffness was specified. For a homogeneous shell made of a linear, orthotropic elastic material, where the stiff material direction aligns with the element's local 1-direction, the transverse shear stiffness for the shell section should be (by ABAQUS/STANDARD user’s manual): 5 5 K = G t, K = G t, K =0.0 (Eq. 5) 11 6 13 22 6 23 12 where K , K and K are transverse shear stiffness values in the local 11, 22 and 12 directions 11 22 12 of the shell section, respectively; G and G are the material’s shear moduli in the out-of-plane 13 23 8 directions; and t is shell thickness. Although shear moduli are concepts in linear elasticity and not strictly applicable for finite elasticity, they can still provide a reasonable starting point. The shear moduli for the Fung-type material were estimated based on “virtual” uniaxial extension tests. Briefly, uniaxial extension of a unit cube (edge length of 1 mm, represented in ABAQUS by a 3-D, 8-node continuum element, C3D8) in the local 1 direction was simulated using the isotropic form (c = c ) of the Fung model (Eq. 1). At different stretch ratios, the stress and strain 1 2 values were used to estimate the Young’s modulus per engineering definition: E = σ /ε . Then 11 11 the shear modulus was estimated using the formula for an incompressible material: G = E/3. Based on (Eq. 5), the transverse shear stiffness values are: K = K = 5Gt/6 = 5E/18 (N/mm), 11 22 where t is wall thickness, which is 1.0 mm for the unit cube. As expected, the computed E (therefore K and K ) increased nonlinearly as the stretch ratio increases. In the current study 11 22 the values at stretch ratio of 1.01 was used (K = K = 2.4 N/mm). 11 22 2.3 Estimation of local material fiber directions In this study, it was assumed that collagen fiber is the main load-bearing structural component of the aneurysm wall, and transmurally there is a single effective collagen layer for the aneurysm wall. The fibers in the layer are oriented to give planar orthotropy with one strong and one week direction. In addition, the orientation of collagen fibers in this effective layer can vary from point to point within the aneurysm surface. For lack of a better alternative, we postulate that effective collagen fiber direction at a given point on the aneurismal surface is largely dependent on local stress [26]. Since wall stress is affected by surface curvature for thin- 9 walled structure, the local material fiber directions are likely related to surface curvature for stable aneurysms. For an axisymmetric membrane under internal pressure P, the stress resultants (tension) in the local 1 (meridional) and local 2 (circumferential) directions can be expressed as [27]: P T = (Eq. 6) 1 2k 2 P k T = (1− 1 ) (Eq. 7) 2 k 2k 2 2 where k and k are local principal curvatures for the meridional and circumferential directions, 1 2 respectively. Therefore, for an axisymmetric membrane under uniform pressure, the stresses and local surface curvatures are inherently related, i.e. the stresses are determined by local surface shape, and the principal stress directions coincide with principal curvature directions. From (Eq. 6) and (Eq. 7), P(k −k ) T −T = 1 2 (Eq. 8) 1 2 2k2 2 From (Eq. 8) it follows that if k > k , then T > T . Therefore, for axisymmetric convex surface, 1 2 1 2 direction of larger stress resultant corresponds to that of the larger surface curvature. Based on the above observations, we propose that a similar principle applies to general thin membranes of arbitrary shape. In this study we assume that: (1) local material fiber directions are orthogonal to each other; (2) stiff fiber aligns with the direction where stress tends to be larger; (3) larger stress corresponds to larger principal surface curvature. According to the above premises, it was assumed that locally at each point, the surface curvatures and effective
Description: