The Hemodynamic Effect of Modified Blalock–Taussig Shunt Morphologies: A Computational Analysis Based on Reduced Order Modeling

Abstract

The Modified Blalock Taussig Shunt (MBTS) is one of the most common palliative operations in case of cyanotic heart diseases. Thus far, the decision on the position, size, and geometry of the implant relies on clinicians’ experience. In this paper, a Medical Digital Twin pipeline based on reduced order modeling is presented for fast and interactive evaluation of the hemodynamic parameters of MBTS. An infant case affected by complete pulmonary atresia was selected for this study. A three-dimensional digital model of the infant’s MBTS morphology was generated. A wide spectrum of MBTS geometries was explored by introducing twelve Radial Basis Function mesh modifiers. The combination of these modifiers allowed for analysis of various MBTS shapes. The final results proved the potential of the proposed approach for the investigation of significant hemodynamic features such as velocity, pressure, and wall shear stress as a function of the shunt’s morphology in real-time. In particular, it was demonstrated that the modifications of the MBTS morphology had a profound effect on the hemodynamic indices. The adoption of reduced models turned out to be a promising path to follow for MBTS numerical evaluation, with the potential to support patient-specific preoperative planning.

1. Introduction

Congenital heart defects are currently present in about 9 of every 1000 live-born children [1,2]. Several of these pathologies, such as pulmonary atresia and Tetralogy of Fallot, are linked with cyanotic heart diseases in which low levels of oxygen in the blood are encountered. Especially in the case of the pulmonary atresia, the pulmonary valve plane is absent, leading to full blockage of the blood flow towards the lungs. Among the most diffused approaches to tackle these complications, the modified Blalock–Taussig shunt (MBTS) is a surgical procedure which consists of implanting a synthetic shunt between the subclavian and pulmonary artery. The MBTS is frequently adopted for the treatment of newborn children suffering from pulmonary hypoperfusion triggered by congenital heart diseases [3]. In this way, sufficient levels of oxygenated blood are delivered to the pulmonary circulation. Although the implanted shunt does not guarantee the holistic repair of pulmonary vascular blockages, it offers a time extension until the child grows in order to subsequently proceed with a permanent surgical treatment [4]. In general, the MBTS is considered as a low-risk management option for the palliation of congenital heart diseases. Nevertheless, the danger associated with this intervention should not be downplayed, as the overall mortality and composite morbidity rates remain relatively high at 7.2% and 13.1%, respectively [4,5].

To date, in clinical practice the appropriate MBTS implant is specified by exclusively considering its cross-sectional diameter. In fact, surgeons select the optimal MBTS implant for each patient according to their experience, their evaluation of the patient’s condition and weight, and other clinical factors [6]. Beyond these aspects, several works from the literature have underlined the importance of a patient-specific approach that personalises treatment for each patient [7,8]. In practice, if the shunt is undersized, the thrombosis risk soars [9]. On the contrary, if the shunt is oversized, a state of pulmonary overflow with consequent diastolic pressure and systemic perfusion decrease is detected. From the clinical state of the art, it appears clear that MBTS sizing is pivotal in managing complex congenital heart diseases [10,11].

In light of this necessity, the exploration of new tools for MBTS fluid dynamic analysis appears to be very challenging. In recent years, Computational Fluid Dynamic (CFD) techniques have been increasingly used to study cardiovascular systems and procedures within both the systemic [12,13,14] and pulmonary circulation [15,16,17]. In this clinical context, the goal is to develop valuable engineering tools for the evaluation of post- and pre-operative information. Nowadays, modern imaging techniques are supported by advanced multidisciplinary tools and are able to fill the gap between the world of engineering and clinics [18,19,20,21] in terms of both morphological and functional information. This trend has been followed for specific studies on MBTS procedures as well. In this context, there are several studies reporting computational methods for the analysis of MBTS implants [7,22]. Nevertheless, only patient-specific cases with central and MBTS configurations have been considered, with limited geometrical exploration.

