Signaling-biophysical modeling unravels mechanistic control of red blood cell phagocytosis by macrophages in sickle cell disease

Abstract

Red blood cell (RBC) aging manifests through progressive changes in cell morphology, rigidity, and expression of membrane proteins. To maintain the quality of circulating blood, splenic macrophages detect the biochemical signals and biophysical changes of RBCs and selectively clear them through erythrophagocytosis. In sickle cell disease (SCD), RBCs display alterations affecting their interaction with macrophages, leading to aberrant phagocytosis that may cause life-threatening spleen sequestration crises. To illuminate the mechanistic control of RBC engulfment by macrophages in SCD, we integrate a system biology model of RBC-macrophage signaling interactions with a biophysical model of macrophage engulfment, as well as in vitro phagocytosis experiments using the spleen-on-a-chip technology. Our modeling framework accurately predicts the phagocytosis dynamics of RBCs under different disease conditions, reveals patterns distinguishing normal and sickle RBCs, and identifies molecular targets including Src homology 2 domain-containing protein tyrosine phosphatase-1 (SHP1) and cluster of differentiation 47 (CD47)/signal regulatory protein α (SIRPα) as therapeutic targets to facilitate the controlled clearance of sickle RBCs in the spleen.

Introduction

More than 2.5 million red blood cells (RBCs) are generated every second by erythropoiesis in the human body, each with an average lifespan of ∼120 days (1). To maintain homeostasis, aged and damaged RBCs are constantly cleared from circulation as new ones are generated (2). The spleen is the primary organ for the filtration and clearance of RBCs (3). In the red pulp of the spleen, traversing RBCs interact with splenic macrophages where normal RBCs are allowed to move toward the interendothelial slits (IESs) for biomechanical filtration while aged or damaged RBCs are retained and eliminated via erythrophagocytosis (4, 5). During the interaction between RBCs and macrophages, ligands exposed on the surface of the RBCs bind to the receptors on the surface of macrophages and modulate downstream signaling that controls phagocytosis (6).

These signaling interactions can be broadly categorized into “eat me” signals that promote phagocytosis and “do not eat me” signals that inhibit it (5). The exposure of phosphatidylserine (PS) on the surface of RBCs increases as the cells age, serving as one of the “eat me” signals (7, 8). Exposed PS molecules on RBCs bind to PS-recognizing receptors on the surface of macrophages, including integrins, PS receptors (PSRs) that include transmembrane immunoglobulin and mucin domain proteins (TIMs), and several other receptor classes, and activate downstream signaling that promotes phagocytosis (5, 7). Band-3 is another highly expressed transmembrane protein on RBCs that was shown to be a target of opsonization by naturally occurring antibodies (NAbs). Band-3 binding with NAbs can activate Fc-gamma receptors (FcγRs) on macrophages and thus serves as another “eat me” signal to stimulate phagocytosis (5, 9). In addition, the rigidity of RBCs also acts as an “eat me” signal and stimulates the activation of myosin in macrophages (10). On the other hand, the cluster of differentiation 47 (CD47) proteins, sometimes referred to as integrin-associated proteins (IAPs), is highly expressed on the surface of RBCs and serves as a self-marker and “do not eat me” signal (11). CD47 binds to signal regulatory protein α (SIRPα) receptors on the surface of macrophages and inhibits the activation of nonmuscle myosin IIa through Src homology 2 domain-containing protein tyrosine phosphatase-1 (SHP1), thereby attenuating phagocytosis signaling (12). Markedly, CD47 in aged RBCs was shown to undergo conformational changes, allowing it to bind to thrombospondin-1 (TSP-1), a matricellular protein. This activates alternate downstream signaling through SIRPα that instead promotes phagocytosis, suggesting that CD47 is both a potential “eat me” signal in aged RBCs and a “do not eat me” signal in healthy RBCs (13).

The signaling control of erythrophagocytosis plays an important role in both physiological homeostasis and many RBC-related pathological conditions. As RBCs age and accumulate physical and molecular alterations, the exposure of PS increases (8), and CD47 proteins undergo conformational changes that allow them to bind TSP-1 and switch to an “eat me” signal (13). In addition, the deformability of RBCs decreases with age (14). These changes allow senescent RBCs to preferentially activate engulfment signaling in macrophages, making them more prone to phagocytosis than young RBCs (13). In pathological conditions that cause excessive oxidative damage to RBCs, such as glucose-6-phosphate dehydrogenase (G6PD) deficiency and malaria, damaged RBCs can stimulate their own phagocytosis by upregulating band-3 (15). In sickle cell disease (SCD), an inherited blood disorder that affects the shape, stiffness, and function of RBCs (16), RBCs display many physical and molecular changes that impact their interaction with splenic macrophages. Their decreased deformability and increased expression of adhesion molecules make them more likely to be retained by macrophages (17). Studies have shown that sickle RBCs are more susceptible to erythrophagocytosis in vivo (16, 18). The surface expression levels of the signaling ligands on sickle RBCs are also changed compared to those on normal RBCs. PS exposure was shown to be elevated in sickle RBCs by up to 15-fold (19, 20). The expression of band-3 slightly decreases in sickle RBCs (21). Compared to healthy RBCs, sickle RBCs have also been demonstrated to express less CD47, which can further enhance their pro-phagocytic signaling when interacting with macrophages (22, 23). Sickle RBCs also have decreased deformability compared to healthy RBCs (24), likely leading to their retention by IESs in the spleen (25, 26), the most difficult challenge to RBC deformability in human circulation. When the retained sickle RBCs cannot be phagocytosed efficiently, a life-threatening acute splenic sequestration crisis (ASSC) may develop and cause splenic dysfunction, severely endangering the life of the patient (27, 28). Computational models that focus on the mechanics of erythrocyte sickling alone and the resulting aggregation in the spleen have been developed to advance our understanding of this pathological condition (25, 29). However, these models have not yet considered erythrophagocytosis or how signaling interactions between RBCs and macrophages control downstream cell fate; therefore, integrative approaches with advanced multiscale physiological details are needed to better address this question.

