# Two-dimensional strain-hardening membrane model for large deformation behavior of multiple red blood cells in high shear conditions

- Swe Soe Ye
^{1}, - Yan Cheng Ng
^{1}, - Justin Tan
^{1}, - Hwa Liang Leo
^{1}and - Sangho Kim
^{1, 2}Email author

**11**:19

https://doi.org/10.1186/1742-4682-11-19

© Ye et al.; licensee BioMed Central Ltd. 2014

**Received: **13 December 2013

**Accepted: **12 March 2014

**Published: **13 May 2014

## Abstract

### Background

Computational modeling of Red Blood Cell (RBC) flow contributes to the fundamental understanding of microhemodynamics and microcirculation. In order to construct theoretical RBC models, experimental studies on single RBC mechanics have presented a material description for RBC membranes based on their membrane shear, bending and area moduli. These properties have been directly employed in 3D continuum models of RBCs but practical flow analysis with 3D models have been limited by their computationally expensive nature. As such, various researchers have employed 2D models to efficiently and qualitatively study microvessel flows. Currently, the representation of RBC dynamics using 2D models is a limited methodology that breaks down at high shear rates due to excessive and unrealistic stretching.

### Methods

We propose a localized scaling of the 2D elastic moduli such that it increases with RBC local membrane strain, thereby accounting for effects such as the Poisson effect and membrane local area incompressibility lost in the 2D simplification. Validation of our 2D Large Deformation (2D-LD) RBC model was achieved by comparing the predicted RBC deformation against the 3D model from literature for the case of a single RBC in simple shear flow under various shear rates (dimensionless shear rate G = 0.05, 0.1, 0.2, 0.5). The multi-cell flow of RBCs (38% Hematocrit) in a 20 μm width microchannel under varying shear rates (50, 150, 150 s^{-1}) was then simulated with our proposed model and the popularly-employed 2D neo-Hookean model in order to evaluate the efficacy of our proposed 2D-LD model.

### Results

The validation set indicated similar RBC deformation for both the 2D-LD and the 3D models across the studied shear rates, highlighting the robustness of our model. The multi-cell simulation indicated that the 2D neo-Hookean model predicts noodle-like RBC shapes at high shear rates (G = 0.5) whereas our 2D-LD model maintains sensible RBC deformations.

### Conclusion

The ability of the 2D-LD model to limit RBC strain even at high shear rates enables this proposed model to be employed in practical simulations of high shear rate microfluidic flows such as blood separation channels.

## Keywords

## Introduction

The transport behavior of blood in microcirculatory flows can be characterized by the mechanical response of the two major constituents of the blood mixture to the fluidic stresses driving the bulk flow. The first major constituent of blood is plasma, which under physiological conditions has Newtonian properties and a viscosity similar to water. The second major component is the red blood cells (RBCs) that make up about 35% to 45% of the systemic blood volume for an average individual [1]. The RBC phase contributes significantly to the complex behavior of blood in micro-flows, such as shear-thinning, the Fahraeus effect and the Fahraeus-Lindqvist effect. Due to these significant contributions to blood microrheology, the mechanical characteristics of RBCs have been studied extensively. The shape of the RBC can be defined by a membrane that separates the inner fluid (cytoplasm) from the suspending fluid (blood plasma). The most notable properties of the RBC membrane are its hyperelastic and viscoelastic response to high shear stress, membrane area incompressibility and the ability to recover its initial shape with the removal of external stress [2, 3].

In accordance with these properties, many previous studies [4–12] have been undertaken to describe the mechanical behavior of RBCs *in silico*. In these studies, the various numerical discretization techniques for the RBC model range from the mesoscale approach with cytoskeletal network models [4, 5] and particle method models [6] to the macroscopic approach of viscoelastic spring network models [7, 8], finite element method (FEM) models [9] and the boundary integral method (BIM) models [10]. With regards to the macroscopic RBC modeling approach, most previous studies have assumed either the neo-Hookean [10–13] or the Skalak constitutive relations [3, 14–16] to describe the non-linear stress-strain shear response of the RBC membrane deformation. In general, the membrane shear response is considered to be the most dominant deformation modality in the RBC membrane mechanical response. Although 3D simulations based on the macroscopic RBC models have previously been performed, the usefulness of such 3D simulations may be very limited due to the extremely high computational cost; consideration of the RBC interactions (aggregation and disaggregation) in high hematocrit flows for a 3D simulation is not feasible without employing sophisticated parallel computing techniques. Consequently, many previous numerical studies have instead utilized 2D RBC models to simulate physiological blood flows [6–8, 13, 17].

