A fibrin enhanced thrombosis model for medical devices operating at low shear regimes or large surface areas

Over the past decade, much of the development of computational models of device-related thrombosis has focused on platelet activity. While those models have been successful in predicting thrombus formation in medical devices operating at high shear rates (> 5000 s−1), they cannot be directly applied to low-shear devices, such as blood oxygenators and catheters, where emerging information suggest that fibrin formation is the predominant mechanism of clotting and platelet activity plays a secondary role. In the current work, we augment an existing platelet-based model of thrombosis with a partial model of the coagulation cascade that includes contact activation of factor XII and fibrin production. To calibrate the model, we simulate a backward-facing-step flow channel that has been extensively characterized in-vitro. Next, we perform blood perfusion experiments through a microfluidic chamber mimicking a hollow fiber membrane oxygenator and validate the model against these observations. The simulation results closely match the time evolution of the thrombus height and length in the backward-facing-step experiment. Application of the model to the microfluidic hollow fiber bundle chamber capture both gross features such as the increasing clotting trend towards the outlet of the chamber, as well as finer local features such as the structure of fibrin around individual hollow fibers. Our results are in line with recent findings that suggest fibrin production, through contact activation of factor XII, drives the thrombus formation in medical devices operating at low shear rates with large surface area to volume ratios.

Hemostasis is the ensemble of biochemical and biophysical processes that control blood 2 clotting when a blood vessel is damaged. These processes work in a delicate balance to 3 stop bleeding yet prevent excessive thrombus formation (thrombosis) [1]. Medical 4 devices that contact blood disrupt the physiological hemostatic balance leading to 5 thrombotic or hemorrhagic complications. This may occur in two ways. First, platelets 6 are susceptible to supra-physiological shear stresses induced by blood flow in devices 7 such as prosthetic heart valves (PHV) or ventricular assist devices (VAD) [2][3][4]. When 8 stressed, platelets become activated releasing platelet agonists that amplify synergistic 9 activity in a positive feed-back loop. Second, when blood is in contact with foreign 10 bodies, proteins such as fibrinogen, factor XII, and albumin adsorb to the surface. 11 Subsequently, factor XII becomes activated by surface contact triggering a series of 12 enzymatic reactions known as the intrinsic coagulation cascade [5]. One product of the 13 coagulation cascade is fibrin monomers, which associate to form protofibrils to form a 14 fibrotic mesh that traps platelet aggregates and other blood cells [6]. These two events 15 can jointly or independently lead to thrombus formation and influence the composition 16 of the thrombus. For example, in low shear stress environments such as found in blood 17 oxygenators and associated circuits, the level of platelet activation is low, and most 18 thrombi are fibrin-rich red clots [7]. On the other hand, the formation of platelet-rich 19 white thrombi has been linked to regions of high shear stress within PHV and 20 VAD [8,9]. In devices operating at low shear rates or with high surface area to volume 21 ratios, selective inhibition of the intrinsic branch of the coagulation cascade has been 22 able to dramatically reduce clot formation [10][11][12], suggesting that this fibrin-dominated 23 pathway is the primary driving force of clot formation in such devices. 24 Computational models have been used to study and predict the complex interactions 25 between hemodynamics and coagulation activation factors, platelets, and platelet 26 agonists [13]. Thrombosis models for medical devices have been derived from 27 mathematical models designed to predict in-vivo clot formation. For example, Taylor et 28 al. [14] adapted five existing models [15][16][17][18][19] to develop a macroscopic-continuum 29 thrombosis model. Another example, is the model of Wu et al. [20] who extended the 30 model of Sorensen et al. [18] and included platelet erosion effects first presented by 31 Goodman et al. [21]. Other models have been developed following the same principle of 32 adapting existing complex descriptions to perform simulations of thrombogenesis in 33 complex geometries [22][23][24][25]. 34 In most device-related thrombosis models, platelet activity plays a central and 35 oftentimes exclusive role. For example Blum et al. [22] and Wu et al. [26] developed 36 models to primarily predict platelet activation and its effects in rotary blood pumps. 37 One of the disadvantages of an overemphasis on platelet activity is that these models 38 cannot capture the clot formation initiated by factor XII adsorption which then 39 activates the intrinsic pathway of the coagulation cascade [12,27,28]. Therefore, this 40 deficit needs to be addressed to capture clotting across a wider range of medical devices, 41 including those devices in which platelet activity plays a minor role, such as blood 42 oxygenators, coils, venous catheters, vena cava filters, hemodialyzers, etc. [29]. 43 In contrast to device-related thrombosis modeling, a vast amount of work has been 44 performed to understand fibrin formation dynamics due to a blood vessel injury. The 45 modeling approaches range from kinetic pathway simulations of fibrin gel formation 46 with ordinary differential equations [30][31][32] to complex gelation models that include 47 blood flow interactions [33][34][35][36], as highlighted in a recent review by Nelson et al. [37]. 48 In the current work, the platelet based thrombosis model of Wu et al. [20] was 49 augmented with a partial model of the intrinsic coagulation cascade that includes fibrin 50 formation and contact activation. We calibrated the resulting fibrin-augmented model 51 by simulating a backward-facing-step flow channel (BFS) that has been extensively 52 characterized in-vitro [38,39]. Finally, to further test the predictive capabilities of the 53 model for a more clinically relevant geometry, we performed blood perfusion 54 experiments through a microfluidic chamber mimicking a hollow fiber membrane 55 oxygenator and validated the in-silico model against these experimental observations. 56 Materials and methods 57 In this section, the biochemical and fluid dynamic governing equations of the 58 fibrin-augmented thrombosis model are presented. In addition, two numerical setups are 59 introduced for the BFS and a hollow fiber bundle chamber chamber. The benchtop 60 experimental setup and blood perfusion experiments of the hollow fiber bundle are also 61 detailed. 62 Thrombosis Model 63 Following the approach of Wu et al. [20], the fibrin-augmented model combines the 64 incompressible Navier-Stokes (NS) equations with a system of 65 convection-diffusion-reaction (CDR) equations for the biochemical blood constituents 66 involved in thrombus formation. As the clot grows, a hindrance term is included in the 67 NS equations to impede blood flow in the regions where clot is present. 68

