Numerical-experimental observation of shape bistability of red blood cells flowing in a microchannel

Red blood cells (RBCs) are the major constituent of mammalian blood and therefore determine the majority of its flow properties. One of the most amazing features of RBCs is their deformability, allowing them to squeeze through channels with diameters much smaller than their own equilibrium size Freund (2013); Picot et al. (2015); Salehyar and Zhu (2016). Another consequence of their deformability is the wide range of stationary and non-stationary shapes assumed by the RBCs in microchannel flows with dimensions similar to or slightly larger than the RBC equilibrium radius Fedosov et al. (2014); Aouane et al. (2014); Tahiri et al. (2013). Understanding and being able to predict these shapes is of high importance for a variety of reasons. From a fundamental point of view, it serves as the foundation in a bottom-up approach to understand the properties of red blood cell suspensions which are chiefly determined by single particle behavior Vitkova et al. (2008); Fedosov et al. (2011); Krüger et al. (2013); Thiébaud et al. (2014); Katanov et al. (2015); Lanotte et al. (2016). From an applied perspective, a series of recent investigations have devised promising approaches for sorting cells based on their mechanical properties either in lateral displacement devices Henry et al. (2016) or using high-speed video microscopy Otto et al. (2015). Finally, knowledge of the precise cell shape is also essential for accurately measuring geometric properties of cells Merola et al. (2017).

The most frequently observed shapes of RBCs in microchannel flows are the so-called “croissant” and “slipper” shapes. Examples are depicted in figure 1. Some researchers refer to croissants also as parachutes, although here we prefer the term croissant since our shapes are not perfectly rotationally symmetric (similar to the ones found by Farutin and Misbah (2014)). Probably one of the earliest experimental study on isolated red blood cells in flow was performed by Gaehtgens et al. Gaehtgens et al. (1980), where slippers as well as parachutes have been found depending on the diameter of the cylindrical channel. Suzuki et al. Suzuki et al. (1996) presented a phase diagram of parachutes and slippers as a function of velocity and confinement in a cylindrical tube. Slippers dominated at smaller diameters and higher velocities. Secomb et al. Secomb et al. (2007) compared experiments with 2D simulations in cylindrical channels of $8\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter for a cell velocity of approximately $1.25\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$. Furthermore, two other publications Faivre (2006); Abkarian et al. (2008) considered the flow of RBCs at very low viscosity ratios of $\lambda\lesssim 0.27$. They presented a phase diagram showing parachutes and slippers, where the velocity was varied in the very high regime of $10$ to $170\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$. Tomaiuolo et al. Tomaiuolo et al. (2009) found parachutes at smaller and slippers at higher velocities in cylindrical channels of $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter. A subsequent study Tomaiuolo and Guido (2011) as well as Prado et al. (2015) considered the transient during start-up of the flow. Cluitmans et al. (2014) detected croissants at lower ($\lesssim$5\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$) and slippers at higher velocities ($\gtrsim$10\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$) in rectangular channels with widths $\leqslant$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. Moreover, Quint et al. (2017) found a stable slipper and a metastable croissant at the same set of parameters in a wider channel of $$25\text{\,}\mathrm{\SIUnitSymbolMicro m}$\times$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. Other publications presenting experiments in channel flow also touch the subject of RBC shapes but focus on other aspects such as the methodology Hochmuth et al. (1970); Seshadri et al. (1970); Zharov et al. (2006); Tomaiuolo et al. (2007); Guido and Tomaiuolo (2009); Gorthi and Schonbrun (2012); Lanotte et al. (2014); Tomaiuolo et al. (2016), dense suspensions and cell interactions Claver´ıa et al. (2016); Guest et al. (1963); Skalak and Branemark (1969); Gaehtgens et al. (1980); Kubota et al. (1996); Abkarian et al. (2008); Tomaiuolo et al. (2012); Wagner et al. (2013); Brust et al. (2014); Tomaiuolo et al. (2016) or use vastly larger channel diameters Goldsmith and Marlow (1972); Lanotte et al. (2016).

