Effects of Patent Ductus Arteriosus on the Hemodynamics of Modified Blalock–Taussig Shunt Based on Patient-Specific Simulation

The question of preserving the patent ductus arteriosus (PDA) during the modified Blalock–Taussig shunt (MBTS) procedure remains controversial. The goal of this study was to investigate the effects of the PDA on the flow features of the MBTS to help with preoperative surgery design and postoperative prediction. In this study, a patient with pulmonary atresia and PDA was included. A patient-specific three-dimensional model was reconstructed, and virtual surgeries of shunt insertion and ductus ligation were performed using computer-aided design. Computational fluid dynamics was utilized to analyze the hemodynamic parameters of varied models based on the patient-specific anatomy and physiological data. The preservation of the PDA competitively reduced the shunt flow but increased total pulmonary perfusion. The shunt flow and ductal flow collided, causing significant and complicated turbulence in the pulmonary artery where low wall shear stress, high oscillatory shear index, and high relative residence time were distributed. The highest energy loss was found when the PDA was preserved. The preservation of PDA is not recommended during MBTS procedures because it negatively influences hemodynamics. This may lead to pulmonary overperfusion, inadequate systemic perfusion, and a heavier cardiac burden, thus increasing the risk of heart failure. Also, it seems to bring no benefit in terms of reducing the risk for thrombosis.


INTRODUCTION
Pulmonary atresia is a complex congenital heart defect characterized by flow interruption between the right ventricle and pulmonary artery. The ductus arteriosus is the fetal blood vessel between the aorta and pulmonary artery, and it can provide pulmonary blood flow postnatally if kept open. Pulmonary perfusion is dependent on the patent ductus arteriosus (PDA) in some patients with pulmonary atresia (Reddy and Hanley, 2012). However, the pulmonary perfusion can still be insufficient due to limited ductal flow. Thus, the modified Blalock-Taussig shunt (MBTS) is required to increase pulmonary blood flow. A Gore-Tex conduit is surgically interposed between the innominate artery (or subclavian artery) and the pulmonary artery to direct blood flow from the systemic circulation to pulmonary circulation (de Leval et al., 1981).
The question of whether to preserve the PDA when performing MBTS procedures in patients with ductal-dependent pulmonary blood flow is still controversial. Surgeons usually make decisions according to experience because there is a lack of studies to ascertain which procedure is better. One reason to preserve the PDA is that it can provide surgeons with additional time to do shunt revision and save lives when the shunt fails due to shunt thrombosis or obstruction. However, the coexisting PDA may be a dangerous flow source for the pulmonary circulation, resulting in pulmonary overperfusion and diastolic runoff of systemic circulation (El-Rassi et al., 2012). These two sources increase the difficulty in regulating the pulmonary blood flow during the perioperative period (Dirks et al., 2013).
Computational fluid dynamics (CFD) and computeraided design (CAD) enable quantitative investigations of hemodynamic characteristics to find the optimal surgery design Celestin et al., 2015;Piskin et al., 2017). Esmaily-Moghadam et al. (2015) utilized the idealized vascular models to explore the hemodynamics of PDA and MBTS and concluded that MBTS with the preservation of PDA did not bring benefits in oxygen delivery but increased the risk of thrombosis. This research is limited in clinical practice due to not considering the patient-specific anatomy. Zhang et al. (2020) also investigated the effects of PDA on MBTS based on the medical images of the patient and the virtually designed PDA. However, the study by Zhang did not give a definite conclusion about the preservation or ligation of PDA but emphasized the importance of considering individual patient conditions in the surgical decision-making process. This study sought to analyze the hemodynamic characteristics of MBTS (with and without the PDA) based on the patient-specific anatomy and physiological conditions. The flow distribution, velocity streamlines, wall shear stress (WSS), oscillatory shear index (OSI), relative residence time (RRT), and energy loss (EL) were analyzed using CFD. This was performed to understand the flow features and to aid in preoperational surgery design and prediction of postoperative prognosis.