In this study, we develop a multiphysics computational framework to formulate a mechanistic and quantitative understanding of the signaling control of erythrophagocytosis and apply it to the context of sickle cell disease (SCD) to characterize the key mechanisms modulating engulfment dynamics. In this hybrid framework, we build a multiscale mechanism-based mathematical model to describe not only the molecular level details of various “eat me” and “do not eat me” signal transduction pathways, but also the signaling-driven macrophage engulfment of RBCs and their individual cell fates under numerous in vitro conditions. In addition, biophysical modeling is coupled to elucidate the effects of flow velocity and contact area, as well as provide validation from the biomechanical standpoint. Simulations from this combined modeling framework are further validated against erythrophagocytosis experimental data collected from our specially designed organ-on-a-chip microfluidic device. We demonstrate that sickle RBCs exhibit a faster and higher degree of engulfment by macrophages compared to normal RBCs due to their altered molecular expression on the cell membrane and increased cell rigidity. Moreover, the uptake process can be further accelerated under hypoxia. Overall, the novel hybrid modeling framework that we develop provides mechanistic insights into the signaling control of erythrophagocytosis and serves as a platform to quantitatively predict the effects of molecular perturbations and lend insight into therapeutic developments in SCD.

Results

A mechanistic computational model of the signaling interaction between erythrocytes and splenic macrophages predicts the activation of myosin under different conditions.

We first construct a computational systems biology model of the signaling interactions between RBCs and splenic macrophages with detailed molecular mechanisms. A diagram of the signaling model is shown in Fig. 1A. The binding of PS to PS-binding proteins, binding of band-3 to NAb-opsonized FcγRs, and rigidity serve as “eat me” signals and promote the activation of nonmuscle myosin IIa in splenic macrophages. The binding of CD47 to SIRPα receptors activates SHP1 and subsequently inhibits the activation of myosin, serving as a self-marking “do not eat me” signal. Activation of SHP1 can also inhibit the activation of SIRPα as negative feedback. In aged or oxidatively damaged RBCs, conformationally altered CD47, labeled CD47_alt, binds to TSP1 and SIRPα, leading to downstream signaling that promotes myosin activation. A detailed list of model reactions and a table of all parameters of the model, along with their descriptions, are available in Tables S1 and S2. The model is also available in the Supplementary Materials as a model file in systems biology markup language (SBML) (File S1). To assess the structural identifiability of the parameters, we use an input-output projection-based method to identify structurally unidentifiable parameters (30). Structurally unidentifiable parameters are fixed to their reference values obtained from the literature, and identifiable parameters are fitted to the experimental data (see Table S3). To estimate the parameter values, the model is calibrated to experimental measurements of SIRPα and myosin activation in the presence and absence of anti-CD47 antibody (31), as well as time-dependent activation of SIRPα during erythrophagocytosis (11). The calibrated model predicts that blocking the action of CD47 with an antibody results in a 3.7-fold decrease in SIRPα activation and a 4.2-fold increase in the activation of nonmuscle myosin IIa (Fig. 1B). The model also demonstrates the time-dependent inhibitory effect of anti-CD47 antibody on both SIRPα activation (Fig. 1Ci) and SHP1 (Fig. 1Cii), as well as the stimulatory effect on the activation of myosin (Fig. 1Ciii).

Integrated model of erythrophagocytosis is calibrated to experimental data and predicts the effects of opsonization, CD47 blockade, and rigidity on RBC engulfment.

To simulate the level of RBC engulfment on a macroscopic level, we generate virtual populations of RBCs and simulate their signaling interactions with macrophages. For each RBC–macrophage interaction, the activation of myosin over their simulated interaction time is calculated. The activation of myosin is then integrated over the contact time to calculate the area under the curve (AUC), or exposure, of activated myosin. In the simulations, the RBC is determined to be engulfed when the exposure of activated myosin reaches a threshold level. The diagram of the modeling framework consisting of the signaling model, population model, and biophysical model is shown in Fig. 2A. The integrated model is referred to as biophysics-informed model (BIM) of signaling and engulfment of populations of RBCs, which is then used to simulate the phagocytosis of RBCs at different levels of opsonization by antiserum. Here we describe the calibration of BIM using experimental measurements and biophysical simulation results of RBC phagocytosis under different opsonization, CD47 blockade, and rigidity conditions to estimate the model parameters.

To calibrate the model, parameters (see Table S1) are fitted to the following datasets: experimental measurement of the engulfment of human RBCs measured by Sosale et al. (10), which demonstrates a dose-dependent increase in the engulfment as the level of opsonization by antiserum increases from 0 to 10 mg/mL (Fig. 2Bi); experimentally measured engulfment of native, opsonized RBCs that is further enhanced by inhibiting their self-marker CD47 by Sosale et al. (10) (Fig. 2Bii); simulations of the phagocytosis of RBCs of varying rigidity, from native to 30-fold increase, predicted by the particle-based biophysical model of macrophages and RBCs (see Materials and Methods and SI Methods for details of the biophysical model development and Fig. S1 for its validation) (Fig. 2Biii). Representative frames of the dataset produced by the biophysical model showing the adhesion and partial engulfment of flexible RBCs and the complete engulfment of rigid RBCs are shown in Fig. 2C. Figure 2C(i) (see also Movie S1) shows that when the adhesive force between macrophages and RBC is low (4 pN/μm²), corresponding to RBCs with a low concentration of ligands expressed on their membrane, the target RBC only adheres to the macrophage without initiating the uptake process. On the other hand, at a higher level of adhesive force (50 pN/μm²), our results show that macrophage can engulf rigid RBCs (see Fig. 2Cii and Movie S2) but fail to do so for flexible RBCs (see Fig. 2Ciii and Movie S3). As shown in Fig. 2Biii, the model predicts that at a moderate level of opsonization, an increase in rigidity promotes the uptake of RBCs, in alignment with the simulation results from the biophysical model. To validate the model, we simulate the engulfment of rigid and opsonized RBCs at varying CD47 inhibition levels and compared the simulation results to an experimental dataset from Ref. (10) which is not used in model calibration. To validate the model, BIM simulations of rigid RBCs are compared to experimental measurements by Sosale et al. (10), both showing a high level of engulfment regardless of CD47 inhibition (Fig. 2Bii, black). To assess the uncertainty of the model, all simulations of erythrophagocytosis are repeated 200 times with randomly generated virtual populations of RBCs, and the 98% confidence intervals of the simulations are visualized in Fig. 2B.

