Numerical simulation of the three‐dimensional dynamics of healthy and hardened red blood cells passing through a stenosed microvessel by immersed boundary‐lattice Boltzmann method

Red blood cells (RBCs) infected by a malaria‐like disease can change their membrane stiffness and deformability drastically. The objective of this article is to study the three‐dimensional dynamics of healthy and hardened red blood cells passing through a stenosed microchannel in a Poiseuille flow by combining the immersed boundary and lattice Boltzmann methods. We numerically simulate the transient motion of cells of different degrees of membrane stiffness and observe that hardened cells require a longer traverse time to squeeze through the constriction. We found that, compared with a healthy cell, a hardened cell shows much less change in its bending behavior. We also analyzed the projected deformation index and the velocity profiles of the cells passing through the constricted vessel. Different values of viscosity contrast are examined and show a minor effect on the passage time and bending energy change. In addition, a previous study of the dynamics of a group of three RBCs passing through a constricted microchannel is extended to three‐dimensional in the current work.

Vahidkhah and Fatouraee 2 study the dynamics of single and a pair of red blood cells in a stenosed arteriole by a combined 2D immersed boundary-lattice Boltzmann method (IB-LBM). They analyzed the sensitivity of elastic moduli, the cell-cell aggregation behavior and concluded that the interactions between fluid viscous forces, the shear and normal stresses in the RBC membrane, and the intercellular adhesion strength are the main factors that influence the motion of RBCs in the channel. Alizadeh and Dadvand 3 simulate in two-dimensional IB-LBM the motion of a single high deformable (hRBC) and a low deformable (iRBC) red blood cell in a stenosed microvessel. By measuring the movement and deformation of an RBC under different Reynolds numbers, they found that a healthy RBC travels faster than a sick one. Hoque et al, 4 utilizing a framework of finite-sized dissipative particle dynamics (FDPD), investigate the 2D deformation dynamics of an RBC flowing through a stenosed microchannel. They observed that, compared with a healthy cell, the deformation of a hardened cell was much less, leading to an increase in passage time. Using 2D IB-LBM, Wang et al 5 study the movement and deformation of a single healthy and an infected red blood cell in a Poiseuille flow through a constricted microchannel. They specifically examined the dynamics and transient behavior of a group of three healthy RBCs and a group of three hardened cells passing through a constricted microchannel. It was found that the RBCs encountered more significant forced deformation with declining constriction ratio as they squeezed through the region of constriction compared with other parts of the microchannel, and the total traverse time of the sick cells is longer than the healthy ones.
Shi et al 6 suggest that studies in two-dimensional are essential but may also be insufficient, as they are not possible to describe the three-dimensional out-of-plane deformation of the cells. By using the smoothed particle hydrodynamics method (SPH), Wu and Feng 7 simulate in 3D the traverse of a single healthy and a malaria-infected red blood cell through a converging microfluidic channel formed by a narrow cubic duct, and investigate the gradual loss of cell deformability. One of their conclusions is that the rigidity of the cell membrane impedes its passage across microfluidic channels so that the transit time increases with the growing membrane modulus. Gallage et al, 8 also utilizing SPH, study the three-dimensional deformation behavior of an RBC when it moves through a stenosed capillary. They analyzed the deformation index, mean velocity, and energy change of RBCs with different membrane stiffness. Soleimani et al, 9 using a combined shell-fluid analysis based on the SPH method, model the deformation of a single healthy and a hardened RBC passing through a stenosed capillary. They concluded that the rheological behavior of an RBC changes because of infection by a disease, and is reflected in its deformability while passing through microvessels. Hashemi and Rahnama 10 investigate the behavior of a single mildly stiffened RBC in a straight microcapillary flow by three-dimensional IB-LBM where only the shear elasticity modulus of the cell is considered. By varying the values of shear modulus, the initial orientation of the cell, and the gradient of flow pressure, they present the statics of deformation index and time to reach steady state, and the transient behavior of RBCs. Three-dimensional numerical simulations specifically for high hematocrit and dense suspension of RBCs in a constricted blood vessel can be found in References 11,12, and for those immersed with other particles in  This work presents a quantitative and qualitative study comparing the three-dimensional motion and deformation of healthy red blood cells (hRBC) and infected red blood cells (iRBC) passing through a stenosed microvessel in a Poiseuille flow by a combined immersed boundary-lattice Boltzmann method (IB-LBM). We consider five different forces/constraints of RBC plus a repulsion force between the wall. Motions of cells with different degrees of hardness are simulated. Statistics between the dynamics of healthy and hardened cells are obtained, including traverse time, changes in bending energy, mean velocity, and projected deformation index. The effects of viscosity contrast for aged cells are also separately studied. We then extend part of the work of Wang et al, 5 which studies the 2D dynamics of a group of three healthy RBCs and a group of hardened cells passing through a constricted microvessel, into a 3D scenario.

