Sickle cell disease is induced by a mutation that converts normal adult hemoglobin to sickle hemoglobin (HbS) and engenders intracellular polymerization of deoxy-HbS and erythrocyte sickling. Development of anti-sickling therapies requires quantitative understanding of HbS polymerization kinetics under organ-specific conditions, which are difficult to assess with existing experimental techniques. Thus, we developed a kinetic model based on the classical nucleation theory to examine the effectiveness of potential anti-sickling drug candidates. We validated this model by comparing its predictability against prior in vivo and in vitro experimental results. We used the model to quantify the efficacy of sickling inhibitors and obtain results consistent with recent screening assays. Global sensitivity analysis on the kinetic parameters in the model revealed that the solubility, nucleation rate prefactor, and oxygen affinity are quantities that dictate HbS polymerization. This finding provides quantitative guidelines for the discovery of intracellular processes to be targeted by sickling inhibitors.

( A ) Schematic of the inputs and outputs of the kinetic model (top), which describes the nucleation, growth, and branching of HbS fibers and RBC sickling (bottom). ( B ) Concise flow chart of prediction of RBC sickling through the kinetic model. The model inputs are temperature T, oxygen tension P O 2 , deoxygenation rate, HbS/HbA/… mole fraction X s /X A /…, total hemoglobin concentration C t , and RBC volume V c . The kinetic model computes the following key quantities: the fraction of oxyhemoglobin Θ, hemoglobin solubility C e and supersaturation Δμ, the nucleation delay time τ d , homogeneous and heterogeneous nucleation rates J and J s , fiber growth rate R, fiber length, and the RBC sickling time. The theories or empirical relations used to compute these kinetic quantities are given between each step. The parameters in the kinetic model are highlighted in the parentheses at each step.

Substantial progress has been made in elucidating the kinetics of HbS polymerization and RBC sickling ( 3 , 9 – 11 ), which provided insights into the development of anti-sickling drugs. However, few of the prior studies focus on testing the efficacy of the SCD drugs. In this work, we connect fundamental physical principles of HbS polymerization at the molecular level with clinical data of RBC sickling by integrating a kinetic model and a mechanistic model, aiming to build a detailed framework that could be used to assess potential SCD drugs and individualized pathology. To the best of our knowledge, this framework uses nucleation theory to explicitly describe the kinetics of homogeneous and heterogeneous nucleation of HbS monomers and chemical rate laws to model the growth dynamics of HbS fibers to capture the complexity of SCD. These quantities are not considered explicitly in a previous model ( 11 ), which was developed on the basis of numerical fitting of the correlation between RBC sickling delay time and supersaturation. In addition, the proposed kinetic model can predict the fraction of sickled RBCs under organ-specific conditions, which are difficult to assess through in vivo studies of animal models or in vitro microfluidic experiments. An important feature of this proposed framework is that each of its components can be refined by future experimental or analytical studies to improve the accuracy and capability of the model predictions. We also perform sensitivity analysis of the kinetic variables implemented in the model, through which we identify the quantities that can be considered as potential druggable targets. Specifically, the new model, as shown in Fig. 1 , can predict the fraction of sickled RBCs according to patient-specific inputs and organ-specific environmental conditions and examine the therapeutic efficacy of anti-sickling agents. Furthermore, it can quantify the impact of multiple polymerization factors in the process of HbS polymerization and identify the key governing parameters that can be targeted by drugs to efficiently prevent RBC sickling.

Therapeutic treatments have been developed on the basis of increased understanding of the pathophysiology of SCD ( 5 ). Transplantation of stem cells from a matched sibling donor has been used to cure SCD patients with high severity. Other curative treatments, such as gene replacement and gene editing, are at the stage of clinical trials. These approaches require advanced medical facilities associated with high costs, which render them impractical for most of the SCD patients. Drug therapy, therefore, is still the most feasible way to treat SCD patients. However, hydroxyurea was the only available SCD drug in the past 60 years before Endari was approved in July 2017. Hydroxyurea inhibits RBC sickling by diluting intracellular HbS through boosting the content of fetal hemoglobin (HbF), a type of hemoglobin that does not assemble into fibers ( 6 ). Endari is thought to reduce the oxidative stress in RBCs, but the underlying mechanism is still not well understood ( 7 ). Currently, most of the new drugs in clinical trials target either the adhesion of blood cells (RBCs, leukocytes, neutrophils, and platelets) to endothelium or vascular inflammation, which are downstream consequences of HbS polymerization ( 8 ). Few drug candidates aim at inhibiting HbS polymerization and sickling of RBCs, the root cause of SCD.

