Fluid dynamic model for extracorporeal membrane oxygenation support and perfusion in cardiogenic shock

The cardiogenic shock (CS) is a pathological condition, usually a consequence of myocardial infarction, which leads to a significant decrease in the left ventricle (LV) function in terms of cardiac output (CO).¹ The venoarterial extracorporeal membrane oxygenation (VA-ECMO) is used as a cardiac support technique to compensate the CO reduction caused by cardiac shocks and to provide new oxygenated blood to the arterial circulation. In practice, the VA-ECMO system draws blood from the venous circulation, forces it into an oxygenator, and pumps into the arterial system. In most cases, the oxygenated blood is inserted in the arterial circulation peripherally by imposing a percutaneous cannulation on the iliac artery.² Even if the VA-ECMO is widely used in clinics, there are many risks and pitfalls linked with its adoption. An example is given by the risk of LV afterload increase caused by the retrograde flow of VA-ECMO.³ For this reason, different LV unloading techniques are usually adopted in clinics.⁴ The right balance between perfusion and LV loading resides mainly in the expertise of the clinician. Another point of insight to be investigated is the interaction between VA-ECMO and native flow from the LV. Given their opposite direction of flow, a mixing zone is formed along the aorta denominated as watershed zone.⁵ The formation of the watershed zone is linked with oxygenation unbalancing between the upper and lower parts of the body, also known as Harlequin syndrome.⁶ Given the context of VA-ECMO, it appears clear that new investigations are still required to better understand its fluid dynamics. To achieve this, computational fluid dynamics (CFD) represent a valid tool to explore different hemodynamic issues.^7,8 In the state of the art, various numerical approaches were proposed for the analysis of the VA-ECMO system.⁹ A single group presented a fluid–structure interaction approach to go beyond the rigid wall hypothesis of CFD.¹⁰ Nevertheless, the a-priori knowledge of aortic stiffness remains a challenging task.^11,12 Some previous studies focused on localized fluid dynamics linked with the cannulation technique and its relative position along the aortic lumen.^13–15 Given the context, it appears that different approaches are limited by the adoption of ideal geometries⁸ or by the choice of nonspecific boundary conditions.¹⁶ While the position of watershed zone is deeply investigated, only a few groups focus their attention on the perfusion of the different aortic vessels.^9,17 Moreover, the effect of ECMO support on afterload is usually neglected or partially considered. The unbalancing of vessel perfusion and the systemic pressure increase linked with ECMO support remains a relevant clinical issue.

The principal aim of this study is to present a numerical method for VA-ECMO analysis. In particular, a CFD approach for a complete patient-specific aortic case is developed by considering anatomical data from computed tomography (CT). The model was finalized to investigate the modification of systemic afterload, the position of the watershed zone, and the perfusion of aortic branches, including supra aortic vessels and renal arteries. To cover a wide range of cases, different VA-ECMO support levels and cardiogenic shock conditions were simulated. The numerical results will be then processed in order to find the minimum ECMO support level to grant the perfusion of all the involved arteries and avoid unnecessary afterload increase. To do this, a novel model to correlate cardiogenic shock severity to the minimum ECMO support level was presented. From a technical point of view, the shock severity was expressed as a percentage reduction of the LV cardiac output with respect to the healthy value.

In this section, the main steps of the numerical workflow are summarized. In particular, the patient-specific data processing in terms of geometry reconstruction is presented first. Then, the numerical setup is described with the definition of the boundary condition to describe the combined scenarios in terms of CO severity and ECMO levels. The different simulated cases are presented, and the different output parameters are defined.

First, the patient-specific data from clinical imaging are processed. The CT dataset of a whole aorta was considered for the study [Fig. 1(a)]. The CT dataset was adopted to obtain the morphological data.¹⁸ A semi-automated segmentation algorithm was adopted to obtain the full aortic geometry, including the ascending and descending sections.¹⁹ To allow for a suitable smoothness, a global Taubin spatial filter, with a weighting factor of 0.5 and the maximum number of iterations of 15, was imposed on the raw aortic surface right after segmentation.²⁰ The arteries considered are: left and right coronary arteries (LCA-RCA), brachiocephalic artery (BCA), left common carotid artery (LCCA), left common subclavian artery (LCSA), and left and right iliac artery (LIA-RIA) [Fig. 1(b)]. To simplify the abdominal portion of model, the left and right renal arteries (LRA-RRA) were segmented and used to represent all the vessels of splanchnic and renal circulation. According to the clinical guidelines,²¹ a 20-French cannula was chosen for the simulation, to impose the VA-ECMO flow. The cannula was designed starting from an actual clinical cannula, in a 3D computer-aided design (CAD) environment. The CAD model was then positioned at the right iliac artery level [Fig. 1(c)]. Its geometry was reproduced from product data sheets of commercially available cannulae, and it was positioned in order to have the tip main axis aligned with the RIA outlet centerline. The cannula dimension as well as the cannulation site and positioning were in agreement with the standard clinical practice.