MODEL AND METHODS
In this section, we briefly summarize the mechanical model and the accompanied elastic forces of an RBC, the lattice Boltzmann method (LBM), the immersed boundary method (IB), and the numerical simulation environment.

RBC model
Blood is essentially a concentrated suspension of deformable RBCs, each of which is considered as an active 2D membrane, composed of an outer lipid bilayer and an underlying spectrin cytoskeleton, and immersed in a surrounding 3D fluid. An undeformed RBC has a biconcave discoid shape with geometry profile given by, 16 where C 0 = 0.207, C 1 = 2.002, C 2 = −1.122, and R 0 is the radius of an RBC at its relaxed state. A biconcave RBC discoid of diameter 7.82 μm and thickness of 2 μm formulated in Equation (1) is modeled by a spring-network with 1280 triangular faces and 642 vertices as shown in Figure 1 where each edge of the triangle represents an elastic spring connecting two adjacent nodes. 17 The elastic forces of an RBC are contributed by five components, 18 which are the stretching force F s , the bending force F b , the local area constraint F a , the global area constraint F A , and the volume constraint F V . The theoretical background can be found in References 17-21. Note that because of the geometric nature, it is common for two-dimensional simulation to include only the stretching and the bending forces.

Stretching force
The stretching force is responsible for preserving the relaxed length of the spring between two adjacent nodes in the triangulated mesh of the cell. Identical to the force of an elastic spring, it pulls the two end nodes back when the spring is elongated and pushes them apart when the spring is compressed. Such a stretching and compressing motion of the spring resembles the cytoskeleton's resistance to the in-plane shearing. For each spring edge, the stretching force acted on one end node can be written as where k s is the stretching modulus, is the stretching ratio, ( ) is the compensation for the nonlinearity of the spring, 22 ΔL is the spring length deviation from its relaxed state, and n 12 is the unit vector pointing from node X 1 to node X 2 .

Bending force
Provided by the lipid bilayer of the cell membrane, the bending force represents the resistance to bending. On each node X k it is written as where and k b is the bending modulus, is the triangle index, X k is the vertex of the tessellated membrane, T(X k ) is the set of triangles that share the same vertex X k , V( ) is the set of vertices of triangle , dA(X k ) is the effective area of X k , H(X p ) is the mean curvature at X p , and n is the unit outward normal vector of triangle . Figure 1 shows an example of triangular patches for the calculation.

Local area conservation
To keep the density of the lipid unchanged in the bilayer, 21 the membrane will shrink or expand to preserve its surface area.
In the model, the function of the local area force is to keep the area of each triangle in the meshed network nearly a constant. For each vertex on a triangle, the local area force is where k a is the local area constraint modulus, X k is the vertex, A 0 is the area of triangle at rest, ΔA is the area dilation of triangle , and w k is the unit vector pointing from the centroid of triangle to vertex X k . The total local area force for node X k is the contribution of all triangles that contain X k as a member.

Global area conservation
Since an RBC is enclosed by a membrane, the cell's volume and the surface stay nearly a constant. The global area force is introduced to conserve strongly the total surface area of the cell membrane. It can be regarded as a penalty for surface deviations while the cell deforms. For each vertex on triangle , the allocated global area force is where k A is the global surface constraint modulus, A 0 T is the total surface area at equilibrium, ΔA T is the total area deviation, and w k is the unit vector pointing from the centroid of triangle to vertex X k . The total global area force for node X k is the contribution of all triangles that contain X k as a member.

Volume conservation
Likewise, the volume force conserves the total volume of the cell strongly. For each vertex on triangle , the allocated volume force is where k V is the volume constraint modulus, V 0 is the total volume of the undeformed cell, ΔV is the volume deviation, A is the area of triangle which contains vertex X k , n is the unit outward normal vector of triangle . Similarly, the total volume force acting on node X k is the sum of contributions of all triangles that contain X k as a member.

