 Research
 Open access
 Published:
A comprehensive MRIbased computational model of blood flow in compliant aorta using radial basis function interpolation
BioMedical Engineering OnLine volume 23, Article number: 69 (2024)
Abstract
Background
Properly understanding the origin and progression of the thoracic aortic aneurysm (TAA) can help prevent its growth and rupture. For a better understanding of this pathogenesis, the aortic blood flow has to be studied and interpreted in great detail. We can obtain detailed aortic blood flow information using magnetic resonance imaging (MRI) based computational fluid dynamics (CFD) with a prescribed motion of the aortic wall.
Methods
We performed two different types of simulations—static (rigid wall) and dynamic (moving wall) for healthy control and a patient with a TAA. For the latter, we have developed a novel morphing approach based on the radial basis function (RBF) interpolation of the segmented 4Dflow MRI geometries at different time instants. Additionally, we have applied reconstructed 4Dflow MRI velocity profiles at the inlet with an automatic registration protocol.
Results
The simulated RBFbased movement of the aorta matched well with the original 4Dflow MRI geometries. The wall movement was most dominant in the ascending aorta, accompanied by the highest variation of the blood flow patterns. The resulting data indicated significant differences between the dynamic and static simulations, with a relative difference for the patient of 7.47±14.18% in timeaveraged wall shear stress and 15.97±43.32% in the oscillatory shear index (for the whole domain).
Conclusions
In conclusion, the RBFbased morphing approach proved to be numerically accurate and computationally efficient in capturing complex kinematics of the aorta, as validated by 4Dflow MRI. We recommend this approach for future use in MRIbased CFD simulations in broad population studies. Performing these would bring a better understanding of the onset and growth of TAA.
Background
Rupture of thoracic aortic aneurysm (TAA) is an acute medical condition, with a fatality rate of almost 95% [1]. Because of the high fatality, properly diagnosing and treating this dangerous condition is of utmost importance. However, the conventional guidelines that focus on the diameter and growth rate of TAA were shown to be inadequate in many cases [2, 3]. This emphasizes the need for new biomarkers that aim for patientspecific prediction of TAA rupture and look beyond the analysis based solely on the aorta geometry [4]. The blood flow information must be assessed to establish new predicting biomarkers. Such information can be obtained from 4D flow magnetic resonance imaging (MRI); however, its spatial and temporal resolution is limited [5]. In recent years, the clinical imagebased computational fluid dynamics (CFD) [6] was successfully applied to provide the patientspecific blood flow features in great detail, for example, flow in aorta [7,8,9] and TAA [10] as well as in wider population studies [11, 12]. However, one important aspect of modeling the aorta or TAA is often omitted in the literature—the movement of the aorta (i.e., aorta kinematic). Because of the beating heart during the cardiac cycle, the aortic root moves downwards during systole and returns to its original position during diastole. This movement was reported to be approximately nine millimeters in the downward direction [13] with a clockwise twist up to 20° [14]. Furthermore, the aortic compliance causes the wall to expand and contract radially during the cardiac cycle due to the changing transmural pressure gradient over time [15]. It was reported that changes in the thoracic aorta diameter were in the 1.7 to 3.6 mm range [16]. These combined effects of the aortic wall kinetics can significantly affect the blood flow simulations, and consequently, they should be included in the CFD simulation [17].
To model the blood vessel movement, two simulation strategies have been applied in previous studies in the literature: (i) the fully coupled fluid–structure interaction (FSI), and (ii) the predefined wall displacement. The FSI studies of aorta hemodynamics were applied in [17,18,19,20,21]. However, the FSI method for patientspecific situations suffers from numerous limitations. These include the lack of detailed information on the aortic wall properties (i.e., nonhomogeneous thickness and elasticity), difficulties with the physiological boundary conditions (for example, pressure), as well as the quite intensive computational costs (for example, iterative prestressing procedure, fluid/structure mechanics coupling). The estimation of the aorta motion was the focus of several studies in the literature [22,23,24]. A simplified method for the aortic wall motion was proposed in [24, 25]. The developed movingboundary method (MBM) tuned with the noninvasive clinical images (2D cineMRI) provided a good agreement with the FSI results [24] as well as with the measurements in terms of the luminal crosssectional area [25]. The MBM method was also less computationally expensive. However, while the MBM methods show excellent agreement with the measurements in terms of the change of luminal radius throughout the cardiac cycle, they cannot capture the rotational and/or longitudinal movement.
To overcome this, the meshmorphing approach based on radial basis function (RBF) was proposed for mimicking the motion of biological tissue, utilizing oneway coupling. Unlike the previously discussed FSI and MBM, the present method directly enforces the movement based on imaging data and therefore has the potential to closely mimic the complete movement of the arterial wall, while also being numerically efficient. While the choice of oneway coupling limits the method in the applications that consider the future progress of diseases, it can be an excellent tool for accurately simulating the present flow in arteries. RBF was successfully implemented in mimicking the motion of the aortic valve [26], left ventricle with mitral valve [27], and thoracic aorta [28, 29]. In the case of RBF application in the aorta, Capellini et al. [28, 29] presented an approach where only the ascending thoracic aorta (excluding root) was considered dynamic, and the rest of the domain was assumed to be rigid. Additionally, only a simplified inlet velocity boundary condition was implemented. These assumptions bring considerable simplifications to the complexity of motion and flow in the aorta.
To bridge these simplifications, we propose a proofofconcept approach for a 4Dflow MRIbased compliant model of the aorta. This approach will be evaluated for the healthy control subject and patientspecific aorta with a large root aneurysm. For both cases, the movement of the aorta and all inlet and outlet boundary conditions will be extracted from the corresponding 4Dflow MRI scans. The dynamic behavior of the aorta will be mimicked by a morphing approach utilizing the radial basis function (RBF) interpolation based on Xu and Kenjereš [27]. We adapt the method to account for the motion of the whole thoracic aorta. To define the motion, we utilize 4Dflow MRI data at several points of the cardiac cycle. In addition, we present an automatic registration protocol of the 4Dflow MRIderived velocity profile at the inlet for the moving aorta.
This article first introduces our main findings within the "Results" section, followed by the "Discussion" of these findings, and "Conclusions". Finally, the last section of our article focuses on a detailed explanation of the "Methods" that are utilized for this study.
Results
MRIbased wall movement
First, we want to assess the quality of the prescribed movement. The imagebased movement of the aortic wall for both studied subjects: the healthy control (HC) and the patient with a TAA (P), is shown in Fig. 1. To validate the results of the RBFbased interpolation, we compare geometries extracted from the 4Dflow MRI (blue isosurface) and RBFbased reconstruction (red isosurface)—both at the midacceleration time instant; Fig. 1a. Furthermore, we also compare characteristic circumferential wall profiles at various crosssections: (1) proximal ascending aorta (pAscAo); (2) distal ascending aorta (dAscAo), (3) proximal descending aorta (pDescAo), and (4) distal descending aorta (dDescAo), respectively. As can be seen in Fig. 1a, the RBF surface matches well the original 4Dflow MRI surface in the majority of the planes, except in the proximity of the root. Here, we can observe more variation. Additionally, in Fig. 1b, we show the timeevolution of the RBFbased circumferential profiles in identical crosssections (i.e., planes (1–4)) at the four keyframes (black—midacceleration, red—peak systole, blue—middeceleration, green—early diastole). These circumferential profiles visualize local change in the area during the dynamic simulations, which is most pronounced close to the aortic root. Finally, the profiles in Fig. 1 give us only a qualitative understanding of local variation between RBF and 4Dflow MRI surfaces. Therefore, to see the variation for the whole surface, we have calculated the absolute Euclidean distance between each point of RBFgenerated surfaces and the MRI segmentations. These data are shown for HC and P (at each keyframe) in Fig. 1c, and the median, mean, and standard deviation of the whole domain are reported in Table 1. Fig. 1c again showcases that the agreement between RBF and MRI surfaces is overall good, except in the proximity of the root. In addition, we can also observe an increasing variation in the agreement further from the peak systole.
To quantify the level of the aorta movement, the normalized timeaveraged aortic wall displacements (magnitude and corresponding coordinate directions) are shown in Fig. 2. The normalization was done using the radius of the inlet plane (i.e., \(r_{\text{HC}}^{\text{in}} = 1.49\cdot 10^{2}\) m, and \(r_\text{P}^\text{in} = 1.58\cdot 10^{2}\) m). The simulated displacement is lower for healthy control in comparison to the patient. In addition, we can observe a clear higher displacement in the ascending aorta for both studied cases.
Computational time
To compare the effect of the wall motion on the computational time for three simulated cycles, we report the wall time for one of the cases (HC). Both static and dynamic simulations (for this comparison) were run on 16 processors of AMD Opteron 6234. The reported wall time for static simulation was 33:07:59 (in [h:min:s]), and for the dynamic simulations, 38:29:31. The main differences in the computational time were related to the I/O intensive tasks (reading the prescribed mesh), rather than the simulations themselves.
Effect of wall movement on blood flow
Contours of the velocity magnitude at the preselected crosssections (dAscAo, pAscAo, and pDescAo) at four time instants of the cardiac cycle (midacceleration, peak systole, mid deceleration, and early diastole), for the healthy and patientspecific cases (both with static and dynamic simulations) are shown in Fig. 3. For MRI, the velocity field was reconstructed from 4Dflow MRI data extracted in planes in the flow direction. Note that for better visualization, used color maps are specifically adjusted for different crosssections and time steps. Overall, the computed profiles for HC resemble well the ones obtained by MRI, more differences can be observed for P, especially in planes 2 and 3. In addition, the results demonstrate a clear influence of movement on the calculated flow field, which is more significant for P.
Effect of wall movement on \(\mathrm {TAWSS}\) and \(\mathrm {OSI}\)
Timeaveraged wall shear stress (TAWSS) and oscillatory shear index (OSI) are two important flowderived quantities often mentioned in the literature as potential biomarkers to evaluate the progression of aortic aneurysms. The contours of the \(\mathrm {TAWSS}\) and \(\mathrm{ OSI}\) for the CFD_{static} and CFD_{dynamic} simulations are shown in Fig. 4. To make the comparison between the static and dynamic simulations easier, we also provided contours of the percentage differences \(\Delta \mathrm {TAWSS}\) and \(\Delta \mathrm {OSI}\). Note that for the dynamic simulation, contours of the TAWSS and OSI are shown for the peaksystole geometry. As can be seen for both subjects, we can observe slight differences between static and dynamic simulations for TAWSS, especially close to the root. For P, the differences in TAWSS are more pronounced in the whole ascending aorta. Dynamic simulations show significantly more differences for OSI; in both cases the differences in this quantity are more pronounced over the whole surface.
Next, we have also calculated the mean values of the absolute difference over the whole aorta surface (without side branches) for TAWSS and OSI. The mean value of absolute \(\Delta \mathrm {TAWSS}\) for the healthy control was \(\Delta \mathrm {TAWSS}_{\text {HC}}=2.72\pm 4.93\%\) Pa, and \(\Delta \mathrm {TAWSS}_{P}=7.47\pm 14.18\%\) Pa for the patientspecific geometry. For the OSI difference, the mean value (of absolute \(\Delta \mathrm {OSI}\) ) was \(\Delta \mathrm {OSI}_{\text{HC}}=12.87\pm 43.92\%\) for the former, and \(\Delta \mathrm {OSI}_{P}=15.97\pm 43.32\%\) for the latter.
Finally, while we can understand the spatial distribution of \(\Delta \mathrm {TAWSS}\) and \(\Delta \mathrm {OSI}\) based on the surface plots, they are unable to show the locality of the highest differences between static and dynamic simulations with respect to the range of \(\mathrm {TAWSS}\) and \(\mathrm {OSI}\). To overcome this, Fig. 5 shows the relationship between the dynamic TAWSS or OSI and the respective percentage difference (\(\Delta \mathrm {TAWSS}\) or \(\Delta \mathrm {OSI}\) [%]) for HC and P. In addition, we have plotted the average value of \(\Delta \phi\) for each of the assessed quantities and a binned average for the \(\mathrm {OSI}\). In the case of the binned average, the data were grouped based on a specific range of \(\mathrm {OSI}\) values (0.01) for the whole domain and average for the outliers with high \(\mathrm {OSI}\) (i.e., when the number of points within a range of OSI was less than 100).
Discussion
In the present work, we proposed an imagebased method for prescribing the motion of the aortic wall for CFD using RBF. The RBF method was chosen since it proved to be a viable approach to represent the complex motion of the aorta. By performing simulations with the predefined motion of the blood vessels, we avoid the necessity of obtaining detailed information regarding the vessel wall (i.e., elasticity and thickness). This information is usually not readily available, making the patientspecific studies challenging for the traditional FSI methods [18,19,20, 22,23,24]. Additionally, the computational time of simulations with prescribed motion is comparable with static simulations, as can be seen from our results. This is an important factor, especially considering the clinical applications with larger population studies. We have performed static and dynamic simulations for two geometries: the healthycontrol (HC) and patientspecific case (P) with a large TAA located close to the aortic root.
The prescribed RBFbased motion of the thoracic aorta matched well with the 4Dflow MRI, as illustrated in Fig. 1a for midacceleration. Some differences could be observed in the proximity of the aortic root. Close to the aortic root, the displacement is also the largest, as visualized in Figs. 1b and 2. In this region (ascending aorta), the motion is rather complex and is a result of the superposition of the axial and radial displacements. The longitudinal displacement (in the feet–head (FH) direction) generated by the physiological strain from the heart, makes an important contribution to the total displacement, as previously reported in several studies on aortic kinematics [13, 14, 30, 31]. It is important to note that this FHcomponent of aortic motion is not included in the FSI studies of the aorta [18,19,20,21], nor in the studies that model the predefined aorta motion based on its wall compliance [22,23,24], which can have a significant impact on the final results. In contrast to the ascending aorta, the descending aorta is less susceptible to motion due to the presence of the spinal column, and its displacement is predominantly in the radial direction [31], which can be also found in our results; Fig. 2.
The discrepancies between RBF and MRI in the proximity of the aortic root can also be found for the other keyframe geometries; Fig. 1c. The agreement between the original segmentation and RBF, in terms of distance between the surface vertices, is good for most of the investigated aortic domain (as seen in Table 1), except for the root. While this could (potentially) be avoided by locally increasing the density of the control points, the final moving geometry does not have to improve with respect to realistic aortic kinematics. This is caused by the accuracy of the original segmentation, which decreases for keyframes further from the peak systole [32], (Fig. 1c). The segmentation variability can be high, especially close to the root, as shown previously for healthy aortas [33]. Additionally, a similar argument can also be made for including more keyframe geometries. Currently, we only considered four geometries for the proofofconcept study. This choice was motivated by the segmentation procedure being very timeconsuming (approximately four hours per subject), with many manual adjustments necessary. Hence, including more details in the RBF procedure can, on the contrary, increase the amount of uncertainty and degrade our simulations.
Although we have only segmented a limited number of geometries to prescribe the movement, the selected keyframes probably contain the most radial expansion during the cycle. This allows us to capture most of the movement with the least amount of information required. Additionally, the ability of the proposed method to capture the movement based on only a limited number of segmentations is an important advantage for clinical use. In this case, the observer does not need to segment all the phases, saving valuable clinical practice time.
While the limitations behind segmentation of 4Dflow MRI are wellknown [32,33,34], our choice to use this method was motivated by the direct availability of measured flow data. This allows us to obtain accurate inlet and outlet boundary conditions and an estimation of the moving domain from a singular measurement. As shown previously by several studies [35,36,37,38] and in Appendix C, velocity profile from measurements should always be imposed as an inlet boundary condition, if available. Additionally, Gallo et al. [39] showed the importance of including as much patientspecific information as possible on all of the outlets of the studied domain to obtain accurate results. By utilizing 4Dflow MRI to obtain all of these, we can create a wellinformed, fully patientspecific model of the moving aorta without the necessity of additional measurements on the subjects.
To validate the proposed model, we can directly utilize the 4D flow MRI. As can be seen in Fig. 3, the computed profiles resemble well the ones acquired by MRI, especially in plane 1. More differences can be observed for the other planes, especially for lower velocity values (i.e., further from peak systole). These discrepancies could originate from the computational model, e.g., inaccurate inlet plane for boundary conditions [40] or the choice of outflow boundary conditions [41], which can have an effect up to \(5D_i\) upstream, what correlates to the locations of the examined planes. On the other hand, in 4Dflow MRI data, a higher noisetosignal ratio is present due to the static VENC, for these phases. This can cause limited velocity field acquisition [42] for these phases. Finally, 4Dflow MRI flow has been shown to be consistently underpredicted due to the temporal averaging of the data [43]. However, the simulations generally predict the blood flow behavior in the studied aortas.
Since an increasing number of studies are investigating the effect of aberrant blood flow on the development and rupture of aortic aneurysms, we highlight here the effects of the aorta motion on changes in blood flow patterns. Based on the presented results, certain differences in blood flow patterns were obtained with the static (a rigid wall assumption) and dynamic (predefined aorta motion) simulations; Fig. 3. The latter showed slightly better agreement, qualitatively, with the 4Dflow MRI measurements, especially during the decelerating part of the systole and early diastole. For example, in plane 2 and plane 3 for HC, during middeceleration and early diastole, the dynamic simulations were more accurate in capturing the velocity profile measured with 4D flow MRI. On the other hand, the differences between dynamic and static simulations for P are more striking, especially for plane 2 and plane 3, and we cannot make similar conclusions for this case. While these observations are only qualitative, they clearly show the effect of the aortic movement on the flow and the need to examine these differences in a larger number of patients and volunteers.
The flowderived variables such as the TAWSS and OSI are often used as potential biomarkers to indicate the onset and growth of aortic aneurysm [5]. Elevated regions of WSS were related to more rapid degradation of the extracellular matrix [44]. This causes weakening of the aortic wall (due to lack of elastin), and it was linked to the growth of the TAA [44]. On the other end of the spectrum, low TAWSS may lead to endothelial dysfunction, correlated with thickening of the aortic wall [9, 45] Furthermore, high OSI affects the response of the endothelial cells, and it was connected to the onset of atherosclerosis [46], a condition with high prevalence in patients with aortic aneurysms [47]. Our simulations revealed significant differences between TAWSS and OSI calculated from the static and dynamic simulations for both geometries; Fig. 4. In static simulations, the mean TAWSS is slightly overpredicted, especially in the aortic arch of P; Fig. 4. Nevertheless, the general trends in TAWSS distribution are similar for both of the simulation methods in each studied case. On the other hand, we can observe significantly more differences in OSI, which is consistently underpredicted by the static simulations, especially in the region of interest—ascending aorta; Fig. 4. These observations, including the maximal differences, are in accordance with other studies in the literature that considered the aorta movement—either by FSI [17] or by prescribed motion [24, 25, 29].
Additionally, we observed that aortic wall movement has a considerable effect on the whole range of TAWSS and OSI; Fig. 5. For OSI, it can be seen that the error introduced by the simulations is highest in the regions with low OSI. On the other hand, for higher values of OSI, which are physiologically important [46], the difference between static and dynamic simulation is lower, yet, still as high as 21% for HC and 37% for P. In the case of TAWSS, for which both low and high TAWSS are physiologically important [9, 44, 45], we could not observe such a clear correlation. Here, the difference between static and dynamic simulations is similarly dominant for the whole range of TAWSS (on average, up to 9% for HC and 15% for P, which are similar to other studies in literature [17, 24]). These findings highlight the importance of including aortic wall motion in the simulations, to prevent misinterpretation of the results.
We also found that the differences between static and dynamic CFD for both TAWSS and OSI are correlated with the magnitude of displacement; Fig. 2. The differences in both variables were most significant in the proximity of the aortic root and the ascending part of the aorta, i.e., in the regions where displacement is most prominent. In the arch and descending aorta, the differences between static and dynamic simulations and resulting TAWSS and OSI were smaller. Capellini et al. [29] presented an approach where only a portion of the aorta (ascending thoracic aorta) is considered moving, and the rest of the domain is static. For this case, they showed that the differences downstream of the moving region are negligible. On the contrary, as shown in the presented study, the movement in the arch and descending aorta still affects the flow considerably and should not be omitted. In conclusion, the aspects of the movement of the whole aorta should be included in the new generation of CFD simulations for accurate modeling of blood flow.
Finally, we need to contextualize our findings with respect to other possible sources of uncertainty in the simulations. As shown in our previous study, WSS is highly affected by the segmentation variability, with a local deviation of up to 50% (at peak systole) [33]. Similar or higher uncertainty as found for OSI and TAWSS in our results, was also reported due to inflow rates [48, 49], outflow boundary conditions [39], or inclusion of turbulence modeling [8, 50]. Nevertheless, due to a lack of data for dynamic simulations, specifically for simulations with prescribed motion, it is not possible to generalize whether rigid assumption for the aorta is sufficient in terms of uncertainty, unlike for other parts of the cardiovascular system [51]. For this, a more thorough followup study is necessary, including a larger number of pathologies.
Next, we address several limitations of the present work. To demonstrate the proofofconcept of the adopted RBFbased morphing approach in mimicking the aortic motion, we have considered two geometries: the healthy control and the patientspecific TAA. Future studies can include significantly larger numbers of both subject and patientspecific cases. Moreover, the patientspecific cases should include additional aortic pathologies such as dissection and coarctation [52]. Nevertheless, our work aimed to investigate the feasibility, accuracy, and numerical efficiency of the proposed method. Since we have selected an advanced stage of TAA as one of the test cases, it is expected that the method will also perform well for lessdeveloped pathologies. We also assumed that there was no aortic movement during the diastole. This assumption was a consequence of the unattainable segmentation of the 4Dflow MRI scans due to very low blood flow intensity during this period of the cardiac cycle. However, this assumption is valid since the aortic motion during diastole is limited [16], and we do not expect significant deviations from our findings. Finally, the presented simulation method with aortic motion was coupled with the 4Dflow MRI clinical data; here, we need to address two points: (1) the 4Dflow MRI acquisition is affected by several acquisition parameters such as efficient respiratory motion compensation, VENC, and Sense factor that reflects the amount of parallel imaging for acceleration. All of these can have an effect on the signaltonoise ratio and hence the segmented data. (2) The current segmentation procedure requires significant manual adjustments to properly capture the exact wall position at particular time instants of the cardiac cycle. We also addressed some of this segmentation variability on the calculated WSS in our previous study [33], where we observed significant variability in WSS due to the segmentation procedure. Since a similar protocol was also used in this study, this could also affect the prescribed wall movement. Additionally, using this technique hinders proper capturing of the aortic dilatation since the absolute difference between the root diameter of the systolic and diastolic phase can be lower than the resolution of 4Dflow MRI, as reported by De Heer et. al. [53]. However, here developed numerical simulation methodology can be directly integrated with other clinical imaging procedures as well (US, MRA, CT), which would improve the segmentation variability and the resolution to capture the motion properly.
Conclusions
In the present work, we showed how the aortic wall motion can be simulated by applying an efficient imagebased geometry morphing approach based on the radial basis function (RBF) interpolation. The simulated aortic motion was in good agreement with the 4Dflow MRI extracted geometries. The developed method proved to be accurate and numerically robust for both considered cases: the healthycontrol and the patientspecific aorta with an aneurysm in the aortic root. The computational time for dynamic simulations (with moving aortic walls) was similar to their static (with rigid wall assumption) counterparts, confirming the numerical efficiency of the proposed method. Effects of wall motion in the dynamic simulations were most prominent in the ascending aorta and this improved agreement with the 4Dflow MRI in comparison to the static simulations. We also report on the largest differences between the calculated TAWSS and OSI for static and dynamic simulations in the ascending part of the aorta. This shows the importance and necessity to include aortic wall motion in the CFD simulations in obtaining more accurate flow and flowderived biomarkers, such as the TAWSS and OSI. Based on here presented proofofconcept study on two geometries and improved agreement with the 4Dflow MRI, we propose to apply the presented moving wall approach on larger cohorts of patientspecific cases with various aortic pathologies.
Methods
Studied cases
Two subjects were included in this study—a healthy control (HC) and a patient (P). The patient had a root aortic aneurysm with a diameter \(D=50\) mm and aortic valve regurgitation of 33%. Additional characteristics for both subjects can be found in Table 2.
MRI acquisition and data processing
For both subjects, 4Dflow MRI was performed on a 3T MRI system (Elition, Philips Healthcare, Best, The Netherlands) using a hemidiaphragm respiratory navigator with retrospective electrocardiogram gating without echoplanar imaging. Additional parameters in the MRI sequence can be found in Table 3.
The acquired 4Dflow MRI data sets were segmented using CAAS MR Solutions v5.2. (Pie Medical Imaging BV, Maastricht, The Netherlands). The protocol for segmentation is identical for both studied subjects. The analysis is initialized by manually placing starting and ending points of the domain at peak systole. The starting point is placed in the aortic root, and the ending points are placed in all major branching arteries of the arch (brachiocephalic trunk, left common carotid artery, and left subclavian artery) and in the abdominal aorta. Subsequently, a 3D volume at peak systole is automatically segmented and manually adjusted (if discrepancies are observed). The manual adjustments for the peaksystolic phase are mostly necessary for the regions with flow recirculation, i.e., in the proximity of the aortic root and downstream the aortic arch. After successful segmentation, the peaksystolic 3D volume is copied to the next phase of interest and manually adjusted for the movement. In this case, the manual interventions are more complex and timeconsuming (up to three hours per phase) due to the arterial movement, both caused by the compliance of the aortic wall as well as the movement of the heart. This process is especially timeconsuming in the ascending aorta due to its complex movement through the cardiac cycle. The segmentation procedure is then repeated for all of the phases of interest (in total, four instants of the cardiac cycle were extracted—midrising systole (point 3 for HC and 2 for P), peak systole (point 6 for HC and 5 for P), middecreasing systole (point 9 for HC and 7 for P), and beginning of diastole (point 12 for HC and 10 for P)).
Geometry preprocessing
The initial surface obtained via segmentation of 4Dflow MRI is not suitable for CFD due to the relative ’roughness’ of the surface mesh (i.e., variation of the normal vector direction of the segmented surface from its ideal form due to segmentation errors) and inconsistent boundary faces of inlets and outlets. To remove these imperfections, we performed preprocessing of the extracted surfaces using Vascular Modelling Toolkit (VMTK) [54]. The initial (4Dflow MRI) and final surface after preprocessing for peak systole are shown in Fig. 6; note that while the whole aorta was extracted for HC, only the thoracic part was considered for the further simulations and analysis.
To obtain the final surface, the following steps are executed: first, the surface inlet and outlets are cut perpendicular to the arterial centerline. Next, the smoothing step is performed. In this step, the most optimal smoothing should account for the regions with high variation of the normal vectors while preserving the total volume. We have utilized the Taubin smoothing, with passband 0.1 and 100 iterations. Compared to other methods, the Taubin smoothing procedure ensures proper smoothing of regions with high variations in surface curvature and avoids extensive shrinkage of the surface [55]. Next, the surface mesh (triangular) is subdivided using the Butterfly method [56] to ensure better surface definition for the computational model. Finally, we added cylindrical extensions on the inlet and outlets in the normal direction of the respective planes. The diameter and the length of the extension are determined based on the diameter of the respective boundary (\(D_i\)). The length was kept constant for all outlets (\(5 \cdot D_i\)), and only a very short flow extension was created for the inlet (\(0.5\cdot D_i\)) to assure the reliability of the applied inlet velocity profile while ensuring the stability of the moving mesh implementation.
Computational model
The casespecific computational model was developed to take into account the detailed aorta geometry (and its movement), as well as the inlet and boundary conditions (BC) from the 4Dflow MRI scans. The entire algorithm is illustrated in the flowchart shown in Fig. 7. We have performed simulations with the rigid (static) and moving (dynamic) aortic wall for both subject and patientspecific geometries. For the latter, additional algorithm details are given in Fig. 7b and will be discussed below.
Fluid dynamics
In the present work, we adopt the ALE (Arbitrary Lagrangian–Eulerian) formulation for conservation of mass and momentum for a moving numerical mesh, for which the following governing equations are solved [57]:
where \(\rho\) is the fluid density, \({\textbf {v}}\) is the fluid velocity, \({\textbf {v}}_g\) is the grid (or mesh) velocity, p is the pressure, and \(\overline{\overline{\tau }}\) is the viscous stress tensor \(\left( \overline{\overline{\tau }}= \mu \left( \nabla {\textbf {v}} + \nabla {\textbf {v}}^T \right) \right)\), with \(\mu\) as the dynamic viscosity of the fluid. Note that for the static simulations (the rigid wall assumption), we have \({\textbf {v}}_g =0\). Additionally, since we employed a moving grid approach for part of our simulations, we need to define the space conservation law as follows:
where dV/dt is the volume derivative of the arbitrary control volume V, \(\partial V\) is the boundary of the arbitrary control volume V, \({\textbf {A}}\) is the face vector area, \(n_f\) is the number of faces j. Finally, the dot product on the righthand side is calculated from
where \(\delta V_j\) denotes the volume swept out by the control volume face j over each time step \(\Delta t\) [58].
Finally, we did not employ any turbulence model and assumed the flow to be laminar. This choice is justified since the mean Reynolds number (Re), was lower than the critical Re reported for aorta [59] for both of the cases (\(Re_{HC}=1890\), \(Re_{P}=1480\)).
Boundary and initial conditions
The inlet plane boundary condition was specified as a velocity inlet where all three velocity components at particular instants of the cardiac cycle were extracted from the reconstructed 4Dflow MRI (similarly to other studies [35,36,37,38]). This was done using an inhouse developed software tool for proper time registration and interpolation of the clinical data. All steps of this procedure are shown in the flowchart diagram shown in Fig. 7a, and can be summarized as:

