1.
Introduction
Cyanotic congenital heart disease (CCHD), or ‘blue baby syndrome’, occurs due to a defect in the heart that prevents separation of oxygenated from deoxygenated blood and/or limits flow to the lungs. As a result, the heart pumps blood that is not fully oxygenated and cannot sustain the demands of the body. Systemic to pulmonary shunt procedures are the mainstay in treatment of congenital heart diseases with limited pulmonary blood flow. The shunt, which in essence is a connecting tube, diverts blood that goes to the body into the lungs, enhancing blood oxygenation. The modified Blalock-Taussig shunt (MBTS; also referred to as modified Blalock-Taussig-Thomas shunt or MBTTS) is the most commonly used shunt for this purpose. Depending on the defect pathology, the shunt constitutes either a sole pulmonary blood provider or an auxiliary pathway that increases pulmonary blood flow to ameliorate oxygenation.
Despite numerous improvements, MBTS overall mortality and composite morbidity rates remain relatively high at 7.2% and 13.1%, respectively [1]. The major causes of morbidity and mortality after MBTS insertion are: 1) MBTS occlusion, which frequently occurs due to coagulation inside the shunt as a result of slow flow conditions; and 2) over-shunting, which basically results in pulmonary over-circulation [2],[3]. The key to avoid occlusion and over-shunting complications is to maintain a properly balanced flow of blood: to the lungs, reaching adequate pulmonary perfusion; and through the shunt, minimizing thrombosis risks. Determining the shunt characteristics that are better suited to achieve this goal could be done through virtual surgery approaches: computer simulations that, given shunt characteristics, baby cardiovascular geometries, and physiological conditions, precisely predict outcomes of interventions [4]. In the context of a balanced circulation, the outcome sought is the blood flow dynamics resulting from placement of the shunt. Current virtual surgery approaches, however, present two main limitations: i) prediction accuracy; and ii) time required (and computational resources needed) to perform the simulation [5]. Our study attempts to overcome these limitations.
Simulations of blood flow in the MBTS have progressively become more refined. Flow simulations are now performed using dynamic multiscale models that solve for flow velocities and pressures over several cardiac cycles, using conditions that approximate the patient response from the systemic circulation [6]–[11]. One of the limitations of this approach is that it requires extensive computational time (on the order of days in a super computing cluster), limiting the number of virtual surgery cases that can be simulated while the patient is waiting for interventions. To overcome this difficulty, in this study we explore the accuracy and efficiency of simply computing average per cycle values of flow variables. As expected, this approach drastically decreases computational time (to a couple of hours at most), while accurately estimating clinically relevant mean flow parameters including cardiac index (CI), the ratio of pulmonary to systemic blood flow (Qp/Qs), and systemic oxygen delivery (OD).
To assess the accuracy of simulations, model predictions were compared to clinical data, both patient specific data and previously published data from a cohort of MBTS patients [12]; as well as a previously developed dynamic multiscale computational model [7]. Moreover, we performed simulations in which a patent ductus arteriosus (PDA) was either left intact or resected after MBTS insertion, to determine the impact of a PDA on blood flow dynamics. The PDA is a fetal vessel that allows blood to bypass the non-functional lungs before birth. This conduit, however, closes due to increased blood oxygenation right after birth, but can be kept open in babies with cyanotic heart defects to improve oxygenation to the lungs (due to pressure differences, after birth blood flow through the PTA in cyanotic babies reverses direction and directs blood from the systemic circulation to the lungs). Herein we describe our proposed approach in details highlighting its advantages, disadvantages, and possible future applications.
2.
Materials and methods
2.1. Patient data
We retrospectively obtained patient data from a neonate with a complex CCHD who underwent MBTS insertion as part of palliation. The study was approved by the Oregon Health & Science University Institutional Review Board (IRB00011857). Patient data (see Table 1), which included a computed tomography angiography (CTA) image, were obtained from postoperative clinical procedures. Systemic vascular resistance (SVR), pulmonary vascular resistance (PVR), mean arterial pressure (MAP), hemoglobin concentration (Hb) and body oxygen consumption (CVO2) from the patient were used as input variables for our model; while other parameters (CI, Qp/Qs, OD, oxygen saturations) were used for model validation and to evaluate model predictions.
2.2. Multiscale model of the MBTS circulation
A multiscale computational model was employed to simulate blood flow through the infant circulation. It consisted of two coupled models: a 3D computational fluid dynamics (CFD) model of blood flow in and around the MBTS; and a lumped parameter network (LPN) model that accounted for the effect of the infant circulation on MBTS blood flow. The integrated model, based on [6], was implemented using a custom MATLAB R2014a (The MathWorks, Inc., Natick, MA) code that interfaced with the CFD solver ADINA-F v9.0.6 (ADINA R&D, Inc., Watertown, MA). The geometry of the 3D CFD model (Figure 1) was derived from the patient CTA image; while the LPN model incorporated patient-specific systemic and pulmonary vascular resistances (SVR and PVR, respectively; Figure 2). The heart was modeled as a constant pressure generator, with aortic blood pressure equal to the patient specific MAP (67 mmHg).
Coupling between the 3D CFD and LPN models was achieved through an iterative algorithm. At each CFD-LPN interface, blood pressure, calculated by the LPN, was imposed to the 3D CFD model as a uniform traction boundary condition. Conversely, flow, the integration of normal velocity over the 3D outlet surfaces, was imposed as an inlet condition to the corresponding LPN. To ensure convergence to equilibrium, iterations were performed until the change in pressure and flow at all 3D-LPN interfaces was less than 0.1% between successive iterative steps.
2.2.1. 3D CFD model: region of interest
The 3D CFD model comprised a portion of the aortic arch, innominate artery, pulmonary arteries, MBTS, and outlets to the carotid, subclavian and vertebral arteries (Figure 1). Surface meshes were imported into ADINA and the model geometry volume was meshed with flow-condition-based interpolation (FCBI) elements, which are stable for a wide range of flows. Simulations assumed that blood flow was incompressible and Newtonian, with a density (ρ) of 1060 kg/m3 and a dynamic viscosity (µ) of 0.004 Pa·s [6]. Thus, neglecting the effects of gravity, the equilibrium equations describing the flow of blood in the 3D CFD model are:
where
v is the flow velocity vector, ∇ is the gradient operator, and ∇
2 the Laplacian operator.
Two different geometries were generated (see Figure 1): 1) MBTS with no patent ductus arteriosus (PDA); and 2) MBTS with PDA, in which an artificial conduit with the patient specific PDA preoperative diameter (4.0 mm) was inserted during mesh generation. These geometries allowed us to perform simulations to determine whether leaving a PDA could be an advantage.
2.2.2. LPN model: systemic and pulmonary systems
A simplified LPN (Figure 2; Table 2), was used to model the connection of the 3D CFD model with the rest of the circulation [6]. Vascular resistance (R) was used to simulate the effect of the systemic and pulmonary circulations. The LPN consisted of two subsystems: 1) the systemic, and 2) the pulmonary circulations. The systemic circulation was divided into the upper body and lower body, each of which was assumed to contribute half the total systemic vascular resistance (SVR). Resistance values were set such that the equivalent resistance across the entire systemic circulation was equal to the patient specific SVR (4.32 mmHg·sec·mL−1; see Table 2), and the equivalent resistance across the entire pulmonary circulation was the patient-specific pulmonary vascular resistance (PVR = 0.816 mmHg·sec·mL−1; Table 2). To couple the 3D model to the LPN model, the following boundary condition equations (at vascular outputs) were employed:
Where
Qi is the volume flow rate of blood exiting the 3D model at output
i,
n is a unit vector normal to the 3D output surface,
Ai is the area of the output surface,
Pi is the pressure imposed by the LPN model onto the 3D model at output
i, and
Ri is the resistance associated with the output
i.
In quantifications, PVR was varied from 20% to 200% of the patient PVR (range of PVRindex: 0.68–6.8 mmHg·m2·L−1·min) to simulate the effects of changing PVR in the circulation.
2.3. Oxygen transport model
An oxygen transport model, based on [12], was implemented to predict oxygen delivery (OD), systemic arterial saturation (SatA), and systemic venous saturation (SatSV). In the shunted configuration (see Figure 3), oxygen mass conservation establishes:
where,
CO is cardiac output (
CO =
QS +
QP),
CA is systemic arterial oxygen concentration, C
SV is systemic venous oxygen concentration,
CPV is pulmonary venous oxygen concentration, and
CVO2 is the patient specific whole body oxygen consumption (134.4 mL·min
−1·m
−2; see
Table 1).
Pulmonary venous oxygen saturation (SatPV = 95%; Table 1) was based on patient specific data and pulmonary venous oxygen concentration computed as follows:
where,
O2Cap is the maximal oxygen carrying capacity (in mL/mL blood,
O2Cap = 0.0134 ·
Hb, where
Hb is the patient specific hemoglobin concentration in g/dL). OD was then defined as:
2.4. Blood flow dynamics and thrombosis
Risk of thrombosis in the MBTS and PDA were estimated based on blood flow dynamics. Like in previous studies, e.g. [7],[12], we used the following surrogates associated with thrombotic risk: 1) wall shear stress (WSS), with low values leading to clot formation, and high values to platelet activation leading to coagulation; and 2) residence time (RT), with higher values indicating greater risk of thrombosis due to flow stagnation. We also computed the RT in the pulmonary arteries (PA) as an indicator of over-shunting.
2.4.1. Wall shear stress
WSS were calculated directly from velocity data via post processing in EnSight v10.1.1 (b) (CEI, Inc., Apex, NC) using a custom Python script with the following equations,
where,
σ(
x) is the stress tensor evaluated at point
x, µ is dynamic viscosity,
v is the velocity vector,
n is the wall surface unit normal vector, ∇ is the gradient operator, and T indicates transposition.
2.4.2. Residence time
RT calculations were performed in EnSight using a particle based method. To measure RT in the MBTS, 1024 evenly spaced particles were ‘emitted’ from a plane intersecting the shunt just below the proximal anastomosis (where flow enters the MBTS). The time for each particle to exit the MBTS, defined as the time it took the particle to reach a plane intersecting just above the distal anastomosis (exiting plane), was recorded. These results were then averaged to calculate mean RT. Similarly, to measure RT in the pulmonary arteries, 1024 particles were emitted from mid-planes intersecting each shunt (i.e. MBTS and PDA) right before flow entered the pulmonary arteries. Particles were then tracked to the boundary of the 3D domain with the time in the 3D modeled portion of the pulmonary arteries recorded.
3.
Results
3.1. Mesh sensitivity
To determine the best mesh for simulations, we first completed a mesh convergence analysis. To this end, we performed CFD simulations of the patient 3D geometry, which we progressively meshed with smaller elements. We computed cardiac index (CI) for each simulation, and plotted it versus number of elements in the mesh to determine convergence. CI convergence was defined as the point at which increasing the number of elements in the mesh changed the resulting CI by <1 %. Mesh convergence optimizes computational cost versus accuracy. Convergence was found when average element length was 0.25 mm (∼1.5 million elements), and this mesh was then used for all simulations.
3.2. Model results and validation
We simulated the patient-specific model using the selected mesh (∼1.5 million elements). An important consideration is that in our model (and the patient simulated) there was no connection between the heart and the pulmonary circulation (the pulmonary artery was blocked in the patient, pulmonary atresia), so that the pulmonary circulation relied on the shunt diverting blood from the systemic circulation. Two cases were considered: 1) the MBTS was the only shunt communicating the systemic and pulmonary circulations; and 2) in addition to the MBTS a PDA was shunting blood to the pulmonary circulation. In addition, we performed simulations in which PVR was changed from 20% to 200% of the patient value, to mimic possible manipulations that surgeons could perform to change PVR in the infant, and determined outcomes. Results obtained are summarized in Table 3.
3.2.1. Validation against previous dynamic model
To validate our results (Table 3), we first compared them to cycle-average data from a previously published dynamic model [7]. Comparisons focused on mean flow parameters describing overall circulation performance. We found that flow parameters were extremely consistent between the models, with percent differences for mean CI, Qp, Qs and Qp/Qs within 1.5% (Table 4). These results suggest our model can accurately capture mean flow characteristics.
3.2.2. Validation against patient data
Predictions from our model were then compared to data from the specific patient modeled as well as clinical data from a cohort of 28 Norwood patients [12], represented as the mean and standard deviation of individual patient data. Clinical data featured only patients who had an MBTS with no PDA. Our model predictions, however, were within one standard deviation of the Norwood cohort averaged data (see Figures 4 and 5). When comparing predicted and patient-specific data, percent differences in CI and Qp/Qs were ∼40% and ∼20%, respectively (Table 5), highlighting uncertainties in both patient-specific measurements and modeling parameters.
3.3. Case results: MBTS alone vs MBTS with PDA
Our results were in general agreement with those of other groups considering similar intervention scenarios with idealized models [6],[7],[13]. Leaving the PDA intact resulted in a higher CI (increased by 57%) and Qp/Qs (increased by 117%) than when considering an MBTS alone (Figure 4a; and Table 3). While Qs remained relatively constant, Qp (lung perfusion) increased, leading to higher blood oxygen saturation and increased OD (Figure 4b), at the expense of over-shunting.
Moreover, parameters including CI, Qp/Qs, and pulmonary arterial pressure displayed a stronger dependency on PVR when the PDA was intact (e.g. steeper slopes in Figures 4a and 5); with pulmonary arterial pressure (PAP) exhibiting the strongest relationship with PVR. The dependency of blood oxygenation and OD on PVR, however, was relatively weak (Figure 4b).
The presence of a PDA drastically changed local hemodynamics around the MBTS, particularly in the pulmonary arteries (Figure 6). Keeping the PDA decreased flow through the MBTS despite an increase in Qp. However, the distribution of flow to the right and left pulmonary arteries was slightly more symmetrical when both the MBTS and PDA were present (see Figure 6b).
3.4. Thrombosis risk and over-shunting indicators
Compared to the case with an MBTS alone, including the PDA increased the spatially averaged WSS over the entire 3D domain by approximately 42% (Figure 7a, b). Because WSS is directly related to velocity, PVR had a greater effect with an intact PDA. On the segment of the pulmonary artery between the MBTS and PDA, where flow competition would be likely to occur, average WSS increased by approximately 116% with the addition of a PDA. Conversely, because flow through the MBTS decreased with PDA, spatially averaged WSS on the MBTS decreased by approximately 25%. The highest WSS (590 dyn·cm−2) was observed on the PDA in the case with the lowest PVR. Over the simulated range of PVR, on average, spatially averaged WSS on the PDA was 32% greater compared to the corresponding MBTS.
Because PDA caused a decreased flow through the MBTS, mean RT in the MBTS increased by approximately 45% (Figure 8). However, because pulmonary perfusion increased with an intact PDA, RT in the pulmonary arteries decreased by approximately 48% (Figure 8).
4.
Discussion
As any model, our model of the MBTS circulation has limitations, advantages and disadvantages.
4.1. Model limitations and potential improvements
While physiological flow is pulsatile, our model (by design) does not capture the dynamics of pulsatile flow. Nevertheless, our results compared well with clinical data (e.g. mean flow, CI, Qp/Qs), and extremely well with predictions from dynamic models (averaged over the cardiac cycle; Table 4). Therefore our computations allow for a quick screening and testing of different scenarios before selecting the most promising ones for a full dynamic simulation.
The assumption that the heart is a constant pressure source is another limitation. Based on clinical data from a previous study [14], however, during the early postoperative period following MBTS insertion (< 18 hrs), mean arterial pressure (MAP) is about 15% lower in patients with PDA than patients without PDA; but this difference becomes negligible after 18 hrs. Thus, our model describes this later dynamics. In any case, implementation of a simplified cardiac model, while outside the scope of this paper, would circumvent the assumption that the heart is a constant pressure source in the future.
Uncertainties in modeling patient-specific cases remain mainly in finding good mathematical descriptions of the circulation model and properly assessing patient data uncertainties [5],[15]. With advances in imaging technologies, uncertainties in delineating and segmenting images for modelling purposes are becoming less of an issue. However, it is still imperative to find accurate parameters for the LPN models that closely represent the infant circulation under consideration. In this regard, our model used published parameters that were adjusted proportionally only to ensure that SVR and PVR were the same as the patient under consideration. This approach was not enough (see Table 5). While outside the scope of this paper, a thorough uncertainty and sensitivity analysis is required to accurately represent patients and properly manage data uncertainty [15],[16]. In particular, resistance values employed are key, and more sensitive than dynamic parameters. Uncertainty and sensitivity analyses will be far easy and effective to perform (at least initially) on our simplified model than in a dynamic model of the infant circulation. To improve our modeling technologies, therefore, future work needs to focus on uncertainty and sensitivity analyses [15].
4.2. Validity of the approach
In clinical practice, accurate and meaningful results are needed in the shortest possible time. This is particularly important if time-sensitive intervention decisions were to rely on simulated scenarios. Our approach was designed to compute clinically relevant mean circulation variables (e.g. CI, Qp/Qs) at a fraction of the time involved in performing simulations accounting for the dynamics of several cardiac cycles, which can easily take several days to compute even in a supercomputer. Our proposed model (with simulations that take only a few hours) seems sufficient for screening scenarios, from which a few could then be chosen for a more in depth analysis if needed.
Our model results compared well against both a fully dynamic model [7] (Table 4) and patient cohort data [12] (Figures 4 and 5). The small difference between our results and mean data from state-of-the-art dynamic models (< 1.5%) suggests that, in computing mean circulation data, our model performs as well as more involved models. Agreement of our results with cohort data further provides confidence in our model predictive capabilities, although there is room for improvement both in our approach as well as that of others. As mentioned above, we believe that this improvement will require a through sensitivity and uncertainty analysis of models.
Differences between patient data and model predictions (Table 4) highlight both measurement errors associated with clinical data and model uncertainties, and were comparable to differences reported when using more involved dynamic models [12]. Future studies addressing these uncertainties are needed to determine the confidence with which we can model a specific patient under consideration, and how this confidence can be improved to increase the accuracy of patient-specific computational predictions [15]. More research is needed in this direction if we are hoping to rely on simulated data for intervention decisions. Our results are beginning to highlight similarities and differences between dynamic and static flow models for evaluating CCHD scenarios and their associated advantages and disadvantages. Simulation of more cases (both dynamic and static) are nevertheless needed to fully develop and validate an effective patient-specific approach.
4.3. Insights gained
4.3.1. PVR control
The advantages of noninvasively managing circulation characteristics are inherently obvious. For example, inhaled nitric oxide (iNO) selectively decreases PVR in neonates, including those with pulmonary hypertension [17], providing PVR control. Indeed, model predictions confirm the effectiveness of managing pulmonary hypertension by decreasing PVR in neonates with an MBTS (Figure 5). With an intact PDA and MBTS the influence of PVR in regulating pulmonary pressure became more pronounced (compared to the case with no PDA; Figure 5). Potential benefits gained from increased PVR response, however, appear to be outweighed by the associated increase in cardiac output and pulmonary perfusion (Figures 4 and 5). Other therapies including sub-atmospheric oxygen and induced respiratory acidosis have been shown to affect PVR as well as other parameters such as SVR [18]. Future studies that model the effects of manipulating multiple variables simultaneously could provide insight into which parameters or combination of parameters result in the most efficient circulatory control techniques.
4.3.2. Thrombosis risk
Previous studies have associated platelet activation and thrombosis with regions of very high and very low WSS [19]. In addition to WSS magnitude, exposure or residence time (i.e. RT) has been implicated as an important factor. PDA resulted in markedly increased WSS in the pulmonary arteries while leading to decreased WSS in the MBTS. Although in each case average WSS in the MBTS was below previously reported values for platelet activation (i.e. < 315 dyn/cm2) [20], WSS near anastomoses, was greater, suggesting possible platelet activation and thrombotic risk. The increase in RT within the MBTS for the case with PDA may increase the risk of shunt thrombosis by allowing for a greater opportunity for platelet deposition and agammaegation following activation. Although flow through both the MBTS and PDA was mostly laminar (Re < 2000 for most cases; exceptions occurred for a few of the MBTS with PDA cases with the lowest PVR index values; Re ∼2200), regions immediately preceding the proximal anastomoses and the pulmonary arteries following the distal anastomoses exhibited recirculation, the later appearing more pronounced. Consequently, platelet activation could lead to the expansion of thrombosis from the MBTS into the pulmonary arteries near the anastomosis. Furthermore, the thrombus formed could conceivably lyse (break off) leading to complication such as pulmonary embolism, heart attack, or stroke.
Although model predictions suggest a possible increase in thrombotic risk following PDA preservation, the association is not immediately clear. As previously explained, MBTS with PDA resulted in a non-optimized shunt configuration which lead to elevated cardiac output and pulmonary perfusion. Because the surrogates used for assessing thrombotic risk are all functions of blood flow, predicted WSS within the MBTS were naturally lower whereas RT was expectedly higher. From our results, we can infer that MBTS with PDA does not only lead to over-shunting but may also increase the risk of thrombosis within the MBTS. This suggests that the two leading causes of morbidity and mortality associated with MBTS insertion, shunt thrombosis and over-shunting, each involve a blood coagulation component, the former being more obvious than the latter.
4.3.3. PDA preservation vs ligation
Thrombosis represents the most likely cause of shunt occlusion. Although preserving the PDA allows for shunt redundancy, model predictions suggest that an MBTS with a PDA lead to pulmonary over-circulation and hypertension, while increasing the risk of thrombosis as flow through the MBTS decreases with respect to a configuration without a PDA. The main benefit of redundancy appears to be greater time for re-intervention in case of complications; but this must be weighed against thrombotic risk.
PDA preservation, however, may be beneficial if establishing sufficient pulmonary circulation renders itself difficult. For example, in young neonates (< 2 wk) when PVR is naturally high, the PDA may be initially beneficial to help establish, together with the MBTS, sufficient pulmonary perfusion. In such situations, the PDA could initially act as a temporary “self-closing” shunt.
4.4. Other considerations and possible model extensions
As presented, our computational model applies for both univentricular and biventricular hearts, when there are no additional sources of blood flow to the pulmonary circulation other than the MBTS or MBTS plus PDA. This is because the heart was modeled as a pressure generator (using mean arterial pressure, MAP). Blood flow in our model adjusts to the patient MAP, reflecting the patient condition. Moreover, our model includes the geometry of the innominate, subclavian and pulmonary arteries (see Figure 1), and therefore variations in the size of these arteries, that occur in clinical practice, are taken into account. In clinical situations in which the subclavian artery is very small and flow limiting, or the arterial tree is anomalous, other arteries, such as the brachiocephalic artery are chosen as the inflow (proximal artery) for the shunt [21],[22]. The computational model presented here could be easily extended to simulate the effect of various shunt sizes in situations in which the brachiocephalic artery, rather than the subclavian artery, is used as the shunt inflow.
The computational model described in this study considered a left aortic arch MBTS. In cases of complex cyanotic heart defects, however, aortic positions are variable, and sometimes the right aortic arch is used for the MBTS. Once again, our computational model could be easily extended for simulating a right aortic arch variant of the MBTS. This is indeed the case since the flow dynamics are similar in both cases and the only difference between these variants is the spatial configuration of the arch in relation to the pulmonary artery, which is taken into account by the modelled portion of the arterial branches (Figure 1). While we presented here a very specific example using the arterial geometry of one patient to model blood flow for that particular case, the model employed is versatile and can be adjusted to diverse anomalous geometries and patient situations.
5.
Conclusions
The approach presented here could provide a step forward in intervention planning. Faster time to results (simulation time) can lead to testing of more scenarios, and to undergo sensitivity studies to fine-tune patient response parameters (e.g. LPN parameters) to more accurately model specific patients. The benefits of reduced computational expense in quickly estimating clinically relevant parameters (e.g. CI and Qp/Qs), enables a more thorough screen of surgical options, opening the doors for enhanced intervention planning through virtual surgery approaches. Patient specific CFD models could prove to be a useful clinical tool in the future by allowing for the preoperative simulation of a multitude of surgical scenarios and help inform clinicians as to the selection of the optimal approach.
Model predictions indicate that PDA ligation/preservation during MBTS insertion has a dramatic impact on postoperative circulation, the efficacy of which may be patient dependent. Our model also predicts that manipulating PVR is an effective means for reducing pulmonary arterial pressure in neonates with MBTS dependent circulations, without drastically altering overall hemodynamics. Similar methodologies could be applied to a range of interventions beyond those intended to treat CCHD and could aid in the creation of new surgical strategies for cardiovascular disease.
DownLoad:
