Numerical modeling of blood damage: current status, challenges and future prospects

Numerical tools including computer-aided design, finite element methods and computational fluid dynamics (CFD) programs are now widely used for the development, optimization and investigation of artificial organs. The impact of these tools has increased, especially in terms of the development of long-term (>5 years) artificial organs, such as prosthetic heart valves or ventricular assist devices.

Hemolysis, thromboembolic complications and internal bleeding as a result of the consumption of anticoagulants continue to present complications in the development of artificial organs [1–3]. Blood damage reduction would result in the decline of anticoagulant consumption, which in turn would improve the performance of artificial organs. This is the main goal in the development of the new generation devices – the functionality of which has already been proven in clinical practice.

Hemolysis and platelet lysis or platelet activation combined with thrombus formation are processes that are directly related to the hydrodynamic situation in artificial organs, as shown by investigations [4,5]. The draft for combined international (ISO 5840) and European standards (EN 12006–1) for cardiovascular implants, such as cardiac valve prostheses, recommends the reduction of turbulence in order to reduce hemolysis and avoid stagnation, and the introduction of recirculation regions in order to reduce thromboembolic complications. The technical committees, ISO/TC 150 ‘Implants for Surgery’ and CEN/TC 285 ‘Nonactive Surgical Implants’, who prepared this draft, also recommend using the experimental methods of fluid mechanics – particle image velocimetry or laser Doppler velocimetry – and/or CFD for the study of the hemodynamic situation. However, they also recommend the validation of the CFD results by experimental methods.

In fact, the numerical study of the 3D flow field using CFD tools is unnecessary, except for a deeper understanding of the hemodynamic situation, for the design optimization of artificial organs. The 3D flow field information will be necessary for numerically predicting the blood damage and thrombus formation in artificial organs when it becomes possible to do so. Normally, artificial organs are tested in vitro before evaluation in animal experiments and clinical tests. These tests also include an estimation of blood damage potency, usually using human or animal blood as a test fluid [6,7]. However, the challenges associated with handling blood and blood variability makes the results of these tests problematic to interpret and to compare [8,9]. For these reasons, it would be advantageous to have a numerical method for estimating the blood damage in a computer model. In order to use the CFD results, we have to know the relationship between the flow field and the parameters describing the blood damage. In other words, some cost function is necessary.

In order to quantify the blood damage caused by an artificial organ, some nondimensional parameters were introduced. The results of the damaging effects on erythrocytes are given by the hemolysis index (HI) [5], which is built with the total hemoglobin concentration (Hb; mg/100 ml) and the increase of the free plasma hemoglobin concentration (dPHb; mg/100 ml) after shear loading: HI = 100·dPHb/Hb (%) (1)

In some works a normalized index of the hemolysis (NIH) is used [10]: NIH = V·(100 – Hkt)·dPHb/100·F·t (2)

where dPHB (mg/100 ml) is the increase of free plasma hemoglobin concentration during time t (min), V is the circuit volume (l), Hkt is the hematocrit (%) and F is the flow rate (l/min).

The results of the damaging effects on platelets are given by the platelet lysis index L_pl: L_pl = 100·dPl/Pl (%]) (3)

where dPl is the change in total platelet concentration (Pl) [11]. The indices, introduced in equations 1–3 may be used as a cost of the blood damage cost function and is generally referred to later as the damage index (DI).

After the problem of blood damage was recognized, attempts were made to elucidate its mechanism [4–8,12,13]. Primarily, the loading stress applied, τ (laminar shear stress and turbulent or Reynolds shear stress) and the exposure time, t, influence the amount of hydrodynamically caused blood damage. These are the main parameters of the cost function DI = f(τ, t), which has to be found.

A large number of blood damage investigations were performed with blood samples or single cells [4,6–8,13–17]. These known experimental results of the shear-induced blood damage allow the definition of the following requirements for the physically consistent mathematical model of blood damage:

•	Principle of causality, which does not allow a reduction of the DI due to decreasing shear stress or, in other words, always requires a positive production of the DI: dDI(t)/dt > 0 (4)
•	Reproduction of known experimental data for uniform mechanical loading, which is described by the power law: DI(t) = A·tα·τβ (5)
•	Take into account the dependence of the DI from the load history for time-dependent loading and reproduce the experimental results of Yeleswarapu [17]: dDI(t)/dt = f(τ[0...t]) (6)