Lattice Boltzmann method
In recent years, the lattice Boltzmann method (LBM) has become a developing alternative to computational fluid dynamics (CFD) in simulating hydrodynamic problems including incompressible fluid flows, porous media flows, multiphase and multicomponent flows, blood flows, and thermal flows. 23 Examples of its attractive features are simplicity, easy implementation, and parallel in nature. 17,19 When Mach number and Knudsen number are small enough, LBM can be shown to recover the macroscopic Navier-Stokes equations by the Chapman-Enskog analysis. 19,24 Based on the Bhatnagar-Gross-Krook (BGK) approximation, 25 the LBM equation with single relaxation time and external forces is written as: where f i (x, t) can be interpreted as the discrete distribution of the probability of having a particle at position x and time t, traveling with velocity c i , Δt is the time step, f eq i (x, t) is the equilibrium distribution function, is the single relaxation time determining the rate of approach to local equilibrium, and S i (x, t) is the sum of external forces.
The function of the equilibrium distribution is defined as follows, where , u are the fluid density and velocity, respectively, i is the lattice weight coefficients, and c s is the lattice speed of sound.
The external force f from an immersed object is incorporated into Equation (7) as, 26 D3Q19, the most commonly employed lattice model in 3D problems, is used in the present work and shown in Figure 2. The lattice weight coefficients i are defined as, The 19 lattice velocities along each direction are, Finally, the macroscopic fluid density and velocity are derived from the first two moments of the populations, where the second term Δt 2 f in Equation (12) is included to keep second-order space-time accuracy of the fluid solution. 26 In the meanwhile, the kinematic viscosity of the fluid is determined as a function of the relaxation time by,

Immersed boundary method
The immersed boundary method is a coupling formulation that obeys the no-slip condition at the interface of the fluid structure. 27 When the LBM grid distance is set to unity, the boundary forces F on the Lagrangian immersed object nodes X are spread into the Eulerian fluid lattice x by, where is the three-dimensional Dirac delta distribution factorized by the relation (x) = (x) (y) (z).
The new velocity U of cell nodes is interpolated back from the fluid velocity u by, While there are various discretized delta functions having different interpolation ranges, we adopt the so-called 4-point stencil as follows,

Computational setup
In the present numerical simulation, a physical space of 30 μm × 15.5 μm × 15.5 μm is set to 60 × 31 × 31 lattices. One lattice unit is equal to 0.5 μm which is selected after a grid-independence study, and one time step is equal to 4.17 × 10 −5 milliseconds. Dimensions hereafter are all in lattice units unless otherwise specified. For a healthy RBC, the parameters k s = 0.03, k b = 0.006, k a = 0.01, k A = 0.1, and k V = 0.3 have been used. The Reynolds number of the fluid is set to 1.0. The mean edge length of the meshed RBC surface is around 1.01 lattice unit so that the node-sticking and fluid-penetration problems of IB are avoided. 28 A constant body-force-driven Poiseuille flow is implemented along the x-direction. Periodic boundary conditions are applied at the channel's inlet and outlet; a full-way bounce back is applied for the vessel wall. The lattice fluid Mach number is verified to be far less than 0.1 during the simulation. Besides the five forces described in Section 2.1, a repulsion force F r 19 between the cell node and the wall node is also introduced and is activated only when the node-to-node distance falls below a cut-off length d cut , where k r is the repulsion modulus, d is the distance between the cell node X k and the in-contact wall node, d cut is set to 0.5 μm, and n r is the unit vector pointing from X k to the wall node. The simulation algorithm consists of the following steps: 1. Compute at the cell nodes the various elastic forces as described in Section 2.1 and the repulsive force if necessary. 2. Spread the total forces at each cell node to fluid lattices by immersed boundary method. 3. Perform the LBM algorithm (ie, compute propagation, equilibrium distribution, collision, and full-way bounce back). 4. Calculate the fluid velocity with force correction as in Equation (12). 5. Interpolate the fluid velocity to the cell nodes. 6. Update the cell position. 7. Go back to step 1 for the next time step.

RESULTS AND DISCUSSIONS
In this section, we first validate our model and the numerical method by the well-known tank-threading motion and the deformation index analysis. We then present the dynamics of healthy and infected red blood cells passing through a stenosed microvessel in a Poiseuille flow.

Validation studies
Tank treading motion in a shear flow Fischer 29 found that an RBC in a shear flow exhibits a steady motion in which the cell deforms into an ellipsoid-like particle at a stationary orientation with the membrane circulating about the cell interior similar to the motion of a tank tread. Kraus et al 28 show that the stationary state of an ellipsoidal vesicle in a linear shear flow is characterized by the inclination angle between its major axis and the plane of shear and by its revolution frequency of a tank treading motion. In addition, both the inclination angle and the revolution frequency depend strongly on the reduced volume of the vesicle. In a cubic computational domain of size 60 × 31 × 31, we induce a simple shear flow by moving the top and bottom walls with equal velocity w x = 0.08 but in opposite x-directions. To get different values of reduced volume, ellipsoidal vesicles are formed by varying the lengths of the three major axes. The stretching modulus is set to zero to reflect the nature of a vesicle. As shown in Figure 3, we obtained the steady-state profiles of the vesicle in good agreement with Kraus

