Assessment of haemolysis models for a positive-displacement total artificial heart

Abstract

The assessment and reduction of haemolysis within mechanical circulatory support (MCS) remains a concern with regard to device safety and regulatory approval. Numerical methods for predicting haemolysis have typically been applied to rotary MCS devices and the extent to which these methods apply to positive-displacement MCS is unclear. The aim of this study was to evaluate the suitability of these methods for assessing haemolysis in positive-displacement blood pumps. Eulerian scalar-transport and Lagrangian particle-tracking approaches derived from the shear-based power-law relationship were used to calculate haemolysis in a computational fluid dynamics model of the Realheart total artificial heart. A range of power-law constants and their effect on simulated haemolysis were also investigated. Both Eulerian and Lagrangian methods identified the same key mechanism of haemolysis: leakage flow through the bileaflet valves. Whilst the magnitude of haemolysis varied with different power-law constants, the method of haemolysis generation remained consistent. The Eulerian method was more robust and reliable at identifying sites of haemolysis generation, as it was able to capture the persistent leakage flow throughout the entire pumping cycle. This study paves the way for different positive-displacement MCS devices to be compared across different operating conditions, enabling the optimisation of these pumps for improved patient outcomes.

Introduction

Heart failure impacts approximately 64 million people globally. ¹ For severe cases, a heart transplant is recommended, however the limited number of donor hearts have led to growing waiting lists. ² Longer wait times increase the risk of transplant failure, ³ so to bridge the gap to transplant, doctors can use mechanical circulatory support (MCS) devices to assist in cardiac function.

Most MCS devices can be separated into two categories based on pumping method: rotary and positive-displacement. Rotary devices use axial or centrifugal rotation to impart energy to the fluid, generating a pressure head, whilst positive-displacement devices repeatedly displace a fixed fluid volume. Ventricular assist devices (VADs) are usually axial ⁴ or centrifugal^5,6 though some are positive-displacement.^7–10 Total artificial hearts (TAHs) are usually positive-displacement pumps,^11–13 however some are rotary. ¹⁴ The Realheart TAH, developed by Scandinavian Real Heart, is a novel four-chamber positive-displacement TAH that has both atria and ventricles, and uses a displacing atrioventricular (AV) plane to produce pulsatile flow governed by a pair of bileaflet mechanical heart valves. ¹⁵

Haemolysis and thromboemboli are noted side effects of chronic MCS. Haemolysis is the damage and subsequent rupture of red blood cells, leading to the release of haemoglobin into the surrounding plasma. ¹⁶ Mechanically induced haemolysis has been shown to occur when red blood cells are exposed to elevated fluid stresses. ¹⁷ The haemolysis index $(HI)$ , or ratio of free haemoglobin to total blood haemoglobin $\left(\frac{\Delta Hb}{Hb}\right)$ , and its relationship with shear stress $(\tau ,\mathrm{Pa})$ , and exposure time $(t,\mathrm{s})$ , can be characterised by a power-law damage function, shown in equation (1)

$$HI=\frac{\Delta Hb}{Hb}=C{t}^{\alpha }{\tau }^{\beta }$$

(1)

where $C$ , $\alpha$ and $\beta$ are empirically derived constants obtained from in vitro experiments. Traditionally, these constants are generated from simplified experiments involving the exposure of blood to a constant shear and exposure time,^18–21 but some studies have attempted to quantify device specific constants using both in vitro and in silico approaches.^22–25 However, there are uncertainties and variabilities in the calculated constants, influenced by device design, operating condition and blood composition.

Various attempts have been made to model haemolysis in blood-contacting medical devices, where damage can either be calculated in an Eulerian or Lagrangian frame of reference. Garon and Farinas ²⁶ developed an Eulerian approach, linearising the damage function (equation (1)) in time and applying it as a source term in a transport equation for the linearised index of haemolysis $(H{I}^{\prime }=H{I}^{\frac{1}{\alpha }})$ shown in equation (2).

$$\frac{\partial H{I}^{\prime }}{\partial t}+\mathbf{u}\cdot \nabla H{I}^{\prime }=(C{\tau }^{\beta }{)}^{\frac{1}{\alpha }}$$

(2)

