# In silico biophysics and rheology of blood and red blood cells in Gaucher Disease

Zhaojie Chai, Guansheng Li, Papa Alioune Ndour, Philippe Connes, Pierre A Buffet, Mélanie Franco, George Em Karniadakis

► To cite this version:

Zhaojie Chai, Guansheng Li, Papa Alioune Ndour, Philippe Connes, Pierre A Buffet, et al. In silico biophysics and rheology of blood and red blood cells in Gaucher Disease. PLoS Computational Biology, 2025, 21 (9), pp.e1012705. 10.1371/journal.pcbi.1012705. hal-05357245

HAL Id: hal-05357245

https://hal.science/hal-05357245v1

Submitted on 10 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# In silico biophysics and rheology of blood and red blood cells in Gaucher Disease

Zhaojie Chai<sup>1*</sup>, Guansheng Li<sup>1</sup>, Papa Alioune Ndour<sup>2</sup>, Philippe Connes<sup>3</sup>, Pierre A. Buffet<sup>2*</sup>, Melanie Franco<sup>2*</sup>, George Em Karniadakis<sup>1*</sup>

1 Division of Applied Mathematics, Brown University, Providence, Rhode Island, United States of America, 2 Université Paris-Cité, INSERM, EFS, BIGR U1134, Paris, France, 3 University Lyon, LIBM EA7424, Vascular Biology and Red Blood Cell Team, Universite Lyon 1, Villeurbanne, France

* zhaojie_chai@brown.edu (ZC); pabuffet@gmail.com (PAB); melanie.franco@inserm.fr (MF); george_karniadakis@brown.edu (GEK)

<figure><img src="image_2.png" /></figure>

## Abstract

Gaucher Disease (GD) is a rare genetic disorder characterized by a deficiency in the enzyme glucocerebrosidase, leading to the accumulation of glucosylceramide in various cells, including red blood cells (RBCs). This accumulation results in altered biomechanical properties and rheological behavior of RBCs, which may play an important role in blood rheology and the development of bone infarcts, avascular necrosis (AVN) and other bone diseases associated with GD. In this study, dissipative particle dynamics (DPD) simulations are employed to investigate the biomechanics and rheology of blood and RBCs in GD under various flow conditions. The model incorporates the unique characteristics of GD RBCs, such as decreased deformability and increased aggregation properties, and aims to capture the resulting changes in RBC biophysics and blood viscosity. This study is the first to explore the Young's modulus and aggregation parameters of GD RBCs by validating simulations with confocal imaging and experimental RBC disaggregation thresholds. Through *in silico* simulations, we examine the impact of hematocrit, RBC disaggregation threshold, and cell stiffness on blood viscosity in GD. The results reveal three distinct domains of GD blood viscosity based on shear rate: the aggregation domain, where the RBC disaggregation threshold predominantly influences blood viscosity; the transition area, where both RBC aggregation and stiffness impact on blood viscosity; and the stiffness domain, where the stiffness of RBCs emerges as the primary determinant of blood viscosity. By analyzing RBC mechanical properties and blood viscosity in relation to bone disease, we find that the RBC aggregation properties, deformability, and blood viscosity, may contribute to its onset. These findings enhance our understanding of how changes in RBC properties impact on blood viscosity and may affect bone health, offering a partial explanation for the bone complications observed in GD patients.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Chai Z, Li G, Ndour PA, Connes P, Buffet PA, Franco M, et al. (2025) In silico biophysics and rheology of blood and red blood cells in Gaucher Disease. PLoS Comput Biol 21(9): e1012705. https://doi.org/10.1371/journal.pcbi.1012705

**Editor:** Feilim Mac Gabhann, Johns Hopkins University, UNITED STATES OF AMERICA

**Received:** December 9, 2024

**Accepted:** August 21, 2025

**Published:** September 10, 2025

**Copyright:** © 2025 Chai et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data availability statement:** All data are in the manuscript and/or Supporting information files.

**Funding:** This study was supported by the National Institutes of Health (Grant No. R01HL154150 (G. E. K.)) and the Laboratory of Excellence GR-Ex (Grant ANR-11-LABX-0051). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.
## Author summary

In Gaucher Disease (GD), a genetic deficiency in the enzyme glucocerebrosidase leads to the accumulation of glucosylceramide in red blood cells (RBCs), resulting in altered biomechanical properties. These changes affect blood flow characteristics, particularly blood viscosity, and may contribute to bone health issues seen in GD patients, including bone infarcts, avascular necrosis (AVN), and other bone diseases. In our study, we apply dissipative particle dynamics (DPD) simulations to explore how GD impacts RBC behavior under various flow conditions. We model GD RBCs with decreased deformability and increased aggregation, examining how these properties influence blood viscosity across three distinct shear rate domains: aggregation, transition, and stiffness. By validating our simulations with confocal imaging data and experimental RBC disaggregation thresholds, we investigate the effects of RBC stiffness, aggregation, and hematocrit levels on blood flow in GD. We find that the RBC aggregation properties, deformability and blood viscosity, may contribute to the onset of bone disease. These findings improve our understanding of how changes in RBC properties influence blood viscosity and may contribute to bone health issues, providing a partial explanation for the bone complications observed in GD patients.

## Introduction

Gaucher Disease (GD) is a rare genetic disorder resulting from a deficiency in glucocerebrosidase (GCase) and stands as the most common autosomal recessive lysosomal storage disorder. The enzyme deficiency is caused by a biallelic mutation in the GBA1 gene [1]. Three types of GD are known. Type 1, which is the most common (90% of cases in Europe and the USA, but less in other regions), is characterized by hepato-splenomegaly, anemia, thrombocytopenia, complex bone disorders such as osteonecrosis with bone micro-infarcts in ischemic regions, and avascular necrosis (AVN) [2,3]. The extremely rare Type 2 (acute neuronopathic) and Type 3 (sub-acute neuronopathic) GD, in addition to the symptoms present in Type 1, are associated with severe neurological disorders in Type 2 and variable neurological symptoms in Type 3 [4]. GCase deficiency leads to the accumulation of its substrates, glucosylceramide and secondary sphingolipids, in macrophages. These cells, as natural phagocytes, are responsible for breaking down membrane glycolipids from red blood cells (RBCs) and leukocytes. The Lipid-laden causes macrophages to transform into Gaucher cells (GC) [5]. These GCs have a distinctive phenotype, characterized by a striated “wrinkled tissue paper” appearance [6]. In GD, these chronically, alternatively activated macrophages infiltrate various organs, including the liver, spleen, bone marrow, and sometimes the lungs [1]. The alterations of the monocytic lineage and the rise of GC are considered primary contributors to the major signs of GD [7]. However, several manifestations, such as hepatosplenomegaly, anemia, and bone disorders, may result, at least in part, from the presence of altered RBCs. RBCs from GD patients often exhibit abnormal morphologies, sphingolipid accumulation, reduced deformability, and increased aggregability [8–16]. RBC progenitors display differentiation defects leading to anemia [15].

In GD, RBCs show reduced deformability and higher aggregation properties [8]. These abnormal properties are attributed to the accumulation of glucosylceramide in RBCs [17]. However, investigations into biomechanical properties of RBCs from GD patients and their consequences on blood circulation and bone pathology remain limited [8, 14, 16, 18, 19].
An increase in hematocrit and in RBC membrane's shear modulus is known to cause a significant increase in blood viscosity [19]. Four main biophysical factors therefore may affect blood viscosity in GD: 1) hematocrit level, 2) RBC deformability, 3) plasma viscosity, and 4) sphingolipid concentration, which may result in attractive interactions between RBCs [8,20–22]. Hematocrit is higher in asplenic Gaucher patients than in healthy controls and is therefore a potentially important determinant of blood viscosity in GD [20]. By contrast, plasma fibrinogen concentration is similar in GD patients and healthy controls, indicating that increased RBC aggregation likely results from alterations in the RBC membrane than from plasma factors [8].

Because experimental *in vitro* measurements rely on rare samples with limited power to validate correlations, we decided to use computational methods to limit experimental variations during sample preparation [23]. Computational methods for studying blood dynamics and rheology include continuum-level models, which describe blood flow using constitutive laws, and discrete particle methods, which reconstruct different blood components. Continuum-level models often use either Newtonian or non-Newtonian representations of blood viscosity but rarely provide systematic sensitivity analyses of biophysical parameters [24]. In contrast, discrete particle methods offer a more accurate model of RBC membranes as coarse-grained particulate material, allowing for precise modeling of cell-cell interactions and accurate recovery of whole blood viscosity [25,26].

