Axial and Nonaxial Migration of Red Blood Cells in a Microtube

Abstract

Human red blood cells (RBCs) are subjected to high viscous shear stress, especially during microcirculation, resulting in stable deformed shapes such as parachute or slipper shape. Those unique deformed RBC shapes, accompanied with axial or nonaxial migration, cannot be fully described according to traditional knowledge about lateral movement of deformable spherical particles. Although several experimental and numerical studies have investigated RBC behavior in microchannels with similar diameters as RBCs, the detailed mechanical characteristics of RBC lateral movement—in particular, regarding the relationship between stable deformed shapes, equilibrium radial RBC position, and membrane load—has not yet been fully described. Thus, we numerically investigated the behavior of single RBCs with radii of 4 $$\mathsf{\mu }$$m in a circular microchannel with diameters of 15 $$\mathsf{\mu }$$m. Flow was assumed to be almost inertialess. The problem was characterized by the capillary number, which is the ratio between fluid viscous force and membrane elastic force. The power (or energy dissipation) associated with membrane deformations was introduced to quantify the state of membrane loads. Simulations were performed with different capillary numbers, viscosity ratios of the internal to external fluids of RBCs, and initial RBC centroid positions. Our numerical results demonstrated that axial or nonaxial migration of RBC depended on the stable deformed RBC shapes, and the equilibrium radial position of the RBC centroid correlated well with energy expenditure associated with membrane deformations.

1. Introduction

The flow behavior of human red blood cells (RBCs) in capillaries has been intensively studied from various points of view, including shape, deformability, and physiological functions such as oxygen transport. Due to the large number of blood cells, the temporal shapes of individual RBCs during flow are of fundamental importance both in micro- and macro-scale hemorheology [1,2,3]. Especially in the microcirculation, where flow can be assumed to be almost inertialess, RBCs are subjected to high viscous shear stress, resulting a bistable shape ( parachute or slipper shape) in capillaries [4,5,6]. Recent numerical analysis further showed that the asymmetrical slipper shape of RBCs was observed not only in circular microchannels, but also in narrow rectangular microchannels, whose width was comparable to the thickness of an RBC [7]. The unique flow behavior of RBCs, which is often accompanied by axial or nonaxial migration, cannot be described using the traditional framework regarding lateral movement of deformable spherical particles, as originally reported in [8]. In this framework, a deformable spherical particle tends to move toward the channel axis and remains there. Assuming that the stable shape of RBCs under flow results from the force balance between internal/external hydrodynamic shear force and inner elastic force on the membrane, model analysis will provide insight into the mechanical background of both stable shape and equilibrium radial position.

So far, the impact on the stable RBC shape of various mechanical conditions, including flow speed, membrane elasticity, channel confinement, and cell volume–surface ratio, has been investigated using a 2D elastic spring model [9], a 2D vesicle model [1,2,10,11], 3D models such as vesicle model [12], a mesoscopic molecular dynamics model [13], and a capsule model [5,6]. These investigations showed that RBCs tend to form a stable parachute shape in a circular microchannel with a diameter comparable to that of an RBC (≤10 $\mathsf{\mu }$m) [6,12,13]. This shape was also observed in a narrow rectangular microchannel with a width comparable to the thickness of an RBC ≤ 3.5 $\mathsf{\mu }$m [7]. In a slightly larger microchannel (≥12 $\mathsf{\mu }$m), shape bistability became significant [5,13]. Guckenberger et al. (2018) further showed the effect of the initial position of the RBC centroid on the stable deformed configuration in Stokes flow for different capillary numbers $Ca$, which defines the ratio between a fluid viscous force and a membrane elastic force [5]. Despite these efforts, the relationships between those stable deformed shapes, the equilibrium radial positions of RBCs, and the membrane load have not yet been fully described.

Along with the aforementioned numerical studies, recent microfluidic techniques have allowed us to conduct high-throughput measurements of single-cell behavior under confined channel flow [14,15,16,17]. Since the stable shape of an RBC under flow is highly reliant on cell mechanical properties, quantifying cell shapes in microfluidic systems will be useful for understanding cell states—including the ability of RBCs to function as oxygen carriers—and might be extended to the diagnosis of blood diseases [18,19]. In particular, patients with sickle cell anemia have a high hemoglobin concentration that results in abnormal rheology [20,21,22]. Hence, alternations in membrane elasticity and the relationship between viscosity ratio and the steady shape of RBCs are clinically important. If the stable shape of RBCs is changed depending on viscosity ratios of the internal to external fluids of RBCs, the shape would be a hallmark to identify the cell state, or be useful in cellular-level diagnoses for blood diseases.

Therefore, the objective of this study is to reveal the relationship between stable deformed RBC shapes, their equilibrium radial positions, and the membrane loads of flowing RBCs in a microtube. We numerically investigated the behavior of a single RBC with a major diameter of 8 $\mathsf{\mu }$m in a straight circular microchannel with 15 $\mathsf{\mu }$m-diameter. The RBC was modeled as a biconcave capsule, whose membranes followed the Skalak constitutive (SK) law [23]. Internal and external fluids were modeled as an incompressible, Newtonian viscous fluid. The problem was characterized by Reynolds number and the capillary number $Ca$. The flow was assumed to be almost inertialess. The power (or energy dissipation) associated with membrane deformations was considered to quantify the state of membrane loads. Simulations were performed for different capillary numbers and viscosity ratios, as well as different initial positions of the RBC centroid.

2. Methods

2.1. Flow and RBC Model

We consider a cellular flow consisting of an external fluid (plasma), internal fluid (cytoplasm), and RBC with radius a in a circular channel of diameter D (2R), with a resolution of 16 fluid lattices per radius of RBC. This resolution was also applied in studies of channel flows [24,25] and the rheology of RBC suspensions [3]. The channel length is set to be 20a. The same computational length was also applied in the numerical analysis in [13]. In this study, we focused on the transition of RBC shape, especially from parachute shape to slipper shape [5,13], which usually occurs for $D\ge$ 12 $\mathsf{\mu }$m assuming a = 4 $\mathsf{\mu }$m [5,13,24]. Thus, we set the channel diameter to be D = 15 $\mathsf{\mu }$m. The RBC is modeled as a biconcave capsule, or a Newtonian fluid enclosed by a thin elastic membrane, with a major diameter d = 8 $\mathsf{\mu }$m (2a) and maximum thickness 2 $\mathsf{\mu }$m (a/2) [26]. The RBC is placed in a computational domain and is shown in Figure 1, where material points at the initial concave node point are represented by green dots, and those at the initial edge node point are indicated by blue dots.