Equation (2) is solved using the velocity field and describes the haemolysis production within the entire domain. This method is commonly employed in simulating haemolysis in MCS devices, particularly rotary flow VADs. Fraser et al. ²² investigated haemolysis in various axial and centrifugal blood pumps, deriving device-specific empirical constants. Gil et al. ²⁷ showed superior haemocompatibility in the HeartMate 3 compared to the HVAD, while Zhang et al. ²⁸ suggested better haemocompatibility of the CH-VAD compared to the HVAD and HeartMate II under normal and hypertensive conditions. However, the Eulerian approach, as highlighted by Faghih and Sharp, ²⁹ has several limitations, most notably that it only holds true for cases with no streamwise velocity variation, something inherent to many blood pumps. ³⁰

In the Lagrangian approach, damage is assessed by integrating the shear stress along particle tracks from the inlet of the device to the outlet. ³¹ Various damage models based on this formulation, from equations (3) to (5), are summarised in Table 1. Excluded from this list is the direct integration method by Yano et al. ³² and Song et al. ³³ which has been deemed inaccurate, and the time-dependent shear loading history method by Goubergrits and Affeld, ³⁴ found to be similar to L2. ³¹

Table 1.

Summary of Lagrangian damage accumulation methods.

Model name	Source	Notes
L1 (equation (3))	Zimmer et al., 35 Lim et al. 36 and Grigioni et al. 37	Gives uniform shear stress but does not account for damage history
L2 (equation (4))	Garon and Farinas 26	Summation of the linearised damage function
L3 (equation (5))	Grigioni et al. 38	Considers a mechanical dosage, holds true for uniform shear stress cases and also addresses damage history

$$L1=\sum _{inlet}^{outlet}\alpha \cdot C\cdot t^{\alpha -1}\cdot {\tau }^{\beta }\Delta t$$

(3)

$$L2=C{\left(\sum _{inlet}^{outlet}\Delta t\cdot {\tau }^{\frac{\beta }{\alpha }}\right)}^{\alpha }$$

(4)

$$\begin{array}{c}L3={\sum }_{inlet}^{outlet}\alpha \cdot C\cdot {\left[{\sum }_{j=1}^i\tau {(t_j)}^{\frac{\beta }{\alpha }}\cdot \Delta t_j+D(t_0)\right]}^{\alpha -1}\\ \cdot \tau {(t_i)}^{\frac{\beta }{\alpha }}\cdot \Delta t_i\end{array}$$

(5)

where equation (5) was derived from the time derivative of a mechanical dose, $D$

$$D=t\cdot {\tau }^{\frac{\beta }{\alpha }}$$

(6)

The Lagrangian method is less commonly applied to MCS, possibly due to additional post-processing requirements or the simplicity of implementing the Eulerian approach. Taskin et al. ³¹ compared both Eulerian and Lagrangian methods in axial shearing devices and VADs, finding neither accurately predicted the magnitude of haemolysis, though the Eulerian method could predict relative haemolysis. Wiegmann et al. ³⁹ used a Lagrangian approach to assess rotary pump design, showing that haemolysis was negatively correlated with hydraulic efficiency.

Modelling haemolysis in positive-displacement devices is less common compared to rotary devices, with several shortcomings in the approaches that have been used. Xu et al. ⁴⁰ used a Lagrangian method to analyse the effect of pusher-plate motion profiles on haemolysis in a positive-displacement VAD, but did not consider the motion and subsequent haemolysis from the bileaflet valves, a site of significant haemolytic potential.^41,42 Richez et al. ⁴³ estimated haemolysis in the Carmat TAH using average shear values during the systolic phase, neglecting the full range of shear stresses. Bierewirtz et al. ¹³ included an Eulerian estimation of haemolysis in the Shuttlepump TAH, although in that work the method for incorporating transient geometry and flow is not investigated.

The aim of this study was to assess the effectiveness of different shear-based haemolysis models in capturing haemolytic effects within a positive-displacement MCS device. Using a computational fluid dynamics (CFD) model of the Realheart TAH, Eulerian and Lagrangian methods were compared. Utilising different empirical constants from the literature for damage estimation, the suitability of these approaches and their applicability to positive-displacement blood pumps were assessed.

Methods

Computational domain

A CFD model of the Realheart TAH V11.2 was created using an approach similar to that used for the V11c. ⁴⁴ This method used fluid-structure interactions (FSI) for the valve leaflets, enabling their movement in response to fluid forces. ⁴⁵ Both V11.2 and V11c are considered legacy devices and differ slightly in shape, but both operate using the same positive-displacement principle. In the V11.2, the outflow valve has been repositioned upstream and rotated by 90°, while the atrial inflow has been centrally relocated to enhance anatomical fit (see Figure 1(a)).

