Computational Fluid Dynamics Support for Fontan Planning in Minutes, Not Hours: The Next Step in Clinical Pre-Interventional Simulations

Computational fluid dynamics (CFD) modeling may aid in planning of invasive interventions in Fontan patients. Clinical application of current CFD techniques is however limited by complexity and long computation times. Therefore, we validated a “lean” CFD method to magnetic resonance imaging (MRI) and an “established” CFD method, ultimately aiming to reduce complexity to enable predictive CFD during ongoing interventions. Fifteen Fontan patients underwent MRI for CFD modeling. The differences between lean and established approach, in hepatic and total flow percentage to the left pulmonary artery (%LPA), power loss and relative wall shear stress area were 1.5 ± 4.0%, -0.17 ± 1.1%, -0.055 ± 0.092 mW and 1.1 ± 1.4%. Compared with MRI, the lean and established method showed a bias in %LPA of -1.9 ± 3.4% and -1.8 ± 3.1%. Computation time was for the lean and established approach 3.0 ± 2.0 min and 7.0 ± 3.4 h, respectively. We conclude that the proposed lean method provides fast and reliable results for future CFD support during interventions. Graphical abstract


Introduction
Congenital heart disease (CHD) is the most common class of major congenital malformations and is the leading cause of mortality from birth defects [1]. The reported prevalence of CHD in the general population varies between 8 and 10 per 1000 live births [2]. Modern surgical and medical care have significantly extended life expectancy and this leads to increasing number of adult CHD patients, thereby increasing the needs for health resources in this patient group [1,3].
Single ventricle (SV) includes a wide range of complex univentricular CHD where children typically require 2 to 3 surgical palliations, aiming to achieve lung perfusion via connection between the systemic caval veins and the pulmonary artery branches, called total cavopulmonary connection (TCPC) or "Fontan" circulation.
The first palliation is often performed neonatally and consists of either surgical or transcatheter shunt between a major systemic central vessel or the right ventricle and the pulmonary artery. At the second stage ("bidirectional Glenn"), the superior vena cava (SVC) is surgically anastomosed to the right pulmonary artery (RPA). At the third stage ("TCPC"), blood in the inferior vena cava (IVC) is directed into the pulmonary arteries via an extracardiac GoreTex™ conduit [4,5]. While early outcomes are generally acceptable, severe adverse events may occur in the long-term including proteinlosing enteropathy, declining functional status, increased Associate Editor Jozine ter Maaten oversaw the review of this article Petter Frieberg and Nicolas Aristokleous contributed equally to this work pulmonary vascular resistance, exercise intolerance, heart failure, and hepatic and renal dysfunction [6,7]. These complications may lead to redo surgery or transcatheter interventions. These procedures are often complex and would therefore greatly benefit from predictive computer simulations of treatment outcome. Ideally simulations should be used earlier to predict the results of Glenn to TCPC surgery and guide the initial surgery to optimize the Fontan circulation.
Currently, simulations using computational fluid dynamics (CFD) offer the possibility to model the Glenn and TCPC circulation and can provide measures of flow and how these would be affected by interventions that change the geometry. Furthermore, advanced hemodynamic measures such as power loss ( Ė loss ) [8] and wall shear stress (WSS) [9] can be obtained from CFD and may further the understanding of pathophysiology in the Fontan circulation.
However, clinical use of CFD is still limited. A major reason is the complex methods, which involves the use of advanced software, powerful computers, and interaction between clinicians, medical physicists, and engineers [10]. Other reasons found are sources of uncertainty that accumulate during the modeling process from image artifacts, noise, and inherent limitations of the image acquisition, to simulation assumptions and simplifications [11]. Trusty et al. provided a detailed description of this time-consuming procedure, with 60 h required for simulation of one surgical option at one physiologic condition [10].
Our goal is to have a clinically integrated and lean CFD framework by which clinicians and engineers can interactively perform interventional planning and get CFD results on pre-prepared models, all within minutes. Ultimately, we aim to realize a framework that could enable near real-time predictive CFD support adjacent to the operating room during ongoing interventions.
We have previously proposed a simplified "lean" method to reduce the complexity in this process so that more centers can adopt CFD in clinical routine [12], but it has not been validated and compared with established CFD approaches. One further simplification to reduce computation time is to assume steady flow and ignore the pulsatile nature of venous flows [13], but how this simplification affects lean CFD results has not been previously examined.
Therefore, we aimed to validate both a clinically feasible lean CFD approach and an established CFD approach using identical boundary conditions, against flow measurements from magnetic resonance imaging (MRI). Next, we aimed to compare the lean and established approach for additional hemodynamic measurements such as hepatic flow distribution, wall shear stress, and power loss. Furthermore, we aimed to verify previous findings [13] on whether steady (time-averaged) flow is an acceptable simplification compared to pulsatile flow computed using both the lean and the established CFD solver. Finally, we aimed to record the required effort in user-and computation time to produce reliable CFD results using the proposed methods.