Calibrated model of erythrophagocytosis can predict the change in phagocytosis of aged RBCs in the presence and absence of TSP1.

In aged RBCs, increased PS exposure and the conformational change of CD47 lead to increased engulfment that is further enhanced by TSP1 peptide stimulation (13, 32). To validate the model, we generate populations of aged RBCs and simulate their engulfment in the absence and presence of TSP1 peptide stimulation using BIM and compare our simulation results with experimental observations not used in model calibration. Based on literature evidence, we generate virtual populations of young and aged RBCs by assuming that PS exposure in aged RBCs increases by 50% and that 20% of CD47 molecules have undergone the conformational change that allows them to bind TSP1 and activate prophagocytic signaling (8, 13, 32). The expression level distributions of PS, CD47 without conformational change, and conformationally changed CD47 that could bind TSP1, noted as CD47_alt in the model, are presented in Fig. 3A. Our model predicts that young, healthy RBCs in control conditions have a low level of engulfment of ∼2%. Inhibiting CD47 with an antibody significantly increases the engulfment level to around 45%. Compared to young RBCs, in the absence of TSP1, aged RBCs are subject to a higher level of engulfment (Fig. 3B). These results agree well with the experimental measurements by Burger et al. (13). The model further predicts that without TSP1, the conformational change of CD47 alone causes a moderate increase in engulfment. Adding TSP1 enhances the pro-phagocytosis effect of the conformational change of CD47, reaching a 50% engulfment when 30% of CD47 is conformationally changed (Fig. 3C). Simulations of the engulfment of healthy and aged RBCs in the absence and presence of TSP1 predict that aged RBCs are subject to faster engulfment that is further accelerated by TSP1 (Fig. 3D). Simulations of erythrophagocytosis are repeated 200 times with randomly generated virtual populations of RBCs, and the 98% confidence intervals of the simulations are visualized in Fig. 3C and D.

BIM predicts that sickle RBCs exhibit stronger “eat me” signals that are further enhanced by hypoxia, resulting in a faster and higher degree of engulfment.

To further demonstrate BIM's ability to predict RBC engulfment dynamics with detailed mechanistic control and validate it, we generate virtual populations of sickle RBCs under hypoxia and normoxia to predict their engulfment over time and compare them to experimental measurements that are not used in model calibration.

Compared to healthy RBCs, sickle RBCs have increased PS exposure (19), making them more susceptible to engulfment by splenic macrophages. The expression of CD47 has been shown to slightly decrease in sickle RBCs (22), while others report the change to be insignificant (13, 23, 33). The conformational change of CD47 in sickle RBCs further contributes to their decreased antiphagocytic signaling (13). The expression of band-3 has also been observed to be slightly decreased (21). In a splenic sequestration crisis, obstruction of RBCs and anemia can induce local hypoxia in the spleen, which has been shown to significantly affect the phagocytosis dynamics of RBCs (17). To investigate the effect of the molecular changes induced by hypoxia on sickle RBCs and their subsequent ingestion by macrophages, we used flow cytometry to measure the expression of the signaling ligands PS, band-3, and CD47 on sickle RBCs under normoxic and hypoxic conditions from two blood samples of SCD patients. Our experimental results show that hypoxic sickle RBCs have an ∼25% increase in band-3 expression and an 84% increase in PS expression compared to normoxic sickle RBCs. However, we did not find any significant changes in the relative expression of CD47 in sickle RBCs under hypoxia compared to that in sickle RBCs under normoxia (Fig. 4A, see also Fig. S2A). The oxygenation level also did not affect the expression of these ligands in RBCs collected from healthy patients (Fig. S2B).

According to the deformability measurements of healthy and sickle RBCs under hypoxia and normoxia using electrodeformation as described in our previous studies (17, 34–36), the rigidity of sickle RBCs under normoxia is measured to be significantly increased by 1.8-fold compared to healthy RBCs, and sickle RBCs under hypoxia show a further increase of 6.2-fold in rigidity compared to healthy RBCs (Fig. 4C). Informed by literature data described above, we generated virtual populations of sickle RBCs by assuming decreased band 3, increased PS, 10% decreased total CD47 with 25% conformational change (Fig. 4B). Note that while Fig. 4B shows the simulated absolute expression levels, and Fig. 4A shows the relative expressions measured using flow cytometry. Figure 4Biii shows the amount of nonconformationally changed CD47 indicative of their antiphagocytic signaling. We then generate virtual populations of sickle RBCs under hypoxia based on our experimental measurements of their ligand expression and rigidity change under hypoxia. The expression levels of band-3, PS, and CD47 in the virtual populations of RBCs (healthy and sickle) are shown in Fig. 4B.