Sickle cell disease (SCD), a genetic blood disorder, originates from the mutation of a single nucleotide in the gene for hemoglobin, a group of oxygen-binding proteins residing in the red blood cells (RBCs; erythrocytes) ( 1 ). This mutation enables polymerization of HbS molecules into polymers under hypoxia. Polymerization of HbS molecules is initiated by deoxygenation and the associated conformational change in hemoglobin from R (relaxed) state and the T (tense) state ( 2 ). HbS molecules rapidly aggregate once a nucleus forms, leading to the growth of fibers. The subsequent branching of fibers has been cast as double nucleation mechanism ( 3 ). The growing HbS fibers distort RBCs and impair cell deformability, potentially contributing to the initiation and propagation of vaso-occlusion (VOC), a hallmark of SCD. Recurrent and unpredictable episodes of VOC lead to stroke, frequent painful crisis, and other serious complications, including acute chest syndrome, splenic sequestration, and aplastic crises ( 4 ).

RESULTS

We propose a predictive model of RBC sickling by integrating a stochastic kinetic model, which describes the homogeneous and heterogeneous nucleation of HbS polymer fibers and their growth, with a mechanistic model, which considers the interaction between HbS fibers and the RBC membrane. The driving force for HbS fiber nucleation and their growth is the supersaturation (12), i.e., the excess of deoxy-HbS chemical potential in the cytosol over that of HbS in the polymer. Prior studies have introduced an all-encompassing delay time between the setting of supersaturation and sickling of RBCs and used an empirical scaling relation to evaluate this time (3, 11). Here, we use a fundamentally different strategy. We represent the nucleation and growth of HbS fibers with established biophysical kinetic laws (9, 13), which allows for a detailed analysis of initiation of RBC sickling. As shown in Fig. 1A, the inputs of the kinetic model include environmental variables, such as temperature, deoxygenation time, and oxygen tension (the partial pressure of oxygen in blood), and patient-specific variables of RBCs, such as mean corpuscular volume (MCV; the average volume of RBCs), hemoglobin distribution, and concentration of other hemoglobin variants. The outputs are numbers of nuclei generated through homogeneous and heterogeneous pathways, respectively, the lengths of HbS fibers, the configuration of fiber domains, and the sickling of RBCs. Detailed information of the new stochastic kinetic model can be found in Materials and Methods.

