An Accelerated Thrombosis Model for Computational Fluid Dynamics Simulations in Rotary Blood Pumps

Purpose Thrombosis ranks among the major complications in blood-carrying medical devices and a better understanding to influence the design related contribution to thrombosis is desirable. Over the past years many computational models of thrombosis have been developed. However, numerically cheap models able to predict localized thrombus risk in complex geometries are still lacking. The aim of the study was to develop and test a computationally efficient model for thrombus risk prediction in rotary blood pumps. Methods We used a two-stage approach to calculate thrombus risk. The first stage involves the computation of velocity and pressure fields by computational fluid dynamic simulations. At the second stage, platelet activation by mechanical and chemical stimuli was determined through species transport with an Eulerian approach. The model was compared with existing clinical data on thrombus deposition within the HeartMate II. Furthermore, an operating point and model parameter sensitivity analysis was performed. Results Our model shows good correlation (R2 > 0.93) with clinical data and identifies the bearing and outlet stator region of the HeartMate II as the location most prone to thrombus formation. The calculation of thrombus risk requires an additional 10–20 core hours of computation time. Conclusion The concentration of activated platelets can be used as a surrogate and computationally low-cost marker to determine potential risk regions of thrombus deposition in a blood pump. Relative comparisons of thrombus risk are possible even considering the intrinsic uncertainty in model parameters and operating conditions. Supplementary Information The online version contains supplementary material available at 10.1007/s13239-021-00606-y.


INTRODUCTION
Rotary blood pumps such as ventricular assist devices (VADs) and extracorporeal pumps are used in long-term and short-term mechanical circulatory support therapies for patients with acute and chronic heart failure. 13,15 Despite tremendous improvements in patient survival and overall quality of life over the last decades, one of the major device related complications is the formation of thrombi (blood clots) and the associated risk for cardioembolic stroke. 7 The causes for thrombus formation are multifactorial involving the molecular interplay of proteins and enzymes in the coagulation cascade, non-physiological blood flow conditions and interactions with foreign surfaces. 9 Ultimately the aforementioned factors enhance the hemostatic response by chronically activating platelets, the initiators of blood clots. 2 Computational modeling can help to understand the underlying complexity of thrombus formation involving a wide range of mathematical models of thrombosis spanned across different spatial and temporal scales. 6 Considering model application in rotary blood pumps in particular poses additional requirements towards thrombosis models. Computational costs are the main challenge as the complex mechanisms of thrombus formation need to be combined with the intricate flow physics within rotary blood pumps comprising turbulent flow structures, complex geometries and high rotational speeds. In addition, thrombus models should be implementable within computational fluid dynamics (CFD) solvers, the gold-standard to determine blood flow characteristics in arbitrary geometries. Optimally, such models should be compatible with commercial CFD software to allow a wider adoption throughout the scientific and engineering community.
Several approaches that serve as good candidates for model development exist in the literature. Many studies in the past adopted the idea of a platelet activation state (PAS) 22,26 that determines stress accumulation of platelet flow paths either through Eulerian or Lagrangian approaches. In consequence, each device exhibits a ''thrombogenicity fingerprint'' that can then be compared between different devices. This approach is similar to the hemolysis models existing in the literature. 34 Despite being computationally cheap, the problem to obtain spatially resolved thrombus deposition within the domain remains, with recent approaches addressing this unmet need. 4 The thrombosis model presented in Taylor et al. 27 employs conventional CFD extended with species transport of non-activated platelets, activated platelets and adenosine diphosphate (ADP) by means of additional convection diffusion equations. In addition, a Brinkman term is added to the momentum equations of the flow solution to assign areas specified as thrombus with a velocity approaching zero, thus enabling a coupling of the flow solution and the thrombus growth. The model has been applied to a simple backward facing step geometry and validated with time-resolved in vitro experiments. 28,32 For this geometry, experimental thresholds of wall shear stresses, which control the growth and detachment of thrombus areas, have been determined. However, these values could be dependent on the geometry and flow conditions. To the best of our knowledge, the model has not been applied in more complex flow geometries such as blood pumps. Furthermore, the model has been developed for laminar flows and resolves the formation time of a macroscopic thrombus which was around 30 min. This proves to be impractical for the simulation of rotating cardiovascular devices where a high time resolution is required for adequate simulation accuracy. Hence, it would be too computationally expensive to simulate the entire physical formation time of a stable thrombus.
Another thrombosis model is presented in Wu et al. 30 This model follows a similar approach as the Taylor model, but includes 10 biological and chemical species to allow a finer description of thrombosis including activation, deposition, accumulation and clearing of platelets. For the correct use of the model, 3 parameters such as the deposition rate of platelets for the materials in direct blood contact must be determined experimentally. The model has already been applied in Ref. 31 to a cardiovascular device (Heart-Mate II VAD). However, the simulations were computationally expensive (4800 core-hours-using 40 cores with 2.4 GHz) and had limited predictive power when compared to clinical studies. More specifically, thrombus formation was predicted only in one region (stator) and not in other known/common regions such as the blade area of the impeller or further downstream locations such as the bearing and diffusor, see Besides these models there exist many studies using the flow field results to identify regions prone to thrombus formation such as recirculation regions, 5,20 low velocities and low washout 25 or high residence times. 18 However, all of these studies do not model thrombosis (or parts of thrombosis) explicitly. In addition, there are other approaches using multiscale simulations 10 and other computational methods besides CFD. 35 The interested reader is referred to the extensive overviews existing in the literature. 6,33 In summary, there is still a lack of computational models capable of predicting locally resolved thrombus formation in complex geometries with reasonable computational effort.
The aim of this study is to develop such a model for rotary blood pumps focusing on low computational costs and ease of implementation in a state-of-the-art CFD solver. Moreover, information about parameter sensitivity of all mechanisms important for the computational modeler and a ready-to-use setup are provided.