The analysis was carried out with ANSYS Fluent software (Canonsburg, Pennsylvania). The domain was discretized with a water-tight polyhedral mesh with 1.1 × 10⁶ elements by also including five layers of inflation. Polyhedral cells were used as they provide more robust and more accurate computational results, they are less sensitive to geometry tortuosity, and they produce computationally less expensive grids compared to standard tetrahedral meshes.²³ The initial computation domain was uniformly and sequentially refined until grid-independent results were attained.

After introducing the numerical setup, a set of case studies was defined. The simulation was set up to emulate a patient case with a CS condition but without respiratory failure. In particular, the case studies were set according to different imposed inlet conditions Q_NAT [Fig. 2(b)] and Q_ECMO. In particular, the different degrees of Q_NAT are expressed to the cardiac output (CO), while for those of Q_ECMO the flow values are reported. First, the healthy conditions (HLT) were simulated. For the healthy conditions, an aortic waveform was considered for Q_NAT with CO^HLT = 4.90 l/min, with no ECMO support (Q_ECMO = 0 l/min). After simulating the healthy baseline, three different pathological cases of cardiogenic shock (CS) were imposed. The CS conditions in terms of Q_NAT were simulated by scaling the healthy aortic flow waveform to impose a 30% (CS_30%), 50% (CS_50%), and 70% (CS_70%) reduction of CO, corresponding to mild (CO = 3.43 l/min), medium (CO = 2.45 l/min), and severe (CO = 1.47 l/min) CS conditions [Fig. 2(b)]. For each CS case, four levels of ECMO support were considered, with Q_ECMO ranging from 0 to 6 l/min with steps of 2 l/min. A total of 13 scenarios were considered.

Concerning the results, the data regarding the afterload are presented first. In particular, Fig. 3 summarizes the maps of $P ¯ sys$ for CS_30%, CS_50%, and CS_70% cases for the different levels of ECMO support. The $P ¯ sys$ map of HLT case is presented as a reference as well.

For each case, the domain-averaged $P ¯ sys$ values were plotted as a function of Q_ECMO, as depicted in Fig. 4.

The position of the watershed zone is presented in Figs. 5 and 6 for the systolic peak and diastole phase, respectively. For both phases, flow trajectories are presented for CS_30%, CS_50%, and CS_70% conditions and for all the different ECMO support levels, excluding the Q_ECMO = 0 l/min case. For each case, a magnification of the mixing zone position is also reported. The results concerning the perfusion of the aortic arteries are then presented in the graphs of Fig. 7. In particular, the $Q ¯$ values for each outlet at different CS and ECMO support conditions are reported: for each outlet, the mean flow level (⁠ $Q ¯$⁠) for the healthy case is reported with dashed lines as baseline reference. The $Q ¯$ trends are reported for CS_30% [Fig. 7(a)], CS_50% [Fig. 7(b)], and CS_70% [Fig. 7(c)] conditions at the different levels of ECMO support. The trends from Fig. 7 were fitted according to the linear model of Eq. (3). The correlation coefficients R² from the fitting are reported in Table II for each outlet at the different CS conditions. From the trends reported in Fig. 7, the values of $Q ECMO HLT , i$ were also calculated by evaluating the intersection between the trends and the healthy flow reference for each outlet. The $Q ECMO HLT , i$ values for all outlets are reported in Table III for each CS case. From the results reported in Table III, by considering the maximum $Q ECMO HLT , i$ according to Eq. (6), the minimum level of ECMO support to grant the overall flow reestablishment was calculated for each CS condition and is reported in Table IV.

Finally, the values of $Q ECMO HLT / C O HLT$ were evaluated. The values of $Q ECMO HLT / C O HLT$ for each cardiogenic shock percentage were calculated by considering the maximum value of $Q ECMO HLT , i$ for each shock percentage value p and by using CO^HLT as a normalization factor. The trend of $Q ECMO HLT / C O HLT$ as a function of p is reported in Fig. 8, together with the corresponding linear fit. The linearity of the trend according to the model of Eq. (6) was checked as well, and it revealed a $R 2 > 0.99$ with $Ω = 1.41$⁠.

