Microfluidic Channel for Simulating Wall Shear Stress Waveforms in Carotid Bifurcation Blood Flow

1. Introduction

The carotid sinus, located at the commencement of the internal carotid artery, is characterized by a slightly thinner vascular wall and an expanded lumen. This anatomical feature induces blood flow separation near the vessel wall and vortex formation, resulting in a region of low oscillatory Wall Shear Stress (WSS). Extensive research [1,2,3,4,5,6,7,8,9] indicates that such low oscillatory WSS can cause arterial endothelial damage and dysfunction, promoting platelet aggregation and the release of growth factors. These processes facilitate the migration of smooth muscle cells from the arterial wall’s medial layer to the intima, where they proliferate. This migration and proliferation lead to the accumulation of cells, contributing to the formation of atherosclerotic plaques and ultimately leading to carotid artery stenosis. When stenosis progresses to occlusion, reduced blood flow can cause cerebral ischemia, leading to prolonged brain hypoperfusion and ischemia, which in turn may result in brain atrophy, cognitive decline [10], neurological dysfunction [11], and reduced cerebrovascular reactivity [12]. Hence, low oscillatory WSS-induced endothelial dysfunction in the carotid sinus is a critical initiating factor in these pathologies.

Understanding the pivotal role of low oscillatory WSS in atherosclerosis development and progression, as well as investigating the relationship between low oscillatory WSS and endothelial dysfunction, including its mechanobiological mechanisms, is at the forefront of cardiovascular biomechanics research [13,14,15]. While animal models and human clinical trials offer direct research approaches, they are influenced by various physiological factors such as respiration and neural activity. Moreover, these methods face challenges like significant individual variability, high costs, extended observation periods, and ethical concerns. In vitro models that simulate endothelial WSS waveforms offer an alternative, free from these physiological influences, and allow for a focused quantitative exploration of the relationship between low oscillatory WSS and endothelial dysfunction and their underlying mechanobiological mechanisms [16].

Parallel-plate flow chamber systems were among the first methodologies for simulating the fluid shear stress microenvironment in vitro to construct endothelial cell culture models. Chiu et al. [17] pioneered the use of a parallel-plate flow chamber with a stepped structure for examining endothelial cells’ mechanobiological responses under oscillatory WSS. Subsequently, researchers have employed such chambers with varied structures and dimensions to investigate the impact of low oscillatory WSS on aspects like endothelial cell signal transduction, gene expression, and structural and functional changes [18,19], as well as the differentiation of endothelial progenitor cells [20]. More recently, Wang et al. [21,22]developed an endothelial cell culture model that simulates the WSS waveform of oscillatory flow in the common carotid artery post-exercise, exploring the mechanobiological effects and mechanisms triggered by exercise-induced oscillatory WSS.

Despite their utility, these systems for simulating oscillatory WSS typically suffer from limitations such as large size and high manufacturing costs. With the advancement of micro- and nano-manufacturing and microfluidic technologies, researchers have devised various microfluidic chip systems for simulating oscillatory flow WSS, enabling in vitro studies of endothelial cell mechanobiology [23,24,25,26]. These microfluidic systems offer advantages over parallel-plate flow chambers in terms of reduced size and cost [27,28,29,30]. However, most of these in vitro models have only qualitatively and coarsely represented the characteristics of oscillatory WSS waveforms, failing to precisely depict the oscillatory flow WSS features at arterial bifurcations and bends. A particular challenge has been the simultaneous representation of both the low oscillatory WSS waveform in the carotid sinus and the high pulsatile flow WSS waveform in the downstream uniform arterial segments within a single endothelial cell culture model.

This study aims to design a microfluidic chip system with a variable cross-section and a stepped structure to concurrently reproduce both the low oscillatory WSS in the carotid sinus and the high pulsatile flow WSS waveforms in the downstream uniform arterial segments. The approach involves designing a stepped, variable cross-section microfluidic channel; simplifying the three-dimensional (3D) fluid dynamics problem into a two-dimensional (2D) Computational Fluid Dynamics (CFD) model, considering the channel’s significantly smaller height compared to its length and width; and finally, using 2D CFD simulations to quantitatively analyze the impact of the microfluidic channel’s step height and width on the oscillatory flow WSS waveforms, leading to the optimization of the channel’s structure and dimensions.

2. Materials and Methods

2.1. Design of the Variable Cross-section Microfluidic Channel