Accelerated Thrombosis Model
The accelerated thrombosis model is closely based on the thrombosis model presented by Taylor et al. 27 with several reductions and adjustments to save computational time. A graphical overview of our model is shown in the Supplementary Fig. S1. Detailed information about the Taylor model can be taken from the original paper.
In brief, our approach uses three convection-diffusion equations describing the species transport of nonactivated platelets (NP), activated plates (AP) and adenosine diphosphate (ADP).
In these equations D x is the respective diffusivity coefficient, A M the mechanical activation, A C the chemical activation and R ADP the amount of ADP contained in a platelet.
Platelet activation (change from NP to AP) occurs either with mechanical cues through a shear stress based power law model, 26 Eq. (4), or with chemical cues through a threshold based activation triggered by the ADP concentration, Eq. (5).
u f is the ratio of activated platelets to total platelet count. a, b and C are the power law model parameters taken from Ref. 26. ADP t is the threshold for chemical platelet activation and t ADP is the characteristic time of the platelet activation rate. Finally, s is the scalar shear stress calculated according to Eq. (6) with r ij as the individual components of the shear stress tensor.
Both activation cues are summed up to determine the rate of platelet activation that is included in the according convection-diffusion equations as a source component. The AP species is considered as the model output and represents the thrombus risk.
In addition, the implementation of the viscous stress tensor presented so far was also compared with an implementation of additional Reynolds stress components. 29 Equations (7) and (8) show the definition of the total shear stress tensor with the viscous stress tensor part on the left side and the Reynolds stress part on the right side of the minus sign. Àq In these equations l is the dynamic viscosity, l t the eddy viscosity, q the density of the fluid, k the turbulent kinetic energy, d ij the Kronecker delta, u i the velocity components and x i the spatial components.
In comparison to the Taylor model we did not include a sink in the momentum equations through a Brinkman term which reduces the velocity of regions with high AP concentrations. Moreover, the accelerated thrombosis model uses an uncoupled approach by first determining the velocity field within the pump through a (steady state or last time step of transient) CFD simulation and then, in a second calculation, resolving the species transport based on the output of the first simulation. For the results shown in this manuscript, steady state simulations were always used as initialization for the species transport simulations. The implemented model in the FDA benchmark blood pump 19 together with a tutorial on how to run the simulations can be found at https://doi.org/10.5281/ze nodo.5116063.
An overview of all parameters used in the accelerated thrombosis model are given in Table 1. In addition to the formula symbols, the definition and the according equation, this table also contains information on how the values of the parameters were determined and whether the parameters were part of the parameter study presented later on. Finally the values of the parameters are also shown. The sources of the parameters can be found in Taylor et al. 27 With the underlying information for R ADP found in Ref. 16 and for ADP t found in Ref. 11. In the column sensitivity assessment, ''not complete'' means that this variable was varied but no complete sensitivity analysis was carried out.

