Computational fluid dynamics-based risk stratification of modified Blalock-Taussig-Thomas shunt thrombogenicity

Study Overview and Model Development

Building on prior identified geometric risk factors for mBTTS thrombosis¹⁶ in a local cohort (weight, 2.8–4.1 kg), we developed idealized computational models to systematically assess how specific geometric features influence thrombogenic flow. This study was approved by Washington University’s institutional review board (ID# 202303162; March 29, 2023).

Using SimVascular (SimVascular), we first generated 3-dimensional mBTTS anatomies from patient computed tomography/magnetic resonance imaging. Insights from these models guided creation of a simplified H model (Figure 1, B) that enabled precise control of key geometric parameters while excluding confounding factors. This approach aligns with established computational hemodynamics methods, where idealized models are used to isolate relationships between flow and geometry before reintroducing complexity.^32–37 Inlet and outlet extensions (10 × subclavian vessel diameter) ensured fully developed flow at the mBTTS region (Figure 1, D). Future validation studies will reintroduce patient-specific anatomy to confirm clinical relevance in more complex geometries.

Geometric Parameter Selection and Model Generation

Four key geometric parameters were identified from our previous clinical study.¹⁶ Three parameters were evaluated across the minimum, median, and maximum values observed in our clinical cohort (Table 1) that can be studied in the idealized H model: the distance between shunt insertion and right carotid artery (RCA) insertion $\left({\text{L}}_1\right)$ (−7 to +7 mm, with negative values indicating insertion distal to the RCA), the insertion angle between the shunt and the pulmonary artery (PA) $(\alpha )$ (60°−120°), and the PA diameter ${(\text{D}}_{\text{PA}})$ (2.5–5.0 mm). The fourth parameter, shunt diameter $\left({\text{D}}_{\text{shunt}}\right)$, was evaluated at 2 clinically standard sizes: 3.5 mm and 4.0 mm. This resulted in 54 unique geometric configurations, each maintaining consistent vessel paths while varying only the specified parameters. In addition, shunt length (10.9–20.4 mm) and shunt tortuosity (1.0–1.12 mm) were evaluated specifically in the safest and most thrombogenic configurations because they are controllable by the surgeon. Tortuous grafts curved 1 mm from the proximal anastomosis and 2 mm from the distal anastomosis to a controlled arc length (Figure E1, A). Tortuosity was calculated as the ratio of the trimmed centerline length (arc length) to the straight-line (Euclidean) distance between its start and end points. Different curvature directions (convex vs concave) were analyzed as well (Figure E1).

Inlet flow was constant (1 L/min), and outflow ratios were calculated using Murray’s Law³⁸ based on derived median vessel radii (Table E1), distributing flow amongst the RCA, right subclavian artery, and PAs. This approach simplifies physiology, neglecting the ratio between pulmonary vascular resistance and systemic vascular resistance that would alter flow distributions in vivo. Future models incorporating patient-specific resistances could improve presurgical planning in patients with pulmonary hypertension or altered vascular resistance.

Mesh Generation and Validation

A mesh was generated in geometry using Ansys Fluent (Ansys). To capture boundary layer effects, 10 prismatic inflation layers were added near the vessel walls, with tetrahedral elements throughout the remaining volume (Figure 1, C). Mesh convergence was reached with ~800,000 cells, and mesh quality was evaluated using standard quality metrics with all values within acceptable thresholds (maximum skewness <0.82, minimum orthogonal quality >0.20, and maximum aspect ratio <40).

Computational Flow Analysis

Blood flow simulations were performed using Ansys Fluent with the steady-state, k-$\omega$ shear stress transport turbulence model, with a low Reynolds number correction. Blood was treated as a Newtonian fluid with constant density $\left(\rho =1060\text{kg}/{\text{m}}^3\right)$ and viscosity (0.004 kg/m/sec). We applied no-slip boundary conditions at rigid walls and specified constant velocity inlet profiles based on median flow velocities measured by Doppler echocardiography in our patient cohort.¹⁶ Steady-state flow was used to minimize confounding variables and isolate geometric effects. The inlet Reynolds number (1830), although indicated laminar flow for a straight vessel, could lead to flow instability in bifurcations. For each configuration, we calculated 4 hemodynamic parameters associated with thrombosis risk. WSR was calculated as:

where $v_t$ is blood velocity parallel to the wall and $n$ is the distance perpendicular to the wall. ESR was quantified as the rate of stretching along the flow streamlines:

where $\tilde{E}=\frac{1}{2}\left(\overrightarrow{\nabla }\overrightarrow{v}+(\overrightarrow{\nabla }\overrightarrow{v}{)}^{\text{T}}\right)$ is the strain rate tensor and $\overrightarrow{v}$ is the blood velocity at a point. TI was quantified as:

where $⟨\overrightarrow{v}⟩$ is the time average of $\overrightarrow{v}$ and $k=\frac{1}{2}(\overrightarrow{v}-⟨\overrightarrow{v}⟩)\cdot (\overrightarrow{v}-⟨\overrightarrow{v}⟩)$ is the turbulent kinetic energy representing the mean kinetic energy per unit mass associated with turbulent velocity fluctuations.

Power loss $(PL)$ was quantified dissipation across the model³⁹:

where $p$ is the static pressure, $v_n$ is the velocity component normal to the inlet/outlet surface (positive outward), and $A$ is cross-sectional area.

A parameter $\lambda$ classified extensional versus rotational flow⁴⁰:

where $\tilde{\Omega }=\frac{1}{2}\left(\overrightarrow{\nabla }\overrightarrow{v}-(\overrightarrow{\nabla }\overrightarrow{v}{)}^{\text{T}}\right)$ is the vorticity tensor and $\Vert \tilde{E}\Vert =\sqrt{\mathit{\text{tr}}\left({\tilde{E}}^T\tilde{E}\right)}$ and $\Vert \tilde{\Omega }\Vert \equiv \sqrt{\mathit{\text{tr}}\left({\tilde{\Omega }}^T\tilde{\Omega }\right)}$ are Frobenius norms, in which $tr$ represents the trace of a tensor. $\lambda =1$ corresponds to pure extensional flow, $\lambda =-1$ represents pure rotational flow, and $\lambda =0$ indicates simple shear.

Optimization

To objectively identify optimal geometric configurations, we applied a composite score balancing ESR, WSR, and TI:

where subscript max and min represent extrema over all geometric parameters, and $i$ indicates a specific parameter combination. Weighting factors $\alpha ,\beta$, and $\gamma$ were determined using principal component analysis. Lower $f_i$ represents more optimal configurations.

Statistics

Data were summarized using frequencies with percentages, means ± SD, or medians with interquartile range, as appropriate. Spearman’s correlation coefficient $(r)$ was used to assess associations between parameters and shunt insertion outcomes. Differences among parameters with 3 levels were compared using analysis of variance, followed by Tukey-adjusted pairwise comparisons when the overall model was significant. If normality was violated, rank-transformed data were analyzed. For 2-level parameters, either a 2-sample $t$ test or Wilcoxon rank-sum test was applied based on data distribution. Analyses were conducted in SAS version 9.4 (SAS Institute Inc).

Spatial Distribution of Peak Flow Characteristics

Analysis of the 54 mBTTS configurations revealed distinct peak locations for each flow metric (Figure 3, C, F, and I, and Table E2). WSR peaks occurred primarily at bifurcation points, with 46.3% at the distal attachment corner of the shunt-subclavian artery junction, 29.6% at the proximal attachment corner, 20.4% (11 models) at the RCA insertion, and 3.7% (2 models) along the shunt wall. ESR peaks were similarly concentrated at junctions, with 66.7% occurring at the distal corner of the shunt-subclavian artery junction, 27.8% at the RCA bifurcation, and 1.8% at the shunt-PA junction. TI peaks differed, with 70.4% occurring within the shunt conduit, 11.1% at the inlet, 7.4% distal to the shunt insertion, and 5.6% split equally between the carotid and shunt-PA junctions. Representative contour plots show the distribution of WSR, ESR, and TI across the model (Figure 3, A, B, D, E, G, and H).

Influence of Geometric Parameters on Flow Characteristics

Statistical analysis revealed correlations between geometric parameters and flow characteristics (Figure E2 and Table E3). Moving the shunt insertion from proximal $\left({\text{L}}_1=-7\text{mm}\right)$ to distal $\left({\text{L}}_1=+7\text{mm}\right)$ moderately decreased ESR and TI, with WSR showing a similar but nonsignificant trend $(r=-0.26;P=.054)$. A larger $\alpha$ reduced ESR $(r=-0.40)$ but increased TI $(r=0.37)$. ${\text{D}}_{\text{PA}}$ strongly correlated with all flow metrics $(r>0.3)$, whereas ${\text{D}}_{\text{shunt}}$ was inversely associated with WSR $(r=-0.33)$ and showed weak, nonsignificant effects on ESR and TI.

Overall and pairwise comparisons confirmed these trends (Figure 4). Significant effects were observed for $\alpha$ and ${\text{D}}_{\text{PA}}$ $(P<.05)$, with ${\text{D}}_{\text{PA}}$ increasing all flow metrics, especially TI $(P<.0001)$. ${\text{L}}_1$ significantly influenced ESR $(P=.028)$ and TI $(P=.0042)$, whereas ${\text{D}}_{\text{shunt}}$ selectively influenced WSR $(P=.016)$. These findings identify ${\text{D}}_{\text{PA}}$ and $\alpha$ as key determinants of all flow metrics, with ${\text{L}}_1$ and ${\text{D}}_{\text{shunt}}$ exerting more targeted effects.

Models exhibiting extreme values for each flow metric demonstrated distinct geometric configurations (Table 2). Maximum WSR (61,400 s⁻¹) occurred in a model with ${\text{D}}_{\text{PA}}=5.0\text{mm}$ and ${\text{D}}_{\text{shunt}}=3.5\text{mm}$, whereas minimum WSR (24,200 s⁻¹) was observed with ${\text{D}}_{\text{PA}}=2.5\text{mm}$ and ${\text{D}}_{\text{shunt}}=4.0\text{mm}$. Peak ESR (7910 s⁻¹) was found in a configuration with ${\text{L}}_1=-7\text{mm}$ and $\alpha =60°$, whereas minimum ESR (2740 s⁻¹) occurred with ${\text{L}}_1=7\text{mm}$ and $\alpha =120°$. Maximum TI (33.9%) was observed in the same configuration as peak WSR, whereas minimum TI (6.8%) occurred in a model with ${\text{L}}_1=7\text{mm},\alpha =60°,{\text{D}}_{\text{PA}}=2.5\text{mm}$, and ${\text{D}}_{\text{shunt}}=4.0\text{mm}$.

Geometric Configurations Associated With Flow Optimization

Analysis of all configurations that consistently minimized or maximized potentially adverse flow characteristics was identified through the optimization method defined by equation 6 (Figure 5). By applying the optimization approach, the weighting factors for WSR, ESR, and TI were determined as 0.38, 0.25, and 0.37, respectively, when calculating the optimization factor, $f$ (Table E2). We ranked all 54 models by $f$ value and identified clear trends. The 6 models with the lowest $f$ values shared ${\text{L}}_1=7\text{mm}$ and ${\text{D}}_{\text{PA}}=2.5\text{mm}$, whereas the 5 highest $f$ value models had ${\text{D}}_{\text{shunt}}=3.5\text{mm}$ and ${\text{D}}_{\text{PA}}=5.0\text{mm}$. These patterns suggest $f$ is particularly sensitive to 1 or more of these parameters. Surface plots of $f$ against geometric parameter pairs (Figure E3) further supports this. Lower $f$ values are observed in models with either ${\text{L}}_1=7\text{mm}$ and the smallest ${\text{D}}_{\text{PA}}$, or ${\text{D}}_{\text{shunt}}=4.0\text{mm}$ and the smallest ${\text{D}}_{\text{PA}}$. Together, these findings indicate that a distal shunt insertion $\left({\text{L}}_1=7\text{mm}\right)$, a 4-mm ${\text{D}}_{\text{shunt}}$, and a smaller ${\text{D}}_{\text{PA}}$ are associated with more favorable hemodynamic characteristics. This choice additionally depends on flow resistance, a surrogate for systemic steal (Table E2), which is largest for small shunt diameters. Increasing shunt angle generally increases resistance. Power loss was highest for smaller $\alpha$ and shunt insertion proximal to the carotid bifurcation. Shunt placement influences both in-graft hemodynamics and the amount of flow expected to pass through the shunt.

The configuration demonstrating most favorable flow ($f=0.14,$ WSR = 32,300 s⁻¹, ESR = 3080 s⁻¹, and TI = 6.8%) combines a distal shunt insertion $\left({\text{L}}_1=7\text{mm}\right)$, acute PA angle $\left(\alpha =60°\right)$, minimal ${\text{D}}_{\text{PA}}(2.5\text{mm})$, and larger ${\text{D}}_{\text{shunt}}(4.0\text{mm})$. $f$ remained uniformly low across the domain (Figure 5, A and B). The configuration associated with the highest-risk flow characteristics ($f=0.90,$ WSR = 60,500 s⁻¹, ESR = 5920 s⁻¹, and TI = 33.3%) combined proximal shunt insertion $\left({\text{L}}_1=-7\text{mm}\right)$, obtuse PA angle $(\alpha =120°)$, maximal ${\text{D}}_{\text{PA}}(5.0\text{mm})$, and smaller ${\text{D}}_{\text{shunt}}(3.5\text{mm})$. High $f$ values were concentrated at the shunt-subclavian junction and shunt-PA connection, especially along the shunt and near bifurcating zones, primarily driven by TI. Intermediate configurations showed consistent trends where distal shunt insertions $\left({\text{L}}_1>0\right)$ combined with acute angles $(\alpha =60°)$ reduced $f$ independent of vessel diameters. Models with ${\text{D}}_{\text{PA}}\le 3.75\text{mm}$ demonstrated more moderate flow characteristics across all measured parameters compared with those with larger PA diameters.

To further evaluate thrombogenicity, $\lambda$, representing elongational flow and linked to VWF extension, was mapped across the symmetry plane. For the lowest $f$ configuration, $\lambda$ was centered around 0, corresponding to simple shear, with small positive values (extensional flow) near bifurcations and stagnation points (Figure 5, C). Negative $\lambda$ was minimal at the shunt-subclavian junction, indicating limited rotational flow. For the greatest $f$ configuration, $\lambda >0$ occurred predominantly along the shunt, at the shunt-subclavian and shunt-PA junctions, aligning with elevated ESR. Negative $\lambda$ remained limited, suggesting extensional flow dominated over rotational flow.

Building on the representative least and most thrombogenic configurations, shunt length and tortuosity were varied across the clinically observed range to evaluate their influence on local flow behavior (Figure E4). Across both configurations, shunt length influenced turbulence, but WSR and ESR remained largely unaffected (Figure E5). In the least thrombogenic configuration, the shortest shunt exhibited elevated turbulence (TI = 15.5% at the shunt-PA junction), which decreased and remained low at longer lengths; turbulence remained high across all shunt lengths in the most thrombogenic configuration.

With shunt tortuosity, both configurations showed direction-dependent effects (Figure E6). In the safest configuration, WSR remained nearly unchanged, whereas ESR rose slightly with increasing shunt tortuosity in a way that depended on the direction of curvature (Figure E1); TI remained low and unaffected. In the most thrombogenic configuration, curvature generally reduced WSR relative to the straight shunt, with WSR increasing with shunt tortuosity. Here, ESR trends were nonmonotonic and strongly dependent on curvature direction, whereas TI remained high regardless of tortuosity.

Mechanistic Implications for Thrombosis

Mechanobiological pathways can trigger thrombosis. High WSR can induce shear-induced platelet activation and aggregation through VWF, known to occur as a function of exposure time and shear stresses, observed in <1 second for shear stress exceeding 400 dynes/cm² (10,000 s⁻¹),^17,18 as observed at bifurcations, especially the shunt-subclavian anastomosis. ESR is theoretically shown to induce VWF extension at lower thresholds (100–500 s⁻¹).^19,20 ESR is highest within the core of the shunt flow and at the stagnation region, seen as flow impinges on the PAs, which could potentially prime platelets as they hit the wall. Although not fully understood, WSR and ESR could induce mechanical tension across platelet glycoprotein Ib-IX-V (>20 pN) when binding to VWF, triggering granule release and calcium signaling, thereby activating integrins and promoting platelet aggregation.^E2 Also, combined ESR and WSR may lead to VWF self-association, potentiating further thrombotic risk.²⁸ More studies are needed to fully understand the role of ESR relative to WSR in causing a thrombotic environment.

TI, which peaked within the shunt conduit, represents a third mechanical pathway to thrombosis. Turbulent flow enhances platelet-platelet collisions and exposure to repeated mechanical stress, potentially amplifying the prothrombotic effects initiated at bifurcations.²³ Multiple studies have suggested that transient changes in stress can increase the likelihood for platelet activation.^26,E3 Turbulence can also lead to a converse situation of bleeding risk due to its potential role in acquired von Willebrand syndrome (loss of high molecular weight VWF multimers.)^E4

Shunt length and tortuosity can modulate WSR, ESR, and TI. Short shunts limit flow development, generating local recirculation and elevated turbulence. Longer shunts allowed laminar flow to develop fully. In most favorable configurations, tortuosity caused minimal change in WSR, but increased ESR due to enhanced elongational stretching along converging streamlines. In most thrombogenic configurations, tortuosity reduced WSR and ESR overall. TI remained largely unchanged by tortuosity. These findings indicate that shunt length shapes the development of flow and alters how blood is distributed into the shunt, whereas mild tortuosity has minimal influence on turbulence but shifts shear and elongational stresses in ways that depend on both the shunt’s geometry and the direction of curvature.

Current antithrombosis strategies target biochemical pathways but have limited efficacy on the mechanical pathways driving thrombosis.^{22,32,E5–E8} Traditional antiplatelet agents like aspirin have minimal effect on shear-induced platelet activation and aggregation, whereas anticoagulants target the coagulation cascade, which depends on low convection/flow that prevents washout of reaction products,^E9 potentially explaining their limited efficacy in preventing mBTTS thrombosis. Caplacizumab, the only drug targeting mechanical pathways, indirectly inhibits VWF-glycoprotein Ib-IX-V binding by stabilizing its autoinhibitory module, but it remains expensive and is primarily used for thrombotic thrombocytopenia purpura rather than surgical thrombosis prevention.^E10 Our findings suggest that geometric optimization, by minimizing adverse flow characteristics, could complement pharmacologic strategies, potentially reducing thrombosis risk without increasing bleeding.

Clinical Translation and Surgical Implications

Our computational findings suggest specific surgical adjustments to improve mBTTS outcomes. These include placing the shunt distal to the RCA bifurcation and using a 60° insertion angle into the PA. Although patient anatomical variations (eg, vessel size and branching) may limit full implementation, these hemodynamically informed engineering principles should guide shunt size and placement decisions. In some cases, these adjustments can be made without added complexity, optimizing outcomes while maintaining the core surgical approach.

The 8% to 12% rate of mBTTS thrombosis, with 50% mortality upon occlusion, underscores the need for better strategies.^1–7 Our findings suggest that smaller PAs (≤3.75 mm) paired with larger ${\text{D}}_{\text{shunt}}(4.0\text{mm})$ improve otherwise prothrombotic hemodynamics, challenging traditional views on vessel size matching as the primary guiding principle. However, these recommendations must be weighed against individualized anatomical constraints and surgical risks. In patients with single-ventricle physiology, shunts with shorter length or larger diameter have lower resistance, which can predispose to pulmonary overcirculation and systemic steal. Therefore, whereas larger shunts may reduce thrombosis risk, their potential to drive excessive pulmonary flow that diverts blood from coronary and systemic circulation must be considered in surgical planning.

Because outlet flows were prescribed as fixed fractions of total flow to each vessel, based on their diameters rather than on distal vascular resistance, absolute pulmonary-to-systemic flow ratios are not physiologically predictive. Moreover, our simulations encompassed only the innominate artery, the shunt, and the PAs, so systemic flow was not directly evaluated. To compensate, we report resistance across the shunt as a comparative metric for all configurations (Table E3), reflecting the pressure drop required to drive blood through the graft: lower resistance favors greater pulmonary flow, whereas higher resistance limits it and reduces systemic steal.

Power loss values calculated for mBTTS are higher than those reported in Fontan studies.^39,E11 Greater power loss indicates less-efficient flow, which is expected due to smaller vessel and shunt sizes when compared with Fontan circulation. Nevertheless, our results help identify designs likely to minimize energy loss and hemodynamic risk in future patient-specific modeling.

In some patients, these findings could be implemented with minimal changes to current surgical practice. Preoperative planning should consider shunt placement relative to the RCA bifurcation, and intraoperative consideration should be placed on maintaining acute insertion angles. Shunt-to-vessel diameter ratio also warrants attention due to its influence on flow. Although exact values may vary, general trends, like favoring acute over obtuse angles, are more important. Further studies are needed to define the geometric parameter space fully.

Although these theoretical recommendations may not apply to every patient, they could help identify those at higher risk for thrombosis, guiding clinical vigilance and management. Geometric optimization complements pharmacologic therapy by targeting mechanical contributors to clot formation without requiring new drugs or protocols. These findings also underscore the value of patient-specific preoperative modeling, enabling surgeons to tailor shunt geometry to individual anatomy.

Limitations and Future Directions

This study has several limitations. Geometric features were derived from a small retrospective cohort, limiting generalizability. Despite this, 54 configurations based on these patients provided statistical strength and give insight into trends. Computationally, steady-state inlet conditions replaced pulsatile flow, and blood was modeled as Newtonian, potentially underestimating low-shear-rate effects. Despite mesh independence, resolution limits may miss fine-scale gradients at sharp transitions. The k-$\omega$ turbulence model with low Reynolds corrections may not fully capture transitional flow. Assumptions of rigid vessel walls and simplified boundary conditions excluded compliance and distal vascular bed resistances. Lastly, only a subset of geometric variables was analyzed.

Validation with in vitro flow models is needed to confirm the predictions. Patient-specific simulations that incorporate complete anatomy and pulsatility could refine geometric recommendations, and large-animal studies may help assess clinical impacts. Development of real-time surgical guidance tools could support intraoperative application of these findings. Ultimately, combining geometric optimization with novel antithrombosis strategies may offer a comprehensive approach to reducing mBTTS thrombosis risk.