Patient-specific modeling of individual sickle cell behavior under transient hypoxia

Sickle cell disease (SCD) is a highly complex genetic blood disorder in which red blood cells (RBC) exhibit heterogeneous morphology changes and decreased deformability. We employ a kinetic model for cell morphological sickling that invokes parameters derived from patient-specific data. This model is used to investigate the dynamics of individual sickle cells in a capillary-like microenvironment in order to address various mechanisms associated with SCD. We show that all RBCs, both hypoxia-unaffected and hypoxia-affected ones, regularly pass through microgates under oxygenated state. However, the hypoxia-affected cells undergo sickling which significantly alters cell dynamics. In particular, the dense and rigid sickle RBCs are obstructed thereby clogging blood flow while the less dense and deformable ones are capable of circumnavigating dead (trapped) cells ahead of them by choosing a serpentine path. Informed by recent experiments involving microfluidics that provide in vitro quantitative information on cell dynamics under transient hypoxia conditions, we have performed detailed computational simulations of alterations to cell behavior in response to morphological changes and membrane stiffening. Our model reveals that SCD exhibits substantial heterogeneity even within a particular density-fractionated subpopulation. These findings provide unique insights into how individual sickle cells move through capillaries under transient hypoxic conditions, and offer novel possibilities for designing effective therapeutic interventions for SCD.

Sickle cell disease (SCD) is a highly complex genetic blood disorder in which red blood cells (RBC) exhibit heterogeneous morphology changes and decreased deformability. We employ a kinetic model for cell morphological sickling that invokes parameters derived from patient-specific data. This model is used to investigate the dynamics of individual sickle cells in a capillary-like microenvironment in order to address various mechanisms associated with SCD. We show that all RBCs, both hypoxia-unaffected and hypoxia-affected ones, regularly pass through microgates under oxygenated state. However, the hypoxia-affected cells undergo sickling which significantly alters cell dynamics. In particular, the dense and rigid sickle RBCs are obstructed thereby clogging blood flow while the less dense and deformable ones are capable of circumnavigating dead (trapped) cells ahead of them by choosing a serpentine path. Informed by recent experiments involving microfluidics that provide in vitro quantitative information on cell dynamics under transient hypoxia conditions, we have performed detailed computational simulations of alterations to cell behavior in response to morphological changes and membrane stiffening. Our model reveals that SCD exhibits substantial heterogeneity even within a particular density-fractionated subpopulation. These findings provide unique insights into how individual sickle cells move through capillaries under transient hypoxic conditions, and offer novel possibilities for designing effective therapeutic interventions for SCD.

Author summary
Sickle cell disease is a genetic blood disease that causes vaso-occlusive pain crises. Here, we investigate the individual sickle cell behavior under controlled hypoxic conditions through patient-specific predictive computational simulations that are informed by companion microfluidic experiments. We identify the different dynamic behavior between individual sickle RBCs and normal ones in microfluidic flow, and analyze the hypoxia-