1.
Using an inhouse developed tool for 4Dflow MRI data analysis, all three velocity components \({\textbf {v}} \left( v_x, v_y, v_z \right)\) are extracted from the reconstructed 4Dflow MRI data at the inlet plane of the studied case for each acquired time step (n = 34 for the healthy control, and n = 32 for the patientspecific acquisitions, respectively).

2.
Velocity data are linearly interpolated for each (n) and (n+1) time step, where timestep size (\(\Delta t\)) is based on the requirements of CFD (in the present work, we have \(\Delta t\) = 1 ms, for both cases).

3.
Inlet (represented by the CFD mesh) is imported from the base mesh (at peak systole), and the interpolated velocity profile is registered on this mesh; for dynamic simulations, the inlet is imported for each time step from the generated moving mesh.

4.
The velocity components are then interpreted in the CFD software as a Profile and interpolated and projected on the inlet mesh using inversedistance interpolation.
An example of the inlet flow rate and the interpolated velocity profiles at peak systole for the HC and P cases can be seen in Fig. 8.
Outlet boundary conditions were treated as outflow with a predefined fraction of mass flow per outlet. The outflow boundary condition assumes zerodiffusive flux for all flow variables and a mass balance correction at the outlet. The flow fractions at each outlet (\(w_N\)) were defined based on the 4Dflow MRI measurements. Since the 4Dflow data in the supraaortic arteries are unreliable due to a low number of voxels, we have exported planes in the upstream and downstream proximity of each bifurcation and, by that, estimated the net flow leaving through each outlet. In addition, due to the presence of bovine arch in P, we have extracted additional planes directly located in the outlets, downstream the brachiocephalic trunk, to account for the flow repartition. Afterward, the fractions at each time step were calculated as follows:
where \(w_n\) and \(Q_n\) are the flow fractions and the net flow of the respective outlets, \(Q_i\) is the net flow at the inlet and m is the total number of outlets. Finally, each outlet flow fraction is scaled by the sum of all of the fractions to satisfy \({ \sum _{n=1}^{m} Q_i\cdot w_N }=Q_i\). Then, the scaled fraction at each outlet (\(w_N\)) is defined as:
Applying the measured data at each timestep proved to be more accurate in the definition of the patientspecific simulations, as shown previously by Gallo et al. [39].
The noslip condition was applied at the wall for both static and dynamic simulations. The definition of wall movement for the dynamic simulations is discussed in detail in the next section. The transient simulations were initialized using the steadystate solution at the peak systole. In total, we have simulated three cardiac cycles to eliminate the influence of initial conditions. We have used only results from the last cardiac cycle for the final analysis.
Moving wall
In the present study, for dynamic simulations, we have adopted a predefined moving wall approach as shown in Fig. 7b. The wall motion was defined from four keyframe geometries extracted from the 4Dflow MRI (as previously described in MRI acquisition and segmentation section). The full process of the moving mesh generation over the entire cardiac cycle can be summarized as:

