Hemodynamic analysis for stenosis microfluidic model of thrombosis with refined computational fluid dynamics simulation

Disturbed blood flow has been increasingly recognized for its critical role in platelet aggregation and thrombosis. Microfluidics with hump shaped contractions have been developed to mimic microvascular stenosis and recapitulate the prothrombotic effect of flow disturbance. However the physical determinants of microfluidic hemodynamics are not completely defined. Here, we report a refined computational fluid dynamics (CFD) simulation approach to map the shear rate (γ) and wall shear stress (τ) distribution in the stenotic region at high accuracy. Using ultra-fine meshing with sensitivity verification, our CFD results show that the stenosis level (S) is dominant over the bulk shear rate (γ0) and contraction angle (α) in determining γ and τ distribution at stenosis. In contrast, α plays a significant role in governing the shear rate gradient (γ′) distribution while it exhibits subtle effects on the peak γ. To investigate the viscosity effect, we employ a Generalized Power-Law model to simulate blood flow as a non-Newtonian fluid, showing negligible difference in the γ distribution when compared with Newtonian simulation with water medium. Together, our refined CFD method represents a comprehensive approach to examine microfluidic hemodynamics in three dimensions and guide microfabrication designs. Combining this with hematological experiments promises to advance understandings of the rheological effect in thrombosis and platelet mechanobiology.

www.nature.com/scientificreports/ concomitant cellular interactions 15,16 . Thus, it is important to map the hemodynamic distribution at both the pre-and post-stenosis regions 5 . In a healthy patent blood vessel (Fig. 1A, top), the blood flow pattern is laminar and steady, yielding a constant velocity profile, γ and τ distribution. However, when a vessel is narrowed concentrically from all directions (e.g. atherosclerotic plaque formation; Fig. 1A, middle) 30 , or eccentrically narrowed from one direction (e.g. medical device insertion; Fig. 1A, bottom) 31,32 , the τ and γ 0 become elevated at the stenosis with blood flow acceleration (Table 1). To mimic the pathological microvascular conditions, stenosis microfluidic models of thrombosis have been designed with a default width of Y 0 = 100 µm, height of Z 0 = 130 µm, α = 85° and S = 80% (Fig. 1B) 13,33 . It is worth noting that the flow profile within such microchannels is considered two-dimensional (2D), as the geometry and flow vary only in the x-y plane at the stenotic region ( Fig. 1B) 7,13,34 . Previous studies have justified these stenosis microfluidic models of thrombosis which have been well characterized and validated experimentally 6,7,13,16,38,40 , demonstrating the shear rate gradient effects on platelet aggregation in blood flow disturbance. The advantages of such microfluidic approaches enable direct visualization of biomechanical platelet aggregation at the downstream face of stenosis (Fig. 1C). Nevertheless, the detailed rheological mechanisms, or specifically the exact τ, γ and γ ′ thresholds that trigger such biomechanical platelet aggregation are incompletely understood. In this context, CFD simulation represents the first step in revealing the hemodynamic profile within disturbed blood flow and correlating with experimental results of thrombotic response.
CFD is the most popular computational method to model hemodynamic parameters i.e. τ, γ and γ ′ and simulate microfluidic outcomes before experiments are done 7,15 . However, rarely do the CFD practices nowadays thoroughly show the mesh sensitivity verification, or systematically benchmark the effects of the microfluidic boundary conditions on the hemodynamic outcome with certain control variables i.e. γ 0 , S, α and fluid medium Blood flow in a healthy vessel without stenosis (top), in a vessel with concentric stenosis due to atherosclerotic plaque (middle) and a vessel with eccentric stenosis due to medical device insertion (bottom). (B) CFD contour maps in the eccentric (top) and concentric (bottom) stenosis microfluidic channels. Note that the channel walls were colored to display the WSS distribution; the representative streamline of a platelet trajectory was colored to display the shear rate γ distribution. (C) Differential interference contrast microscopy image of biomechanical platelet aggregation in an eccentric stenosis microfluidic channel after whole blood perfusion at γ 0 = 1,800 s −1 , mimicking the physiological shear rate in arteries and arterioles 7 . To visualize platelet aggregates in better clarity, the microchannel was washed with Tyrode's buffer after whole blood perfusion. Note that platelets aggregate at the downstream face of the stenosis. Scale bar = 20 µm. (D) ANSYS finite volume meshing scheme for the concentric stenosis model as in panel C. The entire microfluidic channel was meshed into 1,638,336 hexahedral elements for CFD analysis. The coordinate origin is located at the center of the front bottom edge. Note that for illustration purpose, a coarse mesh is shown in the figure. www.nature.com/scientificreports/ ( Table 2). The caveat is that the CFD solutions may have some extents of error with coarse mesh sizes. Besides, although the non-Newtonian property of blood has significant impacts on the viscosity, its influence in the disturbed flow region remains poorly characterized. Hereby, we present an ultra-fine CFD study to map hemodynamic profiles for stenosis microfluidic models of thrombosis, and address the above concerns with mesh sensitivity verification and analytical validation.