To effectively replicate both the low oscillatory WSS at the carotid sinus and the high mean pulsatile flow WSS waveform in the downstream uniform arterial segments, a variable cross-section microfluidic channel structure was developed based on principles of fluid dynamics and microfluidic chip technology. As illustrated in Figure 1(a), this channel primarily comprises a stepped variable cross-section flow chamber and circular inlets and outlets with radius R. The stepped flow chamber is divided into two sections: (I) a narrower rectangular area at the front and (II) a wider rectangular area at the rear. The stepped design in the height direction of the narrower section (I) generates vortices near its rear, simulating the low oscillatory flow WSS waveform at the carotid bifurcation’s sinus [31,32]. The wider rear section (II) exhibits laminar flow, suitable for modeling the laminar WSS in the uniform part of the carotid bifurcation. By adjusting the height of the step in region (I) relative to the width of region (II), a precise simulation of both low oscillatory WSS and high pulsatile flow WSS waveforms in the downstream uniform arterial segment is achieved.

2.2. Governing Equations of Flow in the Microchannel

Assuming that the fluid within the microfluidic channel is an incompressible Newtonian fluid, its flow adheres to the Navier-Stokes equation

$$\frac{\partial \overrightarrow{V}}{\partial t}+(\overrightarrow{V}\cdot \nabla )\overrightarrow{V}=-\frac{1}{\rho }\nabla p+\frac{\mu }{\rho }{\nabla }^2\overrightarrow{V}$$

(1)

and the continuity equation

$$\nabla \overrightarrow{V}=0$$

(2)

where $\overrightarrow{V}$represents the fluid velocity, μ denotes the fluid viscosity, and ρ signifies the liquid density.

Assuming a uniformly distributed velocity profile at the inlet, the boundary conditions are established as

$$\overrightarrow{V_i}=-\frac{Q(t)}{\pi R^2}\overrightarrow{j}$$

(3a)

$$\overrightarrow{V_w}=0$$

(3b)

$$p_0=0$$

(3c)

Where $\overrightarrow{V_i}$ is the inlet velocity, Q(t) represents the volumetric flow rate of fluid entering the inlet, $\overrightarrow{j}$ is the unit vector along the y-axis in the Cartesian coordinate system, $\overrightarrow{V_w}$ is the velocity at the channel wall, and p₀ is the pressure at the outlet.

Obtaining analytical solutions for equations (1) and (2) is typically challenging, so three-dimensional CFD numerical simulation is generally employed. Given that the height (h₁ and h₂) of the microchannel structure in this study is considerably smaller than its width (b₁ and b₂) and length (L₁ and L₂), the influence of side boundaries on fluid flow can be neglected in the microfluidic channel far from the inlet and outlet. Therefore, the 3D flow model is simplified to a 2D model on the oxy plane (as shown in Figure 1(c)). In the rectangular coordinate system oxy, equations (1) and (2) are simplified to

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

(4)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2})$$

(5)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\frac{\mu }{\rho }(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2})$$

(6)

where u and v are the flow speeds in the x-direction and y-direction, respectively.

The corresponding boundary conditions are adapted as

$$u_i=-\frac{Q(t)}{h_1{b}_1}$$

(7a)

$$u_w=v_w=0$$

(7b)

$$p_0=0$$

(7c)

Despite the simplification, obtaining analytical solutions for equations (4) to (6) remains complex. Consequently, this study employs 2D CFD numerical simulation to resolve the flow field.

Upon determining the flow field velocity distribution through numerical simulations, the formula for calculating the WSS at the bottom of the microfluidic channel is given by

$${\tau }_{w1}=\mu \frac{\partial u}{\partial y}{|}_{y=-h_1,0}$$

(8)

Specifically, in the wider uniform rectangular channel section (II), disregarding the effects of the front variable cross-section section, side, and outlet boundaries, the flow is treated as a 2D problem between infinitely large microchannels, further simplifying equations (5) to (7) to

$$\frac{\partial u}{\partial x}=0$$

(9)

$$\frac{\partial u}{\partial t}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial y^2}$$

(10)

$$0=-\frac{1}{\rho }\frac{\partial p}{\partial y}$$

(11)

Here, p represents the pressure, with the boundary conditions further simplified to

$$u{|}_{y=-h_2,0}=0$$

(12)

Under quasi-steady flow conditions, the formula for calculating the WSS at the bottom of the uniform section in the wider rectangular area (II) is [33]

$${\tau }_{w2}(t)=\frac{6\mu Q(t)}{b_2{h_2}^2}$$