69
The pressure and velocity fields p and u are obtained by solving the equations of 70 conservation of mass and linear momentum: where µ is the asymptotic dynamic viscosity of blood and ρ is the density. The 72 scalar field ϕ represents the volume fraction of thrombus which comprises deposited 73 platelets and fibrin. (See Table 2.) To avoid a singularity in the source term, the 74 denominator is set as (1 − ϕ) ∈ [ϵ, 1], where ϵ is a small number. C 2 is the hindrance 75 constant introduced by Wu et al. [20], which assumes that the thrombus is composed of 76 densely compact spherical particles (2.78 µm). f (ϕ) is the hindrance function: that approximates interaction forces between thrombus and blood flow, and its 78 functional form is taken from Batchelor [40].

79
The scalar shear rateγ which is used throughout the thrombosis model is computed 80 via OpenFoam source code as [41]: Biochemical components

82
The transport of the biochemical constituents in time and space is modeled using a 83 system of CDR equations: where: c i is the concentration of species i, D i is the corresponding diffusion 85 coefficient, u is the velocity vector field, and r i is the reaction source term that 86 accounts for biochemical interactions.  [20]. k rpdB and k apdB are the rates at which resting and activated platelets deposit to the surface. Resting platelets can be activated by mechanical shear or the combination of agonists: ADP, TxA 2 and thrombin. Once deposited, AP d can be detached by flow shearing forces. B) Platelet deposition on deposited platelets including cleaning by large shear stresses.