Methods
Stenosis microfluidic models of thrombosis. Microfluidics is a miniaturized approach that has significant advantages in controling hemodynamic parameters, ligand presentation and agonist stimulation at the micro-scale 13,[35][36][37] . Microfluidic channels with hump-like contractions have been developed and utilized to mimic microvessel stenosis and examine the prothrombotic effect of flow disturbance through which whole blood or washed platelets are perfused ( Fig. 1C) 6,7,13,38-41 . Numerical formulation and governing equations. The commercially available software ANSYS FLU-ENT version 2020 R1 is utilized to computationally simulate the flow profile. The flow was assumed as steady, laminar, and incompressible. Under these assumptions, the fluid can be described using the continuity Eq. (1) and Cauchy momentum Eq. (2) as follows: where u is the 3D velocity vector (m s −1 ), P is the pressure (Pa), µ is the dynamic viscosity (Pa s) and ρ is the density (kg m −3 ). This can be simplified to the well-known Navier-Stokes equation for the case of Newtonian fluids (fluids with constant dynamic viscosity): In ANSYS FLUENT, to assess the viscosity effect we simulated 3 working mediums: water, blood as a Newtonian fluid, and blood as a non-Newtonian fluid under the Generalized Power-Law (GPL) viscosity model ( Table 2).
Numerous non-Newtonian models have been used throughout literature for simulating blood flow: Carreau 42 , Modified Cross Law (Carreau-Yasuda) 43 , Power-Law (PL) 44 , non-Newtonian Power-Law 44 , Generalized Power-Law 45 , Casson 46 and Walburn-Schneck Law 47 are among the well-recognized models 48 . In this study, the GPL model is utilized due to its accuracy in calculating the shear rate γ 48 , which is calculated from: where K is the consistency index, and n is the Power-Law index. The two parameters are calculated from 45

Geometrical properties, boundary conditions and Reynolds number.
A schematic diagram of the microfluidic channel is illustrated in Fig. 1B. By default, the axial length is X 0 = 200 µm, the width is Y 0 = 100 µm and the height of the cross-section is Z 0 = 130 µm. A zero-gauge pressure boundary condition is applied at the inlet of the microchannel (Fig. 1B). A no-slip boundary condition is applied to the walls. Finally, a boundary condition of the experimentally implemented flow rate Q (µL min −1 ) is implemented at the outlet of the flow region ( Fig. 1B,C). The flow rate calculation is defined as 49 : www.nature.com/scientificreports/ where γ 0 (s −1 ) is the bulk shear rate, A (µm 2 ) is the cross-sectional area, D h (m) is the hydraulic diameter, and, λ is the shape factor of the microfluidic cross-section 49 . A, D h and are calculated by the following equations: The above also allows for the calculation of the Reynolds number, a dimensionless group used to indicate whether fluid flow in a channel is characterized as laminar or turbulent. The Reynolds number for a Newtonian fluid is given by: Substituting Q from Eq. (7) into Eq. (11) yields the following expression of the Reynolds number for this particular investigation: Using γ 0 = 150 and 3,000 s −1 for Eq. (12), the Reynolds number for this study is found to range from 0.04 to 0.9. Unlike the turbulent state of blood flow in the stenotic region of the common carotid arteries with large diameters shown in some studies [50][51][52] , a Reynolds number with such low magnitude indicates a laminar flow regime, as flow typically does not transition into a turbulent regime before the critical point of Re = 2,300 53 , validating our laminar flow assumption.
Furthermore, our steady flow assumption stems from the Windkessel effect which dampens the pulsatile flow changes far away from the heart, leading to steady flow in the capillary region 54,55 . This steady assumption is further validated by the Womersley number which is a dimensionless group to evaluate the pulsatile flow effect in relation to oscillation frequency and viscosity 56 .
where D is the diameter of the vessel (m); f and ω are the oscillation and angular frequency (both Hz) respectively. We selected f = 2 Hz and D = 130 µm, corresponding to a heart rate of 120 beats per minute to validate our assumption at even extreme conditions. The calculated Womersley number of 0.128 is much less than the critical point of Wo = 1 56 , validating our steady flow assumption.
By default, the shear rate and velocity profile are simulated with a streamline sampled as the single platelet trajectory passing 1 µm (half of the diameter of a single platelet) above the stenosis apex, and z = 30 µm from the bottom (Fig. 1B).
In addition, we set up an ultra-fine CFD mesh with approximately 2,000,000 elements ( Fig. 1D and Supplementary Fig. S1) for all the studies to achieve high accuracy. www.nature.com/scientificreports/ Note that the peak shear rate γ max occurs at the stenosis apex. (C and G) The peak shear rate γ max and peak WSS τ max linearly correlated with the input bulk shear rate γ 0 for both eccentric and concentric stenosis microfluidics respectively. (D and H) The shear rate and shear rate gradient history of a single platelet particle trajectory 30 μm from the bottom wall and 1 μm from the symmetric geometry surface at γ 0 = 2000 s −1 for eccentric and concentric stenosis microfluidic channels. Note that the γ max ≈50,000 s −1 and γ ' max ≈5,000 µm −1 s −1 for eccentric and concentric stenosis respectively.  A and E) The streamlines of blood flow in the eccentric and concentric stenosis microfluidic channels for S = 70% (top) and 90% (bottom) respectively. The representative streamline of a platelet trajectory was colored by shear rate values. (B and F) Peak shear rate γ max exponentially correlates with the input stenosis level S for both eccentric and concentric stenoses respectively. Note that the concentric γ max is higher than that of eccentric stenosis channel. The shear rate γ (C and G) and shear rate gradient γ ′ (D and H) distribution is plotted along a sample streamline 1 µm above the stenosis apex spanning the shear acceleration (x = −100 to 0 µm) and deceleration (x = 0-100 µm) zones. Note that γ max and γ ' max occur at the same location as in Fig. 2