Numerical simulations and semi-analytical calculations of isolated particles in microchannels mostly studied axisymmetric RBCs Secomb et al. (1986); Secomb (1987); Secomb et al. (2001) or 2D vesicles Secomb and Skalak (1982); Kaoui et al. (2009a, 2011, 2012); Shi et al. (2012); Tahiri et al. (2013); Lázaro et al. (2014); Aouane et al. (2014). The numerical work by Aouane et al. (2014), for example, identified a large amount of dynamics including deterministic chaos. The first full 3D simulation of single cells with a realistic RBC model (but with a ratio of inner to outer viscosity of $\lambda=1$) was conducted by Noguchi and Gompper Noguchi and Gompper (2005) who used a cylindrical tube with a diameter of $9.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$. They found the typical discocyte shape below and parachutes above a critical velocity which depends on the elastic parameters. A subsequent study by the same group additionally explored this threshold as a function of confinement McWhirter et al. (2011). Moreover, Fedosov et al. Fedosov et al. (2014) presented very detailed phase-diagrams where the velocity and confinement was varied for three different sets of elastic moduli and a viscosity ratio of $\lambda=1$. They observed four distinct regions where snaking, tumbling, slippers and parachutes occurred. Recently, Ye et al. (2017) considered the shapes of an RBC with $\lambda=1$ in rectangular microchannels (with width $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and aspect ratios $1$ to $2$) for the three cell velocities $4$, $20$ and $100\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ and observation times up to $\approx$0.03\text{\,}\mathrm{s}$$. Snapshots after this short initial transient showed parachutes or slightly slipper-like shapes.

Bistability, i.e. the observation of two different stable shapes depending on the initial condition but at otherwise identical system parameters, was barely considered so far. It was observed only numerically for simpler situations such as close-to-spherical vesicles in unbounded Poiseuille flow Farutin and Misbah (2014) or near a single wall Kaoui et al. (2009b), for a 2D RBC model in bounded Poiseuille flow Secomb et al. (2007), for the initial transient of a red blood cell in a rectangular channel Ye et al. (2017) or for simple shear flows Cordasco et al. (2014); Peng et al. (2014); Sinha and Graham (2015); Lanotte et al. (2016). No systematic experimental investigations exist for cells flowing in microchannels. Moreover, the 3D simulations and experimental investigations that were mentioned above and that consider the RBC shapes in microchannels in more detail all used a viscosity ratio of $\lambda\leqslant 1$, although 2D simulations showed that choosing a physiologically more realistic value of $\lambda\approx 5$ Cokelet and Meiselman (1968) can significantly affect RBC dynamics Kaoui et al. (2012); Tahiri et al. (2013).

Here we present a detailed systematic experimental-numerical study on the steady-state shape of isolated red blood cells in a rectangular microchannel. We use the physiological viscosity ratio of $\lambda=5$ appropriate for healthy human red blood cells in the microcirculation Cokelet and Meiselman (1968). The initial position is varied in the simulations directly, while experimentally we determine it via measurements at the channel entrance. Our central finding is that the initial starting position of the RBC has a decisive influence on the final steady-state shape of the red blood cell.

We begin by outlining our experimental and numerical methods in section II. Afterwards, the results from our experiments (section III) and simulations (section IV) are presented, while section V is dedicated to their detailed comparison. Finally, we conclude our work in section VI.

We classify cells in the experiments either as croissants, slippers or “other” not uniquely identifiable or completely different shapes. Typical slipper and croissant shapes are shown in the photographs (b) and (e) of figure 1. See the supplementary information (SI) for a collection of all images.