Introduction
In research investigations of hematological disorders, most experiments are performed on groups of cells with the underlying assumption that all of the cells in a particular type are identical. However, recent evidence reveals that individual cells within the same population may differ drastically in size, shape, mechanical properties and protein levels, and these variations can have important consequences for the health and biological function of the entire cell population [1]. A representative case is sickle cell disease (SCD), one of the most common inherited genetic blood disorders affecting more than 270,000 new patients each year [2,3]. SCD has been characterized as the first molecular disease [4], being linked to the mutation of a single nucleotide in the hemoglobin molecule. The primary pathophysiological event in SCD is the polymerization of sickle hemoglobin (HbS) into long fibers upon deoxygenation (DeOxy) [5,6]. The fibers distort RBCs into irregular and heterogeneous shapes-e.g. granular, elongated, oval, and crescent (classic sickle) shapes [7,8]. The hypoxia-affected RBCs are also heterogeneous in their cell density in a range of less than 30 g/dL to more than 46 g/dL [9], which are usually fractioned into four arbitrary cell density subpopulations (fractions I-IV) in vitro analysis [7]. Heterogeneous cell fractions engender heterogeneity in cell rigidity [10][11][12][13]. These hypoxia-affected RBCs are more sticky and stiff, causing frequent painful episodes of vaso-occlusion and depriving oxygen from tissues and organs [10,14]. The decrease in RBC deformability contributes to impaired blood flow and other pathophysiological origins of the disease. However, the clinical expression of SCD is heterogeneous, as the hypoxia-affected RBCs do not all behave in the same way all the time, and the variance is considerable even within a same density-fraction, making it hard to predict the risk of a vaso-occlusive crisis [15][16][17]. This poses a serious challenge for disease management.
Precision medicine [18], which accounts for individual variability, is an emerging approach for treatment and prevention of disease [19,20]. Developing such an approach, however, inevitably requires resolving various heterogeneity-related issues, both at whole cell population and single-cell levels [21]. Although developments in quantitative, in vitro microfluidic assays provide greater understanding of cell dynamics under hypoxic conditions in patients' blood samples [16], there is a critical need to develop single-cell level assays to assess sickle RBC behavior under transient hypoxic conditions for better therapeutic interventions [22]. These considerations lead to the motivation for the present work whose aim is to address the following question: To what extent does morphological sickling affect cell dynamics in microcapillaries under transient hypoxia, with consequences for hemodynamics and vaso-occlusion?
In order to gain a better understanding of vaso-occlusive crisis, it is necessary to obtain direct and real-time observations of the traversal individual sickle RBCs through microcapillaries and of alterations in cell biodynamics and biorheology in response to controlled changes in oxygen (O 2 ) concentration. However, owing to the phenotypic heterogeneity of SCD, existing experimental capabilities do not readily provide this information. As a result, there is a critical need to develop predictive, patient-specific cell models to quantify alterations in cell biodynamics under transient hypoxic conditions. Such models could be validated by recourse to a variety of well-controlled and independent in vitro experiments.
In this article, we investigate the dynamic behavior of individual sickle cells in real time, in a capillary-like microenvironment, and perform quantitative analysis of hypoxia-induced alteration in cell behavior and response to obstruction of capillaries by combining predictive simulations with controlled and quantitative information obtained from microfluidic experiments.

Sickle RBC samples
De-identified SCD blood samples from two patients with HU therapy (on-HU) and two patients without HU therapy (off-HU) were selected for this study, following institutional review board (IRB) approvals from the National Institutes of Health (NIH) and Massachusetts Institute of Technology (MIT). All samples were collected into EDTA anticoagulant and stored at 4˚C during shipping and storage. Table 1 shows selected hematologic and hemorheologic parameters in these four blood samples. The experiments were performed exactly in the same way as those in Ref. [23], but with a specific focus on (1) correlations between cell shape, mean corpuscular hemoglobin concentration (MCHC), transit velocity and trajectories through the microgates, and (2) movements of individual cells with and without disturbances from adjacent blockages.
After washing twice with PBS (Thermo Scientific) at 821 × g for 5 min at 21˚C, the RBCs in the blood sample were separated into four density fractions, using a stepwise gradient medium prepared with OptiPrep medium (Sigma Aldrich, Saint Louis, MO, USA) and Dulbecco's PBS (HyClone Laboratories, Inc., South Logan, UT, USA). The estimated MCHC values were 27.3, 30.9, 34.9, and > 45.0 g/dL for fraction I-IV, respectively (Table 2). Fraction I (SS1, reticulocyte rich) and fraction II (SS2, discocyte rich) have moderate MCHC values, which are similar to those of healthy RBCs. Fractions III (SS3) and IV (SS4) mainly comprised rigid discocytes and irreversible sickle cells (ISCs), respectively, with their MCHC values considerably higher than The symbols S-P-I and S-P-II represent two blood samples from SCD patients not treated with HU, whereas S-P-III and S-P-IV represent the other two blood samples from SCD patients treated with HU. Hct, hematocrit; MCV, mean corpuscular volume; MCHC, mean corpuscular hemoglobin concentration; HbA, hemoglobin A (α 2 β 2 ), also known as adult hemoglobin; HbA 2 , a normal variant of hemoglobin A (α 2 δ 2 ); HbF, fetal hemoglobin (α 2 γ 2 ); HbS, sickle hemoglobin (α 2 β S 2 ). those of healthy RBCs. The fractionated RBCs were washed with PBS and re-suspended in RPMI-1640 medium with 1% (wt/vol) Bovine Serum Albumin (Sigma-Aldrich, Saint Louis, MO, USA) for cell sickling measurement in a polydimethylsiloxane (PDMS)-based microfluidic hypoxia assay [23]. This platform provided measurements of cell sickling under controlled O 2 concentrations at 37˚C, including a fully Oxy state (20 vol% O 2 ), and a hypoxic condition, in which the O 2 concentration decreased from 20 vol% to below 5 vol% within 15 s and maintained at 2 vol% for the rest periods. Patient-specific data relevant to the present work along with cell morphological data are provided in Table A in S1 Text for predictive simulations.