We then use our BIM model to simulate the engulfment dynamics of RBCs. Simulations predict that sickle RBCs undergo a much faster and higher degree of ingestion than healthy RBCs, with up to 70% of RBCs being phagocytosed at 30 min. Our model simulations further suggest that sickle RBCs under hypoxia are susceptible to an even higher level of engulfment, with over 90% of RBCs engulfed at 30 min. The engulfment of hypoxic sickle RBCs is faster than that of normoxic RBCs, with a predicted half-engulfment time of ∼13 min compared to 20 min observed for normoxic sickle RBCs (Fig. 4D). Our model predictions agree well with the experimental measurements of the engulfment dynamics of sickle RBCs using an oxygen-mediated spleen-on-a-chip platform, as described in our previous study (17). In the experiments, we also observed that sickle RBCs under hypoxia are generally subject to a higher degree of engulfment and faster ingestion, both at the individual and cell population level (Fig. 4D and E); Fig. 4E shows time-lapse images of representative engulfment processes of a single nonsickle RBC under normoxia (first row of Fig. 4E and Movie S4) and a sickle cell under hypoxia (second row of Fig. 4E and Movie S5).

Sensitivity analysis suggests that myosin activation can be enhanced by targeting the CD47 axis or by increasing opsonization.

Inadequate phagocytosis of retained RBCs can lead to congestion in the spleen, impairing spleen function, and can potentially cause life-threatening spleen complications (17). A mechanistic model allows us to quantitatively characterize the effects of perturbing each molecular mechanism on promoting phagocytosis signaling and accelerating it. To identify the key molecular mechanisms that significantly impact the signaling control of the engulfment of sickle RBCs, we perform a global sensitivity analysis on the signaling model by sampling the parameter space and calculating their partial rank correlation coefficients (PRCCs) with regard to the activation of myosin. The PRCCs of the most sensitive parameters are shown in Fig. 5A. Parameters with positive PRCCs are positively correlated with active myosin IIa. Sensitivity analysis shows that the top sensitive parameters that positively affect myosin activation include those related to the expression of myosin (Myo2a_inact_0), expression of molecules related to “eat me” signals (PSR_0, FCR_0, CD47alt_0), inactivation of myosin-inhibiting SHP1 (SHP_hillcoeff, kf_SHP1inact), activation of “eat me” signals (kf_myo2aact_PSR, kf_myo2aact_FCR), dissociation of CD47 and SIRPα (kr_CD47SIRPα), area of RBC (S_rbc), and contact area (S_contact). Promoting the increase in these parameters can be considered to boost myosin activation (Fig. 5Ai).

On the other hand, parameters with negative PRCCs are negatively correlated with myosin activation. Sensitivity analysis shows that parameters related to the inactivation of myosin by SHP1 (kf_myo2ainact, SHP_kd, kf_SHP1act), expression of SHP1 (SHP1_inact_0), expression of “do not eat me” signal-related CD47 and SIRPα (CD47_0, SIRPα_0), activation of SIRPα by CD47 (kf_CD47SIRPaact, kf_CD47SIRPa), degradation of active “eat me” signal receptors (kdeg_PSRact, kdeg_FCRact), and dissociation of conformationally changed CD47 with SIRPα (kr_CD47altSIRPa). Thus, inhibiting molecular processes related to these parameters can potentially promote myosin activation (Fig. 5Aii).

Based on the sensitivity analysis results, we select parameters that correspond to molecular mechanisms that could be potential therapeutic targets and simulate the effects of modulating these parameters on the engulfment of sickle RBCs using BIM. From the positively correlated parameters, increasing the expression of FCγR (FCR_0) or the dissociation rate of CD47 with SIRPα (kr_CD47SIRPRa) can increase the engulfment of sickle RBCs at 20 min from 5% to over 90%. Increasing the inactivation of SHP1 (kf_SHP1inact) also increases the engulfment from 10 to 75% (Fig. 5B and Di). From the parameters with negative PRCCs, simulations show that increasing the parameter related to myosin inactivation by SHP1 (kf_myo2ainact) can drastically decrease the engulfment of sickle RBCs (from 100 to 0%) and increasing the expression of SHP1 and SIRPα (SHP1_0, SIRPa_0) also results in a decrease in engulfment (Fig. 5C). On the other hand, decreasing these negative parameters can also significantly elevate RBC engulfment (Fig. 5Dii).

To identify key pair interactions of the parameters affecting myosin activation and RBC engulfment, we calculate the pairwise second-order Sobol sensitivity indices of the model parameters by sampling the parameter space (Fig. S3A). The parameter pairs with the highest second-order Sobol indices are shown in Fig. S3B. Our results indicate that the expression rates of SIRPα (SIRP1_0) and CD47 (CD47_0) have pair interactions with the expression of myosin IIa (myo2a_inact_0), suggesting that the effects of single therapies targeting SIRPα or CD47 alone can be prone to resistance through overexpression of myosin IIa. In addition, the inactivation rate of myosin (kf_myo2ainact), also previously identified by partial rank correlation coefficient (PRCC) as a sensitive parameter, exhibits pair interactions with many other parameters, including the expression levels of PSR (PSR_0), Fc Receptor (FCR) (FCR_0), SHP1 (SHP1_inact_0), myosin IIa (myo2a_inact_0), and others, suggesting that modulating the deactivation of myosin IIa by SHP1 alone can be compensated by many other molecular mechanisms including the expression levels of macrophage receptors and SHP1. The total-order Sobol analysis of all parameters identified a set of sensitive parameters consistent with PRCC, and all total-order Sobol indices are shown in Fig. S3C.

Altogether, our integrated model reveals that the activation of myosin IIa and RBC engulfment can be enhanced by strategies including targeting the CD47/SIRPα axis, increasing opsonization, and modulating the inactivation of myosin IIa by SHP1. The second-order Sobol analysis further suggests that single therapies targeting CD47, SIRPα, or inactivation of myosin IIa alone are prone to compensatory mechanisms that may hinder their efficacy.

