Computational Fluid Dynamics Simulations at Micro-Scale Stenosis for Micro ﬂ uidic Thrombosis Model Characterization

Platelet aggregation plays a central role in pathological thrombosis, preventing healthy physiological blood ﬂ ow within the circulatory system. For decades, it was believed that platelet aggregation was primarily driven by solu-ble agonists such as thrombin, adenosine diphosphate and thromboxane A2. However, recent experimental ﬁ nd-ings have unveiled an intriguing but complementary biomechanical mechanism — the shear rate gradients generated from ﬂ ow disturbance occurring at sites of blood vessel narrowing, otherwise known as stenosis, may rapidly trigger platelet recruitment and subsequent aggregation. In our Nature Materials 2019 paper [1], we employed micro ﬂ uidic devices which incorporated micro-scale stenoses to elucidate the molecular insights underlying the prothrombotic effect of blood ﬂ ow disturbance. Nevertheless, the rheological mechanisms associated with this stenotic micro ﬂ uidic device are poorly characterized. To this end, we developed a computational ﬂ uid dynamics (CFD) simulation approach to systematically analyze the hemodynamic in ﬂ uence of bulk ﬂ ow mechanics and ﬂ ow medium. Grid sensitivity studies were performed to ensure accurate and reliable results. Interestingly, the peak shear rate was signi ﬁ cantly reduced with the device thickness, suggesting that fabrication of micro ﬂ uidic devices should retain thicknesses greater than 50 µm to avoid unexpected hemodynamic aberra-tion, despite thicker devices raising the cost of materials and processing time of photolithography. Overall, as many groups in the ﬁ eld have designed micro ﬂ uidic devices to recapitulate the effect of shear rate gradients and investigate platelet aggregation, our numerical simulation study serves as a guideline for rigorous design and fabrication of micro ﬂ uidic thrombosis models. blood ﬂ ow rendering was colored by shear rate at constant viscosity. Note the shear rate maxima occurs at the stenosis apex and near the wall while the center forms a low shear pocket. (B and F) Peak shear rate γ max linearly correlates with input bulk shear rate γ 0 for both eccentric and concentric stenoses respectively. Note the concentric stenosis geometry demonstrates better linearity. The shear rate γ (C and G) and shear rate gradient γ ’ (D and H) distribution were analyzed along a sample streamline 1µm above stenosis apex spanning the shear acceleration ( x = − 100 – 0 µm) and deceleration ( x = 0 – 100 µm) zones. Note the γ and γ ’ present equal-space increment in regard to γ 0 . The γ max is located at x = 0 µm while the γ ’ max is located at x = − 5 µm


Introduction
Cardiovascular diseases have become the world's No.1 killer, as they are responsible for approximately a quarter of modern mortality worldwide [2,3]. Leading cardiovascular diseases such as coronary artery disease, carotid artery disease and deep vein thrombosis, which cause heart attack, stroke and embolism respectively, share a common root of being attributed to platelet aggregation [3]. These diseases together place a significant burden on the healthcare system and individuals, thus thrombosis is a clinically relevant research area of high demand, due to its prevalence and impact in society today.
Although it is well known that platelet adhesion, activation and aggregation play a central role in thrombosis, the precise interplay of mechanisms surrounding biochemical and biomechanical factors are yet to be fully delineated [4][5][6]. Our 2019 Nature Materials publication discovered that an integrin α IIb β 3 intermediate affinity state mediates biomechanical platelet aggregation and subsequent thrombus formation [1], and is induced by platelet mechanosensing pathways [7]. Biomechanical or rheological factors have long been recognized as crucial to thrombus formation [5,6], thus identifying, characterizing, and performing experiments on the entire range of physiologically relevant hemodynamic conditions is desired for further advancing our understanding of flow and shear-dependent thrombosis, and development of better antithrombotic drug therapies and biomaterials [3,8].