Clinical Data Acquisition
For this study, a 2-month-old female child born with pulmonary atresia and PDA was included. Informed consent was obtained from the guardians of the patient, and the local institutional review board together with the regional research ethics committee of the Shanghai Children's Medical Center (SCMC) affiliated to the Shanghai Jiao Tong University School of Medicine approved all associated studies. Patient-specific medical images were acquired by a 64-slice spiral CT scanner (GE Discovery CT750 HD, USA) before surgery. The slice thickness was 0.625 mm, and each view was of a 512 × 512 pixel field. The echocardiographic data were obtained simultaneously.

Model Reconstruction and Virtual Surgery
The three-dimensional (3D) vascular model (Model 1) was reconstructed by Materialise R Mimics Innovation Suite 20.0 (Materialise NV, Leuven, Belgium) based on the patient-specific CT images. The PDA has an irregular shape, as its diameter reduces between the end of the aorta and the end of the pulmonary artery. There is an expanded lumen in the confluent site of the left pulmonary artery (LPA) and right pulmonary artery (RPA). CAD was utilized to perform virtual surgeries.
A 4-mm conduit was selected to connect the innominate artery (IA) and RPA according to the suggestions of the surgeon. The newly created model was named Model 2, and the PDA of Model 2 was preserved. Model 3 was created based on Model 2, with the PDA being wholly removed to mimic the surgical ligation of PDA. Figure 1 displays the three models in detail.

Governing Equations
The blood was regarded as a Newtonian fluid, and the 3D incompressible Navier-Stokes equation was solved to simulate the blood flow of the shunt domain (Qian et al., 2010), shown as follows: where i, j = 1, 2, 3, x 1 , x 2 , and x 3 represent coordinate axes; u i and u j are velocity vectors; p is pressure; µ is viscosity; ρ is density; and t is time. The term f i means the action of body forces. The standard k-ε turbulence model was applied (Qian et al., 2010).

Mesh Generation
Three models were discretized by ANSYS R -ICEM CFD 2020 (ANSYS Inc., USA) for calculation. The tetrahedral mesh was fabricated in the interior, and five layers of prism mesh were generated along the vascular wall in the near-wall region to make calculated parameters more accurate. The inlet was extruded with a length 20 times the inlet diameter to make the inflow fully developed. Each outlet was lengthened 60 times the outlet diameter to simulate the peripheral vasculature for developing antegrade flow. Grid-independent verification was performed using Model 1, and the EL was set as the indicator with the number of grids changing. The results indicated that a grid number of around 1 million would generate the most efficient mesh. Grid numbers of all models ( Table 1) were more than 1 million, which were sufficient to produce robust calculations in this study.

