Hostname: page-component-6898cfcb55-w5hgh Total loading time: 0 Render date: 2026-11-11T05:08:20.247Z Has data issue: false hasContentIssue false
  • English
  • Français

Stability of the tank treading modes of erythrocytes and its dependence on cytoskeleton reference states

Published online by Cambridge University Press:  20 April 2015

Zhangli Peng,
Sara Salehyar and
Qiang Zhu
Zhangli Peng
Affiliation:
Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, IN 46556, USA
Sara Salehyar
Affiliation:
Department of Structural Engineering, University of California San Diego, La Jolla, CA 92093, USA
Qiang Zhu*
Affiliation:
Department of Structural Engineering, University of California San Diego, La Jolla, CA 92093, USA
*
†Email address for correspondence: qizhu@ucsd.edu
Rights & Permissions [Opens in a new window]

Abstract

We studied the tank treading motion of an erythrocyte (red blood cell, or RBC) in linear shear flows by using a boundary-element fluid-dynamics model coupled with a multiscale structural model of the cell. The purpose was to investigate the correlation between the reference (stress-free) state of the cytoskeleton and the cell dynamics in shear flows with relatively high capillary numbers. We discovered that there exist two distinctive modes of tank treading, mode 1 and mode 2. In mode 1 the membrane elements originating from the dimple areas keep close to the central plane, whereas in mode 2 these elements remain near the farthermost locations from the central plane. Mode 1 is also characterized by significantly higher breathing and swinging oscillations. During tank treading one mode may become unstable and switch to the other. Their stability depends on the viscosity ratio and the capillary number. At a fixed viscosity ratio, when the capillary number is increased the cell experiences sequentially a region dominated by mode 2, a mode 1/mode 2 bistable region and a region dominated by mode 1. More profoundly, these regions are highly sensitive to the reference state of the cytoskeleton. For example, compared with a cell with a biconcave reference state, a cell with a spheroidal reference state features a much smaller region dominated by mode 2. This finding may guide experiments to identify the actual reference state of these cells.

JFM classification

Biological Fluid Dynamics: Biological Fluid DynamicsBiological Fluid Dynamics: Capsule/cell dynamics

Information

Type
Papers
Information
Journal of Fluid Mechanics , Volume 771 , 25 May 2015, pp. 449 - 467
Check for updates
Copyright
© 2015 Cambridge University Press 

1. Introduction

With one of the simplest molecular structures among all cells, the red blood cell (RBC) provides an ideal model system to study the mechanical properties and the mechanics-versus-function relation of cells. A mature RBC has no nucleus – it is essentially a droplet of cytosol enclosed within a highly flexible membrane. The membrane itself features a composite structure, consisting of a lipid bilayer on the outside of the cell and a protein skeleton underneath (Mohandas & Evans Reference Mohandas and Evans1994). The skeleton contains approximately 33 000 repeating units called junctional complexes (JCs). Each JC has a central piece of actin protofilament surrounded by (up to) six spectrin dimers, forming a hexagon shape. The skeleton is pinned to (not embedded in) the lipid bilayer. Moreover, the pinning points made of transmembrane proteins are mobile so that the skeleton and the bilayer can slide against each other. Subsequently, these two components of the cell membrane may undergo different local deformations. This has been confirmed via fluorescent measurements of the variation of the skeleton density during micropipette aspirations (Discher, Mohandas & Evans Reference Discher, Mohandas and Evans1994). However, this fact has not been represented in most existing models, which describe the membrane as a single layer of uniform material.

In the cardiovascular system RBCs often travel in non-uniform flows. In microfluidic experiments, theoretical models and numerical simulations these non-uniform flows are usually simplified as steady linear shear flows (see, for example, Fischer, Stohr-Liesen & Schmid-Schonbein Reference Fischer, Stohr-Liesen and Schmid-Schonbein1978; Barthès-Biesel Reference Barthès-Biesel1980; Keller & Skalak Reference Keller and Skalak1982; Eggleton & Popel Reference Eggleton and Popel1998; Abkarian, Faivre & Viallat Reference Abkarian, Faivre and Viallat2007; Noguchi Reference Noguchi2009a ,Reference Noguchi b ; Fedosov, Caswell & Karniadakis Reference Fedosov, Caswell and Karniadakis2011; Vlahovska et al. Reference Vlahovska, Young, Danker and Misbah2011; Yazdani, Kalluri & Bagchi Reference Yazdani, Kalluri and Bagchi2011; Fischer & Korzeniewski Reference Fischer and Korzeniewski2013). It has been illustrated that depending on the shear rate the RBC may display one of the following motions: (a) tumbling, in which the cell somersaults in a rigid-body-like rotation around an axis perpendicular to the shear plane; (b) rolling (Bitbol Reference Bitbol1986), in which the cell rotates like a wheel around its shortest axis which is perpendicular to the shear plane; (c) frisbee-like motion (Dupire, Socol & Viallat Reference Dupire, Socol and Viallat2012), which is similar to the rolling motion except that the axis of rotation is within the shear plane; and (d) tank treading, in which the cell itself remains almost steady (except for oscillations in its orientation and aspect ratio), while its membrane circulates around it. Omori et al. (Reference Omori, Imai, Yamaguchi and Ishikawa2012a ) studied the reorientation of a capsule in shear flow with different capillary numbers $Ca$ . They found that a material point on the capsule moves towards the out-of-shear-plane axis $z$ when $Ca=1.0$ and moves towards $z=0$ when $Ca=0.3$ (see figure 2 in that paper). More interestingly, they found that the final orientation is independent of the initial orientation.

Our recent study indicated that there are in fact two sub-modes in the tank treading motion (Peng & Zhu Reference Peng and Zhu2013). These two modes are distinguished mostly by the trajectories of the membrane elements originating from the dimple areas (hereafter referred to as the dimple elements) during tank treading. To elaborate, in mode 1 these dimple elements remain close to the central plane, whereas in mode 2 these elements are away from the central plane. In addition, mode 1 is also characterized by significant breathing (with changes in aspect ratio) and swinging (with changes in pitching angle) oscillations. Both are negligibly small in mode 2. In that study, however, no systematic simulations were conducted to examine the stability of these tank treading modes and the transition between them. Moreover, all the results were obtained under the assumption that the skeleton is stress-free at its natural biconcave state. The validity of this assumption, however, has been questioned in some recent publications.

Indeed, the prestress inside the protein skeleton at the natural state is considered to be one of the most intriguing mysteries in RBC mechanics (Hoffman Reference Hoffman2001). Since the stress can be decomposed into isotropic stress (pressure) and shear stress, there are actually two questions, both of which remain unanswered. With regard to the isotropic stress, Svoboda et al. (Reference Svoboda, Schmidt, Branton and Block1992) found that the cytoskeleton was prestretched in the natural biconcave state, whereas Discher, Boal & Boey (Reference Discher, Boal and Boey1998) found that it was compressed. The effect of isotropic stress upon the cell mechanics has also been numerically investigated by Omori et al. (Reference Omori, Ishikawa, Barthès-Biesel, Salsac, Imai and Yamaguchi2012b ). The question about the shear stress, on the other hand, is usually expressed as the search for the shear-stress-free state of the skeleton (see, for example, Lim, Wortis & Mukhopadhyay Reference Lim, Wortis and Mukhopadhyay2002). Hereby the stress-free state is also called the reference state. It refers to a three-dimensional geometry of the skeleton in which there is no prestress. The overall surface area of this geometry is the same as that of the cell in its resting state (due to the area conservation imposed by the large area modulus of the lipid bilayer), whereas the internal volume is different. The resting shape of the RBC and its stress-free state can be connected through a deflation process during which the internal volume of the stress-free state is changed gradually until the actual surface area to volume ratio of a cell is reached (Peng, Mashayekh & Zhu Reference Peng, Mashayekh and Zhu2014). To date, most models assume that the reference state is biconcave so that there is no pretension in the skeleton. This is convenient in explaining the resting shape of the cell (which coincides with the minimum bending energy in the lipid bilayer) (see, for example, Mohandas & Evans Reference Mohandas and Evans1994). It is also consistent with the reported shape memory, which suggests that the reference state must be anisotropic so that it cannot be a sphere (Fischer Reference Fischer2004). This reference state, however, creates significant difference in strain energy between dimple and rim elements so that there is a high-energy barrier between these two areas. A direct consequence of this energy barrier is that it discourages tank treading motion since it involves location exchanges between the dimple and rim elements.