The membrane is modeled as an isotropic and hyperelastic material following the SK law [23]. The strain energy w of the SK law is given by

(1)$w=\frac{G_s}{4}\left(I_1^2+2I_1-2I_2+C{I}_2^2\right),$

where $G_s$ is the surface shear elastic modulus; C is a coefficient representing the area incompressibility; $I_1$ (${\lambda }_1^2+{\lambda }_2^2-2$) and $I_2$ (${\lambda }_1^2{\lambda }_2^2-1=J_s^2-1$) are the first and second invariants of the Green–Lagrange strain tensor; ${\lambda }_i$ (i = 1 and 2) are the two principal in-plane stretch ratios; and $J_s={\lambda }_1{\lambda }_2$ is the Jacobian, which expresses the ratio of the deformed to reference surface areas. In this study, we set $C={10}^2$ [27]. Bending resistance was also considered [28], with a bending modulus $k_b=5.0\times {10}^{-19}$ J [29]. By mimicking a previous stretch experiment involving RBCs [30], the surface shear elastic modulus and the area incompressibility coefficient of RBCs at $Re$ = 0.2 were determined to be $G_s$ = 4.0 $\mathsf{\mu }$N/m and C = 10${}^2$, respectively [3,24]. These membrane parameters successfully reproduced the deformation of RBCs in shear flow [3,24] and also the thickness of the cell-depleted peripheral layer in circular channels [24]. We define the initial shape of RBCs as biconcave.

Neglecting inertial effects on membrane deformation, the static local equilibrium equation of the membrane is given by

(2)${\nabla }_s·T+q=0,$

where ${\nabla }_s(\left(I-nn\right)·\nabla )$ is the surface gradient operator, $n$ is the unit normal outward vector in the deformed state, and $T$ is the in-plane elastic tension that is obtained from the SK law (1).

It is known that the usual distribution of hemoglobin concentration in individual RBCs ranges from 27 to 37 g/dL corresponding to the internal fluid viscosity being ${\mu }_1$ = 5–15 cP (5–15 mPa·s) [31], while the normal plasma viscosity is ${\mu }_0$ = 1.1 − 1.3 cP (1.1 − 1.3 mPa·s) for plasma at 37 ℃ [32]. Hence, the physiologically relevant viscosity ratio can be taken as $\lambda ({\mu }_1/{\mu }_0)$ = 4.2 − 12.5 if the plasma viscosity is set to be ${\mu }_0$ = 1.2 cP. Hence, in our study, the physiological relevant viscosity ratio is set to be $\lambda$ = 5−10. The fluids are modeled as an incompressible Navier–Stokes equation, with a governing equation of fluid velocity $v$:

(3)$\begin{array}{cc}\hfill \rho \left(\frac{\partial v}{\partial t}+v·\nabla v\right)& =-\nabla p-\nabla ·\tau +\rho f,\hfill \end{array}$

(4)$\begin{array}{cc}\hfill \nabla ·v& =0,\hfill \end{array}$

where p is the pressure, $\rho$ is the fluid density, $f$ is the body force, and $\tau$ is the stress tensor of liquids and can be expressed by

(5)$\begin{array}{cc}\hfill \tau & =-\mu \left(\nabla v+\nabla v^T\right),\hfill \\ \hfill & =-\left\{\left(1-\alpha \right){\mu }_0+\alpha {\mu }_1\right\}\left(\nabla v+\nabla v^T\right),\hfill \end{array}$

where $\alpha$ is volume fraction of the inner fluid, which is in the range of 0 $\le \alpha \le$ 1. The dynamic condition requires that the load $q$ must be equal to the viscous traction jump across the mechanics:

(6)$\begin{array}{c}\hfill q=\left({\tau }_{out}-{\tau }_{in}\right)·n.\end{array}$

The problem is characterized by Reynolds number $Re$ and the capillary number $Ca$:

(7)$\begin{array}{cc}\hfill & Re=\frac{\rho D{V}_{max}^{\infty }}{{\mu }_0},\hfill \end{array}$

(8)$\begin{array}{cc}\hfill & Ca=\frac{{\mu }_0{\dot{\gamma }}_{m}a}{G_s}=\frac{{\mu }_0{V}_{max}^{\infty }}{G_s}\frac{a}{4R},\hfill \end{array}$

where $V_{max}^{\infty }$ (2$V_m^{\infty }$) is the maximum plasma velocity in the absence of any cells and ${\dot{\gamma }}_m$ ($V_m^{\infty }/D$) is the mean shear rate. Especially for the human microcirculation, flow can be assumed as inertialess, and is represented by $Re$ = 0.2 (corresponding to particle Reynolds number $R{e}_p(\rho {\dot{\gamma }}_m{a}^2/{\mu }_0=(a^2/\left(2D^2\right))Re)\approx$ 0.007) in this study. Although such a finite but low $Re$ accurately represents the capsule dynamics solved by the boundary integral method (BIM) in Stokes flow [3,24,33], we further tested a capsule deformation including large $Ca$ (≥1) (see also Appendix A). The condition defined by $Ca$ = 0.05 (and $Re$ = 0.2) corresponds to a typical venular wall-shear rate of 333 s${}^{-1}$ [34], and $Ca$ = 0.1 corresponds to an arteriolar wall-shear rate of 670 s${}^{-1}$ [35] if the surface shear elastic modulus is considered as $G_s$ = 4 $\mathsf{\mu }$N/m. Increasing $Re$ under constant $Ca$ corresponds to increasing $G_s$, namely, a harder RBC. Unless otherwise specified, we show the results obtained with $Re$ = 0.2 ($R{e}_p\approx$ 0.007).

2.2. Numerical Simulation

We used the D3Q19 LBM [36] coupled with the finite element method (FEM) [37]. Based on the virtual work principle, the above strong form (2) can be rewritten in weak form as

(9)${\int }_S\widehat{u}·qdS={\int }_S\widehat{ϵ}:TdS,$