Study Population and MR Imaging
Patients with SV (n = 15, median age 6.7, range 2.3 to 17 years, 6 females) with either bidirectional Glenn (n = 6) or TCPC (n = 9) were included in the study. The Glenn patients were prospectively included in a research protocol and imaged at a 1.5 T Siemens Aera scanner (Siemens Healthineers, Erlangen, Germany). The TCPC patients were retrospectively included in a previous study and imaged at a 1.5 T Philips Achieva scanner (Philips Healthcare, Best, The Netherlands) [12]. Cine images were acquired using a steady-state free precession (SSFP) sequence (typical parameters TR/TE/flip angle: 2.9 ms/1.5 ms/60°, slice thickness 5 mm, in-plane resolution 1.2 mm × 1.2 mm). Two-dimensional phase contrast MRI (2D PC-MRI) flow measurements were acquired using a velocity encoded fast field echo sequence (typical parameters TR/TE/flip angle: 10 ms/6.5 ms/15° in-plane resolution 1.2 mm × 1.2 mm). Flow measurements of the SVC, the IVC/TCPC conduit, the pulmonary artery branches, and pulmonary veins were typically acquired using velocity encoding 80 cm•s −1 and for the aorta 200 cm•s −1 . Image analysis was done using the freely available software Segment (Medviso AB, Lund, Sweden, http:// segme nt. heibe rg. se) [14]. The study was approved by the Regional Ethical Review Board in Lund, Sweden. Adult patients gave written informed consent. For patients under 18 years of age, written informed consent was given by their parents. Patient characteristics are shown in Table 1.

Surface Model Reconstruction
Patient-specific 3D models of the proximal Glenn/TCPC anastomosis were constructed in Creo Parametric (PTC, Boston, MA, USA). The geometry was created by importing MRI segmentation curves created in Segment into Creo Parametric, where geometry was created on the curve boundaries. Examples of TCPC and Glenn models are shown in Fig. 1.

Numerical Approach
In this study, two different commercial CFD software packages and patient specific results from MRI were used to obtain hemodynamic and time resource parameters. The following comparisons were made under identical boundary conditions in CFD: 709 a) Pulmonary flow distribution to the left pulmonary artery obtained from the lean and established CFD were compared with results from MRI as the reference standard. b) Pulmonary-and hepatic flow distribution to the left pulmonary artery, power loss, and normalized wall shear stress obtained from the lean and established CFD were compared under steady and pulsatile inlet boundary conditions. c) Time required for meshing, simulation set-up, and solving were compared for the lean and established CFD.