Based on this property, Dupire et al. (Reference Dupire, Socol and Viallat2012) conducted a sequence of microfluidic experiments in low-shear-rate flows and documented the transitions from tumbling and rolling to tank treading motions. The results indicate that the critical shear rates between tumbling/rolling and tank treading are too low to be explained by the biconcave reference state. In addition, in these tank treading motions in low-shear-rate flows, the cell keeps its biconcave shape (Abkarian et al. Reference Abkarian, Faivre and Viallat2007; Dupire et al. Reference Dupire, Socol and Viallat2012), whereas in most simulations the cell reaches an elongated shape during tank treading (see, for example, Dodson & Dimitrakopoulos Reference Dodson and Dimitrakopoulos2010). To explain these findings, it was suggested that the reference state may be a spheroid so that the energy barrier between dimple and rim elements is greatly reduced (Dupire et al. Reference Dupire, Socol and Viallat2012). This hypothesis has been confirmed numerically through a boundary-element model (Tsubota, Wada & Liu Reference Tsubota, Wada and Liu2013), a multiscale model (Peng et al. Reference Peng, Mashayekh and Zhu2014) and a front-tracking immersed-boundary model (Cordasco, Yazdani & Bagchi Reference Cordasco, Yazdani and Bagchi2014). Cordasco et al. (Reference Cordasco, Yazdani and Bagchi2014) studied the phase diagram of off-shear plane dynamics, including tank treading, tumbling, rolling, kayaking and hovering, for low-shear-rate flows. It was shown that kayaking motion occurs when the capillary number is just above the critical value while the viscosity ratio ${\it\Lambda}$ is below the physiological range. Furthermore, they observed a hovering motion in which the cell axis first moves towards the shear plane and then aligns at a constant angle with respect to the shear plane, and the cell exhibits a combination of tank treading and a rigid-body-like spinning, when the cell is initially aligned with a large off-shear plane angle. This hovering motion is similar to the frisbee-like hovering motion observed in the experiment of Dupire et al. (Reference Dupire, Socol and Viallat2012). Nevertheless, due to the small range of shear rate considered in these studies the difference between the two reference states is subtle. This may lead to ambiguity in the conclusion. To further test this conclusion, it will be helpful to find additional mechanisms that can be easily tested through experiments.

The purpose of this study is to numerically study the tank treading motion of RBCs, especially the stability of the two modes, over a wide range of parameters (capillary number and viscosity ratio), and identify additional characteristics that can be used as guidance in the search for the reference state in future experiments. Towards this end, we will apply a multiscale fluid–cell interaction model developed earlier and reported in our previous publication (Peng, Asaro & Zhu Reference Peng, Asaro and Zhu2011).

The rest of the paper will be organized as follows. We first present the physical problem associated with an RBC in a linear shear flow. This is followed by a brief description of the numerical method and its validations (the details are included in our previous publications). Numerical results, including characteristics of the two tank treading modes, their stability properties and in particular the effect of the reference state, are included in the results section. These are followed by an intuitive discussion about the underlying physical mechanisms. Finally, conclusions are drawn.

2. Description of the physical problem

The set-up of our numerical simulations is shown in figure 1. We consider the dynamics of an RBC within a linear shear flow in the $x$ direction with the velocity profile $U(y)=\dot{{\it\gamma}}y$ , where $\dot{{\it\gamma}}$ is the shear rate. In this scenario, the dynamics of the cell is determined mostly by the following dimensionless parameters: (a) the capillary number $Ca$ , defined as $Ca\equiv {\it\eta}_{1}R\dot{{\it\gamma}}/{\it\mu}_{s}$ , where ${\it\eta}_{1}$ is the viscosity of the fluid outside of the cell, $R$ is the effective radius of the cell (i.e. the radius of a sphere with the same surface area as the cell – it is chosen to be $3.25~{\rm\mu}\text{m}$ in this study) and ${\it\mu}_{s}$ is the shear stiffness of the skeleton; (b) the viscosity ratio ${\it\Lambda}$ , defined as ${\it\Lambda}\equiv {\it\eta}_{2}/{\it\eta}_{1}$ , where ${\it\eta}_{2}$ is the viscosity of the internal fluid (i.e. the cytosol). The value of ${\it\eta}_{2}$ is chosen to be $0.006~\text{pN}~\text{s}~{\rm\mu}\text{m}^{-2}$ (Chien Reference Chien1987). We note that in our multiscale model (and also in the real cell), ${\it\mu}_{s}$ is a nonlinear function of the skeleton deformation so that it varies both spatially and temporally. To avoid any ambiguity, in calculating the capillary number the value ${\it\mu}_{s}$ is extracted from the undeformed cytoskeleton.

Figure 1.

Numerical set-up of an RBC in a linear shear flow. The membrane-attached bead marks the dimple area.

As mentioned in the introduction, depending on the aforementioned dimensionless parameters ( $Ca$ and ${\it\Lambda}$ ), the cell may display tumbling, wheel-like rolling, frisbee-like and tank treading motions. At certain viscosity ratios, tumbling, rolling and frisbee-like motions occur at relatively low values of the capillary number. These motions, especially their transitions to the tank treading motion, have been investigated extensively in recent numerical studies (Yazdani & Bagchi Reference Yazdani and Bagchi2011; Tsubota et al. Reference Tsubota, Wada and Liu2013; Cordasco et al. Reference Cordasco, Yazdani and Bagchi2014; Peng et al. Reference Peng, Mashayekh and Zhu2014). In the current study, we will concentrate on higher capillary numbers ( $Ca=O(1)$ ), where the tank treading motion dominates. During tank treading motion, in addition to shape changes the cell may also undergo oscillatory rotations around the $x$ , $y$ and $z$ axes, which are depicted as ${\it\theta}(t)$ , ${\it\phi}(t)$ (kayaking) and ${\it\psi}(t)$ (swinging) respectively. In practice, these angles are defined based on the orientation of a unit vector $\boldsymbol{l}=(l_{x},l_{y},l_{z})$ along the longest axis of the deformed cell, so that we have ${\it\theta}=\tan ^{-1}(l_{z}/l_{y})$ , ${\it\phi}=\tan ^{-1}(l_{z}/l_{x})$ and ${\it\psi}=\tan ^{-1}(l_{y}/l_{x})$ . The initial values of these angles, ${\it\theta}_{0}$ , ${\it\phi}_{0}$ and ${\it\psi}_{0}$ , represent the initial orientation of the cell.

The physical properties of the lipid bilayer and the cytoskeleton used in this study are summarized in tables 1 and 2. In table 2 the shear modulus ${\it\mu}_{s}$ is obtained from the undeformed cytoskeleton. The friction coefficient between the lipid bilayer and the skeleton is chosen to be $c_{f}=({\it\rho}/{\it\rho}_{0})144~\text{pN}~\text{s}~{\rm\mu}\text{m}^{-3}$ , where ${\it\rho}$ and ${\it\rho}_{0}$ are the current and initial protein densities of the skeleton respectively (Peng et al. Reference Peng, Asaro and Zhu2011).

Table 1.

Parameters of the lipid bilayer. Here, $h_{b}$ is the bilayer thickness (different from reality due to the homogeneous shell assumption), ${\it\mu}_{b}$ is the bilayer shear stiffness (a very small value to stabilize the numerical algorithm), $K_{b}$ is the bilayer areal stiffness (Mohandas & Evans Reference Mohandas and Evans1994),  $k_{c}$ is the bilayer bending stiffness (Mohandas & Evans Reference Mohandas and Evans1994) and ${\it\nu}_{b}$ is the bilayer viscosity (Otter & Shkulipa Reference Otter and Shkulipa2007).

Table 2.

Parameters of the cytoskeleton. Here, $h_{s}$ is the cytoskeleton thickness, ${\it\nu}_{s}$ is the cytoskeleton viscosity (Tran-Son-Tay, Sutera & Rao Reference Tran-Son-Tay, Sutera and Rao1984),  $p_{f}$ is the persistence length of folded domains in spectrin, $p_{u}$ is the persistence length of unfolded domains (Zhu & Asaro Reference Zhu and Asaro2008),  $L_{f}$ is the contour length of folded domains (Zhu & Asaro Reference Zhu and Asaro2008),  $L_{u}$ is the contour length of unfolded domains (Zhu & Asaro Reference Zhu and Asaro2008),  ${\rm\Delta}{\rm\Delta}x^{\ast }$ is the the difference between the activation length of the unfolding process and that of the refolding process (Zhu & Asaro Reference Zhu and Asaro2008),  $F_{1/2}$ is the the force corresponding to the state when half of the domains are unfolded (Zhu & Asaro Reference Zhu and Asaro2008) and ${\it\mu}_{s}$ is the initial shear modulus of the cytoskeleton (Peterson, Strey & Sackmann Reference Peterson, Strey and Sackmann1992). A spectrin consists of 19 domains in our model.