To systematically investigate the occurrence of the different shapes, we vary the imposed pressure drops from 20 to $1000\text{\,}\mathrm{m}\mathrm{b}\mathrm{a}\mathrm{r}$. The corresponding cell velocities range from $0.14\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ to $10.6\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$, covering the whole physiological range in microchannels Pries et al. (1995); Popel and Johnson (2005); Baskurt et al. (2011). We consider the cells $10\text{\,}\mathrm{m}\mathrm{m}$ away from the channel entrance where most of the cells reached a steady state Claver´ıa et al. (2016). Figure 3(a) depicts the fraction of observed shapes as a function of the measured cell velocities, constituting our central result from the experiments. This distribution was obtained by considering typically more than 100 cells per imposed pressure drop. The average velocities were computed by averaging over all cells at a certain pressure drop, with the horizontal error bars showing the corresponding standard deviations $\sigma_{u}$ in cell velocity. Not all velocities are the same because croissants and slippers have different velocities at otherwise identical flow conditions Quint et al. (2017), and because of the natural variations of cell properties such as elasticity and size, as also noted by Tomaiuolo et al. (2009). See the supplementary information for more details. Considering figure 3(a), high velocities obviously favor slippers while croissants are the most prominent for medium velocities. A pronounced peak exists from around $1$ to $2\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$. Very small velocities produce mostly shapes that fall outside our simple two-state classification.

Figure 3(b) illustrates the corresponding estimated probability density function of the center of mass $y$-position of the cells at the various pressure drops. This estimate was obtained from the measured $y$-positions by using the kernel density estimator as implemented in MATLAB R2017a (ksdensity) with a support of $[-6,6]\,$\mathrm{\SIUnitSymbolMicro m}$$ and otherwise default settings. Thus, croissants and “others” occurring at lower velocities are centered in the channel, while slippers occurring at high velocities show a pronounced off-centered position. The assumed shapes therefore imply a certain $y$-position within the channel with slippers being off-centered and croissants centered. This is confirmed when analyzing the offset distribution separately for each shape class as shown in the supplementary information.

From figure 3(a) it is tempting to conclude that the flow velocity is the major parameter that determines the RBC shape with low velocities favoring centered and high velocities favoring off-centered flow positions. However, looking at the cell positions near the channel entrance (figure 4) we find that already upon entering the channel RBCs are not homogeneously distributed. At low velocities we observe a clear bias towards a centered initial position, with the distribution becoming approximately homogeneous only at the highest measured velocities. These experimental observations allow two distinct parameters as the reason for the dominance of the slipper shapes at high velocities: either the higher flow velocity itself or the more off-centered entry into the channel. To disentangle these two possibilities we now present numerical simulations whose geometry directly corresponds to the experimental setup.

We numerically study the behavior of a single RBC in a rectangular microchannel by varying the imposed flow velocity, the initial shape and the initial offset $r_{\mathrm{init}}$ from the centerline of the tube (see section II.2). After starting the flow, we wait until the RBC reaches the steady state where the shape as well as the radial position does no longer change, or alternatively until periodic motion is observed.

In the majority of cases, we observe two different states: A croissant shape (which moves as a rigid body, figure 1(c)) and a slipper shape (figure 1(f)). The latter exhibits tank-treading (TT) and oscillatory contractions similar to the slippers seen by Fedosov et al. (2014) (see the SI for a movie and the insets in figure 5). Tank-treading refers to the motion of the membrane around a (more or less) static shape. Note that perfectly axisymmetric parachutes are suppressed by the rectangular channel flow, contrary to the situation for cylindrical tubes Fedosov et al. (2014) or unbounded Poiseuille flows Farutin and Misbah (2014).

