- Email: [email protected]
A reduced-order model of a patient-specific cerebral aneurysm for rapid evaluation and treatment planning
Journal Pre-proofs A Reduced-Order Model of a Patient-Specific Cerebral Aneurysm for Rapid Evaluation and Treatment Planning Suyue Han, Clemens M. Schirmer, Yahya Modarres-Sadeghi PII: DOI: Reference:
S0021-9290(20)30060-9 https://doi.org/10.1016/j.jbiomech.2020.109653 BM 109653
To appear in:
Journal of Biomechanics
Accepted Date:
21 January 2020
Please cite this article as: S. Han, C.M. Schirmer, Y. Modarres-Sadeghi, A Reduced-Order Model of a PatientSpecific Cerebral Aneurysm for Rapid Evaluation and Treatment Planning, Journal of Biomechanics (2020), doi: https://doi.org/10.1016/j.jbiomech.2020.109653
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Elsevier Ltd. All rights reserved.
A Reduced-Order Model of a Patient-Specific Cerebral Aneurysm for Rapid Evaluation and Treatment Planning
Suyue Han1, Clemens M. Schirmer2, Yahya Modarres-Sadeghi1* *Corresponding
author: [email protected], Tel:(413)545-2468, Fax:(413)545-1027
1Department
of Mechanical and Industrial Engineering,
University of Massachusetts, Amherst, Massachusetts, 01003
2Department
of Neurosurgery and Neuroscience Institute,
Geisinger, Wilkes-Barre, PA 18711
Abstract We discuss a computationally-efficient numerical Reduced Order Model (ROM) for patientspecific cerebral aneurysms, which can be used to examine the rupture-related hemodynamic parameters over a range of relevant physiological flow parameters, rapidly. The ROMs were derived using a Proper Orthogonal Decomposition (POD) technique, which also took advantage of the method of Snapshot POD. Initially a series of CFD training runs were performed, which were subsequently improved using a QR-decomposition technique to satisfy the various boundary conditions in physiological flow problems. We used the obtained ROMs to study the influence of Pulsatility Index (PI) on a patient-specific aneurysm’s Wall Shear Stress (WSS) and Oscillatory Shear Index (OSI). In addition, we discuss how each of the obtained high-energy
1
POD modes represents a separate significant flow pattern that is believed to influence the aneurysm’s WSS and OSI.
1. Introduction Image-based Computational Fluid Dynamics (CFD) has been widely used in attempts to identify factors leading to cerebral aneurysm rupture, understand flow characteristics of patientspecific aneurysms and support clinical decision-making and physical treatment planning (Cebral et al., 2011; Xiang et al., 2011; Sugiyama et al., 2012; Xiang et al., 2014; Chung and Cebral, 2015). It is well accepted that hemodynamics parameters, such as Wall Shear Stress (WSS) and Oscillatory Shear Index (OSI) are important indicators for aneurysm development, growth and rupture (Xiang et al.; 2014, Sano et al., 2017; Meng et al., 2014; Lu et al., 2011; Valen-Sendstad et al., 2014). Hemodynamic analysis of cerebral aneurysms through CFD simulation has been significantly limited by a variety of barriers such as the absence of patient-specific in-vivo flow input conditions (Chung and Cebral, 2015; Les et al., 2010; Quarteroni et al., 2000), limitations of image-based scanners (van Disseldorp et al., 2016) and high computational cost. Current correlations between hemodynamic parameters and aneurysm growth and rupture are still under debate (Meng et al., 2014; Valen-Sendstad et al., 2014; Zhou et al., 2017), which makes the clinical usage of hemodynamic parameters uncertain. These conflicting findings may come from oversimplified aneurysm geometries, idealized boundary conditions, inconsistent parameter definitions, small datasets, and intrinsic complexities in intracranial aneurysm growth and rupture (Cebral et al., 2011; Zhou et al., 2017). Therefore, larger datasets, better analyses, and increased understanding of hemodynamic-biologic mechanisms are necessary to build more 2
accurate predictive models for intracranial aneurysm risk assessment based on CFD results. However, due to the high computational cost and lack of access to noninvasive patient-specific boundary condition measurements, CFD simulations are difficult to apply in routine clinical procedures. Reduced Order Models (ROMs) have received a lot of attention in recent years, due to their ability to capture the original dynamics of the full-scale simulation with only a few degrees of freedom (Hosain and Fdhila, 2015). Because of the nonlinear convection terms and geometric variability, fluid systems always have strong nonlinearities, which cannot be treated by simple linearization (Lassila et al., 2014), especially when the unsteady flow is being studied and longterm transient behavior needs to be predicted using ROMs. As a model reduction method, Proper Orthogonal Decomposition (POD) is widely used, especially to extract ROMs for unsteady Navier-Stokes equations (Grinberg et al., 2009). The POD method essentially provides optimalordered, orthogonal bases to represent the data in a least-square sense (Algazi and Sakrison, 1969). The optimal low-dimensional approximation of the data is obtained by truncating the bases. There are several benefits in using POD as a model reduction method (Rathinam and Petzold, 2003; Pearson, 1901; Lumley, 1967; Pinnau, 2008): (i) The model decomposition of POD method is purely data-dependent and does not require any prior knowledge of the system behavior. Thus, the POD method is helpful to explore the underlying pattern in the data, even when prior knowledge is insufficient to provide model basis. (ii) The POD method can be easily combined with the Galerkin method and provide a powerful tool for building reduced order models for dynamical systems with very high dimensions or even infinite dimensions. (iii) Using the POD method, the dimensionality of a system is reduced by projecting the original system to a subspace. In the meanwhile, the nonlinearity of the original system is preserved, since the 3
subspace is composed of bases which inherit the special characteristics of the original system. (iv) Only basic matrix computation is required for both linear and nonlinear equations. In real-would problems, the dimensions of the system’s phase space could be very large, and the computational cost of generating the POD bases could be high. The so-called “method of snapshots” is then used (Chatterjee, 2000; Sirovich, 1987). A series of snapshots from the system’s trajectory are chosen with a certain time step. Therefore, instead of solving an eigenvalue problem which has dimensions of the full phase space, an eigenvalue problem with much smaller dimensions (the number of snapshots) will be conducted. The introduction of snapshot POD provides a powerful tool for the computation of POD bases. In a previous study (Chang et al., 2017), we considered a simplified, completely symmetric Abdominal Aortic Aneurysms (AAA) model and discussed a method for efficient reconstruction of flow pattern inside it. The snapshot POD plus Galerkin method was used to build an ROM, and the time variant boundary conditions were handled using QR-decomposition. The resulting ROM reconstructed the flow pattern and the WSS distribution over a range of inflow angles and pulsatility frequencies, efficiently, which allowed us to study the system across a range of different parameters without conducting repetitive CFD simulations. In the present study, we focus on a cerebral aneurysm with a patient-specific geometry and derive an ROM to be used for hemodynamic parametric studies. The physiological flow conditions over a range of different values of Pulsatility Index (PI) are examined. We have chosen the PI as the independent variable because: (i) recent studies have shown that, compared to the mean flow rate and the inflow waveform shape, the PI has a more significant impact on the hemodynamic parameters, especially on WSS-related parameters (Xiang and Siddiqui, 2014), and (ii) in clinical use, besides the fact that non-invasive measurements of patient-specific flow 4
waveform is usually not accessible, the PI can easily change after therapies such as clipping, coil embolization or adding a stent, and even due to a different measurement site. We have used a snapshot POD method to construct the ROMs, and we have applied QR-decomposition in order to handle the time-variant boundary conditions. We have conducted a comparison between CFD simulation and ROM reconstruction of the flow field and WSS distribution to evaluate the efficiency and accuracy of the ROM.
2. Method 2.1.
A Patient-Specific Cerebral Aneurysm Model The vascular modeling toolkit VMTK was used to segment the cerebral aneurysm from CT
angiography (Figure 1). The inlet and outlet surfaces were trimmed to be perpendicular to the centerline of the artery. The cerebral aneurysm’s surface model was then imported to ANSYS ICEM for volume mesh modeling. For the patient-specific geometry, a proper mesh is essential for accurate CFD simulation (Spiegel et al., 2009; Chung and Cebral, 2015; Robertson and Watton, 2012). It has been suggested that each patient-specific model requires its own grid convergence study for CFD accuracy (Spiegel et al., 2011; Hodis et al. 2012). For the present case, three-layer prism layer was used, and mesh convergence test was conducted to specify the critical mesh size of the model. The boundary layer mesh size was chosen to be very small (1×10-4 ~ 3×10-4 mm) to ensure the WSS accuracy, and the inner mesh size was chosen to be relatively large (2×10-3 ~ 6×10-3 mm) compared to the aneurysm wall prism layers in order to save computational cost. Delaunay mesh was used to adapt to the complex geometrical domains of patient-specific cerebral aneurysms (Weatherill, 1992).
5
2.2.
The Training CFD Run The first step toward building an ROM is to create a CFD dataset which contains enough
information and can be used to construct the POD modes. In this study, the fluid was assumed to be laminar and exhibit Newtonian behavior, the viscosity was assumed to be µ = 4×10-3 Pa.s and the density to be ρ = 1050 kg/m3. The CFD simulation was conducted in open-source CFD software OpenFOAM using incompressible transient flow solver pisoFoam on a 24-core workstation. In this case, the mean Reynolds number was Remin = 174, and the Euler scheme was used in the CFD simulation in OpenFOAM with a time step of 5×10-5 s. The convergence criteria were tolerances of 1×10-5 and 1×10-6 for the flow velocity and pressure fields, respectively. Pulsatility Index (PI) was chosen as the independent parameter mainly due to its strong correlation with WSS and OSI, and also because it has been shown to be one of the strongest independent risk factors for anterior communicating artery (ACoA) aneurysm (Kaspera et al., 2014). The PI is defined as: PI u max u min / u average ,
(1)
where umax is the maximum velocity of the waveform, umin is the minimum velocity of the waveform and uaverage is the average velocity of the waveform over one cycle. A typical ICA waveform (Ford et al., 2005) was used as the reference waveform of the inlet flow velocity (Figure 2(a)). During the training simulation, the reference waveform was scaled by a timedependent parameter, e(t):
uwave u t uaverage e(t ) uaverage ,
(2)
where, e(t ) 0.5 0.1 t ,
(3) 6
in which t is the simulation time. A 6-second CFD simulation containing 6 pulsatility cycles with a period of cardiac cycle, T = 1, s was conducted. The continuously changing e(t) resulted in varying PI over a range of 0.5 to 1.1 (Figure 2(b)), where PI is normally measured to be 0.5 to 1.19 in different patients (Gosling and King, 1974). In the CFD simulation, the mean inflow velocity stayed constant (umean = 0.38 m/s) as the PI value varied. A series of snapshots was then selected with a sample frequency of 500 Hz to be used as the training dataset for the ROM. In order to capture the correct temporal dynamics, the characteristics of the dynamical system should be preserved by the ROM. By varying the PI in multiple pulsatile cycles, we obtained information on transitional motion of the flow inside the aneurysm influenced by varying PI, which was needed for constructing an ROM with the POD plus Galerkin method. It ought to be mentioned here that, no converged result was expected within this six-second training CFD simulation, because the purpose of this training CFD dataset was to cover a range of continuously changing flow conditions to gather transitional flow information.
2.3.
Generating POD Modes After the training CFD was conducted, the collected snapshots were used to generate POD
modes using snapshot POD method enhanced by QR-decomposition. The QR decomposition is a method that is used to decompose a matrix A into A=QR, where Q is an orthogonal matrix and R is an upper triangle matrix. The POD modes can also be used to build the local ROMs using the local POD plus Galerkin Method. The accuracy of the ROM depends on the assumption that the trajectories of the original system can be approximated on a much lower-dimensional submanifold. As the first step of constructing a reliable ROM for the full-scale computational result, this assumption needs to be verified after the full-scale computational model has been 7
constructed. By collecting the snapshots of different trajectories of the full-scale system in one matrix, the singular values of this matrix can be found using the singular value decomposition. A proper ROM can be built with a small number of modes if the decay of the singular value is rapid enough.
2.4.
Constructing ROM To decide on the number of bases used to build the ROM, the Relative Information Content
(RIC) is calculated. RIC is defined as (Lassila et al., 2014): n
RIC
i
i
,
i 1 N
(4)
i 1
where n is the number of POD modes used, N is the total number of snapshots and σi are the eigenvalues calculated in the singular value decomposition. As the equation shows, the RIC is the ratio of the summation of the eigenvalues corresponding to the POD-bases selected to the summation of all the eigenvalues of the correlation matrix. By choosing a proper number for RIC, the most energetic POD modes will be retained. Another crucial point in constructing an ROM is to make sure that the Galerkin Method approximates the true dynamics. To do so, the root mean square (RMS) distance of the snapshots from the resulting POD manifold can be used as an error estimator (Rapun and Vega, 2010): RMS
1 N
N
2 i
.
(5)
i n 1
This relies on the assumption that the Galerkin approximation converges to the correct solution if an appropriate number of modes are used. The POD modes and the ROMs were created using an in-house Matlab-based code. 8
3. Results 3.1. The Reduced Order Model Following the steps discussed in the methods section, an ROM was built for the patientspecific aneurysm of Figure 1. Toward making these ROMs, first we plotted the eigenvalues for each POD mode and observed that they drop rapidly within the first 10 modes (Figure 3(a)). This suggests that a proper ROM can be built using these POD modes for this patient-specific aneurysm. Then to decide on the number of POD modes that we needed to construct the ROM, we plotted the RIC for each mode. As shown in Figure 3(b), in this patient-specific case 99.8% energy is captured by the first 5 POD modes, which suggests using only 5 modes in constructing the ROM could be sufficient. As the last test to see if the POD modes are appropriate for our Galerkin approximation, we plotted the ROM error of the Galerkin approximation versus the number of modes used. We found that for the patient-specific case considered here, the RMS error is smaller than 1.45×10-4 when 5 or more modes are used (Figure 3(c)), which again indicates that only 5 modes would be enough to construct the ROM with an acceptable error. We then used this ROM with 5 modes to construct the flow velocity and wall shear stress for this patient-specific aneurysm for different values of PI. Figure 4(a,b) shows the flow velocity and wall shear stress distribution at an arbitrary cross-section of the aneurysm for PI = 0.7, as an example. It is observed that the flow is asymmetric and that a large vortex exists inside the dome. Before more conclusions are drawn from this POD-based ROM, we will fist compare these results with those from an independent series of CFD simulations for validation.
9
3.2. Comparison with Independent CFD Simulations In order to compare the flow field constructed using our ROM and that calculated from traditional CFD simulations, we conducted a series of independent CFD simulations for given values of PI using the same solver that was used for the training CFD. These Independent simulations were run for twelve different values of PI, from PI=0.5 to PI=1.05, with increments of 0.05. For each case, the PI value was introduced to the ROM as a new boundary condition. For each value of PI, the ROM-based flow field was reconstructed and compared with that based on the direct CFD simulations. Results from a sample of these independent runs are shown in Figure 4 (c,d) for PI=0.7. A qualitative comparison reveals that the asymmetric flow structure, including the vortex inside the dome, which was observed in the ROM-based results are also observed in the direct CFD simulations. While qualitative agreements are desirable, a more rigorous quantitative comparison is needed, as we discuss in what follows. In order to conduct a quantitative comparison between the CFD and ROM results, we used four measures: (i)
To compare the flow velocity distribution, the kinetic energy of the flow field was calculated for both the ROM-based and CFD-based results as:
1 2 KE mu i i dV , 2 V
(6)
where mi are the mass of each cell contained in the volume mesh which is related to the cell volume and density, and ui are the flow velocities of each element.
(ii,iii) For the WSS comparison between the ROM-based and CFD-based results, the spatially averaged WSS was calculated as:
10
A dA i
WSS average
i
A
Atotal
,
(7)
where τi are the WSS values of each element on the aneurysm wall, Ai are the surface areas of the wall elements, and Atotal is the total surface area of the aneurysm wall. Since previous studies have also suggested a role for high WSS in aneurysm rupture (Cecchi et al., 2011), we report both Averaged Wall Shear Stress (AWSS, time- and spatiallyaveraged wall shear stress), and Maximum Wall Shear Stress (MWSS, maximum value of time-averaged wall shear stress) as measures for comparison.
(iv)
Oscillation Sheer Index (OSI) is a non-dimensional parameter representing the directional change of WSS during the cardiac cycle (He et al., 1996) and is calculated as: T dt i 0 1 OSI 1 T , 2 i dt 0
(8)
where T is the duration of the cardiac cycle. We report Maximum Oscillatory Shear Index (MOSI, the maximum WSS direction change during one cardiac cycle).
Figure 5 shows how the errors between ROM-based and CFD-based results change versus the number of POD modes for these four quantitative measures for PI=0.7. It was observed that with only 5 POD modes the errors for all the parameters are within 5%, which implies that the first 5 modes could represent the flow behavior inside the aneurysm. This is in agreement with our earlier observations in Section 3.1, where we investigated the validity of the POD modes that we used for constructing these ROM. The error increased slightly as modes higher than 5 are 11
included in the reconstruction. This slight increase (from a maximum of 5% error with 5 modes to a maximum error of 7% with higher modes) is because higher modes might contain more system noises and may behave as redundant information for the ROM. The error for all the other cases (different values of PI) remained within the same range.
3.3. Obtaining hemodynamic information from POD modes The information contained in the POD modes are valuable not only to be able to construct an ROM to represent the flow computationally efficiently, but also to obtain important hemodynamic information of the cerebral aneurysm. Figure 6 shows slices of the first 3 homogenous flow velocity modes of the ROM. The POD method sorts the eigenvalues based on their energy, which is the kinetic energy of the flow field. It is clear from the figure that each mode represents an important aspect of the flow: The first mode (Figure 6(a)) corresponds to the flow at the parent artery, the second mode (Figure 6(b)) corresponds to the major vortex pattern inside the aneurysm dome, and the third mode (Figure 6(c)) stands for the cross-neck flow (CNF). The other modes represent more complex flow behavior and overall contain 0.26% of the total energy. In the present analysis, both CNF and vortex pattern are captured by independent highenergy POD modes. This makes it possible for future studies to investigate their influences on the aneurysm growth and rupture by focusing on their relevant modes, rather than by considering the entire flow field.
4. Discussion The role of hemodynamics parameters in evaluation and treatment planning of cerebral aneurysm is well accepted. However, predictive tools for these parameters are not yet capable of 12
numerically efficient parametric studies, due to the high computational cost of direct CFD simulations. Reliable ROMs can be used as such tools and prepare the ground for clinical use of numerical simulations for patient-specific aneurysms. As an example, for a single cerebral aneurysm such as the one used in the present study, each independent CFD simulation for a single boundary condition takes 20 CPU hrs. That means that if a range of system parameters are to be investigated, for each set of parameters, a 20-CPU-hour run is needed. In the ROM technique, however, the computational cost is divided into 3 parts: (i) the time for getting the training snapshots, which is 20 CPU hrs (this is conducted once only for a parametric study), (ii) the time to use POD plus Galerkin method to construct the ROMs, which is 0.1 CPU hrs, and (iii) the time to use the ROM for each set of parameters, which is 8.3×10-3 CPU hrs. Now, if one decides to conduct a parametric study, the computational cost of examining different sets of parameters using traditional CFD simulation increases linearly with the number of sets (e.g., 2000 CPU hrs for 100 sets), whereas the computational cost of our ROM technique reaches a plateau after the training is conducted and the ROM is constructed, and it only takes 21 CPU hrs to simulate 100 and even more sets, which is considerably lower than the computational cost for a full order CFD simulation. Computational methods are necessary to study the physiological flow field within the aneurysm. Non-invasive measurement techniques which could give detailed flow conditions within the aneurysm are very limited and the resolution of these measurement techniques are too low to provide a precise flow velocity or WSS distribution in vivo. While numerical methods are necessary for these studies, they have several limitations, due to their high-computational times, their inability to handle very large databases of aneurysms, and their simplified geometries, to name a few. Since ROMs are computationally efficient, they have the potential to overcome 13
several of these limitations and make it possible to study large aneurysm databases, time varying and complex boundary conditions, and complex patient-specific aneurysm geometries. ROMs can provide precise WSS distributions with aneurysms under different boundary conditions (e.g., different incoming flow velocity, different geometries), which will make it possible to investigate the influence of variabilities in the system parameters on the biomechanics of aneurysm.
5. Conclusions We have presented computationally efficient ROMs for a patient-specific cerebral aneurysm model using POD plus Galerkin method enhanced by QR-decomposition. The derived ROM makes it possible to explore the role of different system parameters across a range of clinically relevant values on WSS and OSI without repetitive runs of CFD simulations, resulting in significant savings of computational resources. Comparisons with independent CFD simulations showed that a rather low number of POD modes are capable of reconstructing the flow field with errors less than 5%. The POD-based ROMs rely on the training CFD results. If the training CFD is not reliable, then the ROM-based results will not be reliable either, since the ROM gets all its information about the dynamics of the flow from the training CFD. If the CFD method used in the training run can resolve certain features of the flow, then the POD-based ROMs have the possibility of resolving the same features. Therefore, the POD-based ROMs should be thought as complementary methods to the direct CFD methods. The results of the present study show how significant features of the flow inside the aneurysm can be observed in the high-energy (i.e., the first 3) POD modes. The flow in the artery, 14
the vortex in the dome and the flow in the neck were clearly the main features in the first 3 POD modes. The fact that these high-energy POD modes corresponded to the main flow features in the aneurysm, means that each of these POD modes can be used independently to analyze the influence of each of the specific flow features on the resulting WSS and OSI distribution. This can further save both computational and analysis time by making it possible to focus on a specific feature of the flow (i.e., a specific POD mode) for a specific analysis, rather than the entire flow.
Conflict of Interest There is no conflict of interest to report.
References Algazi, V. and Sakrison, D., 1969. On the optimality of the Karhunen-Loève expansion (Corresp.). IEEE Transactions on Information Theory, 15(2), pp.319-321. Babiker, M.H., Gonzalez, L.F., Ryan, J., Albuquerque, F., Collins, D., Elvikis, A. and Frakes, D.H., 2012. Influence of stent configuration on cerebral aneurysm fluid dynamics. Journal of biomechanics, 45(3), pp.440-447. Baharoglu, M.I., Schirmer, C.M., Hoit, D.A., Gao, B.L. and Malek, A.M., 2010. Aneurysm inflow-angle as a discriminant for rupture in sidewall cerebral aneurysms: morphometric and computational fluid dynamic analysis. Stroke, 41(7), pp.1423-1430. Castro, M., Putman, C., Radaelli, A., Frangi, A. and Cebral, J., 2009. Hemodynamics and rupture of terminal cerebral aneurysms. Academic radiology, 16(10), pp.1201-1207. Cebral, J.R., Mut, F., Weir, J. and Putman, C., 2011. Quantitative characterization of the hemodynamic environment in ruptured and unruptured brain aneurysms. American Journal of Neuroradiology, 32(1), pp.145-151. Cecchi, E., Giglioli, C., Valente, S., Lazzeri, C., Gensini, G.F., Abbate, R. and Mannini, L., 2011. Role of hemodynamic shear stress in cardiovascular disease. Atherosclerosis, 214(2), pp.249-256. Chang, G.H., Schirmer, C.M. and Modarres-Sadeghi, Y., 2017. A reduced-order model for wall shear stress in abdominal aortic aneurysms by proper orthogonal decomposition. Journal of biomechanics, 54, pp.33-43. Chatterjee, A., 2000. An introduction to the proper orthogonal decomposition. Current science, pp.808817. Chung, B. and Cebral, J.R., 2015. CFD for evaluation and treatment planning of aneurysms: review of proposed clinical uses and their challenges. Annals of biomedical engineering, 43(1), pp.122-138.
15
Dolan, J.M., Meng, H., Sim, F.J. and Kolega, J., 2013. Differential gene expression by endothelial cells under positive and negative streamwise gradients of high wall shear stress. American Journal of Physiology-Cell Physiology, 305(8), pp.C854-C866. Ford, M.D., Alperin, N., Lee, S.H., Holdsworth, D.W. and Steinman, D.A., 2005. Characterization of volumetric flow rate waveforms in the normal internal carotid and vertebral arteries. Physiological measurement, 26(4), p.477. Fukunaga, K., 2013. Introduction to statistical pattern recognition. Elsevier. Gosling, R.G. and King, D.H., 1974. The Role of Measurement in Peripheral Vascular Surgery: Arterial Assessment by Doppler-shift Ultrasound. Grinberg, L., Yakhot, A. and Karniadakis, G.E., 2009. Analyzing transient turbulence in a stenosed carotid artery by proper orthogonal decomposition. Annals of biomedical engineering, 37(11), pp.2200-2217. Gunzburger, M.D., Peterson, J.S. and Shadid, J.N., 2007. Reduced-order modeling of time-dependent PDEs with multiple parameters in the boundary data. Computer methods in applied mechanics and engineering, 196(4-6), pp.1030-1047. He, X. and Ku, D.N., 1996. Pulsatile flow in the human left coronary artery bifurcation: average conditions. Journal of biomechanical engineering, 118(1), pp.74-82. Hodis, S., Uthamaraj, S., Smith, A.L., Dennis, K.D., Kallmes, D.F. and Dragomir-Daescu, D., 2012. Grid convergence errors in hemodynamic solution of patient-specific cerebral aneurysms. Journal of biomechanics, 45(16), pp.2907-2913. Hoi, Y., Meng, H., Woodward, S.H., Bendok, B.R., Hanel, R.A., Guterman, L.R. and Hopkins, L.N., 2004. Effects of arterial geometry on aneurysm growth: three-dimensional computational fluid dynamics study. Journal of neurosurgery, 101(4), pp.676-681. Hosain, M.L. and Fdhila, R.B., 2015. Literature review of accelerated CFD simulation methods towards online application. In Energy Procedia (Vol. 75, pp. 3307-3314). Kaspera, W., Ładziński, P., Larysz, P., Hebda, A., Ptaszkiewicz, K., Kopera, M. and Larysz, D., 2014. Morphological, hemodynamic, and clinical independent risk factors for anterior communicating artery aneurysms. Stroke, 45(10), pp.2906-2911. Lassila, T., Manzoni, A., Quarteroni, A. and Rozza, G., 2014. Model order reduction in fluid dynamics: challenges and perspectives. In Reduced Order Methods for modeling and computational reduction (pp. 235-273). Springer, Cham. Les, A.S., Shadden, S.C., Figueroa, C.A., Park, J.M., Tedesco, M.M., Herfkens, R.J., Dalman, R.L. and Taylor, C.A., 2010. Quantification of hemodynamics in abdominal aortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics. Annals of biomedical engineering, 38(4), pp.1288-1313. Lu, G., Huang, L., Zhang, X.L., Wang, S.Z., Hong, Y., Hu, Z. and Geng, D.Y., 2011. Influence of hemodynamic factors on rupture of intracranial aneurysms: patient-specific 3D mirror aneurysms model computational fluid dynamics simulation. American Journal of Neuroradiology, 32(7), pp.1255-1261. Lumley, J.L., 1967. The structure of inhomogeneous turbulent flows. Atmospheric turbulence and radio wave propagation. Meng, H., Tutino, V.M., Xiang, J. and Siddiqui, A., 2014. High WSS or low WSS? Complex interactions of hemodynamics with intracranial aneurysm initiation, growth, and rupture: toward a unifying hypothesis. American Journal of Neuroradiology, 35(7), pp.1254-1262. Mut, F., Löhner, R., Chien, A., Tateshima, S., Viñuela, F., Putman, C. and Cebral, J.R., 2011. Computational hemodynamics framework for the analysis of cerebral aneurysms. International journal for numerical methods in biomedical engineering, 27(6), pp.822-839. Pearson, K., 1901. LIII. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), pp.559-572.
16
Pinnau, R., 2008. Model reduction via proper orthogonal decomposition. In Model Order Reduction: Theory, Research Aspects and Applications (pp. 95-109). Springer, Berlin, Heidelberg. Quarteroni, A., Tuveri, M. and Veneziani, A., 2000. Computational vascular fluid dynamics: problems, models and methods. Computing and Visualization in Science, 2(4), pp.163-197. Rapún, M.L. and Vega, J.M., 2010. Reduced order models based on local POD plus Galerkin projection. Journal of Computational Physics, 229(8), pp.3046-3063. Rathinam, M. and Petzold, L.R., 2002. Dynamic iteration using reduced order models: a method for simulation of large scale modular systems. SIAM journal on Numerical Analysis, 40(4), pp.14461474. Rathinam, M. and Petzold, L.R., 2003. A new look at proper orthogonal decomposition. SIAM Journal on Numerical Analysis, 41(5), pp.1893-1925. Robertson, A.M. and Watton, P.N., 2012. Computational fluid dynamics in aneurysm research: critical reflections, future directions. Rosenfeld, A., 1976. Digital picture processing. Academic press. Sano, T., Ishida, F., Tsuji, M., Furukawa, K., Shimosaka, S. and Suzuki, H., 2017. Hemodynamic differences between ruptured and unruptured cerebral aneurysms simultaneously existing in the same location: 2 case reports and proposal of a novel parameter oscillatory velocity index. World neurosurgery, 98, pp.868-e5. Sirovich, L., 1987. Turbulence and the dynamics of coherent structures. I. Coherent structures. Quarterly of applied mathematics, 45(3), pp.561-571. Spiegel, M., Redel, T., Zhang, Y.J., Struffert, T., Hornegger, J., Grossman, R.G., Doerfler, A. and Karmonik, C., 2009, September. Tetrahedral and polyhedral mesh evaluation for cerebral hemodynamic simulation—a comparison. In Engineering in Medicine and Biology Society, 2009. EMBC 2009. Annual International Conference of the IEEE (pp. 2787-2790). IEEE. Spiegel, M., Redel, T., Zhang, Y.J., Struffert, T., Hornegger, J., Grossman, R.G., Doerfler, A. and Karmonik, C., 2011. Tetrahedral vs. polyhedral mesh size evaluation on flow velocity and wall shear stress for cerebral hemodynamic simulation. Computer methods in biomechanics and biomedical engineering, 14(01), pp.9-22. Sugiyama, S.I., Meng, H., Funamoto, K., Inoue, T., Fujimura, M., Nakayama, T., Omodaka, S., Shimizu, H., Takahashi, A. and Tominaga, T., 2012. Hemodynamic analysis of growing intracranial aneurysms arising from a posterior inferior cerebellar artery. World neurosurgery, 78(5), pp.462-468. Valen-Sendstad, K. and Steinman, D.A., 2014. Mind the gap: impact of computational fluid dynamics solution strategy on prediction of intracranial aneurysm hemodynamics and rupture status indicators. American Journal of Neuroradiology, 35(3), pp.536-543. van Disseldorp, E.M., Hobelman, K.H., Petterson, N.J., van de Vosse, F.N., van Sambeek, M.R. and Lopata, R.G., 2016. Influence of limited field-of-view on wall stress analysis in abdominal aortic aneurysms. Journal of biomechanics, 49(12), pp.2405-2412. Weatherill, N.P., 1992. Delaunay triangulation in computational fluid dynamics. Computers & Mathematics with Applications, 24(5-6), pp.129-150. Xiang, J., Natarajan, S.K., Tremmel, M., Ma, D., Mocco, J., Hopkins, L.N., Siddiqui, A.H., Levy, E.I. and Meng, H., 2011. Hemodynamic–morphologic discriminants for intracranial aneurysm rupture. Stroke, 42(1), pp.144-152. Xiang, J., Siddiqui, A.H. and Meng, H., 2014. The effect of inlet waveforms on computational hemodynamics of patient-specific intracranial aneurysms. Journal of biomechanics, 47(16), pp.38823890. Xiang, J., Tutino, V.M., Snyder, K.V. and Meng, H., 2014. CFD: computational fluid dynamics or confounding factor dissemination? The role of hemodynamics in intracranial aneurysm rupture risk assessment. American Journal of Neuroradiology, 35(10), pp.1849-1857. Zhou, G., Zhu, Y., Yin, Y., Su, M. and Li, M., 2017. Association of wall shear stress with intracranial aneurysm rupture: systematic review and meta-analysis. Scientific reports, 7(1), p.5331.
17
Figure 1. Patient-specific cerebral aneurysm model derived from computed tomographic angiography.
18
Figure 2. (a) Reference waveform from an averaged cerebrovascular waveform of 17 patients (Ford et al., 2005), and (b) six cardiac cycles of scaled waveforms with PI
[0.5,1.1].
19
Figure 3. (a) Eigenvalue versus the mode number, and (b) RIC, and (c) RMS error versus the number of modes used in the POD-based reconstruction for the patient-specific model used here.
20
Figure 4. Flow field for PI=0.7 at an arbitrary cross-section (left) and wall shear stress (WSS) (right) of the patientspecific aneurysm studied here based on (a,b) the ROM-based results, and (c,d) direct CFD simulations.
21
Figure 5. The error between the CFD-based and ROM-based results versus total number of modes used for PI=0.7.
22
Figure 6. Slices of homogenous flow velocity modes 1 – 3 of the ROM.
23
S0021-9290(20)30060-9 https://doi.org/10.1016/j.jbiomech.2020.109653 BM 109653
To appear in:
Journal of Biomechanics
Accepted Date:
21 January 2020
Please cite this article as: S. Han, C.M. Schirmer, Y. Modarres-Sadeghi, A Reduced-Order Model of a PatientSpecific Cerebral Aneurysm for Rapid Evaluation and Treatment Planning, Journal of Biomechanics (2020), doi: https://doi.org/10.1016/j.jbiomech.2020.109653
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Elsevier Ltd. All rights reserved.
A Reduced-Order Model of a Patient-Specific Cerebral Aneurysm for Rapid Evaluation and Treatment Planning
Suyue Han1, Clemens M. Schirmer2, Yahya Modarres-Sadeghi1* *Corresponding
author: [email protected], Tel:(413)545-2468, Fax:(413)545-1027
1Department
of Mechanical and Industrial Engineering,
University of Massachusetts, Amherst, Massachusetts, 01003
2Department
of Neurosurgery and Neuroscience Institute,
Geisinger, Wilkes-Barre, PA 18711
Abstract We discuss a computationally-efficient numerical Reduced Order Model (ROM) for patientspecific cerebral aneurysms, which can be used to examine the rupture-related hemodynamic parameters over a range of relevant physiological flow parameters, rapidly. The ROMs were derived using a Proper Orthogonal Decomposition (POD) technique, which also took advantage of the method of Snapshot POD. Initially a series of CFD training runs were performed, which were subsequently improved using a QR-decomposition technique to satisfy the various boundary conditions in physiological flow problems. We used the obtained ROMs to study the influence of Pulsatility Index (PI) on a patient-specific aneurysm’s Wall Shear Stress (WSS) and Oscillatory Shear Index (OSI). In addition, we discuss how each of the obtained high-energy
1
POD modes represents a separate significant flow pattern that is believed to influence the aneurysm’s WSS and OSI.
1. Introduction Image-based Computational Fluid Dynamics (CFD) has been widely used in attempts to identify factors leading to cerebral aneurysm rupture, understand flow characteristics of patientspecific aneurysms and support clinical decision-making and physical treatment planning (Cebral et al., 2011; Xiang et al., 2011; Sugiyama et al., 2012; Xiang et al., 2014; Chung and Cebral, 2015). It is well accepted that hemodynamics parameters, such as Wall Shear Stress (WSS) and Oscillatory Shear Index (OSI) are important indicators for aneurysm development, growth and rupture (Xiang et al.; 2014, Sano et al., 2017; Meng et al., 2014; Lu et al., 2011; Valen-Sendstad et al., 2014). Hemodynamic analysis of cerebral aneurysms through CFD simulation has been significantly limited by a variety of barriers such as the absence of patient-specific in-vivo flow input conditions (Chung and Cebral, 2015; Les et al., 2010; Quarteroni et al., 2000), limitations of image-based scanners (van Disseldorp et al., 2016) and high computational cost. Current correlations between hemodynamic parameters and aneurysm growth and rupture are still under debate (Meng et al., 2014; Valen-Sendstad et al., 2014; Zhou et al., 2017), which makes the clinical usage of hemodynamic parameters uncertain. These conflicting findings may come from oversimplified aneurysm geometries, idealized boundary conditions, inconsistent parameter definitions, small datasets, and intrinsic complexities in intracranial aneurysm growth and rupture (Cebral et al., 2011; Zhou et al., 2017). Therefore, larger datasets, better analyses, and increased understanding of hemodynamic-biologic mechanisms are necessary to build more 2
accurate predictive models for intracranial aneurysm risk assessment based on CFD results. However, due to the high computational cost and lack of access to noninvasive patient-specific boundary condition measurements, CFD simulations are difficult to apply in routine clinical procedures. Reduced Order Models (ROMs) have received a lot of attention in recent years, due to their ability to capture the original dynamics of the full-scale simulation with only a few degrees of freedom (Hosain and Fdhila, 2015). Because of the nonlinear convection terms and geometric variability, fluid systems always have strong nonlinearities, which cannot be treated by simple linearization (Lassila et al., 2014), especially when the unsteady flow is being studied and longterm transient behavior needs to be predicted using ROMs. As a model reduction method, Proper Orthogonal Decomposition (POD) is widely used, especially to extract ROMs for unsteady Navier-Stokes equations (Grinberg et al., 2009). The POD method essentially provides optimalordered, orthogonal bases to represent the data in a least-square sense (Algazi and Sakrison, 1969). The optimal low-dimensional approximation of the data is obtained by truncating the bases. There are several benefits in using POD as a model reduction method (Rathinam and Petzold, 2003; Pearson, 1901; Lumley, 1967; Pinnau, 2008): (i) The model decomposition of POD method is purely data-dependent and does not require any prior knowledge of the system behavior. Thus, the POD method is helpful to explore the underlying pattern in the data, even when prior knowledge is insufficient to provide model basis. (ii) The POD method can be easily combined with the Galerkin method and provide a powerful tool for building reduced order models for dynamical systems with very high dimensions or even infinite dimensions. (iii) Using the POD method, the dimensionality of a system is reduced by projecting the original system to a subspace. In the meanwhile, the nonlinearity of the original system is preserved, since the 3
subspace is composed of bases which inherit the special characteristics of the original system. (iv) Only basic matrix computation is required for both linear and nonlinear equations. In real-would problems, the dimensions of the system’s phase space could be very large, and the computational cost of generating the POD bases could be high. The so-called “method of snapshots” is then used (Chatterjee, 2000; Sirovich, 1987). A series of snapshots from the system’s trajectory are chosen with a certain time step. Therefore, instead of solving an eigenvalue problem which has dimensions of the full phase space, an eigenvalue problem with much smaller dimensions (the number of snapshots) will be conducted. The introduction of snapshot POD provides a powerful tool for the computation of POD bases. In a previous study (Chang et al., 2017), we considered a simplified, completely symmetric Abdominal Aortic Aneurysms (AAA) model and discussed a method for efficient reconstruction of flow pattern inside it. The snapshot POD plus Galerkin method was used to build an ROM, and the time variant boundary conditions were handled using QR-decomposition. The resulting ROM reconstructed the flow pattern and the WSS distribution over a range of inflow angles and pulsatility frequencies, efficiently, which allowed us to study the system across a range of different parameters without conducting repetitive CFD simulations. In the present study, we focus on a cerebral aneurysm with a patient-specific geometry and derive an ROM to be used for hemodynamic parametric studies. The physiological flow conditions over a range of different values of Pulsatility Index (PI) are examined. We have chosen the PI as the independent variable because: (i) recent studies have shown that, compared to the mean flow rate and the inflow waveform shape, the PI has a more significant impact on the hemodynamic parameters, especially on WSS-related parameters (Xiang and Siddiqui, 2014), and (ii) in clinical use, besides the fact that non-invasive measurements of patient-specific flow 4
waveform is usually not accessible, the PI can easily change after therapies such as clipping, coil embolization or adding a stent, and even due to a different measurement site. We have used a snapshot POD method to construct the ROMs, and we have applied QR-decomposition in order to handle the time-variant boundary conditions. We have conducted a comparison between CFD simulation and ROM reconstruction of the flow field and WSS distribution to evaluate the efficiency and accuracy of the ROM.
2. Method 2.1.
A Patient-Specific Cerebral Aneurysm Model The vascular modeling toolkit VMTK was used to segment the cerebral aneurysm from CT
angiography (Figure 1). The inlet and outlet surfaces were trimmed to be perpendicular to the centerline of the artery. The cerebral aneurysm’s surface model was then imported to ANSYS ICEM for volume mesh modeling. For the patient-specific geometry, a proper mesh is essential for accurate CFD simulation (Spiegel et al., 2009; Chung and Cebral, 2015; Robertson and Watton, 2012). It has been suggested that each patient-specific model requires its own grid convergence study for CFD accuracy (Spiegel et al., 2011; Hodis et al. 2012). For the present case, three-layer prism layer was used, and mesh convergence test was conducted to specify the critical mesh size of the model. The boundary layer mesh size was chosen to be very small (1×10-4 ~ 3×10-4 mm) to ensure the WSS accuracy, and the inner mesh size was chosen to be relatively large (2×10-3 ~ 6×10-3 mm) compared to the aneurysm wall prism layers in order to save computational cost. Delaunay mesh was used to adapt to the complex geometrical domains of patient-specific cerebral aneurysms (Weatherill, 1992).
5
2.2.
The Training CFD Run The first step toward building an ROM is to create a CFD dataset which contains enough
information and can be used to construct the POD modes. In this study, the fluid was assumed to be laminar and exhibit Newtonian behavior, the viscosity was assumed to be µ = 4×10-3 Pa.s and the density to be ρ = 1050 kg/m3. The CFD simulation was conducted in open-source CFD software OpenFOAM using incompressible transient flow solver pisoFoam on a 24-core workstation. In this case, the mean Reynolds number was Remin = 174, and the Euler scheme was used in the CFD simulation in OpenFOAM with a time step of 5×10-5 s. The convergence criteria were tolerances of 1×10-5 and 1×10-6 for the flow velocity and pressure fields, respectively. Pulsatility Index (PI) was chosen as the independent parameter mainly due to its strong correlation with WSS and OSI, and also because it has been shown to be one of the strongest independent risk factors for anterior communicating artery (ACoA) aneurysm (Kaspera et al., 2014). The PI is defined as: PI u max u min / u average ,
(1)
where umax is the maximum velocity of the waveform, umin is the minimum velocity of the waveform and uaverage is the average velocity of the waveform over one cycle. A typical ICA waveform (Ford et al., 2005) was used as the reference waveform of the inlet flow velocity (Figure 2(a)). During the training simulation, the reference waveform was scaled by a timedependent parameter, e(t):
uwave u t uaverage e(t ) uaverage ,
(2)
where, e(t ) 0.5 0.1 t ,
(3) 6
in which t is the simulation time. A 6-second CFD simulation containing 6 pulsatility cycles with a period of cardiac cycle, T = 1, s was conducted. The continuously changing e(t) resulted in varying PI over a range of 0.5 to 1.1 (Figure 2(b)), where PI is normally measured to be 0.5 to 1.19 in different patients (Gosling and King, 1974). In the CFD simulation, the mean inflow velocity stayed constant (umean = 0.38 m/s) as the PI value varied. A series of snapshots was then selected with a sample frequency of 500 Hz to be used as the training dataset for the ROM. In order to capture the correct temporal dynamics, the characteristics of the dynamical system should be preserved by the ROM. By varying the PI in multiple pulsatile cycles, we obtained information on transitional motion of the flow inside the aneurysm influenced by varying PI, which was needed for constructing an ROM with the POD plus Galerkin method. It ought to be mentioned here that, no converged result was expected within this six-second training CFD simulation, because the purpose of this training CFD dataset was to cover a range of continuously changing flow conditions to gather transitional flow information.
2.3.
Generating POD Modes After the training CFD was conducted, the collected snapshots were used to generate POD
modes using snapshot POD method enhanced by QR-decomposition. The QR decomposition is a method that is used to decompose a matrix A into A=QR, where Q is an orthogonal matrix and R is an upper triangle matrix. The POD modes can also be used to build the local ROMs using the local POD plus Galerkin Method. The accuracy of the ROM depends on the assumption that the trajectories of the original system can be approximated on a much lower-dimensional submanifold. As the first step of constructing a reliable ROM for the full-scale computational result, this assumption needs to be verified after the full-scale computational model has been 7
constructed. By collecting the snapshots of different trajectories of the full-scale system in one matrix, the singular values of this matrix can be found using the singular value decomposition. A proper ROM can be built with a small number of modes if the decay of the singular value is rapid enough.
2.4.
Constructing ROM To decide on the number of bases used to build the ROM, the Relative Information Content
(RIC) is calculated. RIC is defined as (Lassila et al., 2014): n
RIC
i
i
,
i 1 N
(4)
i 1
where n is the number of POD modes used, N is the total number of snapshots and σi are the eigenvalues calculated in the singular value decomposition. As the equation shows, the RIC is the ratio of the summation of the eigenvalues corresponding to the POD-bases selected to the summation of all the eigenvalues of the correlation matrix. By choosing a proper number for RIC, the most energetic POD modes will be retained. Another crucial point in constructing an ROM is to make sure that the Galerkin Method approximates the true dynamics. To do so, the root mean square (RMS) distance of the snapshots from the resulting POD manifold can be used as an error estimator (Rapun and Vega, 2010): RMS
1 N
N
2 i
.
(5)
i n 1
This relies on the assumption that the Galerkin approximation converges to the correct solution if an appropriate number of modes are used. The POD modes and the ROMs were created using an in-house Matlab-based code. 8
3. Results 3.1. The Reduced Order Model Following the steps discussed in the methods section, an ROM was built for the patientspecific aneurysm of Figure 1. Toward making these ROMs, first we plotted the eigenvalues for each POD mode and observed that they drop rapidly within the first 10 modes (Figure 3(a)). This suggests that a proper ROM can be built using these POD modes for this patient-specific aneurysm. Then to decide on the number of POD modes that we needed to construct the ROM, we plotted the RIC for each mode. As shown in Figure 3(b), in this patient-specific case 99.8% energy is captured by the first 5 POD modes, which suggests using only 5 modes in constructing the ROM could be sufficient. As the last test to see if the POD modes are appropriate for our Galerkin approximation, we plotted the ROM error of the Galerkin approximation versus the number of modes used. We found that for the patient-specific case considered here, the RMS error is smaller than 1.45×10-4 when 5 or more modes are used (Figure 3(c)), which again indicates that only 5 modes would be enough to construct the ROM with an acceptable error. We then used this ROM with 5 modes to construct the flow velocity and wall shear stress for this patient-specific aneurysm for different values of PI. Figure 4(a,b) shows the flow velocity and wall shear stress distribution at an arbitrary cross-section of the aneurysm for PI = 0.7, as an example. It is observed that the flow is asymmetric and that a large vortex exists inside the dome. Before more conclusions are drawn from this POD-based ROM, we will fist compare these results with those from an independent series of CFD simulations for validation.
9
3.2. Comparison with Independent CFD Simulations In order to compare the flow field constructed using our ROM and that calculated from traditional CFD simulations, we conducted a series of independent CFD simulations for given values of PI using the same solver that was used for the training CFD. These Independent simulations were run for twelve different values of PI, from PI=0.5 to PI=1.05, with increments of 0.05. For each case, the PI value was introduced to the ROM as a new boundary condition. For each value of PI, the ROM-based flow field was reconstructed and compared with that based on the direct CFD simulations. Results from a sample of these independent runs are shown in Figure 4 (c,d) for PI=0.7. A qualitative comparison reveals that the asymmetric flow structure, including the vortex inside the dome, which was observed in the ROM-based results are also observed in the direct CFD simulations. While qualitative agreements are desirable, a more rigorous quantitative comparison is needed, as we discuss in what follows. In order to conduct a quantitative comparison between the CFD and ROM results, we used four measures: (i)
To compare the flow velocity distribution, the kinetic energy of the flow field was calculated for both the ROM-based and CFD-based results as:
1 2 KE mu i i dV , 2 V
(6)
where mi are the mass of each cell contained in the volume mesh which is related to the cell volume and density, and ui are the flow velocities of each element.
(ii,iii) For the WSS comparison between the ROM-based and CFD-based results, the spatially averaged WSS was calculated as:
10
A dA i
WSS average
i
A
Atotal
,
(7)
where τi are the WSS values of each element on the aneurysm wall, Ai are the surface areas of the wall elements, and Atotal is the total surface area of the aneurysm wall. Since previous studies have also suggested a role for high WSS in aneurysm rupture (Cecchi et al., 2011), we report both Averaged Wall Shear Stress (AWSS, time- and spatiallyaveraged wall shear stress), and Maximum Wall Shear Stress (MWSS, maximum value of time-averaged wall shear stress) as measures for comparison.
(iv)
Oscillation Sheer Index (OSI) is a non-dimensional parameter representing the directional change of WSS during the cardiac cycle (He et al., 1996) and is calculated as: T dt i 0 1 OSI 1 T , 2 i dt 0
(8)
where T is the duration of the cardiac cycle. We report Maximum Oscillatory Shear Index (MOSI, the maximum WSS direction change during one cardiac cycle).
Figure 5 shows how the errors between ROM-based and CFD-based results change versus the number of POD modes for these four quantitative measures for PI=0.7. It was observed that with only 5 POD modes the errors for all the parameters are within 5%, which implies that the first 5 modes could represent the flow behavior inside the aneurysm. This is in agreement with our earlier observations in Section 3.1, where we investigated the validity of the POD modes that we used for constructing these ROM. The error increased slightly as modes higher than 5 are 11
included in the reconstruction. This slight increase (from a maximum of 5% error with 5 modes to a maximum error of 7% with higher modes) is because higher modes might contain more system noises and may behave as redundant information for the ROM. The error for all the other cases (different values of PI) remained within the same range.
3.3. Obtaining hemodynamic information from POD modes The information contained in the POD modes are valuable not only to be able to construct an ROM to represent the flow computationally efficiently, but also to obtain important hemodynamic information of the cerebral aneurysm. Figure 6 shows slices of the first 3 homogenous flow velocity modes of the ROM. The POD method sorts the eigenvalues based on their energy, which is the kinetic energy of the flow field. It is clear from the figure that each mode represents an important aspect of the flow: The first mode (Figure 6(a)) corresponds to the flow at the parent artery, the second mode (Figure 6(b)) corresponds to the major vortex pattern inside the aneurysm dome, and the third mode (Figure 6(c)) stands for the cross-neck flow (CNF). The other modes represent more complex flow behavior and overall contain 0.26% of the total energy. In the present analysis, both CNF and vortex pattern are captured by independent highenergy POD modes. This makes it possible for future studies to investigate their influences on the aneurysm growth and rupture by focusing on their relevant modes, rather than by considering the entire flow field.
4. Discussion The role of hemodynamics parameters in evaluation and treatment planning of cerebral aneurysm is well accepted. However, predictive tools for these parameters are not yet capable of 12
numerically efficient parametric studies, due to the high computational cost of direct CFD simulations. Reliable ROMs can be used as such tools and prepare the ground for clinical use of numerical simulations for patient-specific aneurysms. As an example, for a single cerebral aneurysm such as the one used in the present study, each independent CFD simulation for a single boundary condition takes 20 CPU hrs. That means that if a range of system parameters are to be investigated, for each set of parameters, a 20-CPU-hour run is needed. In the ROM technique, however, the computational cost is divided into 3 parts: (i) the time for getting the training snapshots, which is 20 CPU hrs (this is conducted once only for a parametric study), (ii) the time to use POD plus Galerkin method to construct the ROMs, which is 0.1 CPU hrs, and (iii) the time to use the ROM for each set of parameters, which is 8.3×10-3 CPU hrs. Now, if one decides to conduct a parametric study, the computational cost of examining different sets of parameters using traditional CFD simulation increases linearly with the number of sets (e.g., 2000 CPU hrs for 100 sets), whereas the computational cost of our ROM technique reaches a plateau after the training is conducted and the ROM is constructed, and it only takes 21 CPU hrs to simulate 100 and even more sets, which is considerably lower than the computational cost for a full order CFD simulation. Computational methods are necessary to study the physiological flow field within the aneurysm. Non-invasive measurement techniques which could give detailed flow conditions within the aneurysm are very limited and the resolution of these measurement techniques are too low to provide a precise flow velocity or WSS distribution in vivo. While numerical methods are necessary for these studies, they have several limitations, due to their high-computational times, their inability to handle very large databases of aneurysms, and their simplified geometries, to name a few. Since ROMs are computationally efficient, they have the potential to overcome 13
several of these limitations and make it possible to study large aneurysm databases, time varying and complex boundary conditions, and complex patient-specific aneurysm geometries. ROMs can provide precise WSS distributions with aneurysms under different boundary conditions (e.g., different incoming flow velocity, different geometries), which will make it possible to investigate the influence of variabilities in the system parameters on the biomechanics of aneurysm.
5. Conclusions We have presented computationally efficient ROMs for a patient-specific cerebral aneurysm model using POD plus Galerkin method enhanced by QR-decomposition. The derived ROM makes it possible to explore the role of different system parameters across a range of clinically relevant values on WSS and OSI without repetitive runs of CFD simulations, resulting in significant savings of computational resources. Comparisons with independent CFD simulations showed that a rather low number of POD modes are capable of reconstructing the flow field with errors less than 5%. The POD-based ROMs rely on the training CFD results. If the training CFD is not reliable, then the ROM-based results will not be reliable either, since the ROM gets all its information about the dynamics of the flow from the training CFD. If the CFD method used in the training run can resolve certain features of the flow, then the POD-based ROMs have the possibility of resolving the same features. Therefore, the POD-based ROMs should be thought as complementary methods to the direct CFD methods. The results of the present study show how significant features of the flow inside the aneurysm can be observed in the high-energy (i.e., the first 3) POD modes. The flow in the artery, 14
the vortex in the dome and the flow in the neck were clearly the main features in the first 3 POD modes. The fact that these high-energy POD modes corresponded to the main flow features in the aneurysm, means that each of these POD modes can be used independently to analyze the influence of each of the specific flow features on the resulting WSS and OSI distribution. This can further save both computational and analysis time by making it possible to focus on a specific feature of the flow (i.e., a specific POD mode) for a specific analysis, rather than the entire flow.
Conflict of Interest There is no conflict of interest to report.
References Algazi, V. and Sakrison, D., 1969. On the optimality of the Karhunen-Loève expansion (Corresp.). IEEE Transactions on Information Theory, 15(2), pp.319-321. Babiker, M.H., Gonzalez, L.F., Ryan, J., Albuquerque, F., Collins, D., Elvikis, A. and Frakes, D.H., 2012. Influence of stent configuration on cerebral aneurysm fluid dynamics. Journal of biomechanics, 45(3), pp.440-447. Baharoglu, M.I., Schirmer, C.M., Hoit, D.A., Gao, B.L. and Malek, A.M., 2010. Aneurysm inflow-angle as a discriminant for rupture in sidewall cerebral aneurysms: morphometric and computational fluid dynamic analysis. Stroke, 41(7), pp.1423-1430. Castro, M., Putman, C., Radaelli, A., Frangi, A. and Cebral, J., 2009. Hemodynamics and rupture of terminal cerebral aneurysms. Academic radiology, 16(10), pp.1201-1207. Cebral, J.R., Mut, F., Weir, J. and Putman, C., 2011. Quantitative characterization of the hemodynamic environment in ruptured and unruptured brain aneurysms. American Journal of Neuroradiology, 32(1), pp.145-151. Cecchi, E., Giglioli, C., Valente, S., Lazzeri, C., Gensini, G.F., Abbate, R. and Mannini, L., 2011. Role of hemodynamic shear stress in cardiovascular disease. Atherosclerosis, 214(2), pp.249-256. Chang, G.H., Schirmer, C.M. and Modarres-Sadeghi, Y., 2017. A reduced-order model for wall shear stress in abdominal aortic aneurysms by proper orthogonal decomposition. Journal of biomechanics, 54, pp.33-43. Chatterjee, A., 2000. An introduction to the proper orthogonal decomposition. Current science, pp.808817. Chung, B. and Cebral, J.R., 2015. CFD for evaluation and treatment planning of aneurysms: review of proposed clinical uses and their challenges. Annals of biomedical engineering, 43(1), pp.122-138.
15
Dolan, J.M., Meng, H., Sim, F.J. and Kolega, J., 2013. Differential gene expression by endothelial cells under positive and negative streamwise gradients of high wall shear stress. American Journal of Physiology-Cell Physiology, 305(8), pp.C854-C866. Ford, M.D., Alperin, N., Lee, S.H., Holdsworth, D.W. and Steinman, D.A., 2005. Characterization of volumetric flow rate waveforms in the normal internal carotid and vertebral arteries. Physiological measurement, 26(4), p.477. Fukunaga, K., 2013. Introduction to statistical pattern recognition. Elsevier. Gosling, R.G. and King, D.H., 1974. The Role of Measurement in Peripheral Vascular Surgery: Arterial Assessment by Doppler-shift Ultrasound. Grinberg, L., Yakhot, A. and Karniadakis, G.E., 2009. Analyzing transient turbulence in a stenosed carotid artery by proper orthogonal decomposition. Annals of biomedical engineering, 37(11), pp.2200-2217. Gunzburger, M.D., Peterson, J.S. and Shadid, J.N., 2007. Reduced-order modeling of time-dependent PDEs with multiple parameters in the boundary data. Computer methods in applied mechanics and engineering, 196(4-6), pp.1030-1047. He, X. and Ku, D.N., 1996. Pulsatile flow in the human left coronary artery bifurcation: average conditions. Journal of biomechanical engineering, 118(1), pp.74-82. Hodis, S., Uthamaraj, S., Smith, A.L., Dennis, K.D., Kallmes, D.F. and Dragomir-Daescu, D., 2012. Grid convergence errors in hemodynamic solution of patient-specific cerebral aneurysms. Journal of biomechanics, 45(16), pp.2907-2913. Hoi, Y., Meng, H., Woodward, S.H., Bendok, B.R., Hanel, R.A., Guterman, L.R. and Hopkins, L.N., 2004. Effects of arterial geometry on aneurysm growth: three-dimensional computational fluid dynamics study. Journal of neurosurgery, 101(4), pp.676-681. Hosain, M.L. and Fdhila, R.B., 2015. Literature review of accelerated CFD simulation methods towards online application. In Energy Procedia (Vol. 75, pp. 3307-3314). Kaspera, W., Ładziński, P., Larysz, P., Hebda, A., Ptaszkiewicz, K., Kopera, M. and Larysz, D., 2014. Morphological, hemodynamic, and clinical independent risk factors for anterior communicating artery aneurysms. Stroke, 45(10), pp.2906-2911. Lassila, T., Manzoni, A., Quarteroni, A. and Rozza, G., 2014. Model order reduction in fluid dynamics: challenges and perspectives. In Reduced Order Methods for modeling and computational reduction (pp. 235-273). Springer, Cham. Les, A.S., Shadden, S.C., Figueroa, C.A., Park, J.M., Tedesco, M.M., Herfkens, R.J., Dalman, R.L. and Taylor, C.A., 2010. Quantification of hemodynamics in abdominal aortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics. Annals of biomedical engineering, 38(4), pp.1288-1313. Lu, G., Huang, L., Zhang, X.L., Wang, S.Z., Hong, Y., Hu, Z. and Geng, D.Y., 2011. Influence of hemodynamic factors on rupture of intracranial aneurysms: patient-specific 3D mirror aneurysms model computational fluid dynamics simulation. American Journal of Neuroradiology, 32(7), pp.1255-1261. Lumley, J.L., 1967. The structure of inhomogeneous turbulent flows. Atmospheric turbulence and radio wave propagation. Meng, H., Tutino, V.M., Xiang, J. and Siddiqui, A., 2014. High WSS or low WSS? Complex interactions of hemodynamics with intracranial aneurysm initiation, growth, and rupture: toward a unifying hypothesis. American Journal of Neuroradiology, 35(7), pp.1254-1262. Mut, F., Löhner, R., Chien, A., Tateshima, S., Viñuela, F., Putman, C. and Cebral, J.R., 2011. Computational hemodynamics framework for the analysis of cerebral aneurysms. International journal for numerical methods in biomedical engineering, 27(6), pp.822-839. Pearson, K., 1901. LIII. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), pp.559-572.
16
Pinnau, R., 2008. Model reduction via proper orthogonal decomposition. In Model Order Reduction: Theory, Research Aspects and Applications (pp. 95-109). Springer, Berlin, Heidelberg. Quarteroni, A., Tuveri, M. and Veneziani, A., 2000. Computational vascular fluid dynamics: problems, models and methods. Computing and Visualization in Science, 2(4), pp.163-197. Rapún, M.L. and Vega, J.M., 2010. Reduced order models based on local POD plus Galerkin projection. Journal of Computational Physics, 229(8), pp.3046-3063. Rathinam, M. and Petzold, L.R., 2002. Dynamic iteration using reduced order models: a method for simulation of large scale modular systems. SIAM journal on Numerical Analysis, 40(4), pp.14461474. Rathinam, M. and Petzold, L.R., 2003. A new look at proper orthogonal decomposition. SIAM Journal on Numerical Analysis, 41(5), pp.1893-1925. Robertson, A.M. and Watton, P.N., 2012. Computational fluid dynamics in aneurysm research: critical reflections, future directions. Rosenfeld, A., 1976. Digital picture processing. Academic press. Sano, T., Ishida, F., Tsuji, M., Furukawa, K., Shimosaka, S. and Suzuki, H., 2017. Hemodynamic differences between ruptured and unruptured cerebral aneurysms simultaneously existing in the same location: 2 case reports and proposal of a novel parameter oscillatory velocity index. World neurosurgery, 98, pp.868-e5. Sirovich, L., 1987. Turbulence and the dynamics of coherent structures. I. Coherent structures. Quarterly of applied mathematics, 45(3), pp.561-571. Spiegel, M., Redel, T., Zhang, Y.J., Struffert, T., Hornegger, J., Grossman, R.G., Doerfler, A. and Karmonik, C., 2009, September. Tetrahedral and polyhedral mesh evaluation for cerebral hemodynamic simulation—a comparison. In Engineering in Medicine and Biology Society, 2009. EMBC 2009. Annual International Conference of the IEEE (pp. 2787-2790). IEEE. Spiegel, M., Redel, T., Zhang, Y.J., Struffert, T., Hornegger, J., Grossman, R.G., Doerfler, A. and Karmonik, C., 2011. Tetrahedral vs. polyhedral mesh size evaluation on flow velocity and wall shear stress for cerebral hemodynamic simulation. Computer methods in biomechanics and biomedical engineering, 14(01), pp.9-22. Sugiyama, S.I., Meng, H., Funamoto, K., Inoue, T., Fujimura, M., Nakayama, T., Omodaka, S., Shimizu, H., Takahashi, A. and Tominaga, T., 2012. Hemodynamic analysis of growing intracranial aneurysms arising from a posterior inferior cerebellar artery. World neurosurgery, 78(5), pp.462-468. Valen-Sendstad, K. and Steinman, D.A., 2014. Mind the gap: impact of computational fluid dynamics solution strategy on prediction of intracranial aneurysm hemodynamics and rupture status indicators. American Journal of Neuroradiology, 35(3), pp.536-543. van Disseldorp, E.M., Hobelman, K.H., Petterson, N.J., van de Vosse, F.N., van Sambeek, M.R. and Lopata, R.G., 2016. Influence of limited field-of-view on wall stress analysis in abdominal aortic aneurysms. Journal of biomechanics, 49(12), pp.2405-2412. Weatherill, N.P., 1992. Delaunay triangulation in computational fluid dynamics. Computers & Mathematics with Applications, 24(5-6), pp.129-150. Xiang, J., Natarajan, S.K., Tremmel, M., Ma, D., Mocco, J., Hopkins, L.N., Siddiqui, A.H., Levy, E.I. and Meng, H., 2011. Hemodynamic–morphologic discriminants for intracranial aneurysm rupture. Stroke, 42(1), pp.144-152. Xiang, J., Siddiqui, A.H. and Meng, H., 2014. The effect of inlet waveforms on computational hemodynamics of patient-specific intracranial aneurysms. Journal of biomechanics, 47(16), pp.38823890. Xiang, J., Tutino, V.M., Snyder, K.V. and Meng, H., 2014. CFD: computational fluid dynamics or confounding factor dissemination? The role of hemodynamics in intracranial aneurysm rupture risk assessment. American Journal of Neuroradiology, 35(10), pp.1849-1857. Zhou, G., Zhu, Y., Yin, Y., Su, M. and Li, M., 2017. Association of wall shear stress with intracranial aneurysm rupture: systematic review and meta-analysis. Scientific reports, 7(1), p.5331.
17
Figure 1. Patient-specific cerebral aneurysm model derived from computed tomographic angiography.
18
Figure 2. (a) Reference waveform from an averaged cerebrovascular waveform of 17 patients (Ford et al., 2005), and (b) six cardiac cycles of scaled waveforms with PI
[0.5,1.1].
19
Figure 3. (a) Eigenvalue versus the mode number, and (b) RIC, and (c) RMS error versus the number of modes used in the POD-based reconstruction for the patient-specific model used here.
20
Figure 4. Flow field for PI=0.7 at an arbitrary cross-section (left) and wall shear stress (WSS) (right) of the patientspecific aneurysm studied here based on (a,b) the ROM-based results, and (c,d) direct CFD simulations.
21
Figure 5. The error between the CFD-based and ROM-based results versus total number of modes used for PI=0.7.
22
Figure 6. Slices of homogenous flow velocity modes 1 – 3 of the ROM.
23