Species i Definition
adenosine diphosphate λ j (k apa RP +k spa RP +k apa RP d +k spa RP d +θk rpd RP) 2.57 × 10 −10 TxA thromboxane A2 s pj (AP+AP d ) −k T xA2 TxA 2.14 × 10 −10 Table 1. Biochemical species involved in platelet activity. reaction rates with the constants k rpdB and k apdB , respectively. Resting platelets 99 can be activated chemically by TxA 2 , ADP, IIa, and mechanically by flow-induced 100 high shear stress. Once deposited, platelets can detach from the surface due to 101 shear cleaning.

102
(B) The artificial surface is covered by an initial layer of platelets followed by 103 additional deposition of platelets mediated by the k ra and k aa constant rates. The 104 source terms representing platelet activation, deposition, and cleaning dynamics 105 are listed in Table 1. The constant parameters involved in platelet activity and 106 hindrance terms are listed in Table 2.

107
Coagulation cascade and fibrin/fibrinogen deposition 108 Fig 2. A) Schematic depiction of fibrin production due to thrombin (Tb) cleavage and deposition of fibrinogen and fibrin at k F gB and k F nB kinetic rates, respectively. Deposited fibrin and fibrinogen (F g d , F n d ) detach from the boundary or thrombus due to the flow shear stress. B) Fibrin deposition to existing thrombus.
Fibrin formation and deposition is depicted in Fig. 2 that depicts the following 109 processes.

110
(A) Figure 2A shows fibrinogen and fibrin depositing to a surface according to k F gdB 111 and k F ndB following the Vroman effect. Fibrin production is mediated by IIa 112 which can be synthesized from activated platelets or via coagulation 113 reactions [18,46]. Coagulation reactions are modeled using the partial model of Thrombin generation rate on the surface of AP 3.69 kapa Biochemical platelet activation rate Ω Agonist activation weight na j=1 wj a j a jcrit [20] kspa Mechanical platelet activation rate   XIIa Activated factor XII k 1c XII (Only applied as reactive boundary surfaces) 5   [49] were simulated using OpenFoam [41]. prior to use, the blood was recalcified. Magnetic resonance imaging (MRI) was used to 136 quantify thrombus formation inside the recirculating region of the BFS in real-time.

137
The evolution of thrombus height and length were used as benchmarks to calibrate our 138 model.  Table 2. pressure value of zero was applied at the outlet, and zero gradient BCs were used at the 150 inlet and walls. Figure 4 shows the steady-state velocity field for the case with human 151 blood.

152
Variable Value/expression Table 5. Model parameters used in the bovine BFS simulation.
For the biochemical species, the physiologic concentration values listed in Table 6 153 were applied at the inlet. The initial condition was set to 0.1 % of the inlet 154 concentration for each species. The walls were prescribed to be reactive by allowing 155 factor XII to be activated due to contact as described in Méndez Rojano et al. [50].

156
Both platelets and fibrinogen deposited to the walls as described previously. The 157 specific model parameters used in the study are listed in Table 6; these values were also 158 used for the hollow fiber bundle simulation, described in the next section. For the 159 bovine simulation, the parameter values are listed in Table 5 160

Species
Value AP Human  A dual-time-step strategy was used in the thrombosis simulation. A large time step 165 was used to advance the biochemical species and thrombus formation (dt = 0.01 s), 166 while the flow equations were solved with a smaller time step to ensure numerical 167 stability (dt = 0.1 µs, CF L << 1). This strategy has been previously used in [51,52] 168 showing minimal impact on the simulation results. To improve the stability of the 169 thrombosis simulation, a steady-state flow solution was used as the initial condition.

Grid
Elements Relative Error (%) Coarse 70,000 NA Medium 158,400 0.21 Fine 280,000 0.001 Table 8. Grids used in the mesh convergence study for the BFS case. The relative error is based on the predicted recirculation length.    Table. 9. 215 Reactive boundary conditions were applied at the walls of the microfluidic chamber and 216 outer surfaces of the hollow fibers following Table 7. As mentioned previously, platelets, 217 fibrinogen, and fibrin can deposit on reactive boundaries. In addition, factor XII 218 activation takes place only on reactive boundaries.