1.
The 4Dflowbased geometries at keyframes are preprocessed using VMTK (as described in Sec. ) yielding the initial surface of the aorta (in.stl format).

2.
For the geometry at the peak systole, various crosssection markers (planes) are introduced to separate the static (branching arteries) and dynamic (the rest of the aorta) segments.

3.
The numerical mesh is created for this aortic geometry (base mesh), with a refinement close to the wall.

4.
The control points are introduced for the peak systole and all keyframes by the following procedure:

(i)
Control points for the inlet and outlet are defined (circumferential equidistant distribution).

(ii)
A finite number of the planes perpendicular to the flow direction with uniform longitudinal distances are selected; in each of these planes, the radial distances are defined similarly to (i);

(iii)
Additional manually adjusted control points are introduced at locations in the proximity of the branching arteries.

(iv)
The final form of the structured control points matrices are established with \(i\times j\) control points (i = number of planes, j = control points per plane), for HC = 19 \(\times\) 6 and P = 18 \(\times\) 6.

(i)

5.
The base mesh (generated in step 3) is morphed using Radial Basis Function (RBF) interpolation of the control points (defined in the previous step), resulting in the morphed surface geometries for all selected keyframes.

6.
The keyframes surface geometries are then interpolated in time over the entire cardiac cycle using spline interpolation with smoothing parameter \(p=0.999\), resulting in a total of n = 1018 and 1212 frames. Note that we assumed no aortic movement during diastole.