Combining the biophysical model and the signaling model predicts engulfment patterns of erythrocytes under different flow conditions

Engulfment of RBCs by macrophages can also be impacted by the speed of fluid flow in the surrounding environment (37, 38), although the blood flow in the red pulp of the spleen is thought to be sluggish (100–500 μm/s) to facilitate the interactions between RBCs and macrophages. To demonstrate how our model can be adapted to investigate in silico how flow condition affects the interaction between RBCs and macrophages, here we demonstrate our integration of our signaling model with biophysical model. Using the same biophysical model as shown in Fig. 2C, we further incorporate the impact of the blood flow on the contact area between RBCs and macrophages. Representative frames of the simulations under normal and high-speed flows are shown in Fig. 6Ai and Bii, respectively, predicting a marked change in the contact area between the macrophages and RBCs with varying stiffness under different flow conditions (Fig. 6B). Movies of biophysical simulations under normal and high-speed flows are shown in Movies S6 and S7. The biophysical model simulations in Fig. 6A show that blood velocity (from 100 to 500 μm/s) does not affect the contact area between RBCs and macrophages in the case of low-rigidity RBCs. For rigid RBCs with a 10-fold increase in their stiffness, a flow velocity of 400 μm/s reduces the contact area from 55 to 15 μm², and a further increase in flow velocity will cause the RBCs to be washed away without contacting the macrophage (Fig. 6A). For RBCs with a 20-fold increase in their stiffness, a flow velocity of 300 μm/s can already reduce the contact area from 51 to 19 μm² (Fig. 6B).

Here, we establish the biophysics-coupled multiscale model (BCMSM) by incorporating the above flow-area relationships from biophysical modeling with our signaling and population model framework. After incorporating the biophysical model, our mechanistic simulations predict that rigid RBCs under slow flow undergo a faster and higher level of engulfment (∼40% at 40 min) compared to healthy controls and that rigid RBCs under fast flow undergo a lower degree of phagocytosis (∼10% at 40 min) compared to the case of slow flow (Fig. 6C). We also investigate the effect of contact area on the engulfment of RBCs with different rigidities. Model simulations predict that for healthy control RBCs, a decrease in contact area to 20 μm² does not affect the level of engulfment. For healthy rigid RBCs with 10- and 20-fold increase in their rigidity, decreasing the contact area until 20 μm² results in a decreased engulfment from 45% to ∼30%, whereas a decrease beyond 20 μm² results in a slight increase in engulfment (Fig. 6Di). Simulation of the engulfment of sickle RBCs demonstrates a similar trend under both hypoxia and normoxia (Fig. 6Dii, note that the phagocytosis of sickle RBCs is higher compared to normal RBCs and the y-axis is from 75 to 100%, different from Fig. 6Di). In biophysical models of RBC phagocytosis, macrophage engulfment requires sufficient contact area with RBCs, and reducing contact area between the two cells is always predicted to reduce phagocytosis. Experimental results, however, have shown that there can be a paradoxical increase in engulfment when the contact area becomes very small while the cell is still rigid and remains in contact with the macrophage (10). Our integrated signaling-biophysical model takes into account the loss of “self/antiphagocytic” signaling through CD47 and is able to predict this phenomenon.

Discussion and summary

In this study, we use a combined approach by integrating a mechanistic systems biology model and a biophysical model of RBC–macrophage interactions during erythrophagocytosis to quantitatively characterize the dynamic control of RBC engulfment by splenic macrophages. The ligand-receptor interactions, feedback, and regulatory mechanisms, signal redundancy, and interplay between the “eat me” and “do not eat me” signals form a complex reaction network that should be studied at the systems level. Systems biology models allow for quantitative characterization of the effects of individual biochemical reactions on a network level, and we have previously developed mechanistic computational models to investigate the complex signaling networks in endothelial cells and macrophages (39–42). Biophysical models, on the other hand, allow for probing the effects of biomechanical properties and have been used to simulate different processes in RBC sickling and retention (25, 26). Our multiphysics and multiscale framework that integrates the systems biology model, biophysical model, and experimental studies can significantly facilitate our mechanistic understanding of the biochemical and biomechanical processes during erythrophagocytosis. In this study, we demonstrate our integrated model calibrated using experimental and biophysical simulation data, can accurately predict the changes in the engulfment dynamics of aged and sickle RBCs in relevant pathophysiological conditions. Our validated model further predicts that targeting SHP1 and SIRPα/CD47 can facilitate RBC clearance during splenic sequestration crisis and provides mechanistic insights into the effect of changing the flow speed in the spleen on phagocytosis, motivating further investigations for the therapeutic development of this life-threatening condition. The model can also be adapted and extended to model other chemical, molecular, and physical processes of RBC clearance such as the digestion of RBCs following internalization to better understand the pathophysiology and simulate potential therapeutic strategies to alleviate spleen sequestration.