Computational Fluid Dynamics
The numerical solution was obtained using two commercial CFD packages. For the lean CFD, we used a solver Simcenter FloEFD for Creo (Siemens EDA, Wilsonville, OR, USA). FloEFD is embedded in the user interface of Creo Parametric and uses the immersed boundary method which is well known and has been previously used in Fontan CFD simulations [15,16] As established CFD we used as solver STAR-CCM + (v2019.1, Siemens PLM Software, Plano, TX, USA), whose numerical core is based on the finite volume method (FVM). Each new patient set-up was based on a typical "Fontan template" with pre-populated fields with non-patient-specific values for inlet flows, outlet boundaries, measurement of hemodynamic parameters, and rheology in the software. These are out-of-the-box capabilities that we found was the most time-efficient approach in both methods.
All simulations and setup thereof were performed by the same user who is an expert operator well acquainted with both softwares.
b) rigid vessel walls c) laminar flow d) patient specific mass flow rates from MRI were applied at each inlet e) zero pressure at the outlets representing the common pressure of the atrium The models included a proxy of pulmonary vascular resistance (PVR) based on MRI flows. This allows the pulmonary artery flow to change following future anatomic interventions, assuming PVR will not change in the short term. The used model of resistance satisfies the fundamental equation: where Δp is the transpulmonary gradient, R is the pulmonary vascular resistance, and F is the pulmonary artery flow. Linear porosity satisfies this equation and is ubiquitous among commercial CFD packages. Porosity was used to model constant values of R, assigned to porous baffles near the outlets (Fig. 1). Single ventricle patients often present with aortopulmonary collaterals (APC), which can be quantified with MRI as the difference between pulmonary vein flow and pulmonary artery flow [18]Compared to the proximal pulmonary artery flow, APC flow increases the total pulmonary flow F with a factor k: which can be obtained from MRI for each patient and for each lung. Instead of physically including APC flow in the simulations, the constant pulmonary vascular resistance R was multiplied with patient specific values of k, as a proxy for the effect of the increased pulmonary flow F due to APC [12]. In the established CFD, the solution was considered converged when continuity and velocity-scaled residuals dropped below < 10 −5 . In pulsatile simulations, 4000 iterations per cycle were used. In the lean CFD, the default internal algorithms were used. User-specified measurements of power loss, total and hepatic mass flow to the left pulmonary artery were monitored for convergence during the solution process.
In all pulsatile simulations inlet velocity waveforms acquired by MRI were decomposed in 15 harmonics using Matlab R2019a (The MathWorks Inc., Natick, MA, USA).
In TCPC patients, results were reported from the 6th cycle to reach a cyclic steady-state in the hepatic flow from IVC to the pulmonary arteries. In Glenn patients, results were reported from the 3d cycle to reach a cyclic steady-state.

Meshing
Systematic mesh convergence studies were performed for both methods. More details on meshing have been described previously [19,20]. In the established CFD we found mesh-independent results with tetrahedral element sizes of approximately 1/15 the average pulmonary artery diameter, corresponding to 2-3% of the IVC diameter. On average, this represented 1.6 million elements in the TCPC patients and 1.1 million elements in the Glenn patients. Elements were constructed semi-automatically in ICEM CFD v19.2 (Ansys Inc., Canonsburg, PA, USA). The established CFD program reported no invalid elements, typical aspect ratios of 0.8 to 1.0 for 85% of the elements (min aspect ratio typically > 0.45) and no skewness angles > 75° (max skewness angle typically < 65°).
In the immersed boundary method used by the lean method, the software automatically generates hexagonal elements within the minimal volume required to completely immerse the model in elements, of which the vast majority will be perfectly hexagonal. The immersed boundary method then "cuts" any corners that protrude outside the computational domain using well-understood correction methods [16]. Mesh-independent results were found with approximately 150,000 elements for TCPC patients and approximately 140,000 for the Glenn patients. This corresponded to a mesh size of approximately 10% of the IVC diameter.

Statistics
Statistical analysis was performed using GraphPad (v8.0, La Jolla, CA, USA). Linear regression was used to analyze correlation between calculated results and the reference method. Accuracy and precision were calculated according to Bland-Altman analysis as mean ± SD of the measured difference against the mean of the results. Results with p < 0.05 were considered significant.