There is a broad agreement among scientists about the requirements mentioned above [18–20]. Considering equation 5, we have to note that the power law relationship between blood damage, exposure time and shear stress may be found only in a few works [11,21]. The major part of the experimental works describe the relationship between the critical value of the shear stress and the exposure time which is a so-called ‘all or nothing’ blood damage process (see Table 1, equations 9 and 10). However, this relationship cannot be used as the desired cost function, since no cost is defined by this approach. Table 1 gives an overview of the known experimental data composed into the relationships for erythrocyte and platelet lysis for both approaches.

Most of the published works [14,17,19–22] about the numerical estimation of blood damage use the power laws 7 and 8 proposed by Giersiepen [11] as the result of a 2D regression analysis based on the data of Wurzinger [4]. This equation is used in spite of the known fact of the overestimated blood damage that was generated in Wurzinger’s experiments owing to heat produced by a used seal [16]. The reason is that the power constants in this equation are nearly correct, since they reflect the mechanical properties of the cell membrane, whereas the constant A is overestimated and should be corrected [19].

Challenge one of blood damage modeling

The main challenge in using the power law for blood damage modeling is a nonlinear dependence of DI on the time (α < 1), which is a direct consequence of load history dependence. Many researchers neglect this important phenomenon and simply integrate the damage index over all cells in the numerical space (Eulers formulation) or integrate along the path lines after their discretization (Lagrangian formulation) [22,23]: DI(t+dt) = DI(t) + A·(dt)^α·τ(t+dt)^β (11)

For the past few years, the use of the Lagrangian formulation has been dominant in the mathematical models of blood damage. One of the possible solutions for considering the blood history was proposed by the author [19]. The new method fulfilled all of the above mentioned requirements. The method is based on the integration of the blood damage along the path lines. In order to consider the time dependent load history, some additional virtual time is introduced for each time step of the discretized path line. The calculation of this virtual time is based on the assumption that the existing damage index is the only parameter that describes the loading history. A similar approach has also been proposed by Grigioni and colleagues [20,24]. This new approach results in the following modification of equation 11: DI(t+dt) = A·(dt+(DI(t)/(A·τ(t+dt)^β))^1/α)^α·τ(t+dt)^β (12)

Once the problem of developing a model to satisfy all the above requirements have been solved, the next problem arises.

Challenge two of blood damage modeling

How can we validate the mathematical model of blood damage? Until now, this has been an unsolved problem. The same problems with blood handling and variability during in vitro tests, which forced us to develop a mathematical model of blood damage, impede the validation of this developed model. In order to solve this problem, we have to consider it from different points of view and develop new approaches.

One approach would be to apply equations 7–10 to the known lifespan of blood cells in the human circulatory system (120 days for erythrocytes and 14 days for platelets). As has already been shown in [19], the correction of coefficient A in equation 7 by a factor of 0.0416 (A_corr = A/24) adopts this equation to the cells’ lifespan. This correction also correlates with the known fact of the overestimated blood damage in Wurzinger’s experiments and also eliminates the discrepancy between power law and threshold approaches. In order to be able to apply this approach fully, new in vivo tests are needed. These experiments have to, for example, define the blood cells’ lifespan in patients with different prosthetic heart valves. Proved techniques for such experiments exist. These may be one of the following: a classical chromium 51 (51Cr)-labeling; a biotin labeling of the cells that does not expose the patient to radiation; or a measurement of the erythrocyte creatine, a marker of erythrocyte age that increases with shortening erythrocyte survival [25,26].

Another approach should combine the threshold with the power law approach, which would allow some cross correlation between the two approaches. Note that the threshold approach represented in equations 9 and 10 is based on a significantly larger experimental data basis than the power law. According to equation 9, the critical value of shear stress is defined as 112 Pa for the exposure time of 0.5 s. Substituting these values into equation 7 would result in the following: the HI may not be higher than 0.07% after one passage of the blood through the artificial organ. This is a very high value, since the estimated HI after one passage through the blood system is only 0.00058% [19].

