QUANTIFICATION OF ANATOMICAL AND FLUID-DYNAMIC ANOMALIES IN FONTAN PATIENTS BASED ON MAGNETIC RESONANCE IMAGING

Background. Univentricular diseases are lethal congenital diseases affecting about 2 % of newborns in the western world. Due to these pathologies, only one ventricle pumps blood into the circulatory bed, and arterial and venous blood are mixed, preventing from properly providing tissues and organs with oxygen. These pathologies are currently treated through the so-called Fontan procedure, which is a multi-step and complex surgical approach. The Fontan procedure aims at obtaining the anatomical separation between the systemic and pulmonary circulations, and hence between oxygenated and non-oxygenated blood. However, the only ventricle present in the heart remains the only pumping organ, and blood flow in the pulmonary circulation is merely passive. Also, and importantly, the postsurgical anatomy of the junction between systemic veins and pulmonary arteries is markedly non-physiological. As such, it is associated with altered blood fluid dynamics, undesired energy losses, and, ultimately, sub-optimal quality of life and short life expectancy. Objective. On this basis, clinicians need tools to 1) quantify the post-surgical anatomical and fluid-dynamic alterations, 2) correlate these anatomies to the patients’ prognosis, and 3) identify criteria to improve Fontan surgery. Methods. In order to support the pursue of these goals, we developed a computational tool for the processing of 4D flow data, i.e., phase contrast magnetic resonance images yielding information on the velocity of tissues within a 3D domain. The tool allows for reconstructing the 3D geometry of the surgically treated anatomical district and, through a semi-automated user-interface, extracting relevant geometrical scores, as well as quantifying flow rates in the different vessels, energy losses, and wall shear stresses. A numerical method based on the finite element approach was implemented to estimate relative pressures. Results. The developed tool was preliminarily applied to the analysis of the datasets of six pediatric patients. The analysis of data obtained by two independent users highlighted a good repeatability of geometrical reconstructions, and hence of the quantification of geometrical scores. The method for the quantification of relative pressures was preliminarily tested in a simplified model of the thoracic aorta, with encouraging results. Conclusions. The developed computational tool, which, to the best of our knowledge, is completely novel, helps clinicians to quantify the post-surgical anatomical and fluid-dynamic alterations. Ongoing activities include its application to the real datasets, and the extension of the analysis to a wider cohort of patients, so to check for correlations between the quantitative geometrical and fluid-dynamic indexes with the patients’ prognosis. Such possible correlations could help identifying criteria to improve Fontan surgery.


Introduction
Univentricular diseases are lethal congenital heart diseases that affect every year about 2 % of newborns in the United States [1], and require surgical treatment. Early surgical approaches consisted in palliative treatments (i.e, BT shunt and BDG shunt) conceived to save the patient's life in the short term. Advances in cardiac surgery have led to the approach that is the current state of the art, i.e., the Fontan procedure ( Fig. 1). Through a multi-step surgery, the Fontan procedure aims to a long term solution by reconfiguring the native structures involved in the pathology. The inferior vena cava (IVC), i.e., the vein returning venous blood from the lower systemic cir-culation is connected to the pulmonary arteries (LPA and RPA, respectively) either through the right atrium (atrium-pulmonary connection or APC) or through a graft bypassing the right atrium (total cavopulmonary connection or TCPC) and through the superior vena cava (SVC). In this way, the blood pumped by the left ventricle feeds the peripheral organs and is directly redirected to the lungs without the action of a proper right ventricle, and successively returns into the left atrium through the pulmonary veins.
Unfortunately, this configuration is far from being physiological [2]: the hydraulic impedance downstream of the left ventricle is much higher than the physiological one, and, due to the abnormal anatomy of the reconstructed APC/TCPC, local fluid dyna-mics is altered and characterized by large eddies and stagnation regions that may lead to thromboembolitic events. Also, sub-optimal distribution of blood flow rates between the left and right pulmonary arteries can be obtained, leading to atrophic remodeling of the branch receiving poor flow rate and hence to an increase in its hydraulic impedance [3].
Based on these evidences, the Fontan procedure should be performed aiming at optimizing the balance between flow rates into the pulmonary arteries, as well as to minimize geometrical distortions that could lead to flow disturbances. Also, flow rate redistribution and energy loss at the APC/TCPC, which is a surrogate measure of flow disturbances, could be used as prognostic indexes to judge the postsurgery evolution of the patient.
In this work we hypothesized that blood fluiddynamics with the APC or TCPC of Fontan patients could be quantified based on 4D flow imaging, i.e., phase-contrast magnetic resonance imaging sequences that yield the 3D velocity components that characterize biological tissues within a volume of interest [4]. In particular, we developed in house software to exploit such information to quantify the 3D geometry of the surgically treated anatomical district and to computing the flow rates in the pulmonary arteries, local viscous energy losses, wall shear stresses, and pressure drops.