Mechanics of RBC sickling in SCD Sickling of RBCs is induced by the growth of intracellular HbS fibers, which distort RBCs into a variety of shapes, including elongated, granular, oval, holly leaf, and crescent (classic sickle) shapes (14). However, a quantitative correlation between the mechanical properties of HbS fibers and the resulting morphological changes of RBCs has not been investigated in sufficient detail. Here, we combine a coarse-grained molecular dynamics (CGMD) HbS fiber model (15) with a CGMD RBC model (16) to simulate the interaction between an RBC and HbS fibers, through which we clarify the relation between the stiffness of HbS fibers and the deformations of RBCs. The RBC model simulates the lipid bilayer and the cytoskeleton as well as transmembrane proteins explicitly by coarse-grained (CG) particles (Fig. 2A). The HbS fiber model is an intracellular fiber bundle with a chain of CG particles. The diameter of the fiber chain can be tuned to mimic HbS fiber bundles that contain various numbers of single HbS fibers. The bending stiffness of the fiber bundle with an increased number of constituent fibers (blue curve in Fig. 2B) is tuned such that it follows the experimental data reported in (17) (black curve in Fig. 2B). Fig. 2 Simulations of RBC sickling with a mechanistic model. (A) Configuration of an HbS fiber bundle in a stretched RBC. Green particles, actin junctions; purple particles, spectrin filaments; black particles, glycophorin proteins; yellow particles, band-3 proteins; red particles, lipid bilayer. Yellow bar inside the RBC represents a chain of CG HbS particles with varying thickness, mimicking an HbS fiber bundle. (B) Bending stiffness of fiber bundles plotted as a function of the number of constituent fibers. Bending stiffness of fiber bundle dictates the variations in the ensuing fiber bundle and RBC shapes, as illustrated in the plot: (I) fiber buckled, (II) fiber buckled with one protrusion on the RBC, (III) fiber buckled with two protrusions on the RBC, (IV) nonbuckled with two protrusions on RBC. Fibers are plotted with an exaggerated thickness for clarity. In the insets, a part of the RBC membrane is not shown such that the HbS fiber configuration inside the RBC can be viewed. To probe the correlation between the mechanical properties of HbS fibers and the resulting morphological changes of RBCs, we initially placed a fiber bundle with a length of 15 μm in a prestretched RBC (Fig. 2A), following the strategy used in (18). Once the simulation begins, the RBC attempts to recover the biconcave shape and interacts with the fiber bundle. The final configurations of the RBC and fiber bundle are dictated by the bending stiffness of the fiber bundle. Figure 2B illustrates that when the fiber bundle contains equal to or less than 13 single HbS fibers, the fiber bundle becomes convoluted after interacting with the RBC membrane. When the number of fibers increases to 14, the fiber bundle still becomes buckled, but it creates one protrusion on the cell surface. As the number of fibers in the bundle increases, the fiber bundle becomes less buckled due to the increased bending stiffness. Meanwhile, two growing protrusions are observed on the surface of the RBC. Notably, these simulation results are consistent with an analytical study, which demonstrated that an intracellular fiber bundle, consisting of 10 or more tightly bonded single HbS fibers, sickles an RBC (19). Therefore, in the kinetic model, RBCs are considered to be “sickled” when an intracellular bundle, consisting of 10 or more HbS fibers, reaches the width of an erythrocyte (6 μm or longer). A parametric study for these criteria is shown in the Supplementary Materials.