Despite the recent advances in terms of computational power and algorithms, high-fidelity CFD simulations remain prohibitively expensive with regard to the CPU time and the memory demand. As an alternative, Reduced Order Models (ROMs) have spread in recent years, although mainly in classical mechanics. Lately, this concept has been successfully applied in the biomechanics field for the evaluation of the hemodynamics of parametrised ascending thoracic aneurysms [23]. To date, the ROM approach has never specifically been used as a tool for the analysis of the hemodynamic parameters of an MBTS case. The adoption of model reduction methods could contribute to a remarkable improvement in terms of the speed and usability of these computational tools in the clinical field [24]. To the best of our knowledge, the potential of using computational methods together with mesh morphing and ROM tools for the analysis of MBTS morphology has not yet been explored. As mentioned, the current literature concerning for numerical analysis of MBTS uses the straightforward approach of full CFD analyses, while the adoption of ROMs for fluid dynamics within the context of cardiovascular shunting has not yet been explored and requires full study and validation. The same holds true for mesh morphing methods. The possibility of modifying MBTS morphology thanks to the mesh morphing approach has not previously been explored. The verification of such approaches by exploiting new cutting-edge computational tools could provide a significant contribution in the context of computational analysis of MBTS.

Considering the current status of the literature, the present study aims to define a new and fast method to provide computational analysis of different MBTS complex configurations. The objective is to combine ROM, CFD, and mesh morphing techniques to overcome the computational expense and obtain high-fidelity analysis of a wide range of MBTS shapes. In the next section, the theoretical foundation of the morphing and ROM approaches is presented. Then, the computational model setup for the MBTS is described in the material and methods section. Afterwards, the results of our analysis are displayed and discussed. Finally, the study’s conclusions are reported.

2. Theoretical Background

2.1. The Mesh Morphing Technique

Radial basis functions consist of a powerful mathematical tool for interpolating data which has further proven to be reliable and accurate in the mesh morphing of a discretised computational model’s domain [25]. The mesh morphing action is performed by imposing a desired displacement on a set of points, named source points. Then, a displacement field is interpolated and applied to the nodes of the computational model, named as target points. The interpolation result depends on the selection of the radial function and on the type of support (local or global) that it offers. The nucleus of the RBF theory is based on the calculation of unknown weight functions and polynomial terms from known data (displacement at source points in case of mesh morphing). The RBF approach can be applied to multidimensional problems as well [26,27].

The general form of the interpolation function consists of a radial basis and a polynomial term, as shown in Equation (1):

$$s\left(\mathbf{x}\right)=\sum _{i=1}^N{\gamma }_i\varphi \left(\right||\mathbf{x}-{\mathbf{x}}_{ki}\left|\right|)+h\left(\mathbf{x}\right).$$

(1)

The scalar function $s\left(\mathbf{x}\right)$ is defined for each arbitrary point in space $\mathbf{x}$ and represents the point’s displacement after the morphing action. The movement of a point can be considered as the aggregation of its interactions with all the source points ${\mathbf{x}}_{ki}$; this is known as the Euclidean distance between source ${\mathbf{x}}_{ki}$ and target points $\mathbf{x}$ multiplied by the radial function $\varphi$ and the coefficient ${\gamma }_i$. The additional polynomial factor $h\left(\mathbf{x}\right)$ is introduced to guarantee the existence and the uniqueness of the solution.

The linear system becomes solvable with the fulfillment of a desired displacement, ${\mathbf{g}}_i$, at the source points and the implementation of the orthogonality conditions:

$${\displaystyle s\left({\mathbf{x}}_{ki}\right)=g_i,1\le i\le Nand\sum _{i=1}^N{\gamma }_{i}h\left({\mathbf{x}}_{ki}\right)=0.}$$

(2)

According to [28], a unique interpolation polynomial exists if the basis function is a conditionally positive definite function. Assuming that the basis functions are conditionally positive definite of order $m\le 2$, the linear polynomial can be expressed in 3D space as follows:

$$h\left(\mathbf{x}\right)={\beta }_1+{\beta }_{2}x+{\beta }_{3}y+{\beta }_{4}z.$$

(3)

The weights of the linear polynomial of Equation (3), ${\beta }_i$, are calculated such that the given value at the source points ${\mathbf{x}}_{ki}$ can be retrieved exactly. It is crucial to underline that mesh morphing is obtained through a vector field displacement, while RBF are scalar functions; therefore, the interpolation of a 3D set of displacement at source points can be performed using the following equations:

$${\displaystyle s_x\left(\mathbf{x}\right)=\sum _{i=1}^N{\gamma }_i^x\varphi (\mathbf{x}-{\mathbf{x}}_{ki})+{\beta }_1^x+{\beta }_2^{x}x+{\beta }_3^{x}y+{\beta }_4^{x}z.}$$

(4)

$$s_y\left(\mathbf{x}\right)=\sum _{i=1}^N{\gamma }_i^y\varphi (\mathbf{x}-{\mathbf{x}}_{ki})+{\beta }_1^y+{\beta }_2^{y}x+{\beta }_3^{x}y+{\beta }_4^{y}z.$$

(5)

$${\displaystyle s_z\left(\mathbf{x}\right)=\sum _{i=1}^N{\gamma }_i^z\varphi (\mathbf{x}-{\mathbf{x}}_{ki})+{\beta }_1^z+{\beta }_2^{z}x+{\beta }_3^{z}y+{\beta }_4^{z}z.}$$

(6)

The morphing action can be restricted in the space by establishing a domain beyond which the mesh remains unaltered.

RBF techniques are gaining more and more ground in the computer-aided engineering (CAE) field. Their application has already spread in the automotive [29], naval [30], aeronautical [31], energy [32], and medical fields [33]. Their main advantage is their potential to morph the mesh without changing its topology or its consistency while retaining disk usage at low levels. Thanks to this feature, the time consuming re-meshing process can be avoided. Currently, RBF-based mesh morphing allows for precise control of the source points’ final position while at the same time preserving a high degree of mesh quality.

2.2. Reduced Order Modeling

ROMs offer a prompt approximation of the exact solution of complex non-linear problems, including fluid dynamics problems, guaranteeing low error levels [34]. There are two main techniques for building ROM, namely, the model-based and data-driven approaches. In this study, the latter method is adopted, and is described in the current subsection. In particular, the ROM workflow can be split into two main sections, the ROM creation and ROM consumption parts, as illustrated in the workflow shown in Figure 1.

The ROM creation section involves a computational solver block and an ROM Builder block. Initially, the determination of the set of the input parameters and the space to be explored in a particular problem is performed. Afterwards, a list of several possible parameter combinations, named scenarios, is built up and utilised for feeding the training of the ROM model. In order to carry out this operation in a systematic way, a Design of Experiments (DOE) method should be introduced to detect the optimal selection of scenarios inside the predefined bounded space. The term optimal implies defining the least number of scenarios that can be distributed in the space in such a manner that the flow behaviour of the computational model is captured with sufficient accuracy [35]. After specifying the DOE, the corresponding simulations and parametrical data are generated in the form of snapshots. These snapshots are sets of matrices in which the parametric fluid dynamic variables of the surface/volume of the domain are stored. For the build-up of the ROM model, snapshots are compressed into a small number of modes via singular value decomposition (SVD) [36]. The SVD approach is based on the generalization of diagonalization for non-square matrixes, which was discovered by Eugenio Beltrami and Camille Jordan in 1873 [37]. The nucleus of the SVD theory is nested in the fact that every rectangular matrix ($\mathbf{A}$) can be decomposed into three matrices: two orthogonal ($\mathbf{P}$, $\mathbf{Q}$) and one diagonal ($\Delta$).

$$\mathbf{A}=\mathbf{P}\Delta {\mathbf{Q}}^{\top }.$$

(7)

Assuming that the dimension of matrix $\mathbf{A}$ is n×d, then, the dimensions of the matrices $\mathbf{P}$, $\Delta$, and ${\mathbf{Q}}^{\top }$ are equal to n×n, n×d, and d× d, respectively. The diagonal values of matrix $\Delta$, known as singular values, are decreasingly classified according to their statistical representativity. By exploiting this characteristic, the SVD algorithm offers an approximation of matrix $\mathbf{A}$ by neglecting part of the matrices $\mathbf{P}$, $\Delta$, and ${\mathbf{Q}}^{\top }$.

