Computational Analysis of the Pulmonary Arteries in Congenital Heart Disease: A Review of the Methods and Results

With the help of computational fluid dynamics (CFD), hemodynamics of the pulmonary arteries (PA's) can be studied in detail and varying physiological circumstances and treatment options can be simulated. This offers the opportunity to improve the diagnostics and treatment of PA stenosis in biventricular congenital heart disease (CHD). The aim of this review was to evaluate the methods of computational studies for PA's in biventricular CHD and the level of validation of the numerical outcomes. A total of 34 original research papers were selected. The literature showed a great variety in the used methods for (re) construction of the geometry as well as definition of the boundary conditions and numerical setup. There were 10 different methods identified to define inlet boundary conditions and 17 for outlet boundary conditions. A total of nine papers verified their CFD outcomes by comparing results to clinical data or by an experimental mock loop. The diversity in used methods and the low level of validation of the outcomes result in uncertainties regarding the reliability of numerical studies. This limits the current clinical utility of CFD for the study of PA flow in CHD. Standardization and validation of the methods are therefore recommended.


Introduction
Stenosis of the pulmonary arteries (PA's) is commonly seen in patients with congenital heart disease (CHD). It can occur as a native substrate or after surgery [1]. Diagnosing patients with PA stenosis can be challenging as cardiac echography may be inconclusive. To confirm the diagnosis, often evaluation with multiple imaging modalities as CTA and CMR is necessary. The decision whether to treat the stenosis is primarily based on pressure gradients which need to be confirmed by cardiac catheterization [2]. However, pressure gradients might resolve under anesthesia resulting in possible under treatment of PA stenosis. In addition, restenosis and intima proliferation can occur after treatment. The mechanisms causing this are still not well understood, but several flow characteristics-i.e., turbulence, wall shear stress (WSS), and the interaction of flow and (stent) geometry-are hypothesized to be of influence [3][4][5]. Computational fluid dynamics (CFD) provides the opportunity to study these factors and enhance our knowledge of hemodynamics in the PA's. It allows for detailed flow visualization and simulation of exercise and treatment outcomes. However, its clinical use is still limited as there is a great variety in used methods and validation of the numerical outcomes is often challenging [6][7][8][9][10][11]. The aim of this review was to analyze the available literature on numerical studies of the pulmonary arteries in biventricular CHD, focusing on the used methods and validation of the results.
Pediatric and adult as well as animal studies were considered. Papers with the focus on pulmonary hypertension were excluded as well as papers from before 2001, reviews, and if no full text was available. All inclusion and exclusion criteria are shown in Figure 1.

Study Selection.
The search in the PubMed and Embase databases resulted in a total of, respectively, 266 and 286 papers. After undoubling and title abstract screening, 45 papers remained. The full-text screening of these papers revealed exclusion criteria for 11 articles. This resulted in a total of 34 eligible papers for this review. The flowchart for the selection process is shown in Figure 1.

Study
Overview. The 34 selected papers included four animal studies and 25 papers presenting human cases. In one paper, both human and animal cases were described [12]. In six studies, PA's were represented by straight or curved tubes [13][14][15][16][17][18]. Altogether, the articles presented a total of 256 geometries based on 126 subjects (Table 1).
In the following sections, the various used methods for numerical analysis of the PA's will be compared. The sections are subdivided into the major steps necessary for solving a CFD case. First is the anatomic reconstruction, then the meshing and setting the boundary conditions followed by the numerical setup and finishing with postprocessing and validation of the results. Table 2 provides a summary of the different strategies used in the selected papers.

Computational and Mathematical Methods in Medicine
A mesh independence test to evaluate mesh quality was performed in 16 of the 34 papers [14-16, 18, 19, 23, 24, 31, 33, 38-44]. In 12 of these, information on the criteria for mesh independency was provided [14,18,19,23,24,31,33,[38][39][40][41]44]. These criteria included velocity profiles on different locations and a difference of <5% in calculations of pressure, velocity, or WSS. In 18 of the 34 papers, a mesh independence test was not performed or at least not mentioned. Three studies applied an element size determined by a mesh independence test performed on a different geometry [17,21,30]. The number of elements for the final mesh varied between 30.000 and 4 million.

Boundary Conditions
3.5.1. Inlet. Table 3 shows all the used inlet and outlet boundary conditions in the included papers. The inlet boundary conditions were pulsatile in 23 and constant in 11 papers. The most common inlet boundary condition was flow rate (l/min) followed by velocity and pressure in 18, 10, and four studies, respectively. In two studies, an electrical system was applied at the inlet [14,26]. The inlet conditions were patient-specific in seven papers [12,16,27,[39][40][41]43]. In these studies, the stroke volume as well as the waveform was patient-specific. These conditions were mostly obtained by MRI or invasive measurements during cardiac catheterization. In five studies, a patient-specific stroke volume was implemented but with a general waveform [25,30,32,37,45]. This waveform was scaled to a cardiac index suitable for the analyzed geometry. In the other 23 studies, a general inlet boundary condition was used.
In 13 papers, the applied velocity profile was specified. This was a flat or plug flow velocity profile in six and a parabolic profile in three studies [12, 16, 17, 19, 22-24, 28, 34]. Three articles implemented a Womersley flow at the inlet, and one study used a patient-specific velocity profile obtained by phase-contrast MRI [32,[39][40][41]. However, most studies did not mention the kind of velocity profile they used on their inlet.
3.5.2. Outlet. The most applied outlet boundary condition was a constant pressure outlet. 14 studies used a variation of this condition, i.e., atmospheric pressure, zero pressure, or the mean PA or aortic pressure obtained by cardiac catheterization [

Computational and Mathematical Methods in Medicine
pulsatile pressure on the LPA outlet was implemented while a Womersley velocity profile was set on the RPA outlet [40]. In another study, the outflow boundary condition was defined by the LPA : RPA flow split [42]. In five papers, the threeelement Windkessel model was imposed to the outlet [32-34, 37, 45]. The pure resistance strategy was used in five studies [23,24,31,41,44]. Spilker et al. described a method for the impedance boundary condition in which they reconstructed a one-dimensional (1D) anatomy and calculated the impedance value for the pulmonary anatomy [12]. One study applied multiple outlet boundary conditions and compared results. These included zero pressure, constant pressure, prescribed flow split, and a lumped parameter model [17]. In two papers, the outlet boundary conditions were not specified [13,38].
Patient-specific information was used to calculate resistance and Windkessel values in four papers [16,27,31,32,37,40,41,45]. For this, i.e., the flow split derived from flow perfusion scans, pressure from catheterization, and cardiac output from catheterization or echocardiography were taken. The other authors estimated values based on more general information.
3.6. Vessel Wall Compliance. In nine of the included papers, FSI was used to simulate deformation of the vessel wall during the cardiac cycle [12,14,15,19,20,32,35,37,45]. One of these groups applied case-specific compliance [20]. Here, the mesh was subdivided in five regions. Young's modulus for each region was obtained by stretch testing the tissue of freshly harvested porcine pulmonary roots. These values were then imposed to the in silico geometry of the porcine pulmonary roots. No studies were available using patientspecific compliance in human cases. The eight other groups assumed a global and constant value for the compliance of the artery wall. The highest Young modulus was 5 * 10e7 Pa, and the lowest 2:6 * 10e5 [19,32]. One group implemented a Young modulus varying between 2.6 and 4:2 * 10e5 Pa [37]. They tuned the value until the computed outcomes matched the desired patient-specific outcomes. The wall thickness was assumed to be between 0.5 and 1.5 mm. In one article, a variable thickness of 10% of the diameter of the vessel was applied [14]. Four articles specified the Poisson ratio used. These were, respectively, 0.42, 0.45, 0.49, and 0.5 [14,19,35,37].

Fluid Characteristics.
In 29 of the 34 included papers, blood was assumed to behave as a Newtonian fluid while in one paper, it was assumed to be a non-Newtonian fluid [40]. Three papers did not specify the assumption they made [13,34,37]. In most studies, the blood density was set to be 1060 kg/m 3 . Only three studies assigned a different density of 1050 kg/m 3 and 1000 kg/m 3 , respectively [33,35,40]. In nine papers, the used density was not described. Viscosity was mostly assumed to be 0.004 kg/ms [12, 14, 15, 17, 21, 22, 28-30, 33, 39, 41, 42, 45]. Other imposed values were 0.0035 kg/ms, 0.003 kg/ms, and 0.00371 kg/ms [25,26,35,38,43,44]. Four authors applied a varying viscosity number [16,18,36,40]. One of them varied the viscosity of blood between 0.003 and 0.008 kg/ms depending on hematocrit levels varying between 30 and 55% [16]. Two authors used the Carreau model to capture the varying viscosity of blood depending on the shear rate [18,40]. The last one analyzed the stability of their algorithm with varying viscosity numbers [36].
3.8. Computational Time. The computational time for the cases was described in six papers. The reported time per case varied between a couple of hours up to 1-2 months [21,23,35,38,39,41]. In one paper, the difference in computational time for the same case with a different number of cores was shown. With the use of "supercomputers," computational time was reduced to a couple of hours for highly complex cases [35].
3.9. Results and Validation. The results presented in the included papers varied according to the research question proposed. In most papers, two or three numerical outcomes were presented, i.e., pressure and WSS or streamlines and flow rates. In the majority of the papers, figures showed the results of the peak systolic and one or two diastolic time steps. Velocity and WSS were the results most reported followed by, respectively, pressure, streamlines, and flow rates. Pressure outcomes were presented either by the peak systolic numbers and energy loss over a stenosis or as the time-pressure curve over the cardiac cycle.
In nine papers, clinical data of the included patients was used to validate the CFD results (Table 4) [12,19,27,32,33,37,39,40,45]. In seven of these papers, one hemodynamic outcome was validated [12,19,27,32,33,37,39,45]. These were flow results in four, pressure in two papers, and regurgitation fraction in one paper. One paper verified pressure as well as flow rate results [40]. In all the papers validating computational pressure outcomes, cardiac catheterization data was used as the golden standard. In three papers, the absolute numbers of diastolic and systolic pressure outcomes were presented and validated [12,27,32]. The invasively measured pressure curve was compared to the computational pressure curve in one study [40]. Other sources for validation were cardiac MRI or lung perfusion scans. With data from these sources, flow rate, flow split, and regurgitation fraction outcomes were validated. In one paper, the source of validation of the flow rates was not specified [33]. One paper verified wall deformation results of a non-patient-specific case using an experimental mock-loop setup [19]. In the 26 other papers, there was no comparison between CFD outcomes and clinical data.

Discussion
The use of advanced imaging modalities to describe the hemodynamic impact of PA stenosis is increasing. CFD is one of these techniques providing detailed visualization of patient-specific hemodynamics. The aim of this paper was to review the numerical methods and-clinical-validation of CFD for evaluation of PA's in biventricular CHD. All of the papers included in this review emphasize the importance of hemodynamic evaluation of the PA's in biventricular CHD. They show the use and feasibility of CFD for this purpose and the wide variety of applications of the technique, i.e., for surgical or interventional treatment planning, research on complications, and exercise simulation. However, this review also shows limitations of the current available literature.
The literature reveals a large diversity in the setup for the numerical analysis of PA stenosis. This heterogeneity is important as variations in the numerical case setup significantly influence the outcomes. Results of patient-specific analysis are highly dependable on the source and quality of the anatomic reconstruction. In addition, small differences in the applied inlet or outlet boundary conditions can have a major impact. Outcomes of WSS and velocity can differ up to 30% with different boundary conditions [46][47][48][49]. In the included papers for this review, 10 different sources for anatomic reconstruction were used and 10 different methods were identified for the definition of the inlet boundary conditions. The largest variety however is seen in the definition of outlet boundary conditions. The 34 papers described 17 different approaches to assign outlet boundary conditions with very limited use of complete patient-specific boundary conditions. In the majority of cases, assumptions or generalizations defined inlet and or outlet boundary conditions. In several papers, key methodological information was missing. This included missing information on mesh size (23/34), number of cardiac cycles simulated (15/34), and convergence criteria (21/34).
The heterogeneity, assumptions, and generalizations in the computational setup result in uncertainties regarding the outcomes. The validation of methods and results is therefore of major importance. It provides direct feedback on the used methods and increases confidence in the reliability of the technique. This review shows that the level of validation of the CFD outcomes is very low. Studies with the main aim to validate CFD outcomes were completely missing, and only nine of the 31 papers compared their outcomes to clinical data. This lowers the translational value of the studies.
Another important limitation for the clinical utility of CFD is the computational time. This was reported to be as long as several days to even months per case. However, more and more progress is made in speeding up the computational process. By use of "super computers," improved algorithms, and cloud-based environments, the simulation time can be significantly reduced. Great examples of these efforts are shown by the two papers of Kong et al. included in this review. They show how the use of multiple cores and adjustment of algorithms can decrease computational time with several hours [35,36]. This can be expected to further decrease in the coming years.

Conclusion
The aim of this review was to evaluate the available literature on numerical analysis of the PA's in biventricular CHD. The 8 Computational and Mathematical Methods in Medicine focus was on the used methods and the rate of-clinical-validation of the outcomes. To the best of our knowledge, this is the first review evaluating the different strategies for numerical studies of the PA's. The included literature shows the wide variety of applications of CFD and emphasizes the added value of numerical studies for hemodynamic assessment of the PA's. However, this review also shows the large heterogeneity in used methods in all parts of the numerical setup and little validation of the results. This limits the current clinical utility of CFD. To increase the translation towards clinical use, standardization of the methodologies is desirable. Future research should therefore be pointed towards validation of methods of numerical studies.