Figure 1.

Computational domain and mesh of the Realheart V112 showing (a) location of the inlet (blue), outlet (red), mitral and aortic valves (dark grey). (b) 3D rendering from micro-CT scan of ON-X bileaflet mechanical heart valve from which valve geometry was derived. (c) Valve gap sizes implemented into the CFD model of the ON-X valve, with the gap model implemented between the leaflets approximating full closure at this location. (d) Mid-plane mesh of the CFD model, showing location of the six overset zones: background, AV cylinder, two mitral valve leaflets and two aortic valve leaflets. Left view is mid-plane mesh of the closed aortic valve, whilst right view is mid-plane of closed mitral valve. The same set of views are presented for contour plots in the results.

Valve model

A 3D model of a closed ON-X bileaflet mechanical heart valve was generated using a Nikon XT 225 micro-CT scanner (Nikon Metrology, Inc. MI, USA), shown in Figure 1(b). This enabled the size of the leaflets, housing and gaps around the edges of the leaflets to be accurately measured. The average peripheral tip gap was 50 µm and the gap near the hinge mechanism was 140 µm. These measurements were incorporated into an existing CAD model of the ON-X valve using Ansys DesignModeler V2022R2 (Ansys Inc. PA, USA), as shown in Figure 1(c). This model excluded the hinge geometry, instead using a constant 140 µm side gap and limiting the leaflets to rotate between 0° (fully closed) and 44° (fully open) around a single axis. ⁴⁴ Micro-CT post-processing was achieved using Avizo 9.0 (Thermo Fisher Scientific, MA, USA).

Meshing

An overset meshing approach was used to combine the moving components into a single continuous mesh which contained a background zone, a translating AV cylinder zone, two mesh zones for the translating and rotating mitral valve leaflets and two mesh zones for the rotating aortic valve leaflets (Figure 1(d)). Ansys Fluent Meshing V2022 R2 (Ansys Inc. PA, USA) was used to mesh each individual overset zone before later combining to form the full device.

Regurgitant leakage flow through a closed mechanical heart valve has been shown to play a major role in the haemolytic potential of MCS devices.^46–49 Therefore, to create the appropriate mesh for the full device, a mesh sensitivity study was conducted in which a pressure gradient was applied to a static closed valve within a fluid cylinder, inducing leakage flow. Four meshes were created by reducing both global and local mesh sizing around the valve leaflets. Between the two leaflets, where no gap existed and both leaflets touched, the gap model function, as previously used in the V11c model, ⁴⁴ was applied to approximate a perfect seal. Bulk flow mesh parameters for the rest of the pump were the same from the validated model of v11c in our previous study. ⁴⁴ The most suitable mesh, resulting in a total of 5.1 million elements, was chosen for the full device. Further details on the mesh sensitivity study are available in the Supplemental Material.

Haemolysis models

Haemolysis was solved using Lagrangian and Eulerian approaches. For both model types, three sets of empirically derived damage function constants, $C,\alpha$ and $\beta$ were used, outlined in Table 2. Additionally, both approaches employed Bludszuweit’s scalar shear stress formulation for a single shear stress value, as shown in equation (7)

$$\tau ={\left[\frac{1}{12}{({\tau }_{ii}-{\tau }_{jj})}^2+\frac{1}{2}{\tau }_{ij}^2\right]}^{\frac{1}{2}}$$

(7)

noting that $i\ne j$

Table 2.

Damage function constants.

Model	Source	$C$	$\alpha$	$\beta$
HO	Heuser and Opitz 18 (Porcine Blood)	1.8e-8	0.765	1.991
ZT	Zhang et al. 20 (Ovine Blood)	1.228e-7	0.6606	1.9918
FZ	Fraser et al. 22 (Ovine Blood)	1.745e-8	0.7762	1.963

In the Eulerian approach, the transient form of the passive scalar transport model was solved, as shown in equation (2), to calculate $H{I}^{\prime }$ . The source term of equation (2) ( $C{\tau }^{\beta }{)}^{\frac{1}{\alpha }}$ , was applied in each of the six overset zones. The mass flow weighted average of $HI$ at the outlet, $H{I}_{\mathrm{outlet}}$ , was calculated using equation (8)

$$H{I}_{\mathrm{outlet}}=\frac{{\sum }_{i=1}^n{(H{I^{\prime }}_i)}^{\alpha }{Q}_i}{{\sum }_{i=1}^n{Q}_i}$$