where $\widehat{u}$ and $\widehat{ϵ}=\left({\nabla }_s\widehat{u}+{\nabla }_s{\widehat{u}}^T\right)$ are the virtual displacement and virtual strain, respectively. The FEM is used to solve Equation (9) and obtain the load $q$ acting on the membrane [37]. FEM and LBM were coupled by the immersed boundary method [38]. All procedures were fully implemented on a GPU to accelerate the numerical simulation. The flow was driven by a pressure gradient. Periodic boundary conditions were imposed on the inlet and outlet (z-direction). No-slip conditions were employed for the walls (x- and y-directions). The mesh size of the LBM for the fluid solution was set to be 250 nm, and that of the finite elements describing the membrane was approximately 250 nm (an unstructured mesh with 5120 elements was used for the FEM). This resolution has been shown to successfully represent single- and multicellular dynamics [24]; further, the results of multicellular dynamics are not changed by using twice the resolution for both the fluid and membrane meshes [24].

2.3. Analysis

To quantify the effect of the radial position of the RBC centroid and the deformed cell shape on fluid flow, the power (or energy dissipation) associated with membrane deformations is considered, and is given by

(10)$\begin{array}{cc}\hfill \delta W_{mem}& =\int \widehat{q}·\left(v^{\left(m\right)}-V^{\infty }\left(r\right)\right)dS,\hfill \\ \hfill \to \delta W_{mem}^{*}& =\frac{2}{{\mu }_{0}D{V}_{max}^{\infty 2}}\int {\widehat{q}}^{*}·\left(v^{\left(m\right)*}-V^{\infty *}\left(r\right)\right)d{S}^{*},\hfill \end{array}$

(11)$\begin{array}{cc}\hfill & =\frac{2}{{\mu }_{0}D{V}_{max}^{\infty 2}}\int \left[{\widehat{q}}_x^{*}{v}_x^{\left(m\right)*}+{\widehat{q}}_y^{*}{v}_y^{\left(m\right)*}+{\widehat{q}}_z^{*}\left(v_z^{\left(m\right)*}-V_z^{\infty *}\left(r\right)\right)\right]d{S}^{*},\hfill \end{array}$

where $V^{\infty }\left(r\right)=\left(0,0,V_{max}^{\infty }\left[1-{(r/R)}^2\right]\right)$ is the fluid flow velocity without cells, $\widehat{q}$ is the load acting on the membrane and includes the contribution of bending rigidity, r is the membrane distance from the channel center, $v^{\left(m\right)}$ is the interfacial velocity of the membrane, and S is the membrane surface area. Here, nondimensional variables are defined as ${\widehat{q}}^{*}=\widehat{q}/\left({\mu }_0{\dot{\gamma }}_m\right)$, $v^{\left(m\right)*}=v^{\left(m\right)}/V_{max}^{\infty }$, $V^{\infty *}=V^{\infty }/V_{max}^{\infty }$ and $S^{*}=S/D^2$.

3. Results

3.1. Effect of Capillary Number $Ca$ on RBC Shapes

First, we investigated the behavior of RBCs with a viscosity ratio $\lambda$ = 5 for different $Ca$ and different initial radial positions of the RBC centroid $r_0$. Figure 2a–d show examples of snapshots of flowing RBCs in a steady state (${\dot{\gamma }}_{m}t$ = 800, see also Videos S1–S4). Initially off-centered RBCs ($r_0/R$ = 0.4) subjected to the lowest $Ca$ (0.05) gradually migrated toward the channel center and exhibited a non-tank-treading (non-TT) discoidal shape, as shown in Figure 2a (see also Video S1). This occurred even when the initial position was set to be the channel center (Figure 2c, see also Video S3). Initially off-centered RBCs that were subjected to the highest $Ca$ (1.2) exhibited a TT slipper shape [5,13], as shown in Figure 2b (see also Video S2). Since initially centered RBCs exhibited a non-TT parachute shape (Figure 2d, see also Video S4), the stable shape of RBCs subjected to a higher $Ca$ depends on the initial position of the RBC centroid, which qualitatively agrees with previous experimental and numerical results in the rectangular microchannel [5]. Our numerical results further showed that a non-TT discoidal shape was observed for $Ca\le$ 0.1, which shifted to a non-TT parachute shape for 0.4 $\le Ca\le$ 0.8, regardless of initial position. Bistable shapes (non-TT parachute shape and TT slipper shape) were only observed for $Ca$ = 1.2. Note that at $Ca$ = 1.2, a non-TT parachute shape was observed for $r_0\le$ 1.5 $\mathsf{\mu }$m (i.e., $r_0/R\le$ 0.2), while a TT slipper shape was observed for $r_0\ge$ 2 $\mathsf{\mu }$m (i.e., $r_0/R\ge$ 0.27) (see also Figure A2 in Appendix B).

Figure 2e shows the time history of the radial position of each RBC centroid r, normalized by the channel radius R. Spherical capsules uniformly exhibited axial migration independent of $Ca$ (Figure A3 in Appendix C, see also Videos S5 and S6), which is consistent with classical principles regarding the axial migration of deformable spherical particles [8]; however, RBCs exhibited well-centered and off-centered migration depending on $Ca$ and the initial radial position. In this study, we hereafter define axial migration as having an order of magnitude $O(⟨r⟩/R)\le$ 10${}^{-2}$, and nonaxial migration as having an order of magnitude $O(⟨r⟩/R)>$ 10${}^{-2}$, where $⟨·⟩$ denotes the time average. The time-averaging was performed for data after ${\dot{\gamma }}_{m}t$ = 200 based on Figure 2e. The non-TT discoidal shape, obtained with $Ca$ = 0.05 (Figure 2a,c), flowed near the channel axis (blue and black lines in Figure 2e). The TT slipper shape, obtained with $Ca$ = 1.2 and $r_0/R$ = 0.4 (Figure 2b), exhibited nonaxial migration with temporal fluctuations (red line in Figure 2e). The non-TT parachute shape, obtained with $Ca$ = 1.2 and $r_0/R$ = 0 (Figure 2d), showed stable flow near the channel axis (orange line in Figure 2e). Note that the equilibrium radial position of each RBC was independent of its initial orientation, which was randomly determined (data is not shown).