Stenosis-Induced Hemodynamic Effects of Flow Disturbance
Hemodynamic force largely exists with a wide range of physiological (Fig. 1A, left) and pathological (Fig. 1A, middle and right) shear conditions present in circulating blood flow, especially at a stenosis (Tab. 1). Stenosis describes a vessel or channel narrowing, which leads to flow velocity increase. This in turn also increases the shear rate γ (s −1 )-the tangential strain rate acting on vessel walls and blood components due to fluid flow. Clinically, there are two types of vessel stenoses typically observed: firstly, 'concentric stenosis' (Fig. 1A, middle) which usually occurs due to atherosclerotic plaque formation [9] and symmetrically narrows the vessel from two directions; and secondly 'eccentric stenosis' (Fig. 1A, right) which usually occurs due to medical device intervention [10,11] or existing thrombi formation [12], and asymmetrically narrows the vessel from one direction.
In the context of blood rheology, recent studies have highlighted the intriguing association between enhanced platelet adhesion, activation, aggregation and blood flow disturbance, specifically the shear rate gradient γ' as implicated in elongational flow, eddy currents, recirculation and reciprocation [5]. It is also hypothesized that the primed platelet adhesion is not simply a localized response to elevated shear and extensional force at the stenosis, but may be a result of the cumulative effects of the entire 'shear history' and concomitant cellular (particle) interactions experienced by platelets and red blood cells as they enter the stenosis contraction, pass through the stenosis apex and exit the stenosis expansion [18,19]. Thus, a CFD approach that can map 'shear distribution and rheological history' would be instrumental in investigating biomechanical thrombosis using microfluidics [6].

Application of Microfluidics to Thrombosis Researches
Microfluidics, or 'lab on a chip' systems, are popular tools for micro-scale investigation of hematological processes, by miniaturized control of agonist stimulation and ligand presentation [1,[20][21][22]. Moreover, taking advantage of their fluidic nature, microfluidic devices are increasingly invented and employed to examine the hemodynamic force effect on vascular and hematological cell behaviors [23][24][25], particularly the unique mechanosensitive features of platelet activation and thrombus formation [6]. In this context, microfluidic channels may be designed and fabricated with specific materials to recapitulate specific cardiovascular or medical device shear microenvironments. For example, a microcontraction or 'hump' may be designed to mimic concentric stenosis and flow disturbance [1,26] (Fig. 1B). Furthermore, while in vivo stenosis may have a complex 3D geometry, the hemodynamic shear effects of such 3D geometries may be simply mimicked using a precisely constructed 2D geometry, to mirror the same hemodynamic environment in vitro in an easily fabricated manner. The challenge in these microfluidic studies is to accurately recapitulate physiologically relevant flow conditions, which is mostly determined by flow medium selection.

Computational Fluid Dynamic Modelling for Interpreting Microfluidic Shear Experiments
While previous experimental studies have used CFD to demonstrate stenosis induced shear gradients [18,26], none of them have conducted whole channel numerical simulation nor systematically modeled the effect of different bulk flow mechanisms and fluid media on blood rheology. Furthermore, the employed verification schemes are usually unclear or void. Hereby, this CFD study comprehensively mapped the shear rate distribution of various flow conditions and examined the geometrical influence on microfluidic hemodynamics. Our efforts of numerical verification established CFD characterization principles for microfluidic models of thrombosis in the field.

Numerical Formulation and Governing Equations
In the present study, a finite volume approach using ANSYS FLUENT version 2020 R1 is utilized to solve the governing equations of fluid flow dynamics and obtain the flow field. It is assumed that the flow is incompressible, laminar and steady. Under these assumptions, the well-known Navier-Stokes equation can be expressed as where u; P; q; l denote the velocity vector, the pressure of the fluid, the density and the dynamic viscosity, respectively. In ANSYS FLUENT, the standard second order scheme as well as the second order upwind scheme were applied to discretize the pressure and momentum equations, with the Coupled algorithm utilized for pressure-velocity coupling.
To assess the effect of various physical properties of fluid on the shear rate distribution, two working fluids were examined; water and blood. The physical properties of the investigated working fluids are listed in Tab. 2.

Geometrical Properties and Boundary Conditions
A schematic diagram of the microfluidic device is illustrated in Fig. 1C. The axial length (X 0 = 200 µm), the width (Y 0 = 100 µm) and the height of the microfluidic device's cross-section (Z 0 = 130 µm) were kept as default geometric constants for all simulations.
A zero-gauge pressure boundary condition was applied at the inlet of microchannel (Fig. 1B). The noslip boundary condition was applied to the walls. Finally, a boundary condition of the experimentally implemented mass flow rate Q (kg/s) was implemented at the outlet of the flow region (Figs. 1B and 1C) as ANSYS inputs. The mass flow rate calculation is defined as: where c 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 device's cross-section [27]. A, D h and are calculated by the following equations:

Control Variables Selection
To investigate the flow rate influences on shear rate distribution, two studies were conducted with variables defined as per Tab. 3. Study 1) The bulk (wall) shear rate was varied from γ 0 = 150 to 3,000 s −1 by adjusting outlet mass flow rate from 32.5 to 649.6 µg s −1 as calculated from Eq. (3); Study 2) We assessed water and blood scenarios by adjusting the density and viscosity of fluid medium.
The shear rate and velocity profile were exported with a default streamline sampled as the single platelet trajectory passing 1 µm (half of a typical platelet diameter) above the stenosis apex, and z = 30 µm from the bottom (Fig. 1B, bottom).