(8)

where $n$ is the total elements in the outlet plane, $i$ is the current element and $Q$ is the flow rate through the element. To arrive at a single $HI$ value for comparison to the Lagrangian approaches, the time average value of $H{I}_{\mathrm{outlet}}$ was taken over the last (ninth) cycle to obtain $\overline{E}$ .

For the Lagrangian approach, 100 particles were seeded at a random location on the inlet face at every time step during the second cycle, totalling approximately 26,000 particles. After the second cycle, the injection was disabled and the particles were tracked until leaving the domain through either the inlet or outlet. $\tau$ was calculated for each particle at every time step and exported alongside the $x,y,z$ coordinate data.

The three Lagrangian formulations outlined in equations (3–5) were then calculated for each of the empirical constants shown in Table 2 for only those particles that exited the domain through the outlet. Particle post-processing was performed using MATLAB 2022b (Mathworks, MA, USA).

To calculate an average Lagrangian haemolysis from all the particles, $\overline{L}$ , each particle was weighted by its initial velocity, $v_i$ , shown in equation (9).

$$\overline{L}=\frac{\sum L_i{v}_i}{\sum v_i}$$

(9)

Collectively, three Eulerian and nine Lagrangian values of $HI$ were calculated and compared.

Boundary conditions

A constant pressure of 10 mmHg was assigned at the inlet. At the outlet, a two-element Windkessel model was applied to approximate downstream compliance, $C$ (mL/mmHg), and resistance, $R$ (Wood Units). Values of $R=22$ Wood Units and $C=0.8$  mL/mmHg were used, generating a mean pressure head of 100 mmHg. No slip conditions were used on all walls. A specified value of $H{I}^{\prime }=0$ was assigned on the inlet boundary for all Eulerian haemolysis models, and a specified flux of 0 was applied at the outlet. For the particle tracking, the inlet and outlet were set to escape boundaries, whilst all walls were set to reflect.

Numerical procedure

Ansys Fluent V2022 R2 (Ansys Inc. PA, USA) was used to solve: the Unsteady Reynolds Averaged Navier Stokes equations, with a $K-\omega$ SST model for turbulence closure, and a pressure-based coupled solution approach; the FSI motion of the valves, as detailed in Bornoff et al.^44,45 with small modifications to improve solution stability; and the Eulerian scalar transport and Lagrangian particle tracking haemolysis models. Stability improvements included reducing pressure and momentum relaxation factors to 0.35, reducing motion relaxation to 0.25 when the angle of the valve leaflet was less than 10°, and reducing the maximum time step to achieve a maximum angle change of 1° during this same phase. At all other times, pressure and momentum relaxation was 0.75, motion relaxation was 0.4 and maximum leaflet angle change per time step was 2°.

Second-order upwind discretisation schemes were used for all equations utilising a least-squares gradient method. High resolution tracking was used for the particles with a maximum number of steps of 5000 and a step length factor of 5.

Heart rate and stroke length were set to 105 bpm and 23.7 mm, respectively, based on unpublished in vitro estimates for producing a cardiac output of approximately 5 L/min.

The model was solved on a Microsoft Azure cloud HPC system operated by the University of Bath, ⁵⁰ using 120 AMD EPYC 7V73X (Milan-X) CPU cores, taking approximately 40,000 core hours to complete. Cyclic convergence was assessed by calculating the relative error between cyclic averages of $H{I}_{\mathrm{outlet}}$ for the three empirical constants, finding that nine cycles were sufficient for convergence. Shear stress data was obtained during the ninth cycle.

Results

Flow field

Figure 2 shows the key flow field results from the last (ninth) cycle. Average flow rate was 4.95 L/min, peaking at 18 L/min at mid systole. Mean pressure head was 107 mmHg, with a pulse pressure of 45 mmHg and a peak of 140 mmHg shortly before end systole. Leakage flow, driven by pressure gradients across the valves, was observed during systole (around the mitral valve) and diastole (around the aortic valve). During diastole the mean aortic leakage flow rate was 0.5 L/min, accumulating 3 mL of leakage volume per cycle. Peak leakage velocity was in the side gaps, with an approximate magnitude of 4.75 and 3.75 m/s during systole (mitral valve) and diastole (aortic valve) respectively. Peak leakage flow through the peripheral tip was 1.5 and 1.1 m/s during systole and diastole respectively. Leakage flow during systole correlated with increasing pressure gradients and mitral valve translation velocity, peaking around mid-systole, while pressure gradients peaked just before end-systole. The ninth cycle velocity and pressure contours can be seen in Supplemental Videos S1 and S2.

