By He Li, Lu Lu, Peter Vekilov, and George Em Karniadakis
Human red blood cells (RBCs) carry oxygen from the lungs to peripheral tissues and in turn transport carbon dioxide from tissues to the lungs for exhalation. RBCs, which are approximately 8 μm wide, undergo drastic deformations during their circulation in the human vasculature as they repeatedly pass through arterioles and capillaries as small as 2-3 μm.
Individuals with sickle cell disease (SCD) carry a genetic mutation in the gene for hemoglobin, a group of oxygen-binding proteins residing inside RBCs [9]. Under hypoxia, the mutated hemoglobin (HbS) can polymerize into stiff HbS fibers that distort the RBCs into various pathological shapes, ranging from elongated, granular, oval, and holly-leaf to crescent — the classic sickle shape [2]. In addition to the heterogeneous morphological changes, sickle RBCs also exhibit increased stiffness and adhesion; these features potentially contribute to the initiation and propagation of vaso-occlusion crises (VOCs), a hallmark of SCD. Recurrent and unpredictable episodes of VOC lead to stroke, frequent painful crises, and other serious complications like acute chest syndrome, splenic sequestration crisis, and acute liver crisis [10].
SCD affects roughly 100,000 people in the U.S. and millions worldwide. The disease is most prevalent among African Americans, afflicting one out of every 365 babies. Drug therapy is still the most affordable way to treat the majority of SCD patients. For 60 years, the cancer drug Hydroxyurea was the only available medication for SCD. The U.S. Food and Drug Administration (FDA) recently approved two new drugs: Endari, which is thought to reduce oxidative stress in RBCs, and Adakveo, which seems to prevent VOCs by reducing blood cell adhesion (a downstream target in SCD). However, no molecular-level mechanisms of these two drugs have been documented in the biomedical literature.
There is a pressing need to develop new drug therapies that target the sickling of RBCs — the root cause of SCD. Discovery and development of a novel drug is a long process that normally requires four to six years of laboratory experiments on human cells and animal models before it can enter clinical trials. Given recent advances in quantitative multiscale modeling for biological systems, it is now possible to integrate traditional drug development approaches with in silico methods to accelerate laboratory studies and expeditiously explore new SCD drug targets.
Mathematical Models
Our group has generated particle-based microscopic and mesoscopic models to simulate biological problems in SCD and improve our understanding of the disease’s pathogenesis and pathophysiology. At the microscopic level, we created coarse-grained molecular dynamics (CGMD) HbS fiber models [6, 7] to simulate HbS molecules’ self-assembly into fibers and the consequent fiber-fiber interaction. In parallel work, we developed OpenRBC [13], a CGMD code that uses a single shared memory commodity workstation to simulate an entire RBC at the molecular level. The simulation includes multiple millions of particles that explicitly represent the lipid bilayer and cytoskeleton. At the mesoscopic level, we produced dissipative particle dynamics (DPD) models of RBCs based on the rigorous theory of the Mori-Zwanzig formalism from statistical physics [5]. DPD-based RBC models enable us to simulate biological processes beyond the cellular level, such as the initiation of VOC by sickle RBCs in capillaries and venules [3].
We recently put forth a kinetic model [8] based on classical nucleation theory to explicitly describe the kinetics of HbS molecules’ homogeneous and heterogeneous nucleation, as well as the growth dynamics of HbS fibers. Our team built a computational framework by integrating this kinetic model with our mechanistic model to assess potential SCD drugs and individualized pathology (see Figure 1).
Figure 1. A computational framework that integrates a kinetic model with a mechanistic model can predict the efficacy of anti-sickling drugs and discover new drug targets. 1a. The kinetic model predicts the number of fibers inside red blood cells (RBCs) and the rate at which they grow under patient-specific and organ-specific conditions. 1b. The mechanistic model determines the sickling of RBCs based on the kinetic model’s outputs. Figure adapted from [8]. Refer to [8] for definitions of kinetic quantities in 1a.
In Silico Analysis for Enhanced Drug Development
We validated this new framework by demonstrating that our model predictions are consistent with prior in vivo [11] and in vitro studies [1, 4]. In particular, our model can predict the SCD drug efficacy in agreement with the drug screening assay in Figures 2a and 2b, where each point of the experimental data (black symbols) represents 12 measurements with durations of either four or 24 hours [4] — not including experimental setup time. On the other hand, it only takes about one minute to obtain the corresponding simulation results (blue symbols) by running the computational model on a laptop computer. Given this high efficiency, researchers can use computational modeling as an extra drug pre-screening modality to reduce the number of experimental tests without compromising measurement reliability. This accelerates laboratory studies before new drugs enter expensive clinical trials.
Figure 2. Selection of critical parameters via global sensitivity analysis. 2a and 2b. The fraction of sickled red blood cells (RBCs) decreases with increased dosage of anti-sickling agents monensin A (2a) and gramicidin A (2b). Blue curves represent our simulation results and black curves denote the experimental results reported in [4]. Squares and diamonds correspond to four-hour incubations, and triangles correspond to 24-hour incubations. 2c and 2d. Global parametric sensitivity analysis on the kinetic quantities with model inputs following an in vitro study (2c) [4] and an in vivo study (2d) [11]. The diameter of the green circles is the first-order sensitivity index of each parameter. The blue ring surrounding the green circles represents the sensitivity index’s 95 percent confidence interval. The line thickness between two parameters indicates the second-order sensitivity index of each two-parameter pair. Figure adapted from [8].
We performed a global sensitivity analysis of sickled fractions of RBCs to the model’s kinetic quantities to identify key quantities for consideration as potential druggable targets. We employed variance functional decomposition [12]—in which we simultaneously varied 13 quantities—to conduct this analysis. As shown in Figures 2c and 2d, our global sensitivity analysis indicates that RBC sickling is sensitive to hemoglobin solubility \((C_{e0})\), nucleation rate pre-factor \((J_0)\), oxygen-binding affinity \((K_R)\), and hemoglobin activity \((B_2)\). These findings provide therapeutic rationales for anti-sickling targets, which are consistent with previous findings but also identify new druggable targets. For example, a new anti-sickling drug called Voxelotor—recently approved by the FDA for an accelerated clinical trial—functions by increasing the oxygen-binding affinity, one of four targets identified by our model. Furthermore, our model can provide quantitative guidance on the required modification levels of these four kinetic quantities to achieve therapeutic efficacy.
The clinical outcomes of SCD drugs are heterogeneous [14]. This heterogeneity may be due to the variation in patient-specific RBC variables—like hemoglobin distribution and composition, and mean cell volume—as well as organ-specific environmental variables such as temperature, deoxygenation time, and initial and final oxygen pressure. Our computational framework enables predictions of RBC sickling using both patient- and organ-specific data, which opens the path to monitoring disease progression for individual SCD patients with various degrees of severity. Such model predictions can help clinicians make patient-specific prognoses and offer guidance on the optimal time for medical intervention. Our kinetic model can also provide direction for drug dosage based on patient-specific data, thus facilitating the development of precision medicine to eliminate the side effects of SCD drugs.
Acknowledgments: This work was supported by National Institutes of Health grant U01HL114476. The authors thank Ming Dao at the Massachusetts Institute of Technology; Pankaj Qasba at the National Heart, Lung, and Blood Institute; and Grace Peng at the National Institute of Biomedical Imaging and Bioengineering for helpful discussions.
References
[1] Du, E., Diez-Silva, M., Kato, G.J., Dao, M., & Suresh, S. (2015). Kinetics of sickle cell biorheology and implications for painful vasoocclusive crisis. Proc. Nat. Acad. Sci., 112(5), 1422-1427.
[2] Ferrone, F.A. (2004). Polymerization and sickle cell disease: A molecular view. Microcirc., 11(2), 115-128.
[3] Lei, H., & Karniadakis, G.E. (2013). Probing vasoocclusion phenomena in sickle cell anemia via mesoscopic simulations. Proc. Nat. Acad. Sci., 110(28), 11326-11330.
[4] Li, Q., Henry, E.R., Hofrichter, J., Smith, J.F., Cellmer, T., Dunkelberger, E.B., …, Eaton, W.A. (2017). Kinetic assay shows that increasing red cell volume could be a treatment for sickle cell disease. Proc. Nat. Acad. Sci., 114(5), E689-E696.
[5] Li, Z., Bian, X., Caswell, B., & Karniadakis, G.E. (2014). Construction of dissipative particle dynamics models for complex fluids via the Mori-Zwanzig formulation. Soft Matt., 10(43), 8659-8672.
[6] Lu, L., Li, H., Bian, X., Li, X., & Karniadakis, G.E. (2017). Mesoscopic adaptive resolution scheme toward understanding of interactions between sickle cell fibers. Biophys. J., 113(1), 48-59.
[7] Lu, L., Li, X., Vekilov, P.G., & Karniadakis, G.E. (2016). Probing the twisted structure of sickle hemoglobin fibers via particle simulations. Biophys. J., 110(9), 2085-2093.
[8] Lu, L., Li, Z., Li, H., Li, X., Vekilov, P.G., & Karniadakis, G.E. (2019). Quantitative prediction of erythrocyte sickling for the development of advanced sickle cell therapies. Sci. Advan., 5(8), eaax3905.
[9] Pauling, L., Itano, H.A., Singer, S.J., & Wells, I.C. (1949). Sickle cell anemia, a molecular disease. Science, 110(2865), 543-548.
[10] Piel, F.B., Steinberg, M.H., & Rees, D.C. (2017). Sickle cell disease. New Eng. J. Med., 376(16), 1561-1573.
[11] Serjeant, G.R., Petch, M.C., & Serjeant, B.E. (1973). The in vivo sickle phenomenon: A reappraisal. J. Lab. Clin. Med., 81(6), 850-856.
[12] Sobol, I.M. (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Sim., 55(1-3), 271-280.
[13] Tang, Y.-H., Lu, L., Li, H., Evangelinos, C., Grinberg, L., Sachdeva, V., & Karniadakis, G.E. (2017). OpenRBC: A fast simulator of red blood cells at protein resolution. Biophys. J., 112(10), 2030-2037.
[14] Ware, R.E. (2015). Optimizing hydroxyurea therapy for sickle cell anemia. ASH Edu. Prog. Book, 2015(1), 436-443.
He Li is a research assistant professor in the Division of Applied Mathematics at Brown University. Lu Lu is a senior Ph.D. and master’s student in computer science in the Division of Applied Mathematics at Brown. Peter Vekilov is a John and Rebecca Moores Professor of Chemical and Biomolecular Engineering and of Chemistry at the University of Houston. He is a fellow of the American Physical Society. George Em Karniadakis is the Charles Pitts Robinson and John Palmer Barstow Professor of Applied Mathematics at Brown. He is a SIAM Fellow.