Current clinical management of ASSC commonly involves prophylactic blood transfusion or splenectomy in severe cases (28). Here, we use the integrated mechanistic modeling framework to identify targets for medical treatments that can serve as alternatives to more invasive options. Using global sensitivity analysis, we identified the TSP1/CD47/SIRPα axis as a key modulator of erythrophagocytosis of sickle RBCs. The effect of TSP1 on phagocytosis is well characterized in aged RBCs (13), but its exact role in the phagocytosis of sickle RBCs is unclear, although the conformational change of CD47 that allows for TSP1 binding has been observed in sickle RBCs (33). TSP1 is also observed to induce the adhesion of sickle RBCs to macrophages (43) and its plasma expression is elevated in SCD (44). The role of CD47 is more apparent in addition to the various in vitro findings discussed in this work; in vivo in a murine model of SCD, silencing of CD47 has been observed to correlate with significantly less vascular congestion in the lung (44), suggesting the value of CD47 blockade to ameliorate SCD symptoms which is in alignment with our model predictions. On the other hand, therapeutically targeting CD47 to promote immune activation has been a potential strategy for immunotherapy in many cancer types with numerous drug candidates in clinical trials (45). However, given the high expression of CD47 on RBCs, anemia due to undesired RBC killing is a common and serious complication of anti-CD47 therapies (46). Antibodies targeting SIRPα, such as BI-765063 (47) and ADU-1805 (48), are also under development as anticancer therapies with supposedly no anemia-related concerns (46). In this sense, our model also provides a unique computational framework to facilitate the investigation of possible strategies to circumvent and mitigate treatment-related adverse events (e.g. anemia) during anti-CD47/SIRPα drug development. Inhibiting SHP1 is another potential strategy to accelerate phagocytosis. The loss of SHP1 has been observed to enhance macrophage activity and is linked to antitumor activity (49), and SHP-1 inhibitors have been therapeutically explored in the context of cancer (50). The second-order Sobol sensitivity analysis further allows us to identify key pair interactions and preemptively predict the potential resistance mechanisms and strategies for combination therapies including simultaneously targeting the CD47/SIRPα axis and SHP1-inactivation of myosin to overcome resistance.

We note that our model aims mainly to capture the effects of molecular, biophysical, and environmental changes on phagocytosis, and some ligand-receptor interactions are simplified representations of their actions. In addition to binding to FcγR after opsonization by NAbs, band-3 is also known to activate prophagocytic signaling by activating the complement system (15). Macrophages are also known to express other PS-recognizing receptors including TIMs, integrins, and Stabilins, to recognize the exposed PS on erythrocytes and promote the “eat me” signal (7). In addition, both PS exposure and the rigidity of RBCs are implicated in changes in intracellular calcium levels (51). Therefore, we model the effect of increased rigidity by increased activation of PSRs and subsequent myosin activation. We also note that the total number of targets internalized by the macrophage may be limited by its surface membrane and is difficult to count due to the vesicle fusion process. Our model assumes that the simulated interaction between RBCs and macrophages is the internalization during the initial engulfment of the first RBC by the macrophage, consistent with the experimental measurements of RBC internalization dynamics, where we only quantified the internalization time during the initial engulfment of the first RBC as the phagocytic index. These simplifications of the mechanistic details of our model are made carefully in consideration of the limited data and model complexity without affecting the predictive power of the model. We also note that our integrated model focuses on the engulfment and uptake of RBCs by splenic macrophages, which is a key step in the initiation of erythrophagocytosis. The digestion of RBCs by macrophages following internalization is not considered in the current model as the detailed mechanisms that regulate digestion are still not clear. New qualitative experiments are required to build a kinetic model to describe the digestion process (17).

In summary, we present a framework for integrating systems biology modeling, biophysical modeling, and experimental studies to elucidate the mechanisms underlying erythrophagocytosis in the context of SCD. Our established model captures and reproduces the effects of opsonization, CD47 inhibition, and rigidity on the phagocytosis levels of RBCs consistent with experimental data and biophysical simulations and predicts the increased phagocytosis of aged RBCs that is further increased by TSP1 stimulation. Using experimental measurements of the ligand expression and rigidity changes of sickle RBCs under normoxia and hypoxia as inputs, the model predicts a faster and higher degree of erythrophagocytosis of sickle RBCs that is further accelerated by hypoxia. For model validation, simulations of the phagocytosis dynamics of sickle and aged RBCs show good consistency with experimental measurements. Using global sensitivity analysis, we identify the CD47/SIRPα axis and macrophage receptor expression as potential targets for therapeutic interventions to modulate erythrophagocytosis. Informed by the biophysical simulations, the model also predicts the effect of flow velocity on the phagocytosis of sickle RBCs. Overall, our model serves as a platform to quantitatively characterize the signaling control of erythrophagocytosis and to simulate the effects of molecular perturbations to inform experimental design and drug development.

Materials and methods

Systems biology modeling of the signaling interactions between RBCs and macrophages during erythrophagocytosis

We use a mass action and ordinary differential equation (ODE)-based computational model of the signaling interactions between RBCs and macrophages during erythrophagocytosis to simulate the activation of phagocytosis signaling in macrophages by RBCs. The model is constructed using the SimBiology toolbox and MATLAB, release R2022a (MathWorks Inc., Natick, MA, USA). A diagram of the signaling model is shown in Fig. 1A. Ligands PS, band-3, and CD47 are expressed at the surface of RBCs. Receptors FcγR, SIRPα, and PSR are expressed at the surface of macrophages. During erythrophagocytosis, PS binds and leads to the activation of PSR. The activation of PSR by PS is also enhanced by the rigidity of RBCs. Band-3 on the surface of RBCs, once opsonized by NAbs, binds, and activates FcγR on macrophages. Activated PSR and FcγR stimulate the downstream activation of nonmuscle myosin IIa, acting as major “eat me” signals to promote the phagocytosis of RBCs. CD47 on RBCs binds to and activates SIRPα receptors on the macrophages, leading to the downstream activation of SHP1. SHP1 inhibits the activation of nonmuscle myosin IIa, acting as a “do not eat me” signal. Conformationally changed CD47 can bind to TSP1 and activate an alternate downstream signaling pathway of SIRPα, causing it to act as an “eat me” signal to activate myosin instead of the “do not eat me” signal associated with CD47. Initial estimations of the ligand and receptor expression are derived from literature data (52–56). The expression levels of SHP1 and nonmuscle myosin IIa in splenic macrophages are also estimated based on literature data (57). The initial estimations of the binding, phosphorylation, and dephosphorylation rates are based on reported values in the literature (58–63). The surface areas of RBCs and macrophages are also derived from literature data (64, 65). A detailed list of the reactions of the model along with a table of all parameters in the model are included in Table S2. The model file used to generate the model network is also available in the Supplementary Materials in the SBML. All subsequent model analyses are performed using MATLAB R2022a (MathWorks Inc.). The model is calibrated to experimentally measured activation of SIRPα and myosin IIa in the presence and absence of anti-CD47 antibody from two independent sources (11, 31), and the patternsearch algorithm from the global optimization toolbox in MATLAB (MathWorks, Inc.) is used for parameterization.