Deformation index
In their experimental study, 30 Tsukada et al report that an RBC deforms into a steady parachute configuration in a pressure-driven microchannel flow. To quantify the degree of deformation, they defined a deformation index (DI = l/d) as shown in Figure 4 (left). In our numerical simulation, a healthy RBC positioned with an initial inclination of = ∕2 at the centerline of a Poiseuille flow is progressed until a steady state. Together with the experimental data of Tsukada, 30 the deformation index is plotted as a function of mean cell velocity in Figure 4 (right), which shows that our simulation result falls well in the range of the experiment.

Effect of the membrane stiffness
When an RBC is infected by a disease like malaria, it gets hardened, loses its deformability, and tends to block up the microvessel. In this section, we compare the deformation behavior and the passage time between a regular healthy RBC (hRBC) and hardened ones (iRBC) with different degrees of stiffness. Literature that studies the dynamics of a biconcave RBC passing through a microvessel used to match the initial position of the cell centroid to the centerline of a symmetric flow domain. As a result, the cell marches and deforms symmetrically, and the transient three-dimensional out-of-plane motion, which differentiates mostly from 2D computation, of the cell is not recorded. In this study, we intentionally position an RBC slightly off the centerline with a bias toward the lower vessel wall so that the 3D asymmetrical out-of-plane deformation of the cell can be easily observed. Figure 5 depicts the motion and deformation of a healthy RBC passing through a constricted vessel where the narrowest diameter is of 8.5 μm at x = 30 in a steady Poiseuille flow. The tapered stenosis is geometrically constructed by fully rotating a Gaussian distribution function perpendicular to the flow direction. When passing through the stenosis, the cell gradually changes its shape to a slipper, Fosbury-jumps over the stenosis, and recovers its biconcave shape. The total traverse time is 1.65 milliseconds (39 685 time steps), for the mass center of the cell starting at x = 10 and reaching x = 50. Figure 6A shows the motion and deformation of a healthy RBC passing through a constricted vessel where the narrowest diameter is of 5.5 μm at x = 30. The total traverse time is 15.03 milliseconds (3.61 × 10 5 time steps). Since a healthy RBC becomes hardener when infected by a disease, its rigidity can be reflected in the increase of the bending and shear moduli. We define the hardened ratio (iRatio) as the increase of the moduli k b and k s . For example, when iRatio equals 3, the bending modulus k b and the shear modulus are 3-fold of a healthy RBC. In Figure 6B-D, the motion and deformation of iRBCs with different degrees of hardness passing a stenosed vessel of a diameter of 5.5 μm is shown, where the iRatio is 3, 10, and 20, respectively. The total passage time is 27.57, 31.56, and 44.84 milliseconds, respectively, and the corresponding temporal variation is presented in Figure 8A. The variations of shapes and paths in each stage are clearly observable; each cell's end trace and orientation appear quite different. Figure 7 shows the total traverse time spent for eight cases of different iRatio. The R-square between the linear regression line and the simulation value is 0.97.

Deformation index
Deformation index is one of the measures to describe the shape changes of cells in transition. For examples, the Tsukuda index in Section 3 and the frequently seen Taylor index, both are suitable for automatic 2D image processing. For tracking asymmetrical 3D deformation, we define more appropriately the projected deformation index as, DI i = L i /D , where D is the diameter of a healthy biconcave RBC at rest, and L i is the projected length of RBC in the axial directions (i = x, y, z ). When an RBC is positioned initially perpendicular to the flow, its DI x is about 0.31, and both DI y and DI z are equal to one. As the cell approaches the narrow segment of the constricted channel, an abrupt change in all DI's can be noticed. When entering the stenosis, the cell gradually changes its posture from vertical to horizontal. DI x rises to 1.0 and keeps almost constant afterward. DI y and DI z are relatively steep at both sides of the narrowest point of the microvessel at x = 30, since the RBC turns from a biconcave discoid to a slipper marching in the flow direction. Projected deformation indices reveal qualitative information on how the cell changes its three-dimensional shape during the simulation. The significant difference in the evolution of DI y and DI z signals a three-dimensional out-of-plan motion that does not appear in a two-dimensional simulation.