Geometries, Mesh and Simulation Parameters
A reverse engineered version of the HeartMate II (HM2) geometry, an axial flow ventricular assist device (VAD), was used. The pump geometry was meshed with ANSYS 2021 R1 Meshing tool (ANSYS Inc., Canonsburg, USA) using unstructured tetrahedral elements and prism layers to assure a y+ value < 1. The bearings were modelled as hemispherical domains of 0.05 mm thickness and a bearing shell diameter of 3 mm, see Supplementary Fig. S3, and meshed using the sweep mesh method with a segmentation of 20. As can be seen in Fig. S2, a mesh independence study resulted in mesh element numbers of approx. 12 million elements for the HM2 geometry. The non-Newtonian blood model by Ballyk et al. 1 with a density of q = 1056.4 kg m À3 (hematocrit level of 44%) was used. The k-x Shear Stress Transport model was taken as the turbulence model. For the steady-state initialization simulations, the auto-time step method was used with a time step of around 0.1 milliseconds, and the stabilization of the pressure difference between the inlet and outlet of the simulation was used as the termination criterion. This took about 0.15 simulated seconds which correspond to 2 times of the pass-through time from inlet to outlet. For the simulations of the accelerated thrombosis model, a physical time step of one second was used and the stabilization of the volume-averaged AP concentration in each simulation domain was used as the termination criterion.
The computational mesh of the pump is shown in Fig. 1 with blue as the parts with a rotating-wall boundary condition and grey as the parts with a nonrotating-wall boundary condition.The baseline simulations for the HM2 were performed at 9000 rpm and a flow of 2 L/min. ANSYS CFX was used as the CFD solver.

Comparison with Clinical Data
In order to compare the computational model with in vivo data on thrombus occurrence in a blood pump, the study by Rowlands et al. 23 was used. In this study, a retrospective analysis of 29 explanted HM2 with suspected pump thrombosis was performed to identify the frequency, composition and localization of macroscopic thrombi. The locations used in Ref. 23 are named Inflow Straightener, Fore Bearing, Impeller Blade-Blade and Outlet Bearing/Stator. The corresponding regions in the simulation are named Re-gion1-4. It has to be noted, that in Ref. 23 another region upstream of the inflow straightener, the fore connector seam, was considered. As this region was not present in our HM2 geometry, it could not be included in the analysis. Most thrombi were found in the Outlet bearing/stator region (region 4). Furthermore, due to the different locations and frequencies of the thrombi, a hypothesis was made that thrombi were probably transported downstream and did not necessarily originate where they were found. Because there was no information on the respective operating points from the explanted pumps, the baseline 9000 rpm at 2 L/min according to Ref. 31 was assumed and later expanded for more speeds and pump flows to cover a broader operating range. Statistical comparison between in vivo and in silico data was performed by linear regression with Python. The output was the coefficient of determination (R 2 ) and the according p value.

