The buckling instability of aggregating red blood cells

TL;DR

This study investigates the buckling instability in aggregating red blood cells, revealing how contact zone deformation occurs and identifying specific cell configurations through experiments and simulations, with implications for blood flow and tissue mechanics.

Contribution

It provides the first detailed experimental and theoretical analysis of RBC aggregation morphology and the buckling instability at the contact zones.

Findings

Identified three RBC contact zone shapes: parachute, male-female, and sigmoid.

Demonstrated nonlinear increase in deformation with interaction energy.

Predicted a buckling bifurcation leading to contact zone deformation.

Abstract

Plasma proteins such as fibrinogen induce the aggregation of red blood cells (RBC) into rouleaux, which are responsible for the pronounced shear thinning behavior of blood, control the erythro- cyte sedimentation rate (ESR) a common hematological test and are involved in many situations of physiological relevance such as structuration of blood in the microcirculation or clot formation in pathological situations. Confocal microscopy is used to characterize the shape of RBCs within rouleaux at equilibrium as a function of macromolecular concentration, revealing the diversity of contact zone morphology. Three different configurations that have only been partly predicted before are identified, namely parachute, male-female and sigmoid shapes, and quantitatively recovered by numerical simulations. A detailed experimental and theoretical analysis of clusters of two cells shows that the deformation increases nonlinearly with the interaction energy. Models indicate a forward bifurcation in which the contacting membrane undergoes a buckling instability from a flat to a de- formed contact zone at a critical value of the interaction energy. These results are not only relevant for the understanding of the morphology and stability of RBC aggregates, but also for a whole class of interacting soft deformable objects such as vesicles, capsules or cells in tissues.

Equations

$$y(x)=A{L}_1\sin \frac{2\pi x}{L_1}$$

$$E_i=-{ϵ}_{ad}{L}_1\left(1+A^2{\pi }^2-\frac{3}{4}{\pi }^4{A}^4\right).$$

$$E_{b1}$$

$$E_{b2}$$

$${\beta }_c=\frac{8{\pi }^2\left(9\sqrt{1-\tau }-8\right)}{-2\sqrt{1-\tau }\tau +3\tau +3\sqrt{1-\tau }-3}$$

$$E_B=\oint \frac{{\kappa }_B}{2}(c-c_0{)}^{2}dl$$

$$\mathbf{f}(\mathbf{X})={\kappa }_B[\frac{{\partial }^{2}c}{\partial l^2}+\frac{c^3}{2}]\mathbf{n}-c\zeta \mathbf{n}+\frac{\partial \zeta }{\partial l}\mathbf{t}$$