(13)

In this formula, Q(t) is the volumetric flow rate entering the microfluidic channel, and b₂ and h₂ respectively represent the width and height of the wider rectangular channel in the latter part of the microchannel.

2.3. Carotid Bifurcation Blood Flow WSS Waveform Extraction

Figure 2. (a) 3D model of the carotid bifurcation with the stenosis at the internal carotid sinus; (b) Waveform of flow rate at the common carotid artery measured by ultrasound Doppler; (c) Low oscillatory WSS waveform at the site of carotid sinus stenosis; (d) High pulsatile WSS waveform in the downstream segment of the internal carotid artery.

Figure 2. (a) 3D model of the carotid bifurcation with the stenosis at the internal carotid sinus; (b) Waveform of flow rate at the common carotid artery measured by ultrasound Doppler; (c) Low oscillatory WSS waveform at the site of carotid sinus stenosis; (d) High pulsatile WSS waveform in the downstream segment of the internal carotid artery.

To capture WSS waveforms at various positions within the carotid bifurcation, numerical simulations of the carotid artery’s flow field were performed using COMSOL 6.0, a CFD simulation software. The 3D carotid artery model, reconstructed from actual human medical imaging, includes the common carotid artery, the internal carotid artery with an enlarged sinus, and the external carotid artery. A protruding plaque at the internal carotid sinus was integrated into the model to simulate in vivo stenosis conditions. The model’s total length is 52 mm, with diameters of approximately 6.4 mm for the common carotid artery, 2.6 mm for the internal carotid artery, and 2.2 mm for the external carotid artery. The common carotid artery’s inlet boundary condition was derived from in vivo ultrasound Doppler-measured axial blood flow velocities, processed through the Womersley formula [34] (see Figure 2(b)). The outlets of the internal and external carotid arteries were set to a zero-pressure boundary condition (default reference pressure is one standard atmosphere) [35,36]. Under low Reynolds number laminar flow and no-slip boundary conditions, the fluid velocity and WSS distribution at the carotid bifurcation were computationally determined. This study focuses on the low oscillatory flow WSS at the disturbed flow site of carotid sinus stenosis (Figure 2(c)) and the high pulsatile WSS in the downstream uniform arterial segment (Figure 2(d)) as the target WSS waveform features for simulation in the microfluidic channel.

2.4. CFD Simulation and Optimization of WSS in the Variable Cross-section Microchannel

Figure 3. Schematic of the optimization procedures for CFD simulation and structure optimization of the microfluidic channel.

Figure 3. Schematic of the optimization procedures for CFD simulation and structure optimization of the microfluidic channel.

To replicate the carotid artery’s in vivo WSS waveform characteristics within the microfluidic channel, fluid simulation software was utilized for WSS modeling and microchannel dimension optimization. As depicted in Figure 3, the process commenced with calculating the fluid volumetric flow rate Q(t) in the microchannel using Equation (13), guided by the high pulsatile shear stress waveform in the downstream uniform arterial segment (shown in Figure 2(c)), and under the predefined channel dimensions h₂ = 0.15mm and b₂ = 5mm. Subsequently, the velocity waveform at the inlet boundary was deduced from the microchannel inlet’s cross-sectional area. CFD simulation was then employed to obtain the flow field distribution and target WSS waveform at the step, followed by a similarity analysis with the intended low oscillatory WSS waveform. If the detected WSS waveform did not meet the desired similarity, the step dimensions b₁ and h₃ in the microchannel were adjusted, and the process iterated until the step’s WSS waveform achieved satisfactory similarity with the expected pattern. The optimized microchannel dimensions were finalized upon reaching this criterion.

In optimizing the step dimensions within the microchannel, the influence of step height h3 and width b1 on the WSS waveform was initially analyzed, setting the heights and widths to h₃ = 0.03 mm, 0.06mm, 0.09mm and b₁ = 2mm, 3mm, 4mm, respectively. The step was subsequently optimized to attain the desired fitting effect. During optimization, the step height varied in 0.03 mm increments, and the width in 1 mm increments. Furthermore, altering the width of the rear wider area to b₂ = 4mm, the impact of Q(t) variations on the step’s WSS waveform was examined under the same optimization protocol. The best-sized 3D model post-optimization was chosen to assess the WSS waveform consistency in the wider rear area and at the step. Using the coordinate system from Figure 1, five points at different x-positions were selected in both areas along the microchannel’s central axis (z = 0mm, y = -0.15mm) to examine the WSS waveform. Additionally, vortex changes at the moments of maximum (t_p = 0.11s) and average (t_m = 0.16s) flow velocities were observed at different cross-sectional areas. Lastly, to more directly observe how the low oscillatory WSS waveform varied with changes in the cross-sectional area of region (I), the waveform was extracted at a distance of l = 0.025mm from the step.