Over the past few years, we have developed a fluid–structure interaction model to simulate the dynamics of RBCs in shear flows. In the structural part, we employ a multiscale framework to represent the viscoelastic responses of the cell membrane at different length scales. In the fluid part, we apply a boundary-element method to simulate the dynamics of the fluids inside and outside of the cell based on the Stokes flow assumption. A brief review of the primary characteristics of these models is included in appendix A.

3. Results

In the following simulations, we use ABAQUS (ABAQUS Inc., Providence, RI) to generate the computational mesh for both the finite-element solver (the level III model) and the boundary-element solver. On each layer, 1000 elements are applied. The number of boundary elements on the outer surface is also 1000. For numerical stability, a time step of $2\times 10^{-5}$ s is used.

3.1. Resting shape

When the reference shape of the cytoskeleton is biconcave, it is straightforward to set the resting geometry of the cell as (Evans & Fung Reference Evans and Fung1972)

(3.1) $$\begin{eqnarray}y=\pm 0.5R_{0}\sqrt{1-\frac{x^{2}+z^{2}}{R_{0}^{2}}}\left[C_{1}+C_{2}\frac{x^{2}+z^{2}}{R_{0}^{2}}+C_{3}\left(\frac{x^{2}+z^{2}}{R_{0}^{2}}\right)^{2}\right],\end{eqnarray}$$

where $C_{1}=0.21$ , $C_{2}=2.03$ , $C_{3}=-1.12$ and $R_{0}=3.95~{\rm\mu}\text{m}$ . The internal volume is $0.65V_{0}$ , where $V_{0}$ is the volume of a sphere with the same surface area (i.e. a sphere whose radius is $R=3.25~{\rm\mu}\text{m}$ ). The skeleton is prescribed to be stress-free at this state.

The determination of the resting state of a cell with a spheroidal reference shape requires a deflation process (Peng et al. Reference Peng, Mashayekh and Zhu2014). During this process, we start from a spheroid-shaped cell whose internal volume is $V_{s}$ and surface area is $4{\rm\pi}R^{2}$ . It is easy to see that the eccentricity of this spheroid can be uniquely extracted from $V_{s}/V_{0}$ , the ratio between its internal volume and the volume of a sphere with radius $R$ . Hereafter this volume ratio is used as the parameter to characterize a spheroidal shape. At this initial state there is no stress in the skeleton so that this spheroid coincides with the reference shape. The internal volume of the cell is then gradually reduced from $V_{s}$ to 0.65 $V_{0}$ while its surface area does not change. Thus obtained, the resting cell shape is shown to be dependent on the initial stress-free shape and the spontaneous curvature of the lipid bilayer $C_{0}$ (which is usually expressed in a dimensionless form as $c_{0}\equiv C_{0}R$ ).

Systematic studies have been conducted using a Monte Carlo approach to document the resting shapes of the cell at different values of $V_{s}/V_{0}$ and $c_{0}$ (Lim et al. Reference Lim, Wortis and Mukhopadhyay2002; Lim, Wortis & Mukhopadhyay Reference Lim, Wortis and Mukhopadhyay2008). According to these simulations, to create the biconcave resting shape of RBCs, the reference shapes are most likely within the range $0.925<V_{s}/V_{0}<0.976$ . At each $V_{s}/V_{0}$ , as $c_{0}$ increases beyond a threshold value ( $c_{0}^{\ast }$ ) the resting shape switches from a stomatocyte (bowl shape) to a discocyte (biconcave shape). However, if the $V_{s}/V_{0}$ is too close to 1, no discocyte can be reached at any value of $c_{0}$ . These tendencies have been quantitatively confirmed in our previous study using the multiscale model with the aforementioned deflation method (Peng et al. Reference Peng, Mashayekh and Zhu2014). Furthermore, our simulations show that to duplicate the experimentally observed transition from tumbling/rolling motions to tank treading motion at low capillary numbers, the value of $V_{s}/V_{0}$ may be a little higher than the range proposed by Lim et al. (Reference Lim, Wortis and Mukhopadhyay2002) (a similar conclusion was reached by Cordasco et al. (Reference Cordasco, Yazdani and Bagchi2014) using a front-tracking immersed-boundary model).

To study cells with spheroidal reference states we will first concentrate on a case with $V_{s}/V_{0}=0.975$ and $c_{0}=8$ (although other values of $V_{s}/V_{0}$ and $c_{0}$ will also be explored later) and compare it with the case in which the reference shape is biconcave (in which $V_{s}/V_{0}=0.65$ and $c_{0}=0$ ). The resting shape of the cell and the distributions of skeleton deformations in this case are shown in figure 2. Hereby the in-plane deformations of the skeleton are quantified through two independent parameters, ${\it\alpha}\equiv {\it\lambda}_{1}{\it\lambda}_{2}$ and ${\it\beta}\equiv {\it\lambda}_{1}/{\it\lambda}_{2}$ , where ${\it\lambda}_{1}$ and ${\it\lambda}_{2}$ are principal in-plane stretches and by definition ${\it\lambda}_{1}\geqslant {\it\lambda}_{2}$ . The parameter ${\it\alpha}$ represents area change ( ${\it\alpha}>1$ corresponds to area expansion and ${\it\alpha}<1$ to area compression) and ${\it\beta}$ represents shear deformation. In the case shown in figure 2, it is seen that in the resting state the cytoskeleton is compressed in the area surrounding the dimple and expanded near the rim. The maximum area expansion ${\it\alpha}_{max}$ (associated with reduced skeleton density) is approximately 1.07, occurring on the rim. The maximum shear deformation ${\it\beta}_{max}$ also occurs near the rim. Its value is approximately 1.22.

Figure 2.

(a) The shape, (b) the area deformation of the skeleton and (c) the shear deformation of the skeleton at the resting configuration. Here, $V_{s}/V_{0}=0.975$ and $c_{0}=8$ .

Figure 3.

Snapshots of the cell in tank treading motion in (a) mode 1 and (b) mode 2. Two beads are attached to the membrane to mark the dimple elements. Here, $Ca=3.51$ and ${\it\Lambda}=1.0$ . Mode 1 is achieved with ${\it\theta}_{0}=0^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$ . Mode 2 is achieved with ${\it\theta}_{0}=90^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$ .

3.2. Characteristics of the two tank treading modes

As reported in Peng & Zhu (Reference Peng and Zhu2013), through numerical simulations we have shown that the tank treading motion actually includes two distinguishable modes. As demonstrated in figure 3, the most pronounced difference between these two modes is in the location and motion of the membrane elements originating from the dimple area (i.e. the area marked by the bead in figure 1). In the first mode (referred to as mode 1), these dimple elements remain close to the central plane (the $x{-}y$ plane) so that they travel the longest distance during tank treading (see figure 3 a). In the second mode (mode 2), the dimple elements move to the farthermost points from the central plane and their translational motion during tank treading is minimized (figure 3 b). Incidentally, similar features (including the existence of two distinguishable tank treading modes of erythrocytes in shear flow, the dependence upon initial orientations and the key characteristics of these modes) have also been observed by Cordasco et al. (Reference Cordasco, Yazdani and Bagchi2014) in their simulations using a front-tracking immersed-boundary model. This phenomenon is also similar to findings by Dupont, Salsac & Barthès-Biesel (Reference Dupont, Salsac and Barthès-Biesel2013) in the off-plane motion of a prolate capsule in shear flow.

Another important difference between the two tank treading modes lies in the magnitudes of the breathing and swinging motions. Hereby the breathing motion refers to the periodic oscillation in the aspect ratio of the cell geometry (usually characterized by variations of the Taylor deformation parameter $D_{xy}$ ) and the swinging motion refers to the oscillatory pitch around the $z$ axis (the ${\it\psi}(t)$ motion). Both the breathing and the swinging motions have twice the frequency of the tank treading motion. From figure 4 it is seen that mode 1 is characterized by significant breathing and swinging. In mode 2, on the other hand, both breathing and swinging are usually negligibly small.

Figure 4.

Time histories of the inclination angle ${\it\Psi}$ (a) and Taylor deformation parameter $D_{xy}$ (b) of mode 1 (solid lines) and mode 2 (dashed lines). Here, $Ca=3.51$ and ${\it\lambda}=1.0$ . Mode 1 is achieved with ${\it\theta}_{0}=0^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$ . Mode 2 is achieved with ${\it\theta}_{0}=90^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$ . The reference state is biconcave.