Hardware
A workstation with an 8-core Intel i7700k processor @ 5.0 GHz was used for all steady-state simulations (lean and established CFD), and for the pulsatile simulations using the lean method.
The pulsatile CFD simulations in the established CFD were performed on a computation server with 44 compute cores (2 × 22-core Intel Xeon Gold 6152, Intel Corporation, Santa Clara, CA, USA), of which 32 cores were used for the CFD simulations. where p is the static pressure, is the blood density, v is the velocity, and A is the cross-section area of all inlets and outlets. 4. Wall shear stress, WSS ( w ): where μ is the dynamic viscosity and u y =̇ is the rate of shear.

Time average wall shear stress (TAWSS):
where T is the physical time of one heartbeat and τ w is the local wall shear stress. 6. Pulsatile index (PI) [21]: where Q max , Q min , and Q avg are the maximum, minimum, and average flow rates. 7. Weighted pulsatile index (wPI) to represent an overall pulsatility level of the Fontan connection [21]: where C i is the relative flow split of vessel i and Normalized wall shear stress area was calculated as the ratio of the area with WSS or TAWSS < 0.4 Pa, compared to the total vessel wall area. This measure has been suggested to be associated with increased risk of atherosclerosis [22] and has been previously reported in Fontan patients [23].

Results
Qualitative comparison of velocity-coded flow patterns from the lean, steady established, and pulsatile established CFD showed essentially identical flow patterns in all patients (Fig. 2). Visual comparison of WSS from the lean, steady established, and pulsatile established CFD also showed very similar patterns (Fig. 3). Simulations showed a posterior vortex in the SVC/IVC junction in six of the nine TCPC patients.
Linear regression analysis of %PFD LPA from the steady lean and pulsatile established CFD both showed a correlation coefficient of 0.94 with the MRI data. Bland-Altman analysis of %PFD LPA from the steady lean and pulsatile established CFD compared to MRI data showed a bias of -1.9 ± 3.4% for the lean CFD and -1.8 ± 3.1% for the established CFD (Fig. 4).
The correlation of weighted pulsatility index (wPI) with the differences in %HFD LPA and Ė loss between lean steadystate and pulsatile simulations are shown in Fig. 6.
In the lean CFD, meshing was part of the automatic solver process and took on average 10 ± 2.8 s per model. In the established CFD, semi-manual meshing took on average 21 ± 4.1 min per model. The process of setting up the CFD model, such as adding boundary conditions, preparation of measurements and convergence criteria took 15 ± 3 and 33 ± 5.6 min for the lean and established CFD, respectively ( Table 2). Among Glenn patients, solver times were for the steadystate lean and pulsatile established CFD 3.0 ± 2.0 min and 203 ± 65 min (3.4 ± 1.1 h), respectively. Among TCPC patients, solver times were for the steady-state lean and pulsatile established CFD 3.7 ± 2.3 min and 570 ± 115 min (9.5 ± 1.9 h), respectively ( Table 2).
Additional results in terms of user-and solver time, as well as separately reported hemodynamic results for Glenn and TCPC patients using the lean and established steady-state and pulsatile CFD are shown in Tables 2, 3, and 4, respectively.

Discussion
This study shows that the proposed lean CFD approach and the established CFD approach produced very similar results, but with vastly shorter solution times for the lean, steady-state simulations. Our findings indicate that steadystate simulations using the lean method can be performed in approximately 4% of the time required for achieving nearly identical results using the established method with pulsatile inflows, whereas previous findings of Khiabani et al. showed savings of approximately 50% [21].
While the results show that the lean CFD approach can reproduce in vivo measurements from MRI and compares well with established CFD, clinical research is needed before broad application in clinical routine to show that the CFD simulations are accurate in their prediction of outcome after surgery and catheter interventions. The first study comparing CFD with post-operative outcome was recently published by Trusty et al. and showed in a retrospective analysis of 12 patients (whereof 7 Fontan completion surgeries) a fairly large bias between predicted flows and outcome [24]. This may be explained by changes in the patient's physiology, thereby changing the boundary conditions. Examples of physiological changes are patient's growth as well as changes in oxygenation, collateral flow and cardiac output. Validation of the possibility to predict intervention outcome using the lean CFD approach is thus needed.