Operating Point and Sensitivity Analysis
In order to judge the usability of the thrombosis model, the influence of operating point variation and the model sensitivity to changes in model parameters need to be assessed.
An operating point analysis for the HM2 was performed to determine the influence of speed variation (7000-11,000 rpm) and the flow rate (1-4 L/min) on the model output. Sensitivity analysis was performed for the shear stress power law parameters C, a and b (± 10% variation and ± 20% for the most sensitive parameter) and the initial values of background activation u a and u b (ratio between activated and nonactivated platelets of 1, 5, 10 and 20%). As preliminary simulations did not deliver a high contribution of chemical activation, the parameters ADP t and t ADP were only adjusted once to identify their influence on the results. Lastly, the influence of the Reynolds stresses (latter part of the right side of Eq. (7) on the platelet activation was evaluated.

Model Results and Comparison with Clinical Data
The model results with two visualizations of the AP variable are shown in Fig. 2. For better readability the AP species was scaled such that the per mille increase in relation to the initial value of 25e12 (m À3 ) was determined and this scaling was kept consistent throughout the manuscript. In Fig. 2a the geometry is described by blue dots and blood volumes of scaled AP concentration above 5.7 are shown in red. Small volumes of elevated scaled AP concentration are located behind the front stator and along the blade edges. Larger volumes can be seen in the rear stator. In Fig. 2b the scaled value of the AP variable is displayed on the surface by a contour plot. The scale is chosen such that differences on the surface are clearly visible.
To ensure comparability of the results, the same scale is used throughout the manuscript unless otherwise specified. On the surface, the areas of the trailing edges of the front stator, parts of the blade passage and especially the areas after the bladed area can be identified with high AP values. Figure 3a shows the volume averaged scaled AP concentration and Fig. 3b shows the frequency of the thrombi as identified from the study 23 in the respective fluid regions. Both data suggest region 4, the rear stator and rear bearing as the location most prone for thrombus formation. An excellent correlation between model output and observed thrombi exists, (R 2 = 0.99, p value = 5.4eÀ5).

Operating Point Analysis
The influence of speed variation at constant flow rate of 2 L/min on the AP concentration is shown in Fig. 4a). The data shown in black corresponds to the reference operating point of 9000 revolutions per minute and 2 L/min. As the speed increases, the concentration of AP in each region also increases. The level of AP concentration in region 3 approaches the concentration level of region 2 as the speed increases, meaning that the platelet activation potential of the bladed impeller area approaches the potential of the bearing region. In general, there is an approximately linear relationship between speed and AP concentration.     of AP can already be seen in region 1 (inflow straightener) and in region 4 (outlet bearing/stator), which exhibits the highest AP concentration of all regions. Figure 5a) shows the surface distribution of the AP variable for 1, 2 and 4 L/min flow at 9000 rpm and Fig. 5b) the volume renderings of scaled AP concentrations higher than 5.7. The different flow conditions, which result from the different operating points, also result in different distributions of the AP concentration on the surface. In the case of 1 L/min, almost the entire wake region of the impeller has high AP concentration values. A comparison between the 2 L/min scenario and the 4 L/min scenario exhibits differences in AP concentration at the trailing edge of the first stator. In addition, the shape of the AP distribution in the bladed passage changes and the concentration of AP in the rear area is lower in the 4 L/min than in the 2 L/min case. A common feature of all simulations is that the two bearing areas are always indicated with a high AP concentration.
The correlation between the different operating points and the clinical data is shown in the Supplementary Fig. S4. The lowest R 2 for all scenarios was 0.93, p value = 0.035.