3.3. Stability of the modes and its relation with the reference state

At the initial stage of its tank treading motion, the cell motion is determined mostly by its original orientation (especially the roll angle ${\it\theta}$ ). Generally speaking, when the cell is initially oriented towards the $x{-}z$ plane (i.e. the value of ${\it\theta}_{0}$ is close to $0^{\circ }$ ), it displays mode 1 motion when the shear flow is turned on. On the other hand, when the initial orientation of the cell is more towards the $x{-}y$ plane (i.e. the value of ${\it\theta}_{0}$ is close to $90^{\circ }$ ), it will start with mode 2 motion. After long-term evolutions, however, the cell motion may eventually switch to a different mode. Indeed, in our simulations it has been discovered that both mode 1 and mode 2 can become unstable and switch to the other mode.

Figure 5 shows a typical scenario when mode 1 becomes unstable and switches to mode 2. The most pronounced change is that the dimple elements, represented by the beads, drift from the central plane to the two sides of the cell. This is clearly seen in the trajectories $z_{b}$ of the two beads in the $z$ direction (figure 5 a – each curve in that figure corresponds to the trajectory of one bead). Meanwhile, there are dramatic reductions in the oscillations of both $D_{xy}$ (breathing) and ${\it\psi}$ (swinging) (figure 5 b,c). According to figure 5(d), this transition process is also accompanied by an increase in the kayaking motion ${\it\phi}(t)$ , which eventually dies out when the cell settles down at mode 2. To make sure that the observed instability was not caused by accumulated error due to time integration, we conducted sensitivity tests by varying the time step ${\rm\Delta}t$ and found that the stability/instability behaviour was not sensitive to ${\rm\Delta}t$ .

Figure 5.

A transition event from mode 1 to mode 2. (a) Trajectories of the dimple elements in the $z$ direction. (b) The breathing motion characterized by variations of the Taylor parameter $D_{xy}(t)$ . (c) The swinging motion ${\it\Psi}(t)$ . (d) The kayaking motion ${\it\phi}(t)$ . Here, $Ca=2.34$ ,  ${\it\Lambda}=1.0$ and ${\it\theta}_{0}=0^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$ . The reference state is biconcave.

Figure 6 demonstrates the opposite process of transition from mode 2 to mode 1. As shown in the figure, the instability of mode 2 is characterized by the relocation of the dimple elements to the central plane (as shown by the trajectories $z_{b}$ of the dimple elements in figure 6 a). This is accompanied by large increases in breathing (figure 6 b) and swinging (figure 6 c). Large kayaking motion is also observed during this process (figure 6 d).

Figure 6.

A transition event from mode 2 to mode 1. (a) Trajectories of the dimple elements in the $z$ direction. (b) The breathing motion characterized by variations of the Taylor parameter $D_{xy}(t)$ . (c) The swinging motion ${\it\Psi}(t)$ . (d) The kayaking motion ${\it\phi}(t)$ . Here, $Ca=3.90$ ,  ${\it\Lambda}=0.5$ and ${\it\theta}_{0}=90^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$ . The reference state is biconcave.

We conducted systematic simulations to document the stability of the two tank treading modes on the ${\it\Lambda}{-}Ca$ plane. Specifically, we focused on the range of ${\it\Lambda}$ from 0.5 to 1.0. At each viscosity ratio, the lower limit of the stable region of mode 1 was determined (through a bisection approach) as $Ca^{-}$ , and the upper limit of the stable region of mode 2 was determined as $Ca^{+}$ . Typical stability diagrams including both the $Ca^{-}$ versus ${\it\Lambda}$ curves and the $Ca^{+}$ versus ${\it\Lambda}$ curves are plotted in figure 7. In the region below the $Ca^{-}$ curve (marked as region III in the figure), the only stable mode is mode 2. Within this region, even if the cell starts with a mode 1 motion at certain initial conditions, it will eventually transform into mode 2 via the process illustrated in figure 5. Between the $Ca^{-}$ and $Ca^{+}$ curves (region II) both mode 1 and mode 2 are stable so that the system displays a bifurcation behaviour. The exact response depends on the initial condition. In fact this is the behaviour we concentrated on in our previous work (Peng & Zhu Reference Peng and Zhu2013). Above the $Ca^{+}$ curve (region I) the system response is dominated by mode 1. This is the region where the transition process illustrated in figure 6 may occur.

Figure 7.

Stability diagrams in the ${\it\Lambda}{-}Ca$ plane of (a) an RBC with biconcave reference state and (b) an RBC with spheroidal reference state ( $V_{s}/V_{0}=0.975$ , $c_{0}=8$ ).

An enormous difference exists between a cell with a biconcave reference state and one with a spheroidal reference state in terms of the stability properties of the two tank treading modes. This is clearly demonstrated in the comparison between figures 7(a) and 7(b) (note that the vertical scales in these two figures are different). Specifically, for a cell with spheroidal reference state (figure 7 b) region III (the region in which only mode 2 is stable) is greatly diminished, and region II (the bistable region) almost disappears. The implication is that the spheroidal reference state discourages the occurrence of mode 2.

It is more illustrative to express this difference in terms of shear rates instead of capillary numbers. For example, with the physical parameters chosen in this study (see tables 1 and 2), at ${\it\Lambda}=1.0$ for a cell with biconcave reference state ( $V_{s}/V_{0}=0.65$ and $c_{0}=0$ ) mode 1 is stable when the shear rate is above $425~\text{s}^{-1}$ . If the cell has a spheroidal reference state with $V_{s}/V_{0}=0.975$ and $c_{0}=8$ , this value is greatly reduced to approximately $165~\text{s}^{-1}$ . Similarly, the highest shear rate for mode 2 to be stable is reduced from approximately $595~\text{s}^{-1}$ for the cell with biconcave reference state to approximately $168~\text{s}^{-1}$ for the cell with spheroidal reference state. Such large differences will not be difficult to detect in experiments. Our results on multiple tank treading modes and their stability properties will provide a valuable measure in future experiments to determine whether the reference state is biconcave or spheroidal.

Another interesting question is whether or not this phenomenon can be applied to pinpoint the exact reference state. To shed light on this question, we have conducted additional numerical tests to show the sensitivity of $Ca^{-}$ and $Ca^{+}$ to the detailed eccentricity of the spheroidal reference state. Towards this end, in figure 8 we plot the variations of $Ca^{-}$ and $Ca^{+}$ over the range $0.90\leqslant V_{s}/V_{0}\leqslant 0.995$ ( $c_{0}$ is fixed as 8 and ${\it\Lambda}$ is fixed as 0.75). It is seen that both $Ca^{-}$ and $Ca^{+}$ decrease quickly with $V_{s}/V_{0}$ , especially when it approaches 1 (i.e. the reference state approaches a sphere). Therefore, carefully measured values of these critical capillary numbers may indeed be used as criteria to determine the exact value of $V_{s}/V_{0}$ .

Figure 8.

Stability diagram in the $V_{s}/V_{0}\sim Ca$ plane for an RBC with spheroidal reference state. Here, ${\it\Lambda}=0.75$ and $c_{0}=8$ .

Figure 9.

Schematic illustration of the membrane contour within the $x{-}y$ plane during tank treading motions for (a) mode 1 and (b) mode 2 over half a tank treading period (the time increases sequentially from $t_{1}$ to $t_{5}$ ). The dimple elements are shown as circles and the rim elements are shown as bullets.

A remaining concern is that the spontaneous curvature of the lipid bilayer characterized by the parameter $c_{0}$ , which is still an unknown factor, may also affect the cell dynamics during tank treading motion. According to the Monte Carlo simulation by Lim et al. (Reference Lim, Wortis and Mukhopadhyay2002) and our own study (Peng et al. Reference Peng, Mashayekh and Zhu2014), at each value of $V_{s}/V_{0}$ there exists a critical value of $c_{0}$ ( $c_{0}^{\ast }$ ), beyond which the resting shape of the cell switches from a stomatocyte to a discocyte. Although $c_{0}^{\ast }$ , which depicts the lower bound of the possible range of $c_{0}$ , can be determined from simulations, there is no existing knowledge about the exact value of $c_{0}$ itself. To address this effect, we conducted a sensitivity test using the case when $V_{s}/V_{0}=0.975$ , ${\it\Lambda}=0.75$ and $8\leqslant c_{0}\leqslant 12$ . According to our results, within this range of $c_{0}$ the variations of both $Ca^{-}$ and $Ca^{+}$ are non-significant.