$$\varphi (r)=4ϵ[(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]$$

$${\mathbf{f}}^{\varphi }(\mathbf{X})=-{\int }_{{\sum }_{j\ne i}\partial {\Omega }_j}\frac{\partial \varphi (r)}{\partial r}\frac{\mathbf{r}}{r}dl(\mathbf{Y})$$

$$f(d)=\frac{3\pi ϵ}{64}{\left(\frac{\sigma }{d}\right)}^6\left(160-231{\left(\frac{\sigma }{d}\right)}^6\right)$$

$$\begin{gathered} {ϵ}_{ad}=-{\int }_{-\infty }^{infty}\varphi \left(\sqrt{x^2+d_0^2}\right)dx \\ =\frac{48\times 5^{5/6}\times 6^{1/6}}{11\times 77^{5/6}}\pi ϵ\sigma \\ \simeq 1.89273ϵ\sigma \end{gathered}$$

$$\begin{gathered} \mathbf{u}({\mathbf{X}}_0)= \\ \frac{1}{2\pi {\mu }_1(1+\lambda )}\underset{{\sum }_i\partial {\Omega }_i}{\int }\mathbf{G}(\mathbf{X},{\mathbf{X}}_0)[\mathbf{f}(\mathbf{X})+{\mathbf{f}}^{\varphi }(\mathbf{X})]dl(\mathbf{X}) \\ +\frac{(1-\lambda )}{2\pi (1+\lambda )}\underset{{\sum }_i\partial {\Omega }_i}{\int }\mathbf{u}(\mathbf{X})\cdot \mathbf{T}(\mathbf{X},{\mathbf{X}}_0)\cdot \mathbf{n}(\mathbf{X})dl(\mathbf{X}) \end{gathered}$$

$$y(x)=A{L}_1\sin \frac{2\pi x}{L_1}$$

$$L_2={\int }_{-L_1/2}^{L_1/2}\left(1+2{\left(\pi A\right)}^2{\cos }^2\frac{2\pi x}{L_1}\right)dx=L_1\left(1+A^2{\pi }^2-\frac{3}{4}{\pi }^4{A}^4\right)$$

$$E_i=-{ϵ}_{ad}{L}_1\left(1+A^2{\pi }^2-\frac{3}{4}{\pi }^4{A}^4\right).$$

$$E_{b1}=\frac{\kappa }{R^2}2\pi R=\frac{2\pi \kappa }{R}$$

$$C(x)=-\frac{y^{\prime \prime }(x)}{(1+(y^{\prime }(x){)}^2{)}^{3/2}}.$$

$$E_{b2}$$

$$=\kappa {\int }_{-L_1/2}^{L_1/2}\frac{16{\pi }^4{A}^2{\sin }^2\left(\frac{2\pi x}{L_1}\right)}{L_1^2{\left(4{\pi }^2{A}^2{\cos }^2\left(\frac{2\pi x}{L_1}\right)+1\right)}^{5/2}}dx$$

$$\simeq \kappa \frac{8{\pi }^4{A}^2-20{\pi }^6{A}^4}{L_1}$$

$$L$$

$$S$$

$$1$$

$$\tau /4\pi$$

$$R$$

$$L_1$$

$$E$$

$$+A^2\left(-\frac{1}{2}{\pi }^2\beta \sqrt{1-\tau }+\frac{{\pi }^2\beta }{4}+\frac{16{\pi }^4}{\sqrt{1-\tau }}+\frac{2{\pi }^4}{\sqrt{1-\tau }-1}\right)$$

$$\displaystyle+A^{4}\biggl{(}-\frac{5\pi^{4}\beta\tau}{16\sqrt{1-\tau}}+\frac{\pi^{4}\beta}{4\sqrt{1-\tau}}+\frac{\pi^{4}\beta}{16}-\frac{7\pi^{6}\tau}{2\left((\sqrt{1-\tau}-1\right)^{3}\sqrt{1-\tau}}$$

$$\displaystyle\hskip 28.45274pt+\frac{2\pi^{6}\tau}{\left(\sqrt{1-\tau}-1\right)^{3}}\frac{3\pi^{6}}{\left(\sqrt{1-\tau}-1\right)^{3}\sqrt{1-\tau}}-\frac{40\pi^{6}}{\sqrt{1-\tau}}+\frac{8\pi^{6}}{1-\tau}-\frac{3\pi^{6}}{\left(\sqrt{1-\tau}-1\right)^{3}}\biggr{)}$$

$$A=\{\begin{array}{l}0\text{ for }\beta <{\beta }_c\\ \pm \frac{\sqrt{\beta \left(4\sqrt{1-\tau }-2\right)+16{\pi }^2\left(\frac{1}{1-\sqrt{1-\tau }}-\frac{8}{\sqrt{1-\tau }}\right)}}{\pi \sqrt{\frac{\beta \tau \left(-5\sqrt{1-\tau }\tau -\tau +4\sqrt{1-\tau }+1\right)+8{\pi }^2\left(5\left(\sqrt{1-\tau }+1\right)+\left(11-84\sqrt{1-\tau }\right)\tau \right)}{(1-\tau )\tau }}}\text{ for }\beta \ge {\beta }_c\end{array}$$

$${\beta }_c=\frac{8{\pi }^2\left(9\sqrt{1-\tau }-8\right)}{-2\sqrt{1-\tau }\tau +3\tau +3\sqrt{1-\tau }-3}$$

Full text

Introduction

Red blood cells (RBCs) have been a model system and an inspiration for the biophysics of the membrane for decades[1]. Among the numerous works in the literature, many studies have been devoted to adhesion of cells or vesicles and capsules on flat and solid surfaces [2, 3], where, in the case of elastic capsules only, a buckling instability may occur due to the interplay between bending and stretching energy[4]. Buckling has also been observed in freestanding capsules due to osmotic pressure[5, 6] or shear forces by external flow [7], and in vesicles due to asymmetric lipid distributions [8] or shear forces by external flow [9]. The influence of intercellular adhesion [10] is less known despite its relevance in dense environments such as blood at physiological hematocrit or even tissues. Theoretical studies on aggregation between two RBCs based on minimizing the free energy of the membrane[11, 12] indicate a significant change of the geometry of contact zones between cells when varying the interaction energy. Some of these shapes have been observed experimentally [13] but not quantified as a function of the dextran or fibrinogen concentrations. A variety of membrane shapes has been predicted theoretically for RBC doublets[12], depending mainly on the nondimensionalized reduced adhesion strength $\gamma$ and the reduced volume $\nu=3V\sqrt{4\pi}S^{-3/2}$ where $V$ and $S$ are the enclosed volume and surface area of the membrane, a parameter that is indeed physiologically related to the hydration state and age of RBCs and can vary over a significant range. Examples include male-female (convex-concave) shapes and more or less pronounced sigmoid shapes obtained by varying $\nu$ and $\gamma$ as depicted in Fig. 1(a). The 3D reconstructions of our experimental confocal images show RBC doublets that look very similar to the theoretical predictions, too (Fig.1(b)). Experimental images also show that the amplitude of the deformation of contact zones depends non monotonically on macromolecular (dextran) concentration.

There are direct physiological implications of aggregation since RBCs in blood in stasis are known to form aggregates in the form of so called rouleaux. This clustering process is reversible and typical hydrodynamic shear forces in physiological flow are usually sufficient to break rouleaux, at least in larger vessels. In the 1960’s, fibrinogen could be identified as one of the main plasma proteins causing RBC aggregation or clustering [14, 15, 16]. It is known that the number and size of aggregates increase with fibrinogen concentration in the physiological range (from approximately $3mg/ml$ for healthy adults up to $10mg/ml$ in acute inflammatory phases). The level of RBC aggregation is used in one of the most fundamental standard hematological blood tests, the erythrocyte sedimentation rate (ESR) which is carried out worldwide many thousand times each day. The more and the larger the aggregates, the faster the sedimentation of the RBCs: by simply measuring visually the sedimentation front in a standardized glass capillary one gets a robust and quick, if unspecific, indication on the inflammatory state of the patient. Of course, the tendency of the RBCs to form aggregates may also increase the risk of thrombosis and cardiovascular diseases, especially in combination with stenosis.

The molecular mechanisms of macromolecular induced RBC aggregation have been the subject of controversial studies [17, 18, 19, 20, 21, 22]. Two models have been proposed: One is based on the physical effect of depletion [23, 24, 25, 26, 27] and the other on physisorption or bridging [28, 29, 14, 30, 31]. In addition to the specific adhesion mechanism, the strength of RBC aggregation also depends on their physical properties, such as deformability [32], surface charge [33] and reduced volume [11, 12, 34]. It can also be induced by other macromolecules such as dextran, a widely used molecule in laboratory experiments or as a plasma expander and in veterinary medicine. AFM based single cell force spectroscopy, rheology and sedimentation of RBCs [27, 35, 36] as well as theoretical models [22] depict that an increasing dextran concentration leads first to an increase of the interaction strength among RBCs up to a specified macromolecule concentration. Beyond that maximum the interaction strength decreases again, and cell–cell aggregation vanishes. This leads to a characteristic bell–shaped adhesion-energy versus concentration curve. For the case of fibrinogen, it was found the interaction energy increases with the concentration, but data from the literature are typically limited to lower concentrations than for the case of dextran [27, 35, 36, 37].

We present here a systematic experimental study in order to analyze quantitatively the evolution of the contact zone as a function of aggregating molecule concentration, in a large domain of parameter space by using either dextran or fibrinogen. We analyze the case of RBC doublets as well as that of larger rouleaux. A systematic numerical investigation supports the experimental results. The numerical results show that in the case of two adhering cells, the interface remains flat up to a critical value of interaction energy above which it bifurcates into a deformed, buckled state. Our analytical model reveals a forward bifurcation with critical parameters close to those of the full numerical problem.

Analysis of RBC doublet shapes

Experiment

The analysis of the interfacial shape of two cells rouleaux, here called doublets, was performed by confocal imaging of doublets during free sedimentation, i.e. while not in contact with surfaces. Typical images are shown in Fig. 2 for doublets in dextran 70kDa solutions (see S. I. for the fibrinogen case). There is a significant evolution of the interfacial shape with polymer concentration: At low and high concentrations the interface is only weakly deformed while in the intermediate range where interaction energy is higher, the interface is strongly bent. Qualitatively, three different kinds of interfacial shapes were identified: parachute, male-female and sigmoid. In the sigmoid (S-shape interface), cells were mostly laterally displaced, with a non-monotonous curvature of the contact zone, while in other shapes the two cells had a common symmetry axis. The parachute shape was defined by at least one of the cells having a convex and a concave side. The male-female shape was defined when the interface had a bulge, with rather flat non-interacting large membrane parts. Cells with a flat interface were also classified under this category. Similar shapes were found for fibrinogen solutions. As can be seen in Fig. 2(a), all shapes can be found for all polymer concentrations, probably due to dispersity of cell properties in the sample, with a probability that depends on concentration.

A quantitative analysis of the distribution of the different interfacial shapes for the different concentrations of dextran is shown in Fig. 2(b), revealing an increase of the occurrence of sigmoid shapes in the intermediate range of polymer concentrations, i.e. at larger interaction energy favoring slightly larger contact areas. At low and high concentrations of dextran, doublets are mostly in the parachute (about 60 %) and male-female shapes. For fibrinogen (Fig. 2(c)), mostly parachute shapes were observed in the physiological range of the protein’s concentration. The case of echinocyte RBC doublets, obtained after sedimentation and contact with glass, was also studied (see S. I.).

Quantification of the interfacial shape deformation was done by fitting it with a heuristically chosen sinus function (see S.I.), whose amplitude was taken as a measure of the interfacial curvature (Fig. S2). A bell shape dependency of the deformation amplitude on polymer concentration is seen (Fig. 3), similar to the dependency found in AFM single cell force spectroscopy [27] or aggregation index (see S.I.). For the case of dextran, both doublets with discocytic and with echinocytic cells were analysed. The amplitude for discocytic cells can reach much larger values than echinocytes that have a close-to-spherical shape and probably a higher membrane stiffness [34].

Simulation and theory

To compute the interaction energies corresponding to the experimental doublets and gain further insight on the variability of shapes we performed numerical simulations (see S.I. for details of the method). Based on several successful studies (especially equilibrium shapes and shape dynamics under external flow) which have shown that 2D and 3D models capture the same essential features [38, 39, 40], we restrict ourselves here to a 2D model. In 2D the notion of in-plane shear elasticity looses its meaning and the vesicles and inextensible capsule model (often evoked to model RBC) are equivalent.

For reduced areas $\tau=4\pi S/L^{2}$ (where $S$ is the area of the 2D vesicle and $L$ its perimeter) in the range ($0.40,0.65$) with different combinations (cells with same or different reduced areas) and a membrane stiffness of $\kappa_{B}=4\times 10^{-19}$ J, we varied the interaction energy $\epsilon_{ad}$ in the range 0-8.23 $\mu$J/m2. The resulting interfacial shapes for different combinations of reduced area $\tau$, interaction energy $\epsilon_{ad}$ and curvature energy $\kappa$ (Fig. 4) were in qualitative agreement with the experimentally obtained shapes, namely male-female, S-shapes and parachute shapes.

Figure 4(a) shows the evolution of the deformation as a function of the interaction energy. As for the analysis of experimental data we fit the interfacial shape with a sine function. The mean amplitude of the fitted curve was compared to experimental results (Fig. 3), by heuristically choosing a correspondence between concentration and interaction energy: We assume a linear increase of the interaction energy with concentration up to 6.9 $\mu$J/m2 at 60 mg/ml, followed by a symmetric linear decrease, in order to get a good compromise between an agreement on the width and the maximum value of bell shaped curves. This piecewise linear variation of the interaction energy is a simplification which nevertheless yields a qualitative agreement for the overall behavior. Interestingly, the numerical results also show that the transition from a flat to a deformed interface only occurs at some finite interaction energy $\epsilon_{ad}$ (Fig. 4(b)). In principle, this method could allow the determination of the experimental interaction energy for a given concentration. However, the numerical results depend quite sensitively on the reduced areas $\tau$ of the two cells (Fig. 4(b)), a parameter that varies slightly between cells in experiments, even in the same sample. This dispersity is probably also the reason why the bell shape in Fig. 3(a) of the experimental data does not show this sharp transition.

The evolution of the amplitude of contact zone’s deformation in Fig. 4(b) suggests that a bifurcation akin to buckling occurs at a critical value of the interaction energy, triggering a change of shape from flat to S-shape. To get some insight into the involved mechanism, this behavior can actually be recovered by an analytical model through a minimization of the energy of the system. Let us suppose that the shape of an RBC in a symmetric doublet (in 2D) is approximated as shown in Fig. 12(a), with two semi-circular caps of radius $R$, a roughly straight portion where curvature is negligible (length $L_{1}$) and a contact zone with a sinusoidal shape of curvilinear length $L_{2}$ and amplitude $A$ described by the equation:

[TABLE]

for $-L_{1}/2<x<L_{1}/2$. The corresponding interaction energy (per unit length in the third dimension) is, at order 4 in $A$:

[TABLE]

where $\epsilon_{ad}$ is the interaction energy per unit area.

The curvature energy of the system is concentrated in the circular caps $E_{b1}$ and the deformed contact zone $E_{b2}$ given by:

[TABLE]

The total energy $E=E_{i}+E_{b1}+E_{b2}$ therefore depends on $A$, $R$ and $L_{1}$ which are related through the constant area $S$ and constant perimeter $L$ constraints.

In the following we define a dimensionless energy $E^{*}=EL/\kappa$, a rescaled interaction energy $\beta=\epsilon L^{2}/\kappa$, and rescale all lengths by $L$: $R*=R/L$ and so on, and introduce the reduced area $\tau=4\pi S/L^{2}$. The total energy $E^{*}$ can then be derived and expanded at order 4 in $A$ (see supplemental material). Figure 12(b) shows the energy due to deformation for $\tau=0.6$ and different interaction energies $\beta$. Below a critical value $\beta_{c}$, the minimum of energy corresponds to $A=0$ (flat interface), while for $\beta>\beta_{c}$, the energy is minimal for a finite value of $A=A_{eq}$ (see S. I.).

This defines the bifurcation threshold $\beta_{c}$ corresponding to $A_{eq}=0$, which is an increasing function of $\tau$ (Fig. 12(c)):

[TABLE]

The bifurcation diagram of the deformation amplitude (Fig. 12(d)) reveals a classical supercritical bifurcation, owing to the competition between interaction energy and curvature energy through the constant perimeter and area constraints. Orders of magnitude of the bifurcation threshold and amplitude of the deformation are in good agreement with experiments and numerical simulations. We have $\epsilon=\kappa\beta/L^{2}$. If we take $L\sim 24\mu$m and $\kappa\sim 4\times 10^{-19}$ J as in simulations, the critical interaction energy is $\epsilon_{c}\sim 2\mu$J/m2 for $\tau=0.6$ and $\epsilon_{c}\sim 3.2\mu$J/m2 for $\tau=0.65$, which is in reasonable agreement with the simulations of Fig. 4(b), while the amplitude of the deformation $A$ saturates at $\sim 0.13L=3.1\mu$m in this model, rather close to the value $\sim 2~{}\mu$m found in the simulations given the simplicity of the model.

Interestingly, the energy of the equilibrium configuration decreases linearly with $\beta$ as shown in Fig. 12(d), with no sign of the shape bifurcation due to the smallness of the non linear terms arising from the deformation in the curvature and interaction energy terms. This result can be of interest in the modeling of clustering dynamics in flow for instance, where the dissociation rate due to hydrodynamic stresses can be assumed to simply scale like $1/\epsilon_{ad}$ despite the possibly complex shapes of contact zones.

Analysis of membrane shapes in larger RBC clusters

A characterization of larger clusters was performed with rouleaux of at least 7 cells. The fast sedimentation and arbitrary orientation of the rouleaux only allowed us to image them when settled on the (BSA treated) cover slip. Since fibrinogen led to significant inter-connections between rouleaux which eventually led to a percolated gel like structure at 45% hematocrit, we only focused on rouleaux in linear configuration in dextran solutions (Fig. 6a). At 10mg/ml the cells in a rouleau were typically laterally shifted and we did not observe large clusters. At higher concentrations of dextran larger rouleaux were observed and as shown in Fig. 6, only the first three interfaces exhibited significant deformation (Fig. 6a) while the interfaces far from the tip were always in a male-female shape with very low deformation, mostly below the experimental resolution limit. Again, a large variety of shapes was observed for these first three interfaces (Fig. 6b).

As for doublets, these interfaces can be sorted into different categories and the distribution of shapes is shown separately for the first, second and third interface starting from the doublet tip in Fig. 7a. Unlike the case of doublets, we observed mostly male-female shapes and no parachute shapes at all, likely due to the constraints induced by the next attaching cell. For concentrations in the lower to medium-range, S-shape contact zones exist more frequently for the first contact zone than for the next ones, while male-female shapes are indeed prevalent everywhere except at very low dextran concentrations where RBCs tend to stay very close to their free discoid shape due to the low interaction energy, in agreement with the theoretical results on low-aggregation behavior mentioned in [11].

The amplitudes of the interfacial deformations (determined by the same method as for doublets) are shown in Fig. 7a. Again, a bell shape behavior for all three interfaces was observed. Interestingly, the amplitude of the first contact zone was very similar to the values measured for echinocyte doublets (see S.I.), too. Indeed, the aggregation into larger clusters gives less degree of freedom to RBCs due to higher constraints on the membrane, leading to a higher effective stiffness of RBCs.

Numerical simulations were performed for 7 cells and different interaction energies $\epsilon_{ad}$, a constant reduced area $\tau=0.65$ and membrane stiffness $\kappa_{B}=4\times 10^{-19}$ J (Fig. 8). The results agree qualitatively with the experimental results and the amplitudes of the interfacial deformations of the first and second interface were in good agreement with the experimental data using the same energy scale as for doublets (see Fig. 7). However, the amplitudes between experiment and theory differed significantly for the third shape. Here the variability of the experimentally observed shapes plays a role, as well as the blurring of the contact zone due to resolution in experiments, which tends to increase the amplitude of the fitted sine function.

Discussion

Red blood cells aggregate due to the presence of plasma proteins or biomimicking polymers. This is an example of adhesion between two very soft interfaces, different from the case of adhesion on a solid substrate, where we observe a buckling of the cell–cell interface. The nature of this instability is revealed by our numerical simulations and explained by an analytical model. The critical threshold of buckling and the final geometry of the aggregate depends strongly on the physical parameters, i.e. on the reduced volume of the cells. For large reduced volumes that correspond to more spherical cells, the critical interaction energy for buckling increases. In the case of doublets of two cells having identical reduced volumes close to that of a typical RBC, a sigmoid deformation of the interface takes place, while for two different reduced volumes male-female shapes can be obtained due to the asymmetry of properties that forces a different mode. We also characterized long linear clusters experimentally and numerically and have found that only the three interfaces at the two ends are considerably deformed. It would undoubtedly be interesting to check how these contact shapes obtained in quasi-static conditions evolve in pathological conditions, or under flow with the additional constraint of hydrodynamic stresses.

Methods

Human blood withdrawal from healthy donors as well as blood preparation and manipulation were performed according to regulations and protocols that were approved by the ethic commission of the Ärztekammer des Saarlandes (reference No 24/12). We obtained informed consent from the donors after the nature and possible consequences of the studies were explained.

RBC preparation

Venous blood of three healthy donors (one female, 26 years, blood group 0 negative; one male, 43 years, blood group 0 negative; one male, 27 years, blood group B positive) was drawn into conventional EDTA tubes (S-Monovette; Sarstedt, Nümbrecht, Germany) and all measurements were performed within 4 hours after the blood had been drawn. After washing the RBCs three times (704 g, 3 min) with phosphate-buffered solution (PBS, Life Technologies, Waltham, MA, USA) 0.1 $\mu$l CellMask (Life Technologies, Waltham, MA, USA) was added for fluorescence labeling to 200 $\mu$l RBCs in 1ml PBS for 10min. After washing two times (704 g, 3min) the supernatant was removed and the RBCs were resuspended in solutions of dextran (Dextran 70 from Leuconostoc mesenteroides, Sigma-Aldrich, St Louis, USA) or (washed) fibrinogen (fibrinogen from human plasma, Sigma-Aldrich, St Louis, USA) in PBS. The fibrinogen as delivered contains approximately 60% Fibrinogen, 25% sodium chloride and 15% of sodium citrate. In order to remove the salt the solutions were ultrafiltrated with a 10kDa filter (Vivaspin Turbo 4, Satorius AG, Goettingen, Germany) at 3500g in a centrifuge (Hermle Z326K, Hermle Labortechnik GmbH, Wehingen, Germany) for 10 hours by adding 2ml of PBS every hour. The remaining protein solution was freeze-dried and the proteins were resuspended in PBS. Fibrinogen was used in the measurements within 24 hours after freeze-drying and the osmolality of this solutions was verified to be in the order of PBS (+/- 5%).

Well preparation

Cylindrical wells of diameter 5 mm, height 1.6 mm ($\mu$-Slide 18 Well-Flat from Ibidi, Munich, Germany) were used as sample chambers. RBCs sitting on the bottom slide of the well changed to echinocytes due to the glass effect. For the measurements in which we wanted to preserve the discoyte shape of sedimented cells, the wells were covered with 1 mg/ml bovine serum albumine (BSA, Sigma-Aldrich, St Louis, USA) in PBS for 30 min at 37∘C. The BSA-PBS solution was removed and after a drying time of 10 min at 37∘C the RBC suspensions were added.

Imaging of the interfacial shape

All observations were performed with the confocal microscope described in the previous section at a temperature of 23 ∘C with an hematocrit of 0.2 % for the imaging of RBC doublets and 0.4 % for clusters larger than 6 RBCs. Images of doublets of RBC in their normal discoid shape could be recorded in free volume during the slow sedimentation process which prevented any contact with the glass surface of the cover slip. This procedure did allow to study the aggregates without any other influence but we should also note that it was also not possible to investigate the interfacial shape of sedimented cells in discocytic state. In this configuration, they were horizontally oriented with their interface and the RBC membrane was very thin. The optical 2-D sections showed a good signal only when the membrane was oriented along the axis of the point spread function of the laser focus. However, images of sedimented RBCs were also taken on the non-treated cover slips, i.e. for the case of echinocytes which were often oriented with their interface perpendicular to the cover slip, most likely due to their rounder shape. In total, 30 cells for discocytes and 15 for echinocytes were evaluated. For clusters that were larger than 6 cells, measurements could be only be performed on sedimented cells. In this case the cover slip was always treated with BSA to preserve the discocytic shape, and contact zones were always oriented perpendicular to the cover slip, thereby allowing a good characterization. The raw images were contrast enhanced and the position of the interfacial membrane contours was manually determined by visual inspection. The point cloud that indicated the interface was extracted with Matlab (The MathWorks GmbH, Ismaning, Germany) and stored into a one dimensional array to allow for a fitting procedure with Origin (OriginLab Corporation, Northampton, MA, USA). The interface was fitted with a heuristically chosen sinus function which yielded reasonable agreement. For sigmoid shapes the wave length of the sinus was left as a free parameter and for the male-female and parachute shape the wave length was limited to a range of $\pm 20\%$ around the double of the cell diameter.

3-D imaging

This accounts only for the images shown in Fig.1. In the preparation step for the wells CellTak (CellTak, BC Bioscience, San Jose, USA) was added as specified by the manufacturer in order to immobilize the cells on the BSA treated cover slip. A stack of 50 images was taken at a distance of 0.1 $\mu$m in z-direction, while the x-y area was 12 x 12 $\mu$m. A surface rendering with Imaris (Bitplane AG, Zuerich, Switzerland) allowed to look through the cell membranes attached to the bottom onto the contact zones.

Acknowledgements

The research leading to this results has received fundings from the German Research Foundation (DFG, SFB 1027), the Centre National de la Recherche Scientifique (CNRS) and the German French University (DFH/UFA). T. P. and C. M. acknowledge support from CNES and LabEx Tec21.

Author contributions statement

D. F. performed experiments and data analysis, O. A. ran numerical simulations, L. K. supervised the experiments, C.R. programmed the analysis software, S.S. compared 3D analysis, C. M. supervised numerical simulations, T. P. designed the research, made the analytical model, supervised data analysis, C. W. designed the research, supervised experiments and numerical simulations.

Supplementary Information

Microscopic aggregation index

A simple way to quantitatively characterize the clustering of the RBCs is based on the quantification of sedimented aggregates within a certain area on a microscope slide. Chien and Jan [33] defined a microscopic aggregation index (MAI) as the ratio of the total number of cells to the total number of ”particles” (individual cells and aggregates), representing the average number of cells per particle. Here, we investigate RBC aggregation by confocal microscopy and define the MAI as the the ratio of the number of aggregated cells at the bottom of the well after sedimentation to the total number of cells, i.e. the fraction of aggregated cells.

All observations were performed at a temperature of 23 ∘C with a hematocrit of 10 %. This concentration ensured an almost fully covered bottom surface after 5 min of sedimentation for most sample solutions, with an average coverage of $2.9\times 10^{4}$ cells/mm2. Imaging was performed using a confocal microscope (TSC SP5 II; Leica Microsystems, Mannheim, Germany) with a laser of 633nm as the excitation wavelength. RBCs were classified as part of a cluster if the membranes of two cells were optically indistinguishable over a distance of at least 4 pixels, which corresponds to 1 $\mu$m in our image with a size of 144 x 144 $\mu$m2. White-light microscopy was used as a control to exclude any possible influence of CellMask on aggregation (data not shown). Fibrinogen and dextran solutions were compared to PBS as negative control. Due to the increase of viscosity with increase in dextran or fibrinogen concentration, the sedimentation slows down, which leads to a lower total number of cells per image at higher viscosities in comparison to RBCs in PBS.

Figure 9(a) shows representative images for various concentrations of dextran and fibrinogen.

In dextran solutions, the number of rouleaux and the MAI exhibit a bell shaped curve with a maximum around a concentration of 50 mg/ml and almost no rouleaux at 120 mg/ml. (Fig. 9(b)). For fibrinogen, a monotonous increase of the MAI is seen up to a concentration of 12mg/ml, which corresponds to the physiological range. Both results are reminiscent of the dependencies of interaction energy, in-flow clustering and rheology on dextran concentration obtained in previous studies [27, 35]. Another qualitative observation is that the dextran and fibrinogen induced rouleaux look slightly different with more interconnections in the fibrinogen case.

Contact zone contour extraction

All measurements of contact zone shapes of RBC clusters were performed at 0.2% v/v hematrocrit. The procedure for extracting interfacial shapes for further analysis is exemplified in Fig. 10. The contact zones of RBC clusters are first isolated. The images are then thresholded and skeletonized to extract the interface shape on which a sine function can be fitted to measure the amplitude of deformation.

Echinocyte clusters

When the doublets were allowed to sediment on the BSA treated cover slip, the cells remained in discocyte state for several minutes but they were mostly oriented with their interface parallel to the coverslip which did prevent a characterization of the thin interfacial membrane that was then oriented perpendicular to the z-axis of the point spread function. This was different for the non-treated cover slips: due to the so called ”glass effect”, the RBCs transformed into the echinocyte state and we observed only male-female and sigmoid shapes (Fig. 11(a)). A quantitative analysis of the distribution of the different interfacial shapes for the different concentrations of dextran is shown in Fig. 11(b). For all different cases we observed an increase in the number of sigmoid interfaces in the intermediate range of polymer concentrations, i.e. at larger interaction strengths. The sigmoid shape allows for the largest contact zone, favorized by large interaction strength. At low and high concentrations of dextran we observe only few sigmoid shapes and more parachute than male-female shapes for the cells in discocytic state. With fibrinogen, mostly parachute shapes were observed. This is most likely related to the low (but physiological) concentrations of the protein. For cells in echinocyte state no parachute shape was observes which is to be expected in view of the spherical shape of a single echinocyte.

Modelling and numerical methods

Based on several successful studies (especially equilibrium shapes and shape dynamics under external flow) that have shown that 2D and 3D models capture the same essential features [38, 39, 40], have restricted ourselves here to a 2D model. In 2D the notion of in-plane shear elasticity looses its meaning and the vesicles and inextensible capsule model (often evoked to model RBC) are equivalent. Vesicles consist of closed bilayer membranes enclosing an inner fluid and suspended in an outer fluid. Using the typical reduced area (equivalent of the surface-to-volume ratio in three dimensions) of the RBC, the minimum energy shape of an unstressed vesicle is the same as for an RBC despite the lack of a cytoskeleton. Indeed, the stress free shape is dictated by the curvature energy and the rigidity of the cell. To describe the mechanical behavior of the vesicle, we employ (for each cell) the two dimensional Helfrich free energy[41, 42] (per unit length in the 3rd dimension):

[TABLE]

where $c$ is the mean curvature, $dl$ is the arclength element along the membrane contour, $\kappa_{B}$ is the membrane bending rigidity modulus, and $c_{0}$ is the spontaneous curvature that can be neglected in case of symmetrical membranes; actually in 2D the spontaneous curvature plays no role. In order to fulfill the condition of membrane inextensibility we supplemented the curvature energy with a term $\oint\zeta dl$, where $\zeta$ is a Lagrange multiplier that depends on the position along the membrane . The fluid inside and outside the membrane is described by the Stokes equations (the inertial effect is small and is neglected) with boundary conditions at the membrane and far away from the membrane. Far away from the membrane the flow vanishes (in the absence of any imposed external flow). At the membrane we require velocity continuity, and force balance. In addition, the velocity field along the membrane has to be divergence-free in order to enforce membrane incompressibility. This condition represents an implicit equation that allows one to determine the still unknown Lagrange multiplier $\zeta$; the details of how this is achieved in practice can be found in [43]. The hydrodynamics stress jump across is counterbalanced by the membrane force. This force is composed of the bending force obtained as a functional derivative of the Helfrich free energy supplied with a tension-like force to fulfill the condition of inextensibility of the membrane.

[TABLE]

where $\mathbf{n}$ and $\mathbf{t}$ are the normal and tangent unit vectors. $\mathbf{X}$ designates the vector position of a given point on a given cell; the curvature as well as the Lagrange multiplier are defined for each cell. In addition to this force we had to consider the force due to cell-cell interaction (due to the presence of proteins or dextran). For this purpose we refer to the work of Neu and Meiselman [44]. Neu and Meiselman have presented a model that takes into account the strong electrostatic repulsion due to the negative charge of RBC membranes, and the weak osmotic forces of attraction due to the depletion effect induced by macromolecules surrounding RBCs [44]. This model contains several parameters depending on the physiochemical properties and the concentration of the polymers. Often, the interaction energy between cells is described by a Morse potential (or alternatively a Mie type potential), having an exponential form, which has proven to adequately describe experimental data [45]. As is the case in many physical problems, a Morse potential can be well approximated by a Lennard-Jones potential, and we adopted this form in the present study for practical reasons. The interaction potential per unit (length of the 1st membrane $\times$ length of the 2nd membrane $\times$ length in the 3rd (invariant) dimension) is

[TABLE]

where $\epsilon$ is a microscopic parameter with dimension [energy/length3].

This potential is minimum at distance $r_{*}=\sigma 2^{1/6}$ and the minimal interaction potential is equal to $\phi_{*}=-\epsilon$.

The cell-cell force is given by

[TABLE]

$\mathbf{r}=\mathbf{X}-\mathbf{Y}$, $r=\|\mathbf{X}-\mathbf{Y}\|$, and $X$ and $Y$ are two position vectors belonging to the $ith$ and $jth$ membrane ($\partial\Omega$) respectively.

In this formulation, $\epsilon$ is a microscopic parameter of interaction between surface elements. Given the range of the interaction potential (typically a few times $\sigma$), an integration of the potential over planar parallel surfaces of extension much larger than this range must be performed in order to derive the macroscopic interaction energy per unit area $\epsilon_{ad}$.

The force per unit surface between two infinite parallel membranes at distance $d$ is:

[TABLE]

The zero force position is therefore $d_{0}=(231/160)^{1/6}\sigma$, which is smaller than $r_{*}$. The interaction energy per unit surface between infinite plates (i.e. whose extension is much larger than $d_{0}$) is then:

[TABLE]

This gives the relation between the microscopic parameters of the potential $\epsilon$ and $\sigma$ and the macroscopic interaction energy per unit area $\epsilon_{ad}$.

Recently the RBC-RBC adhesion energy has been quantified using the atomic force microscopy-based single cell force spectroscopy [27]. From these data, in principle, we can determine the adhesion energy $\epsilon_{ad}$ and equilibrium distance $d_{0}$ from which the two parameters entering the Lennard-Jones potential can be determined, namely the minimal energy $\epsilon$ and the range of the potential $\sigma$. However, here we chose to adapt the results from the numerical simulations, namely the amplitude of the interfacial deformation, directly to the experimental data and use this to deduce the interaction energy. In this procedure, due to numerical constraints, the zero force length (i.e. the length corresponding to minimal potential) was fixed in our simulations to 489 nm. This was probably significantly larger than the experimental value but still in the order of the optical resolution limit. In addition, note that as long as $\sigma$ is small compared to the lateral extension of the membrane, the product $\epsilon\sigma$ is the most quantity here. The Stokes equations were solved using the boundary integral technique as presented in  [46, 40]

[TABLE]

where $\mathbf{u}(\mathbf{X_{0}})$ is the velocity of a vector position $\mathbf{X_{0}}$ belonging to the membrane $\partial\Omega_{i}$, $\lambda$ is the viscosity contrast between the inner and outer fluids with the respective viscosities $\mu_{2}$ and $\mu_{1}$, and $\mathbf{G}$ and $\mathbf{T}$ are the Green’s functions associated with the velocity and the stress fields respectively.

To match the average values of a human RBC, the reduced area ($\tau$) of the vesicle typically lies in the range ($0.60,0.65$), the typical radius is $R_{0}=3\mu m$, the membrane stiffness is $\kappa_{B}=4\times 10^{-19}J$. Here we did not investigate any dynamic effects and we therefore fixed the viscosity ratio simply to unity, even if typically a factor of 5 was assumed. We considered doublets and sets of $7$ cells, separated by an initial surface-to-surface interdistance of 600 nm.

Analytical model

We provide here the details of the analytical model showing buckling of contact zones in symmetric RBC doublets. We suppose that the shape of an RBC in a symmetric doublet (in 2D) is approximated as shown in Fig. 12, with two semi-circular caps of radius $R$, a roughly straight portion where curvature is negligible (length $L_{1}$) and a contact zone with a sinusoidal shape of curvilinear length $L_{2}$ and amplitude $A$ described by the equation:

[TABLE]

for $-L_{1}/2<x<L_{1}/2$. The corresponding interaction energy of this contact zone is in expansion at order 4 in $A$:

[TABLE]

This corresponds to an interaction energy $E_{i}$:

[TABLE]

where $\epsilon_{ad}$ is the interaction energy per unit surface.

The curvature energy of the system is concentrated in the circular caps $E_{b1}$ and the deformed contact zone $E_{b2}$. The circular cap term is:

[TABLE]

The curvature of the contact zone is:

[TABLE]

This leads to the energy term:

[TABLE]

The total energy $E=E_{i}+E_{b1}+E_{b2}$ therefore depends on $A$, $R$ and $L_{1}$ which are related through the constant area $S$ and constant perimeter $L$ constraints

[TABLE]

In the following we define a dimensionless energy $E^{*}=EL/\kappa$, a rescaled interaction energy $\beta=\epsilon_{ad}L^{2}/\kappa$, and rescale all lengths by $L$: $R*=R/L$ and so on, and use the definition of the reduced area $\tau=4\pi S/L^{2}$.

The constant perimeter and area constraints can then be written as the following (dropping the * for simplicity):

[TABLE]

and can be used to eliminate $R$ and $L_{1}$ as in equations 22 and 23.

[TABLE]

The total energy $E^{*}$ can then be derived and expanded at order 4 in $A$ as in equation LABEL:energy, dropping the ∗ for simplicity.

[TABLE]

This energy has only a minimum for a finite, critical value of the rescaled interaction energy $\beta$.

If we solve $\partial E/\partial A=0$ to find the equilibrium configurations we get the two solutions whose expressions are given by 24, which defines the bifurcation threshold $\beta_{c}$ corresponding to $A_{eq}=0$ (Eq. 25)

[TABLE]

[TABLE]