The tolerance error is determined by the number of selected singular values, r. The result of the approximation is the modes, which are used to represent the fluid flow data. Each dataset case can be linearly expressed through a set of scalar coefficients by exploiting the genetic aggregation response surface (GARS) technique [38]. GARS is a method that aims at selecting the best response surface for a particular design of experiments. It solves several surrogate models, performs aggregation, and selects the best response surface according to the cross-validation errors. The calculation of the variable v on a selected point x(x,y,z) is considered as the aggregation of the product of the response surface factor and the mode data for all of the modes r:

$$v\left(\mathbf{x}\right)=\sum _{i=1}^{r}GAR{S}_i\left(\mathbf{x}\right){mode}_i\left(\mathbf{x}\right).$$

(8)

Completing the above-described ROM creation procedure can require hours or even days.

Finally, the ROM consumption section involves the use of solvers and standalone systems. Having acquired the GARS of all the modes, the calculation of the solution of each new set of input parameters is computed almost in real-time. The building of a ROM is considered successful when it estimates the flow parameters of input combinations that are not included in the training set with adequate precision. The quality of the built ROM is highly associated with the number of training scenarios, the number of selected modes, and the accuracy of the snapshots, meaning the convergence of the fluid dynamic problem of interest.

3. Materials and Methods

3.1. Patient-Specific Data Pre-Processing

For the purpose of this study, an infant case affected by complete pulmonary atresia who underwent the MBTS implant operation was selected from the clinical database of the Heart Hospital (Massa, Italy) of Fondazione Toscana Gabriele Monasterio. The research activity was approved by the local institutional review board. In accordance with the ethical guidelines of Declaration of Helsinki (1975), informed consent for personal data management was obtained for the case. The patient-specific data were taken from a Computed Tomography (CT) dataset. In particular, a conventional clinical CT 640-detector scanner (Toshiba Aquilion One, Toshiba, Japan) was employed. The 3D patient-specific model was reconstructed with a semi-automated segmentation algorithm (Figure 2a,b) [39]. The dataset was segmented to obtain the systemic vessels, including the ascending/descending aorta sections, the aortic arch and the epiaortic vessels, and the pulmonary branch, including the left and the right pulmonary arteries and the pulmonary valve plane affected by atresia.

Using the discretised surface model available after segmentation (i.e., STL file), the Computer-Aided Design (CAD) model of the vessels of interest was generated by reconstructing the shapes. This operation was accomplished using the standard features of the ANSYS^® SpaceClaim^® software. Exploiting the advanced features of this software, a semi-automated procedure was implemented to create the CAD model of MBTS. In order to set up this parametric-based procedure, the characteristic geometrical quantities of the MBTS were assigned. These tools allowed us to create the geometry of the MBTS with a high level of flexibility. The obtained aortic CAD model with the different components highlighted is shown in Figure 3.

3.2. Mesh Morphing Set-Up

The RBF Morph (www.rbf-morph.com) software program was employed to mold the configuration of the shunt, using the previous baseline as a starting point. The CAD baseline geometry from the previous step was meshed using the meshing tool of ANSYS^® Fluent^® in order to obtain an initial computational mesh to be morphed. Source points organized into four sets were selected on the implant’s district, as depicted in Figure 4a,b. The morphing action was constrained inside the borders of the blue cylinder indicated in Figure 4a. The four different sets of source points were generated and placed at:

• the inferior boundary of the shunt’s geometry (ISP)
• the central segment of the shunt, corresponding
  to the maximum cross-sectional diameter (DSP)
• the cylindrical periphery of the shunt’s geometry (CSP)
• the superior boundary of the shunt’s geometry (SSP)

Adhering to the clinicians’ needs, twelve RBF transformations were introduced to cover the wide spectrum of possible shunt configurations. The shunt’s geometry was modified by imposing linear translation and rotation to its upper, central and lower section. The combination of all the RBF modifiers enables the creation of specific shunt configurations. For their determination, a local coordinate reference system (LCRS) was defined. The LCRS was placed on the shunt’s cross-section centerline (Figure 4b).

The center line of the implant was defined according to [23]. The setup of each modifier was built by activating a group of source points subjected to a predefined motion by following the description reported below:

• dl-1-vol—Rigid motion of ISP along the ±x direction: sliding of the shunt’s root along the length of the pulmonary artery (Figure 5a).
• dl-2-vol—Rigid motion of SSP along the ±x direction: sliding of the shunt’s top segment along the length of the right subclavian artery (Figure 5b).
• dr-1-vol—Rigid motion of ISP along the ±y direction: sliding of the shunt’s root along the width of the pulmonary artery (Figure 5c).
• dr-2-vol—Rigid motion of SSP along the ±y direction: sliding of the shunt’s top segment along the width of the right subclavian artery (Figure 5d).
• mid-dl-vol2—Rigid motion of DSP along ±x direction: inflation or deflation towards the ±x axis (Figure 5e).
• mid-dl-vol—Rigid motion of CSP along ±x direction: inflation or deflation towards the ±x axis (Figure 5f).
• mid-dr-vol2—Rigid motion of DSP along ±y direction: inflation or deflation towards the ±y axis (Figure 5g).
• mid-dr-vol—Rigid motion of CSP along ±y direction: inflation or deflation towards the ±y axis (Figure 5h).
• rl-1-vol—Rotation of ISP with respect to the y axis of LCRS: the shunt’s root, which is located on pulmonary boundary, is revealed or hidden on the zx plane (Figure 5i).
• rl-2-vol—Rotation of SSP with respect to the y axis of LCRS: the shunt’s upper segment which is located on the right subclavian aortic boundary, is revealed or hidden on the zx plane (Figure 5j).
• rr-1-vol—Rotation of ISP with respect to the x axis of LCRS: the shunt’s root is revealed or hidden on the yx plane (Figure 5k).
• rr-2-vol—Rotation of SSP with respect to the x axis of LCRS: the shunt’s upper segment is revealed or hidden on the yx plane (Figure 5l).

The shape factors were employed by using amplification factors to determine the extent of their influence on the morphology of thw MBTS. The amplification factor of the dl-1-vol, dl-2-vol, mid-dl-vol, mid-dl-vol2, mid-dr-vol, and mid-dr-vol2 ranged between −0.5 and 0.5, whereas the range of the rest of the shape parameters was modified between −5 and 5. The action of each shape parameter is displayed in Figure 5, where the maximum and minimum effect of each modifier is displayed.

For every scenario, the mesh was modified successively according to the amplification factor that was imposed for each shape modifier. This sequence of mesh morphing actions induced the risk of mesh quality deterioration. In order to control this hazard, the maximum mesh cell squish of each scenario was monitored throughout the DOE study.

3.3. CFD Set-Up

ANSYS^® Fluent^® - Release 21R1 was used to solve the governing flow equations. A steady state flow regime was adopted and the k-$\omega$ SST model was employed to integrate the turbulence [40,41]. An unstructured mesh consisting of 1.79 million polyhedral elements was built. In order to capture the laminar ongoing phenomena close to the boundaries, ten inflation layers with a growth rate of 1.05 and a total thickness of 2.5 mm were introduced. The blood was considered as a Newtonian fluid with a density and viscosity equal to 1060 kg/m ${}^3$ and 3.5 $·{10}^{-3}$ Pa s, respectively. The systolic peak flow data were used as the inlet boundary condition by applying a velocity condition equal to 0.24 m/s at the aortic valve section. For the aortic branches, the Right/Left Common Carotid Arteries (RCCA/LCCA), Left Subclavian Artery (LSA), and descending aorta (DA) were handled with a constant pressure condition of 55 mmHg. The Right Subclavian Artery (RSA) was treated with a 48 mmHg condition. For the pulmonary branch, the right and the left pulmonary artery (RPA/LPA) boundary conditions were set to a constant pressure condition of 8 mmHg and 7 mmHg, respectively. These pressure values correspond to the physiological pressure conditions of an infant. A wall condition was imposed on the pulmonary valve boundary to simulate the pulmonary atresia condition.

3.4. ROM Set-Up