Figure 2.

(a) Pressure and (b) velocity flow field contour plots at mid diastole, end diastole, mid systole and end systole. Pressure contours are clipped between 0 and 150 Pa. Velocity contours clipped between 0 and 2 m/s. Black arrows indicate direction of AV plane motion, downwards during systole and upwards during diastole. (c) Volume flow rate and pressure at the outlet, showing the time points equating to mid diastole, end diastole, mid systole and end systole.

Shear stress

Figure 3(a) shows the volume-averaged scalar shear stress, with peaks labelled 1–4 corresponding to key events in the cardiac cycle. At these times, the volumetric distribution of shear stress was investigated and presented in Figure 3(b) alongside the time-averaged mean value. Higher stresses were present at valve openings (2 and 4), but over a small volume of fluid. About 99% of the volume was below 3.9, 2.9, 3.5 and 2.9 Pa for points 1–4 respectively, and 3.1 Pa for the mean value.

Figure 3.

(a) Volume average shear stress over one cycle with time points 1: mid diastole, 2: aortic valve, opening 3: mid systole, 4: mitral valve opening. (b) Percentage volume above a threshold shear stress for the four defined time points, including the mean percentage volume above a threshold shear in green. (c) Wall shear stress on the valve leaflets and pump housing, and scalar shear stress on the mid-plane of the valves at mid systole and mid diastole. Closed valves are clipped to between 0 and 1000 Pa, whilst open valves and pump housing are between 0 and 150 Pa. Scalar shear stress clipped to 0 and 10 Pa. Note the aortic valve at mid systole is viewed upside down.

Figure 3(c) shows wall shear stress on the valve leaflets and pump housing at time points 1 and 3, as well as the scalar shear stress in the bulk flow. The highest wall shear on the valves occurs in the side gaps (140 µm width) of both aortic and mitral valves, peaking at 530 Pa during diastole (aortic valve) and 540 Pa during systole (mitral valve). Shear in the peripheral tip gaps (50 µm width) were lower, with values of 430 Pa during diastole (aortic) and 440 Pa during systole (mitral). The highest stresses observed on the pump housing was around the closed valve leaflets as expected, with some elevated wall shear, approximately 30–50 Pa, in the housing around the open valve at points 1 and 3. Scalar shear stress in the bulk flow was elevated downstream of the open valves at each time point where the three-jet structure was present, with approximate values between 5 and 7 Pa.

For the open valves, the highest shear was on the leading edge of the leaflets: flow in the $+y$ direction during systole around the aortic valve and $+y$ motion of mitral leaflets during diastole. Low wall shear was observed on the inner faces of the open valves and on the upper and lower faces of the closed valves.

Eulerian haemolysis

In Figure 4(a), $H{I}_{outlet}$ increased non-linearly with time towards an asymptotic cycle average for all sets of empirical constants. $H{I}_{outlet}$ increased sharply as the aortic valve opened at the beginning of systole and fell as it closed at the beginning of diastole. Throughout systole, $H{I}_{outlet}$ remained relatively stable and there was some variation during diastole. Despite similar transient profiles for the three sets of constants, their magnitudes differed. $\overline{E}$ was $1.20\times {10}^{-6}$ , $2.28\times {10}^{-5}$ and $8.93\times {10}^{-7}$ for HO, ZT and FZ respectively, with a small error between eighth and ninth cycles (0.09%, 0.96% and 0.11% respectively), indicating convergence. The $HI$ contour for the FZ constants for the ninth cycle can be seen in Supplemental Video S3.

Figure 4.

(a) Transient Eulerian mass-weighted haemolysis index results for HO, ZT and FZ constants at the outlet. Contour plots of the Eulerian haemolysis index at (b) mid systole and (c) mid diastole for HO, ZT and FZ constants. Note that HO and FZ are clipped between 0 and 1e-5, whilst ZT is clipped between 0 and 1e-4.

The $HI$ scalar, reflecting the concentration of plasma free haemoglobin, was calculated using the scalar field of $H{I}^{\prime }$ raised to the power of $\alpha$ for each set of empirical constants and plotted for mid-systole and mid-diastole in Figure 4(b) and (c). $HI$ was greatest on the upstream side of the closed valves due to high shear leakage flow generated by pressure gradients across the valves (Figure 3(a)). During systole (Figure 4(b)), $HI$ accumulated upstream of the mitral valve and moved towards the apex of the ventricle during diastole when the valve opened. During diastole (Figure 4(c)), $HI$ collected downstream of the aortic valve in the ventricle, which was transported past the aortic valve and out of the domain during systole.