We then investigated GD RBC properties using dissipative particle dynamics (DPD) simulations. DPD accurately modeled blood viscosity by representing RBC membranes as coarse-grained particles [26–31]. This is the first time the DPD method was introduced to study the biophysics and rheology of blood and RBCs in GD. We validated the RBC disaggregation threshold and stiffness with experimental data. The computational model incorporated hemorheological parameters to simulate the behavior of RBCs under different flow conditions. By quantitatively assessing RBC deformability, aggregation, and flow characteristics, we examined the relationship between these parameters and GD pathophysiology, particularly its association with bone infarcts and other bone diseases.

# Materials and methods

## Ethics statement

We used data from different experiments previously done in Franco et al. [8], approved by the local ethics committee in accordance with the principles of the Declaration of Helsinki (CPP [Personal Protection Committee] agreement – Ile de France VI of the Hospital Group, obtained on 16.12.2010). These experimental data form the basis of this study: adhesion assay on laminin at 1 dyne/cm²; RBC aggregation analysis by syllectometry; and blood viscosity measured by cone/plate viscometer.

## GD RBC model

To investigate the role of GD RBC biophysical factors on blood viscosity, we needed to model the individual RBCs accurately, without compromising the essential hydrodynamics of the system and the convective transport processes that govern blood flow. RBC is a highly deformable biconcave nucleus-free membrane cell with a diameter of ~8 µm. Known as a viscoelastic object, RBC shows both liquid-like (viscous) and solid-like (elastic) responses to an applied deformation at the same time. This allows RBCs to be very deformable to squeeze through capillaries as small as 3 µm, preserve their shapes at small deformation rates, and orient in the flow direction at larger deformation rates in wider arteries [32–34]. The membrane
of RBC is modeled as a set of $N_v$ DPD particles with the three-dimensional coordinates $X_i$ ($i \in 1, ..., N_v$) in a triangulated network of springs with a dashpot on the surface. To model the incompressibility of the RBC membrane, area and volume constraints are applied. Also, considering the bending resistance between all neighboring triangles, the bending rigidity of the membrane can be mimicked.

In this work, we considered the aggregation of normal and GD RBCs under different disaggregation thresholds and stiffness, with the shear modulus and disaggregation threshold of RBCs under these conditions summarized in Table 1. The shear modulus and bending modulus of normal RBCs in our work were selected to be $Es_0 = 4.73\mu N/m$, $Eb_0 = 2.4 \times 10^{-19} J$ [35]. In this study, RBC stiffness is modulated by adjusting the shear elastic modulus ($E_s$). To model GD-RBCs with increased stiffness, we varied $E_s$ relative to control RBCs, while keeping the bending modulus ($Eb$) constant. We further validated our comprehensive RBC model by comparing it with experimental data, ensuring accurate simulations that capture the changes in biomechanical, rheological, and dynamic behaviors of GD RBCs in response to morphological alterations and membrane stiffening. For detailed formulas and parameters to maintain RBC morphology and mechanical properties, please refer to the Supporting Information (as shown in S1 Text).

## Cell-Cell aggregation models

In DPD method, particles interact through a series of pairwise interactions that represent the chemical and physical forces acting between different components. Details of the DPD method and the parameters adopted in this work are shown in the supplementary material.

RBC aggregation is primarily driven by two mechanisms: depletion and bridging [36]. In the depletion mechanism, molecules in the plasma create an osmotic pressure that induces an attractive force between adjacent cells, while in the bridging mechanism, specific macromolecules physically crosslink two cell surfaces [36,37]. In computational models, depletion-induced aggregation is typically represented using a Morse potential [38–40], or alternatively, a Lennard-Jones potential [41,42], while bridging-mediated aggregation is often modeled through reaction kinetics [37]. In this study, RBC aggregation is modeled using a Morse potential to approximate depletion-mediated interactions. The Morse potential captures both the attractive forces at intermediate distances and the short-range repulsive forces due to electrostatic effects between cell membranes. This approach aligns with previous modeling studies on RBC aggregation [38–40]. The Morse potential is defined as

$$V(r) = D_e \left( e^{-2\beta(r-r_0)} - 2e^{-\beta(r-r_0)} \right), \quad (1)$$

where $r$ denotes the distance between two particles, $D_e$ is the depth of the potential well, $\beta$ controls the interaction range, and $r_0$ represents the zero-force distance.

The Morse potential is applied to the RBC membrane to model cell-cell adhesion, with blue vertices (as shown in Fig 1Bi and 1Bii) representing three severity levels—severe, moderate, and mild—as described in previous work [33]. We emphasize that the blue vertices in the figure represent levels of adhesion strength, rather than specific spatial locations of adhesion on the RBC membrane. To construct RBC models that reflect varying levels of cell-cell adhesion strength, we adopt the hypothesis that sphingolipid molecules—accumulating in both the plasma and the RBC membrane in Gaucher Disease—can promote intercellular adhesion and enhance aggregation behavior [8,33,43].
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Overview of simulation setup and Cell-Cell interaction models. A) The front view of a typical simulation of blood viscosity at $Ht = 45\%$ (corresponding to 698 RBCs) and the size of the simulation RBC populations ($60\mu m \times 50\mu m \times 50\mu m$) employed in this study. B) Vertex-vertex interactions for the aggregation model. Bii) Representation of three different disaggregation thresholds through the Morse potential between constituents of different RBCs. Here, the presence of more blue vertices corresponds to a stronger RBC disaggregation threshold, which is commonly observed in blood with higher glucocerebroside accumulation levels. The severe condition corresponds to the highest RBC disaggregation threshold of $680 s^{-1}$ in the experiment; the moderate condition corresponds to a threshold of $250 s^{-1}$; and the mild condition corresponds to a threshold of $150 s^{-1}$. Detailed model parameters are provided in the Supporting Information and Methods section.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g001

In Fig 1 Bii, the presence of more blue vertices corresponds to a stronger RBC disaggregation threshold, which is commonly observed in blood with higher glucocerebroside accumulation levels. The severe condition corresponds to the highest GD RBC disaggregation threshold of $680 s^{-1}$ in the experiment; the moderate condition corresponds to a mean threshold of $250 s^{-1}$; and the mild condition corresponds to a lowest threshold of $150 s^{-1}$.

The adhesion strength is mainly correlated with the parameter $D_0$ in Eq 1. Detailed model parameters are provided in the Supporting Information (as shown in S2 Text) and Methods section.

## Computational rheometry

In a DPD simulation, each particle represents the center of mass for a small group of atoms or molecules, with its position and momentum updated in a continuous phase but at discrete time intervals. Particles i and j, located at positions $r_i$ and $r_j$ respectively, occur through three types of pairwise forces: conservative, dissipative, and random forces. The viscosity is determined by dividing the calculated shear stress by the shear rate. Having individual particle velocities and each of the pairwise interactions between particles, the stress tensor, S, can be directly calculated through the Irving-Kirkwood formalism [44] as follows:

$$S = -\frac{1}{V} \left\{ \sum_{i=1}^{N} m_i (\mathbf{v}_i - \mathbf{u}(\mathbf{r}_i)) \otimes (\mathbf{v}_i - \mathbf{u}(\mathbf{r}_i)) + \sum_{j>i} \sum_{i=1}^{N-1} \mathbf{r}_{ij} \otimes \mathbf{F}_{ij} \right\}. \quad (2)$$

Here, V is the volume of the entire fluid and $u(\mathbf{r}_i)$ is the streaming velocity imposed by the shearing flow. The $\otimes$ is the dyadic product of the two vectors. The shear viscosity of the DPD fluid is subsequently defined as the shear component of the stress tensor divided by the
<table><caption>Table 1. Model Parameters Extracted from Available Measured Data for RBCs in Normal and GD Conditions.</caption><thead><tr><th></th><th>CTR-RBC</th><th>GD-RBC1</th><th>GD-RBC2</th><th>GD-RBC3</th></tr></thead><tbody><tr><td>Shear modulus (μN/m)</td><td>4.73</td><td>42.57</td><td>94.60</td><td>18.92</td></tr><tr><td>RBC disaggregation threshold (s<sup>-1</sup>)</td><td>80</td><td>250</td><td>680</td><td>150</td></tr></tbody></table>

Notes: The mechanical properties of CTR-RBC are within the range of previous measurements [45–48]. The membrane properties for GD-RBCs are based on the experimental measurements by Franco et al. [8]. CTR-RBC represents RBCs from healthy controls, while GD-RBC1, GD-RBC2, and GD-RBC3 represent RBCs from GD patients with moderate, severe, and mild conditions, respectively.

https://doi.org/10.1371/journal.pcbi.1012705.t001

imposed deformation rate by the wall velocities

$$ (\eta = \langle S_{ij} \rangle / \dot{\gamma}), \qquad (3) $$