The DOE definition for the ROM was performed on ANSYS^® DesignXplorer^®-Release 22R1. In order to ensure the optimal spatial allocation of the investigated scenarios, an Optimal Space-Filling algorithm was employed. Starting from the baseline configuration, a total of 150 runs were performed by combining the shape modifiers described in Section 3.2. The associated snapshots were imported in ROM Builder (ANSYS^® Electronics Desktop ${}^{\mathrm{TM}}$ - Release 22R1) and used to feed the SVD algorithm. The number of modes chosen for the Digital Twin was 27 for the pressure, 18 for the velocity, and 22 for the wall shear stress (WSS) evaluation. The decision regarding the number of selected modes was carried out with the aim of reaching the best compromise between the accuracy and responsiveness of ROM. For evaluation of the accuracy of the generated ROM, the original full CFD and the reduced computational models were compared. The error between the two models was calculated according to Equation (9)

$$Error=\parallel CF{D}_{solution}-RO{M}_{prediction}\parallel .$$

(9)

In order to simulate the application of the workflow, 30 snapshots subtracted from the available dataset were used for verification of the ROM. For these cases, the evaluation was conducted in terms of velocity, pressure, and WSS maps. The entire computational part was performed on a workstation equipped with two Intel Xeon 3.00 GHz processors, exploiting 14 cores from each of them and 256 GB of shared RAM.

4. Results

In this section, the results of the adopted workflow are presented. These results are organized into three subsections: (i) verification of the RBF mesh morphing; (ii) validation of ROM; and (iii) hands-on demonstration of and evaluation of ROM.

4.1. Mesh Morphing Verification

Concerning the performance of the mesh morphing verification,

Table 1. Four worst-case scenarios presenting the greatest mesh quality degradation. The amplification shape factors close to the range boundaries are shown in bold font.

Shape factor	Scenario 23	Scenario 43	Scenario 123	Scenario 142
dl-1-vol	−0.26	−0.49	−0.36	−0.39
dl-2-vol	−0.03	−0.27	0.48	−0.46
dr-1-vol	−0.17	3.30	0.70	−2.23
dr-2-vol	−3.30	0.90	4.23	−0.70
mid-dl-vol2	−0.04	0.30	−0.08	0.34
mid-dl-vol	−0.02	0.18	−0.18	0.14
mid-dr-vol2	0.48	−0.36	-0.08	0.34
mid-dr-vol	−0.25	0.24	−0.34	0.05
rl-1-vol	−4.63	4.43	−1.83	3.30
rl-2-vol	−2.83	0.57	3.90	2.43
rr-1-vol	−0.57	2.50	−1.37	−3.10
rr-2-vol	4.77	−0.37	3.17	−2.37
Cell Squish	0.982	0.994	0.906	0.973

depicts the setup of the scenarios linked with the highest mesh degradation. The highest cell squish values and corresponding amplification factors of the shape parameters are reported. The shape amplification factors close to the range boundaries are shown in bold font within the Table. As presented, the mesh quality is preserved within the studied range of possible shunt configurations, reaching the peak cell squish index of 0.994 on three mesh elements in scenario 43. It is worth mentioning that although the initial mesh process took about ten minutes, the morphing procedure for each scenario required only one minute. The display of the morphed mesh preview was executed almost in real-time (i.e., less than 1 s).

4.2. ROM Verification

Figure 6 illustrates the velocity flow field evaluated at the MBTS middle cross-section for the following selected scenarios: No. 135 (Figure 6a–e), No.1 (Figure 6b–f), No. 34 (Figure 6c–g), and No. 146 (Figure 6d–h). The full CFD cases and the corresponding ROM solutions are displayed on the upper and the lower level of the Figure, respectively. The velocity estimation of the ROM adequately approaches the accurate CFD solution, achieving an average error of 0.03 m/s. The highest discrepancies are marked locally on the segment of the pulmonary artery (Figure 6c–g,d–h).

In Figure 7 surface pressure ROM is compared with the CFD model for four different shunt morphologies. The local maps for scenarios No. 76 (Figure 7a–e), No. 87 (Figure 7b–f), No. 99 (Figure 7c–g), and No. 101 (Figure 7d–h) are represented at the MBTS level. As can be seen, the ROM approach captures the pressure distribution along the shunt’s geometry. The ROM’s average pressure error is 1.92 mmHg.