Bending energy profile
When an RBC travels through the stenosis, the most significant change is the increase of its bending energy. 8 In the present work, the bending energy E b of an RBC is given by, 20 where H is the mean curvature vector, N v is the total number of vertices in the surface triangulation of the cell, and dA(X k ) is the effective local area at the vertex X k . Figure 9A plots changes in the bending energy of four RBCs with different bending moduli. The most diseased RBC, which is also the most rigid, displays the least change of bending energy. Figure 9B shows the evolution of the x-component of the mean velocity of RBCs in the flow. Inside the constricted zone, the healthy cell deforms into a more curled slipper-shape accepting larger pushing force from the flow and shows higher moving speed. After passing through the stenosis, it is pushed downward, resulting in a smaller x-component of its velocity. In general, the higher the stiffness, the slower the RBC moves inside the constricted segment.

Effects of the viscosity contrast
A diseased or aged RBC may also show an increase in its interior viscosity. We define the viscosity contrast as = interior cell viscosity / fluid viscosity. In Figure 10, The difference in traverse time caused by viscosity contrast is much lower than that of the hardened ratio in Figure 6. Note that for a healthy red blood cell, the actual viscosity contrast equals 5. In the previous sections, we set = 1 to keep the simulation simple.

Multiple cells passing through a constricted microchannel
By IB-LBM, Wang et al 5 study the two-dimensional movement and deformation of three healthy and infected red blood cells in a Poiseuille flow through a constricted microchannel. Here we extend part of Wang's study into the three-dimensional scenario.
The simulation domain is set as 200 × 50 × 50 lattice units equivalent to a physical space of 100 × 25 × 25 μm. Three healthy (or sick cells) are initially positioned vertically in the flow instead of being placed horizontally like those in Wang's. To avoid overlapping the animation images, we set the center cell to start six lattice units ahead of the other two cells. The total traverse time is determined when the mass center of the center cell reaches x = 110. The sick cells are hardened by 10-fold increase in the stretching and bending moduli than the healthy cell. An additional cell-to-cell repulsion force is introduced with the same functional form as Equation (17). Figure 11A plots the total traverse time vs the constriction ratio 1 − d/D where, as defined in Wang's, d is the narrowest diameter at the constriction center, and D is the outer channel diameter. More specifically, D is 48 lattice units, and d is tested at 24, 20, 16, 12, and 10 lattice units, respectively. Given an example, when d equals 12, the constriction ratio is 0.75. As before, one lattice unit is equal to 0.5 μm, which is chosen after a grid-independence study, and one time-step is equal to 4.17 × 10 −5 milliseconds. Figure 11B shows the movement of the mass center of a group of three cells along the x-direction. In the case for sick cells when d equals 10 lattice units (5 μm), we can see that the mass center stays stagnant around the  Figure 11A. While Figure 12A-E shows the three-dimensional evolutions of the deformations of three healthy RBCs in a constricted microchannel where the narrowest diameter is of 5 μm, Figure 12F prints three infected cells stuck at the constriction segment of the same width. After comparing with figures 10 and 11 of Wang's, we can see that, although the simulation setup is different in the starting gesture of the cells and the shape of channel constriction, the behavior of cell dynamics and the degree of passage time in 2D and 3D are qualitatively similar. Narrower constriction of the channel impedes the cell motion as being expected biophysically. Hardened cells take longer time to proceed and get blocked at a width that healthy ones can push through.

CONCLUSION
By three-dimensional IB-LBM, we succeeded in simulating healthy and hardened red blood cells passing through a stenosed microvessel in a Poiseuille flow where the out-of-plane transient motions of the cells are stepwise followed.
We show that, when an RBC gets hardened because of the infection by a disease like malaria, it offers less change in its bending energy, and takes much more time to overcome the stenosis. On the other hand, deterioration of the viscosity contrast has a minor impact on the transverse time. The different degree of deformability of RBCs can act as a metric to differentiate or filter hardened cells from healthy ones. In extending part of Wang's work on a group of three RBCs passing through stenosed microchannel into a three-dimensional scenario, we find that the a group of three hardened cells become stagnant when the narrow width of the stenosis equals to 5 μm comparably. The transient out-of-plane motion is, however, closer to the reality where biconcave, slipper, parachute, and other transient shapes are present.
Our next plan is to simulate the dynamics of dense RBCs passing through a stenosed vessel in a high hematocrit blood flow. When the computation domain gets more extensive in the current work, each simulation loop takes a considerable amount of time, barely affordable on a multicore Intel i7 desktop. Since LBM has the natural advantage of parallelization, it offers an anticipated opportunity to deploy the corresponding capability of modern GPUs to reduce high computational demand.

PEER REVIEW INFORMATION
Engineering Reports thanks the anonymous reviewers for their contribution to the peer review of this work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon request.