Boundary Conditions
The ascending aorta (AAO) was the only inlet, and the inflow was the pulsatile velocity of AAO measured by the real-time  echocardiography ( Figure 2). A reference pressure of 50 mmHg (arterial blood pressure at end-diastole) was imposed on the aorta domain. The outlets included the descending aorta (DA), right subclavian artery (RSA), right common carotid artery (RCCA), left common carotid artery (LCCA), left subclavian artery (LSA), LPA, and RPA. Previous literature suggested that the measured pressure wave is composed of the forward pressure wave and the backward pressure wave, to explain the pressure wave propagation (Parker and Jones, 1990;Khir et al., 2001;Trachet et al., 2010). The forward pressure wave is derived from the inflow ejected by ventricle contraction. The backward pressure wave is the reflected wave formed by the downstream artery branches. To make the simulation close to the physiology, the backward pressure wave was considered and applied in DA and other aortic branches (Liu et al., 2011(Liu et al., , 2015, as is shown in Figure 2. Because the diameters of aortic branched vessels were similar, the backward pressure waves of four outlets (RCCA, RSA, LCCA, and LSA) were assumed the same. A constant pressure of 10 mmHg was simultaneously utilized in the LPA and RPA (Sundareswaran et al., 2006;Sun et al., 2014). Rigid and impermeable vascular walls with no slip were assumed. The same boundary conditions were set in the three models to compare the hemodynamic characteristics of the preservation of PDA or ligation during MBTS procedures.
Simulations ANSYS R -Fluent 2020 (ANSYS Inc., USA) was used to finish the transient simulations of three models. The semi-implicit (SIMPLE) method and the second-order upwind scheme were applied. Because the blood flow in large arteries is pulsatile, high-speed, and high-shear and the diameter of the aorta is much larger than the diameter of a single red blood cell, the blood could be simplified as a Newtonian fluid with a constant density of 1,060 kg·m −3 and viscosity of 4.0 × 10 −3 Pa·s (Liu et al., 2015;Ceballos et al., 2019). The time step was 10 −5 s, and the convergence criterion was 10 −5 . Three cardiac cycles were calculated, and the results of the third cycle were analyzed. The computational analysis system (currently used) was described in detail and validated in previous studies (Qian et al., 2010;Liu et al., 2011Liu et al., , 2015.

Hemodynamic Parameters
The flow rates, pulmonary/systemic ratio (Q P /Q S ), and the LPA/RPA ratio (Q LPA /Q RPA ) were calculated to evaluate the flow interactions and distribution. The WSS is considered to be related to vascular endothelial injuries and thrombosis (Koskinas et al., 2012;Kumar et al., 2018), calculated according to the equation as follows: where u x is the velocity of fluid near the vessel wall and y is the distance to the vessel wall. The time-averaged WSS (TAWSS) is defined as follows: where T is the time of a cardiac cycle, equal to 0.5 s. The OSI is a non-dimensional metric to describe the WSS vector deflection from the predominant direction of blood flow FIGURE 2 | The computational system and boundary conditions. The purple arrows represent the direction of inflow, and the black arrows indicate the direction of outflow. The inlet pulsatile velocity and the backward pressure waves in different outlets are shown at the bottom. PDA, patent ductus arteriosus; MBTS, modified Blalock-Taussig shunt; AAO, ascending aorta; DA, descending aorta; RSA, right subclavian artery; RCCA, right common carotid artery; RPA, right pulmonary artery; LCCA, left common carotid artery; LSA, left subclavian artery; LPA, left pulmonary artery. over a cardiac cycle (Ku et al., 1985;He and Ku, 1996) and is defined as follows: The RRT is the residence time that the particles spent around the endothelium, introduced as inversely proportional to OSI and TAWSS (Himburg et al., 2004), The EL is a quantitative indicator to evaluate hemodynamic efficiency and cardiac workload, given by the equation as follows: (7) where V and Q are velocity and flow rate, respectively, in indicates the inlet, and out indicates the outlet.

Flow Distribution
Blood flow from the right ventricle to the pulmonary artery was completely obstructed in the subject of the study. Therefore, pulmonary perfusion was provided by the flow through the MBTS and/or the PDA. Model 2 with coexisting MBTS and PDA was associated with the highest pulmonary perfusion and Q P /Q S , as listed in Table 2. The pulmonary perfusion of Model 1 was the lowest. The shunt flow of Model 2 was lower than that of Model 3, and the ductal flow of Model 2 was lower than that of Model 1. It turns out that competitive flow was present when the PDA was preserved. The Q LPA /Q RPA of Model 2 was located between Model 1 and Model 3, and the difference of Q LPA /Q RPA was insignificant among the three models. The flow distribution between the LPA and RPA was asymmetric and tended to RPA in all models.

Velocity Streamlines
The velocity streamlines at the systolic peak are shown in Figure 3. Shunt flow was associated with high velocity no matter whether the PDA was ligated or not. This shunt flow ran through the vertical conduit and impacted the pulmonary artery wall, causing turbulence in Models 2 and 3. More turbulence was formed when the shunt flow and ductal flow rushed against each other in the pulmonary artery between the MBTS and PDA. Fewer streamlines through the PDA went into the RPA in Model 2 compared with Model 1. Similarly, fewer streamlines through the MBTS went into the LPA in Model 2 compared with Model 3. These results suggested that the flow split between the LPA and RPA was affected by flow interactions between the shunt and ductus.

Wall Shear Stress
The WSS distribution at the systolic peak is shown in Figure 4. The shunts of Model 2 and Model 3 were associated with high WSS. The pulmonary anastomotic site of the shunt was found with sharp WSS distribution compared with other sites of the pulmonary artery. High WSS was also found in the confluent site of the PDA and pulmonary artery, especially evident in Model 1. The TAWSS was also higher in the shunt and the anastomosed site of the pulmonary artery ( Figure 5). In the pulmonary artery between the MBTS and PDA, the TAWSS was notably lower in Model 2.

Oscillatory Shear Index
The OSI value ranged from 0 to 0.5, i.e., 0 indicated noncyclic variation and 0.5 indicated 180.0 • deflection from the WSS direction. A high OSI was found in the shunt of Models 2 and 3 (Figure 6). The area of high OSI distributed in the shunt of Model 2 was larger. The pulmonary artery between the MBTS and PDA was also associated with higher OSI. The maximum OSI reached 0.49, which was more profound in Model 2.

Relative Residence Time
The RRT distribution varied among the three models. The inner side of the shunt was associated with high RRT (Figure 7). High RRT mainly appeared downstream of the pulmonary artery in Model 1. As for Model 2, the pulmonary artery between the MBTS and PDA was associated with high RRT. The area of high RRT distribution was smaller in Model 3 compared with Model 2.

Energy Loss
The EL is used to evaluate energy efficiency and ventricular workload quantitatively. Figure 8 shows the time-averaged EL through a cardiac cycle. The EL of Model 2 was the highest, and the preservation of PDA increased energy dissipation of flow and cardiac workload.

DISCUSSION
In the patient with pulmonary atresia, the pulmonary blood flow comes from the systemic circulation through the PDA and/or surgically created MBTS. Based on the patient-specific anatomy and physiological data, the simulation results demonstrated the existence of flow competition between the MBTS and PDA and the influences of flow interactions. This provided quantitative  and qualitative hemodynamic parameters to help with surgery design and postoperative prediction. Balanced flow distribution between the pulmonary circulation and systemic circulation is intimately related to the prognosis of MBTS, which can be quantitatively evaluated by Q P /Q S . Excessive pulmonary blood flow may lead to pulmonary overperfusion and diastolic runoff of the systemic circulation. Insufficient pulmonary perfusion would cause inadequate oxygen delivery, thrombosis, shunt occlusion, and an underdeveloped pulmonary artery (Fenton et al., 2003;Dirks et al., 2013). The ideal value of Q P /Q S is slightly <1, which helps provide adequate pulmonary blood flow and oxygen delivery without reducing systemic perfusion and myocardial blood flow (Barnea et al., 1994(Barnea et al., , 1998. However, if the value of Q P /Q S is <0.7, a further decrease of Q P /Q S reduces the oxygen delivery and negatively affects hemodynamics (Riordan et al., 1996). Since the same boundary conditions were applied in simulations, the values of Q P /Q S were mainly affected by the anatomy. The Q P /Q S of Model 1 was <0.7, suggesting that the PDA alone cannot provide enough pulmonary perfusion, and this was the reason for cyanosis and corresponded with the diagnosis of the patient. The Q P /Q S of Model 3 was between 0.7 and 1, providing sufficient pulmonary perfusion without resulting in insufficient systemic blood flow. The Q P /Q S of Model 2 was more than 1,  increasing the risk of pulmonary overperfusion and congestive heart failure.
Balanced flow distribution between the LPA and RPA benefits symmetric pulmonary artery growth in patients with bilateral hypoplastic pulmonary arteries. Coexisting MBTS and PDA produced a value of Q LPA /Q RPA located between Model 1 and Model 3. The Q LPA /Q RPA differences were insignificant among the three models, and the pulmonary flow favored RPA in all models. The unimportant differences could be explained by the patient-specific morphology of the pulmonary artery because the boundary conditions of LPA and RPA were assumed to be similar in this study. The RPA is associated with a larger diameter, and the confluent site of LPA and RPA has an expanded lumen. The swirling flow was developed in the confluent pulmonary artery despite the PDA ligation or MBTS insertion, which may be the reason for the decreased left pulmonary perfusion. These results are not in line with the study by Zhang using the ideal vascular model (Zhang et al., 2020). This study reflected the concerning influence of the individualized pulmonary vascular geometry on the flow field and shed light on the importance of numerical simulation based on the patient-specific anatomy.
The hemodynamic EL is an effective indicator caused by blood movement, flow collision, and turbulence to evaluate energy efficiency and ventricular workload quantitatively. The MBTS with PDA was associated with the largest EL, which could be explained by augmented ventricular volume load, flow interactions, and significant turbulence. As time accumulates, coexisting MBTS and PDA will gradually increase cardiac workload. The EL of Model 3 was higher than Model 1, suggesting that the MBTS increased the cardiac burden. Adequate coronary perfusion is critical for heart pump with an augmented cardiac burden. Given that coronary perfusion decreases as pulmonary perfusion increases, the preservation of PDA may augment the risk of heart failure in neonates with hypoplastic myocardium.
Shunt-related complications such as thrombosis and occlusion remain unresolved. Very low or very high WSS is likely to cause endothelial injury and platelet activation, increasing the risk of thrombosis and occlusion (Holme et al., 1997;Kumar et al., 2018). Low WSS and high OSI are associated with endothelial injury and atherosclerosis, and OSI > 0.2 is a risk factor for atherosclerosis (Frydrychowicz et al., 2009). RRT is used to assess the time that the particles spend in the endothelium, and high RRT indicates a higher possibility of particle deposition. The outer sides of the shunt and anastomosed sites in the pulmonary artery were associated with higher WSS. Low WSS and high OSI were developed in the inner side of the shunt. Competitive flow also yielded a region with lower WSS, higher OSI, and RRT inside the pulmonary artery. These results suggest that the preservation of PDA brings no benefits in potential risks of shunt thrombosis and occlusion, which still need to be carefully monitored during the perioperative period and the long-term follow-up.
There were several limitations in this study. First, the assumption of a rigid wall without considering the elasticity and compliance of the great vessels. Although fluid-structure interaction can be applied to consider the elasticity, the rigid wall can achieve comparative results with minor errors (Esmaily-Moghadam et al., 2015). The pressure wave reflection was set in all outlets to reduce the influence of the rigid wall assumption, which was validated in previous studies (Liu et al., 2011(Liu et al., , 2015. Second, only one type of shunt design was discussed in this study. Shunt diameter, shunt position, anastomosis angle, and shunt shape will affect the flow features (Piskin et al., 2017). Nevertheless, the shunt design in this study was appropriate to be applied. According to the advice of the surgeon, a shunt size of 4 mm was selected, and the MBTS was placed in the contralateral side of the PDA. The anastomotic sites were decided based on the spatial position and geometry of the aorta and pulmonary artery branches. Third, the position and size of PDA affect the flow features of MBTS. This study preliminarily investigated the effects of the existence of PDA on the shunt flow. A further study focusing on the effects of the position and morphology of PDA on hemodynamics will be added in the future.

CONCLUSION
The MBTS with PDA is accompanied by more pulmonary perfusion, flow competition, large-scale and complicated turbulence, and higher EL. The ligation of PDA is recommended during the MBTS procedure since the preservation of PDA is riskier, and it may lead to pulmonary overperfusion, insufficient systemic blood flow, and a heavier cardiac burden. Low WSS, high OSI, and high RRT are found in the pulmonary artery between the MBTS and PDA when the PDA is preserved, which does not reduce the risk of thrombosis. The risks of shunt thrombosis and occlusion are still worth noting in the perioperative period and the long-term follow-up. Vascular morphology plays a vital role in hemodynamics, which should be considered in surgery design.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
Written informed consent was obtained from the guardians of the patient. All associated studies were approved by the local institutional review board and regional research ethics committee of Shanghai Children's Medical Center (SCMC) affiliated to Shanghai Jiao Tong University School of Medicine.

AUTHOR CONTRIBUTIONS
JX was responsible for numerical simulation, data analysis, and manuscript preparation. QS was responsible for the clinical data acquirement and virtual surgery. YQ developed the idea of the study and assisted in the manuscript review. LH assisted in acquiring clinical images. ZT assisted in the reconstruction of the 3D vascular model. JinfL was responsible for data analysis and manuscript review. JinlL conceived the study and helped with manuscript revision. All authors contributed to the manuscript and approved the submitted version.