Similarly, in Figure 8 the WSS distribution between the CFD and the ROM model is displayed for four different investigated scenarios: No. 65 (Figure 8a–e), No. 52 (Figure 8b–f), No. 63 (Figure 8c–g), and No. 87 (Figure 8d–h), proving that the latter achieves precise predictions. The average deviation between the WSS ROM estimation and the corresponding CFD solution is 3.83 Pa.

In order to quantify the accuracy of the ROM models, the maximum detected absolute errors over all the investigated scenarios in terms of velocity, pressure, and WSS are depicted on Figure 9. The greatest velocity, pressure, and WSS errors occur in scenarios No. 34, No. 101, and No. 87, respectively.

4.3. ROM Consumption

Having successfully generated the ROM model, the results regarding the ROM consumption are presented in this subsection. The built velocity field ROM file occupies a memory storage of 3 GB, considering its volumetric distribution, while those relative to pressure and WSS fields require 31.2 MB and 35.6 MB, respectively. Concerning the computational requirements, the ROM loading requires 2–3 min. After their mounting is completed, the user is allowed to interactively set up the MBTS geometry of interest by tuning the available sliders which correspond to the amplification factor of the RBF shape modifiers.

A demonstration of the resulting Digital Twin with its interactive environment for estimation of the velocity field is depicted in Figure 10. Additional examples are provided in the Supplementary Materials section. Moreover, a timeframe comparison for the time performance assessment of the ROM approach against the full CFD is reported in Figure 11.

5. Discussion

In the literature, several studies have been presented to investigate the effect of the implant geometry of the MBTS on the cardiovascular hemodynamic fluid flow [7,42,43]. Although most of these agree on the pivotal importance of the shunt’s morphology for the success of an MBTS operation, no procedure has yet been established to perform preoperative patient-specific studies by exploring multiple geometrical configurations. A complete morphological analysis reaching beyond the simple consideration of the implant’s diameter and counting in the unique anatomy of each patient is highly required. In order to meet this need, a comprehensive framework which fuses the ROM approach and the RBF mesh morphing technique is presented in this paper, with the added value of providing fast and interactive prediction of the hemodynamic parameters. An infant case affected by complete pulmonary atresia with an MBTS implant was selected for demonstration of the pipeline’s soundness. Mesh morphing enabled the modification of the position, size, and geometry of the implanted shunt by means of RBF shape modifiers. These modifiers were used to generate possible shunt configurations, i.e., scenarios, for the building up of the ROM. Both morphing and ROM analysis were verified in terms of mesh degeneration and full CFD comparison, respectively. Finally, the ROM consumption allowed the definition of an effective Digital Twin for the evaluation of hemodynamic flow variables in real-time. The developed tool was validated correctly by demonstrating the ROM accuracy. The current work constitutes a significant starting point for a future decisional tool that can be used in actual clinical contexts.

As mentioned in Section 3.2, it is crucial to inspect and control the mesh quality when employing a mesh morphing method. Thereby, the adopted RBF mesh morphing technique was checked. The main outcome flowing from this verification is that mesh quality degradation was demonstrated to be negligible for almost all of the explored scenarios. As reported in Table 1, which illustrates the four scenarios linked with the greatest mesh degradation, the mesh quality is considered acceptable because the maximum cell squish of the mesh remains lower than the critical value of 1. In the reported scenarios, the amplification factors of the rr-2-vol, rl-1-vol, dl-1-vol, dl-2-vol, and mid-dr-vo2l shape parameters lie close to the range extremes. This suggests that these shape factors had the most influence on the meshing morphing performance of the MBTS structure. Even though this study is, to the best of our knowledge, the first to report the adoption of morphing techniques on MBTS implants, this approach is in accordance with previous computational analyses in the cardiovascular context [44].

Moving on to the building of the ROM, the high accuracy of the ROM predictions were verified with success. More precisely, validation was carried out with regard to the velocity (Figure 6), pressure (Figure 7), and WSS (Figure 8).