The CFD simulation process began by importing the microchannel model, designed in SolidWorks 2021, into COMSOL software for mesh generation. The inlet velocity, determined based on the high WSS waveform anticipated in the downstream microchannel, was set as the inlet condition. The channel walls were defined as no-slip surfaces, and the outlet was set to zero pressure. The study focused on a flow chamber with a narrow section width of b₁ = 4mm and step height of h₃ = 0.09 mm for mesh independence verification. It was found that with over 60201 mesh elements, the mesh size had a negligible impact on the results. Therefore, for a balance of accuracy and computational efficiency, CFD simulations were conducted with a mesh count exceeding 60201.

3. Results

Employing a 3D reconstructed carotid artery model from clinical imaging and CFD numerical simulations, we obtained the low oscillatory WSS waveform at the carotid sinus and the high pulsatile WSS waveform at the uniform downstream section of the carotid bifurcation. Utilizing the design of a variable cross-section microchannel and controlled inlet flow, we systematically analyzed the influences of channel width (b₁), step height (h₃), and input flow rate (Q(t)) on the WSS within the microchannel. Through optimization of the microchannel dimensions, we achieved a more precise simulation of both low oscillatory and high pulsatile WSS waveforms in different regions of the microchannel.

3.1. Carotid Artery WSS Waveform Simulation in the Variable Cross-section Microchannel

Figure 4. (a) WSS waveforms in the wider rectangular region of the microfluidic channel; (b) Low oscillatory WSS waveforms at the stepped site.

Figure 4. (a) WSS waveforms in the wider rectangular region of the microfluidic channel; (b) Low oscillatory WSS waveforms at the stepped site.

Based on the target WSS waveforms at the carotid bifurcation (Figure 2c-d), and following the CFD simulation and optimization process for the microchannel (Figure 3), we obtained simulation results for WSS in different regions of the optimized microchannel (Figure 4). In the wider rectangular section at the rear of the microchannel (h₂ = 0.15mm, b₂ = 5mm), the WSS waveform closely matched the targeted high pulsatile WSS waveform in the internal carotid artery. This consistency demonstrates that the microchannel structure designed can accurately replicate the high pulsatile WSS waveform in the uniform segment downstream of the carotid bifurcation. Along the fluid flow direction, WSS changes remained consistent, allowing for precise simulation of high pulsatile WSS in this laminar flow region. In the stepped region, after optimizing the channel dimensions (h₃ = 0.09mm, b₁ = 4mm), the produced low oscillatory WSS waveform at the step showed similarity to the WSS changes at the carotid sinus stenosis in terms of oscillation pattern and amplitude. The subtle variations in WSS waveform along the fluid flow direction were attributed to the vortices formed in the stepped area, where differences in velocity gradients at various positions cause WSS variations. Thus, the optimized variable cross-section microchannel successfully simulates both low oscillatory WSS and high pulsatile WSS waveforms of the carotid artery. Moreover, the WSS waveforms from both 3D and 2D simulations at the same location in the stepped area showed good agreement. To optimize computational resources and time, 2D numerical simulation was employed for flow field modeling and dimension optimization in the stepped area.

3.2. Impact of Step Size on Low Oscillatory WSS

Figure 5. Changes of streamline distributions with step height and width at the stepped site in the time interval of the maximum flow velocity (t_p = 0.11s) for b₂ = 5 mm.

Figure 5. Changes of streamline distributions with step height and width at the stepped site in the time interval of the maximum flow velocity (t_p = 0.11s) for b₂ = 5 mm.