see [26]. The volume fraction of RBCs in the suspension, namely the hematocrit ($Ht$), can be calculated using the volume $V'$ and the number of RBCs $N$ in a system with the volume $V$ as follows: $Ht = N \times \frac{V'}{V}$. The hematocrit was set to 25 % to 45 %, matching the hematocrit observed in healthy subjects and GD patients [8]. We performed the simulations over long timescales, and the shear stress was then averaged over the entire steady-state flow. Plasma viscosity was excluded as relevant protein levels were normal [8]. We use relative viscosity (RBC suspension viscosity normalized by the viscosity of the plasma) to assess the effect of RBCs on blood rheology. Dimensionless parameters are provided in the Supporting Information (as shown in S3 Text).

## Simulation setup

In the viscosity simulations, illustrated in Fig 1A, we modeled a rectangular box with dimensions of 60 µm in the X-direction, 60 µm in the Y-direction, and 50 µm in the Z-direction (perpendicular to the XY plane) and bounded by two walls of width 5 µm. The simulation box considered in the current study (Fig 1A) has been shown to be sufficiently large for accurate rheological characterization [27]. The shear flow was imposed by the bounding walls moving in opposite directions at a constant velocity. The boundary conditions were periodic in the velocity (X) and the vorticity (Z) directions, and no-slip boundary condition was applied near the walls. Each individual RBC was composed of 500 DPD particles, with the number of RBCs modeled ranging from 155 to 698 to reproduce hematocrits of $Ht = 10–45\%$. Accordingly, 600,000 particles were modeled to represent the plasma, resulting in a total number of particles ranging from 677,500 to 949,000.

All simulations were computed by using the extended version of code developed based on LAMMPS. Each simulation takes about 2,000,000 time steps to 4,000,000 time steps. A typical simulation requires 2400 CPU core hours to 4800 CPU core hours by using the computational resources (Intel Xeon E5-2670 2.6 GHz 24-core processors) at the Center for Computation and Visualization at Brown University.

## Experiment setup

For GD RBC deformation during the adhesion assay, the experiment was performed using a plate flow chamber and an ibiTreat-coated or uncoated microslide 10.2 Luer (Ibidi) [8]. Ten µg of laminin α5 chain was immobilized on an uncoated microslide and RBCs, resuspended at 0.5% hematocrit (Hct), were perfused through the microslide under a 0.2 dyn/cm² flow for 5 minutes at 37°C. After cells were washed out, RBCs were fixed at 1 dyn/cm², and RBC
morphology and Lu/BCAM pattern were visualized by immunofluorescence using a specific antibody and confocal microscopy [8]. We qualified the shear modulus for GD RBCs by determining the X-axis length of CTR and GD RBCs under a shear stress of 1 dyn/cm² from the confocal images of adherent CTR and GD RBCs on laminin α5. The X-axis length measurements were performed on several images obtained from adhesion assay experiments conducted at a shear stress of 1 dyn/cm², using CTR (n = 10) or GD (n = 11) RBCs. The X-axis corresponds to the flow direction. The mean X-axis length of CTR RBCs is 11.0 µm (n = 10), while the mean X-axis length of GD RBCs is 8.9 µm (n = 11) (Fig 2B).

The RBC aggregation experiment was conducted at 37°C by syllectometry (i.e., laser backscatter versus time, using the Laser Optical Rotational Red Cell Analyzer; LORRCA, Mechatronics, The Netherland) after adjustment of the hematocrit to 40%. The RBC disaggregation threshold (γ), that is, the minimal shear rate needed to prevent RBC aggregation or to breakdown existing RBC aggregates, was determined using a reiteration procedure [8]. RBC deformability was measured by ektacytometry (LORRCA, Mechatronics, The Netherland) at shear stresses ranging from 0.53 to 3.00 Pa, and at 37°C [8].

Blood viscosity experiment was setup at native Hct, at 37°C, and at various shear rates using a cone/plate viscometer (Brookfield DVIII+ with CPE40 spindle, Brookfield Engineering Labs, USA). All hemorheological measurements were carried out within 4 hours of venipuncture to avoid blood rheological alterations and after complete oxygenation of the blood. We followed the guidelines for international standardization of blood rheology techniques/measurements and interpretation [49].

## Results

In this section, we first validated the simulation parameters against experimental data. Then, we simulated suspensions of RBCs that replicate the experimental conditions of GD erythrocyte suspensions as conducted by Franco et al. [8]. Last, we investigated the effects of hematocrit level, RBC deformability, and RBC disaggregation threshold on blood viscosity by DPD simulation. Detailed simulation value are provided in the Supporting Information (as shown in S1 Data).

## Clinical and hemorheological profile of GD patients with bone disease

This study focused on acute or subacute manifestations of GD such as bone infarcts and AVN. We conducted an analysis on 22 untreated GD patients, among whom 7 exhibited various forms of bone disease, including bone infarcts and vertebral body collapse, indicating a prevalence of bone disease in 31.8% of GD patients. Among these 22 patients, hemorheology experiments were performed on 9 patients [8]. Table 2 provided a detailed overview of the clinical and biological parameters of 3 GD patients who have been diagnosed with bone disease. Each RBC parameter from these 3 patients was compared with mean values calculated from 11 healthy volunteers (M-CTR) and from 9 GD patients (M-GAU) (Table 2). Notably, both Patient 1 and Patient 3 demonstrated substantially lower elongation indices (calculated as

$$EI = \frac{L - W}{L + W}$$

where L is the cell's length along the flow axis and W its width perpendicular to the flow) compared to the M-CTR and M-GAU across shear stress levels ranging from 0.3 Pa to 3.0 Pa. Patient 2 exhibited lower elongation index values compared to the M-GAU values at shear stress level ranging from 0.53 to 0.3 Pa, indicating a significant decrease in RBC deformability.
<table><caption>Table 2. Clinical and biological parameters of the studied patients with hemorheology experiment.</caption><thead><tr><th>Patient No.</th><th>1</th><th>2</th><th>3</th><th>M-CTR</th><th>M-GAU</th></tr></thead><tbody><tr><td>Clinical skeletal manifestations</td><td>-</td><td>-</td><td>bone pains</td><td>-</td><td>-</td></tr><tr><td>Radiological (MRI) skeletal manifestations</td><td>VBC</td><td>BI</td><td>BI</td><td>-</td><td>-</td></tr><tr><td>Hb (g/dl)</td><td>13.9</td><td>12.5</td><td>12.0</td><td>-</td><td>12.6</td></tr><tr><td>EI at 0.3 Pa</td><td>0.004</td><td>0.043</td><td>0.002</td><td>0.055</td><td>0.021</td></tr><tr><td>EI at 0.53 Pa</td><td>0.055</td><td>0.055</td><td>0.053</td><td>0.084</td><td>0.062</td></tr><tr><td>EI at 0.95 Pa</td><td>0.113</td><td>0.114</td><td>0.123</td><td>0.165</td><td>0.131</td></tr><tr><td>EI at 1.69 Pa</td><td>0.202</td><td>0.216</td><td>0.226</td><td>0.262</td><td>0.229</td></tr><tr><td>EI at 3.0 Pa</td><td>0.303</td><td>0.317</td><td>0.331</td><td>0.356</td><td>0.329</td></tr><tr><td>RBC disaggregation threshold (s<sup>-1</sup>)</td><td>800</td><td>190</td><td>215</td><td>111</td><td>245</td></tr><tr><td>Blood viscosity at 22.5 s<sup>-1</sup></td><td>11.85</td><td>7.46</td><td>8.17</td><td>8.03</td><td>9.01</td></tr><tr><td>Blood viscosity at 45 s<sup>-1</sup></td><td>9.6</td><td>6.39</td><td>6.9</td><td>6.62</td><td>7.43</td></tr><tr><td>Blood viscosity at 90 s<sup>-1</sup></td><td>7.89</td><td>5.39</td><td>5.93</td><td>5.74</td><td>6.32</td></tr><tr><td>Blood viscosity at 225 s<sup>-1</sup></td><td>6.26</td><td>4.9</td><td>4.86</td><td>4.98</td><td>5.40</td></tr></tbody></table>

Notes: M-CTR: Mean value of 11 healthy volunteers; M-GAU: Mean value of 9 GD patients; BI: patients whose skeletal MRI revealed at least 1 bone infarct, either painful or not in the disease history; VBC: vertebral body collapse; Hb: hemoglobin; EI: Elongation Index.

https://doi.org/10.1371/journal.pcbi.1012705.t002

Furthermore, all three patients display higher RBC disaggregation thresholds and blood viscosities compared to the M-CTR, with only Patient 1 exceeding the M-GAU threshold. These data were insufficient to conduct robust statistical analysis regarding potential links between hemorheological parameters and acute/subacute bone disease in GD. However, the markedly decreased elongation indices, elevated RBC disaggregation thresholds, and increased blood viscosity suggest underlying biophysical mechanisms in these complications encouraging a deeper exploration of biophysics and rheology of blood and RBCs in GD.

## Calibration of the GD RBC deformability during adhesion assay

Membrane deformation of GD RBCs during flow adhesion experiment was visualized by confocal microscopy. At 1 dyn/cm², CTR RBCs were elongated and underwent a deformation of their membrane, with RBCs exhibiting a typical racket shape described by Franco et al. [8]. In contrast, GD RBCs did not change shape, suggesting reduced membrane deformability (left column of Fig 2B). X-axis lengths were measured using CTR and GD RBC images obtained from adhesion assay (Fig 2A).

In analogy with the experimental setups [8], simulations were conducted where a single RBC was initially attached to a substrate, and a pressure gradient was imposed on the liquid, inducing shear flow on the RBC. The direction of the shear flow followed the X-axis, from left to right, as shown in Fig 2B. Subsequently, a stepwise increasing pressure gradient was applied to the liquid until the X-axis length of the CTR RBC in the simulation matched that of the CTR RBC in the experiment, i.e., 11.6 µm. At this point, the simulated RBC and the experimental RBC exhibited similar morphology (top column of Fig 2B). We determined that the pressure gradient parameter in this simulation corresponded to a shear stress of 1 dyn/cm². After establishing this shear stress, we maintained it at 1 dyn/cm² while varying the shear modulus of the RBC in simulations from 1Es₀ to 25Es₀. After completing the simulations for each shear modulus of the GD RBC, we compared the X-axis length of each RBC under shear stress to the X-axis length of the GD RBC in the experiment (S1 Fig). We found that the RBC with a shear modulus of 9Es₀ closely matched the experimental GD RBC in terms of X-axis length and morphology. The morphology of the RBC in the simulation closely resembled that
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Validation of shear modulus of GD RBC with adhesion assay. A) X-axis length of CTR and GD RBCs from experiments at a shear stress of 1 dyn/cm². The X-axis length measurements were performed on several images obtained from adhesion assay experiments conducted at a shear stress of 1 dyn/cm², using CTR (n = 10) or GD (n = 11) RBCs. B) Morphology of CTR and GD RBC under shear flow. The left two figures are confocal images of adherent CTR and GD RBCs to laminin α5 at a shear stress of 1 dyn/cm². The right two figures are simulation images for CTR and GD RBCs at different shear modulus conditions at a shear stress of 1 dyn/cm². Bar scale represents 1 µm.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g002