Grid Sensitivity and Convergence Analysis
Grid sensitivity study is the key calibration step, which should be conducted before finalizing the results. Hereby, we analyzed blood shear rate and velocity profiles along a 20 μm vertical sample line, which is located 30 µm away from the channel wall in the z-direction and connects between the apex of the stenosis and the ceiling of the microfluidic channel ( Fig. 2A). To assure the independence of meshing size, three grids with 1,005,040 (Fine), 2,010,840 (Finer) and 2,868,187 (Finest) structured hexahedral meshed elements were examined (Fig. 1C) and six simulations for the eccentric stenosis geometry were performed for the lowest γ 0 = 50 s −1 (Figs. 2B and 2E) and the highest γ 0 = 3,000 s −1 (Figs. 2C and 2F). Moreover, during the iterative finite volume analysis, the computation is considered completed when the residuals (the root mean square error of approximation between two successive iterations) of the continuity equation are less than a certain value, defined as the "convergence criteria". To assure the results are also independent of the applied convergence criteria, three simulations were performed with the finer grid and γ 0 = 1,000 s −1 with residuals less than 1e −3 , 1e −6 and 1e −9 (Fig. 2F). We found the velocity and shear rate profiles are superimposed, demonstrating that the results acquired from our simulations are independent of the selected grid and convergence criteria. It is worth noting the following simulations are conservatively performed utilizing the finer grid and 1e −6 convergence criteria. Bulk shear rate γ 0 (s −1 ) γ 0 = 150, 600, 1,000, 1,500, 2,000, 2,500, 3,000 s −1 Fluid medium property Water; Blood; γ 0 = 50-1,050 s −1

Flow Rate Effect on Shear Rate Distribution
Over the past decades, extensive studies using parallel plate flow chambers and capillary slides have demonstrated shear-dependent platelet adhesive behaviors at various shear rates of laminar flow [25,28,29]. Nevertheless, the effect of shear rates in disturbed flow remains elusive. Thus, we first examined the effect of various bulk shear rates γ 0 = 150-3,000 s −1 on the disturbed flow region of eccentric (Fig. 3A) and concentric (Fig. 3E) stenoses. As observed in both types of stenoses, the peak shear rate γ max occurs at the middle of stenosis area while the peak shear rate gradient γ' max occurs at the upstream corner of stenosis geometry (Figs. 3D and 3H; 5 µm from the x-direction midpoint). Blood flow experienced linear increases in γ max and γ' max with respect to γ 0 (Figs. 3B and 3F). The γ max for eccentric and concentric stenoses increased 19.4 and 20.0 folds respectively as γ 0 increases from 150 to 3,000 s −1 (Figs. 3C and 3G). Similarly, the γ' max for both geometries increased 18.8 and 20.5 folds respectively (Figs. 3D and 3E). A close look at the eccentric (Figs. 3A-3D) and concentric (Figs. 3E-3H) stenoses with their γ and γ' distribution in all cases shows an equal-space increasing trend. Similar to eccentric stenosis, concentric stenosis geometry presents identical shear rate distributions in most of the cases. Notably, we found concentric stenosis has a more organized streamline which leads to smaller peak shear rate than eccentric stenosis (compare γ max = 24,300 s −1 in Fig. 3A vs. γ max = 26,900 s −1 in Fig. 3E). Furthermore, the concentric stenosis γ 0 -γ max linearity is stronger than that of eccentric stenosis (compare Figs. 3F vs. 3B). In addition, a more extensive high γ region at the ceiling is observed in the eccentric stenosis geometry (Figs. 3A and 3E), which could lead to higher probabilities of platelet aggregation. These results suggest that increased hemodynamic force is not purely a function of the stenosis level, but that the stenosis topology also plays a role. As a result, further investigation into the effects of geometrical shape and location of atherosclerotic plaque, thrombi, or medical device intervention, may be of clinical implication to controlling shear-dependent platelet activation.