Sickling of RBCs in the hepatic vein The hepatic veins are the location with the lowest oxygen saturation of 51% in the venous circulation, corresponding to P O 2 = 28.1 mmHg (20). To validate the model prediction of RBC sickling in hepatic vein, we compare the simulation results with experimental measurements conducted in (20). Since the distribution of hemoglobin concentration and composition was not reported in the original work, we assume that hemoglobin distribution follows a Gaussian distribution with a mean of 32.7 g/dl and an SD of 2.3 g/dl (Fig. 3B) (21, 22). We also assume that the examined RBCs contain only one hemoglobin variant, HbS. In addition, P O 2 is assumed to decrease linearly from 100 to 28.1 mmHg. The deoxygenation time is selected to be 1 s, comparable to the average capillary transit time of RBCs (23). Other model inputs can be found in table S3 (condition I). Fig. 3 Deoxygenation and fraction of sickled RBCs in hepatic vein. (A) The linear decrease of P O 2 from 100 to 28.1 mmHg in 1 s (dashed line) enforces deoxygenation of hemoglobin from 4.3 to 49% correspondingly (blue curve). Fractions of RBCs classified according their sickling propensities in the hepatic vein of (B) SCD patients and (C) SCT individuals correspondingly. Areas A, B, C, and D denote fractions of RBCs with distinct polymerization propensities at P O 2 = 28.1 mmHg. Area A: Deoxy HbS is undersaturated. Area B: The nucleation delay time is longer than 10 s. Area C: Less than 1 polymer is generated via homogeneous nucleation within 10 s. Area D: RBCs are likely to become sickled within 10 s. (D to H) Variation of supersaturation Δμ/k B T (D), delay time τ d (E), homogeneous nucleus size n* (F), homogeneous nucleation barrier ΔG* (G), and homogeneous nucleation rate J (H) with respect to the total hemoglobin concentration C t at the end of deoxygenation. Vertical dash line highlights the beginning of supersaturation. As shown in Fig. 3A, the concentration of fully deoxygenated hemoglobin in T-state conformation, which is determined by the total hemoglobin concentration C t and oxygen tension P O 2 , increases to 0.49C t at the end of deoxygenation period. The propensity of the deoxy-HbS to polymerize and induce RBC sickling is dictated by the value of C t , e.g., RBCs with higher C t are more likely to sickle. Figure 3 (D to H) present the dependence of the supersaturation, nucleation delay time, nucleus size, nucleation barrier, and nucleation rate on the hemoglobin concentration in RBCs during their passage through hepatic vein. On the basis of these results, the RBCs are divided into fractions with different sickling probabilities, as shown in Fig. 3B. When C t is less than a threshold value of 26.4 g/dl (area A in Fig. 3B), the red cell cytosol is undersaturated after deoxygenation (Δμ/k B T < 0) (see Fig. 3D), and polymerization is thermodynamically prohibited. For RBCs with C t larger than 26.4 g/dl, the intracellular hemoglobin becomes supersaturated, and nucleation proceeds after a delay time of τ d , a fundamental characteristic of nucleation, defined as the length of the transition from a stable homogeneous solution to the moment when nucleation proceeds at a steady rate (9, 13). As shown in Fig. 3E, this delay time could range from hundreds of seconds to 0 s, depending on C t . RBCs with C t less than 30.0 g/dl (area B in Fig. 3B) are likely to return to the lungs without sickling, due to the long delay time (>10 s) (see Fig. 3D). When C t is larger than 30.0 g/dl, the delay time is shortened to several seconds, comparable to RBC capillary transit time (23). However, the probability of sickling is small for RBCs with C t < 34.3 g/dl (area C in Fig. 3B) because the corresponding homogeneous nucleation rate J is smaller than 1.07 × 109 cm−3 s−1, under which fewer than one nucleus is generated in 10 s (see Fig. 3H). Hence, sickling of RBCs is more likely to occur for those RBCs with C t > 34.3 g/dl (area D in Fig. 3B), which accounts for ~24.3% of RBCs. This result is close to the fraction of sickled RBCs (22%) observed in the hepatic vein of the SCD patients examined in (20), demonstrating the power of the kinetic model to predict the fraction of sickled RBCs in vivo. To demonstrate the robustness of the model, we also computed the fraction of sickled RBCs in hepatic vein of sickle cell trait (SCT) individuals. Unlike SCD patients who have two mutated genes that cause the production of HbS, SCT individuals carry only one mutated gene, and thus have a lower percentage of HbS in RBCs compared to SCD patients. People with SCT typically live a normal life and only experience symptoms of SCD under extreme conditions such as severe dehydration and high-intensity physical activities. In our simulation, the distribution of hemoglobin in SCT RBCs is selected to be identical to that of SCD patients, but the hemoglobin composition is 38% HbS and 62% Hemoglobin A (HbA), consistent with the data reported in (22). The reduced HbS percentage increases the solubility of hemoglobin and thus reduces the propensity for polymerization. As shown in Fig. 3C, the threshold values of C t for Δμ/k B T < 0, τ d > 10 s, and JV c < 0.1 s−1 increased. Note that the population of RBCs with the largest propensity to sickle (area D in Fig. 3C) only accounts for less than 0.03% of the RBCs. This finding provides a rationale for the clinical observation that, under normal conditions, individuals with SCT do not experience symptoms or complications similar to those of SCD patients.