of the RBC observed in the experiment, indicating that the simulation parameters accurately reflected the experimental conditions (bottom row of Fig 2B).

## Calibration of the aggregation model based on syllectometry

In analogy with experimental setups [8], simulations were conducted with two RBCs initially in contact one with another forming a doublet. After sufficient contact time and the cell-cell interaction reaching an equilibrium, we imposed the bounding walls moving in opposite directions at a constant velocity. The shear flow was induced by the moving of bounding walls. Subsequently, we imposed a stepwise increasing shear rate on the bounding walls until the doublet broke up, i.e., these two RBCs were fully detached from each other. In the simulation, cell-cell aggregation between the two stacked RBCs occurred after they approached each other. The setup involved two parallel surfaces moving in opposite directions to generate shear flow along the horizontal direction. As a result, the two aggregated RBCs were sheared apart by sliding (Fig 3). By using this method, we could calibrate the pairwise interaction strength in cell-cell aggregation models with the RBC disaggregation threshold observed in the experiment.

To simulate the cell-cell detachment of GD RBC doublets under different shear rates, we also carried out the simulations on the breakup of GD RBC doublets and recorded the minimal shear rate that allowed the full separation of these two RBCs in doublets (Fig 3). Generally, the cell-cell disaggregation threshold increased as the interactive vertices ratio (aggregate severity) increased in both viscometer experiments and simulations. Also, larger $D_e$ introduced the condition that a higher shear rate was required to break up the doublets, i.e., The cases with $D_e = 4.24 \times 10^{-25}$ J required higher detaching shear rates compared with those with $D_e = 8.48 \times 10^{-26}$ J. The increasing RBC disaggregation threshold required to fully detach RBC doublets in the simulation closely resembled the results from cone/plate viscometer experiments, indicating that simulation parameters accurately captured the disaggregation dynamics observed experimentally (Fig 3). Effect of RBC Stiffness on Disaggregation Threshold are provided in the Supporting Information (as shown in S4 Text).
<figure><img src="image_6.png" /><figcaption>Fig 3. Sequential snapshots of the disaggregation dynamics of GD RBCs (GD-RBC1) for a range of varied flow shear rates.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g003

## Modeling of the rheology of GD RBC suspensions

Franco et al. demonstrated increased bulk viscosity in GD RBC suspensions under physi-
logical shear rates (22.5–225 s⁻¹) using a cone/plate viscometer, [8]. For the *in silico* viscosity
measurement, we used our model, taking into account non-Newtonian behavior such as shear
thinning and the Fahraeus–Lindqvist effect [50].

To evaluate the accuracy of our *in silico* viscosity measurements, we compared our values
generated with viscosities measured by Franco et al. [8]. Fig 4 shows the steady shear blood
viscosity as a function of shear rate, which shows an excellent agreement between the sim-
ulation results and the experimental measurements for normal and GD RBC suspensions
[8].

We then extended our analysis to examine the relative blood viscosity in GD under shear
flows with larger shear rates, ranging from 1 to 1000 s⁻¹. This four-orders-of-magnitude
range covers most of the shear levels present in the micro- and macrocirculation, from

<figure><img src="image_7.png" /><figcaption>Fig 4. Blood relative viscosity for GD evaluated in viscometric flows under physiologically relevant shear rates ranging from 22.5 to 225 s⁻¹. Experimental viscosity measurements for GD and CTR RBC suspension by Franco et al. [8]. The viscosity measurements were based on a hematocrit of ~36% for GD RBCs and ~40% for CTR RBCs. Open symbols represent experimental results, whereas solid symbols represent simulation results.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g004
low/moderate rates in veins (around $1 \text{ s}^{-1}$) and arteries (around $80-250 \text{ s}^{-1}$) to high rates (up to $1000 \text{ s}^{-1}$) in capillaries and arterioles [51–53]. In Fig 5, viscosity curves were plotted for normal RBC (CTR-RBC) and GD-RBC1, GD-RBC2, and GD-RBC3, which represent RBCs from GD patients with moderate, severe, and mild conditions, respectively, as shown in Table 1. These curves revealed three distinct regions that we have categorized as separate domains. Under physiologically high shear rates ($>\sim 100 \text{ s}^{-1}$), the elevated blood viscosity of GD patients seems to be representative of the blood rheological anomalies of various vascular diseases such as atherosclerosis [54], and malaria-infected blood [55]. Under physiologically lower shear rates ($\le 8 \text{ s}^{-1}$), the blood viscosity of GD patients could also be representative of these diseases, particularly in diabetic mellitus because of the RBC hyperaggregation [56].

To investigate whether different underlying mechanisms dominate across varying flow domains in GD blood, similar to other blood rheological anomalies, we defined three distinct domains of GD blood viscosity based on shear rate. Then, we investigated whether different underlying mechanisms dominate across varying flow domains in blood from GD patients.

**Impact of the RBC hyperaggregation in GD patients on blood viscosity.** Here, we studied the role of the RBC disaggregation threshold on blood viscosity by increasing the number of RBC constituents that are attractive to one another which will induce inter RBCs attachments. Fig 6A shows the schematic representation of different RBC disaggregation thresholds through the Morse potential between constituents of different RBCs. This methodology accurately represents the increased effective interactions without prohibitive computational costs, which is consistent with the results of Deng et al. [57]. Fig 6B shows the relative viscosity of blood for different hematocrit levels (Ht = 25%, 36%, 45%) over a range of shear rates (1 to $1000 \text{ s}^{-1}$) for different RBC disaggregation thresholds. The shear rates and hematocrits examined in our study were consistent with the physiological ranges observed in the circulatory system of GD patients [8,58]. Here, we added an intermediate threshold of $350 \text{ s}^{-1}$ to reduce the large gap between $250 \text{ s}^{-1}$ and $680 \text{ s}^{-1}$. The relative viscosity plotted revealed a stark