In efforts to precisely simulate the low oscillatory WSS waveform, we analyzed the effects of step height (h₁) and width (b₁) on the low oscillatory WSS generated in the stepped area. Observing vortex changes at the moments of maximum (t_p = 0.11s) and average (t_m = 0.16s) flow velocity (Figure 5 and Figure 6), it was found that as the microchannel step width (b₁) increased, the vortex area on the oxy plane gradually diminished. At a height of h₃ = 0.03mm, increasing width b₁ even resulted in the disappearance of the vortex area, making low oscillatory WSS simulation challenging. The size of the vortex area varied consistently with step width b₁ across different cardiac cycle moments, with the vortex area size being directly proportional to the flow rate. As the flow rate decreased, the vortex in the stepped area vanished under microchannel dimensions of b₁ = 2mm, h₃ = 0.03mm. Correspondingly, within the vortex area, the negative oscillation amplitude of the WSS waveform decreased as step width b₁ increased (Figure 7). This decrease is attributed to reduced flow speed and subsequently lower velocity gradients when inlet flow rate is constant and width increases.

With a fixed microchannel step width (b₁), increasing step height (h₃) led to an enlargement of the vortex area on the oxy plane. The vortex area size correlated directly with the average flow rate at different cardiac cycle moments (Figure 5 and Figure 6). Correspondingly, within the vortex area, the negative oscillation amplitude of the WSS waveform increased as step height h₃ increased (Figure 7).

When keeping the width of region (I) constant and altering step height (h₃), the size and position of the vortex area changed. As the step height (h₃) decreased, the height and width of the vortex area reduced, with its center shifting towards the step. This shift caused the low oscillatory WSS waveform at the same position (Figure 7) to first increase and then decrease in negative oscillation amplitude as the vortex area’s center moved. When width b₁ decreased, the vortex area remained largely unchanged, but the negative oscillation amplitude of the low oscillatory WSS waveform increased.

3.3. Influence of Wide Area Channel Size on Low Oscillatory WSS

Figure 8. Changes of streamline distributions with step height and width at the stepped site in the time interval of the maximum flow velocity (t_p = 0.11s) for b₂ = 4 mm.

Figure 8. Changes of streamline distributions with step height and width at the stepped site in the time interval of the maximum flow velocity (t_p = 0.11s) for b₂ = 4 mm.

Figure 9. Changes of streamline distributions with step height and width at the stepped site in the time interval of the average flow velocity (t_m = 0.16s) for b₂ = 4 mm.

Figure 9. Changes of streamline distributions with step height and width at the stepped site in the time interval of the average flow velocity (t_m = 0.16s) for b₂ = 4 mm.

With a fixed target WSS for simulation in the wide area of the microchannel, the channel size in this area determines the internal flow rate (Equation (12)). As shown in Section 2.2, whether vortices form in the stepped area and how WSS changes under various flow rates within a cardiac cycle depend on the inlet flow rate. Therefore, we investigated the effect of flow rate variations on the vortices and WSS simulation at the step by altering the channel width in the wide area. The simulations (Figure 8 and Figure 9) indicated that changing the channel width in the wide area (b₂ = 4mm) resulted in low oscillatory WSS at the step, with variations in step width b₁ and height h₃ following a trend consistent with when b₂ = 5mm (Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10). Compared to b₂ = 5mm, the reduced flow rate in the microchannel at b₂ = 4mm diminished the vortex area, leading to the disappearance of vortices at b₁ = 2mm, h₃ = 0.03mm (Figure 8 and Figure 9). Simultaneously, the reduced flow rate led to the decrease in the negative oscillation amplitude of the WSS waveform (Figure 10).

4. Discussion

In this study, leveraging microfluidic technology and fluid dynamics principles, we developed a variable cross-section microfluidic cell culture chamber with a stepped design. This chamber effectively replicates both the low oscillatory WSS found at the carotid sinus and the high pulsatile WSS waveform characteristics in the downstream uniform arterial segment as observed in vivo. Simulation results showed that with a step height of 0.09 mm, a width of 4 mm, and a wider channel width of 5 mm in the variable cross-section microchannel, the stepped area efficiently simulates the low oscillatory WSS features associated with carotid sinus stenosis. The wider portion of the microchannel accurately models the high pulsatile WSS waveform in the uniform arterial section downstream of the carotid sinus. Additionally, the formation of vortices and variations in low oscillatory WSS within the stepped region of the microchannel correlate with the step height, width, and dimensions of the wider channel area. An increase in step width b₁ leads to a decrease in both the vortex area and the magnitude of negative WSS oscillations. When the step height h₃ falls below 0.03 mm, the vortex area vanishes, challenging the simulation of low oscillatory WSS patterns. The variable cross-section microfluidic channel designed in this research lays a foundational basis for simulating carotid artery WSS waveforms and for conducting in vitro studies on the interplay between low oscillatory WSS at the site of carotid sinus stenosis and arterial endothelial dysfunction.