4. Discussion

To physically understand the difference between mode 1 and mode 2, we concentrate on the contour of the cell membrane within the central plane (the $x{-}y$ plane). The membrane elements within this plane travel the longest distances during tank treading. As illustrated in figure 9(a), in mode 1 these elements have different origins. Some of them are from the dimple area (depicted as circles), whereas some are from the rim area (represented by bullets). There are significant differences between the strain energy densities stored in these elements from different origins (Dupire et al. Reference Dupire, Socol and Viallat2012). Specifically, the dimple elements contain less energy than the rim elements. The identities of these elements will be preserved even in a deformed cell during tank treading (Peng & Zhu Reference Peng and Zhu2013). This is epitomized by the fact that the cytoskeleton associated with dimple elements and that with rim elements will have different deformations during tank treading. As a consequence, following the circulation of these elements around the cell there is a periodic variation in the structural configuration of the cell. The period of this variation is half that of the tank treading motion itself. For example, in figure 9(a) the cell configurations at $t_{1}$ and $t_{3}$ are clearly different, while at $t_{5}$ (after half a tank treading period) the cell returns to the initial configuration at $t_{1}$ . This variation contributes to the breathing and swinging oscillations of the cell. It also explains why the frequencies of these oscillations are exactly twice that of the tank treading motion. In mode 2, on the other hand, all the elements on this contour are from the rim area so that the structural configuration of the cell remains unchanged during tank treading (figure 9 b). This explains the significant reduction in breathing and swinging oscillations. In terms of the motions of the rim elements, mode 2 appears to be an ‘extension’ of the wheel-like rolling mode into the high-capillary-number regime. The only difference is that mode 2 is associated with significant cell deformation, while in the rolling mode the cell maintains its original biconcave shape.

Figure 10.

Distributions of the strain energy density $E$ (normalized by ${\it\eta}_{1}\dot{{\it\gamma}}R$ ) along the contour within the $x{-}y$ plane in the membrane as the cell switches from (a) mode 1 (within half of a breathing/swinging period starting at $\dot{{\it\gamma}}t=67.5$ ) to (b) mode 2 (at $\dot{{\it\gamma}}t=1440$ ). This is the same case as defined in figure 5, which lies in region III in the stability diagrams. The curves in (a) show four evenly distributed snapshots over time (the variations are caused by the breathing/swinging oscillations in mode 1). In mode 2 the cell shape and configuration remain almost steady so that there is no variation of the strain energy density over time. The inset in (a) describes the definition of the angle ${\it\Theta}$ .

Figure 11.

Distributions of the strain energy density $E$ (normalized by ${\it\eta}_{1}\dot{{\it\gamma}}R$ ) along the contour within the $x{-}y$ plane in the membrane as the cell switches from (a) mode 2 (at $\dot{{\it\gamma}}t=125$ ) to (b) mode 1 (within half of a breathing/swinging period starting at $\dot{{\it\gamma}}t=2000$ ). This is the same case as defined in figure 6, which lies in region I in the stability diagrams.

Furthermore, in mode 1 since the elements within the $x{-}y$ plane are from different origins, along the contour there will be larger differences (compared with mode 2) in terms of the distribution of strain energy (for example, see figures 10 and 11). This may create higher energy barriers to achieve tank treading. Meanwhile, the breathing/swinging motions in mode 1 also consume significant energy in this highly dissipative flow regime. As a result, mode 1 can be considered as a high-energy mode while mode 2 is a low-energy mode. This may be the reason why mode 2 is favoured when the energy input is low (hereby in flows with relatively low shear rates) and mode 1 dominates when the energy input is sufficiently high (i.e. in flows with relatively high shear rates).

The power input from the external fluid has to be balanced by the internal dissipation rate and the rate of elastic energy storage (Skotheim & Secomb Reference Skotheim and Secomb2007). As the membrane elements circulate around the cell, their strain energy increases (i.e. energy extraction) when they go from low-energy regions to high-energy regions, and decreases (energy release) when they return. Therefore, the system contains numerous spring–damper oscillators (each membrane element is an oscillator) that store and release potential energy, and dissipate energy all the time. When the other conditions are the same, mode 1 provides a higher rate of elastic energy storage (i.e. it contains ‘stiffer springs’) than mode 2. In addition, due to the additional breathing motions, mode 1 may also have a higher rate of internal dissipation. When the power input from the external fluid is increased (with increased shear rate), it becomes hard for mode 2 to provide enough ‘room’ for energy dissipation and elastic energy storage (as demonstrated in figure 6, at the same shear rate mode 2 involves larger cell deformation), so that it may become unstable to disturbances and switch to mode 1. This is somewhat similar to the transition from tumbling to tank treading. Since tumbling is close to rigid-body motion, it shows smaller elastic energy and smaller dissipation than tank treading. When the shear rate is increased, the motion switches from the lower-elastic-energy tumbling mode to the higher-elastic-energy tank treading mode. To a certain extent, mode 1 and mode 2 are also analogous to the ‘orbits’ of electrons of an atom. With high energy input more electrons will be excited to higher-level orbits.

Compared with a cell with a biconcave reference state, the one with a spheroidal reference state is characterized by greatly reduced anisotropy on its membrane so that the differences between dimple and rim elements are also mitigated (Dupire et al. Reference Dupire, Socol and Viallat2012; Cordasco et al. Reference Cordasco, Yazdani and Bagchi2014). It therefore takes much less energy to achieve mode 1. This may provide an explanation for the enlarged region dominated by mode 1 in the stability diagrams when the reference state is spheroidal rather than biconcave.

5. Conclusions

By using a fluid–cell interaction model based on the coupling of a boundary-element fluid solver with a multiscale structural model that includes the molecular details, we numerically investigated the tank treading motions of an RBC in linear shear flows. The results confirmed that there existed two different tank treading modes. In mode 1 the dimple elements circulate around the cell near the central plane, whereas in mode 2 these elements go to the farthest points from the plane so that they undergo minimum circulation motions. Compared with mode 2, mode 1 is also characterized by much more significant oscillations in its shape (breathing) and pitching angle (swinging).

Although the occurrence of mode 1 or mode 2 at the initial stage of tank treading is dependent on the initial orientation of the cell, during evolution the cell motion may become unstable and eventually switch to the other mode. When other parameters (e.g. the viscosity ratio) are fixed, there exist two critical capillary numbers, $Ca^{-}$ and $Ca^{+}$ . At low capillary numbers ( $Ca<Ca^{-}$ ), only mode 2 is stable no matter what initial conditions are applied. When $Ca$ is between $Ca^{-}$ and $Ca^{+}$ , both mode 1 and mode 2 are stable. Finally, when $Ca>Ca^{+}$ , mode 1 becomes the only stable mode. Through systematic simulations, the stability diagrams in the ${\it\Lambda}{-}Ca$ plane have been obtained. Specifically, in that plane we have identified a mode 2-dominated region (region III), a bistable region (region II) and a mode 1-dominated region (region I).

An important finding in this study is that the stability diagram of the tank treading modes is highly sensitive to the reference state (the stress-free state) of the cytoskeleton. For example, compared with a cell with a biconcave reference state, a cell with a spheroidal reference state displays much smaller regions II and III. Moreover, as the reference state approaches a sphere, both $Ca^{-}$ and $Ca^{+}$ decrease dramatically so that mode 1 dominates the cell dynamics. These results suggest that the occurrence and stability of the tank treading modes may provide a measure to experimentally determine the exact reference state of the cytoskeleton of RBCs.

It is necessary to point out that the observed tank treading modes and their stability may also be obtained through a simplified RBC model in which the membrane is represented by a single layer. However, there will be differences between predictions made by a single-layer model and those by our two-layer model due to three major factors. First, the local area incompressibility of the lipid bilayer is missing in the single-layer models, although the global area incompressibility is usually considered in single-layer models. For instance, there are two kinds of micropipette experiments: one is to measure the shear modulus (Waugh & Evans Reference Waugh and Evans1979) and the other is to measure the local area modulus of the bilayer (Evans, Waugh & Melnik Reference Evans, Waugh and Melnik1976). Only the two-layer model can reproduce both micropipette experiment results with the same set of parameters (Peng et al. Reference Peng, Li, Pivkin, Dao, Karniadakis and Suresh2013). Second, the mechanical property of the cytoskeleton in our model is derived from a multiscale model, including nonlinear constitutive behaviour with spectrin unfolding considered, which is extensively validated against experiments, so that it is more accurate. Third, the single-layer model does not correctly represent the viscoelasticity of the cell since it does not consider the viscous sliding between the lipid bilayer and the skeleton. This sliding mechanism may have important effects in dynamical processes. For example, without considering this effect it is not possible to study the transient process of the build-up of area deformations in the skeleton during tank treading motion (Peng & Zhu Reference Peng and Zhu2013).