In the 2D modeling approach, the neo-Hookean and Skalak constitutive relations have been reformulated for 2D by removing a principle strain direction from the original 3D formulation [18]. However, employing the 2D formulations without modifying the effective moduli can overpredict the extent of deformation in the RBC membrane due to the disappearance of the Poisson’s effect contributed by the second principle strain direction. Furthermore, one important membrane feature that has been considered in the 3D simulation but not in the 2D simulation is the surface area incompressibility of the RBC membrane attributed by its incompressible lipid bilayer. Essentially, the extensional resistance of the RBC cross-section is not entirely a result of the membrane’s shear resistance. Consequently, the significant loss of these two deformation modalities in the 2D simulation severely limits the accurate prediction of the RBC 2D cross-sectional profiles under complex flow conditions such as high shear rates, crowded cell-cell interactions under high hematocrits and multi-directional RBC strain. This inaccuracy has limited 2D studies in the literature to the low flow regime models where shear rates are typically less than 300 s^{-1}[7, 8].

In the present study, we propose a modification of the 2D neo-Hookean relation to compensate for the apparent softening of the RBC membrane in 2D. The modified membrane model is coupled with the lattice Boltzmann method (LBM) flow solver through the immersed boundary method (IBM) [19]. For the membrane model development, a large deformation scaling coefficient is applied to the neo-Hookean membrane model and its several complementary constitutive relations to introduce a strain hardening effect on the RBC membrane at high strain rates.

## Numerical model

### Lattice Boltzmann method (LBM)

LBM simplifies the original Boltzmann equation by discretizing time, space and momentum [20, 21] through the employment of a lattice grid. Its mesoscopic nature arises from the fact that it considers microscale kinetic conditions of the fluid particles in relation to the macroscopic variables such as continuum mass and momentum. In a 2D discretization of space, the microfluxes are quantified on the square lattice in 9 directions under the D2Q9 approach [22]. The key concept of the LBM approach is that the microstates in the mesosystem evolve as a result of the macroscale (continuum) conditions and the evolution is conducted in key stages known as the streaming and collision stages. To represent the statistical contribution of the 9 directions, the objective quantifiers of the microstates are given by the density distribution function in the LBM formulation, summation of the 9 density distribution functions at a lattice grid point gives the local fluid density of the macroscopic continuum.

*t*, Δt is the lattice time step,

*τ*is the relaxation time and

*B*

_{ i }is the body force term discretized in the 9 lattice directions denoted by the subscript

*i*. In Eq. (1), left-hand side terms represent the streaming stage for the 9 density distribution functions

*f*denoted by the subscript

*i*. The first right-hand side term represents the collision contribution to the distribution functions [24] whereby microstates are disturbed from their equilibrium states which is given by the equilibrium density distribution function ${\mathit{f}}_{\mathit{i}}^{\mathit{eq}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)$:

*ω*

_{ i }is the weight factor, taking the value of 4/9 for

*ω*

_{0}, 1/9 for

*ω*

_{1-4}, and 1/36 for

*ω*

_{5-8}. The lattice speed of sound is given by the form ${\mathit{c}}_{\mathit{s}}=\mathit{h}/\left(\sqrt{3}\mathit{\Delta t}\right)$ where

*h*is the lattice cell (space step) size. The second right-hand side term, which is the aforementioned body force term, includes all the fluid-structure interaction (FSI) forces between the suspending fluid and the RBC membrane. Calculation of the body force term is given as:

where ${\overrightarrow{\mathit{F}}}_{\mathit{f}}$ is the total FSI body force acting on the lattice (fluid) grid node due to on-membrane forces.

*ρ*and continuum velocity $\overrightarrow{\mathit{u}}$ can be obtained from the fundamental equivalency between the macrostate quantities and the summation of microstate quantities (conservation of mass and momentum):

where $\mathit{v}=\left(\mathit{\tau}-\frac{1}{2}\right){\mathit{c}}_{\mathit{s}}^{2}\mathrm{\Delta t}$ and $\mathit{P}={\mathit{c}}_{\mathit{s}}^{2}\mathit{\rho}$.

### Immersed boundary method (IBM)

The IBM [26] was employed in our simulation to account for the acceleration effect of a moving boundary on the fluid through the application of a distributed force density evaluated from the boundary’s constitutive laws. In this method, the fluid domain is represented by an Eulerian mesh where the globally referenced coordinates of the fluid grid point are given by $\overrightarrow{\mathit{x}}$. The RBC membrane boundary is represented by a Lagrangian mesh with body-fixed coordinates $\overrightarrow{\mathit{s}}$ used for the membrane force calculations. The globally referenced coordinate location of a node on the moving membrane can be given by $\overrightarrow{\mathit{X}}\left(\overrightarrow{\mathit{s}},\mathit{t}\right)$.

*m*and

*f*denote the membrane node and fluid grid point indices, respectively. It is given by the discrete delta function:

where ${\mathit{\varphi}}_{\mathit{m}}\left({\overrightarrow{\mathit{x}}}_{\mathit{f}}-{\overrightarrow{\mathit{X}}}_{\mathit{m}}\right)$ is the spread function, given by the same discrete delta function in Eq. (8).

The force density spreading and membrane node velocity interpolation are performed on a 4*h* × 4*h* region [26]. To describe the different properties of blood plasma and cytoplasm within the RBCs, an indicator field approach was employed to update the moving topology of the plasma domain and RBC interior cytoplasm domain at every time-step. By utilizing our recently developed method (flood-fill method) [27], we have assigned a viscosity of 6.0 and 1.2 cP to the cytoplasm and plasma respectively, thereby capturing the viscoelastic response of the membrane deformation due to the RBC interior-exterior viscosity ratio. The details of the flood-fill algorithm used to update the fluid properties during the simulations can be found in our earlier work [27].