Computational model and method
In order to investigate the behavior of sickle RBCs in each density-fractionated subpopulation under transient hypoxic conditions, we have developed a unified modeling framework based on dissipative particle dynamics (DPD). Our computational framework is predicated on recent developments [23] in microfluidics that quantify in vitro single RBC dynamic response and its heterogeneity under controlled oxygen partial pressures in blood samples from SCD patients. Specifically, the geometry of microfluidic channel (Fig 1) and profile of transient hypoxia ( Fig  1 and Figure A in S1 Text) are the same as in the experimental setup [23]. The solid walls of  the microfluidic channel in DPD simulations are modeled by layers of frozen DPD particles, and a force-adaptive is employed for fluid particles to control their density fluctuations [24]. An external body force is exerted on each fluid particle to generate a flow in the microfluidic channel. In this study, the externally applied body force works in the direction of flow (x-direction) on the fluid particles. We studied the dynamics of sickle RBCs under transient hypoxia using a multiscale RBC model [25]. For completeness, the model is briefly summarized below, whereas details of this MS-RBC model are available in Refs. [25,26]. Here, the membrane of MS-RBC is modeled by a two-dimensional triangulated network with N v vertices, which are connected by N s viscoelastic bonds to impose proper membrane mechanics. In this study, the elastic energy of the bond is taken as where l j and l m are the equilibrium (initial edge) length and maximum extension of spring j, p is the persistence length, k p is a constant factor of the spring, and n is an exponent. The bending energy of the RBC membrane is written as where k b is the bending constant, θ j and θ 0 are the instantaneous and spontaneous angles between two adjacent triangles with a shared common edge j. In addition, the volume and surface area of RBC are controlled to mimic the incompressible interior fluid and the area-preserving cell membrane. The corresponding energy is taken as where k a and k d are the global and local area constraint coefficients, k v is the volume constraint coefficients, A 0 is the triangle area, and the terms A 0 tot and V 0 tot are the total area and volume in equilibrium, respectively. The RBC membrane interacts with the fluid and wall particles through DPD forces, and the system temperature is maintained by the DPD thermostat. The surrounding external fluids and cytosol are modeled by collections of free coarse-grained particles and their separation is enforced through bounce-back reflections at the moving surface of RBC membrane. In addition, we consider a positive correlation between the number of intracellular coarse-grained particles and the MCHC value, i.e., for a denser RBC with a higher MCHC value, we include more coarse-grained particles inside the cell. The MS-RBC model has been validated by comparisons with a number of available experiments that examine the mechanics and biorheology of healthy and diseased RBCs [25][26][27][28][29][30][31]. Such connections between simulations and experiments have also indicated that the results developed here are not dependent on the details of refinement of the simulation in that they are independent of the level of coarse-graining for N v ! 500 [25,26]. In this study, we employ a MS-RBC model with a level of coarse-graining (N v = 500) that facilitates computationally efficient simulations of RBCs in microfluidic channel.
Sickle RBCs undergo various morphological transitions, due to the polymerization of intracellular HbS molecules, from the normal biconcave shape to an irregular sickled shape, in the form of granular, elongated and crescent shapes. Instead of directly modeling the formation of HbS polymer fibers and the resulting cell morphology, we consider an "effective surface tension" (stretching force) applied on the cell membrane to mimic the cell distortion exerted by the growing HbS polymer domain [31,32]. We develop a two-step RBC morphological sickling model, as described below.
1) We first apply surface tension on the cell membrane to transform the discocyte RBC into a granular, elongated or crescent shape. Here, we define the directions along the two long axes of RBC as x-and y-directions, and the third direction, perpendicular to the long axes, as the zdirection. We choose four anchor points located at the maximum/minimum values in the x/y directions (Fig 1), where the intracellular polymer fibers can potentially interact with the red cell membrane. Different combinations of the stretching force applied on the anchor points cause different sickle cell shapes: when the stretching force is exerted on two diametrically opposite anchor points (for example, points A and C), an elongated or a crescent cell shape is obtained. Similarly, if the stretching force is exerted on all the four anchor points, a granular shape results. The values of stretching forces for obtaining different shapes of sickle RBCs are shown in Table B in S1 Text. We define the distorted shape as the equilibrium state of the sickle RBC with minimum free energy. Specifically, we adjust the equilibrium length, l j , of each spring j to the edge length, l j,REF , of the final distorted state to eliminate local stress on the cell membrane.
2) Subsequently, we employ the kinetic description for sickling RBCs [33] to model the cell morphological sickling process under transient hypoxia. Here, we define the two parameters, l j and l j,REF , as the reference length of the spring j, respectively, in the RBC under normal Oxy state (~20 vol% O 2 ) and fully DeOxy state (~2 vol% O 2 ). Once the oxygen partial pressure, p O2 , is lower than a critical value p O2,c for cell sickling, we assume that a sufficiently large stretching force, f j (t), is applied on each spring so that its instantaneous edge length, l j (t), is essentially equal to the expected edge length, where t D is delay time for cell sickling. In the present study, the parameters of p O2,c , p O2 (t), and t D are directly obtained or calculated from experiments [23]. Through this approach, we are able to dynamically adjust the bond lengths of the elastic springs to affect the sickling and unsickling processes of RBCs while maintaining the cell membrane free of any total external force. Membrane stiffening of sickle RBCs is closely related to their morphologic change [11]. It seems likely that a step increase in rigidity occurs at the same time as morphological sickling. Similarly, the shear modulus, μ ref (t), of the hypoxia-affected cell membrane is taken as where μ REF is the shear modulus of sickle RBC, described below. The cell unsickling process is analogous to the cell sickling process, by considering the delay time (t D,R ) of cell unsickling. This kinetic cell sickling model, which accounts for (1) cell morphological change, (2) cell sickling and unsickling delay times, and (3) cell membrane stiffening under transient hypoxia, offers an effective method of investigating the sickling and unsickling processes of RBCs in response to changes in O 2 concentration under both static and flow conditions (S1 Video).
In this study, we modeled a healthy RBC using the multiscale RBC model with the following parameters: number of RBC vertices N v = 500; RBC area A 0 = 135.2 μm 2 and volume V 0 = 92.4 μm 3 ; RBC membrane bending modulus k c,0 = 2.4 × 10 −19 J; shear modulus μ 0 = 4.7 pNÁμm -1 ; cytosol viscosity η 0 = 1.2 cp. In contrast to a healthy RBC, the SCD RBC is characterized by extensive impairment in deformability, the extent of which is dependent on the cell density and oxygen partial pressure. Sickle RBCs in different cell density subpopulations exhibit different membrane elasticity in the Oxy state (O 2 concentration~20 vol%) and the DeOxy state (O 2 concentration < 5 vol%). Previous studies have shown that the intracellular polymerization of HbS upon DeOxy leads to significant increases in cytosol viscosity and membrane elasticity [34][35][36][37][38][39]. The cytosol viscosity, η cytosol , under DeOxy state could be two orders of magnitude greater than η 0 [34,39]. In our previous sensitivity study, we have shown that the blood dynamics is nearly independent of cytosol viscosity if η cytosol > 50η 0 [32]. We would expect that the cytosol viscosity also plays a major role in determining sickle cell behavior in a capillary-like microenvironment under transient hypoxia. However, we find the cytosol viscosity in the densest cell fractions IV is only around 4.5 η 0 , because our model does not explicitly include the intracellular HbS polymer fibers. Here, we used an effective shear modulus of RBCs, which accounts for the effects of (1) membrane stiffening of sickle RBCs, (2) high membrane tension induced by the growth of intracellular HbS polymers, and (3) elevated cytosol viscosity. In the Oxy state, we select experimentally determined shear modulus data [35][36][37], and set the effective shear modulus to μ REF = 1.0μ 0 , 1.2μ 0 , 1.5μ 0 and 3.0μ 0 for RBCs in fractions I, II, III and IV, respectively. In the DeOxy state, we consider four distinct types of sickle RBCs with different cell membrane mechanical properties based on experimental observations [11,40]. For an SS1 deformable cell with a low MCHC value, the measured shear modulus increased by up to 10 times after the occurrence of sickling [11], and hence we set the effective shear modulus μ REF = (5-10) μ 0 . For an SS2 cell with a mild MCHC value, in which the effective shear modulus is increased by one or two orders of magnitude compared to that of healthy RBCs [11,40], we set μ REF = (50-100) μ 0 . For an ISC of the SS4 type, its effective shear modulus is at least two or three orders of magnitude greater than the value of healthy RBCs [11], so we set μ REF = (1000-2000) μ 0 . For a SS3 rigid cell with a higher MCHC value, we set μ REF = (250-500) μ 0 . In addition, for the hypoxia-unaffected RBCs in each cell density fraction, we assume comparable deformability with the population in the same cell density fraction under the Oxy state. Similar parameters have been used as inputs for DPD simulations to quantify the shearindependent rheological behavior of blood flow in SCD [32]. We also investigated the functional dependence of shear viscosity of sickle RBC suspensions on the effective shear modulus, and demonstrated that the sickle RBC suspension behaves as a Newtonian fluid only if the effective shear modulus of sickle RBCs increases by two or three orders of magnitude [31]. The bending rigidity k c , following prior work [31,32], is kept constant in the model as the hemoglobin concentration increased.
The cell sickling delay time due to intracellular HbS polymerization is demonstrated to be the primary determinant of clinical severity in SCD [41]. The delay time of cell sickling is extraordinarily sensitive to solution conditions, particularly to HbS concentration. From prior results obtained from in vitro microfluidics experiment on sickle RBCs [23] and in vivo studies [41], we set the mean t D % 8.7 s, 19.8 s and 23.6 s for granular, elongated and classic sickle shaped RBCs, respectively, for the off-HU group blood samples. In the microfluidic experiments, the delay time for sickling for the on-HU group blood samples varied from 28 to 100 s [23], which is much longer than that for the off-HU group. Based on the experimental results, here we set the mean t D to be in the same range for the on-HU group. The cell unsickling delay time, t D,R , is also directly obtained from the microfluidic experiments. According to Du et al. [23], the cell unsickling process after reoxygenation (ReOxy) was much faster (< 20 s) than the cell sickling process, and the delay time distribution of cell unsickling was not significantly different between the two different groups. Based on the experimental results, here we set the mean t D,R to be in the range of 10-15 s.
All simulations were carried out with a representative volume of 180 μm × 60 μm × 5 μm that comprised a total of 30 periodic obstacles with a fluid particle number density of 3 obtained by dividing the fluid particle number of the model system by its volume. Capillary obstruction statistics were collected by running 200 independent simulations and the collective behavior is found to be independent of the initial position and initial orientation of RBCs. For the case of flow under Oxy state (without obstruction), steady state is reached with an average velocity of~120 μm/s. The simulations are performed using the USER MESO functional package that was written based on LAMMPS [42]. The time integration of the motion equations is computed through a modified velocity-Verlet algorithm with λ = 0.5 and time step Δt = 0.001 τ where τ is a characteristic time in DPD units. A typical simulation performed in the current study involves one million time steps and a computing time, on average, of about 20,000 CPU core hours on the Blue Gene/P system at the Argonne Leadership Computing Facility (ALCF).