Shear Stress Power Law Model
The percentage deviation of the AP concentration after parameter variation of a, b and C for the baseline operating condition is shown in Fig. 6a). The front and back bearings are presented as separate regions. A change in the parameter a results in a high percentage deviation of the AP concentration, especially in the bearing areas. The parameters b and C do not seem to have a major influence on the results. In Fig. 6b), the correlation between simulation and clinical data for the different perturbations of a is shown. Although the lines differ in intercept and slope, they all show a good correlation accuracy with the study. 23 As observable in the Supplementary Fig. S5, changes of a influence the absolute values of AP concentration, however the spatial distribution of high and low AP concentrations remains the same.

Influence of Turbulence Modeling
Whether the scalar shear stress is obtained with the viscous stress tensor or with the additional Reynolds stress tensor also has an influence on the mechanical activation and the results are shown in Fig. 7. Reynolds stresses increase platelet activation, especially in the region downstream of the impeller blades (region 4). Contour plots of the AP concentration throughout the pump are found in the Supplementary Fig. S6. The qualitative behavior of the AP concentration throughout the different regions does not change.

Influence of Different Background Activation Levels
The influence of the ratio of activated platelets to non-activated platelets (ratios of 1, 5, 10 and 20%-total concentration of 5e8 platelets/ml) as an initial condition was investigated. The variations delivered AP concentration levels spanned over two orders of magnitude ranging from 5.02e6 platelets/mL to 1e8 platelets/mL, however, did not change the relative AP concentration differences between the regions1-4, as can be seen in Fig. S7.

Influence of Chemical Activation
A closer look at the baseline simulation with 9000 rpm and 2 L/min has shown that the ADP threshold value responsible for the chemical platelet activation was never reached. The highest value within the domain was 1.53eÀ4 (mol m À3 ) which is far below the ADP t value of 2eÀ3 (mol m À3 ). Running the model without any chemical activation thus produced the exact same results. Reducing the ADP t to 1eÀ4 (mol m À3 ) led to model collapse due to very high source terms of the convection-diffusion equation and a negative volume fraction u f = u a /(u a + u n ). To circumvent this problem, we increased the characteristic time of chemical activation t ADP from 1s to 100s. The simulation produced an extreme increase of activated platelets in region 2 (fore bearing) and did not affect the other regions, see Supplementary Fig. S8a and b. A separate analysis of the source terms A M and A C from Eq. (1), as shown in Fig. S8c, reveals a particularly high value for chemical activation in the fore bearing for this case.