Appendix A. A brief review of the numerical model

The mechanical responses of an RBC involve physical mechanisms at different length scales, ranging from dynamics at the complete-cell level (on the micron scale) down to local dynamics of the skeleton and tension-induced structural remodelling such as domain unfolding in proteins (which occur on the nanometre scale). To numerically duplicate these mechanisms, we developed a multiscale approach including three models at different length scales. These models are as follows.

  1. (a) A force-elongation model of spectrin (Sp) (level I). Spectrin is the protein that determines the constitutive behaviour of the cytoskeleton. As a highly flexible biopolymer, it is usually modelled by using the worm-like-chain (WLC) depiction (Discher et al. Reference Discher, Boal and Boey1998; Li et al. Reference Li, Dao, Lim and Suresh2005). This depiction, however, does not take into account the fact that Sp can be overstretched, i.e. it can be stretched beyond its original contour length due to the existence of multiple domains that unfold under sufficient loading. Without proper representation of this process it is impossible to accurately describe the mechanics of the cytoskeleton during large deformations. Towards this end, we have developed a force-elongation model of Sp by representing its folding/unfolding reactions as thermally activated processes (Zhu & Asaro Reference Zhu and Asaro2008).

  2. (b) A molecular-based model of the JC (level II). Consisting of a central piece of actin protofilament connected with (up to) six Sp, the JC is the basic repeating unit in the cytoskeleton of an erythrocyte. A molecular-based model that incorporates the state-of-the-art understanding of the exact molecular architecture of a JC (including, e.g., the exact configuration of actin–spectrin connectivity) has been developed. With this model we are able to predict the area and shear moduli of the skeleton, as well as simulate mesoscale mechanics such as the bistability of a JC (Zhu et al. Reference Zhu, Vera, Asaro, Sche and Sung2007).

  3. (c) A complete-cell model (level III). A key structural feature of the RBC is the mobile connectivity between its lipid bilayer and its cytoskeleton. This is caused by the fact that these two components are connected through pinning points made of transmembrane proteins, which can drift within the fluid-like lipid bilayer. To numerically duplicate this characteristic, in the level III model we depict the cell membrane as two separate layers. The outer layer corresponds to the lipid bilayer, which has large area stiffness, negligible shear stiffness and finite bending stiffness. The inner layer, with finite area and shear moduli, represents the cytoskeleton. The two layers remain in contact, but they are allowed to slide against each other with viscous friction. The viscous coefficient is evaluated by considering the distribution of transmembrane proteins (band 3 and glycophorin C) as well as their mobilities inside the lipid bilayer. Numerically, a finite-element approach based on the Hughes–Liu shell element is employed to solve the mechanics of both the outer and the inner layers (Peng, Asaro & Zhu Reference Peng, Asaro and Zhu2010).

In our multiscale framework, these three models are connected through an information-passing technique, in which the force-elongation relation of Sp predicted by the level I model is used by the level II model, and the constitutive properties of the skeleton (the area and shear moduli) predicted by the level II model are used by the level III model for the inner layer.

An essential mechanical property of the RBC is its viscoelasticity. This property is important for its capacity to sustain large dynamic loads. The elasticity is attributed to the elastic behaviours of the lipid bilayer and the cytoskeleton, characterized by the area stiffness, the shear stiffness and the bending stiffness of each component. These two structural members also possess viscosity, which is modelled by using the generalized Voigt–Kelvin stress–strain relation (Evans & Skalak Reference Evans and Skalak1980). In addition, the sliding between the bilayer and the skeleton also contributes to the viscosity of the system. This effect is modelled as a viscous friction force (see the description of the level III model). The viscosities of the internal and external fluids are included in the fluid-dynamics model described below.

The effects of the fluids outside and inside the cell are simulated by coupling the finite-element model (i.e. the level III model) with a boundary-element approach for the Stokes flow fields inside and outside of the cell. The finite-element model–boundary-element model coupling is achieved through a staggered algorithm with explicit time integration (Peng et al. Reference Peng, Asaro and Zhu2011). Under the Stokes flow assumption, we employ the boundary integral equation of interface dynamics (Pozrikidis Reference Pozrikidis1992) so that the velocity field $\boldsymbol{v\,}^{f}$ satisfies

(A 1) $$\begin{eqnarray}\displaystyle \boldsymbol{v\,}^{f}(\boldsymbol{x}_{0}) & = & \displaystyle \frac{2}{1+{\it\Lambda}}\bar{\boldsymbol{v}}^{f}(\boldsymbol{x}_{0})-\frac{1}{4{\rm\pi}{\it\eta}_{1}({\it\Lambda}+1)}\iint _{{\it\Gamma}^{fb}}\boldsymbol{G}(\boldsymbol{x},\boldsymbol{x}_{0})\boldsymbol{\cdot }{\rm\Delta}\boldsymbol{t}^{f}(\boldsymbol{x})\text{d}{\it\Gamma}(\boldsymbol{x})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1-{\it\Lambda}}{4{\rm\pi}(1+{\it\Lambda})}\int \hspace{-9.70001pt}-\!\!\int _{{\it\Gamma}^{fb}}\boldsymbol{v\,}^{f}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D64F}(\boldsymbol{x},\boldsymbol{x}_{0})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x})\text{d}{\it\Gamma}(\boldsymbol{x}),\end{eqnarray}$$

where $\bar{\boldsymbol{v}}^{f}$ is the prescribed undisturbed velocity field of the shear flow in the absence of the cell, ${\it\Gamma}^{fb}$ is the fluid–solid boundary, the vector ${\rm\Delta}\boldsymbol{t}^{f}=\boldsymbol{t}^{f,1}-\boldsymbol{t}^{f,2}$ is the discontinuity in the interfacial surface traction, where $\boldsymbol{t}^{f,1}$ is the traction in the outside surface ${\it\Gamma}^{fb,1}$ of the interface and $\boldsymbol{t}^{f,2}$ is the traction in the inside surface ${\it\Gamma}^{fb,2}$ of the interface. The surface traction is related to the nodal force through the principle of virtual work (Walter et al. Reference Walter, Salsac, Barthès-Biesel and Le Tallec2010). Here, $\int \hspace{-9.70001pt}-\!\!\int$ denotes the principal value integration. The second term on the right-hand side of (A 1) is the single-layer potential, representing the contribution from the distribution of point forces associated with the Green’s function for velocity. The third term is the double-layer potential, representing the contributions from point sources and point force dipoles.

The matrix $\unicode[STIX]{x1D642}$ contains the free-space Green’s function for velocity $\unicode[STIX]{x1D60E}_{ij}$ . We have

(A 2) $$\begin{eqnarray}\unicode[STIX]{x1D60E}_{ij}(\boldsymbol{x},\boldsymbol{x}_{0})=\frac{{\it\delta}_{ij}}{|\boldsymbol{x}-\boldsymbol{x}_{0}|}+\frac{(x_{i}-x_{0_{i}})(x_{j}-x_{0_{j}})}{|\boldsymbol{x}-\boldsymbol{x}_{0}|^{3}},\end{eqnarray}$$

where ${\it\delta}_{ij}$ is Kronecker’s delta. The matrix $\unicode[STIX]{x1D64F}$ is the Green’s function for stress. Its components are

(A 3) $$\begin{eqnarray}\unicode[STIX]{x1D61B}_{ijk}(\boldsymbol{x},\boldsymbol{x}_{0})=-6\frac{(x_{i}-x_{0_{i}})(x_{j}-x_{0_{j}})(x_{k}-x_{0_{k}})}{|\boldsymbol{x}-\boldsymbol{x}_{0}|^{5}}.\end{eqnarray}$$

