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

Yunduo Charles Zhao,Parham Vatankhah,Tiffany Goh,Jiaqiu Wang,Xuanyi Valeria Chen,Moein Navvab Kashani,Keke Zheng,Zhiyong Li and Lining Arnold Ju,＊

1School of Biomedical Engineering, Faculty of Engineering, The University of Sydney, Darlington, NSW 2008, Australia

2Charles Perkins Centre,The University of Sydney, Camperdown, NSW 2006, Australia

3Heart Research Institute, Newtown,NSW 2042, Australia

4School of Mechanical, Medical and Process Engineering, Queensland University of Technology, Brisbane, 4000, Australia

5Future Industries Institute, University of South Australia, Mawson Lakes, SA 5095, Australia

6South Australian Node,Australian National Fabrication Facility,University of South Australia,Mawson Lakes,SA 5095,Australia

7School of Aerospace, Mechanical and Mechatronic Engineering, Faculty of Engineering, The University of Sydney, Darlington,NSW 2008, Australia

ABSTRACT Platelet aggregation plays a central role in pathological thrombosis, preventing healthy physiological blood flow within the circulatory system.For decades,it was believed that platelet aggregation was primarily driven by soluble agonists such as thrombin,adenosine diphosphate and thromboxane A2.However,recent experimental findings have unveiled an intriguing but complementary biomechanical mechanism—the shear rate gradients generated from flow 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 microfluidic devices which incorporated micro-scale stenoses to elucidate the molecular insights underlying the prothrombotic effect of blood flow disturbance.Nevertheless, the rheological mechanisms associated with this stenotic microfluidic device are poorly characterized.To this end,we developed a computational fluid dynamics (CFD) simulation approach to systematically analyze the hemodynamic influence of bulk flow mechanics and flow medium.Grid sensitivity studies were performed to ensure accurate and reliable results.Interestingly, the peak shear rate was significantly reduced with the device thickness, suggesting that fabrication of microfluidic devices should retain thicknesses greater than 50 µm to avoid unexpected hemodynamic aberration, despite thicker devices raising the cost of materials and processing time of photolithography.Overall, as many groups in the field have designed microfluidic 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 microfluidic thrombosis models.

KEYWORDS Computational fluid dynamics; thrombosis; stenosis; hemodynamics; microfluidics

1 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-6].Our 2019Nature Materialspublication discovered that an integrin αIIbβ3intermediate 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].

1.1 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.

Table 1:Physiological and pathological shear rate ranges of in vivo circulating blood

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].

1.2 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-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-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, whilein vivostenosis 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 environmentin vitroin 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.

Figure 1:Schematics of vessel stenosis and the mimicking microfluidic models of thrombosis.(A)Blood flow in a healthy vessel without stenosis(left),in a vessel of concentric stenosis due to atherosclerotic plaque(middle)and a vessel of eccentric stenosis due to medical device intervention(right).(B)Eccentric stenosis mimicking microfluidic model.Note that platelet aggregation occurs at the downstream of the stenosis(top).Representative streamline of a platelet trajectory flowing through the microfluidic device(bottom).(C) 3D geometry and coordinate system of the eccentric stenosis model.Note the origin is located at the center of the front bottom line.(D) Computational mesh details for grid sensitivity analysis.Note the entire structure was meshed into 2,868,187 hexahetral elements for CFD analysis

1.3 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.

2 Method

2.1 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

whereu，P， ρ， μ 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.

Table 2:Working fluids modeled and physical properties

2.2 Geometrical Properties and Boundary Conditions