Behavior of individual sickle cells under transient hypoxia
Individual sickle cells show marked heterogeneity, and the variance is considerable even within the same density-fraction. This leads to the question: how do individual sickle RBCs from different density-fraction behave differently from healthy ones when they travel though microcapillaries under transient hypoxic conditions? Here, we investigate the motion of individual RBCs flowing in a capillary-inspired microchannel that consists of parallel periodic obstacles forming 15-μm-long, 4-μm-wide and 5-μm-high microgates. Fig 2 and S2 Video in supplementary material, show typical dynamic motion of RBCs travelling through the microfluidic channel under cyclic hypoxia. This figure and the accompanying video reveal that all RBCs are easily deformed and pass readily through the microgates. Following a decrease in O 2 concentration, some RBCs become sickled and get trapped at the microgates; other RBCs can still deform and squeeze through the microgates (Fig 2). When hypoxia continues for a prolonged period, more RBCs become sickled and are obstructed at the microgates. After the O 2 levels are restored, these rigid RBCs remain trapped for a few seconds due to a delay in the cell unsickling process before regaining their deformability; after this delay period, they once again easily traverse the microgates (Fig 2). The simulated individual cell behavior under controlled hypoxia conditions is in qualitative agreement with experimental observations [23]. In addition, we calculated the transit velocity of individual sickle RBCs under transient hypoxia and compared them to experimentally measured data, see Fig 2. We found that the transit velocities obtained from experiments and simulations are mutually consistent, under both Oxy and DeOxy states.
A direct means of triggering the early stages of a vaso-occlusive event is possible by following trajectories of sickle RBCs under controlled transient hypoxia at the single-cell level. For the purpose of illustration and to help with the analysis, we track the trajectories of individual sickle RBCs flowing in the capillary-inspired microchannel, both from patient-specific predictive simulations and companion microfluidic experiments, see Fig 3. These figures show that the sickle RBCs have two different types of motion, i.e., translational and flipping motion, before they arrest at the microgates. These different types of motion may lead to transient or even permanent occlusion of the microcapillary. We observe that a sickle RBC favors moving along the flow direction (translational motion) when it undergoes cell sickling (Fig 3 and corresponding S3-S5 Videos). Specifically, for elongated and crescent-shaped (classic sickle) cells, more than three quarters of the cells follow the translational motion. The flow of an RBC approaching the entrance of a microgate is disturbed because of the size mismatch between the cell and microgate. While a deformable cell is able to squeeze through the microgate, a stiffened (sickled) cell loses its ability to deform dynamically. Instead, it undergoes a continuous rotation until its long axis is nearly parallel to the flow, allowing a part of the cell to enter into the microgate. These cells are usually arrested at the microgates in a parallel manner, i.e., the sickle RBC tends to align with flow streamlines, as shown in Fig 3. Our observations show that such a blockage is initially transient (Fig 3; S3-S5 Videos), due to a relatively smaller contact area of the trapped sickle RBC than the cross-sectional area of the microgates. The trapped sickle RBC may move slowly through the microgates (S3 Video), or eventually stop at the microgates, causing persistent obstruction to RBC flow (S4 and S5 Videos), if the sickle RBCs become stiff, e.g., an ISC in fraction IV. Interestingly, if a sickle RBC moves alongside with other RBCs, while travelling or being trapped, it may also experience considerable rotation because it is subject to a velocity gradient even when the sickle RBC aligns with flow streamlines. Such a cell is more likely to flip, leading to occasionally longitudinal or vertical blockage (sickled RBC aligns perpendicular to flow direction) (Fig 3 and corresponding S3 Video and S6 Video). It should be noted that the vertical blockage is much more likely to become persistent blockage, with more serious implications for vaso-occlusion in SCD patients.
We also track the dynamic behavior of individual sickle RBCs at different density-fractionated subpopulations under controlled transient hypoxic conditions. Deformable sickle RBCs in fractions I and II appear to take a preferred path (if the adjacent microgates in the flow direction are not fully blocked): they twist and turn along a serpentine path once they spot trapped cells ahead of them, see Fig 4. However, the stiff sickle RBCs in fractions III and IV just flow toward the trapped sickle RBCs and eventually stop nearby (S7 and S8 Videos). The difference in cell behavior between deformable and stiff sickle RBCs is probably due to the chaotic motion of RBCs caused by cell-cell interactions and flow perturbations near the obstructions. For a deformable sickle RBC, fluid flow can change the cell shape substantially and drift it along the flow direction; however, a stiff sickle RBC is also disturbed by the fluid flow near the obstruction.