The details of our numerical model are presented in our previous publications (Zhu et al. Reference Zhu, Vera, Asaro, Sche and Sung2007; Zhu & Asaro Reference Zhu and Asaro2008; Peng et al. Reference Peng, Asaro and Zhu2010, Reference Peng, Asaro and Zhu2011; Peng & Zhu Reference Peng and Zhu2013; Peng et al. Reference Peng, Mashayekh and Zhu2014). The multiscale model, together with the fluid–cell interaction algorithm, has been created step by step. For each step we have conducted extensive tests, calibrations and validations. For example, to test the capability of the level I model in simulating the rate-dependent stretching of Sp, we conducted systematic simulations to obtain the force-elongation curves of Sp at different loading rates. The results are consistent with measurements from atomic force microscopy (Zhu & Asaro Reference Zhu and Asaro2008). At level II, we conducted simulations to record the constitutive properties (area and shear moduli) of the skeleton as functions of the persistence length of Sp. These were compared with results from various experiments and simulations (Zhu et al. Reference Zhu, Vera, Asaro, Sche and Sung2007). At level III, to validate the finite-element method we used the model to predict shapes of a vesicle with different volume to surface area ratios. The results match well with the theoretical predictions by Seifert, Berndl & Lipowsky (Reference Seifert, Berndl and Lipowsky1991). We also modelled cell deformations created by optical tweezers and micropipette aspirations. In both cases our results are consistent with measurements. Specifically, we demonstrated that the area deformation of the cytoskeleton predicted by our model was close to the measurements by Discher et al. (Reference Discher, Mohandas and Evans1994) using florescent microscopy (Peng et al. Reference Peng, Asaro and Zhu2010). Finally, the accuracy of the fluid–cell interaction model has also been tested in numerous scenarios, including cell motion inside a pipe, tumbling motion and tank treading. In all the tests good agreement with benchmark results in the literature has been obtained (Peng et al. Reference Peng, Asaro and Zhu2011; Peng & Zhu Reference Peng and Zhu2013).

References

Abkarian, M., Faivre, M. & Viallat, A. 2007 Swinging of red blood cells under shear flow. Phys. Rev. Lett. 98, 188302.Google Scholar
Barthès-Biesel, D. 1980 Motion of a spherical microcapsule freely suspended in a linear shear flow. J. Fluid Mech. 100, 831–853.Google Scholar
Bitbol, M. 1986 Red blood cell orientation in orbit $c=0$ . Biophys. J. 49, 1055–1068.CrossRefGoogle ScholarPubMed
Chien, S. 1987 Red cell deformability and its relevance to blood flow. Annu. Rev. Physiol. 49, 177.CrossRefGoogle ScholarPubMed
Cordasco, D., Yazdani, A. & Bagchi, P. 2014 Comparison of erythrocyte dynamics in shear flow under different stress-free configurations. Phys. Fluids 26, 041902.CrossRefGoogle Scholar
Discher, D. E., Boal, D. H. & Boey, S. K. 1998 Simulations of the erythrocyte cytoskeleton at large deformation. II. Micropipette aspiration. Biophys. J. 75, 1584–1597.Google Scholar
Discher, D. E., Mohandas, N. & Evans, E. A. 1994 Molecular maps of red cell deformation: hidden elasticity and in situ connectivity. Science 266, 1032–1035.Google Scholar
Dodson, W. R. & Dimitrakopoulos, P. 2010 Tank-treading of erythrocytes in strong shear flows via a nonstiff cytoskeleton-based continuum computational modeling. Biophys. J. 99, 2906–2916.CrossRefGoogle Scholar
Dupire, J., Socol, M. & Viallat, A. 2012 Full dynamics of a red blood cell in shear flow. Proc. Natl Acad. Sci. USA 109, 20808–20813.Google Scholar
Dupont, C., Salsac, A. & Barthès-Biesel, D. 2013 Off-plane motion of a prolate capsule in shear flow. J. Fluid Mech. 721, 180–198.Google Scholar
Eggleton, C. D. & Popel, A. S. 1998 Large deformation of red blood cell ghosts in a simple shear flow. Phys. Fluids 10, 1834–1845.Google Scholar
Evans, E. A. & Fung, Y. C. 1972 Improved measurements of the erythrocyte geometry. Microvasc. Res. 4, 335–347.Google Scholar
Evans, E. A. & Skalak, P. 1980 Mechanics and Thermodynamics of Biomembranes. CRC Press.Google Scholar
Evans, E. A., Waugh, R. E. & Melnik, L. 1976 Elastic area compressibility modulus of red cell membrane. Biophys. J. 16, 585–595.CrossRefGoogle ScholarPubMed
Fedosov, D. A., Caswell, B. & Karniadakis, G. E. 2011 Wall shear stress-based model for adhesive dynamics of red blood cells in malaria. Biophys. J. 100, 2084–2093.Google Scholar
Fischer, T. M. 2004 Shape memory of human red blood cells. Biophys. J. 86, 3304–3313.Google Scholar
Fischer, T. M. & Korzeniewski, R. 2013 Threshold shear stress for the transition between tumbling and tank-treading of red blood cells in shear flow dependence on the viscosity of the suspending medium. J. Fluid Mech. 736, 351–365.Google Scholar
Fischer, T. M., Stohr-Liesen, M. & Schmid-Schonbein, H. 1978 The red cell as a fluid droplet: tank tread-like motion of the human erythrocyte membrane in shear flow. Science 202, 894–896.CrossRefGoogle ScholarPubMed
Hoffman, J. F. 2001 Questions for red blood cell physiologists to ponder in this millennium. Blood Cells Mol. Dis. 27, 57–61.Google Scholar
Keller, S. & Skalak, R. 1982 Motion of a tank-treading ellipsoidal particle in a shear flow. J. Fluid Mech. 120, 27–47.CrossRefGoogle Scholar
Li, J., Dao, M., Lim, C. T. & Suresh, S. 2005 Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte. Biophys. J. 88, 3707.Google Scholar
Lim, G., Wortis, M. & Mukhopadhyay, R. 2002 Stomatocyte–discocyte–echinocyte sequence of the human red blood cell: evidence for the bilayer couple hypothesis from membrane mechanics. Proc. Natl Acad. Sci. USA 99, 16766–16769.Google Scholar
Lim, G., Wortis, M. & Mukhopadhyay, R. 2008 Soft Matter, Vol 4: Lipid Bilayers and Red Blood Cells. Wiley-VCH.Google Scholar
Mohandas, N. & Evans, E. A. 1994 Mechanical properties of the red cell membrane in relation to molecular structure and genetic defects. Annu. Rev. Biophys. Biomol. Struct. 23, 787–818.Google Scholar
Noguchi, H. 2009a Dynamic modes of red blood cells in oscillatory shear flow. Phys. Rev. E 81, 061920.Google Scholar
Noguchi, H. 2009b Swinging and synchronized rotations of red blood cells in simple shear flow. Phys. Rev. E 80, 021902.Google Scholar
Omori, T., Imai, Y., Yamaguchi, T. & Ishikawa, T. 2012a Reorientation of a nonspherical capsule in creeping shear flow. Phys. Rev. Lett. 108, 138102.Google Scholar
Omori, T., Ishikawa, T., Barthès-Biesel, D., Salsac, A., Imai, Y. & Yamaguchi, T. 2012b Tension of red blood cell membrane in simple shear flow. Phys. Rev. E 86, 056321.Google Scholar
Otter, W. K. & Shkulipa, S. A. 2007 Intermonolayer friction and surface shear viscosity of lipid bilayer membranes. Biophys. J. 93, 423–433.Google Scholar
Peng, Z., Asaro, R. & Zhu, Q. 2010 Multiscale simulation of erythrocyte membranes. Phys. Rev. E 81, 031904.Google Scholar
Peng, Z., Asaro, R. & Zhu, Q. 2011 Multiscale modeling of erythrocytes in Stokes flow. J. Fluid Mech. 686, 299–337.Google Scholar
Peng, Z., Li, X., Pivkin, I. V., Dao, M., Karniadakis, G. E. & Suresh, S. 2013 Lipid bilayer and cytoskeletal interactions in a red blood cell. Proc. Natl Acad. Sci. USA 110, 13356–13361.Google Scholar
Peng, Z., Mashayekh, A. & Zhu, Q. 2014 Erythrocyte responses in low shear rate flows – effects of non-biconcave stress-free state in cytoskeleton. J. Fluid Mech. 742, 96–118.CrossRefGoogle Scholar
Peng, Z. & Zhu, Q. 2013 Deformation of the erythrocyte cytoskeleton in tank treading motions. Soft Matt. 9, 7617–7627.Google Scholar
Peterson, M. A., Strey, H. & Sackmann, E. 1992 Theoretical and phase contrast microscopic eigenmode analysis of erythrocyte flicker: amplitudes. J. Phys. II France 2, 1273–1285.Google Scholar
Pozrikidis, C. 1992 Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press.Google Scholar
Seifert, U., Berndl, K. & Lipowsky, R. 1991 Shape transformations of vesicles: phase diagram for spontaneous-curvature and bilayer-coupling models. Phys. Rev. A 44, 1182–1202.CrossRefGoogle ScholarPubMed
Skotheim, J. M. & Secomb, T. W. 2007 Red blood cells and other nonspherical capsules in shear flow: oscillatory dynamics and the tank-treading-to-tumbling transition. Phys. Rev. Lett. 98, 078301.Google Scholar
Svoboda, K., Schmidt, C. F., Branton, D. & Block, S. M. 1992 Conformation and elasticity of the isolated red blood cell membrane skeleton. Biophys. J. 63, 784–793.Google Scholar
Tran-Son-Tay, R., Sutera, S. & Rao, P. 1984 Determination of red blood cell membrane viscosity from rheoscopic observations of tank-treading motion. Biophys. J. 46, 65–72.Google Scholar
Tsubota, K., Wada, S. & Liu, H. 2014 Elastic behavior of a red blood cell with the membrane’s nonuniform natural state: equilibrium shape, motion transition under shear flow, and elongation during tank-treading motion. Biomech. Model Mechanobiol. 13, 735–746.Google Scholar
Vlahovska, P., Young, Y. N., Danker, G. & Misbah, C. 2011 Dynamics of a non-spherical microcapsule with incompressible interface in shear flow. J. Fluid Mech. 678, 221–247.Google Scholar
Walter, J., Salsac, A. V., Barthès-Biesel, D. & Le Tallec, P. 2010 Coupling of finite element and boundary integral methods for a capsule in a Stokes flow. Intl J. Numer. Meth. Engng 83, 829.Google Scholar
Waugh, R. E. & Evans, E. A. 1979 Thermoelasticity of red blood cell membrane. Biophys. J. 26, 115–131.Google Scholar
Yazdani, A. & Bagchi, P. 2011 Phase diagram and breathing dynamics of a single red blood cell and a biconcave capsule in dilute shear flow. Phys. Rev. E 84, 026314.Google Scholar
Yazdani, A., Kalluri, R. & Bagchi, P. 2011 Tank-treading and tumbling frequencies of capsules and red blood cells. Phys. Rev. E 83, 046305.Google Scholar
Zhu, Q. & Asaro, R. 2008 Spectrin folding versus unfolding reactions and RBC membrane stiffness. Biophys. J. 94, 2529–2545.Google Scholar
Zhu, Q., Vera, C., Asaro, R., Sche, P. & Sung, L. A. 2007 A hybrid model for erythrocyte membrane: a single unit of protein network coupled with lipid bilayer. Biophys. J. 93, 386–400.Google Scholar
Figure 0