<figure><img src="image_8.png" /><figcaption>**Fig 5.** Extension of blood relative viscosity for GD in shear flows under large shear rates ranging from 1 to $1000 \text{ s}^{-1}$. There are three distinct domains for blood viscosity based on the shear rate: the aggregation domain, ranging up to approximately $8 \text{ s}^{-1}$; the transition area, spanning from around $8$ to $100 \text{ s}^{-1}$; and the stiffness domain, extending from approximately $100 \text{ s}^{-1}$.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g005
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Schematic diagram of different RBC disaggregation thresholds (A) and corresponding simulation results (B and C) of the relative viscosity of blood for different RBC disaggregation thresholds and hematocrit levels with shear modulus of $9E_8$. A) Representation of three different glucocerebroside accumulation through the Morse potential between constituents of different RBCs. Here, the presence of more blue vertices corresponds to a stronger RBC disaggregation threshold. The severe condition corresponds to the highest RBC disaggregation threshold of $680 \text{ s}^{-1}$ in the experiment; the moderate condition corresponds to a threshold of $250 \text{ s}^{-1}$; and the mild condition corresponds to a threshold of $150 \text{ s}^{-1}$, as defined in Table 1; B) The steady shear relative viscosity versus imposed shear rate for three different hematocrit levels and four different RBC disaggregation thresholds at each hematocrit. The viscosities are normalized by the plasma viscosity; the effect of RBC disaggregation thresholds disappears at high shear rates when it causes no stable cell-cell adhesion; C) Relative viscosity against hematocrit for different shear rates and four different RBC disaggregation thresholds of 150, 250, 350 and $680 \text{ s}^{-1}$. In the simulation, we added an intermediate threshold of $350 \text{ s}^{-1}$ to reduce the large gap between $250 \text{ s}^{-1}$ and $680 \text{ s}^{-1}$.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g006

increase of viscosity at low shear rates for both hematocrit levels and RBC disaggregation thresholds studied. The observed variation in viscosity between the different colored curves is due to the different hematocrit levels tested, which shows that a high hematocrit level would contribute to high blood viscosity (Fig 6B). The discrepancy in the same color is caused by the RBC disaggregation threshold, which is much more significant at low shear rates. In contrast, virtually no change is observed at larger shear rates. Thus, the first domain of the curve is dependent on the aggregation properties of the GD RBCs.

Additionally, the increase in viscosity at low shear rates was significantly more pronounced at higher hematocrit values. This is clearly illustrated by the red lines in Fig 6C, where the relative viscosity is plotted against hematocrit for various shear rates and RBC disaggregation thresholds. In contrast, when the relative viscosity is plotted at high shear rates, the resulting lines overlap, indicating that the disaggregation threshold has no significant influence on viscosity at higher shear rates, regardless of the hematocrit level (blue and grey lines in Fig 6C).

**Impact of the cell deformability on blood viscosity in GD.** Here, we directly manipulated the deformability of RBCs in isolation by adjusting the shear modulus that maintains RBC mechanical integrity. The viscosity is reported for three distinct cases of RBC deformability/shear moduli in Fig 7. Fig 7A depicts the representation of three different levels of cell stiffness through different shear modulus conditions at a shear stress of 1 dyn/cm$^2$. Fig 7B and 7C illustrate an overall increase in viscosity across different hematocrits and shear rates, with this augmentation notably more pronounced under physiologically higher shear rates (100–1000 s$^{-1}$, see red lines in Fig 7B and grey lines in Fig 7C). We also noted that blood viscosity in
<figure><img src="image_10.png" /><figcaption>Fig 7. Schematic diagram of different RBC stiffness (A) and corresponding simulation results (B and C) of the relative viscosity of blood for different RBC stiffness and hematocrit levels with disaggregation threshold of 250s⁻¹. A) Representation of three different levels of cell stiffness through different shear modulus conditions at a shear stress of 1 dyn/cm². B) The steady shear relative viscosity versus imposed shear rate for three different hematocrit levels and three different cell stiffness at each hematocrit. The viscosities are normalized by the plasma viscosity; the effect of cell stiffness increases at high shear rates, and the deformation of cells becomes increasingly more dominant C) Relative viscosity against hematocrit for different shear rates and three different RBC stiffness of 20, 9, and 4 E₅₀.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012705.g007

GD, as predicted by the model, rises by approximately 1.4-fold compared with normal values when hematocrit increases from 25% to 45% (Fig 7B and 7C). Then, the third domain of the curve is dependent on the deformability of the RBCs.

At very low shear rates, blood rheology and flow dynamics are predominantly influenced by fluid flow and stress fluctuations between particles, explainable through suspension mechanics. Conversely, at higher shear rates, cell deformation becomes increasingly dominant. Thus, blood containing RBC with decreased deformability exhibits heightened viscosities at larger shear rates compared to blood containing normally deformable RBC. This effect is stronger at higher hematocrit levels.

**Impact of the hematocrit levels on blood viscosity in GD.** In this section, we investigated the impact of hematocrit levels on blood viscosity, in combination with the deformability and aggregation properties of GD RBCs. Thus, we changed the hematocrit from 25% to 45% and explored the non-Newtonian viscosity over a range of shear rates in GD. Figs 6 and 7 present the steady shear relative blood viscosity as a function of shear rate and hematocrit. As generally observed in physiology and disease, at all shear rates blood viscosity increased dramatically with rising hematocrit, as shown in Figs 6C and 7C; however, this increase also depended on the shear rate at which the viscosity was measured. For instance, in Fig 7C, increasing the hematocrit from 25 to 45% at the shear rate of 1 s⁻¹ resulted in a 1.8-fold higher value of blood viscosity, whereas at a shear rate of 1000 s⁻¹, it only increased 1.4-fold the value of blood viscosity. Thus, the effect of hematocrit on blood viscosity is greater at low shear rates.

In GD patients, hematocrit levels are often below the normal range, so we expected decreased blood viscosity. The average hematocrit level in blood samples from 11 untreated GD patients was 36.42% ± 3% (n = 11) [8]. This anemia can be explained by the dyserythropoiesis observed in GD and by the splenic retention of rigid GD RBCs, which results in a drop in hematocrit. However, treated patients may have a hematocrit level of 40%–45%. At low
shear rates, this increase in hematocrit to 45% may result in blood viscosity roughly doubling compared to that of healthy individuals. Similarly, at high shear rates, the viscosity may rise to approximately 1.8 times that of healthy individuals (Fig 6B). The hematocrit levels in Figs 6 and 7 (25%-45%) are within the physiological range for GD. While hematocrit significantly impacts blood viscosity, other factors also contribute to the overall viscosity profile.

Discussion and summary

Bone complications in GD are of multiple types and mechanisms, including bone fragility marked by osteopenia/osteoporosis, osseous deformities, and cortical thinning in long bones [9,10,59,60]. Additionally, acute or subacute occurrences such as bone infarcts and AVN of the head of long bones or vertebral bodies contribute to its complexity. Bone infarcts and AVN occur when bones do not get enough blood due to decreased intraosseous blood flow, causing the tissue to deteriorate and die, resulting in tiny breaks in the bone and eventual bone collapse [10,11]. Although reminiscent of a similar complication in sickle cell disease, where microcirculation impairment is central [12], the mechanisms that contribute to bone involvement in GD are not clearly established. Several studies have shown that bone pathology in GD is attributed to infiltration of the bone marrow and trabecular bone by Gaucher cells [9,13]. On the other hand, it has been shown in different pathologies that blood viscosity abnormalities can induce ischemic processes [61,62]. Detailed analysis of rheological data from 3 Gaucher patients with bone damage suggests that the RBC aggregation properties, as well as their deformability and blood viscosity, could play a role in the onset of bone damage. Furthermore, a large study from the International Gaucher Disease Registry (ICGG) indicates that patients with GD who have undergone surgical removal of the spleen (splenectomy) have a higher incidence of AVN [18]. Given that the spleen is the organ that removes biophysically altered RBCs from circulation, these findings strongly suggest a potential link between AVN occurrence and altered RBC biophysics and blood viscosity [63,64].

Based on the rheological alterations of RBCs in GD, we hypothesized that the pathogenesis of bone infarcts and AVN in GD may involve rigid erythrocytes and elevated RBC aggregation, leading to increased blood viscosity and microcirculatory obstruction [65-69]. It would therefore seem interesting to study more precisely the impact of abnormal Gaucher RBC properties, such as aggregation and deformability, on blood viscosity. Due to the limitations of experimental *in vitro* measurements mentioned earlier, we chose to use computational methods, as advances in computational power over the past two decades have stimulated the development and applications of multiscale biophysical models [23,70-75]. Based on previous results obtained *ex vivo* [8], in this study, we investigated the biomechanics and rheology of blood and RBCs in GD under various flow conditions via *in silico* models.