Patient-specific studies of RBC sickling under hypoxia in microfluidic channel In this section, we compute the fractions of sickled RBCs from seven SCD patients under the experimental conditions of Du et al. (24). In this in vitro study, P O 2 in the microfluidic channel was reduced from 160 mmHg (O 2 concentration of 21%) to 15.2 to 38 mmHg (O 2 concentration of 2 to 5%) in 15 s. The compositions of hemoglobin, MCV, mean corpuscular hemoglobin concentration (MCHC) for the seven SCD patients are listed in table S1. In our simulations, we assume that the distribution of C t in the examined RBCs follows a Gaussian distribution with a mean of MCHC (listed in table S1) and a universal variance of 2.3 g/dl, as the actual variances for each patient were not provided in the original work (24). Other model inputs can be found in table S3 (condition II). First, we compute the fractions of sickled RBCs, assuming that P O 2 is reduced to 15.2 mmHg (2%) after deoxygenation, the lower bound of O 2 concentration range reported in the microfluidic experiments (24). The fractions of sickled RBCs are recorded in the subsequent 240 s after deoxygenation begins. As shown in Fig. 4, significant discrepancies are observed between our simulation results (red dashed curve) and experiential data (black curve) (24) for patient VII, which probably results from the underestimated P O 2 and application of the universal variance for hemoglobin distribution. To confirm our hypothesis, we vary O 2 concentration within 2 to 5%, the range reported in the original experiments (24). Meanwhile, we tune the variance of hemoglobin distribution such that our simulation results gradually approach the experimental results (black curve in Fig. 4). These tests indicate that when O 2 concentration is 3.75% and the variance is 0.7 g/dl, our simulation results (blue dashed curve in Fig. 4) are in agreement with the experimental data. Similar computations are performed for the other six SCD patients, and the results can be found in fig. S5. Fig. 4 Evolution of fraction of sickled RBCs under hypoxia. Deoxygenation completes in 15 s after the simulations begin. Black line represents the experimental results reported in (24). Red dashed line denotes model prediction with O 2 concentration of 2% and variance of HbS distribution of 2.3 g/dl. Blue dashed curve denotes model prediction with O 2 concentration of 3.75% and variance of HbS distribution of 0.7 g/dl. Gray lines represent 20 of 5000 sensitivity tests with O 2 concentration varying at 3.75 ± 0.01% and the variance of hemoglobin distribution varying between 0.7 ± 0.35 g/dl. Blue area represents SD. An additional 5000 simulations for each SCD patient (20 of 5000 sensitivity tests for patient VII are plotted as gray lines in Fig. 4) are performed to examine the sensitivity of the model prediction to O 2 concentration and the variance of hemoglobin distribution. We find that a fraction of sickled RBCs is more sensitive to variation of O 2 concentration. However, we note that the agreement between the simulation and experimental results cannot be achieved by varying the P O 2 alone, suggesting that patient-specific hemoglobin distribution is an essential determinant of the transient sickling of RBCs. The above findings demonstrate the capability of the kinetic model to predict the RBC sickling based on patient-specific inputs.