### RBC model

#### Membrane in-plane shear

where *E*_{s} is the shear elastic modulus of the membrane and *λ* is the stretch ratio on the local membrane segment given by the ratio of the current membrane segment length *l*_{
m
} over the initial membrane segment length *l*_{0}. The shear elastic modulus in this study is set at 6 × 10^{-3} dyn cm^{-1}, which is within the physiological range for RBC elastic properties [30].

As discussed earlier, the limitation of applying the 2D neo-Hookean model is the unrestricted stretching of membrane perimeter (circumference). Previous studies [31, 32] have shown that use of the neo-Hookean model alone is unable to restrict the membrane surface area changes in 3D capsules from exceeding 7.8% at shear rates above 300 s^{-1}. Conversely, the changes in global RBC membrane surface area should be less than 5% under physiological conditions due to the incompressibility of the lipid bilayer in the membrane [31]. The stretching of the RBC membrane is expected to be overpredicted for the case of 2D simulations where the extension of a 1D surface (line) is unrestricted due to the lack of the Poisson’s effect from the second principle strain direction.

*α*, which is a function of the local RBC membrane stretch ratio

*λ*. The 2D-LD neo-Hookean model developed in this study can be presented as follows:

*λ*) is given by:

where *D*_{
LD
} and *β*_{
LD
} are constants. It is of note that the value of *α* from Eq. (12) approaches unity at very low strains and Eq. (11) reverts back to its original form in Eq. (10) under such conditions.

#### Membrane bending

where *E*_{b} is the bending modulus of the RBC membrane, *κ* is the current membrane curvature and *κ* is the spontaneous curvature of the un-deformed RBC. The bending force is similarly scaled by the large deformation coefficient *α* presented in Eq. (12). The scaling of the membrane flexural resistance is necessary to prevent membrane buckling under high compression which instigates numerical instabilities.

#### Cytoplasmic volume conservation

*p*

_{ int }is introduced to act on the RBC membrane, thereby strictly imposing the conservation of the RBC internal volume and mass. The pressure penalty model for a 2D capsule can be expressed as follows:

*k*

_{ p }is the incompressibility coefficient and the argument

*λ*

_{ C }is the RBC perimeter extension ratio given by ratio of the current RBC circumference

*L*

_{ C }over the initial circumference of the circular RBC

*L*

_{0}.

*A*is the internal area of the deformed RBC, and

*A*

_{ref}is the initial internal area of the RBC. The internal area of the RBC is calculated using Green’s theorem:

where *x*_{
m
} and *y*_{
m
} are the coordinates of the points on the RBC membrane curve C.

By taking a sufficiently large incompressibility coefficient *k*_{
p
} and by considering the growing restriction under large deformation using the α(*λ*_{
C
}) term, we can engage a sufficiently large internal pressure *p*_{
int
} to maintain the constant RBC size in the simulation. The maintenance of a constant RBC area is a necessary constraint in order to satisfy the conservation of cytoplasmic mass in the channel flow. Accordingly, the cytoplasmic mass is not allowed to arbitrarily swell or disappear from the movement of the RBC membrane. Consequently, by including the cytoplasmic conservation in our model, the RBC area and 2D hematocrit in the channel can be maintained at a constant value throughout the entire simulation.

#### RBC-RBC interaction

*φ*can be expressed using the Morse-type potential energy function [34]:

*r*[μm] is the separation distance between the pairing membrane nodes and

*r*

_{0}[μm] is the zero force distance specified in the model.

*D*

_{ e }[μJ μm

^{-2}] is the surface energy and

*β*[μm

^{-1}] is the scaling factor that determines the rate of interaction energy decay with distance. In this study,

*r*

_{0},

*β*and

*D*

_{ e }were set with the values of 0.49, 3.84 and 1.3 × 10

^{-7}respectively as reported in previous studies [27, 35]. The total interaction force between the membrane nodes is expressed as the negative derivative of the interaction potential from Eq. (16):

In Eq. (17), a negative *F*_{
αgg
} value when *r* > *r*_{0} indicates an aggregating (attraction) force while a positive value when *r* ≥ *r*_{0} represents a repulsion force. Additionally, the LD coefficient scales only the repulsion force between pairing RBC membranes to prevent cell to cell overlap from the increase in internal forces from the shear, bending and dilatory modalities.

*r*

_{0}, the node-pair experiences a repulsion force. However, when the distance between the two nodes is within the

*R*to

*r*

_{0}range, the node-pair experiences an attraction force.

### Simulation setup

*u*=

*ky*is imposed where the strain rate

*k*can be obtained by the simple relation

*k = U/Y. Y*is the half-height of the simulation domain and

*U*is the maximum magnitude of the velocity at the top and bottom of the simulation domain as presented in Figure 3. The deformation of the capsule is described by the Taylor deformation index

*D*