In our study, we identified three factors that affect blood viscosity in GD: RBC deformability, RBC disaggregation threshold, and hematocrit level. Firstly, experimental findings have shown that hardened RBCs are present in GD patients [8], and these cells exhibit increased resistance to flow compared to healthy ones, leading to elevated blood viscosity [19]. Reduced RBC deformability is also observed in conditions like diabetes mellitus, where shear moduli are twice those of healthy cells, contributing to higher viscosity across various shear rates [35]. Similarly, malaria-infected RBCs lose deformability, further increasing viscosity [76], while in sickle cell anemia, deoxygenated sickled RBCs can become 100-2000 times more rigid, significantly raising blood viscosity [77]. Secondly, the mechanisms underlying RBC hyperaggregation in GD patients remain unclear but likely involve cellular factors [17,78]. Glucocerebroside, a sphingolipid that accumulates in RBC membranes, may play a key role by raising the RBC disaggregation threshold and promoting rouleaux formation, where cells
stack in chain-like structures [17]. Lastly, hematocrit levels are lower in GD patients com-
pared to healthy controls [20], which highlights the importance of considering hematocrit
variations when assessing splenectomy in these patients. Post-splenectomy hemorheological
changes have been linked to bone complications [79–81].

To investigate the effects of RBC deformability, RBC disaggregation threshold, and hemat-
ocrit level on blood viscosity by DPD simulation, we first validated the simulation parameters
against experimental data. This study is the first to investigate the Young's modulus and aggre-
gation parameters of GD RBCs by integrating validated simulations with confocal imaging
and experimental RBC disaggregation thresholds. Our simulation parameters accurately cap-
tured the shear modulus of GD RBCs and the RBC disaggregation dynamics observed exper-
imentally. This approach enables a comprehensive assessment of the mechanical properties of
GD RBCs.

To validate the *in silico* viscosity model against experimental measurements, we simulated suspensions of RBCs that replicate the experimental conditions of GD erythrocyte suspensions as conducted by Franco et al. [8]. Our simulation shows an excellent agreement with the experimental measurements. Experimental studies have shown glucosylceramide accumulation in both the RBC membrane and plasma in GD-RBCs, which may enhance aggregation [8,33]. While our model focuses on membrane-mediated adhesion, we acknowledge that plasma factors also contribute. Due to RBC compatibility and handling challenges, we do not experimentally separate these effects in the current study. After that, we extended our analysis to examine the relative blood viscosity in GD with larger shear rates ranging from 1 to 1000 s⁻¹.

Furthermore, to explore whether different underlying mechanisms dominate across vary-
ing flow domains in GD blood, we defined three distinct domains of GD blood viscosity based
on shear rate. These domains include the aggregation domain, which spanned shear rates
up to approximately 8 s⁻¹; the transition area, which covered shear rates from around 8 to
100 s⁻¹; and the stiffness domain, which extended from approximately 100 s⁻¹ and beyond.
The stiffness domain may correspond to the regions where RBC deformability becomes a
major factor.

Our findings regarding the effect of RBC stiffening on blood viscosity show notable sim-
ilarities to the observations of Lanotte et al. [82]. Specifically, both studies demonstrate that
reduced deformability of RBCs leads to an increase in blood viscosity and modulates shear-
thinning behavior under microcirculatory flow conditions. In agreement with Lanotte et al.,
we observe that hardening of RBCs alters their ability to deform dynamically in response to
shear rate, resulting in elevated viscosity [82]. These consistent trends strengthen the criti-
cal role of RBC mechanical properties in governing hemorheological behavior under phys-
iological shear rates. In addition, at low shear rates (<10 s⁻¹), RBC aggregates remain intact
but deform and align with the flow, leading to increased viscosity due to microstructural
rearrangements [50,83]. In the intermediate range (10–100 s⁻¹), aggregates partially elongate
without full breakup, resulting in a dynamic equilibrium. Some studies report that in this
range, microstructural differences between aggregating and non-aggregating RBC suspen-
sions become negligible, and viscosity converges across media [34,50]. In GD, altered RBCs
could lead to reduced deformability and the ability to navigate through small capillaries, wors-
ening the microcirculatory problems and contributing to increased vascular resistance and
even organ dysfunction, such as bone infarction. Consequently, we can imagine that initially,
reduced deformability of RBCs could induce a slowing of blood flow, which would then be
even more impacted by the high RBC disaggregation threshold.

Finally, we investigated the effects of RBC deformability, RBC disaggregation threshold,
and hematocrit level on blood viscosity. Our study demonstrates that all three mechanisms
can cause blood viscosity to rise many folds, with manifestation in different shear rate con-
ditions and hematocrit levels. An increased RBC disaggregation threshold leads to higher
blood viscosity in GD patients at low shear rates. Consequently, this phenomenon is less likely
to occur in arterioles and capillaries, which experience high shear rates. In terms of RBC
deformability, decreased RBC deformability can cause elevation of blood viscosity, this can
be captured under large shear rates, and observed no change at low shear rates. However, this
phenomenon is less common in veins, which typically experience lower shear rates.

In summary, we identified three viscosity domains based on shear rate, showing that ele-
vated hematocrit, reduced RBC deformability, and increased RBC disaggregation threshold
contribute to increased blood viscosity with manifestation in different shear rate conditions
in GD. In GD patients, lowering hematocrit helps protect the circulatory system from compli-
cations associated with high blood viscosity. Treated patients may reach hematocrit levels of
40%-45%, which could result in high blood viscosity if RBC anomalies were maintained. For-
tunately, treatment significantly reduces RBC abnormalities, and as splenomegaly decreases,
the spleen regains functionality, aiding in the return of blood viscosity to normal [84]. How-
ever, splenectomy can negatively impact this process by increasing hematocrit and releasing
hyperaggregated, rigid GD RBCs into circulation. This underscores the role of "saving the
spleen" objective in maintaining optimal blood rheology in GD patients. This study enhances
our understanding of the biophysics and rheology of blood and RBCs in GD, offering insights
into the broader implications of these factors in pathological conditions.

Our study of blood biophysics and rheology in GD has some limitations. While RBC
membrane parameters were calibrated with experimental data, assuming fixed shear modulus
and RBC disaggregation thresholds may not account for variability across individuals or GD
stages. Additionally, the model reflects RBC dynamics in controlled environments, which may
not fully capture *in vivo* conditions where factors like immune cell interactions affect blood
rheology. Further mapping these simulations to real physiological contexts could provide
deeper insights into GD-related hemorheological phenomena.

While our study focuses on how RBC aggregation and stiffness affect blood viscosity, the
resulting impact on microvascular flow and obstruction remains an important open ques-
tion. Whether microcirculatory blockage arises primarily from RBC clustering or reduced
deformability is not addressed here, but this distinction may play a critical role in the patho-
genesis of bone infarcts and AVN in Gaucher Disease. Future studies integrating detailed flow
simulations and experimental validation are needed to further explore this mechanism. In
our current study, direct quantification of RBC morphology across shear rates is not included.
In future work, a detailed study of cell alignment, the relation between EI and RBC stiff-
ness *in silico*, different shapes at varying hematocrits and deformation rates, and mapping of
such structures with respect to different biophysical factors at play will provide an invaluable
insight into the groundwork of hemorheological phenomena of GD.

**Supporting information**

**S1 Text. Dissipative particle dynamics (DPD) method and particle-based blood cell models.**

(PDF)

**S2 Text. Cell-cell interaction potentials and measurements of membrane stiffness.**

(PDF)
S3 Text. Dimensionless parameters.

(PDF)

S4 Text. Effect of RBC stiffness on disaggregation threshold.

(PDF)

S1 Data. GD biophysics and rheology data set.

(ZIP)

S1 Fig. Morphology of RBC in simulations under different shear modulus conditions.

(TIF)

S2 Fig. Influence of RBC stiffness on disaggregation threshold.

(TIF)

Acknowledgments

Simulations were carried out at the Center for Computation and Visualization of Brown University.

References

1. Watek M, Piktel E, Wollny T, Durnas B, Fiedoruk K, Lech-Maranda E, et al. Defective sphingolipids metabolism and tumor associated macrophages as the possible links between gaucher disease and blood cancer development. Int J Mol Sci. 2019;20(843):10–3390.

2. Grabowski GA. Phenotype, diagnosis, and treatment of Gaucher’s disease. The Lancet. 2008;372(9645):1263–1271.

3. Baris HN, Cohen IJ, Mistry PK. Gaucher disease: the metabolic defect, pathophysiology, phenotypes and natural history. Pediatric endocrinology reviews: PER. 2014;12(0 1):72.

4. Stirnemann J, Belmatoug N, Camou F, Serratrice C, Froissart R, Caillaud C, et al. A Review of Gaucher Disease Pathophysiology, Clinical Presentation and Treatments. International Journal of molecular sciences. 2017;18(2):441.