Examination of sickling inhibitors Although the molecular basis of SCD was discovered about 70 years ago (1), only a single drug, hydroxyurea, was available to SCD patients until Endari was approved in 2017. Hydroxyurea inhibits HbS polymerization by reducing the fraction of HbS in the red cell cytosol through inducing the synthesis of HbF (6). However, clinical evidence indicates that hydroxyurea exhibits heterogeneous outcomes and that the causes of these findings are still under investigation (25). An alternative approach to inhibit RBC sickling is to increase cell volume using compounds that interact with the cell membrane. Two compounds of this class, monensin A and gramicidin A, have demonstrated strong therapeutic effects in recent screening assays (22). In this section, we apply our model to test the suggested mechanism of action of sickling inhibitors. First, we use the kinetic model to evaluate the sickled fraction of SCT cells before application of a drug. We select a hemoglobin composition of 38% HbS and 62% HbA, comparable to the composition of the SCT subjects (38% HbS, 58% HbA, 3.7% HbA2, and 0.3% HbF) reported in (22). We use the hemoglobin concentration distribution for normal RBCs reported in (21). Following the original experiments in (22), the initial P O 2 is set to 0 mmHg. Other model inputs can be found in table S3 (condition III). Our simulation predicts a wide distribution of nucleation delay times primarily due to variation in MCHC. As illustrated in Fig. 5A (blue dashed curve), this wide distribution of nucleation delay time leads to a gradually increased fraction of sickled RBCs with time. Notably, the prediction from our kinetic model is consistent with experimental data reported in (black curve in (22) (black curve in Fig. 5A). Moreover, when the deoxygenation time is shorter than 20 s, our model (blue dashed curve) provides a more accurate prediction of the fraction of sickled RBCs than a prior model (red dashed curve) based on the universal curve of delay time versus supersaturation (11). This discrepancy is likely attributed to the fact that the present model considers the number of HbS molecules in the nucleus n* as discrete integer numbers, whereas the model applied in (22) represents n* as a continuous variable. Continuous representation of n* leads to infinitely high HbS polymerization rates when the delay time is zero, and thus causes an overestimation of the fraction of sickled RBCs at the early stage of measurements (see red dashed curve in Fig. 5A). Fig. 5 Predictions of the efficacy of sickling inhibitors with the kinetic model. (A) Evolution of the fraction of sickled RBCs in the absence of sickling inhibitors. Black line, experimental results reported in (22); red dashed line, the prediction based on the universal curve of delay time versus supersaturation in (22); blue dashed line, prediction of the kinetic model. (B and C) The fraction of sickled RBCs decreases with increased MCV at 60 s after photolysis. The swelling of RBCs is induced by (B) monensin A and (C) gramicidin A. Blue curves represent our simulation results, and black dashed curves denote the experimental results reported in (22). Squares and diamonds denote 4-hour incubations. Triangles denote 24-hour incubations. (D and E) Sickled fractions at 60 after photolysis of SCT RBCs treated with (D) AES103 and (E) GBT440 at a variety of concentrations. Blue curves represent our simulation results, and black curves denote the experimental results reported in (22). Square and diamond symbols denote 4-hour incubations. Triangle symbols denote 24-hour incubations. logarithm bas 10, Log 10 . Next, we examined the efficacy of two anti-sickling agents, monensin A and gramicidin A, both of which operate by reducing the concentration of HbS through increasing the cell volume. To evaluate the increased RBC volume with various dosages of these two drugs, we followed the Coulter counter measurements performed in (22). As shown in Fig. 5B, the kinetic model predicts that as the concentration of the monensin A increases, the sickled fraction decreases in the same pattern as the results of the screening assay for 4- and 24-hour incubations, respectively. This finding demonstrates the capability of the kinetic model in predicting the efficacy of SCD drugs. For gramicidin A, our predictions for 4-hour incubation are close to the results of the screening assay (Fig. 5C). However, the simulation results for 24-hour incubation correspond to a stronger efficacy than the experimental results at the lower end of the examined concentration range (Fig. 5C). This discrepancy suggests that the correlation between the increased cell volume and decreased sickled fraction for gramicidin A–treated RBCs is nonmonotonic (22). Apparently, besides inhibiting RBC sickling through increasing RBC volume, gramicidin A has additional mechanisms of promoting RBC sickling, which is overwhelmed by the effect of increased cell volume as the concentration of gramicidin A increases. Last, we examine the efficacy of another two SCD drugs, AES103 and GBT440, both of which are proposed to inhibit RBC sickling by increasing the oxygen affinity. Following the screening assay in (22), the efficacy of these two drugs is tested in the absence of oxygen. We use condition III in table S3 for model inputs. Since the quantitative correlations between the concentrations of these two drugs and their effects on the oxygen binding parameters, such as c, K R , and L, are not available, we only consider their effects on changing the cell volume, as reported in (22). As shown in Fig. 5D, AES103 exhibits no anti-sickling effect in the absence of oxygen because it does not increase cell volume, consistent with the finding of the screening assay (22). For GBT440, our simulations in Fig. 5E show that after a 4-hour incubation, no considerable anti-sickling effect is observed at a low concentration, but the effect becomes significant at a high concentration, in agreement with the screening assay (22). However, in the case of 24-hour incubation, our simulation predicts notable sickling inhibition by GBT440 at low concentration due to the cell volume change, which was not observed from the screening assay (22). This discrepancy implies that there are the additional mechanisms of GBT440 that cancel out the effect of cell volume changes on inhibiting RBC sickling at low concentration. As the concentration of GBT440 is increased, both our simulation and the screening assay predict a therapeutic inhibition due to the increased cell volume.