A schematic diagram of the microfluidic device is illustrated in Fig.1C.The axial length(X0=200µm),the width(Y0=100µm)and the height of the microfluidic device’s cross-section(Z0=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 rateQ(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 γ0(s−1)is the bulk shear rate,A(m2)is the cross-sectional area,Dh(m)is the hydraulic diameter,and,λ is the shape factor of the microfluidic device’s cross-section [27].A,Dhand λ are calculated by the following equations:

2.3 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−1by adjusting outlet mass flow rate from 32.5 to 649.6 µg s−1as calculated from Eq.(3);Study 2)We assessed water and blood scenarios by adjusting the density and viscosity of fluid medium.

Table 3:Control variable selection

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,andz=30µm from the bottom(Fig.1B,bottom).

3 Results and Discussion

3.1 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 thez-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−1with residuals less than 1e−3, 1e−6and 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−6convergence criteria.

Figure 2:Grid sensitivity and convergence analysis result under blood condition.(A)The vertical sample line of grid sensitivity study,which locates 30µm height to the bottom in the z-direction and connects 20µm spacing between the ceiling and the surface of stenosis apex along y-direction.(B-E)Shear rate and velocity profiles along the sample line under γ0 = 50 s−1 (B and D) and γ0 = 3,000 s−1 (C and E) with 3 different meshing sizes.Note the 3 lines are indistinguishable with each other.(F) The velocity distribution along the sample line under γ0= 1,000 s−1 with 3 convergence criteria of ＜1e−3, ＜1e−6 and ＜1e−9

3.2 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−1on the disturbed flow region of eccentric (Fig.3A) and concentric (Fig.3E) stenoses.As observed in both types of stenoses, the peak shear rate γmaxoccurs at the middle of stenosis area while the peak shear rate gradient γ’maxoccurs at the upstream corner of stenosis geometry (Figs.3D and 3H; 5 µm from thex-direction midpoint).Blood flow experienced linear increases in γmaxand γ’maxwith respect to γ0(Figs.3B and 3F).The γmaxfor eccentric and concentric stenoses increased 19.4 and 20.0 folds respectively as γ0increases from 150 to 3,000 s−1(Figs.3C and 3G).Similarly, the γ’maxfor 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−1in Fig.3Avs.γmax= 26,900 s−1in Fig.3E).Furthermore, the concentric stenosis γ0-γmaxlinearity is stronger than that of eccentric stenosis (compare Figs.3Fvs.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.

Figure 3:CFD simulated shear rates and shear gradients in eccentric(A-D)and concentric(E-H)stenoses geometry under various bulk shear rates of blood condition.(A and E) 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

3.3 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 γmaxfor 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.

3.4 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 optimalz-axis thickness with the dilemma of compromising between smooth blood flow with a largeZ0＞100 µm, while saving time and material costs of microfluidic fabrication photolithography with a smallZ0＜50 µm [31].Hereby, we constructed three eccentric stenoses microchannels of thicknessesZ0= 10, 50, and 130 µm then performed CFD simulation at γ0=1,000 s−1(Fig.4B).Interestingly,γmaxwas significantly reduced from 28,000 s−1at Z0=130µm down to 4,180 s−1at Z0=10µm(Fig.4C).This result suggests that fabrication of microfluidic devices should retainZ0＞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.

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)

4 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 toZ0＜50 µm incurs a caveat that compromises the accuracy of rheological modeling in the disturbed flow regions.

Acknowledgement:We thank Qing Li, Chi Wu, Lawrence Hao Huang and Sarah Keogh for their helpful discussion, and Sydney Manufacturing Hub, Gregg Suaning,Simon Ringer, Kevin Doan and Nadia Court for support of our lab startup.This work was supported by Australian Research Council Discovery Project (DP200101970-L.A.J.and Z.L.), NSW Cardiovascular Capacity Building Program (Early-Mid Career Researcher Grant-L.A.J.) and Sydney Research Accelerator prize (SOAR-L.A.J.), The University of Sydney Faculty of Engineering Startup Fund and Major Equipment Scheme (L.A.J.).Lining Arnold Ju is an Australian Research Council DECRA fellow (DE190100609), an honorary National Heart Foundation Future Leader fellow (102532), Faculty of Engineering Dean’s Research Award 2020 recipient.This work was conducted (in part) using the ‘Research and prototype foundry’ facilities at the New South Whales node and the ‘Design House’ facilities at the South Australia node of the Australian National Fabrication Facility (ANFF), a company established under the National Collaborative Research Infrastructure Strategy to provide nano-and micro-fabrication facilities for Australia’s researchers.

Author Contributions:Y.C.Z., P.V.and T.G.performed the research, analyzed the data and co-wrote the paper.L.Z., J.W., M.K.and K.Z.provided critical advice.X.V.C.helped revise the manuscript.L.A.J.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.