The results presented in Sec. III demonstrate the effects of ECMO support on the different fluid dynamic parameters like afterload, watershed zone position, and arteries perfusion. In particular, different cases of shock and ECMO level support were simulated on a patient-specific geometry of a whole aorta. The obtained outcomes were successfully used as a baseline to develop a preliminary tool for defining the minimum ECMO level support necessary for a satisfactory perfusion.

The results in terms of afterload distribution are presented in Fig. 3. Indeed, as expected, the presence of ECMO support strongly influences the afterload, as the systemic pressure monotonically rises within the whole aorta with the increasing levels of ECMO (Fig. 4). This trend is in accordance with the clinical evidence encountered in the literature.^32,33 It is important to notice that the worst condition in terms of afterload is obtained in CS_30% condition for $Q ECMO =$ 6 l/min, which corresponds to a mild shock with the maximum level of ECMO support. In this worst afterload scenario, the systemic pressure reaches values of 175 mm Hg. This aspect confirms that excessive ECMO support in non-severe conditions of heart shock might lead to unwanted afterload increases. For this reason, the definition of a minimum level of ECMO level to grant the perfusion of arteries becomes fundamental to avoid unnecessary afterload soaring.

Concerning the watershed zone position, the results are presented in Figs. 5 and 6. It appears evident that, regardless of the CS and ECMO support conditions, the position of the mixing zone between native and ECMO flow changes significantly during the cardiac cycle. It is interesting to note that the ECMO flow always reaches the ascending aorta during the diastolic phase (Fig. 6). On the contrary, the position of the watershed zone is strongly dependent on CS and ECMO support conditions during the systolic phase (Fig. 5). In particular, from Fig. 5, it appears that, for the CS_30% shock condition with Q_ECMO levels greater than 2 l/min, the renal region is mostly perfused by ECMO flow. The same trend was encountered also for the CS_50% shock condition. The only configuration providing ECMO flow in systole up to the aortic arch results to be the CS_70% with $Q ECMO =$ 6 l/min. Moreover, as we can notice from Figs. 5 and 6, in all the simulated cases, the coronary arteries are not reached by the ECMO contribution and they are always perfused by native flow only.

Concerning the perfusion results, the different trends of $Q ¯$ in the simulated conditions are reported for each outlet in Fig. 7. The results demonstrate that it was possible to reestablish the healthy flow in each artery with the adequate level of ECMO support, regardless of the CS condition. The results presented in Table II show that the R² values for the linear model of Eq. (3) are always greater than 0.99. This confirms that the perfusion of the different arteries follows a linear behavior with the ECMO support level, regardless of the vessel position and shock conditions. It appears that the linear trend slope of each vessel output is instead depending on the vessel itself.

By inspecting the values of $Q ECMO HLT , i$ from Table III, it is possible to evaluate which are the most difficult arteries to perfuse through ECMO support. It appears that the LIA vessel is the artery requiring the highest level of ECMO support to reach the healthy flow in each CS case. The resulting minimum level of ECMO for flow reestablishment is, in fact, corresponding directly to LIA, as reported in Table IV. This behavior denotes that iliac artery without the cannula is actually the hardest artery to perfuse with the VA-ECMO configuration. The difficulty in perfusion is determined by flow split–curve slope ratio (⁠ $f LIA / m LIA$⁠), as determined in Eq. (6). According to this line of thought, it is possible to assume that LIA requires a significant flow to be perfused (high f^LIA) and, at the same time, it exhibits a poor response to VA-ECMO support given its vessel position and morphology (low m^LIA).

Finally, the $Q ECMO HLT / C O HLT$ trend as a function of shock percentage is reported in Fig. 8. The data were fitted to the model proposed in Eq. (6). The trend was confirmed to be linear through the evaluation of the R² coefficient, which was greater than 0.99. The slope resulting from the fitting was Ω = 1.41. These findings demonstrate the possibility to have a preliminary tool to establish the minimum level of ECMO support for overall perfusion as a function of CS percentage. This concept means that, given a certain percentage of CS, if the $Q ECMO ≥ Q ECMO HLT$ condition is met, the ECMO support grants the overall perfusion. In addition, from the results, it also emerged that the systemic pressure, linked to afterload, monotonically rises with the level of ECMO (see Fig. 4). On the basis of these two phenomena, the previous condition (⁠ $Q ECMO ≥ Q ECMO HLT$⁠) becomes $Q ECMO = Q ECMO HLT$⁠. For this new condition, the perfusion is granted and unnecessary afterload increases are avoided. It is worth underlining also that this model could play a significant role in VA-ECMO weaning process, which, in the current state of the art, still lacks a precise planning.³⁴ The correlation between ECMO minimum level and CS percentage could be the basis for the definition of a successful weaning strategy, in which the support is diminished gradually in correspondence with the heart function.

Numerical simulations of VA-ECMO was investigated by previous authors in the literature. Their main focus was put on flow mixing dynamics on single shock cases,^8,35 idealized geometries⁹ or cannula positioning effects.¹⁴ The main clinical issue of VA-ECMO support, linked with the impossibility to establish the suitable support level to avoid afterload increases and provide a satisfactory perfusion,^34,36 seems to remain open. Key findings of this work are in line with these clinical issues. Indeed, the presented workflow is, to our knowledge, the first one trying to explore the perfusion behavior as a function of ECMO support level and different cardiogenic shock conditions from an overall point of view. Additionally, for the first time, a model to correlate the minimum level of ECMO support to the severity of cardiogenic shock was presented. The minimum level of ECMO support would allow an effective perfusion of the outlet arteries and a minimum afterload increase. The discussed results represent a significant first step in developing a predictive tool to be used in the ECMO support evaluation.

The current work presents some limitations and points of future development. The current model, like already proposed literature contributions,^8,10 does not take into account for the oxygen saturation in the blood supplied by ECMO and native circulation. In this study, a case of a patient with cardiac shock and no respiratory failure was simulated. By assuming this, the native and ECMO flows were qualitatively assumed as blood fluids with the same oxygen saturation and with no actual overall quantification. For future work, it would be interesting to study cases with concomitant cardiac shock and respiratory failure, in which blood saturation plays a key role together with arteries perfusion. In the proposed computational model, no turbulence was considered and the stroke volume during the cardiac cycle is not directly effected by the systemic pressure increase. It would be interesting in the future to account also for wall motion^37–40 and LV distensibility and consequently correlate the stroke volume with the afterload also including a more realistic aortic valve inflow profile.⁴¹ Nevertheless, these considerations would require assumptions on the LV elasticity and the current state of the art heavily relies on rigid wall hypotheses for computational evaluation of VA-ECMO systems. The lumped model for the coronary arteries used in the presented analysis does not account for the left ventricle pressure of the microvasculature during systole. However, the cardiogenic shock conditions provide less ventricular compression. Moreover, the three-element Windkessel model was tuned to maintain the systemic pressure and flow split and, consequently, the average fluid dynamics during the cardiac cycle are still correctly coped. An additional point of development would be to explore how different aorta morphologies, cannula sizes, and orientations influence the watershed zone position and the perfusion factor Ω. This could be done by coupling the simulations to statistic shape model techniques.⁴² Different effects on the perfusion factor Ω are yet to be explored, such as the influence of boundary conditions modifications. Finally, the current work considers a single position of the VA-ECMO cannula. Moreover, the cannula size and orientation effect have been widely studied in the state of the art and their influence on the velocity distribution and the flow within the iliac artery was established.^43,44 Nevertheless, the main focus of the study was to evaluate the perfusion at different ECMO support levels. It would be interesting to also analyze the effect of cannula position, size, and orientation on the perfusion and how this parameter affects the proposed model.⁴⁵ 

In this manuscript, a CFD approach was presented to characterize systemic pressure, watershed zone position, and perfusion of arteries along the aorta. The work demonstrated that it is possible to define a minimum level of ECMO support to grant perfusion and to avoid unnecessary afterload increases. Moreover, a model correlating this minimum level and the cardiogenic shock percentage was proposed and validated with the numerical data. These important achievements represent a first significant step toward the introduction of the concept that ECMO support should match the severity of cardiogenic shock and give a tool to indicate the level of ECMO support for that given shock.

The authors have no conflicts to disclose.

E. Vignali and E. Gasparotti contributed equally to this work.

Emanuele Vignali: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Emanuele Gasparotti: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Writing – original draft (equal). Dorela Haxhiademi: Conceptualization (equal); Validation (lead); Writing – original draft (equal). Simona Celi: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Resources (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.