Lagrangian haemolysis

Of the total 25,800 particles that were injected, 12,172 (47.2%) escaped through the inlet, 11,027 (42.7%) escaped through the outlet and 2,601 (10.1%) remained after nine cycles, predominately in the ventricle and outflow. Only particles that travelled from the inlet to the outlet were considered for further analysis. The paths that particles took can be seen in Supplemental Video S4. Additional information regarding convergence of haemolysis results with regards to number of particles and particle time step can be found in the Supplemental Materials.

Figure 5 a plots the final $L1$ value using FZ constants for each particle against its residence time. Particles experiencing one or two instances of leakage flow are highlighted in orange and red respectively, while those not experiencing it are blue. Particles that experienced leakage flow showed the highest haemolysis levels, with 506 particles (4.6%) encountering at least one leakage event and 0.25% experiencing two events. No particles were exposed three times. The spread of $L1$ was large, with peak values more than 104 times greater than the smallest, with a similar spread observed for $L2$ and $L3$ .

Figure 5.

(a) Particle residence time against Haemolysis index using L1 formulation and FZ constants, showing no leakage exposure (blue), single leakage exposure (yellow) and two leakage exposures (red). (b) Shear stress, approximate exposure time and haemolysis index gained from leakage flow events. (c) Comparison of last cycle mean Eulerian haemolysis index and mean Lagrangian formulations for HO, ZT and FZ constants. Note the HI scale is logarithmic.

As residence time increased, $L1$ values formed vertical groups due to the pulsatile flow, where particles exited in clusters during systole. The average $L1$ of each group increased over time, but the average $L1$ for those particles that experienced leakage flow was relatively consistent from cycle to cycle. The cumulative average $\overline{L1}$ , cumulatively calculated using equation (9), converged towards a steady value as residence time increased.

Figure 5(b) displays the shear stress, exposure time and resulting accumulation of haemolysis for particles experiencing leakage flow. There was relatively consistent exposure around the edges of the aortic valve during diastole (203 Pa ± 87 Pa) and mitral valve during systole (214 Pa ± 99 Pa). Exposure time was calculated using particle velocity and gap distance of 0.85 mm, the thickness of the leaflet when closed. Results aligned with observations of the velocity field in Figure 2(b), indicating side leakage flow had higher velocity, resulting in 4–5× smaller exposure times compared to the peripheral tip gap. Results show that the mitral valve was slightly more haemolytic $(L=11.18\times {10}^{-6}\pm 1.21\times {10}^{-6})$ than the aortic valve $(L=19.59\times {10}^{-7}\pm 8.16\times {10}^{-6})$ , with relatively equal haemolysis between side and peripheral leakage for both valves.

Figure 5(c) presents $\overline{L1}$ , $\overline{L2}$ and $\overline{L3}$ from equation (9) alongside $\overline{E}$ for each set of constants. $\overline{E}$ was consistently greater than all Lagrangian formulations, followed by $\overline{L2}$ , $\overline{L3}$ then $\overline{L1}$ across all constants. The greatest gap between Eulerian and Lagrangian was with ZT constants, where $\overline{E}$ was between 50 and 4 times larger than each Lagrangian formulation. The smallest difference was with FZ constants, where $\overline{E}$ was between 12.5 and 2.7 times larger. Notably, $\overline{E}$ with ZT constants was the largest Eulerian result, 20 times larger than HO and 25 times larger than FZ.

Regarding the remaining particles within the domain, the majority were located in the ventricle and outflow. Using FZ constants and the $L1$ formulation, $\overline{L1}$ for the remaining particles was $7.9\times {10}^{-8}$ , slightly larger than the value of the escaped particles of $7.8\times {10}^{-8}$ , indicating that once these particles leave the domain, they will serve to increase the overall haemolysis, but not by a significant amount. About 74 particles (2.85%) experienced at least one leakage flow event whilst 0.35% experienced two. No remaining particles had experienced three leakage flow events by the end of the simulation.

Discussion

In this study, one Eulerian and three Lagrangian approaches to simulating haemolysis were used in conjunction with three sets of empirical damage function constants from the literature to assess haemolysis generation in a positive-displacement TAH using a CFD and FSI approach.