Microfluidic hollow fiber bundle blood oxygenator chamber
219 Figure 6B shows a portion of the mesh that consisted of 2,756,704 elements with a 220 boundary layer thickness resolution of 5 µm. A mesh sensitivity analysis was performed 221 using the pressure drop across the microfluidic chamber (prior to thrombus growth) as 222 the comparison metric. The mesh sizes used in this study are listed in  Table 9. Inlet concentration for biochemical species for the microfluidic chamber case. The platelet count was taken from Lai [49] Grid Elements ( Table 10. Mesh sizes used in mesh sensitivity analysis for the microfluidic chamber case. The relative error was computed using the pressure drop across the chamber. To demonstrate the impact of fibrin formation, coagulation cascade, and contact 248 activation in the BFS configuration, the above results were compared to a simulation 249 with the original platelet-based model of Wu et al. [20]. Figure 9 shows the platelet Simulated and experimental thrombus formation patterns in the microfluidic hollow fiber bundle chamber at 15 min. Thrombus height was used as a surrogate for thrombus density to compare against experimental deposition maps computed from multiple µCT averaged scans from Lai to create a clot probability map [49] .

Fig 12.
Fibrin concentration field on the middle plane at 15 min. Insets: local fibrin structure around individual hollow fibers compared against experimental images from Lai [49] .

281
Thrombosis in low shear conditions is characterized by a predominance of fibrin, which 282 is the result of factor XII production stimulating the coagulation cascade [13,53]. where coagulopathies suggest increased thrombotic complications [59][60][61][62][63][64]. An additional 325 future application of the current fibrin model could be to simulate fibrinolysis strategies. 326 As highlighted by He et al. [65], another important challenge in blood-contacting 327 medical devices is the controlled formation of neointimal tissue while preventing 328 abnormal fibrovascular tissue known as pannus [9]. The model can be used to provide 329 the fibrin bed field required for neointimal tissue simulations.

330
Several limitations are acknowledged in the current work. Although our model was 331 able to reproduce certain time courses of thrombus growth, it predicted a slower growth 332 rates and different deposition dynamics compared to experiments. This could be 333 explained by the lack of accounting for red blood cells (RBC) in the model of 334 thrombosis, which in reality are trapped within fibrin clots and significantly increases 335 its volume and decreases its flow permeability. Another disparity was the prediction of 336 thrombus formation at the flow reattachment point of the BFS, which, although is in 337 line with the model of Taylor et al. [14], was not observed experimentally. A possible 338 explanation for this is that in the experiment, contrary to the simulation, the 339 reattachment region is not spatially fixed, but varies its location enough to promote 340 shear-induced washout of agonists. The list of parameters currently included in the 341 model is not necessarily universal or independent, and hence should be more carefully 342 and extensively studied coupled with experimental validation [52]. We have assumed 343 that phospholipids and calcium are in excess which translates to an overestimation of 344 thrombin and fibrin production. Additionally, the activation role of polyphosphate in 345 the intrinsic pathway was not considered [66]. The influence of RBCs in thrombus 346 composition, near-wall platelet margination, and non-Newtonian viscosity effects were 347 neglected as their inclusion are anticipated to be computationally prohibitive. Another 348 omission from the current model is any consideration of sudden embolization of clots 349 which may overestimate the final thrombus volume. Approaches for simulating 350 embolization such as the one presented by Zheng et al. [67] could be considered in the 351 future. Fibrin polymerization and subsequent fibrotic structure formation were not 352 explicitly modeled and should be considered if sparse fiber structures were to be 353 simulated in more detail. Finally, fibrin mechanics, deformation, and contraction and its 354 interaction with platelets and other cells [68] were not considered as this would render 355 the simulations excessively computationally intensive.

357
This report introduced an enhanced thrombosis model that incorporates a portion of 358 the coagulation cascade, specifically, the contact activation system and fibrin formation. 359 The importance of these mechanisms was illustrated in the simulation of thrombus