_{xy}which is given by:

*L*is the major diameter of the RBC and

*B*is the minor diameter as shown in Figure 3. Notably, this characterization of the RBC deformation only works for RBCs that adopt a 2D ellipse profile and the value of

*D*

_{ xy }is highly sensitive to the major and minor diameters at low deformation states.

**Multi-cell channel flow simulations**

Case | LD scaling applied | |||
---|---|---|---|---|

Shear | Bending | Vol. | Cell-cell | |

conservation | interaction | |||

I | X | X | X | X |

II | √ | X | X | X |

III | √ | √ | √ | √ |

In the initial condition, twelve circular RBCs were suspended in a periodic arrangement inside a channel of 80 μm by 20 μm to achieve a physiological hematocrit level (38%). The circular RBC cross-sectional profile was chosen for simulation as this 2D profile represents the most extreme shearing orientation for the RBCs in a narrow channel. Pressure boundary conditions were prescribed for the pressure-driven flow to obtain the pseudoshear rates (mean velocity/channel width) of 50, 150, and 500 s^{-1} for each of the three cases. Periodic translations were implemented on the RBCs at the inlet and outlets such that RBCs leaving the simulation domain re-enter from the inlet, thereby maintaining the same number of 12 RBCs for the entirety of the simulation.

## Results and discussion

### 2D large deformation (2D-LD) model validation

*G*was used.

*G*provides a normalized indication of the stress condition on the cell by comparing the estimated fluidic shear stress applied on the RBC membrane (numerator) to the inherent elastic property of the membrane (denominator) [12]:

*μ*is the dynamic fluid viscosity,

*k*is the shear rate, and a is the equivalent radius of the RBC. Breyiannis and Pozrikidis [11] have compared the deformation of 2D solitary capsules against the deformation of 3D spherical capsules and have reported a good correlation. They established a

*D*

_{xy}correspondence between the

*G*values for circular capsules and 3D spherical capsules by using their cross-sectional profiles. Consequently, the empirical equation relating the 3D

*G*to 2D

*G*was reported to be:

*D*

_{xy}values obtained for the dimensionless shear rates

*G*of 0.05, 0.1, 0.2 and 0.5, we observed a reasonable correspondence between our 2D results and the 3D model results of Eggelton and Popel. While the discrepancy is close to 50% at the lowest shear rate, the 2D-LD model can sufficiently limit the RBC deformation to agree with the 3D model data at higher shear rates. It is likely that the low shear rate discrepancy arises as a result of Eq. (20) presenting non-sensible G values for the 2D equivalent at very low shears. For example, conversion of the 3D G at a value of 0.01 using Eq. (20) results in a 2D G value of -0.00383. Thus, this conversion may not be accurate under very low shear conditions.

### RBC deformation in a multi-cell channel flow

*ε*and the earlier introduced perimeter extension ratio

*λ*

_{ C }to describe the overall deformation of the cells in the channel flow;

*ε*is given by:

#### Model comparison, 2D-LD against 2D-neo-Hookean

*λ*

_{ C }. At 50 s

^{-1}, there was no statistical difference among the three cases (

*λ*

_{ C }= 1.049 ± 0.013 for Case I, 1.048 ± 0.002 for Case II, and 1.044 ± 0.043 for Case III). However, at 150 s

^{-1}, there was a ~13% difference (

*P*< 0.001) in the average extension between Case I (1.225 ± 0.070) and Case III (1.083 ± 0.033), but no significant difference between Case II (1.088 ± 0.037) and Case III. Similarly at 500 s

^{-1}, there was a ~135% difference (

*P*< 0.001) between Case I (2.688 ± 0.835) and Case III (1.125 ± 0.040), but no statistical difference between Case II (1.125 ± 0.031) and Case III. The pronounced difference in perimeter extension between Case I and the other two cases was expected since the LD model imposes a larger restrictive force on the membrane when it stretches beyond a stipulated limit. Thus, even in a very high shear condition of 500 s

^{-1}, the RBC perimeter does not extend by more than 12% of the original length for the 2D-LD model in Case III whereas Case I’s RBCs have stretched in length by more than two times of their original perimeter. Figure 6 shows the deformation profiles of the RBCs for Case I – III in the channel flow at a particular instant in time. As observed in Figure 6a, all three cases were initialized from the same symmetrical arrangement but the RBC flow developed differently with time (Figures 6b–d). Due to the over-extension of RBCs in the simulation, the RBC flow for Case I never reached a developed flow condition for the simulation conducted at the highest shear condition of 500 s

^{-1}(Figure 6b). Subsequently, simulation failure occurred before the RBC flow structure could break its initial symmetric arrangement which typically occurs within 0.1 s of RBC flow as observed for Cases II and III in Figures 6c and d. Interestingly, the deformation profiles of cells observed in Case I for Figure 6b indicate that the extensive stretching of RBCs into “noodle-like” profiles occurs predominately for cells located in the high shear rate regions near the channel walls. From this evaluation of Case I’s result, it can be concluded that the 2D neo-Hookean model has a limitation in performing RBCs flows at high shear rates typical to microfluidic devices (> 1000 s

^{-1}). Through a comparative investigation of multi-cell simulations with (Cases II and III) and without (Case I) LD augmentation, we have established that LD augmentation is required for 2D RBC models to maintain physiological 2D RBC deformations, particularly for the cells travelling in close proximity to or impinged against the channel wall (see Figure 6b). A very recent study [13] on 2D multiple RBC flow simulations in a bifurcation also showed this limitation of the 2D RBC deformation simulation. They have illustrated that even at 100 s