Computational modeling of phagocytosis of virtual populations of RBCs

To simulate the level of phagocytosis on a macroscopic level, we generate virtual populations of RBCs and simulate their signaling interactions with macrophages. Virtual populations of RBCs are generated assuming that the expression of PS, band-3, and CD47 follows log-normal distributions (19, 53, 54, 66). The interactions between the RBCs and macrophages are simulated for 45 minutes, consistent with the incubation time in the experiments by Sosale et al. (10). During their interaction time, the prolonged activation of myosin in macrophages by the “eat me” signals leads to phagocytosis (5). To quantitatively simulate the accumulation of the myosin activation over their interaction time, we integrate the activated nonmuscle myosin IIa over the simulation timespan to calculate the exposure, or AUC of activated myosin. Our simulations show the increase of myosin exposure over time during their contact, and the RBC is determined to be engulfed when myosin exposure reaches a threshold level. The threshold for RBC engulfment is calibrated to experimental data as described below. A diagram of the phagocytosis model informed by the signaling interactions between RBCs and macrophages is shown in Fig. 2A. The model implementations and subsequent analyses are performed in MATLAB release 2022a (MathWorks Inc.). The model is used to simulate the phagocytosis of RBCs at different levels of opsonization, and the phagocytosis criterion is calibrated to the experimental data reported in Sosale et al. (10), where at an NAb concentration of 10 mg/mL, ∼33% of RBCs are phagocytosed. The model is also calibrated to experimental data measuring the phagocytosis of native RBCs under varying levels of CD47 inhibition and simulations of RBC phagocytosis at different rigidity by the biophysical model (10). For validation, the simulation of the phagocytosis of rigid RBCs with increasing degrees of CD47 inhibition using the calibrated model is compared to experimental measurement (10). For further validation, the calibrated model is then used to simulate the engulfment of aged and sickle RBCs under different conditions and compare them to experimental measurements, as described in Figs. 3 and 4.

Identifying key parameters using global sensitivity analysis

To identify the key parameters that impact the activation of myosin IIa and phagocytosis level, we perform a global sensitivity analysis on the signaling model by using Latin hypercube sampling (LHS) to sample the parameter space and calculate the PRCCs of each parameter with regard to the signaling outcome activated myosin. This method of characterizing the parameter sensitivity is efficient and well established for systems biology models and is described in detail in Marino et al. (67). The parameters and initial values are sampled with LHS and constrained in the range of 1/10th to 10-fold of the baseline values and their PRCCs are calculated. Parameters with positive PRCCs, indicating that increasing these parameters results in increased outcomes, are shown in Fig. 5Ai. The top parameters with negative PRCCs are shown in Fig. 5Aii. We further select the parameters that are most influential on the activation of myosin to simulate the effect of modulating them using the multiscale BIM to quantitatively characterize how modulating them affects RBC engulfment dynamics (Fig. 5B–D). Second-order and total sobol sensitivity analyses are also performed by sampling the parameter space in the range of 1/10th to 10-fold of the baseline values using the UQLab framework for uncertainty quantification (68). The second-order sobol indices for all parameter pairs, top parameter pairs ranked by their sobol indices, and total-order sobol indices for all model parameters are shown in Fig. S3.

Measurements of oxidized band-3, PS, and CD47 expression on sickle RBCs using flow cytometry

We obtained healthy blood samples from a local blood bank and collected sickle blood samples from homozygous SCD patients in a local clinic under an excess human material protocol approved by the Partners Healthcare Institutional Review Board (IRB) with a waiver of consent. We collected additional sickle blood samples from patients at the University of Pittsburgh under University of Pittsburgh IRB protocol PRO08110422. All blood samples were de-identified before experiments. The samples were stored in an anticoagulant tube containing ethylenediaminetetraacetic acid before the experiment. Red blood cells (RBCs) were washed twice with phosphate buffered saline (PBS) solution (Sigma-Aldrich, St. Louis, MO, USA) at 2,000 rpm for 2 min at room temperature (20 °C). RBCs were opsonized with 0.5 µM IgG (Rockland, Limerick, PA, USA) in PBS solution at room temperature (20 °C) for 30 min. THP-1 cells (American Type Culture Collection [ATCC], Manassas, VA, USA) were cultured in microfluidic chambers containing Roswell Park Memorial Institute (RPMI) media (StemCell Technologies, Cambridge, MA, USA) supplemented with 10% fetal bonvine serum (FBS) (Sigma-Aldrich) and then differentiated for 2 days using 100 ng/mL phorbol myristate acetate (PMA) (StemCell Technologies). Microscopic imaging was performed using Olympus IX71 inverted microscopes (Olympus, Tokyo, Japan) and an Olympus DP72 camera. All the phagocytosis experiments were performed at 37 °C using an incubation chamber (ibidi heating system; ibidi USA, Fitchburg, WI, USA). More details can be found in our recent work (17).

Isolated sickle RBCs were washed and diluted with fluorescence-activated cell sorting (FACS) buffer and 1 × 10⁶ cells/mL were used per assay. Cells were stained with 3 µL of the Pacific Blue antihuman CD47 antibody (BioLegend). Flow cytometry analysis was performed using a BD FACSAria III flow cytometer and FACS DIVA software.