Problem Statement
Our main objective is to develop numerical approaches and corresponding computational tools for the processing of 4D flow data aimed to reconstructing the 3D geometry of the surgically treated anatomical district and to extracting relevant geometrical scores, as well as quantifying flow rates in the different vessels, energy losses, wall shear stresses and relative pressure distributions. It helps clinicians to quantify the post-surgical anatomical and fluid-dynamic alterations, correlate these anatomies to the patients' prognosis, and identify criteria to improve Fontan surgery.

Patients Population and Image Acquisitions.
Cardiac magnetic resonance imaging was performed on 6 pediatric patients (age 11-42 years, 2 males) previously treated through TCPC Fontan surgery. 4D flow acquisitions were performed using a 1.5 T Siemens AERA System, during free breathing, with prospective ECG gating and respiratory navigator. Images were acquired on a stack of sagittal planes to scan the volume of interest (VOI). Images in-plane resolution and thickness were equal to 1.9-2.3 mm and to 2-2.4 mm, respectively. The dimensions of the resulting voxels hence range from 1.91.92.0 to 2.32.32.4 mm 3 . For each plane, three images corresponding to the three velocity components of tissues were obtained (Fig. 2). Velocity-encoding (VENC) was set to 100-150 cm/s, depending on the specific patient, in order to optimize the signal-tonoise ratio. All data were exported in DICOM format.
Image Pre-Processing. As for the subsequent operations, image pre-processing was carried out through in house software implanted in Matlab® (Matlab, The Mathworks, Inc.) and embedded within an ad hoc graphical user interface.  1. Physiological anatomy of the heart (a); heart anatomy in case of univentricular disease (b); heart anatomy following APC (c); heart anatomy following TCPC (d) Images were processed by three consecutive filtering operations, aimed to preventing aliasing, to removing possible off-sets in velocity values and background noise effects, and to smooth velocity data within the VOI [5]. The images related to the single velocity components were processed so to obtain images whose greyscale represented the velocity magnitude at each voxel. Such velocity magnitude images were used for the subsequent segmentation of anatomical structures.
Segmentation of Anatomical Structures. On each relevant plane of the n velocity magnitude images, the contour of the connection between the pulmonary arteries, the inferior vena cava and the superior vena cava was manually traced. Typically, 20-30 points were manually positioned on the image in this process. A polygonal geometry was automatically obtained by connecting the selected points, and every voxel whose center fell within the polygon was selected as part of the fluid-dynamic domain to be analyzed. The whole manual tracing procedure typically required about 10 minutes for one subject.
The entire set of selected voxels was automatically assembled, its boundary was extracted in the form of an isosurface and stored as a triangulated surface mesh. The latter was smoothed through a Laplacian filter to obtain the final 3D geometry of the vessels' wall encompassing the final fluid domain. After manually labeling each of the four vessels of interest, i.e., IVC, SVC, LPA, and RPA, skeletonization of the fluid domain was automatically performed to obtain the corresponding centerlines [6].
Computation of Geometrical Indexes. Two geometrical indexes were computed: the caval off-set, i.e., the distance from the junction between the IVC and the pulmonary arteries to the junction between the SVC and the pulmonary arteries, and the angles characterizing the connection of the four vessels involved in the Fontan procedure. This computation was carried out through the following automated steps: 1) for each labeled vessel, the centerline was scanned moving from the vessel's distal end towards the junction. At each voxel of position 2) the centerline point h X characterized by an abrupt change in cross-sectional hydraulic diameter was identified; 3) the positions of the last voxels, within 10 % of the centerline length, scanned along the centerline throughout steps 1) and 2) were approximated by cubic splines. The local tangent h  t to the spline at the point h  X representing the approximation of h X was computed; 4) given the tangent vectors h  t and k  t associated to the h-th and k-th vessel, the angle of their connection was simply computed as: The caval distance was computed as the distance between IVC X and SVC X (Fig. 3).