The scenarios reported in Figure 6 show a velocity range and distribution in line with those reported in the previous literature [45,46]. It appears that higher velocity discrepancies occur at the MBTS junction with the pulmonary branch. This is particularly evident in scenarios with low-curvature MBTS implants (Figure 6c–g,d–h). This is further confirmed by the error map of Figure 9a, where maximum absolute errors of 1 m/s (about 25% of the full velocity range) are reached in the pulmonary branch zone for snapshot No. 34. However, it is important to highlight that inaccuracies remain highly localized on several mesh cells and involve a limited number of scenarios. Additionally, the area of interest, the MBTS, continues to maintain a lower level of errors. With respect to pressure and WSS, it can be confirmed that the ROM models correctly represent the full CFD flow pattern even in the worst-case scenarios. Figure 9b refers to the highest recorded error of the pressure ROM, which was identified in snapshot No. 101. The greatest pressure error (≈1.5 mmHg) is in the middle part of the shunt. Finally, the prediction of WSS was characterised by localised discrepancies (≈48 Pa) on the upper partition of the implant, as illustrated in snapshot No. 87 (Figure 9c). It is important to underline the fact that these errors constitute about 2.5% and 6% of the full range of values of pressure and WSS, respectively.

Finally, the ROM consumption allowed for the generation of a Digital Twin, as shown in the example in Figure 10. In order to appreciate the exhibited results, the timeframe needed to perform MBTS morphology exploration with the CFD and ROM approaches is presented. It is of paramount importance to point out that the solution of the complete model for every new MBTS configuration requires 45 min on an 8-core laptop. On the other hand, the ROM build-up takes 20 h and 10 min on a 28-core workstation and only has to be calculated once. Afterwards, an estimation of the hemodynamic flow as a function of the selected shape factors of the MBTS implant (ROM consumption part) is provided in several seconds (near real-time) (Figure 11). Moreover, the reduced memory space required by the Digital Twin represents an additional positive feature. Consequently, the ROM can be considered as a rather competitive surrogate model that has the potential to be inserted into demanding clinical practice contexts.

In synthesis, it is worth highlighting the weak and strong points of the presented work pipeline. The weak points to be improved include: (i) the possibiliy of increasing the number of SVD modes/ROM snapshots in order to further reduce the discrepancy; (ii) the necessity of including cardiac motion within the computational model; (iii) the lack of time-varying fluid dynamic conditions; and (iv) the need for future clincal validation. On the other hand, it is important to underline a series of strong points. In particular, the obtained accuracy of the pipeline was acceptable in the case tested here, and our approach proved powerful enough to be adopted for the exploration of different MBTS configurations. In addition, the execution times and required memory space needed to obtain a fast and efficient analysis were very satisfactory.

6. Conclusions

The creation of a patient-specific Digital Twin for preoperative inspection of cardiovascular operations is on the cutting edge of today’s research interests. New multidisciplinary computational methods are currently under development for the prediction of the fluid flow of several pathologies, such as cyanotic heart diseases. Espousing this tendency, a computational framework is proposed within this paper for the presurgical inspection of an MBTS operation. The realization of a fast and interactive pipeline for the exploration of the effect of the shunt’s geometry on the fluid flow was performed employing the reduced order model and radial basis function mesh morphing techniques. The methods and procedure we followed are discussed in this work, and the results appear rather promising. The reported workflow aims to move a step further in the world of cardiovascular pediatrics by contributing to a deeper understanding of the correlation between MBTS geometry and the hemodynamic flow parameters. This study has certain limitations in terms of the accuracy, the adopted assumptions and the feasibility of automating the exhibited pipeline which ought to be addressed. For this reason, possible future developments are suggested. To begin with, the local errors detected on few scenarios between the reduced and full order models can be further minimised by increasing the number of SVD modes and ROM snapshots, or alternatively by reducing the range of exploration of the different shunt geometries. Additionally, the adoption of time-variant inlet flow conditions and lumped-parameter boundary conditions could bring us one step closer to dynamic visualization of the full flow field during the heart cycle. In the current study, we demonstrate has been proven that the way has been paved for computational patient-specific MBTS preoperative planning. It is worth noting that the current study was not aimed at individuating the best MBTS configuration in terms of clinical outcomes, and it requires clinical evaluation in order to become an actual decisional tool. Nevertheless, expanding the strategy of this study for the numerical simulation of more complex flow conditions appears to be an interesting upcoming challenge, with potential impacts on the clinical reality.