Morphological and mechanical factors involved in vaso-occlusive crisis
SCD exhibits heterogeneous morphologies, which depend on the DeOxy rate [10]. Gradual DeOxy is known to result in predominantly elongated-and crescent-shaped RBCs, whereas rapid DeOxy results in less distorted granular-shaped RBCs. Cell morphological sickling is thus identified by visibly changes in cell shape and texture associated with transient hypoxic conditions. We performed computational simulations of patient-based samples with different cell shape count information under the Oxy and DeOxy states. We calculated the transit velocity of sickle RBCs in the microfluidic channel and compared the results to those measured in experiments (Fig 5). The results show that the RBC shape plays an important role in cell traversal through microgates: with the same cell rigidity values under the Oxy state, sickle RBCs with a disc-shape have the lowest transit velocity (v~98 μm/s). By contrast, sickle RBCs with a crescent shape have the highest transit velocity (v~136 μm/s). Therefore, the transit velocity of sickle cells with different shapes is statistically significantly different. The reason for this might stem from the different cross-sectional areas of sickle RBCs with different shapes. In our simulations, the cross-sectional area of sickle RBCs with disc, oval, and crescent shapes are about 32 μm 2 , 26 μm 2 , 23 μm 2 , respectively, which are greater than that of the microgate opening (20 μm 2 ). Thus, the sickle RBC has to deform during its traversal through the microgates. When we applied a fixed pressure gradient in each simulation to force the cells to cross the microgate, the granular-shaped ones suffer the largest deformation, resulting in the lowest transit velocity. The simulation results, albeit counter-intuitive, are quantitatively consistent with the experimental observations of sickle RBCs transiting in the microfluidic channel.
It is known that the loss of deformability of RBCs play a crucial role in vaso-occlusion in SCD [11,14,40]. We perform simulations of sickle RBCs in different density-fractionated subpopulations in order to verify the significant role of cell deformability in determining the dynamic and rheological characteristics of individual sickle RBCs. Here, we consider three distinct types of sickle RBCs under channel flow, i.e., the hypoxia-affected RBCs in fractions I (SS1), II (SS2), and IV (SS4). Our simulations indicate that sickle RBCs show increased flow resistance under the DeOxy state, especially for the sickle RBCs in fraction IV; some sickle RBCs (for example, the granular and disc-shaped ones) in this fraction always are obstructed at the microgates (Fig 5). In addition, considering the sickle RBCs in the same density-fractioned subpopulation, (hence with nearly the same shear modulus), the present results also show an obvious difference in transit velocity for sickle RBCs with different shapes. This confirms the hypothesis that SCD exhibits substantial heterogeneity even within the same densityfractionated subpopulation.
There is currently no universal cure for SCD patients, but symptoms can be managed with fluids, oxygen and medication. HU has a number of characteristics of an ideal drug for SCD patients [43,44]. It increases HbF production and reduces the occurrence of sickling-related complications [45], partially attributed to improved cell hydration [46] and cell deformability [47,48]. However, HU treatment is also associated with an elevated MCV [49,50]. This dual influence of HU on RBC structural and mechanical properties may significantly affect cell traversal through the microgates. In order to develop a quantitative assessment of the efficacy of We study the dynamic behavior of individual sickle RBCs using patient-specific hematological values for individual patients with SCD after treated with HU. When the MCV values are compared between the off-HU and on-HU groups, we find that they are always higher in the on-HU groups (see the MCV values in Table 1), consistent with previous observations [49,50]. Our simulations indicate that sickle RBCs travelling through the microgates are sensitive to MCV, and that the correlation between cell transit velocity and MCV is enhanced with the decrease in the size of the microgates. As shown in Fig 5,