7.
The generated surface geometries at each time step alongside the base mesh are then used as input for the RBFbased meshmorphing during the simulations.
Figure 9 depicts the surface points (both HC and P) for the reference phase (peak systole—red) and for one of the RBFgenerated frames (mid acceleration—blue) together with the control points for the two respective phases.
Physical and solver setup
The blood rheology was accounted for by applying the Carreau–Yasuda model:
where \(\mu _{\text{app}}\) is the apparent viscosity, \(\mu _{\infty }\) the viscosity at infinite shear, \(\mu _0\) the viscosity at zero shear, \(\lambda\) the relaxation time, \(\dot{\gamma }\) the shear rate, \(\alpha\) a shape parameter, and n the powerlaw index. The values for these parameters are adopted from [60], and are \(\mu _{\infty } = 2.2\) mPa\(\cdot\)s, \(\mu _0 = 22\) mPa s, \(\lambda = 0.110\) s, \(\alpha = 0.644\), and \(n = 0.392\). The blood density was kept constant (\(\rho = 1060\) kg/m\(^3\)).
The initial mesh was identical for static and dynamic simulations and consisted of tetrahedral elements with refinement close to the wall. We have performed a mesh dependency study for the peaksystolic flow conditions (all details shown in the Appendix). Based on the mesh dependency study, the final mesh consisted of \(n=1.58\cdot 10^6\) elements for the HC case and \(n=1.47\cdot 10^6\) elements for the P case. Specifically for the dynamic simulations, the smoothing and remeshing of the 3D mesh were conducted if element skewness was higher than 0.9. We have used the springbased smoothing with the spring constant factor of one and a maximum of 250 iterations allowed. For remeshing, the minimal and maximal allowed cell size for the whole domain varied between 1.76 \(\times 10^{4}\) m and 5.76 \(\times 10^{3}\) m.
The simulations were performed using Ansys Fluent 2019 R3 (Ansys, Canonsburg, Pennsylvania, USA). The main computational settings used in this study were: the pressurebased solver, PISO for pressure–velocity coupling, the secondorder upwind scheme used for the discretization of convective terms, the secondorder central differencing scheme (CDS) used for the discretization of diffusive terms, the time integration was performed by the secondorder fully implicit scheme, and the convergence criterion per time step of \(10^{5}\) was used for all quantities.
Postprocessing
The nearwall hemodynamic effects were studied by introducing several quantities averaged over the entire cardiac cycle:
where \(\mathrm {TAWSS}\) is the timeaveraged wall shear stress, T is the length of a cardiac cycle, and \(\overrightarrow{\tau _w}\) is the wall shear stress,
where \(\mathrm {OSI}\) is the oscillatory shear index. For the dynamic simulations, the values of TAWSS and OSI were projected and visualized on the surface geometry at the peak systole. Additionally, we have calculated the percentage difference (\(\Delta \phi\)) between CFD_{static} and CFD_{dynamic} for abovedefined quantities as:
where \(\phi _{\text{stat}}\) and \(\phi _{\text{dyn}}\) are the TAWSS or OSI for the static and dynamic CFD simulations, respectively.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.
References
Johansson G, Markström U, Swedenborg J. Ruptured thoracic aortic aneurysms: a study of incidence and mortality rates. J Vasc Surg. 1995;21(6):985–8. https://doi.org/10.1016/S07415214(95)70227X.
Pape LA, Tsai TT, Isselbacher EM, Oh JK, O’Gara PT, Evangelista A, Fattori R, Meinhardt G, Trimarchi S, Bossone E, Suzuki T, Cooper JV, Froehlich JB, Nienaber CA, Eagle KA. Aortic diameter ≥ 5.5 cm is not a good predictor of type a aortic dissection. Circulation. 2007;116(10):1120–7. https://doi.org/10.1161/CIRCULATIONAHA.107.702720.
Sweeting MJ, Thompson SG, Brown LC, Powell JT. Metaanalysis of individual patient data to examine factors affecting growth and rupture of small abdominal aortic aneurysms. Br J Surg. 2012;99(5):655–65. https://doi.org/10.1002/bjs.8707.
Pasta S, Agnese V, Gallo A, Cosentino F, Di Giuseppe M, Gentile G, Raffa GM, Maalouf JF, Michelena HI, Bellavia D, Conaldi PG, Pilato M. Shear stress and aortic strain associations with biomarkers of ascending thoracic aortic aneurysm. Ann Thorac Surg. 2020;110(5):1595–604. https://doi.org/10.1016/j.athoracsur.2020.03.017.
Cibis M, Potters WV, Gijsen FJ, Marquering H, Van Ooij P, vanBavel E, Wentzel JJ, Nederveen AJ. The effect of spatial and temporal resolution of cine phase contrast MRI on wall shear stress and oscillatory shear index assessment. PloS one. 2016;11(9):e0163316.
Kenjereš S. On recent progress in modelling and simulations of multiscale transfer of mass, momentum and particles in biomedical applications. Flow Turbul Combust. 2016;96:837–60. https://doi.org/10.1007/s1049401596692.
Tse KM, Chiu P, Lee HP, Ho P. Investigation of hemodynamics in the development of dissecting aneurysm within patientspecific dissecting aneurismal aortas using computational fluid dynamics (cfd) simulations. J Biomech. 2011;44(5):827–36. https://doi.org/10.1016/j.jbiomech.2010.12.014.
Perinajová R, Juffermans JF, Mercado JL, Aben JP, Ledoux L, Westenberg JJM, Lamb HJ, Kenjereš S. Assessment of turbulent blood flow and wall shear stress in aortic coarctation using imagebased simulations. Biomed Eng Online. 2021;20(1):84. https://doi.org/10.1186/s12938021009214.
Jamaleddin Mousavi S, Jayendiran R, Farzaneh S, Campisi S, Viallon M, Croisille P, Avril S. Coupling hemodynamics with mechanobiology in patientspecific computational models of ascending thoracic aortic aneurysms. Comput Methods Progr Biomed. 2021;205:106107. https://doi.org/10.1016/j.cmpb.2021.106107.
Jayendiran R, Campisi S, Viallon M, Croisille P, Avril S. Hemodynamics alteration in patientspecific dilated ascending thoracic aortas with tricuspid and bicuspid aortic valves. J Biomech. 2020;110:109954. https://doi.org/10.1016/j.jbiomech.2020.109954.
Zhu Y, Chen R, Juan YH, Li H, Wang J, Yu Z, Liu H. Clinical validation and assessment of aortic hemodynamics using computational fluid dynamics simulations from computed tomography angiography. Biomed Eng Online. 2018;17(1):53. https://doi.org/10.1186/s1293801804855.
Lu Q, Lin W, Zhang R, Chen R, Wei X, Li T, Du Z, Xie Z, Yu Z, Xie X, Liu H. Validation and diagnostic performance of a cfdbased noninvasive method for the diagnosis of aortic coarctation. Front Neuroinform. 2020;14:59. https://doi.org/10.3389/fninf.2020.613666.
Burggraf GW, Mathew MT, Parker JO. Aortic root motion determined by ultrasound: relation to cardiac performance in man. Cathet Cardiovasc Diagn. 1978;4(1):29–41. https://doi.org/10.1002/ccd.1810040104.
Stuber M, Scheidegger MB, Fischer SE, Nagel E, Steinemann F, Hess OM, Boesiger P. Alterations in the local myocardial motion pattern in patients suffering from pressure overload due to aortic stenosis. Circulation. 1999;100(4):361–8. https://doi.org/10.1161/01.CIR.100.4.361.
Nosek T. Essentials of human physiology. Gold Standard Multimedia Incorporated. 1998.
Guo JP, Jia X, Sai Z, Ge YY, Wang S, Guo W. Thoracic aorta dimension changes during systole and diastole: evaluation with ecggated computed tomography. Ann Vasc Surg. 2016;35:168–73. https://doi.org/10.1016/j.avsg.2016.01.050.
Alimohammadi M, Sherwood JM, Karimpour M, Agu O, Balabani S, DíazZuccarini V. Aortic dissection simulation models for clinical support: fluidstructure interaction vs. rigid wall models. Biomed Eng Online. 2015;14(1):1–16.
Reymond P, Crosetto P, Deparis S, Quarteroni A, Stergiopulos N. Physiological simulation of blood flow in the aorta: comparison of hemodynamic indices as predicted by 3D FSI, 3D rigid wall and 1D models. Med Eng Phys. 2013;35(6):784–91. https://doi.org/10.1016/j.medengphy.2012.08.009.
Qiao Y, Zeng Y, Ding Y, Fan J, Luo K, Zhu T. Numerical simulation of twophase nonNewtonian blood flow with fluidstructure interaction in aortic dissection. Comput Methods Biomech Biomed Eng. 2019;22(6):620–30. https://doi.org/10.1080/10255842.2019.1577398.
Pons R, Guala A, RodríguezPalomares JF, Cajas JC, DuxSantoy L, TeixidóTura G, Molins JJ, Vázquez M, Evangelista A, Martorell J. Fluidstructure interaction simulations outperform computational fluid dynamics in the description of thoracic aorta haemodynamics and in the differentiation of progressive dilation in marfan syndrome patients. Roy Soc Open Sci. 2020;7(2):191752. https://doi.org/10.1098/rsos.191752.
Vignali E, Gasparotti E, Celi S, Avril S. Fullycoupled fsi computational analyses in the ascending thoracic aorta using patientspecific conditions and anisotropic material properties. Front Physiol. 2021;12:732561. https://doi.org/10.3389/fphys.2021.732561.
Ganten MK, Weber TF, von TenggKobligk H, Böckler D, Stiller W, Geisbüsch P, Kauffmann GW, Delorme S, Bock M, Kauczor HU. Motion characterization of aortic wall and intimal flap by ecggated ct in patients with chronic bdissection. Eur J Radiol. 2009;72(1):146–53. https://doi.org/10.1016/j.ejrad.2008.06.024.
Bonfanti M, Balabani S, Greenwood JP, Puppala S, HomerVanniasinkam S, DíazZuccarini V. Computational tools for clinical support: a multiscale compliant model for haemodynamic simulations in an aortic dissection based on multimodal imaging data. J R Soc Interface. 2017;14(136):20170632. https://doi.org/10.1098/rsif.2017.0632.
Bonfanti M, Balabani S, Alimohammadi M, Agu O, HomerVanniasinkam S, DíazZuccarini V. A simplified method to account for wall motion in patientspecific blood flow simulations of aortic dissection: comparison with fluidstructure interaction. Med Eng Phys. 2018;58:72–9. https://doi.org/10.1016/j.medengphy.2018.04.014.
Stokes C, Bonfanti M, Li Z, Xiong J, Chen D, Balabani S, DíazZuccarini V. A novel mribased data fusion methodology for efficient, personalised, compliant simulations of aortic haemodynamics. J Biomech. 2021;129:110793. https://doi.org/10.1016/j.jbiomech.2021.110793.
Geronzi L, Gasparotti E, Capellini K, Cella U, Groth C, Porziani S, Chiappa A, Celi S, Biancolini ME. Advanced radial basis functions mesh morphing for high fidelity fluidstructure interaction with known movement of the walls: simulation of an aortic valve. In: Krzhizhanovskaya VV, Závodszky G, Lees MH, Dongarra JJ, Sloot PMA, Brissos S, Teixeira J, editors. Comput Sci ICCS 2020. Cham: Springer International Publishing; 2020. p. 280–93.
Xu F, Kenjereš S. Numerical simulations of flow patterns in the human left ventricle model with a novel dynamic mesh morphing approach based on radial basis function. Comput Biol Med. 2021;130:104184. https://doi.org/10.1016/j.compbiomed.2020.104184.
Capellini K, Vignali E, Costa E, Gasparotti E, Biancolini ME, Landini L, Positano V, Celi S. Computational fluid dynamic study for aTAA hemodynamics: an integrated imagebased and radial basis functions mesh morphing approach. J Biomech Eng. 2018;140(11):1110007. https://doi.org/10.1115/1.4040940.
Capellini K, Gasparotti E, Cella U, Costa E, Fanni BM, Groth C, Porziani S, Biancolini ME, Celi S. A novel formulation for the study of the ascending aortic fluid dynamics with in vivo data. Med Eng Phys. 2021;91:68–78. https://doi.org/10.1016/j.medengphy.2020.09.005.
Pagoulatou SZ, Ferraro M, Trachet B, Bikia V, Rovas G, Crowe LA, Vallée JP, Adamopoulos D, Stergiopulos N. The effect of the elongation of the proximal aorta on the estimation of the aortic wall distensibility. Biomech Model Mechanobiol. 2021;20(1):107–19. https://doi.org/10.1007/s1023702001371y.
Sailer AM, Wagemans BAJM, Das M, de Haan MW, Nelemans PJ, Wildberger JE, Schurink GWH. Quantification of respiratory movement of the aorta and side branches. J Endovasc Ther. 2015;22(6):905–11. https://doi.org/10.1177/1526602815605325.
Juffermans JF, Westenberg JJ, van den Boogaard PJ, Roest AA, van Assen HC, van der Palen RL, Lamb HJ. Reproducibility of aorta segmentation on 4d flow mri in healthy volunteers. J Magn Reson Imaging. 2021;53(4):1268–79. https://doi.org/10.1002/jmri.27431.
Perinajová R, Juffermans JF, Westenberg JJ, van der Palen RL, van den Boogaard PJ, Lamb HJ, Kenjereš S. Geometrically induced wall shear stress variability in cfdmri coupled simulations of blood flow in the thoracic aortas. Comput Biol Med. 2021;133:104385. https://doi.org/10.1016/j.compbiomed.2021.104385.
Corrado PA, Wentland AL, Starekova J, Dhyani A, Goss KN, Wieben O. Fully automated intracardiac 4d flow mri postprocessing using deep learning for biventricular segmentation. Eur Radiol. 2022;32(8):5669–78. https://doi.org/10.1007/s00330022086167.
Morbiducci U, Ponzini R, Gallo D, Bignardi C, Rizzo G. Inflow boundary conditions for imagebased computational hemodynamics: impact of idealized versus measured velocity profiles in the human aorta. J Biomech. 2013;46(1):102–9. https://doi.org/10.1016/j.jbiomech.2012.10.012.
Youssefi P, Gomez A, Arthurs C, Sharma R, Jahangiri M, Alberto Figueroa C. Impact of patientspecific inflow velocity profile on hemodynamics of the thoracic aorta. J Biomech Eng. 2017;140 (1). https://doi.org/10.1115/1.4037857.
Pirola S, Jarral OA, O’Regan DP, Asimakopoulos G, Anderson JR, Pepper JR, Athanasiou T, Xu XY. Computational study of aortic hemodynamics for patients with an abnormal aortic valve: the importance of secondary flow at the ascending aorta inlet. APL Bioeng. 2018;2(2):026101. https://doi.org/10.1063/1.5011960.
Armour CH, Guo B, Pirola S, Saitta S, Liu Y, Dong Z, Xu XY. The influence of inlet velocity profile on predicted flow in type b aortic dissection. Biomech Model Mechanobiol. 2021;20(2):481–90. https://doi.org/10.1007/s10237020013954.
Gallo D, De Santis G, Negri F, Tresoldi D, Ponzini R, Massai D, Deriu MA, Segers P, Verhegghe B, Rizzo G, Morbiducci U. On the use of in vivo measured flow rates as boundary conditions for imagebased hemodynamic models of the human aorta: Implications for indicators of abnormal flow. Ann Biomed Eng. 2012;40(3):729–41. https://doi.org/10.1007/s1043901104311.
Casciaro ME, Pascaner AF, Guilenea FN, Alcibar J, Gencer U, Soulat G, Mousseaux E, Craiem D. 4D flow MRI: impact of region of interest size, angulation and spatial resolution on aortic flow assessment. Physiol Meas. 2021;42(3).
Madhavan S, Kemmerling EMC. The effect of inlet and outlet boundary conditions in imagebased cfd modeling of aortic flow. Biomed Eng Online. 2018;17(1):66.
Ma LE, Markl M, Chow K, Vali A, Wu C, Schnell S. Efficient triplevenc phasecontrast mri for improved velocity dynamic range. Magn Reson Med. 2020;83(2):505–20. https://doi.org/10.1002/mrm.27943.
Markl M, Frydrychowicz A, Kozerke S, Hope M, Wieben O. 4d flow mri. J Magn Reson Imaging. 2012;36(5):1015–36. https://doi.org/10.1002/jmri.23632.
Guzzardi DG, Barker AJ, van Ooij P, Malaisrie SC, Puthumana JJ, Belke DD, Mewhort HE, Svystonyuk DA, Kang S, Verma S, Collins J, Carr J, Bonow RO, Markl M, Thomas JD, McCarthy PM, Fedak PW. Valverelated hemodynamics mediate human bicuspid aortopathy: insights from wall shear stress mapping. J Am Coll Cardiol. 2015;66(8):892–900. https://doi.org/10.1016/j.jacc.2015.06.1310.
Carroll JE, Colley ES, Thomas SD, Varcoe RL, Simmons A, Barber TJ. Tracking geometric and hemodynamic alterations of an arteriovenous fistula through patientspecific modelling. Comput Methods Progr Biomed. 2020;186:105203. https://doi.org/10.1016/j.cmpb.2019.105203.
Chakraborty A, Chakraborty S, Jala VR, Haribabu B, Sharp MK, Berson RE. Effects of biaxial oscillatory shear stress on endothelial cell proliferation and morphology. Biotechnol Bioeng. 2012;109(3):695–707. https://doi.org/10.1002/bit.24352.
Reed D, Reed C, Stemmermann G, Hayashi T. Are aortic aneurysms caused by atherosclerosis? Circulation. 1992;85(1):205–11. https://doi.org/10.1161/01.CIR.85.1.205.
Bozzi S, Morbiducci U, Gallo D, Ponzini R, Rizzo G, Bignardi C, Passoni G. Uncertainty propagation of phase contrastmri derived inlet boundary conditions in computational hemodynamics models of thoracic aorta. Comput Methods Biomech Biomed Engin. 2017;20(10):1104–12. https://doi.org/10.1080/10255842.2017.1334770.
Marzo A, Singh P, Larrabide I, Radaelli A, Coley S, Gwilliam M, Wilkinson ID, Lawford P, Reymond P, Patel U, Frangi A, Hose DR. Computational hemodynamics in cerebral aneurysms: the effects of modeled versus measured boundary conditions. Ann Biomed Eng. 2011;39(2):884–96. https://doi.org/10.1007/s104390100187z.
Manchester EL, Pirola S, Salmasi MY, O’Regan DP, Athanasiou T, Xu XY. Analysis of turbulence effects in a patientspecific aorta with aortic valve stenosis. Cardiovasc Eng Technol. 2021;12(4):438–53. https://doi.org/10.1007/s13239021005369.
Steinman DA, Pereira VM. How patient specific are patientspecific computational models of cerebral aneurysms? An overview of sources of error and variability. Neurosurgical Focus FOC. 2019;47(1):E14. https://doi.org/10.3171/2019.4.FOCUS19123.
Erbel R. Diseases of the thoracic aorta. Heart. 2001;86(2):227–34. https://doi.org/10.1136/hrt.86.2.227.
de Heer LM, Budde RPJ, Mali WPTM, de Vos AM, van Herwerden LA, Kluin J. Aortic root dimension changes during systole and diastole: evaluation with ecggated multidetector row computed tomography. Int J Cardiovasc Imaging. 2011;27(8):1195–204. https://doi.org/10.1007/s105540119838x.
Antiga L, Piccinelli M, Botti L, EneIordache B, Remuzzi A, Steinman DA. An imagebased modeling framework for patientspecific computational hemodynamics. Med Biol Eng Comput. 2008;46(11):1097. https://doi.org/10.1007/s1151700804201.
Taubin G. Curve and surface smoothing without shrinkage. In: Proceedings of IEEE International Conference on Computer Vision. 1995. pp. 852–7. https://doi.org/10.1109/ICCV.1995.466848.
Dyn N, Levine D, Gregory JA. A butterfly subdivision scheme for surface interpolation with tension control. ACM Trans Graph. 1990;9(2):160–9. https://doi.org/10.1145/78956.78958.
Khalafvand SS, Xu F, Westenberg J, Gijsen F, Kenjeres S. Intraventricular blood flow with a fully dynamic mitral valvle model. Comput Biol Med. 2019;104:197–204. https://doi.org/10.1016/j.compbiomed.2018.11.024.
Ansys I, editor. Ansys fluent theory guide. vol. 18.2, Ch. 4, 2017. pp. 39–91.
Stalder A, Russe M, Frydrychowicz A, Bock J, Hennig J, Markl M. Quantitative 2d and 3D phase contrast MRI: optimized analysis of blood flow and vessel wall parameters. Magn Reson Med. 2008;60(5):1218–31. https://doi.org/10.1002/mrm.21778.
Gijsen F, van de Vosse F, Janssen J. The influence of the nonNewtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model. J Biomech. 1999;32(6):601–8. https://doi.org/10.1016/S00219290(99)000159.
Buhmann M. Radial basis functions: theory and implementations, Cambridge monographs on applied and computational mathematics. Cambridge University Press. 2003. https://books.google.co.nz/books?id=TRMf53opzlsC
Biancolini M. Fast radial basis functions for engineering applications. Springer International Publishing. 2018. https://books.google.co.nz/books?id=8rRTDwAAQBAJ
Roache PJ, Ghia KN, White FM. Editorial policy statement on the control of numerical accuracy. J Fluids Eng. 1986;108(1):2. https://doi.org/10.1115/1.3242537.
Acknowledgements
Not applicable.
Funding
This research was financially supported by a grant from the Dutch Heart Foundation, The Netherlands (Grant Number CVON201808RADAR).
Author information
Authors and Affiliations
Contributions
R.P. developed and implemented the models, conducted the simulations, analyzed and interpreted the CFD results and MRI data, and was a major contributor to writing the manuscript, T.vdV. developed the models and was a contributor to writing this article, E.R. and F.X. developed the models, J.J. acquired and analyzed the MRI data, J.W. acquired and analyzed the MRI data and conceived the study, H.L. conceived the study, S.K. conceived the study, interpreted the results, and was a major contributor in writing the manuscript. All authors reviewed the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
This study protocol was approved by the Medical Ethics Committee (METC) of the Leiden University Medical Center (number P18.034 for healthy control and G20.149 for patient) and both subjects signed a consent to participate.
Consent for publication
Consent for publication has been obtained from both studied subjects.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix A: RBF interpolation—mathematical view
RBF interpolation is based on source points, i.e., the nodal points of the original geometry, and target points which are the nodal points of the deformed geometry. It assumes an unknown smooth function f that is given via the set of source and target points. This function is then approximated by an interpolant s(x) of which the general form is defined as [61, 62]:
where s(x) is the approximating smooth function, N the source points, \(\gamma _i\) the weights of the radial basis, \(\varphi\) the radial basis function, \(x_{k_i} = [x_{k_i}, y_{k_i}, z_{k_i}]\) the coordinates of the source points, and h(x) a polynomial part of which the degree depends on the type of radial basis function used. The polynomial part is added to guarantee the existence and uniqueness of the solution. The parameters \(\gamma _i\) and the polynomial coefficients are determined by solving a linear system of equations, with the order equal to N. This is defined by the passage and orthogonality condition [61, 62]:
where \(g_i\) are known discrete values of displacement of the source points \(x_{k_i}\), and q polynomials with a degree less than or equal to the degree of polynomial h. The first criterion ensures that s(x) goes through the given values \(g_{k_i}\). Secondly, the summation of the product of the weights and the polynomial at the source points \(x_{k_i}\) equals zero, which satisfies the orthogonality condition. The values of \(\gamma _i\) and polynomial coefficients \(\beta _i\) are found by solving [61, 62]
where \(M_{k_i}\) is the interpolation matrix containing the evaluation of the radial function based on the source points, and \(H_{k_i}\) the coordinate matrix in which the ith row is \(\big [{\begin{matrix} 1&x_{k_i}&y_{k_i}&z_{k_i} \end{matrix}}\big ]\). Finally, the displacement of the nonsource points \(s_m\) is calculated using [61, 62]:
where \(M_m\) is the evaluated matrix based on basis function \(M_{mij} = \varphi (x_{mj}x_{k_i})\), and \(H_m\) is the coordinate matrix with ith row \([1 \, x_{m_i} \, y_{m_i} \, z_{m_i} ]\).
We used the multiquadratics method for the radial basis function, which is defined as:
where r is the euclidean distance between the source (x) and nonsource (\(x_{k_i}\)) points (\(xx_{k_i}\)) and c is the shape parameter. The shape parameter was estimated based on the mean distance between the source points and their farthest neighbor normalized by the distance to their nearest neighbor. The estimated shape parameters for this study were \(c_{HC}=3.2\times 10^{3}\) for HC and \(c_{P}=4.0\times 10^{3}\) for P.
Appendix B: Mesh dependency
To investigate the mesh dependency of the simulations, we created three different meshes for HC (coarse—\(0.58\times 10^6\), medium—\(1.47\times 10^6\), and fine—\(4.35 \times 10^6\) control volumes) and performed simulations at the peaksystolic flow conditions. The mesh refinement ratio was \(r_{1,2} \approx r_{2,3} \approx 1.40\) for the fine/medium and medium/coarse meshes, respectively. Figure 10 shows WSS for all three meshes and data extracted alongside a line following the aorta for all three meshes.
We have also estimated the Grid Convergence Index (GCI) [63]. The results of the GCI analysis are given in Table 4. Finally, the estimated theoretical value of WSS at zero grid spacing using the Richardson extrapolation to be \(f_{h=0} = 7.79\) Pa.
Appendix C: Effect of inlet boundary conditions
To study the possible effects of various inlet boundary conditions (BC) on the flow and WSS distribution we have performed three simulations (for the static aortic wall) with the following velocity profiles in the inlet plane: the MRIbased, parabolic, and plug. The MRIbased inlet was defined using a reconstructed plane at the inlet from the 4Dflow MRI data. The uniform (plug) and parabolic velocity profiles were reconstructed such that their averaged velocity profiles give the corresponding MRIbased inlet flow rate. Using these settings, we have performed simulations for both the healthy control (Fig. 11) and the patientspecific (Fig. 12) geometries and obtained WSS distributions (and differences between the inlet BCs) are shown in Figs. 11 and 12. To make an easier distinction, the range of the WSS contours and particular WSS differences (the color map values) were adjusted per the considered case. Additionally, we have extracted two characteristic profiles along the ascending (line A) and descending (line B) parts of the aortic wall, to analyze the results in detail (Figs. 11 and 12).
Due to the complexity of the moving wall model, the inlet boundary conditions used in the literature are often simplified and defined as either plug or parabolic profiles. In the present work, we have demonstrated that the effects of the various specifications of the inlet velocity profiles were significant in the ascending aorta. In contrast, in descending part of the aorta, the impact of the various inlet velocity profiles was much less significant. Additionally, the quality of the segmentation in the proximity of the aortic inlet can have a significant impact on the calculated WSS distribution.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Perinajová, R., van de Ven, T., Roelse, E. et al. A comprehensive MRIbased computational model of blood flow in compliant aorta using radial basis function interpolation. BioMed Eng OnLine 23, 69 (2024). https://doi.org/10.1186/s1293802401251x
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1293802401251x