Results
Bulk shear rate impact on WSS and shear rate distribution at stenosis. It has been well recognized that platelets activate and aggregate in response to flow disturbance [57][58][59] . In microcirculation, blood flow is considered laminar and steady as the change in geometry does not transfer the flow into a turbulent regime. Here we investigated a wide range of γ 0 = 150-3,000 s −1 and the effects on the τ and γ distribution using our refined CFD method. Although τ max is similar on both eccentric ( Fig. 2A) and concentric (Fig. 2E) stenosis microfluidic wall with 7% difference, a higher τ distribution is observed at the ceiling of the eccentric stenosis (compare Fig. 2A vs. E), which could lead to higher probabilities of vessel damage 48 .
Then we compared a sample platelet trajectory streamline across the stenotic region (Fig. 2B,F), displaying 10% higher γ max in the eccentric stenosis than the concentric stenosis. Interestingly, as γ 0 increases from 150 to 3,000 s −1 , the τ max and γ max both increased linearly with respect to γ 0 about 20.0-and 19.4-fold respectively for eccentric stenosis (Fig. 2C). Similarly, the τ max and γ max increased by 16.7-and 20.0-fold respectively for Note that the eccentric stenosis upstream shear rate is higher than that of the concentric stenosis channel. The shear rate γ (C and G) and shear rate gradient γ ′ (D and H) distribution is plotted along a sample streamline 1 µm above stenosis apex spanning the shear acceleration (x = −100 to 0 µm) and deceleration (x = 0-100 µm) zones. Note that the shear rates in the post-stenosis zone remain unchanged with respect to α, while the upstream of the pre-stenosis area displays faster acceleration. www.nature.com/scientificreports/ concentric stenosis (Fig. 2G). Notably, we found eccentric stenosis has larger γ max and γ ' max than concentric stenosis (compare γ max = 53,875 s −1 , γ ' max = 25,701 µm −1 s −1 in Fig. 2D vs. γ max = 48,691 s −1 , γ ' max = 23,295 µm −1 s −1 in Fig. 2H). Of note, our simulated τ max displays strong correlation to γ max . Thus, we focused on the shear rate γ as the examining parameter in the following studies.
Relation between micro-contraction geometry and hemodynamic profile at stenosis. Flow disturbance caused by vessel lumen constriction has a major prothrombotic effect via platelet mechanosensing 7 . Hereby, we examined γ and γ ′ distributions as functions of S for degree (Fig. 3) and α for topology (Fig. 4) of vessel narrowing. Notably, we observed that S plays a significant role in the shear rate acceleration towards γ max (Fig. 3A,E). The γ max for eccentric and concentric geometries increased 34.6-and 41.0-fold respectively as S increased from 30 to 95% (Fig. 3C,G). Similarly, the γ ' max for both geometries increased 28.3-and 35.7-fold respectively (Fig. 3D,E). Notably, we identified an exponential γ max increment once the stenosis level exceeded 70% (Fig. 3B,F). Trajectory analyses demonstrated similar increases in γ and γ ′ in the S > 70% regimes for both concentric (Fig. 3C,D) and eccentric (Fig. 3G,H) stenoses.
In contrast, α did not significantly affect the hemodynamics within the stenotic region, and subsequently had little effect on the γ and γ ′ distribution (Fig. 4B,F). It is worth noting that with increasing α, the upstream face or pre-stenosis area displayed rapid acceleration in γ (Fig. 4C,G) and γ ′ (Fig. 4D,H), whilst the post-stenosis downstream area was not affected. Our findings demonstrated that the shear rate profile changes most significantly when S is varied, suggesting that (1) although the development of vessel narrowing might be negligible at earlier stages of development, its hemodynamic impacts are rapidly enhanced once severe pro-occlusive conditions are reached; (2) the impact of the stenosis contraction angle is negligible due to its limited influence on γ max , however, its significant impact on the shear acceleration zone should be further investigated, especially its elongational effects acting on blood cells and plasma proteins 4,9 . Hemodynamic benchmarks of viscosity models. The selection criteria for viscosity models on microfluidic simulation remain incompletely standardized. For example, it is fundamental to understand whether the hemodynamic forces experienced by washed platelets in Tyrode's buffer is the same as those in whole blood 41,60 . Accurate computational modelling of blood is a complicated task due to multiplexed components in whole blood. RBCs consist a large amount, 45% of whole blood by volume or hematocrit 61 . Considering that RBCs have viscoelastic properties, blood exhibits shear-thinning behavior (decreasing viscosity with increasing shear rates), or non-Newtonian feature, at low shear rate γ < 1000 s −1 (Fig. 5A) 62 . While at high shear rates γ > 1000 s −1 , blood behaves as a Newtonian fluid with a linear relationship between γ and τ, in another words nearly constant viscosity (Fig. 5A) 62,63 .
To benchmark the viscosity model choice, we mapped the viscosity changes with respect to γ 0 by applying five different constitutive viscosity models: one Newtonian and four non-Newtonian 62 , where a wide range of bulk shear rates γ 0 = 50-1050 s −1 were examined. All the non-Newtonian models have similar viscosity predictions at low γ 0 , while the viscosity of the PL model becomes inaccurate at high γ 0 region. Besides, we also found negligible difference in γ max for water, blood as a Newtonian fluid, and blood as a non-Newtonian fluid (Fig. 5B). Interestingly, qualitative similarities between water and blood suggest that blood can be considered Newtonian in CFD analyses of shear rate distribution, even in low γ 0 conditions (Fig. 5B).
This finding supports the feasibility of CFD practices with the water Newtonian model for microfluidic characterization as a reductionist approach with much reduced computational costs. Having said that, the water Newtonian model does not fully recapitulate the behavior of whole blood 61 . Viscoelastic models, such as the simplified Phan-Thien-Tanner or Gisekus 61 , are good alternatives to capture the blood rheology in future studies.