The treatment of sickle RBCs by the normoxic and hypoxic conditions was obtained by incubating the RBC suspension in PBS with 1% (w/v) bovine serum albumin (BSA) (EMD Millipore, Billerica, MA, USA) by injecting gas of different mixtures at different oxygen levels for 30 min in a sealed chamber, respectively: 20% oxygen (O₂), 5% carbon dioxide (CO₂) with the balance of nitrogen (N₂), and an oxygen-poor gas mixture: 2% O₂, 5% CO₂ with the balance of N₂. After the normoxic/hypoxic treatments, RBCs were fixed using 4% formaldehyde. RBCs were then stained with the desired antibodies for 30 min on ice and washed twice with PBS with 1% (w/v) BSA for 5 min at 2,000 rpm. The stained RBCs were resuspended in PBS at a desired concentration of 10⁵∼10⁶ cells/mL and analyzed using a FACS Canto II HTS-1 (BD Biosciences) flow cytometer. Data were processed using FlowJo software (BD Bioscience). Antibodies used for staining were allophycocyanin (APC)-conjugated anti-CD233 (1:400, American Research Products), Brilliant Violet 421-conjugated anti-CD47 (1:400, BioLegend), and pSIVA-IANBD (1:100, Novus Biologicals).

Mechanical testing of sickle RBCs and macrophages

The mechanical properties of sickle RBCs were characterized by the shear modulus values of the cell membranes using the electrodeformation method under controllable oxygen tension. The mechanical properties of the macrophages were evaluated by Young's modulus using micropipette aspiration, where a portion of the cell was sucked into a micropipette with a tip diameter of 2 µm by applying negative pressures from a water column with adjustable height. See more detailed information in Ref. (17).

Biophysical modeling and phagocytosis under quiescent and flow conditions

We employ the model introduced in the work of Li et al. (69) for simulating macrophages and the model proposed by Fedosov et al. (70) to represent normal and sickle RBCs. The membrane of macrophages and RBCs is constructed by a 2D triangulated network with N_v vertices. The vertices are connected by N_s elastic bonds to impose proper membrane elasticity. The energies of the macrophage and RBC models (V_cell) are given by V_cell = V_s + V_b + V_{a  +  v}, where the elastic energy V_s represents the elastic interactions of the cell membrane and the bending energy V_b denotes the bending resistance of the cell membrane (71–75). The area and volume constraints V_{a  +  v} are imposed to mimic the area-preserving lipid bilayer and the incompressible interior fluid for macrophages and RBCs. Upon macrophage activation, we apply an additional active force F_a to the activated membrane (highlighted as yellow in Fig. 2C) to mimic the protrusive force induced by actin polymerization. This protrusive force is applied along the normal outward direction of the activated membrane to enable membrane extension toward the targets. Once the protrusions from the macrophages make contact with the surface of the target, they become adhesive (red) to it, representing the collective effects of the engagement of receptors on macrophages to their ligands on the targets. The detailed formulation of the potentials implemented in the blood cell models is introduced in Li et al. and Fedosov et al. (70).

The adhesion between macrophages and their phagocytic targets is controlled by various ligand–receptor interactions on the cell membrane. Since the length scales of these adhesive molecules (∼10 nm) are commonly much smaller than the size of the cells (∼10 µm), it is feasible to coarse-grain their effects into an effective interaction potential that closely mimics the microscopic forces and reproduces realistic dynamical behavior rather than explicitly modeling all constituents, which could be computationally prohibitive. Along this line, we apply a Morse potential between the particles on macrophages and the targets when they are in the vicinity, expressed as

where r denotes the distance between two particles. β denotes the interaction range and r₀ represents the zero-force distance. D_e is the depth of the potential well and it controls the strength of the adhesion between macrophages and their targets. In this study, D_e is tuned such that the adhesive stress (adhesive force per unit area) in the model, which represents the extent of ligand-receptor binding, varies between 0 and 375 pN/µm², a range that is consistent with the adhesive stress reported from phagocytosis experiments (76) and analytical estimation (77).

For the phagocytosis simulation under flow, we place the macrophages on the bottom of a microchannel with a length of l = 50 µm, a width of w = 25 µm, and a height of h = 20 µm, as illustrated in Fig. 6B. We assume firm adhesion between the macrophages and the channel wall. We approximate the flow field in the simulation domain as a Couette flow. The target RBCs flow into the channel from the inlet. The flow velocities are selected to be consistent with the blood flow conditions in the red pulp of the spleen. A detailed description of the biophysical model development and validation is included in SI Methods.

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

Part of this work was carried out at the Advanced Research Computing at Hopkins (ARCH) core facility (rockfish.jhu.edu), which is supported by the National Science Foundation (NSF) grant number OAC1920103. This work is supported by National Institutes of Health Grants R01HL154150 and R01HL101200 and Nanjing Medical University Career Development Fund NMUR20210006.

Author Contributions

Y.Z., H.L, C.Z., M.D., G.E.K., and A.S.P. conceptualized and designed the study. Y.Z. and C.Z. constructed the biochemical signaling model, and phagocytosis model and performed all model analyses. Y.Q. and M.D. designed and performed flow cytometry and phagocytosis experiments. H.L., G.L., and G.E.K. designed and performed biophysical simulations of erythrophagocytosis. L.L., M.D., G.E.K., and A.S.P. provided critical input and jointly supervised the study. Y.Z., C.Z., Y.Q., and H.L. drafted the manuscript. All authors reviewed and edited the manuscript.

Data Availability

The computational model used to generate all simulations in this study is available in the Supplementary file as a .sbml file in systems biology markup language.

References

Author notes