Benefits of the Lean Method
The proposed "lean" numerical pipeline uses less computational resources and less user input and thus use less resources to create value compared to the established methods. The three major differences that we identified making the proposed lean CFD faster and "lean" compared to established CFD are: 1. In the immersed boundary method used by the lean method, the software automatically generates hexagonal elements and benefits from the numerical advantages  that follow, since hexagonal elements are more efficient compared to tetrahedral elements. The established method uses tetrahedral elements which require greater user attention to the mesh density, mesh quality, and the distribution of element sizes near the center and walls of the model. Thus, the lean method can obtain a correct solution faster, with less user interaction and with fewer elements than the established method. 2. The lean method integrates CFD with the geometric editing tool, requiring 0 transfer of data between softwares. This significantly saves time when iterating designs. 3. The lean CFD method has a process-oriented workflow making the use of software easier.
While Glenn and Fontan procedures are not emergency procedures, performing faster CFD for interventional planning with less resources are inherently valuable even for elective surgery.
Regionalization and centralization are increasingly adopted by western countries to promote efficiency, reduce excess mortality, and health care expenses. Thus, a tertiary referral center will cover a large geographic area and reducing the number of outpatient visits is important. With lean CFD, there is enough time to perform predictive simulations even when the pre-interventional MRI is performed after the patient has arrived for elective surgery. This may be a time-efficient way to obtain updated information prior to the surgical intervention and can save time and cost of travel for the family.
In addition, some patients require catheter-based interventions between and after the major Fontan surgical stages. Lean CFD and its short solution times may aid in interventional decision-making based on acute findings during such procedures. If MRI and catheterization is performed under the same general anesthesia, time is short to update the calculations.
Our reported simulation times of less than 5 min mean that such simulations can be performed well within the timeframe of ongoing invasive interventions. These reported time savings are likely conservative since the lean simulations ran on an 8-core computer, whereas the pulsatile simulations of the established method ran using 32 cores on a dedicated server.
While it was shown that the lean CFD approach has computation times as short as 3.0 min in average, another benefit for interventional work is the integration of the lean Fig. 6 a Hepatic flow percentage (HFD) difference between lean steady-state and pulsatile simulations vs. weighted Pulsatile Index (wPI). b Relative power loss error between lean steady-state and pulsatile simulations vs. weighted pulsatile index (wPI) In the lean CFD, the workflow in the user interface involves graphical interaction with the model and does not expose the user to complex details of CFD simulation, thus making it easier to use in a clinical setup. Many established CFD settings are available in the lean approach but are not mandatory for the user to interact with.
The lean CFD approach can thus be used in clinical studies aiming to demonstrate the effect of proposed anatomical interventions on blood flow and offers the possibility to do clinically integrated simulations of surgical or catheterbased interventions in patients with Glenn and Fontan circulations at the hospital without use of advanced computational facilities.
A user of predictive CFD in a clinical setting should still be an expert user with demonstrated skills in the tools and procedures involved and this does not differ in the investigated methods. We believe that it will be easier to have more expert users with the productivity benefits of the lean method and the rapid iterations it enables with the