^{-1}, non-physiologically over-stretched RBC shapes (“noodles”) were obtained in the model simulation due to wall impingement and multi-directional strain near flow bifurcation corners, similar “noodling” of RBCs

*in vivo*has not been reported.

While the presentation of deformation data in Figure 5 indicates that the overall perimeter remains statistically the same in Case II and Case III, the shapes of the RBCs were considerably different in these two LD-applied cases as evidenced by the images in Figures 6c and d. This may be due to the lack of strain-hardening on the remaining three constitutive models for the RBCs in Case II (as summarized in Table 1). The implication of this omission in Case II will be discussed in the following subsection.

#### Significance of bending resistance and contact forces for large deformations

*α*used in the 2D-LD model is a simple multiplicative operator that can be used to apply the same order of strain hardening to all elastic moduli involved in the RBC membrane’s constitutive response to deformation.

#### Cell-free layer width and relative apparent viscosity

The cell-free layer (CFL) and its role in influencing the apparent viscosity of blood is an important characteristic in quantifying microvessel and microchannel flows. Due to the shear-induced migration of RBCs towards the center of the vessel, the formation of a CFL along the vessel wall significantly lowers the apparent viscosity of blood in microvessel flows when compared against the uniform bulk viscosity of blood [1]. Accordingly, we validated our channel flow simulations by comparing the CFL width and the apparent viscosity predictions of our 2D-LD model in Case III (for the pseudoshear rate of 50 s^{-1}) against the earlier work of Zhang and coworkers [29, 30]. Our predicted CFL width was ~26% of the total channel width which is in good agreement with the value (27% – 32% at 58 s^{-1}) reported in their study [29] where the simulation condition was similar to ours.

*μ*

_{ app }of blood in our channel flow simulations was calculated using the Poiseuille formula:

*P*is the pressure difference applied across the channel length

*L*

_{ channel },

*H*is the channel width and

*Q*is the resulting flow rate. For comparison against the literature, the apparent viscosity was normalized by the plasma viscosity

*μ*

_{ plasma }to provide the relative apparent viscosity

*μ*

_{ rel }:

Our result (1.10) falls within the range of the simulated results by Zhang and colleagues: 1.05 in a 20-μm channel [30] and 1.29 in a 12-μm channel [29]. Although the comparisons of the CFL width and the relative apparent viscosity indicate reasonable agreement between our results and theirs, it should be noted that circular RBC profiles were considered for our 2D flow model while they have represented the 2D flow of RBCs using biconcave RBC profiles. Subsequently, even though we have a higher 2D hematocrit of 38% in comparison to their 30.5% [29], our actual number of RBCs in the simulation is much fewer (12 circular RBCs vs. 27 biconcave RBCs). As a result of this, it may be limited to directly compare our relative apparent viscosity and cell-free layer width with the values reported by them.

*λ*

_{ C }. This is in accordance with the result shown in Figure 9 where the RBC deformation increases when the distance between the RBC and the channel wall is reduced.

#### RBC stress and strain

*G*using the pseudoshear rate for

*k*in Eq. (19) could be a poor representation of the actual shear condition applied on the RBCs in a channel flow. Therefore, in this study, we introduce a new dimensionless shear rate parameter

*G**to provide a better approximation of the local shear condition in the channel flow as follows:

where *τ*_{max} is the time-averaged maximum shear stress on the RBC membrane and it is a function of the the minimum RBC membrane to channel wall distance ${\overrightarrow{\mathit{r}}}_{\text{min}}$.*τ*_{max} was obtained by recording the maximum shear stress exerted on the RBC membrane for each successive time-step and thereafter performing a time-averaging calculation on the maximum stress over the period of analysis. In the present study, we have analyzed the RBC mechanics over 10 material transit cycles. The material transit time represents the average time it takes for a cell entering the simulation domain to exit at the channel outlet and this was estimated using the bulk flow velocity and the channel length. Since both *τ*_{max} and the cell perimeter *L*_{
C
} are dependent on the RBC trajectory in the channel, *G** can provide a closer approximation of the spatially and temporally varying shear condition than the traditional *G* for a channel flow study.

*G**, the deformation state of RBCs in a general flow condition can be characterized and compared by using the relation between

*G**and the membrane circumferential strain

*ε*. To prove the validity of the

*G**and membrane strain relation, both parameters were calculated for the two different simulation sets performed in this study. In Figure 10, the relation between the two parameters is compared for the single-cell in simple shear and the multi-cells in channel flow conditions. An exponential curve (

*ε*= – 0.1771

*e*

^{-0.0611G*}+ 0.1771) was fitted by non-linear regression to both the single cell deformation results and multi-cell deformation results to obtain the empirical expression for circumferential strain against varying