Figure 3 shows the time average of the radial position of the RBC centroid as a function of $Ca$, where the error bars represent standard deviations of the time axis. RBCs exhibited axial or nonaxial migration depending on $Ca$ and the initial radial position. For relatively low $Ca\le$ 0.1, RBCs with a non-TT discoidal shape were located slightly away from the channel center regardless of initial position. For 0.4 $\le Ca\le$ 0.8, RBCs with a non-TT parachute shape were much closer to the channel axis. For larger $Ca$ (1.2), RBCs exhibited axial migration with a non-TT parachute shape and nonaxial migration with a TT slipper shape. Hence, the equilibrium radial position depends on the stable deformed RBC shape, which in turn, is related to $Ca$ and the initial position. Note that aforementioned stable deformed shapes (i.e., non-TT discoidal, non-TT parachute, and TT slipper shapes) remained consistent even at high $Re$ = 10, as did the equilibrium position (orange/black dots at $Ca$ = 0.05 and $Ca$ = 1.2).

Figure 4a shows the time average of the volumetric flow rate $⟨Q⟩$ with a RBC measured in simulations and normalized by the flow rate $Q^{\infty }$ ($\pi R^2{V}_m^{\infty }$) without a RBC as a function of $Ca$. The result also represents the change in flow resistance because the apparent viscosity is inversely proportional to the volumetric flow rate. Thus, the decrease of $⟨Q⟩$ means the increase of the apparent viscosity. $⟨Q⟩/Q^{\infty }$ was the largest with the TT slipper shape (red dot at $Ca$ = 1.2 in Figure 4a)—i.e., the flow resistance was the smallest—because the projected area of the deformed RBC to the cross-sectional area of the channel (x-y plane) $A_{xy}$ was the smallest, as shown in Figure 4b, where $A_{xy}$ is normalized by the initial projected RBC area with maximal length, i.e., $\pi a^2$. Since fluid drag force can be described as a proportion of the projected area, the TT slipper shape with small projected area $A_{xy}$ leads to a smaller flow resistance than the other two shapes (non-TT discoidal/parachute shape) (Figure 4a). A reduction in the flow resistance, accompanied with the transition from non-TT discoidal/parachute shapes to TT slipper shapes with the increases of $Ca$, is consistent with [13].

The powers associated with membrane deformations $⟨\delta W_{mem}^{*}⟩$ are shown in Figure 4c. The results of $⟨\delta W_{mem}^{*}⟩$ were not always correlated with those of $⟨Q⟩$ and $⟨A_{xy}⟩$. $⟨\delta W_{mem}^{*}⟩$ obtained with $r_0/R=0.4$ tended to increase for $Ca$ = 0.8, but abruptly decreased for further large $Ca$ (>0.8) (red dots in Figure 4c). $⟨\delta W_{mem}^{*}⟩$ obtained with $r_0/R=0$ tended to increase, at least for $Ca\le$ 1.2 (blue dots in Figure 4c). Since it was expected that $⟨\delta W_{mem}^{*}⟩$ would represent a more precise membrane load state than $⟨Q⟩$ or $⟨A_{xy}⟩$, we replotted $⟨\delta W_{mem}^{*}⟩$ as a function of $⟨r⟩/R$, as shown in Figure 4d. The result reflected the finding that $⟨\delta W_{mem}^{*}⟩$ was not affected by the initial position of $Ca$, but instead, was large for the near-center position with $⟨r⟩/R\le$ 10${}^{-2}$ and small for the off-centered position with $⟨r⟩/R>$ 10${}^{-2}$ (Figure 4d).

3.2. Effect of Viscosity Ratio $\lambda$ on RBC Shapes

Next, we investigated the effect of viscosity ratio $\lambda$ on the behavior of RBCs. Figure 5a shows snapshots of flowing RBCs subjected to low and high $Ca$ in steady state for different $\lambda$. The simulations were started from the off-centered position at $r_0/R$ = 0.4. RBCs subjected to low $Ca$ (0.05) uniformly exhibited the non-TT discoid shape regardless of $\lambda$ (bottom row in Figure 5a). In contrast, the stable shape of RBCs subjected to high $Ca$ (1.2) changed from non-TT parachute shape for $\lambda \le$ 2 to TT slipper shape for $\lambda \ge$ 5 (top row in Figure 5a).

Figure 5b shows the time average of the radial position of the RBC centroid $⟨r⟩/R$ as a function of $\lambda$. As described previously in Figure 3, non-TT parachute shapes, obtained with high $Ca$ (1.2) and $\lambda \le$ 2, approached the channel axis much more closely than the non-TT discoidal shape obtained with low $Ca$ (0.05). The TT slipper shape, which was only found for high $Ca$ (1.2) and $\lambda \ge$ 5, clearly exhibited nonaxial migration (i.e., $O(⟨r⟩/R)>{10}^{-2}$). Furthermore, comparing the results of $⟨r⟩/R$ between $\lambda$ = 5 and $\lambda$ = 10, a higher viscosity ratio allowed RBCs to be positioned away from the channel center (Figure 5b).

Figure 6a shows the time average of the volumetric flow rate $⟨Q⟩$ as a function of $\lambda$. The non-TT discoidal shape, obtained with $Ca$ = 0.05, was associated with a near-center position (Figure 5b) and a relatively large projected area $A_{xy}$ (Figure 6b), resulting in large flow resistance, i.e., small $⟨Q⟩/Q^{\infty }$ (Figure 6a). Although the non-TT parachute shape obtained with $Ca$ = 1.2 and $\lambda \le$ 2 was also associated with a large $A_{xy}$, $⟨Q⟩/Q^{\infty }$ was larger than with the non-TT discoidal shape (Figure 6a). The values of $A_{xy}$ for the TT slipper shape, which was obtained with $Ca$ = 1.2 and $\lambda \ge$ 5, were smaller than for the other two shapes (Figure 6b). Further, with the TT slipper shape, the value of $⟨Q⟩/Q^{\infty }$ was relatively large except for $\lambda$ = 10 (Figure 6a). These results suggest that the projected area $A_{xy}$ is not always correlated with flow resistance represented by $⟨Q⟩/Q^{\infty }$, and hence, a more precise description of membrane dynamics is necessary to better understand the relationship between stable shapes and equilibrium radial positions.