A third approach should include clinical practice. For example, the efficiency of the antithrombotic therapy with anticoagulants is measured now by an international normalized ratio (INR). According to the formula, INR = (patient prothrombin time/mean normal prothrombin time)(^ISI). The international sensitivity index (ISI) defines the responsiveness of each laboratorys’ thromboplastin reagent to that of a reference reagent from the WHO, which corresponds to ISI = 1 [27]. Clinical experience showed that patients with either a caged-ball valve or more than one mechanical prosthetic valve have the lowest incidence of thromboembolic and hemorrhagic events when the INR is between 4.0 and 4.9. For patients with mono-leaflet tilting disk valves the optimal INR is between 3.0 and 3.9, and 2.0–2.9 in those patients with bileaflet disk valves [28,29]. The same parameter is optimal from 2.8 to 3.0 in patients with the rotary blood pump INCOR^® (Berlin Heart AG, Germany). It is imaginable to use the INR as a cost for the blood damage cost function in order to validate mathematical models of the blood damage instead of using in vitro tests with blood.

In addition to the above-mentioned new approaches for blood damage model validation, some new experimental data are necessary.

Challenge three of blood damage modeling

The next challenge of the application of the power law for blood damage modeling is the unknown difference between the impact of the laminar shear stress and the turbulent or Reynolds shear stress on the DI. In fact, equations 7 and 8 are the result of the regression analysis of tests with only laminar shear stress, whereas equations 9 and 10 were based on the regression analysis of the experiments with laminar and turbulent shear stresses. Note that the results for critical turbulent shear stresses were very controversial. The works of Forstrom, Indeglia and Blackshear defined the critical Reynolds shear stress as 4000 Pa or ten-times higher than for laminar shear stress [15,21,30,31]. Sallam found the critical value to be approximately 400 Pa or the same as for laminar shear stress [32], this is also supported by Sutera [33]. Recently published data from Kameneva also does not improve the situation, since she found that turbulent wall shear stress produces approximately five-times more hemolysis than laminar wall shear stress [34].

The future prospects of the blood damage mathematical model of becoming a broadly used, capable tool for the artificial organs’ development are directly dependent on our ability to overcome the above-mentioned challenges. Recapitulated, they are the following: modifying the power law model in order to fulfill all requirements of blood damage modeling, elaborating new approaches for the blood model validation and acquiring, experimentally, some additional knowledge.

Furthermore, the Eulers’ formulation should be used instead of the Lagrangian in future models, since Eulers’ approach is more convenient for the visualization and analysis of the results, has the same resolution as the CFD solution and equally considers all regions of the numerical space. This automatically means that we have to solve the mass transport equation for the cost, including convection, diffusion and the nonlinear time-dependent source term S. In this case, some physical variable P, such as free released hemoglobin or activated platelets, should be used as a cost. This should be done in order to correctly define the diffusion constant D. Equation 13 shows the desired mass transport equation:

Challenge four of blood damage modeling

This challenge is also related to the problem of blood damage mathematical model validation and is a consequence of blood damage loading history dependence. It is not possible to validate the numerical modeling of blood damage through tests with blood without knowledge of the DI (HI or L_pl) before the in vitro test has begun. Currently, we cannot define the damage index of blood at the beginning of the in vitro test. Consequently, we begin the numerical modeling of blood damage with a zero DI. As a result, this prevents us from comparing experiments with the numerical results since the initial conditions are not known. Only a tendency may be predicted using this approach, which seems to be the maximum that may be achieved now from the numerical modeling of blood damage. However, this is not enough for numerical blood damage modeling to become really successful at this point in its development.

Table 1. Summarized experimental data of the blood damage for two main approaches: power law and threshold.

Approach	HI	Lpl
Power law approach	HI = 3.62·10^-5·t^0.785·τ^2.416 (7)	Lpl = 3.66·10^-6·t^0.77·τ^3.075 (8)
Threshold approach	τcrit = 88.905·t^-0.3372 (9)	τcrit = 54.986·t^-0.3585 (10)

HI: Hemolysis index; L_pl: Platelet lysis index.