*G*.*The regression curve produced high

*R*

^{ 2 }values of 0.996 and 0.939 for the simple shear and channel flow data sets respectively, indicating that the two different shear conditions can be described by the same relation between

*G**and

*ε*. This implies that by using

*G**and

*ε*, we can relate the single-cell simple shear flow results to the multi-cell channel flow results.

On the other hand, using the original dimensionless shear rate *G* to relate between the two sets of flow conditions demonstrates the limitation of the *G* parameter as a universal tool for relating the applied shear to the RBC elastic resistance. RBCs in a channel flow under the pseudoshear rate of 150 s^{-1} (*G* = 0.125) exhibited a wide range of deformation *ε* (0.05 – 0.11) that is dependent on the cell location as seen in Figure 9 where the RBC perimeter extension ratio *λ*_{
C
} was plotted against RBC center of mass location. Corresponding to this range of *ε*, the RBC in a simple shear flow exhibited similar values of 0.05 and 0.11 when the applied *G* was 0.2 and 1.0, respectively. Hence, it can be concluded that while the stress-strain relation between *G* and *ε* is valid for the simple shear flow, it is weak for the case of channel flows since a single value of *G* presents multiple strain possibilities depending on the cell location in the channel (Figure 9).

The case for using *G** presented in Figure 10 to signify RBC stress-strain behavior becomes particularly favorable when we wish to provide a more accurate estimation of shear stress acting on the RBC at high shear rates. Using Eq. (24), we can obtain *τ*_{max} acting on the RBC in experiments by obtaining circumferential strain *ε* and RBC perimeter *L*_{c} through imaging and measurement techniques and *G** by reading off its corresponding value against *ε* in Figure 10. Therefore, Figure 10 and Eq. (24) can be used to predict the maximum shear stress on an RBC for a given profile in the channel flow for experimental studies where shear stress cannot be obtained through direct measurement of the shear stress quantity. For our particular study, the maximum pseudoshear rate considered was 500 s^{-1} for a channel of 20-μm width. The corresponding time-averaged maximum shear stress on the RBCs in the flow was approximately 10 Pa (*G** = 47). This is well below the reported shear stress (300 Pa for 120 s of shear) for RBC lysis [40] or even sub-hemolytic damage to the RBC membrane [41]. Hence, even though strain hardening is expected of the RBC at shear rates > 500 s^{-1}, this is by no means suggestive of mechanical damage to the RBC.

### Potential limitations of the 2D-LD model

One possible limitation of our 2D-LD model can arise from the chosen 2D RBC profile for the channel flow simulation. Firstly, a similar simulation on the 2D biconcave RBC profile may give a different set of results due to its higher bending and flexing capability than the circular RBC. While the LD scaling model will not be different in form for biconcave cells, a calibration of the model terms and coefficients would be required. Furthermore, we have assumed that the cross-sectional area remains constant in the 2D plane of investigation. This is necessary for the 2D model to maintain the channel hematocrit and to enforce the 2D conservation of mass. This model is therefore limited to flow situations where 2D RBCs remain in the plane of observation. However, such flow conditions can easily be found in most microfluidic flows where the Reynolds number is very low.

Additionally, our present 2D-LD model lacks the inclusion of the membrane viscosity and its effect on the dynamic deformation behavior of RBCs. It is likely that the membrane viscosity will affect the dynamic behavior of RBCs that are subjected to ever-evolving shear rates due to the variation in their transverse location as they travel along their respective trajectories within the channel. Membrane viscosity is likely to delay the deformation response of RBCs to fluctuations in the shearing condition as a result of the changes in RBC position and orientation. However, our present model does include the effect of cytoplasmic viscosity (6.0 cP) and the plasma viscosity (1.2 cP), and our earlier work [27] has shown that the inclusion of the difference in viscosity between the two fluids can similarly delay the RBC membrane deformation response.

## Conclusion

In the present study, we have presented a 2D large deformation (2D-LD) model to augment the elastic moduli of the RBC membrane in the high shear rate flow regimes. The efficacy of the model and the predictive accuracy of the resulting 2D deformation states were tested on a single circular RBC profile under a simple shear condition and the results were found to be in good agreement with the 3D data. Furthermore, this study highlights the importance of sufficiently scaling the various membrane mechanics models to prevent numerical instabilities in the simulation. In our analysis of the stress-strain relation for the membrane, we have also proposed a new dimensionless shear rate term *G** to generalize the shear condition on a RBC so as to predict the extent of deformation regardless of flow conditions.

Therefore, our 2D-LD model can be applied to blood flows in practical microfluidic studies involving channel bifurcations and cell mechanical partitioning [37] where high and multi-direction strain can be applied on the RBCs at flow dividing locations. These studies would need robust mechanical models to predict the RBC deformation without incurring the high computational cost that is generally required for 3D simulations. Thus, by utilizing our LD model, it would be possible to simulate a blood flow in a microfluidic system, and such a model would enable us to optimize a microfluidic channel structure for biomedical applications at relatively low computing cost.