Computation of Fluid-dynamic Indexes.
The following fluid-dynamic indices were computed: 1) Flow rates into each vessel. To this aim, for each vessel of label h the plane h   normal to the tangent vector h  t was identified, and the corresponding contour of the vessel wall was identified as previously explained. The flow rate through h   was computed as: where A is the area encompassed by the contour, N is the number of voxels falling within it, and v i is the velocity vector at the i.th voxel of interest.
The ratio between flow rates in the four vessels was then computed.
2) The 3D distribution of the wall shear stresses (WSS) acting on the vessels' wall, which was computed following a recently published method, which exploited the use of 3D sobel kernels to obtain the position-dependent velocity gradient [7].
3) The power wasted due to viscous dissipation at the TCPC (W w ), which was computed following the approach recently proposed by [8].
4) Pressure distribution within the lumen of the vessels. We developed a special computational finite-element-based module in our tool for determining of pressure distributions. For this purpose, the module solves the equation p b    (where the pressure gradient b  is derived directly from the Navier -Stokes equation) and accomplishes it in the following steps: 1) partitioning of the actual computational aortal domain into finite elements with consequential computing of the pressure gradient; 2) finite-element solving the equation mentioned above. This method was preliminary tested on CFDderived velocity fields, following a previously published sub-sampling procedure to discretize CFD data [7]. In particular, it was tested within a simplified model of aorta with a grid discretization equal to 2 mm, to estimate time-dependent pressure distributions, and within a more realistic model of the aorta with a grid discretization equal to 1 mm and 2 mm, to analyze the influence of the domain discretization on pressure estimations.

Results and Discussion
Geometrical Indexes. The isosurface and the skeletonization of the wall of the four vessels are reported in Fig. 4, highlighting the geometrical complexity and heterogeneity of the anatomical district of interest of the Fontan connection. The automatic computation of the caval off-set, and of the angles characterizing the connection of the four vessels in- volved in the Fontan procedure, yielded the data reported in Table. In particular, a maximum caval offset equal to 3.7 cm was found in patient 3, being almost double with respect to the corresponding datum obtained for the other patients; values lower than 1.3 cm characterized patients 1, 2 and 6, indicating a low-distorted caval anastomoses. As regards the angles between the four vessels of the connection, these were characterized by high inter-patient variability: IVC-LPA SVC-LPA IVC-RPA SVC-LPA , , ,     ranged from 69 to 109, from 74 to 131, from 46 to 121 and from 85 to 134, respectively. Of note, 33 % of the measured angles fell in the specific 90  10 range.
Fluid-dynamic Indexes. The high variability of the geometrical indexes reflected into a notable variability in flow rates flowing through the connection and into the vessels. For instance, peak flow rates of the IVC and SVC ranged from 0.35 and 2.17 L/min and from 0.17 L/min to 1.16 L/min, respectively. Computing the ratio between peak IVC Q or SVC Q and the global inlet flow rate entering the connection, i.e., IVC+SVC , Q highlighted that contribution of flow rate from the IVC was predominant for every analyzed patient but patient #3 (Table). Also, flow rate was unevenly redistributed to the pulmonary arteries for every analyzed patient but patient#5, who was characterized by a ratio between peak LPA Q or RPA Q and the global flow rate flowing out of the connection, i.e., LPA+RPA , Q equal to 0.49. As regards the computation of 3D WSS distributions, the highest WSS value (1.1 Pa) was found in the lower venous duct of patient# 4, while in patients #1 and #5 WSS reached values up to 0.65 and 0.60 Pa, respectively (Fig. 5). In the other patients, peak WSS values were lower than 0.45 Pa.
Moreover, the computation within the TCPC of the temporal mean of W w highlighted that the conditions of patients #3 and #5 were characterized by viscous dissipations equal to 7e-05 mW and 3.2e-05 mW, nearly one order of magnitude lower with respect to the other patients, who were characterized by an average value of 21.3 e-05 mW.
Pressure Distribution within the Lumen of the Vessels. For the simplified model of the aorta, results estimated in five equidistant consecutive time frames, centered at peak systole, are shown in Fig. 6 (from left to right). The pressure distributions well resembled, both qualitatively and quantitatively, those obtained by Krittian et al. [9]. The results obtained for the more realistic model of aorta, characterized by a sub-sampling discretization equal to 1 mm and 2 mm for isotropic voxels, corresponding to about 270'000 and 39'000 nodes, respectively, are shown in Fig. 7. In both cases, the approach demonstrated to achieve a good qualitative correlation with the original CFD simulations.

Conclusions
At the current state of the art, the developed Matlab tool represents the first investigation platform including all the fundamental steps of 4D-flow MRI data analysis: image opening and viewing, segmentation and post-processing analysis (both geometrical and hemodynamical). Moreover, the developed software offers a user-friendly graphical interface (GUI) for a simple and intuitive managing of the different implemented functions, thus provi- ding a ready-to-use platform for the daily clinical practice. Therefore, it is the first release of a comprehensive analysis tool able to promote a larger clinical use of the 4D-flow MRI technique while studying and featuring the fluid-dynamics in Fontan patients' TCPC; thanks to an enhanced usability of data which are difficult to get by non-expert operators, the developed software may elucidate the comprehension of the effectiveness of the Fontan circuit and help in detecting any pathophysiologic situation (i.e. arte-riovenous malformations and pulmonary arteries stenosis).