To start the systematic study, we take a red blood cell that is initially in the typical discocyte shape with its rotation axis aligned along the tube’s axis (cf. fig. 2). We then vary the radial offset $r_{\mathrm{init}}$ from the center line as described in section II.2 and record the final radial position as well as the shape. The mean of the radial position is extracted by a temporal average once the cell is in the steady state (see the supplementary information for more details). Figure 5 shows the result for a cell velocity of $\approx$6.5\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$. A single sharp transition at $r_{\mathrm{init}}\approx$0.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ from centered croissants to off-centered slippers is observed. The final position of the slippers is mostly offset only along the wider width of the channel ($y$-direction), but not along the smaller height ($z$-direction). Hence we find pronounced bistability: The result is significantly determined by the initial condition and two different shapes coexist. This is consistent with the 2D simulations by Secomb et al. Secomb et al. (2007) and Tahiri et al. Tahiri et al. (2013). It also agrees qualitatively with observations by Farutin and Misbah for 3D simulations of vesicles in unbounded Poiseuille flow Farutin and Misbah (2014).

To study the bistability in more detail, we vary the imposed flow velocity as well as the initial offset $r_{\mathrm{init}}$ and characterize the behavior in the steady state. This yields the shape phase diagram depicted in figure 6(a). The cell velocity is extracted in the steady state via a temporal average. For slippers the velocity varies periodically (similar to the radial position): the minimum and maximum in one period is indicated by the horizontal error bars. Overall, the mean cell velocity $u$ ranges from $0.132\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ to $10.4\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$, matching with the experimentally covered range. The corresponding shear capillary number $\mathrm{Ca}_{\mathrm{S}}:=\mu u/\kappa_{\mathrm{S}}$ varies therefore in the interval $\mathrm{Ca}_{\mathrm{S}}\in[0.0317,2.50]$, while the bending capillary number $\mathrm{Ca}_{\mathrm{B}}:=\mu uR_{\mathrm{RBC}}^{2}/\kappa_{\mathrm{B}}$ lies in the range $\mathrm{Ca}_{\mathrm{B}}\in[8.45,666]$. The reddish area illustrates the approximated region where croissants exist. Furthermore, there is a maximal initial offset $r_{\mathrm{init}}$ above which overlapping with the vessel wall would occur.

The shape phase diagram in figure 6(a) (together with (b) and (c) explained below) constitutes our main result from the simulations. Starting near the channel center (in the reddish region) results in croissants, whereas higher initial offsets lead to slippers. The transition is found to be sharp, and depends significantly on the velocity. Croissants are the only stable steady state in a small region ranging from around $2$ to $3\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$, independently of the initial radial position. Smaller and larger velocities tend to favor slippers. Stable croissants do not appear below $0.25\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$. In the case of the slippers, the final periodic state is usually reached after roughly $2$ to $10\text{\,}\mathrm{s}$. In contrast, the final croissant state is sometimes achieved only after more than $10$ to $30\text{\,}\mathrm{s}$, possibly after an intermediate slipper state that can last several seconds (see figure S4 and the movie in the supplementary information). Hence, shapes observed after less than one second often turn out to be transient, contrary to the interpretation of Ye et al. (2017) but in agreement with Prado et al. (2015).

Considering our results in figure 6(a) in more detail, we find that two different types of croissants and slippers are possible. On the one hand, at very low velocities ($\lesssim$0.7\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$) the slippers no longer exhibit tank-treading motion of the membrane and instead show tumbling behavior: The cell rotates around the $z$-axis while approximately preserving its shape (similar to a rigid-body, see the SI for a movie). The difference compared to the tumbling motion observed by Fedosov et al. (2014) is that the cell still exhibits a clear slipper-like instead of a proper discocyte shape. Hence, we classify this mode still as slipper. On the other hand, at very high velocities ($\gtrsim$7\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$) slightly asymmetric shapes strongly reminiscent of croissants with a distinct tank-treading motion can sometimes be observed (see the inset in figure 6(c) for an example). As the shape itself is very close to a croissant, we will nevertheless consider it to be a croissant below.