Figure 6c shows the powers associated with membrane deformations $⟨\delta W_{mem}^{*}⟩$ as a function of $\lambda$. The results of $⟨\delta W_{mem}^{*}⟩$ at $Ca$ = 0.05 consistently had a relatively small order of magnitude, while those at $Ca$ = 1.2 started to decrease as $\lambda$ increased to be >1 (Figure 6c). These trends are consistent with those in the equilibrium position $⟨r⟩/R$ described in Figure 5b. Although the results of $⟨\delta W_{mem}^{*}⟩$ did not correlate with those of $⟨Q⟩/Q^{\infty }$, as previously described in Figure 4a,c, they did correlate well with the equilibrium position $⟨r⟩/R$ regardless of $\lambda$ and $Ca$, as shown in Figure 6d.

4. Discussion

Various numerical models have been used to systematically investigate the behavior of a single RBC at low $Re$ in microchannels whose scale is comparable to the cell size [5,12,13]. For instance, Fedosov et al., (2014) performed simulations of RBC behavior for a fixed viscosity ratio $\lambda$ (1) in circular microchannels and showed a phase diagram of stable RBCs as a function of shear rates and size ratios in the range of 0.3 $<d/D<$ 0.8 [13]. Guckenberger et al., (2018) also performed simulations of RBCs for a viscosity ratio $\lambda$ = 5 in a rectangular microchannel of width W = 12 $\mathsf{\mu }$m and height H = 10 $\mathsf{\mu }$m (i.e., $d/H$ = 0.8 and $d/W$ = 0.67), and showed a phase diagram of stable RBCs as a function of cell velocity and initial position [5]. However, it has not yet been fully determined how these shapes at equilibrium position contribute to fluid flow. In this study, we further investigated energy expenditure due to membrane deformation of RBCs in a circular microchannel with a diameter of D = 15 $\mathsf{\mu }$m, i.e., $d/D$ = 0.53, and found that it correlated well with the equilibrium position of the RBC regardless of $Ca$, $\lambda$, and the initial position of the RBC centroid (Figure 4d and Figure 6d). The results are summarized in Figure 7, which consists of replotted data from Figure 4d and Figure 6d. The results suggest that the equilibrium radial position of the RBC centroid is determined by the stable deformed shape, due to different energy expenditures associated with various membrane deformations. The non-TT parachute shape allows the RBC to approach the channel axis with an order of magnitude $O(⟨r⟩/R)\sim {10}^{-2}$, while the radial position of the non-TT discoidal shape is shifted slightly away from the channel center with ${10}^{-2}<⟨r⟩/R<{10}^{-1}$ (Figure 7). The TT slipper shape always demonstrates nonaxial migration with $⟨r⟩/R>{10}^{-1}$ (Figure 7). This shape was only observed in limited conditions: $Ca$ = 1.2, $\lambda \ge$ 5, and an initial off-centered position $r_0/R$ = 0.4. Such high $Ca$ (1.2) corresponds to a higher wall shear rate in a circular channel ${\dot{\gamma }}_{wall}\left(8{\dot{\gamma }}_m\right)$ = 8 × 10${}^3$ s${}^{-1}$, which is over 10 times greater than human arterial wall shear rates [35]. Therefore, it is expected that the TT slipper shape with an off-centered position will be found in in vitro systems with artificially high shear rates except for pathological vascular regions, e.g., arterial stenosis [39]; hence, this shape may be a hallmark in cell sorting techniques using microfluidics. Furthermore, at high $Ca$ (1.2), the parachute shape was found instead of the TT slipper shape for $\lambda \le$ 2 (Figure 5a). Thus, the parachute shape may also be useful as an indicator for identifying cytoplasmic viscosity. Since we set the channel diameter to be 15 $\mathsf{\mu }$m, it would be interesting to study how the off-centered position $⟨r⟩/R$ depends on stable deformed shape changes in larger microtubes.

The motion of a TT slipper shape in this study may be redefined as a snaking motion (periodic oscillation of the shape in the form of a snake motion), according to previous numerical analyses using a 2D vesicle model [40,41,42]. However, in this study, we focused on stable deformed RBC shapes, and did not rigorously differentiate between TT slipper shape and snaking motion, following from previous numerical studies [5,43]. Although we have shown that our numerical models successfully reproduce not only the deformation of single RBC but also the thickness of the cell-depleted peripheral layer due to multicellular interactions [24], comparisons in flowing RBC shapes between numerical results and experimental observations in microtubes have not yet been conducted, which is our future study. As in many numerical studies, e.g., [5,43], we also neglected membrane viscosity, which has been introduced only in a few continuum model analyses [44,45,46]. Hence, it would be also interesting to study how such fluid deformable surfaces changes the stable deformed RBC shape relative to a pure elastic membrane. Although we performed simulations for a wide range of $Ca$ and $\lambda$, we are unsure what factors cause RBCs to adopt a stable shape under specific equilibrium positions. Considering the finding that flow resistance characterized by the volumetric flow rate $⟨Q⟩/Q^{\infty }$ is not always described by the projected area of the RBC to the cross-sectional area of the channel $A_{xy}$ (Figure 4b and Figure 6b), membrane dynamics need to be more precisely investigated to clarify this problem. In the future, we will report the precise mechanical characteristics of the stable deformed RBC shapes, as well as the relationship between these stable shapes, the equilibrium radial positions of RBCs, and the membrane load.

5. Conclusions

We numerically investigated the dynamics of translating RBCs in a circular microchannel with a diameter of 15 $\mathsf{\mu }$m for different capillary numbers $Ca$ and viscosity ratios $\lambda$. The flow was assumed to be almost inertialess. Our results demonstrated that the presence of axial or nonaxial migration depends on the stable deformed RBC shapes, and that the equilibrium radial position of the RBC centroid correlated well with the energy expenditure associated with different membrane deformations. The non-TT parachute shape, obtained with high $Ca$ and low $\lambda$, allowed RBCs to approach the channel axis, while the non-TT discoidal shape, obtained with low $Ca$, shifted the radial position of RBCs slightly away from the channel center. The TT slipper shape, obtained with high $Ca$ and high $\lambda$, was always accompanied by obvious nonaxial migration. The energy expenditure decreased in the following order: Non-TT parachute shape, non-TT discoidal shape, and TT slipper shape. In the near future, we will examine the shape stability of deformed RBCs in more details to clarify precise mechanical characteristics of the stable shapes, and report the relationship between these stable shapes, the equilibrium radial position of RBCs, and the membrane load.

Author Contributions