## Declarations

### Acknowledgements

The authors would like to thank the kind reviewers and the editors at *TBioMed* who have invested their time and energy to assist us in improving the quality of this paper with their invaluable comments. This research was supported and made possible by FRC R-397-000-134-112.

## Authors’ Affiliations

## References

- Popel AS, Johnson PC: Microcirculation and hemorheology. Annu Rev Fluid Mech. 2005, 37: 43-69. 10.1146/annurev.fluid.37.042604.133933.PubMed CentralView ArticlePubMedGoogle Scholar
- Hochmuth RM, Waugh RE: Erythrocyte membrane elasticity and viscosity. Annu Rev Physiol. 1987, 49: 209-219. 10.1146/annurev.ph.49.030187.001233.View ArticlePubMedGoogle Scholar
- Evans EA, Hochmuth RM: Membrane viscoelasticity. Biophys J. 1976, 16: 1-11.PubMed CentralView ArticlePubMedGoogle Scholar
- Gov NS: Active elastic network: cytoskeleton of the red blood cell. Phys Rev E Stat Nonlin Soft Matter Phys. 2007, 75: 011921-View ArticlePubMedGoogle Scholar
- Hansen JC, Skalak R, Chien S, Hoger A: An elastic network model based on the structure of the red blood cell membrane skeleton. Biophys J. 1996, 70: 146-166. 10.1016/S0006-3495(96)79556-5.PubMed CentralView ArticlePubMedGoogle Scholar
- Tsubota K, Wada S, Yamaguchi T: Particle method for computer simulation of red blood cell motion in blood flow. Comput Meth Programs Biomed. 2006, 83: 139-146. 10.1016/j.cmpb.2006.06.005.View ArticleGoogle Scholar
- Secomb T, Styp-Rekowska B, Pries A: Two-dimensional simulation of red blood cell deformation and lateral migration in microvessels. Ann Biomed Eng. 2007, 35: 755-765. 10.1007/s10439-007-9275-0.View ArticlePubMedGoogle Scholar
- Secomb TW: Mechanics and computational simulation of blood flow in microvessels. Med Eng Phys. 2011, 33: 800-804. 10.1016/j.medengphy.2010.09.016.PubMed CentralView ArticlePubMedGoogle Scholar
- Dao M, Lim CT, Suresh S: Mechanics of the human red blood cell deformed by optical tweezers [Journal of the Mechanics and Physics of Solids, 51 (2003) 2259-2280]. J Mech Phys Solids. 2003, 53: 493-494.View ArticleGoogle Scholar
- Pozrikidis C: Numerical simulation of the flow-induced deformation of red blood cells. Ann Biomed Eng. 2003, 31: 1194-1205.View ArticlePubMedGoogle Scholar
- Breyiannis G, Pozrikidis C: Simple shear flow of suspensions of elastic capsules. Theor Comput Fluid Dyn. 2000, 13: 327-347.View ArticleGoogle Scholar
- Ramanujan S, Pozrikidis C: Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of fluid viscosities. J Fluid Mech. 1998, 361: 117-143.View ArticleGoogle Scholar
- Yin XW, Thomas T, Zhang JF: Multiple red blood cell flows through microvascular bifurcations: cell free layer, cell trajectory, and hematocrit separation. Microvasc Res. 2013, 89: 47-56.View ArticlePubMedGoogle Scholar
- Evans EA: A new material concept for the red cell membrane. Biophys J. 1973, 13: 926-940. 10.1016/S0006-3495(73)86035-7.PubMed CentralView ArticlePubMedGoogle Scholar
- Evans EA: New membrane concept applied to the analysis of fluid shear- and micropipette-deformed red blood cells. Biophys J. 1973, 13: 941-954. 10.1016/S0006-3495(73)86036-9.PubMed CentralView ArticlePubMedGoogle Scholar
- Skalak R, Tozeren A, Zarda RP, Chien S: Strain energy function of red blood cell membranes. Biophys J. 1973, 13: 245-264. 10.1016/S0006-3495(73)85983-1.PubMed CentralView ArticlePubMedGoogle Scholar
- Barber J, Alberding J, Restrepo J, Secomb T: Simulated two-dimensional red blood cell motion, deformation, and partitioning in microvessel bifurcations. Ann Biomed Eng. 2008, 36: 1690-1698. 10.1007/s10439-008-9546-4.PubMed CentralView ArticlePubMedGoogle Scholar
- Sui Y, Chew YT, Low HT: A lattice Boltzmann study on the large deformation of red blood cells in shear flow. Int J Modern Phys C. 2007, 18: 993-1011. 10.1142/S012918310701108X.View ArticleGoogle Scholar
- Feng Z-G, Michaelides EE: The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems. J Comput Phys. 2004, 195: 602-628. 10.1016/j.jcp.2003.10.013.View ArticleGoogle Scholar
- He X, Luo L-S: A priori derivation of the lattice Boltzmann equation. Phys Rev E. 1997, 55: R6333-R6336. 10.1103/PhysRevE.55.R6333.View ArticleGoogle Scholar
- Abe T: Derivation of the lattice Boltzmann method by means of the discrete ordinate method for the Boltzmann equation. J Comput Phys. 1997, 131: 241-246. 10.1006/jcph.1996.5595.View ArticleGoogle Scholar
- Chopard B, Masselot A: Cellular automata and lattice Boltzmann methods: a new approach to computational fluid dynamics and particle transport. Futur Gener Comput Syst. 1999, 16: 249-257. 10.1016/S0167-739X(99)00050-3.View ArticleGoogle Scholar
- Guo Z, Zheng C, Shi B: Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys Rev E Stat Nonlin Soft Matter Phys. 2002, 65: 046308-View ArticlePubMedGoogle Scholar
- Bhatnagar PL, Gross EP, Krook M: A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev. 1954, 94: 511-10.1103/PhysRev.94.511.View ArticleGoogle Scholar
- Buick JM, Martin AJ, Cosgrove JA, Greated CA, Easson WJ: Comparison of a lattice Boltzmann simulation of steep internal waves and laboratory measurements using particle image velocimetry. Eur J Mech B-Fluids. 2003, 22: 27-38. 10.1016/S0997-7546(02)00002-X.View ArticleGoogle Scholar
- Peskin CS: Numerical analysis of blood flow in the heart. J Comput Phys. 1977, 25: 220-252. 10.1016/0021-9991(77)90100-0.View ArticleGoogle Scholar
- Ju M, Ye SS, Low HT, Zhang J, Cabrales P, Leo HL, Kim S: Effect of deformability difference between two erythrocytes on their aggregation. Phys Biol. 2013, 10: 036001-10.1088/1478-3975/10/3/036001.PubMed CentralView ArticlePubMedGoogle Scholar
- Bagchi P, Johnson PC, Popel AS: Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow. J Biomech Eng. 2005, 127: 1070-1080. 10.1115/1.2112907.View ArticlePubMedGoogle Scholar
- Zhang J, Johnson PC, Popel AS: Effects of erythrocyte deformability and aggregation on the cell free layer and apparent viscosity of microscopic blood flows. Microvasc Res. 2009, 77: 265-272. 10.1016/j.mvr.2009.01.010.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang J, Johnson PC, Popel AS: An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows. Phys Biol. 2007, 4: 285-295. 10.1088/1478-3975/4/4/005.View ArticlePubMedGoogle Scholar
- Eggleton CD, Popel AS: Large deformation of red blood cell ghosts in a simple shear flow. Phys Fluids. 1998, 10: 1834-1845. 10.1063/1.869703.View ArticleGoogle Scholar
- Pozrikidis C: Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow. J Fluid Mech. 1995, 297: 123-152. 10.1017/S002211209500303X.View ArticleGoogle Scholar
- Neu B, Meiselman HJ: Depletion-mediated red blood cell aggregation in polymer solutions. Biophys J. 2002, 83: 2482-2490. 10.1016/S0006-3495(02)75259-4.PubMed CentralView ArticlePubMedGoogle Scholar
- Liu Y, Liu WK: Rheology of red blood cell aggregation by computer simulation. J Comput Phys. 2006, 220: 139-154. 10.1016/j.jcp.2006.05.010.View ArticleGoogle Scholar
- Zhang J, Johnson PC, Popel AS: Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method. J Biomech. 2008, 41: 47-55. 10.1016/j.jbiomech.2007.07.020.PubMed CentralView ArticlePubMedGoogle Scholar
- Kim S, Kong RL, Popel AS, Intaglietta M, Johnson PC: Temporal and spatial variations of cell-free layer width in arterioles. Am J Physiol Heart Circ Physiol. 2007, 293: H1526-H1535. 10.1152/ajpheart.01090.2006.View ArticlePubMedGoogle Scholar
- Bhagat A, Bow H, Hou H, Tan S, Han J, Lim C: Microfluidics for cell separation. Med Biol Eng Comput. 2010, 48: 999-1014. 10.1007/s11517-010-0611-4.View ArticlePubMedGoogle Scholar
- Yang S, Undar A, Zahn JD: A microfluidic device for continuous, real time blood plasma separation. Lab Chip. 2006, 6: 871-880. 10.1039/b516401j.View ArticlePubMedGoogle Scholar
- Alizadehrad D, Imai Y, Nakaaki K, Ishikawa T, Yamaguchi T: Parallel simulation of cellular flow in microvessels using a particle method. J Biomech Sci Eng. 2012, 7: 57-71. 10.1299/jbse.7.57.View ArticleGoogle Scholar
- Nevaril CG, Lynch EC, Alfrey CP, Hellums JD: Erythrocyte damage and destruction induced by shearing stress. J Lab Clin Med. 1968, 71: 784-790.PubMedGoogle Scholar
- Baskurt OK, Uyuklu M, Meiselman HJ: Protection of erythrocytes from sub-hemolytic mechanical damage by nitric oxide mediated inhibition of potassium leakage. Biorheology. 2004, 41: 79-89.PubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.