Edited by William R. Schowalter, Princeton University, Princeton, NJ, and approved May 23, 2013 (received for review December 6, 2012)
Vasoocclusion crisis is a key hallmark of sickle cell anemia. Although early studies suggest that this crisis is caused by blockage of a single elongated cell, recent experiments have revealed that vasoocclusion is a complex process triggered by adhesive interactions among different cell groups in multiple stages. However, the quantification of the biophysical characteristics of sickle cell anemia remains an open issue. Based on dissipative particle dynamics, we develop a multiscale model for the sickle red blood cells (SS-RBCs), accounting for diversity in both shapes and cell rigidities, to investigate the precise mechanism of vasoocclusion. First, we investigate the adhesive dynamics of a single SS-RBC in shear flow and static conditions, and find that the different cell groups (SS2: young-deformable SS-RBCs, ISCs: rigid-irreversible SS-RBCs) exhibit heterogeneous adhesive behavior due to the diverse cell morphologies and membrane rigidities. We quantify the observed adhesion behavior (in static conditions) in terms of a balance of free energies due to cell adhesion and deformation, and propose a power law that relates the free-energy increase as a function of the contact area. We further simulate postcapillary flow of SS-RBC suspensions with different cell fractions. The more adhesive SS2 cells interact with the vascular endothelium and trap ISC cells, resulting in vasoocclusion in vessels less than 
depending on the hematocrit. Under inflammation, adherent leukocytes may also trap ISC cells, resulting in vasoocclusion in even larger vessels.
Sickle cell anemia (SCA) is a genetic disease originating from the abnormal sickle hemoglobin molecules (HbS). In hypoxic conditions, the intracellular HbS solution transitions into a polymerized state, resulting in a series of alterations in the cell membrane functions. According to the seminal study by Kaul et al. (1), suspensions of sickle red blood cells (SS-RBCs) contain heterogeneous cell density groups, which can be roughly divided into four fractions according to the intracellular mean corpuscular hemoglobin concentration (MCHC). Fractions I (SS1) and II (SS2) with moderate MCHC are mainly composed of reticulocytes and discocytes, respectively, with MCHC similar to healthy cells. On the other hand, fractions III (SS3) and IV (SS4) with high MCHC are mainly composed of rigid discocytes and irreversible sickle cells (ISC). Associated with the heterogeneous cell groups is the abnormal rheology of the SS-RBC suspensions. Compared with healthy blood, the SS-RBC suspensions exhibit higher flow resistance in shear flow (2, 3), as confirmed in microfluidic experiments (4) and in isolated vascular systems (5).
Early studies of SCA suggested that the major cause of vasoocclusion is the sickling process of the dense SS4 cells during their circulation in capillaries. However, subsequent investigations (5⇓⇓–8) have shown that vasoocclusion mainly occurs in postcapillaries. Moreover, it has been found that there is no direct correlation between the percentage of the dense SS-RBCs and the disease severity (9). Recent studies have revealed that vasoocclusion is a complex process triggered by the interactions between multiple density groups (10, 11), where each group contributes differently to the occlusion crisis. In addition, in an inflammation-activated state, leukocytes adherent to the vessel wall may further trap the SS-RBCs, causing occlusion in venular flow (12).
Despite these findings, understanding the precise mechanism of vasoocclusion and quantifying cell interactions, including interactions with the vascular endothelium, are still open research questions. Several numerical investigations (13, 14) have focused on the effect of membrane viscoelasticity and morphology of SS-RBCs on blood flow resistance. However, the heterogeneous SS-RBC adhesive properties, the physiological contribution to vasoocclusion for individual SS-RBC groups, and the hemodynamic effect of the leukocytes have not been fully understood yet. Numerical investigation of these issues may enhance our understanding of vasoocclusion and potentially reveal new paradigms for therapeutic treatments by targeting specific adhesive conditions. In this paper, we present a systematic simulation study of the biophysical characteristics of vasoocclusion in SCA.
The entire simulation framework is based on dissipative particle dynamics (DPD), including models for RBCs, leukocytes, plasma, and the wall boundaries (Figs. S1–S3, SI Text). Unlike the healthy SS-RBCs, the SS-RBC membrane expresses multiple types of adhesive receptors (15), which interact with ligands on the endothelial cells through different pathways. We use a stochastic coarse-grained model to represent the interaction between the SS-RBC and the endothelial cells, where the diverse types of adhesive proteins are modeled by effective coarse-grained particles expressed on the cell membrane and vascular wall (Fig. S4 and SI Text).
First, we examine the adhesive dynamics of single SS-RBCs of different density groups; our simulation results validate the aforementioned density-class dependence of the adhesive properties. Subsequently, we explore the adhesion mechanism by computing the free-energy increase as a response of adhesive interaction and quantify the magnitude of the force between different cell groups and the vascular wall. We then use the same sets of parameters to investigate the hemodynamics of SS-RBC suspensions mixed with different cell groups. By analyzing the distinct contribution of individual cell groups, we identify the specific physiological conditions that trigger vasoocclusion.
Experiments in ref. 16 have used an in vitro shear flow system to show that the light density cells exhibit larger adhesion, whereas the densest ISCs correspond to the least adhesive group. However, the following fundamental question remains unresolved: Are the different adhesive responses a consequence of the different adhesive proteins on the cell membranes, or rather of the biomechanical properties of the different cell groups? Kaul et al. (6) explored this issue by examining the alteration of SS-RBC adhesivity after a dehydration/rehydration treatment on individual cell groups and found that the cell adhesivity of the deformable SS2 and the dense SS4 cells can be reversed after controlled treatment. They hypothesized that both light- and high-density cells have similar “adhesion potential,” whereas the heterogeneous cell adhesivity is mainly due to the differences in cell deformability and shape peculiarities among the multiple cell groups. To investigate the validity of this hypothesis, we simulated the cell adhesive dynamics of SS-RBCs with different cell rigidity and morphologic characteristics.
We consider three distinct SS-RBCs under shear flow, as shown in Fig. 1. Following ref. 6, we assume that the three cells share similar “adhesive potential,” and hence we set identical adhesive parameters in all three simulations. Cell I represents a SS2 deformable discocyte of cell rigidity similar to the healthy RBC. We set the shear modulus 
and bending rigidity 
, i.e., similar to healthy cell rigidity (17⇓–19). Cell III represents an ISC of the SS4 type. Experimental measurements (20, 21) of cell rigidity report scattered values, which further depend on the deoxygenation–oxygenation cycle. For example, an ektacytometer study (21) showed that the cell rigidity of oxygenated ISC is ∼12 times larger than a healthy cell, whereas a micropipette study (20) showed that the cell rigidity of the SS4 group is about 2–3 times larger than the healthy cell; however, this value increases to 20 times the healthy value after one deoxygenation–reoxygenation cycle. In the present work, we set the shear modulus 
, because an ISC undergoes at least one deoxygenation–reoxygenation cycle. The modification to cell bending rigidity by the deoxygenation–reoxygenation process is unknown; here we set 
and also conduct sensitivity studies. Cell II represents a rigid discocyte cell (22) of SS3 type with medium MCHC value; we set 
and 
for the comparative study.
Sickle cells in shear flow: (A) Successive snapshots of SS-RBCs in shear flow, Movie S1. Labels I, II, and III correspond to a deformable SS2 cell, rigid SS3 cell, and ISC, respectively. Arrow indicates the flow direction. (B and C) Instantaneous contact area and velocity for SS-RBC in shear flow conditions.
Fig. 1A shows successive snapshots of the three cells. For the same shear rate 
and same membrane adhesive characteristics, the three cells exhibit substantially different adhesive dynamics (Movie S1). Cell I shows firm adhesion to the lower plate after the bond formation stage; the corresponding contact area (see SI Text for definition) is about 
, as shown in Fig. 1B. Cell II shows transient adhesion to the lower plate initially, but then undergoes a periodic flip movement along the flow direction, characterized by the peak values of the instantaneous cell velocity at 0.30, 0.53 and 
s. This indicates a weaker cell adhesivity than cell I. Accordingly, the contact area achieves minimum values at those times. The cell eventually detaches from the plate after two to three flips. Differently from cells I and II, cell III does not show firm/transient adhesion to the plate; instead, it directly detaches from the lower plate and moves freely without adhesive bonds established thereafter. Given the similar adhesive potential, the present results validate the hypothesis of heterogeneous cell adhesive dynamics under shear flow conditions. To further quantify the effect of the cell rigidity and morphology on adhesive properties, next we investigate cell adhesion in static conditions by using free-energy analysis.
Fig. 2 shows the instantaneous contact area between the cell and the plate for the three SS-RBCs. As in the shear flow conditions, the contact area shows an inverse relationship with the cell rigidity. Whereas the contact area for all of the three cells increases to 
within the initial stage, the contact areas between the plate and cells I and II further increase to 
and 
, respectively, at the final stage. These results can be understood by analyzing the change of free energy 
due to the adhesive interaction, i.e.,
where 
represents the free-energy change due to the cell deformation, and 
represents the free-energy change due to the adhesive bond formation between the cell and adhesive ligands. In the present system 
can be approximated by 
, where 
represents the change of contact area and β is the energy coefficient determined by the adhesive affinity. Fig. 2B shows 
for different adhesive affinities (therefore yielding different contact areas). The numerical results are fitted by 
, where α depends on the cell rigidity of individual cells. The equilibrium state is determined by the balance between the two energy terms, e.g., 
, leading to different contact areas. For SS2 cells (corresponding to small α), the adhesive interaction plays a dominant role in the initial stage, resulting in cell deformation and further increase of the contact area. In contrast, a rigid SS4 cell exhibits “solid”-like properties with a larger energy barrier for cell deformation. After the initial stage, the adhesive interaction driven by the bond formation cannot overcome the free-energy increase induced by the cell deformation and therefore prohibits further increase of contact area.
Adhesion in static conditions: (A) Instantaneous contact area between the SS-RBC and the plate coated with adhesive ligands. Black curve represents the contact area of a discocyte with cell rigidity similar to the ISC (cell III). (B) Increase of cell free energy as a function of the contact area. (C) Adhesive force between the cell and the plate as a function of the membrane rigidities for cell morphologies. Error bar represents the adhesive force computed from four independent simulations. (Inset) Sketch of the simulation setup, where a uniform lift force is applied on the upper part of an ISC.
Compared with the shear flow system, we note that the equilibrium contact area obtained in the static condition is smaller. This discrepancy is mainly due to the greater cell deformation under the shear flow condition, i.e., part of the deformation free energy 
is balanced by the hydrodynamic force exerted on the cell membrane under the shear flow condition. The extended cell membrane facilitates the bond formation and results in a larger contact area. Moreover, compared with the ISC, a discocyte with similar cell rigidity yields a larger contact area in equilibrium state, as shown in Fig. 2A. This result demonstrates that the peculiar elongated and curved characteristics of ISC may further prevent the membrane receptors from interacting with the ligands.
To quantify the effects of shape and rigidity on cell adhesion, we compute the adhesive force between the cell and plate for different cases. Starting from the steady state obtained from the static condition, we apply uniformly a lift force to the upper part of the cell membrane (Fig. 2C, Inset). The adhesive force is determined as the lift force that drives the cell detaching from the plate in a quasistatic way. The deformable SS2 cells exhibit the largest adhesive force of 
, whereas this value decreases to 34 pN in the case of rigid cells with 
and 
. Compared with the discocyte, the ISC with similar cell rigidity yields an adhesive force of 17 pN. For smaller cell rigidity (corresponding to the rehydration treatment of ISC in ref. 6) with 
and 
, the adhesive force is approximately 24 pN, compared with 43 pN for a discocyte with similar cell rigidity. In a recent experimental study (23), force spectroscopy was used to quantify the cytoadherence of RBC invaded by a Plasmodium falciparum parasite. It was reported that the rupture force between the invaded RBC and the CHO cells is within the range 
. This result is similar to the adhesive force magnitude for the RBC–endothelium interaction reported in the present work.
Single ISCs may occasionally result in occlusion at precapillary junctions (1, 5, 24); however, most of the occlusion sites are in postcapillaries and exhibit a specific cell pattern, where the adherent deformable cells form a sieve-like configuration and selectively trap the ISCs. We model blood circulation by SS-RBC suspensions first in a tube of diameter 
, as shown in Fig. 3; small green particles represent the adhesive ligands coated on the tube wall. To quantify the distinct role of different cell groups, we infuse different cell groups into the tube and apply a pressure gradient 
.
Vasoocclusion in postcapillaries: instantaneous mean velocity of blood flow in a cylindrical tube of 
infused with different SS-RBC suspensions. Red curve represents the resultant velocity infused with SS2 and ISC cell groups. (Inset) Instantaneous snapshots where SS2 cells adhere to vessel wall, consequently trapping the ISCs and resulting in cell blockage. Green curve represents the blood velocity infused with SS2 and ISC cell groups, where adhesive interaction is only applied to the ISC group. The time axis is scaled by 0.5 for better visualization. The Inset represents a snapshot where transient adhesion is established between ISC and the tube wall. Steady flow is recovered as the cell detaches from the tube wall. Blue curve represents the instantaneous velocity of blood flow infused with SS2 and healthy cell groups. Blood flow exhibits a slowdown but not a full occlusion; Movies S2, S3, and S4.
First, we consider suspensions composed of SS2 cells (labeled by blue) and ISC cells (labeled by red) with hematocrit 
, similar to experiments (5, 6). Steady flow is achieved by turning off the adhesive interaction during the initial stage, yielding a mean flow velocity of 
. Starting from the steady state 
, we turn on the adhesive interaction between the SS-RBCs and the ligand particles and compute the instantaneous flow velocity across the tube, as represented by the red curve in Fig. 3. Steady flow is maintained until one of the SS2 cells adheres to the vessel wall, triggering a sharp decrease of blood flow around 
. As a positive feedback, the decreased blood flow induces more SS2 cells adherent to the vessel wall, leading to a further decrease of flow rate at 
and 
. Moreover, these adherent cells decrease the effective vessel lumen near the adherent region, resulting in a secondary trapping of the ISC groups. The final occlusion state is achieved around 
with cell patterns similar to the experimental observations (5, 6). We have also performed another three simulations starting from different initial conditions, and in all cases we obtained the same final pattern of full occlusion.
To explore the unique contribution of the deformable cell group, we performed a similar simulation of blood flow where the SS2 cell–ligand interaction is turned off. Starting from the steady flow at 
, blood flow is simulated for 6 s, as represented by the green line in Fig. 3. Due to the ISC–ligand interaction, transient adhesive bonds can be formed, resulting in the decrease of blood flow at 
. However, firm adhesion cannot be established: blood flow can recover the initial flow rate when the adherent ISCs detach from the tube wall. Blood occlusion is not observed within the simulation time, in agreement with the experimental observation (6) that the ISCs do not result in consistent microvascular blockage when infused alone. This finding reveals the distinct role of the SS2 cell in the vasoocclusion crisis and is consistent with the positive correlation between the SS-RBC deformability index and disease severity reported in clinical observations (9). The ISC cell group, in contrast, contributes differently to the occlusion crisis. We quantify its unique contribution by simulating the blood suspension mixed with the deformable SS2 and healthy cells. Starting from the steady flow at 
, we compute the mean velocity with the adhesive interaction between the SS2 and ligands, as represented by the blue line in Fig. 3. The flow velocity decreases due to the cell adhesion at 
and 
. Although blood flow exhibits a substantial slowdown, full occlusion is not realized due to the high deformability of the healthy RBCs, which enables the cells to squeeze through the sieve-like region formed by adherent SS2 cells. The present result, in turn, reveals the distinct role of the ISC group in the vasoocclusion process. Although the least adhesive, the ISC group, due to its high cell rigidity and peculiar shape, serves as the particular cell group trapped by adherent cells in the postcapillaries. This finding explains the experimental observation that a large number of dense cells accumulate in the occlusion region but are not observed in the peripheral blood circulation (5).
To quantify the effect of the vascular size on occlusion, we have performed similar simulations but in tubes of different diameters. For 
, occlusion does not occur with tube diameter larger than 
(Fig. S5); the sieve-like pattern formed by SS2 cells cannot fully trap the ISCs due to the larger spatial accommodation near the adhesion region. However, for 
, occlusion is observed with 
but occlusion is prevented with tube diameter larger than 
(see Fig. S6 and SI Text for more details).
Recent studies (25) have shown that SCA is often accompanied by an inflammatory endothelial phenotype, resulting in elevated leukocyte recruitment in blood circulation. Moreover, studies by Turhan et al. (12) in transgenic-knockout mice have shown that the inflammation-stimulated (by cytokine TNF-α) adherent leukocytes, upon interactions with SS-RBCs, lead to occlusion in venular flow. We use our computational framework to investigate the vasoocclusion induced by inflammation-stimulated leukocytes (Methods). We first consider blood flow of SS-RBC suspension 
with a single leukocyte in a tube of diameter 
, as shown in Fig. 4. Due to cell migration, the leukocyte touches the wall at 
and firm adhesion is established thereafter. The mean flow velocity drops from 
to 
during this stage. During the inflammation-stimulated stage, we turn on the cell aggregation interaction between the SS-RBCs and leukocytes (see SI Text for details). Multiple SS-RBCs get trapped by the leukocyte, resulting in full occlusion at 
. For venular flow with larger diameter, multiple leukocytes may accumulate in the inflammation region. Similar to the experiment (12), we simulate a SS-RBC suspension 
with three leukocytes in a tube of diameter 
and compute the mean flow velocity during different time intervals, as shown in Fig. 4. Stage I represents the blood flow during the initial stage of the inflammatory response. Starting from the steady flow at 
, the mean flow velocity drops to 
due to the attachment of leukoctyes on the vessel wall at 
s. Stage II represents the blood flow with moderate RBC–leukocyte interaction; the blood flow exhibits a slowdown due to the adherent leukocytes. Stage III represents the late stage of the inflammatory response, where the RBC–leukoctye interaction is further intensified. Multiple SS-RBCs are trapped on the adherent leukocytes, leading to full occlusion at 
(Movies S5, S6, S7, S8, S9, and S10).
Effect of leukocytes: instantaneous mean velocity of the blood flow in a tube of 
, one leukocyte (A) and 
, 
, three leukocytes (B). (Inset) Snapshots represent blood cells in free motion, leukocyte adhesion, and blood occlusion states. For 
, the Inset plot represents the blood flow velocity of the present study (blue) and the experimental results (red) (12), where measurements are taken on 
venules with average diameter 
and 
before and after inflammation–stimulation; Movies S5, S6, S7, S8, S9, and S10.
Our simulations of individual SS-RBCs in shear flow validated the hypothesis of the importance of cell rigidity and morphology (5, 6, 16), and elucidated the distinct behavior of each individual cell group under adhesive conditions. This behavior was further quantified by computing the adhesive force and free -response for all cell groups. In simulations of postcapillary SS-RBC flow, we identified which combination of cell groups and at what adhesion strengths lead to full or transient occlusion states. Under typical physiological conditions, our sensitivity studies with respect to the tube size and hematocrit show that blood occlusion mainly occurs in postcapillaries with diameter smaller than 
, a finding consistent with the experimental results reported in ref. 5. Further free energy analysis on cell adhesion (Fig. S7) and discussion on the important dimensionless groups is included in SI Text. Moreover, under inflammation-stimulated conditions, the present work quantified the influence of the adherent leukocytes, which may further trap the SS-RBC and lead to vasoocclusion in venular flow, especially in larger vessels. In particular, our simulations revealed and quantified the multiinteractional and multistage nature of vasoocclusion, which may be useful in developing effective biomarkers (26) and therapeutic treatments for SCA. In general, our computational framework provides a capability, which, in conjunction with microfluidic experiments, can be used in the preliminary screening of proposed drugs and the interpretation of clinical outcomes.
The DPD method is a Lagrangian-based coarse-grained particle method (27, 28) used to simulate complex fluids (29⇓–31). Each DPD particle can be viewed as a coarse-grained virtual cluster of multiple atomistic particles (32), where the particle motion is governed by pairwise force imposed between the DPD particles (Table S1 and SI Text).
A healthy RBC (33, 34) is modeled by a triangulated surface (e.g., 500 vertices) of the biconcave RBC shape (35) in equilibrium state. Each vertex on the membrane surface is represented by a DPD particle. Membrane elasticity is imposed by viscoelastic bonds between the DPD particles. Bending rigidity is modeled by the bending resistance between the adjacent triangulated surface. In addition, area and volume constraints are imposed to represent the incompressibility of the lipid bilayer and RBC cytosol. The model has been validated in a number of studies, e.g., for the mechanical properties (33, 34), as well as rheological and hemodynamic properties of blood flow in healthy (36), malaria-infected (37), and SCA conditions (38). The ISC model is constructed by applying surface tension on a healthy RBC membrane (38), representing the distortion of the cell membrane due to the growth of the intracellular sickle hemoglobin polymers. The distorted shape is redefined as the equilibrium state of the ISC with minimum free energy. The leukocyte is modeled by a triangulated surface (e.g., 1,000 vertices) of a sphere of radius 
. Similar to the RBC model, proper cell rigidity is imposed by the viscoelastic bonds and bending resistance; total area and volume constraints are imposed for individual cells. More details on the cell models are presented in SI Text.
The adhesive interaction between the SS-RBC/leukocyte and the adhesive ligands on vascular endothelium is modeled by the stochastic bond formation/dissociation according to refs. 39, 40. Cell vertices can interact with the endothelial ligands within interaction distance 
. During each time step 
, transient bonds can be formed between cell vertices and endothelial ligands with probability 
, whereas existing bonds can be ruptured with probability 
, where 
are, respectively, bond formation and dissociation rate depending on the adhesive affinity between the blood cell and vascular endothelium (Table S2). The interaction between the leukocyte and SS-RBC is modeled by the Morse potential, 
, where r is the distance between cell-membrane vertices of adjacent cells, β determines the range of interaction, and 
and 
define the equilibrium distance and attractive strength, respectively. See SI Text for more details.
H.L. acknowledges Dr. Xin Yi (Brown University), Dr. Petia Vlahovska (Brown University), and Dr. Dmitry A. Fedosov (Forschungszentrum Jülich) for useful discussions. This work was sponsored by National Institutes of Health Grant R01HL094270 and the new Collaboratory on Mathematics for Mesocopic Modeling of Materials (CM4) supported by the Department of Energy (DOE). Computations were possible due to a DOE/Innovative and Novel Computational Impact on Theory and Experiment Award.
Author contributions: H.L. and G.E.K. designed research, performed research, analyzed data, and wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1221297110/-/DCSupplemental.