N.T., analyzed data; N.T. performed experiments; N.T., H.Y., T.O., N.Y. and M.S.-S. interpreted simulation results; N.T. prepared figures; N.T. drafted the manuscript; N.T. edited and revised the manuscript; N.T., H.Y., T.O., N.Y. and M.S.-S. approved the final version of the manuscript; N.T. and T.O. contributed to the research conception and design. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by JSPS KAKENHI Grant Number JP20H04504, JP20H02072, and by the Keihanshin Consortium for Fostering the Next Generation of Global Leaders in Research (K-CONNEX), established by Human Resource Development Program for Science and Technology. N.T. is grateful for the financial support of UCL-Osaka Partner Funds.

Acknowledgments

N.T. thanks Hiroaki Ito, and Marco E. Rosti for their helpful discussions.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

RBC	Red blood cell
LBM	Lattice-Boltzmann method
FEM	Finite element method
IBM	Immersed boundary method
GPU	Graphics processing unit
(non-)TT motion	(non-)tank-treading motion

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Simulation setup: A single RBC with radius of 4 [Forumla omitted. See PDF.]m is placed with random orientation in a circular channel with diameter of 15 [Forumla omitted. See PDF.]m and length of 80 [Forumla omitted. See PDF.]m. Periodic boundary conditions are imposed on the flow (z-direction) and no-slip conditions are employed for the wall (x- and y-direction). Green dots represent material points at the initial concave node point, and blue dots at the initial edge node point.