Hemodynamic Measurements
We used several hemodynamic parameters for comparison between the techniques, some used in clinical routine imaging, such as pulmonary flow distribution. Other measures, such as hepatic flow distribution, are clinically relevant but only available through more advanced techniques, e.g., 4D-flow MRI or CFD. We also compared more advanced hemodynamic measures ( Ė loss , WSS) that have been used in previous studies of single ventricles and shown to be of interest but have not reached clinical practice [17,23,25]. One reason for the lack of clinical use is the complexity in obtaining these values. Thus, our lean approach to CFD in single ventricle patients may lead to increased availability and the possibility to use some of these advanced hemodynamic measures in larger cohorts and demonstrating the potential clinical value. We found a pulmonary distribution (%PFD LPA ) of 40.1%, which is similar to 43% measured by Tang et al. in a previous study with a notably large cohort (n = 108) [26]. Our measurements for the hepatic distribution (%HFD LPA ) were 43% on average, similar to 50% calculated by Wei et al. [23]. Moreover, we calculated a power loss ( Ė loss ) of 1.3 ± 1.4 mW and 1.4 ± 1.4 mW on average for steady and pulsatile simulations, respectively. These values are similar to previous studies, e.g., Baretta [28]; and Marsden et al. found 6.7 mW and 13.9 mW in two patients [29]. Normalized WSS area, which is the area of WSS below 0.4 Pa divided by the total surface area was similar for the lean CFD (22 ± 18%) and established CFD (21 ± 17%). These results were similar to a study of Wei et al. [23] that reported normalized WSS area of 22%.

Pulsatile vs. Steady-State Simulations
We calculated time averaged venous flows as the mean flows from the flow curves obtained in MRI during the cardiac cycle at supine rest during free breathing. However, the flow patterns in our Glenn and Fontan patients are much less pulsatile compared to venous flow in healthy subjects (illustrated in the pulsatility index). As seen in Table 2, computation times are significantly shorter for time-averaged simulations compared to pulsatile flow simulation.
With respect to hepatic flow and pulsatility, Wei et al. [13] showed that "a large portion of the cases can use timeaveraged BCs to save computational cost," and we found that differences in our patients fell within the lower range of this study (Fig. 6a). Wei et al. showed that such differences are mainly determined by weighted pulsatility index (wPI) and the IVC angle. We speculate that the found differences between the studies could be due to low patient number and the relatively perpendicular IVC insertion angle relative to the pulmonary arteries in the current study (Fig. 2), whereas Wei et al. had a wider range of angles. If required, the lean approach can significantly shorten pulsatile TCPC simulations which were shown to take 1.3 ± 0.5 h whereas the established CFD required 9.5 ± 1.9 h on a much faster computer (Table 2). Additionally, our findings of the relationship between power loss and pulsatility were similar to the findings of Khiabani et al. [21] (Fig. 6b).

Limitations
User dependence and limited spatial resolution with possible inaccuracies in the identification of the vessel anatomy and blood flow rate will affect the results. However, our validation used MRI which is the clinical reference standard for measurement of intrathoracic flow, and this showed low bias and variability. Further limitations are simplifications such as assumed rigid walls and inlet plug flow profiles, but these are established assumptions in research studies of patient-specific CFD simulations [8,17,[30][31][32].
Aiming to predict changed pulmonary artery flows following interventions, previous studies have simulated more advanced patient specific physiology at the outlets in different ways. Multiple studies have used compliant Windkessel models of PVR and "lumped parameter" models of the complete cardiovascular system in the CFD loop [26,33]. These methods require more patient data as input and are more resource intensive.
The demonstrated lean method does not predict pulse pressures since it relies on non-invasive MRI data. If catheterized pressures are available, such data could potentially be used with the lean method to include patient-specific linear PVR and atrial reference pressure, but this has not been validated with patient data or compared with results from advanced models such as proposed by Ahmed et al. [33].
Given the limitations of using fixed PVR, our simplified pulmonary circuit including the effects of precapillary aortopulmonary collateral (APC) flow correlated well with the pulmonary flow split measured in patients with a correlation coefficient of 0.94 for both the lean and established CFD (Fig. 4) and has been successfully used to predict outcome of interventions [12].
While this model can be considered a fully non-invasive and lean approach, it should be considered a proxy for resistance and cannot predict absolute patient-specific values of Fontan pressures. However, for the purposes of calculating power loss, wall shear stress, and total and hepatic pulmonary flow distribution, this is not a limitation since they are not based on absolute pressures.