DISCUSSION
The aim of this work was to develop a computationally efficient thrombus model in rotary blood pumps that can be easily implemented and adapted with current CFD solvers.
The model is based on the work of Taylor et al. 27 considering mechanical and chemical activation of platelets and convection-diffusion transport of three species: non-activated platelets, activated platelets and ADP. There are three main differences distinguishing our approach from the model of Taylor   solve the convection-diffusion equations resulting in low additional computational costs between 10 and 20 core-hours. Second, the formation of a clot and its influence on the velocity field is excluded for our model and thrombus risk is considered as the concentration of the activated platelet species AP. Lastly, in our model the shear stress tensor including viscous stresses and Reynolds stresses is considered for the mechanical activation (Eq. 3). Our model correlated well (R 2 > 0.93, p < 0.035) with the study by Rowlands et al. 23 even considering the fact that the operating conditions were not exactly known in the clinical study or were probably different for each patient.
Hence, analyses of different speeds between 7000 and 11,000 rpm with constant flow of 2 L/min and different flows between 1 and 4 L/min with constant speed of 9000 rpm were carried out. The correlation between pressure, flow and speed that occurs in reality is not taken into account in order to better investigate the influence of the individual parameters.
Higher rotational speeds cause higher concentrations of AP. The increase in regions downstream of the impeller was greater than in regions upstream of the impeller. This observation can be explained by the fact that the majority of mechanical activation takes place within the impeller region, with higher shear stresses at higher speeds and the activated platelets accumulate downstream in recirculation regions around the diffusor. In low flow simulations, considerably higher concentrations of activated platelets were predicted than in high flow simulations. This is also consistent with the observations of Wu et al. 31 Low flow operating conditions have also shown to lead to increased hemolysis both in silico 14 and in vitro 24 due to higher recirculation and longer exposure times of red blood cells. Probably a similar effect takes place for the platelets.
Influencing factors of the different components of the accelerated thrombosis model were investigated. The sensitivity analysis of the parameters included in the shear stress power law model showed that only the parameter a had a large influence on the activated platelet concentration. This influence, however, applied equally to the entire domain, thus raising the overall activated platelets concentration level to a similar extent. Hence, the relative contribution of different pump compartments towards thrombus formation stayed the same. In addition, the influence of turbulence effects on the scalar shear stress was investigated. Compared to the viscous stress approach, the activated platelet concentration increased especially downstream of the impeller blades, a region where one would expect a high contribution of turbulent stresses. Of note, a Reynolds-averaged Navier-Stokes equations (RANS) turbulence model was used. Scale resolving models such as Large Eddy or Detached Eddy Simulations might produce different contributions of turbulent stresses, potentially higher than the ones of the RANS model. Another possibility would be to resolve the unresolved turbulence by a turbulent energy dissipation approach. What effect such a small scale resolution of turbulence has on our Eulerian approach could be investigated in further studies. Changes of the initial ratio of activated to nonactivated platelets raised the absolute value of AP but did not affect the relative contribution of different pump components. Lastly, it was shown that chemical activation did not have any influence on the results as the activation threshold value ADP t was never reached. Lowering this threshold value (and increasing the t ADP from 1 to 100 s to achieve stable simulations) led to a drastic increase of activated platelets in the fore bearing region (region 2) but did not affect other parts of the pump. Such a ''switch-like'' behavior is one of the main weaknesses of the model and a continuous relationship between ADP concentration and activated platelets, e.g., through a Hill-like function, should be implemented in the future. A separate investigation of the source terms Am and Ac from Eq. (1) has demonstrated a high proportion of chemical activation in the front bearing in this case. From investigation of the source terms Am and Ac, it can be concluded where and in what intensity activation occurs. This approach provides the user of the model not only with information about the location of thrombus deposition but also with information about the location of activation.
Our approach compares to the existing thrombogenicity models in the following ways. As outlined in the introduction, PAS, which is easy to implement, is often taken as a global measure to compare different designs with each other and does not focus on particular regions of the devices. 12 A recent implementation by Chiu et al. 4 used a particle approach and computed a probability density function from the stress accumulation of particles, compared two VAD designs with each other and also focused on the probability density functions in the different regions of the pumps. The global in silico comparison (device 1 vs. device 2) was confirmed with in vitro and in vivo experiments, however it remains unclear to what degree the localized thrombus risk can be trusted upon. The particle based approach has, in general, limitations regarding correct stress accumulation in narrow gaps of the pump, an observation that has also been made for hemolysis models. 34 The model by Wu et al. 30 provides more information on thrombus mechanisms by consideration of 10 instead of just 3 biochemical species. It clearly gives a more thorough description of thrombosis than our model at a cost of higher parameter uncertainty (3 parameters such as the deposition rate of platelets for the materials in blood contact must be determined experimentally) and higher computational costs of 200 core-hours vs. 20 core-hours.
A recent third modelling approach by Dai et al. 8 employed a two-phase approach of two fluids with identical rheological properties and defined thrombosis risk regions as regions of low washout. Although computationally cheap, the approach is purely phenomenological and depends on definitions of ''washout'' regions. Even if such thresholds are identified for one pump geometry (with existing in vitro data), there is no guarantee that these values will hold for other device designs. Lastly, the model of Taylor et al. 27 was developed and experimentally validated in a backward facing step and, as already explained in the introduction, would exhibit extremely high computational costs when implemented in a rotary blood pump environment.
All things considered, we believe that our model provides a promising trade-off between complexity and computational costs. It includes a higher degree of mechanistic detail compared to the PAS approach but does not include as many mechanisms as the model by Wu et al. 29 The approach of uncoupling the flow and the species transport together with neglecting a physical presence of a thrombus drastically reduces computational costs but still allows a reliable comparison with clinical data. Using a steady state solution of the velocity field while neglecting dynamic flow features potentially works for continuous flow pumps. New generations of VADs include additional pulsatility, such as the Lavare cycle of the HVAD and the artificial pulse of the HeartMate III. 3 To overcome this limitation of our model and to consider the transient influence of pulsatility one could neglect the uncoupling and combine flow field calculations with the species transport. However, this would again lead to long computation times. An approach in which different representative operating points of a pulsatile cycle are calculated in a stationary and uncoupled manner (as shown in this manuscript) would also be conceivable when the results of AP concentration are evaluated individually per operating point or superimposed by means of a heat map. As mentioned above, a more detailed representation of chemical activation (Hill-like activation and consideration of further chemical species) should be a next possible step. In addition, the influence of the coagulation cascade and the influence of artificial materials on platelet activation could be included, like in the model by Mendez et al. 21 Validation of a numerical thrombosis model remains a tremendous challenge, as it is extremely difficult to obtain experimental data on thrombus formation in medical devices. Hence, computational studies about thrombus formation should be handled with care. To the best of our knowledge, the study by Rowlands et al. 23 is the only study that provides quantitatively reliable data on thrombus deposition. Despite its low sample number of n = 29, the observed distribution of thrombi is significantly different from a distribution by chance (p = 0.005, Pearson's v 2 test). Thus, we believe that this data can be used as ''ground truth''. It still has to be tested if the model performs equally well in other pump geometries and if comparable risk thresholds are identified. Our exploration of the Medos DP3 centrifugal blood pump, see Supplementary Fig. S10, showed increased concentration of activated platelets in the bearing region and in the secondary gap which seems to correspond with clinical observations of pump head thrombosis of the DP3. 17 Nevertheless, it is not appropriate to consider this observation a validation of the computational model and more research is needed.

CONCLUSION
Within this study we presented a computationally efficient model to predict thrombosis risk for CFD simulations in rotary blood pumps. The model exhibits equal or better accuracy and provides lower computational costs than common state of the art approaches. The presented model was investigated on the one hand for operating point uncertainties and on the other hand for the following model parameter uncertainties: Shear stress power law model Influence of turbulence modeling Influence of different background activation levels

Influence of chemical activation
The main conclusion that can be drawn from the investigations is that relative comparisons of thrombus risk regions are feasible even considering the uncertainty of operating points or variations of the model parameters. Due to the lack of experimental data for thrombus deposition in blood pumps, the model could only be compared on a geometry with 4 data points. Therefore, we encourage the scientific community to implement our model and compare it on additional data points. The implementation of the model in the FDA Round Robin geometry within the commercial solver ANSYS CFX has been in-cluded under (https://doi.org/10.5281/zenodo.511606 3) for further testing and expansion.

ACKNOWLEDGMENTS
The authors thank Dr. Bente Thamsen for the reconstructed HeartMate II geometry.

AUTHOR CONTRIBUTIONS
All authors contributed to the study conception and design. CB developed the numerical model, performed the simulations, gathered, analysed and discussed the results. SGH, MN and US were involved in the analysis and discussion of the results. MN supervised the project. MN and CB wrote the manuscript based on the input of all co-authors. All co-authors read and approved the final version of the manuscript.

FUNDING
Open Access funding enabled and organized by Projekt DEAL. This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors. Open access funding enabled and organized by Projekt DEAL.

DATA AVAILABILITY
The raw data can be retrieved by request from the authors.

CODE AVAILABILITY
The implementation of the thrombus model in the FDA benchmark blood pump geometry is available on https://doi.org/10.5281/zenodo.5116063.