Figure 2. Snapshots of flowing RBCs in steady state ([Forumla omitted. See PDF.] = 800) for (a,c) [Forumla omitted. See PDF.] = 0.05 (see also Video S1 and S3) and (b,d) [Forumla omitted. See PDF.] = 1.2 (see also Video S2 and S4). The right side is the axial view and the left side is the lateral view; the flow direction is from left to right. For each [Forumla omitted. See PDF.], the initial position of the RBC centroid is set to be (c,d) [Forumla omitted. See PDF.] = 0 and (a,b) [Forumla omitted. See PDF.] = 0.4. (e) Time history of the radial position of the RBC centroid [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The results were obtained with [Forumla omitted. See PDF.] = 5.

Figure 3. Time average of the radial position of the RBC centroid [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for initial position [Forumla omitted. See PDF.] = 0 (triangles) and [Forumla omitted. See PDF.] = 0.4 (inverse triangles), where [Forumla omitted. See PDF.] denotes time average. The error bars represent standard deviations on the time axis. The results at [Forumla omitted. See PDF.] = 10 for low [Forumla omitted. See PDF.] (0.05) and high [Forumla omitted. See PDF.] (1.2) are also plotted, with black squares for [Forumla omitted. See PDF.] = 0.4 and orange circles for [Forumla omitted. See PDF.] = 0. The results were obtained with [Forumla omitted. See PDF.] = 5.

Figure 4. (a) Time average of the volumetric flow rate [Forumla omitted. See PDF.], (b) projected area of the RBC to the cross-sectional area of the channel (x-y plane) [Forumla omitted. See PDF.], and (c) powers associated with membrane deformations [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for different initial positions [Forumla omitted. See PDF.] (0 and 0.4). (d) Replotted data of [Forumla omitted. See PDF.] as a function of equilibrium radial position [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.]. The error bars represent standard deviations on the time axis. The results were obtained with [Forumla omitted. See PDF.] = 5.

Figure 5. (a) Snapshots of flowing RBCs in steady state at [Forumla omitted. See PDF.] = 0.05 (bottom row) and [Forumla omitted. See PDF.] = 1.2 (top row) for different [Forumla omitted. See PDF.]. (b) Time average of the radial position of the RBC centroid [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] = 0.05 (blue inverse triangles) and [Forumla omitted. See PDF.] = 1.2 (red triangles). The results were obtained with an initial off-centered position at [Forumla omitted. See PDF.] = 0.4.

Figure 6. (a) Time average of the volumetric flow rate [Forumla omitted. See PDF.], (b) projected area of the RBC [Forumla omitted. See PDF.], and (c) powers associated with membrane deformations [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.] (0.05 and 1.2). (d) Replotted data of [Forumla omitted. See PDF.] as a function of equilibrium radial position [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.]. The results were obtained with an initial off-centered position [Forumla omitted. See PDF.] = 0.4.

Figure 7. Powers associated with membrane deformations [Forumla omitted. See PDF.] as a function of equilibrium radial position [Forumla omitted. See PDF.] for different shapes. The results are replotted from Figure 4d and Figure 6d.

Supplementary Materials

The following materials are available online at https://www.mdpi.com/article/10.3390/mi12101162/s1. Video S1: Non-TT discoidal shape, obtained with $$Ca$$ = 0.05 and $$r_0/R$$ = 0.4. Video S2: TT slipper shape, obtained with $$Ca$$ = 1.2 and $$r_0/R$$ = 0.4. Video S3: Non-TT discoidal shape, obtained with $$Ca$$ = 0.05 and $$r_0/R$$ = 0. Video S4: Non-TT parachute shape, obtained with $$Ca$$ = 1.2 and $$r_0/R$$ = 0. Videos S5 and S6: Numerical results of a spherical capsule initially placed at $$r_0/R$$ = 0.4 for $$Ca$$ = 0.05 and $$Ca$$ = 1.2, respectively. These numerical results were obtained with $$Re$$ = 0.2 and $$\lambda$$ = 5.

References

1. Lázaro, G.R.; Hernández-Machadoa, A.; Pagonabarraga, I. Rheology of red blood cells under flow in highly confined microchannels. I. Effect of elasticity. Soft Matter; 2014; 10, pp. 7195-7206. [DOI: https://dx.doi.org/10.1039/C4SM00894D]

2. Lázaro, G.R.; Hernández-Machadoa, A.; Pagonabarraga, I. Rheology of red blood cells under flow in highly confined microchannels. II. Effect of focusing and confinement. Soft Matter; 2014; 10, pp. 7207-7217. [DOI: https://dx.doi.org/10.1039/C4SM01382D] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25068313]

3. Takeishi, N.; Rosti, M.E.; Imai, Y.; Wada, S.; Brandt, L. Haemorheology in dilute, semi-dilute and dense suspensions of red blood cells. J. Fluid Mech.; 2019; 872, pp. 818-848. [DOI: https://dx.doi.org/10.1017/jfm.2019.393]

4. Skalak, R.; Branemark, P.I. Deformation of red blood cells in capillaries. Science; 1969; 164, pp. 717-719. [DOI: https://dx.doi.org/10.1126/science.164.3880.717] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/5778020]

5. Guckenberger, A.; Kihm, A.; John, T.; Wagner, C.; Gekle, S. Numerical-experimental observation of shape bistability of red blood cells flowing in a microchannel. Soft Matter; 2018; 14, pp. 2032-2043. [DOI: https://dx.doi.org/10.1039/C7SM02272G] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29473072]

6. Yaya, F.; Römer, J.; Guckenberger, A.; John, T.; Gekle, S.; Podgorski, T.; Wagner, C. Vortical flow structures induced by red blood cells in capillaries. Microcirculation; 2021; 28, e12693. [DOI: https://dx.doi.org/10.1111/micc.12693]

7. Takeishi, N.; Ito, H.; Kaneko, M.; Wada, S. Deformation of a red blood cell in a narrow rectangular microchannel. Micromachines; 2019; 10, 199. [DOI: https://dx.doi.org/10.3390/mi10030199] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30901883]

8. Karnis, A.; Goldsmith, H.L.; Mason, S.G. Axial migration of particles in Poiseuille flow. Nature; 1963; 14, pp. 284-304. [DOI: https://dx.doi.org/10.1038/200159a0]

9. Shi, L.; Pan, T.-W.; Glowinski, R. Lateral migration and equilibrium shape and position of a single red blood cell in bounded Poiseuille flows. Phys. Rev. E; 2012; 86, 056306. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.056308]

10. Kaoui, B.; Biros, G.; Misbah, C. Why do red blood cells have asymmetric shapes even in a symmetric flow?. Phys. Rev. Lett.; 2009; 103, 188101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.103.188101]

11. Hogan, B.; Shen, Z.; Zhang, H.; Misbah, C.; Barakat, A.I. Shear stress in the microvasculature: Influence of red blood cell morphology and endothelial wall undulation. Biomech. Model. Mechanobiol.; 2019; 18, pp. 1095-1109. [DOI: https://dx.doi.org/10.1007/s10237-019-01130-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30840162]

12. Noguchi, H.; Gompper, G. Shape transitions of fluid vesicles and red blood cells in capillary flows. Proc. Natl. Acad. Sci. USA; 2005; 102, pp. 14159-14164. [DOI: https://dx.doi.org/10.1073/pnas.0504243102]

13. Fedosov, D.A.; Peltomäki, M.; Gompper, G. Deformation and dynamics of red blood cells in flow through cylindrical microchannels. Soft Matter; 2014; 10, pp. 4258-4267. [DOI: https://dx.doi.org/10.1039/C4SM00248B] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24752231]

14. Ciftlik, A.T.; Ettori, M.; Gijs, M.A.M. High throughput-per-footprint inertial focusing. Small; 2013; 9, pp. 2764-2773. [DOI: https://dx.doi.org/10.1002/smll.201201770] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23420756]

15. Fregin, B.; Czerwinski, F.; Biedenweg, D.; Girardo, S.; Gross, S.; Aurich, K.; Otto, O. High-throughput single-cell rheology in complex samples by dynamic real-time deformability cytometry. Nat. Commun.; 2019; 10, 415. [DOI: https://dx.doi.org/10.1038/s41467-019-08370-3]

16. Ito, H.; Murakami, R.; Sakuma, S.; Tsai, C.-H.D.; Gutsmann, T.; Brandenburg, K.; Poöschl, J.M.B.; Arai, F.; Kaneko, M.; Tanaka, M. Mechanical diagnosis of human eryhrocytes by ultra-high speed manipulation unraveled critical time window for global cytoskeletal remodeling. Sci. Rep.; 2017; 7, 43134. [DOI: https://dx.doi.org/10.1038/srep43134] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28233788]

17. Kihm, A.; Kaestner, L.; Wagner1, C.; Quint, S. Classification of red blood cell shapes in flow using outlier tolerant machine learning. PLoS Comput. Biol.; 2019; 14, e1006278. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1006278]

18. Lu, X.; Wood, D.K.; Higgins, J.M. Deoxygenation reduces sickle cell blood flow at arterial oxygen tension. Biophys. J.; 2016; 110, pp. 2751-2758. [DOI: https://dx.doi.org/10.1016/j.bpj.2016.04.050] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27332133]

19. McMahon, T.J. Red blood cell deformability, vasoactive mediators, and adhesion. Front. Physiol.; 2019; 10, 1417. [DOI: https://dx.doi.org/10.3389/fphys.2019.01417] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31803068]

20. Chien, S.; Usami, S.; Bertles, J.F. Abnormal rheology of oxygenated blood in sickle cell anemia. J. Clin. Investig.; 1970; 49, pp. 623-634. [DOI: https://dx.doi.org/10.1172/JCI106273]

21. Usami, S.; Chien, S.; Scholtz, P.M.; Bertles, J.F. Effect of deoxygenation on blood rheology in sickle cell disease. Microvasc. Res.; 1975; 9, pp. 324-334. [DOI: https://dx.doi.org/10.1016/0026-2862(75)90069-2]

22. Kaul, D.K.; Xue, H. Rate of deoxygenation and rheologic behavior of blood in sickle cell anemia. Blood; 1991; 77, pp. 1353-1361. [DOI: https://dx.doi.org/10.1182/blood.V77.6.1353.1353] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/2001458]

23. Skalak, R.; Tozeren, A.; Zarda, R.P.; Chien, S. Strain energy function of red blood cell membranes. Biophys. J.; 1973; 13, pp. 245-264. [DOI: https://dx.doi.org/10.1016/S0006-3495(73)85983-1]

24. Takeishi, N.; Imai, Y.; Nakaaki, K.; Yamaguchi, T.; Ishikawa, T. Leukocyte margination at arteriole shear rate. Physiol. Rep.; 2014; 2, e12037. [DOI: https://dx.doi.org/10.14814/phy2.12037] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24907300]

25. Takeishi, N.; Imai, Y.; Yamaguchi, T.; Ishikawa, T. Flow of a circulating tumor cell and red blood cells in microvessels. Phys. Rev. E; 2015; 92, 063011. [DOI: https://dx.doi.org/10.1103/PhysRevE.92.063011] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26764808]

26. Evans, E.; Fung, Y.-C. Improved measurements of the erythrocyte geometry. Microvasc. Res.; 1972; 4, pp. 335-347. [DOI: https://dx.doi.org/10.1016/0026-2862(72)90069-6]

27. Barthés-Biesel, D.; Diaz, A.; Dheni, E. Effect of constitutive laws for two-dimensional membranes on flow-induced capsule deformation. J. Fluid. Mech.; 2002; 460, pp. 211-222. [DOI: https://dx.doi.org/10.1017/S0022112002008352]

28. Li, J.; Dao, M.; Lim, C.T.; Suresh, S. Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte. Phys. Fluid; 2005; 88, pp. 3707-6719. [DOI: https://dx.doi.org/10.1529/biophysj.104.047332] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15749778]

29. de Morales Marinkovic, M.P.; Turner, K.T.; Butler, J.P.; Fredberg, J.J.; Suresh, S. Viscoelasticity of the human red blood cell. Am. J. Physiol. Cell Physiol.; 2007; 293, pp. C597-C605. [DOI: https://dx.doi.org/10.1152/ajpcell.00562.2006]

30. Suresh, S.; Spatz, J.; Mills, J.P.; Micoulet, A.; Dao, M.; Lim, C.T.; Beil, M.; Seufferlein, T. Connections between single-cell biomechanics and human disease states: Gastrointestinal cancer and malaria. Acta Biomater.; 2005; 1, pp. 15-30. [DOI: https://dx.doi.org/10.1016/j.actbio.2004.09.001]

31. Mohandas, N.; Gallagher, P.G. Red cell membrane: Past, present, and future. Blood; 2008; 112, pp. 3939-3948. [DOI: https://dx.doi.org/10.1182/blood-2008-07-161166] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18988878]

32. Harkness, J.; Whittington, R.B. Blood-plasma viscosity: An approximate temperature-invariant arising from generalised concepts. Biorheology; 1970; 6, pp. 169-187. [DOI: https://dx.doi.org/10.3233/BIR-1970-6304] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/5418877]

33. Takeishi, N.; Imai, Y.; Ishida, S.; Omori, T.; Kamm, R.D.; Ishikawa, T. Cell adhesion during bullet motion in capillaries. Am. J. Physiol. Heart Circ. Physiol.; 2016; 311, pp. H395-H403. [DOI: https://dx.doi.org/10.1152/ajpheart.00241.2016]

34. Koutsiaris, A.G.; Tachmitzi, S.V.; Batis, N. Wall shear stress quantification in the human conjunctival pre-capillary arterioles in vivo. Microvasc. Res.; 2013; 85, pp. 34-39. [DOI: https://dx.doi.org/10.1016/j.mvr.2012.11.003]

35. Koutsiaris, A.G.; Tachmitzi, S.V.; Batis, N.; Kotoula, M.G.; Karabatsas, C.H.; Tsironi, E.; Chatzoulis, D.Z. Volume flow and wall shear stress quantification in the human conjunctival capillaries and post-capillary venules in vivo. Biorheology; 2007; 44, pp. 375-386.

36. Chen, S.; Doolen, G.D. Lattice boltzmann method for fluid flow. Annu. Rev. Fluid Mech.; 1998; 30, pp. 329-364. [DOI: https://dx.doi.org/10.1146/annurev.fluid.30.1.329]

37. Walter, J.; Salsac, A.V.; Barthés-Biesel, D.; Le Tallec, P. Coupling of finite element and boundary integral methods for a capsule in a stokes flow. Int. J. Numer. Meth. Eng.; 2010; 83, pp. 829-850. [DOI: https://dx.doi.org/10.1002/nme.2859]

38. Peskin, C.S. The immersed boundary method. Acta Numer.; 2002; 11, pp. 479-517.

39. Case, L.D.C.; Ku, D.N. Thrombus formation at high shear rates. Annu. Rev. Biomed. Eng.; 2017; 19, pp. 413-415. [DOI: https://dx.doi.org/10.1146/annurev-bioeng-071516-044539] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28441034]

40. Aouane, O.; Thiébaud, M.; Benyoussef, A.; Wagner, C.; Misbah, C. Vesicle dynamics in a confined Poiseuille flow: From steady state to chaos. Phys. Rev. E; 2014; 90, 033011. [DOI: https://dx.doi.org/10.1103/PhysRevE.90.033011] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25314533]

41. Kaoui, B.; Tahiri, N.; Biben, T.; Ez-Zahraouy, H.; Benyoussef, A.; Biros, G.; Misbah, C. Complexity of vesicle microcirculation Phys. Rev. E; 2011; 84, 041906.

42. Tahiri, N.; Biben, T.; Ez-Zahraouy, H.; Benyoussef, A.; Misbah, C. On the problem of slipper shapes of red blood cells in the microvasculature. Microvasc. Res.; 2013; 85, pp. 40-45. [DOI: https://dx.doi.org/10.1016/j.mvr.2012.10.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23063869]

43. Ye, T.; Shi, H.; Peng, L.; Li, Y. Numerical studies of a red blood cell in rectangular microchannels. J. Appl. Phys.; 2017; 122, 084701. [DOI: https://dx.doi.org/10.1063/1.5000357]

44. Peng, Z.; Asaro, R.J.; Zhu, Q. Multiscale modelling of erythrocytes in Stokes flow. J. Fluid Mech.; 2011; 686, pp. 299-337. [DOI: https://dx.doi.org/10.1017/jfm.2011.332]

45. Torres-Sánchez, A.; Millán, D.; Arroyo, M. Modelling fluid deformable surfaces with an emphasis on biological interfaces. J. Fluid Mech.; 2019; 872, pp. 271-281. [DOI: https://dx.doi.org/10.1017/jfm.2019.341]

46. Yazdani, A.; Bagchi, P. Influence of membrane viscosity on capsule dynamics in shear flow. J. Fluid Mech.; 2013; 718, pp. 569-595. [DOI: https://dx.doi.org/10.1017/jfm.2012.637]

47. 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, pp. 117-143. [DOI: https://dx.doi.org/10.1017/S0022112098008714]

48. Foessel, É.; Walter, J.; Salsac, A.-V.; Barthés-Biesel, D. Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow. J. Fluid Mech.; 2011; 672, pp. 477-486. [DOI: https://dx.doi.org/10.1017/S0022112011000280]