5. Grabowski GA, Antommaria AHM, Kolodny EH, Mistry PK. Gaucher disease: Basic and translational science needs for more complete therapy and management. Molecular Genetics and Metabolism. 2021;132(2):59–75.

6. Boven LA, van Meurs M, Boot RG, Mehta A, Boon L, Aerts JM, et al. Gaucher cells demonstrate a distinct macrophage phenotype and resemble alternatively activated macrophages. American journal of clinical pathology. 2004;122(3):359–369.

7. Dupuis L, Chauvet M, Bourdelier E, Dussiot M, Belmatoug N, Le Van Kim C, et al. Phagocytosis of erythrocytes from gaucher patients induces phenotypic modifications in macrophages, driving them toward gaucher cells. International Journal of Molecular Sciences. 2022;23(14):7640.

8. Franco M, Collec E, Connes P, van den Akker E, Billette de Villemeur T, Belmatoug N, et al. Abnormal properties of red blood cells suggest a role in the pathophysiology of Gaucher disease. Blood, The Journal of the American Society of Hematology. 2013;121(3):546–555.

9. Hughes D, Mikosch P, Belmatoug N, Carubbi F, Cox T, Goker-Alpan O, et al. Gaucher disease in bone: from pathophysiology to practice. Journal of Bone and Mineral Research. 2019;34(6):996–1013.

10. Lafforgue P. Pathophysiology and natural history of avascular necrosis of bone. Joint Bone Spine. 2006;73(5):500–507.

11. Binnetoglu E, Komurcu E, Sen H, Kizildag B. Gaucher disease with pathological femoral neck fracture. Case Reports. 2013;2013:bcr2013200260.

12. Itzep NP, Jadhav SP, Kanne CK, Sheehan VA. Spontaneous healing of avascular necrosis of the femoral head in sickle cell disease. Am J Hematol. 2019;94(6):E160–E162.

13. Schiffmann R, Mankin H, Dambrosia JM, Xavier RJ, Kreps C, Hill SC, et al. Decreased bone density in splenectomized Gaucher patients receiving enzyme replacement therapy. Blood Cells, Molecules, and Diseases. 2002;28(2):288–296.

14. Adar T, Ben-Ami R, Elstein D, Zimran P, Berliner S, Yedgar S, et al. Increased red blood cell aggregation in patients with Gaucher disease is non-inflammatory. Clinical hemorheology and microcirculation. 2008;40(2):113–118.
15. Reihani N, Arlet JB, Dussiot M, de Villemeur TB, Belmatoug N, Rose C, et al. Unexpected macrophage-independent dyserythropoiesis in Gaucher disease. Haematologica. 2016;101(12):1489.

16. Chipeaux C, de Person M, Burguet N, de Villemeur TB, Rose C, Belmatoug N, et al. Optimization of ultra-high pressure liquid chromatography-tandem mass spectrometry determination in plasma and red blood cells of four sphingolipids and their evaluation as biomarker candidates of Gaucher's disease. Journal of Chromatography a. 2017;1525:116-125.

17. Adar T, Ben-Ami R, Elstein D, Zimran A, Berliner S, Yedgar S, et al. Aggregation of red blood cells in patients with Gaucher disease. British journal of haematology. 2006;134(4):432-437.

18. Khan A, Hangartner T, Weinreb NJ, Taylor JS, Mistry PK. Risk factors for fractures and avascular osteonecrosis in type 1 Gaucher disease: A study from the International Collaborative Gaucher Group (ICGG) Gaucher Registry. Journal of Bone and Mineral Research. 2012;27(8):1839-1848.

19. Shin S, Ku Y, Park MS, Suh JS. Deformability of red blood cells: A determinant of blood viscosity. Journal of Mechanical Science and Technology. 2005;19(1):216-223.

20. Bax BE, Richfield L, Bain MD, Mehta AB, Chalmers RA, Rampling MW. Haemorheology in Gaucher disease. European Journal of Haematology. 2005;75(3):252-258.

21. Paul A, Jacoby G, Laor Bar-Yosef D, Beck R, Gazit E, Segal D. Glucosylceramide Associated with Gaucher Disease Forms Amyloid-like Twisted Ribbon Fibrils That Induce -Synuclein Aggregation. ACS Nano. 2021;15(7):11854-11868.

22. Hellstrand E, Nowacka A, Topgaard D, Linse S, Sparr E. Membrane lipid co-aggregation with -synuclein fibrils. PloS one. 2013;8(10):e77235.

23. Rajabi A, Javadi E, Pejman Sereshkeh SR, Morad MR, Kebriaee A, Nasiri H, et al. Experimental characterization of an extended electrohydrodynamic cone-jet with a hemispherical nozzle. Physics of Fluids. 2018;30(11):114108.

24. Li X, Vlahovska PM, Karniadakis GE. Continuum-and particle-based modeling of shapes and dynamics of red blood cells in health and disease. Soft matter. 2013;9(1):28-37.

25. Yazdani A, Deng Y, Li H, Javadi E, Li Z, Jamali S, et al. Integrating blood cell mechanics, platelet adhesive dynamics and coagulation cascade for modelling thrombus formation in normal and diabetic blood. Journal of The Royal Society Interface. 2021;18(175):20200834.

26. Javadi E, Deng Y, Karniadakis GE, Jamali S. In silico biophysics and hemorheology of blood hyperviscosity syndrome. Biophysical Journal. 2021;120(13):2723-2733.

27. Fedosov DA, Pan W, Caswell B, Gompper G, Karniadakis GE. Predicting human blood viscosity in silico. Proceedings of the National Academy of Sciences. 2011;108(29):11772-11777.

28. Fedosov DA, Noguchi H, Gompper G. Multiscale modeling of blood flow: From single cells to blood rheology. Biomechanics and Modeling in Mechanobiology. 2014;13(2):239-258.

29. Li G, Ye T, Li X. Parallel modeling of cell suspension flow in complex micro-networks with inflow/outflow boundary conditions. Journal of Computational Physics. 2020;401:109031.

30. Li G, Qiang Y, Li H, Li X, Dao M, Karniadakis GE. In silico and in vitro study of the adhesion dynamics of erythrophagocytosis in sickle cell disease. Biophysical Journal. 2023;122:2590-2604.

31. Liu ZL, Li H, Qiang Y, Buffet P, Dao M, Karniadakis GE. Computational modeling of biomechanics and biorheology of heated red blood cells. Biophysical Journal. 2021;120(21):4663-4671.

32. Diez-Silva M, Dao M, Han J, Lim CT, Suresh S. Shape and biomechanical characteristics of human red blood cells in health and disease. MRS bulletin. 2010;35(5):382-388.

33. Deng Y, Papageorgiou DP, Li X, Perakakis N, Mantzoros CS, Dao M, et al. Quantifying Fibrinogen-Dependent Aggregation of Red Blood Cells in Type 2 Diabetes Mellitus. Biophysical Journal. 2020;119(5):900-912.

34. Fedosov DA, Lei H, Caswell B, Suresh S, Karniadakis GE. Multiscale Modeling of Red Blood Cell Mechanics and Blood Flow in Malaria. PLOS Computational Biology. 2011;7:1-13.

35. Chang HY, Li X, Karniadakis GE. Modeling of Biomechanics and Biorheology of Red Blood Cells in Type 2 Diabetes Mellitus. Biophysical Journal. 2017;113(2):481-490.

36. Baskurt O, Neu B, Meiselman HJ. Red blood cell aggregation. 2011.

37. Bagchi P, Johnson PC, Popel AS. Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow. Journal of biomechanical engineering. 2005;127(7):1070-1080.

38. Liu Y, Zhang L, Wang X, Liu WK. Coupling of Navier-Stokes equations with protein molecular dynamics and its application to hemodynamics. International Journal for Numerical Methods in Fluids. 2004;46(12):1237-1252.

39. Liu Y, Liu WK. Rheology of red blood cell aggregation by computer simulation. Journal of Computational Physics. 2006;220(1):139-154.

40. Zhang J, Johnson PC, Popel AS. Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method. Journal of biomechanics. 2008;41(1):47-55.
41. Flormann D, Aouane O, Kaestner L, Ruloff C, Misbah C, Podgorski T, et al. The buckling instability of aggregating red blood cells. Scientific reports. 2017;7(1):7928.

42. Brust M, Aouane O, Thiébaud M, Flormann D, Verdier C, Kaestner L, et al. The plasma protein fibrinogen stabilizes clusters of red blood cells in microcapillary flows. Scientific reports. 2014;4(1):4348.