Discussion
Our numerical studies present a refined CFD approach to map hemodynamic parameters for both concentric and eccentric stenosis microfluidic models of thrombosis, providing unprecedented rheological insights underlying biomechanical platelet aggregation and thrombosis. The results present that (1) the stenosis level S is the major determinant of the shear rate γ and shear rate gradient γ ′ within disturbed flow at the micro-contraction; (2) the contraction angle α plays a significant role in governing γ ′ , while having negligible influence on γ max ; (3) the shear rate γ experienced by washed platelets in Tyrode's buffer is similar to that in whole blood; (4) water as a Newtonian fluid can be applied to microfluidic CFD characterization as a reductionist approach with reduced computational costs.
Previous experimental studies have demonstrated that in vivo aggregation is sensitive to γ ′ at stenosis 7 . The follow-up studies suggested that γ 0 , S and α affect platelet aggregation in terms of kinetics, stability and sizes of aggregates 38 . The underlying platelet mechanobiology at molecular and cellular levels is exciting yet incompletely understood 5 . One possibility is that a γ ′ threshold exists for VWF elongation, and its subsequent conformational activation 6,13,60,64 . The other working mechanism is due to the collision and compression between RBCs and platelets at the stenosis and surface of developing aggregates 16,41 . To investigate these mechanobiological mechanisms, future experiments are required to correlate our CFD results with molecular and cellular behaviors of biomechanical platelet aggregation in these stenotic microfluidic devices 13 .

Data availability
The data that supports the findings of this study are available from the corresponding author upon reasonable request. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.