Fluid Medium Effect on Shear Rate Distribution
The effects of the selected fluid medium on microfluidic experiment models are incompletely defined. For example, it is fundamental to understand whether the flow-induced shear force experienced by purified platelets in water based buffer is the same as that in whole blood [25,29]. It is commonly known that water is a Newtonian fluid, while whole blood displays Newtonian features at γ > 1,000 s −1 [30]. To examine the viscosity influence on shear rate distribution, we mapped the shear rate distribution of the entire microfluidic device. In these studies, water and blood were examined at a wide range of bulk shear rates (γ 0 = 50 to 1,050 s −1 ). The results presented in Fig. 4A illustrate negligible differences in γ max for the two fluid medium. As a result, water can be used to qualitatively map shear rate distribution as a pilot study, and thereby reduce the computational power and time required for microfluidic characterizations in our future practice.

Relationship between the Device Thickness and Shear Rate Distribution
Previous experimental observations have shown that for the same stenosis level, small variations in the height of the microfluidic channel may compromise the reproducibility of thrombosis observation [1]. With respect to the microfabrication principle, it is difficult to ascertain and design the optimal z-axis thickness with the dilemma of compromising between smooth blood flow with a large Z 0 > 100 µm, while saving time and material costs of microfluidic fabrication photolithography with a small Z 0 < 50 µm [31]. Hereby, we constructed three eccentric stenoses microchannels of thicknesses Z 0 = 10, 50, and 130 µm The 3D colormap of shear rate distribution in the eccentric stenosis and concentric stenosis microfluidic device for γ 0 = 1,000 (top) and 2,000 (bottom) s −1 respectively. The interstitial blood flow rendering was colored by shear rate at constant viscosity. Note the shear rate maxima occurs at the stenosis apex and near the wall while the center forms a low shear pocket. (B and F) Peak shear rate γ max linearly correlates with input bulk shear rate γ 0 for both eccentric and concentric stenoses respectively. Note the concentric stenosis geometry demonstrates better linearity. The shear rate γ (C and G) and shear rate gradient γ' (D and H) distribution were analyzed along a sample streamline 1µm above stenosis apex spanning the shear acceleration (x = −100-0 µm) and deceleration (x = 0-100 µm) zones. Note the γ and γ' present equal-space increment in regard to γ 0 . The γ max is located at x = 0 µm while the γ' max is located at x = −5 µm then performed CFD simulation at γ 0 = 1,000 s −1 (Fig. 4B). Interestingly, γ max was significantly reduced from 28,000 s −1 at Z 0 = 130 µm down to 4,180 s −1 at Z 0 = 10 µm (Fig. 4C). This result suggests that fabrication of microfluidic devices should retain Z 0 > 50 µm to avoid unexpected hemodynamic changes, which causes significant shear rate decreases between wall and medium. In addition, 'cutting-corners' by reducing the thickness of the microfluidic device (and thus cutting down on the raw materials required to fabricate the device, as well as the blood volume required per experiment), may compromise the accuracy of rheological results, leading to questionable biological conclusions in experiments.

Conclusion
Our systematic CFD study provides rheological characterization of microfluidics with stenosis under a 2D flow assumption, with respect to bulk shear rates and flow medium effects. The numerical results provide rheological insights that influence biomechanical platelet aggregation [1,26]. Overall, our results present 1) Both the peak shear rate and shear rate gradient linearly correlate with the bulk flow rate as the boundary condition; 2) To simulate microfluidic perfusion with isolated blood components, for example wih washed platelets, the fluid medium can be assumed as water so as to qualitatively compute the shear rate distribution; 3) Reduction of the microfluidic channel thickness to Z 0 < 50 µm incurs a caveat that compromises the accuracy of rheological modeling in the disturbed flow regions. Figure 4: Numerical calculation and modeling of flow medium and device thickness. A) Peak shear rate of water (blue) and blood (red) with respect to bulk shear rate. Note the blue and red lines are superimposed. B) The geometry modeling (left) and meshing results (right) of 10, 50 and 130 µm thickness devices. C) The shear rate distribution along a sample streamline 1 µm above stenosis apex spanning the shear acceleration (x = −100-0 µm) and deceleration (x = 0-100 µm) zones in the mid-plane of z-axis (z = 5, 25, 65 µm in height respectively) supervised the study, designed the research and wrote the paper. Research activities related to this work complied with relevant ethical regulations at the Faculty of Engineering, The University of Sydney.
Funding Statement: The author(s) received no specific funding for this study.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.