Single-cell capillary obstruction
The hypoxia-affected RBCs in SCD cause blockages at the microgates under the DeOxy state. Hence, we performed simulations of patient-based samples to predict the single-cell capillary obstruction and compare the results against experimental data [23].
The capillary obstruction ratio, is defined as the ratio of the total number of trapped sickle RBCs at the microgates to the total number of RBCs in the microfluidic channel during the DeOxy state. This ratio is now examined from the four SCD blood samples including on-HU and off-HU groups. Fig 6 shows the values of capillary obstruction ratio obtained from DPD simulations and experiments. In general, the off-HU group exhibits a significantly higher capillary obstruction ratio than the on-HU group. Such increases are likely caused primarily by changes in sickled fraction. As shown in Table A in S1 Text, 20.2% and 39.0% of RBCs sickle in the off-HU/S-P-I group and off-HU/S-P-II group, respectively, among which 15.8% and 39.0% of hypoxia-affected RBCs fall within cell density fractions III and IV, respectively. In the on-HU groups, fewer RBCs become sickled than those in the on-HU groups, i.e., there are only 9.4% and 3.1% of sickled RBCs in the on-HU/S-P-III group and the on-HU/S-P-IV group, respectively, among which only 6.3% and 3.1% of hypoxia-affected RBCs fall within cell Direct observation of individual sickle cell behavior density fractions III and IV, respectively. As we demonstrated in an earlier section, the sickled RBCs show increased flow resistance under the DeOxy state and the densest ones are usually unable to traverse the microgates. Thus, a relatively large number of sickle RBCs, especially in the denser cell fractions, can cause a marked increase in single-cell capillary obstruction. The results confirm that cell deformability plays a key role in RBC traversal in small microfluidic channels. In addition, we find that the predictions of all DPD simulation cases fall within the range of experimental data. The discrepancy between predictions and experiments could arise from an underestimation of the cell morphological sickling count of individual sickle RBCs.
In summary, in this paper, we demonstrate the unique capabilities benefits of combining dynamic microfluidic experiments with multiscale simulations for characterizing the complex behavior of individual sickle RBCs in a capillary-like microenvironment under transient hypoxia. We show that hypoxia-affected RBCs undergoing sickling significantly alter cell behavior. We also monitor the dynamic behavior of both hypoxia-affected and hypoxia-unaffected RBCs as they travel through the capillary-like microenvironment under cyclic hypoxia. Taken together, these experiments and corresponding systematic particle-based simulations elucidate the effects of irregular geometry, decreased cell deformability, and elevated cell volume on the biodynamic behavior of individual sickle RBCs and their roles in single-cell capillary obstruction. They provide a quantitative measure of the heterogeneity associated with SCD at the single-cell level. Hence, the particle-based simulations and comparisons with available independent experiments offer a powerful means for real-time monitoring of in vitro behavior of individual sickle RBCs under controlled transient hypoxic conditions, providing an objective way of assessing the effectiveness of targeted drug therapy aimed at easing or preventing vasoocclusive crisis associated with SCD.
Our model does not explicitly include the intracellular HbS polymer fibers. Hence, it does not model the interaction between cell membrane and polymer fibers, and the potential influences on morphological distortion of RBCs and the attendant alteration in their mechanical properties. This deficiency could be addressed in future work by recourse to a hybrid model that encompasses the molecular and cellular scales by combing the MS-RBC model with particle-based HbS polymer models developed recently [51][52][53]. In addition, it would require further computational validation and extensive testing of these patient-specific models against clinical and experimental studies to make future predictions more reliable. Such simulations from these patient-specific predictive models would be useful for testing known biomarkers and discovering new biomarkers.