43. Bäumler H, Neu B, Donath E, Kiesewetter H. Basic phenomena of red blood cell rouleaux formation. Biorheology. 1999;36:439-42.

44. Irving JH, Kirkwood JG. The Statistical Mechanical Theory of Transport Processes. IV. The Equations of Hydrodynamics. The Journal of Chemical Physics. 1950;18(6):817-829.

45. Suresh S, Spatz J, Mills JP, Micoulet A, Dao M, Lim C, et al. Connections between single-cell biomechanics and human disease states: gastrointestinal cancer and malaria. Acta biomaterialia. 2005;1(1):15-30.

46. Peng Z, Li X, Pivkin IV, Dao M, Karniadakis GE, Suresh S. Lipid bilayer and cytoskeletal interactions in a red blood cell. Proceedings of the National Academy of Sciences. 2013;110(33):13356-13361.

47. Li X, Li H, Chang HY, Lykotrafitis G, Karniadakis GE. Computational biomechanics of human red blood cells in hematological disorders. Journal of biomechanical engineering. 2017;139(2):0210081.

48. Skalak R, Keller SR, Secomb TW. ASME centennial historical perspective paper: mechanics of blood flow. 1981.

49. Baskurt O, Boynard M, Cokelet G, Connes P, Cooke BM, Forconi S, et al. New guidelines for hemorheological laboratory techniques. Clinical Hemorheology and Microcirculation. 2009;42(2):75-97.

50. Chien S. Shear Dependence of Effective Cell Volume as a Determinant of Blood Viscosity. Science. 1970;168(3934):977-979.

51. Secomb TW. Mechanics of blood flow in the microcirculation. Symposia of the Society for Experimental Biology. 1995;49:305-321.

52. Ku DN. Blood flow in arteries. Annual Review of Fluid Mechanics. 1997;29:399-434.

53. Baskurt OK, Meiselman HJ. Blood rheology and hemodynamics. In: Seminars in thrombosis and hemostasis. vol. 29. Thieme Medical Publishers, Inc.; 2003. p. 435-450.

54. Woo HG, Kim HG, Lee KM, Ha SH, Jo H, Heo SH, et al. Blood viscosity associated with stroke mechanism and early neurological deterioration in middle cerebral artery atherosclerosis. Scientific Reports. 2023;13(1):9384.

55. Scherer EF, Cantarini DG, Siqueira R, Ribeiro EB, Braga ÉM, Honório-França AC, et al. Cytokine modulation of human blood viscosity from vivax malaria patients. Acta Tropica. 2016;158:139-147.

56. Peduzzi M, Melli M, Fonda S, Codeluppi L, Guerrieri F. Comparative evaluation of blood viscosity in diabetic retinopathy. International Ophthalmology. 1984;7(1):15-19.

57. Deng Y, Papageorgiou DP, Chang H, Abidi SZ, Li X, Dao M, et al. Quantifying Shear-Induced Deformation and Detachment of Individual Adherent Sickle Red Blood Cells. Biophyscal Journal. 2019;116:360-371.

58. Nader E, Skinner S, Romana M, Fort R, Lemonne N, Guillot N, et al. Blood rheology: key parameters, impact on blood flow, role in sickle cell disease and effects of exercise. Frontiers in physiology. 2019;10:1329.

59. Mikosch P. Gaucher disease and bone. Best Practice & Research Clinical Rheumatology. 2011;25(5):665-681.

60. Wenstrup RJ, Roca-Espiau M, Weinreb NJ, Bembi B. Skeletal aspects of Gaucher disease: a review. The British Journal of Radiology. 2002;75(suppl_1):A2-A12.

61. Kang J, Oh JS, Kim BJ, Kim JY, Kim DY, Yun SY, et al. High blood viscosity in acute ischemic stroke. Frontiers in Neurology. 2023;14:1320773.

62. Akinola N, Stevens S, Franklin I, Nash G, Stuart J. Rheological changes in the prodromal and established phases of sickle cell vaso-occlusive crisis. British journal of haematology. 1992;81(4):598-602.

63. Pivkin I, Peng Z, Karniadakis GE, Buffet PA, Dao M, Suresh S. Biomechanics of red blood cells in human spleen and consequences for physiology and disease. Proceedings of the National Academy of Sciences. 2016;113:7804-7809.

64. Sloop GD, De Mast Q, Pop G, Weidman JJ, Cyr JAS. The role of blood viscosity in infectious diseases. Cureus. 2020;12(2).

65. Ballas SK, Smith ED. Red blood cell changes during the evolution of the sickle cell painful crisis. Blood. 1992;79(8):2154-2163. https://doi.org/10.1182/blood.V79.8.2154.2154
66. Sanders WJ. A rare case of avascular necrosis in sickle cell trait: a case report. BMC Hematology. 2018;18.

67. Kuznetsova PI, Raskurazhev AA, Shabalina AA, Melikhyan AL, Subortseva IN, Tanashyan MM. Red blood cell morphodynamics in patients with polycythemia Vera and stroke. International Journal of Molecular Sciences. 2022;23(4):2247.

68. Man Y, Kucukal E, An R, Watson QD, Bosch J, Zimmerman PA, et al. Microfluidic assessment of red blood cell mediated microvascular occlusion. Lab on a Chip. 2020;20(12):2086-2099.

69. Mohanty JG, Nagababu E, Rifkind JM. Red blood cell oxidative stress impairs oxygen delivery and induces red blood cell aging. Frontiers in physiology. 2014;5:84.

70. Tang YH, Lu L, Li H, Evangelinos C, Grinberg L, Sachdeva V, et al. OpenRBC: A fast simulator of red blood cells at protein resolution. Biophysical journal. 2017;112(10):2030-2037.

71. Zhang P, Zhang L, Slepian MJ, Deng Y, Bluestein D. A multiscale biomechanical model of platelets: Correlating with in-vitro results. Journal of biomechanics. 2017;50:26-33.

72. Chai Z, Gu S, Lykotrafitis G. Dynamics of the axon plasma membrane skeleton. Soft Matter. 2023;19(14):2514-2528.

73. Chai Z, Tzingounis AV, Lykotrafitis G. The periodic axon membrane skeleton leads to Na nanodomains but does not impact action potentials. Biophysical Journal. 2022;121(18):3334-3344.

74. Chang HY, Li X, Li H, Karniadakis GE. MD/DPD multiscale framework for predicting morphology and stresses of red blood cells in health and disease. PLoS computational biology. 2016;12(10):e1005173.

75. Li G, Ye T, Wang S, Li X, Haq UR. Numerical design of a highly efficient microfluidic chip for blood plasma separation. Physics of Fluids. 2020;32:031903.

76. Fedosov DA, Caswell B, Suresh S, Karniadakis GE. Quantifying the biophysical characteristics of <i>Plasmodium-falciparum</i>–<i>parasitized</i> red blood cells in microcirculation. Proceedings of the National Academy of Sciences. 2011;108(1):35-39.

77. Lei H, Karniadakis GE. Quantifying the Rheological and Hemodynamic Characteristics of Sickle Cell Anemia. Biophysical Journal. 2012;102(2):185-194.

78. Zimran A, Bashkin A, Elstein D, Rudensky B, Rotstein R, Rozenblat M, et al. Rheological determinants in patients with Gaucher disease and internal inflammation. American Journal of Hematology. 2004;75(4):190-194.

79. Englum BR, Rothman J, Leonard S, Reiter A, Thornburg C, Brindle M, et al. Hematologic outcomes after total splenectomy and partial splenectomy for congenital hemolytic anemia. Journal of pediatric surgery. 2016;51(1):122-127.

80. Ciftciler R, Pasayeva A, Aksu S, Ozcebe O, Sayinalp N, Malkan UY, et al. Indications and outcomes of splenectomy for hematological disorders. Open Medicine. 2019;14(1):491-496.

81. Baskurt OK. The role of spleen in suppressing the rheological alterations in circulating blood. Clinical Hemorheology and Microcirculation. 1999;20(3):181-188.

82. Lanotte L, Mauer J, Mendez S, Fedosov DA, Fromental JM, Claveria V, et al. Red cells' dynamic morphologies govern blood shear thinning under microcirculatory flow conditions. Proceedings of the National Academy of Sciences. 2016;113(47):13289-13294.

83. Bishop JJ, Popel AS, Intaglietta M, Johnson PC. Rheological effects of red blood cell aggregation in the venous network: a review of recent studies. Biorheology. 2001;38(2-3):263-274.

84. Dupuis L, Chipeaux C, Bourdelier E, Martino S, Reihani N, Belmatoug N, et al. Effects of sphingolipids overload on red blood cell properties in Gaucher disease. J Cell Molecul Med. 2020;24(17):9726-9736.