Most haemolysis simulations focus on continuous flow rotary blood pumps which have statistically steady temporal behaviour, although the instantaneous flow field is unsteady. This allows for steady-state assumptions which enables the simulation of haemolysis as a post-processing step. However, positive-displacement pumps produce inherently unsteady pulsatile flow, thus requiring a transient simulation and complicating the post-processing step of estimating haemolysis. In this study, Eulerian and Lagrangian haemolysis was simulated at the same time as the flow field, yielding results that asymptotically tended toward comparable single-pass cycle-steady solutions that could be compared with experimental measurements of the modified index of haemolysis: $MIH=HI\times (V/QT)\times (100-Hct)/100\times {10}^6$

Gao et al. ⁵¹ used the same approaches as Taskin et al. ³¹ to simulate haemolysis in a roller pump. Both studies showed the same order for the Lagrangian methods that was found in this study $(L>2L>3L1)$ . Additionally, Taskin et al. noted significant differences between Lagrangian and Eulerian approaches, something also observed in this study. This variance could be due to limitations highlighted by Faghih and Sharp ²⁹ and Faghih et al., ³⁰ that Eulerian models are only valid for simple steady-flow conditions and do not hold true in complex flow fields.

The variation in haemolysis with choice of empirical constants broadly followed the differences in the constants themselves. ZT’s $C$ value was roughly ten times greater than HO and FZ and was reflected in both Eulerian and Lagrangian results. Whilst the three $\beta$ constants were approximately equal, any remaining differences could be attributed to variations in $\alpha$ . A decrease in $\alpha$ appeared to increase the difference between Eulerian and Lagrangian outcomes, as HO and FZ returned similar Eulerian to Lagrangian ratios, but for ZT it was approximately twice as large. Such discrepancies could be attributed to the issue that the linearised damage function only holds true for $\alpha =1$ . ²⁹ Greater deviations from $\alpha =1$ , such as is the case with ZT constants, may increase Eulerian inaccuracies compared to the Lagrangian results. Constants derived from human blood, such as those from Ding et al., ²¹ have very small $\alpha$ values $(\alpha =0.2777)$ , indicating that Lagrangian and Eulerian results may differ even further. The constants considered in this study enabled comparisons with studies in the literature, so constants for human blood were not considered but could be used in future studies.

Bludszuweit’s equivalent scalar shear stress combined the nine components of the stress tensor into a scalar for the haemolysis source term. ⁵² This scalar shear stress was hypothesised to correlate with membrane failure in red blood cells, making it a suitable choice for using in haemolysis simulations. ⁵³ There are several other equivalent stress formulations in use, such as the von Mises, which is based on the second invariant of the viscous stress tensor, or the Tresca criterion, which considers maximum stresses.^54,55 Yu et al. ⁵⁵ investigated haemolysis using various equivalent stress formulations, and showed that the choice of empirical damage constants had a much greater impact on solution outcome than choice of equivalent stress, a result echoed by others. ⁵⁶ Additionally, whilst this study employed a turbulent approach, turbulent stresses were not considered due to the uncertainty surrounding their impact on the solution.^57,58

Another explanation for Eulerian and Lagrangian differences may come from how they captured the key haemolysis mechanism: leakage flow past closed valve leaflets. Micro-CT scans of the ON-X valve showed no b-datum gap, with 140 µm hinge side gaps and 50 µm peripheral tip gaps. This resulted in sustained larger leakage jets around the hinge region, noting that this region was not explicitly modelled, and smaller jets on the periphery. Manning et al. ⁵⁹ observed similar behaviour in vitro, supporting this modelling approach. Leakage flow was persistent throughout the cycle, generating haemolysis consistently in the Eulerian approach. However, the Lagrangian approach indicated that while leakage flow contributed most to haemolysis, high shear exposure occurred only for a few particles, resulting in discontinuous exposure events. The remaining particles after nine cycles were shown to be more haemolytic than the current average of escaped particles. They would ultimately become more haemolysed as they traversed the rest of the pump towards the outlet. However, since the majority were in the ventricle or the outflow further haemolysis would be small. Convergence with regards to the number of seeded particles was observed, as shown in the Supplemental Materials, indicating that seeding additional particles would not significantly change the calculated haemolysis value.