Figure 1.Numerical set-up of an RBC in a linear shear flow. The membrane-attached bead marks the dimple area.

Figure 1

Table 1.Parameters of the lipid bilayer. Here, $h_{b}$ is the bilayer thickness (different from reality due to the homogeneous shell assumption), ${\it\mu}_{b}$ is the bilayer shear stiffness (a very small value to stabilize the numerical algorithm), $K_{b}$ is the bilayer areal stiffness (Mohandas & Evans 1994), $k_{c}$ is the bilayer bending stiffness (Mohandas & Evans 1994) and ${\it\nu}_{b}$ is the bilayer viscosity (Otter & Shkulipa 2007).

Figure 2

Table 2.Parameters of the cytoskeleton. Here, $h_{s}$ is the cytoskeleton thickness, ${\it\nu}_{s}$ is the cytoskeleton viscosity (Tran-Son-Tay, Sutera & Rao 1984), $p_{f}$ is the persistence length of folded domains in spectrin, $p_{u}$ is the persistence length of unfolded domains (Zhu & Asaro 2008), $L_{f}$ is the contour length of folded domains (Zhu & Asaro 2008), $L_{u}$ is the contour length of unfolded domains (Zhu & Asaro 2008), ${\rm\Delta}{\rm\Delta}x^{\ast }$ is the the difference between the activation length of the unfolding process and that of the refolding process (Zhu & Asaro 2008), $F_{1/2}$ is the the force corresponding to the state when half of the domains are unfolded (Zhu & Asaro 2008) and ${\it\mu}_{s}$ is the initial shear modulus of the cytoskeleton (Peterson, Strey & Sackmann 1992). A spectrin consists of 19 domains in our model.

Figure 3

Figure 2.(a) The shape, (b) the area deformation of the skeleton and (c) the shear deformation of the skeleton at the resting configuration. Here, $V_{s}/V_{0}=0.975$ and $c_{0}=8$.

Figure 4

Figure 3.Snapshots of the cell in tank treading motion in (a) mode 1 and (b) mode 2. Two beads are attached to the membrane to mark the dimple elements. Here, $Ca=3.51$ and ${\it\Lambda}=1.0$. Mode 1 is achieved with ${\it\theta}_{0}=0^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$. Mode 2 is achieved with ${\it\theta}_{0}=90^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$.

Figure 5

Figure 4.Time histories of the inclination angle ${\it\Psi}$ (a) and Taylor deformation parameter $D_{xy}$ (b) of mode 1 (solid lines) and mode 2 (dashed lines). Here, $Ca=3.51$ and ${\it\lambda}=1.0$. Mode 1 is achieved with ${\it\theta}_{0}=0^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$. Mode 2 is achieved with ${\it\theta}_{0}=90^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$. The reference state is biconcave.

Figure 6

Figure 5.A transition event from mode 1 to mode 2. (a) Trajectories of the dimple elements in the $z$ direction. (b) The breathing motion characterized by variations of the Taylor parameter $D_{xy}(t)$. (c) The swinging motion ${\it\Psi}(t)$. (d) The kayaking motion ${\it\phi}(t)$. Here, $Ca=2.34$, ${\it\Lambda}=1.0$ and ${\it\theta}_{0}=0^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$. The reference state is biconcave.

Figure 7

Figure 6.A transition event from mode 2 to mode 1. (a) Trajectories of the dimple elements in the $z$ direction. (b) The breathing motion characterized by variations of the Taylor parameter $D_{xy}(t)$. (c) The swinging motion ${\it\Psi}(t)$. (d) The kayaking motion ${\it\phi}(t)$. Here, $Ca=3.90$, ${\it\Lambda}=0.5$ and ${\it\theta}_{0}=90^{\circ }\rightarrow {\it\phi}_{0}=0^{\circ }\rightarrow {\it\psi}_{0}=0^{\circ }$. The reference state is biconcave.

Figure 8

Figure 7.Stability diagrams in the ${\it\Lambda}{-}Ca$ plane of (a) an RBC with biconcave reference state and (b) an RBC with spheroidal reference state ($V_{s}/V_{0}=0.975$, $c_{0}=8$).

Figure 9

Figure 8.Stability diagram in the $V_{s}/V_{0}\sim Ca$ plane for an RBC with spheroidal reference state. Here, ${\it\Lambda}=0.75$ and $c_{0}=8$.

Figure 10

Figure 9.Schematic illustration of the membrane contour within the $x{-}y$ plane during tank treading motions for (a) mode 1 and (b) mode 2 over half a tank treading period (the time increases sequentially from $t_{1}$ to $t_{5}$). The dimple elements are shown as circles and the rim elements are shown as bullets.

Figure 11

Figure 10.Distributions of the strain energy density $E$ (normalized by ${\it\eta}_{1}\dot{{\it\gamma}}R$) along the contour within the $x{-}y$ plane in the membrane as the cell switches from (a) mode 1 (within half of a breathing/swinging period starting at $\dot{{\it\gamma}}t=67.5$) to (b) mode 2 (at $\dot{{\it\gamma}}t=1440$). This is the same case as defined in figure 5, which lies in region III in the stability diagrams. The curves in (a) show four evenly distributed snapshots over time (the variations are caused by the breathing/swinging oscillations in mode 1). In mode 2 the cell shape and configuration remain almost steady so that there is no variation of the strain energy density over time. The inset in (a) describes the definition of the angle ${\it\Theta}$.

Figure 12

Figure 11.Distributions of the strain energy density $E$ (normalized by ${\it\eta}_{1}\dot{{\it\gamma}}R$) along the contour within the $x{-}y$ plane in the membrane as the cell switches from (a) mode 2 (at $\dot{{\it\gamma}}t=125$) to (b) mode 1 (within half of a breathing/swinging period starting at $\dot{{\it\gamma}}t=2000$). This is the same case as defined in figure 6, which lies in region I in the stability diagrams.