A natural question that occurs in light of the profound bistability is the influence of other initial shapes on the result. To this end, we consider a typical croissant as well as a typical slipper as the starting shape. Both were obtained from previous simulations that started with the discocyte form and are depicted in the supplementary information. We once again construct the shape phase diagram as before and display the results in figures 6(b) and (c). Note that the different starting shapes admit a larger initial radial position $r_{\mathrm{init}}$ of the centroid. In short, starting with a croissant favors croissants in the steady state (the reddish area is larger than in figure 6(a)). For slippers it is the other way around: Starting with a slipper tends to produce more slippers (reddish area smaller than in figure 6(a)). Despite this, the croissant-only region from around $2$ to $3\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ still exists unscathed. Overall, only two qualitative differences occur between the phase diagrams of different initial shapes, both at lower velocity when starting with the croissant shape (figure 6(b)): First, stable croissants emerge at very low velocities ($\lesssim$0.7\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$) and second, the croissant-only peak exhibits a “protrusion” into the slipper space. This observation suggests that slippers and croissants can be stable below $2\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ for most $r_{\mathrm{init}}$ values, although in some cases a very precise croissant configuration is required in order to actually get a croissant in the steady state.

Another interesting aspect concerns the radial positions of the centroids in the final steady states. The average values are obtained by computing the temporal average in the steady state first for each simulation, and then combining the results for identical shapes via a weighted arithmetic mean. We use the observation time in the steady state as the weight. This procedure leads to figure 7(a). Obviously, the final radial positions are independent of the initial starting shape, i.e. a particular steady state shape at a certain velocity is always located at the same position. Furthermore, non-tank-treading croissants are always almost centered, with only minor deviations away from zero. These slight deviations in the range from 2 to $4\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ are mainly due to some croissants exhibiting minuscule periodic shape deformations. Moreover, the centroids of tank-treading croissants occurring at velocities $\gtrsim$8\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$ are located near but not directly in the center. Their slight off-centered position is a result of their asymmetry.

In contrast to croissants, slippers are located $0.8$ to $1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ away from the channel’s axis. The minimum position is attained for velocities near the border of the croissant-only region in the phase diagram (at around $2$ and $3\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$, compare figure 6). Above, the off-center position increases and seems to converge to a value of around $1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The reason for this increase is that slippers become more elongated and thinner at higher velocities (up to a certain degree), as shown in figure 7(b) and also observed in previous experiments Tomaiuolo et al. (2009). Thus, they effectively become smaller in the radial direction and their centroids can move closer to the wall. We note that the distance between the wall and the upper side of the slipper approximately remains the same for all velocities. This also hints at that the “optimal” off-center position for the slippers is more than $1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ away from the center, and that this particular value is due to the smallness of the channel.

To summarize, we have performed in-vitro experiments and 3D simulations of healthy red blood cells flowing in a microchannel. The viscosity ratio was approximately $5$ and the flow velocities ranged from around $0.1\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ to $10\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ in both methodologies, corresponding to the typical conditions prevailing in the microvascular system. We found that both the flow velocity as well as the initial starting configuration (offset from channel center, shape) have a major impact on the final steady state of the cells. Using three different starting shapes (discocyte, croissant, slipper), we constructed the corresponding phase diagrams via simulations. In most cases the cells assumed one out of two different forms: either a centered croissant or an off-centered slipper. Interestingly, for most velocities bistability, i.e. a dependence of the final shape on the initial position, was observed. Only in a small range of velocities (at around $\approx$1\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$) was the final shape found to be always a croissant. The experimental diagram showed very good agreement with the numerical result, especially when considering the simulations that used the rather asymmetric slipper as starting shape.

We thus conclude that the employed numerical RBC model can sensibly describe the cell behavior in the presented setup. Moreover, since we used physiological viscosity ratios and flow velocities, we speculate that croissants and slippers can occur in the microvasculature at the same set of system parameters not just as transients but rather that both are states which are intrinsically assumed by the cells. Our results are important for applications where the cells should be in a specific state (e.g. in lab-on-a-chip devices) and allow for a comprehensive validation of numerical models.

There are no conflicts to declare.