Despite differences between the two approaches, Taskin et al. ³¹ highlighted the ability to effectively calculate relative differences in haemolysis in different fluid-flow environments and operating condition using the shear-based Eulerian method, albeit with significant error compared to experimental values. These errors originate from the simplification of the mechanical and biological processes involved in haemolysis that are required to compute it in a simplified way. ¹⁶ Improvements to the capture of these mechanisms, such as strain-based or explicit cell modelling, will ultimately lead to more accurate haemolysis predictions, but at greater complexity and computational cost. The empirical constants in this study were chosen for their extensive use in the literature and the possibility of comparison to other blood pumps. Taking the Eulerian results with FZ constants, and comparing this study to Fraser et al., ²² the Realheart TAH potentially generates lower relative haemolysis than the four simulated VADs from that study (approximately $1\times {10}^{-4}$ in this study compared to $1\times {10}^{-4}$ – $6\times {10}^{-4}$ in Fraser et al. ²² ), noting that operating conditions may not align exactly. Differences may stem from the smaller volumes exposed to elevated shear in the Realheart TAH, limited to leakage flow gaps only, unlike VADs where elevated shear may exist throughout the pump domain. Additionally, shear stress magnitudes in the Realheart TAH vary with time, with average values (Figure 3(a)) within the same range as other studies of steady-state rotary flow blood pumps, indicating that comparisons in the Eulerian frame are possible.

Finally, the Eulerian method is better suited to visualising and understanding the underlying mechanisms of haemolysis generation and ways of optimising the design of the device. When particles are post-processed to understand their shear stress history, the results are viewed in an Eulerian frame anyway (Figure 5(b)). It is uncertain whether a particle, or enough particles, will experience the primary mechanisms for haemolysis, making optimisations of the pump challenging. Thus, the Eulerian method can be the preferred approach for tackling optimisation.

Limitations

Particles were generally exposed to leakage flow shear stress during the middle of systole and diastole when the time step was at its largest. As particle shear stress data was sampled at the same time as the flow field, there may have been insufficient temporal accuracy to capture the short leakage flow exposure of a particle. To overcome this, smaller solution time steps, or local time-stepping for each particle could be used in future studies, but at an increased computational cost. Additionally, micro-CT capture of the valve was not under pressurised conditions, possibly leading to an overestimation of leaflet gap size and subsequent amount of leakage flow. This could be rectified in future studies by imaging the valve under pressure.

Conclusion

In conclusion, a transient CFD model, utilising a previously validated FSI valve motion methodology, was employed to evaluate shear-based haemolysis models in the Realheart TAH V112. The study used Eulerian and Lagrangian methods with different empirical damage constants, finding that both methods captured the primary driver of haemolysis: leakage flow around closed valve leaflets, and that the choice of constants did not impact this capture. However, the Eulerian method, with ease of implementation and improved visualisation of damage, was more suitable for modelling haemolysis in this application. Future research will focus on validating simulated results against experiments and exploring haemolysis variations with operating conditions and device design changes.

Supplemental Material

sj-docx-1-jao-10.1177_03913988241267797 – Supplemental material for Assessment of haemolysis models for a positive-displacement total artificial heart

Supplemental material, sj-docx-1-jao-10.1177_03913988241267797 for Assessment of haemolysis models for a positive-displacement total artificial heart by Joseph Bornoff, Shaikh Faisal Zaman, Azad Najar, Thomas Finocchiaro, Ina Laura Perkins, Andrew N Cookson and Katharine H Fraser in The International Journal of Artificial Organs

Supplemental Material

sj-pdf-2-jao-10.1177_03913988241267797 – Supplemental material for Assessment of haemolysis models for a positive-displacement total artificial heart

Supplemental material, sj-pdf-2-jao-10.1177_03913988241267797 for Assessment of haemolysis models for a positive-displacement total artificial heart by Joseph Bornoff, Shaikh Faisal Zaman, Azad Najar, Thomas Finocchiaro, Ina Laura Perkins, Andrew N Cookson and Katharine H Fraser in The International Journal of Artificial Organs

Supplemental Material

sj-pdf-3-jao-10.1177_03913988241267797 – Supplemental material for Assessment of haemolysis models for a positive-displacement total artificial heart

Supplemental material, sj-pdf-3-jao-10.1177_03913988241267797 for Assessment of haemolysis models for a positive-displacement total artificial heart by Joseph Bornoff, Shaikh Faisal Zaman, Azad Najar, Thomas Finocchiaro, Ina Laura Perkins, Andrew N Cookson and Katharine H Fraser in The International Journal of Artificial Organs