Marginated aberrant red blood cells induce pathologic vascular stress fluctuations in a computational model of hematologic disorders

Red blood cell (RBC) disorders such as sickle cell disease affect billions worldwide. While much attention focuses on altered properties of aberrant RBCs and corresponding hemodynamic changes, RBC disorders are also associated with vascular dysfunction, whose origin remains unclear and which provoke severe consequences including stroke. Little research has explored whether biophysical alterations of RBCs affect vascular function. We use a detailed computational model of blood that enables characterization of cell distributions and vascular stresses in blood disorders and compare simulation results with experimental observations. Aberrant RBCs, with their smaller size and higher stiffness, concentrate near vessel walls (marginate) because of contrasts in physical properties relative to normal cells. In a curved channel exemplifying the geometric complexity of the microcirculation, these cells distribute heterogeneously, indicating the importance of geometry. Marginated cells generate large transient stress fluctuations on vessel walls, indicating a mechanism for the observed vascular inflammation.

The PDF file includes:

Supplementary Text Figs. S1 to S14 Legends for movies S1 to S8 References
Other Supplementary Material for this manuscript includes the following:

Materials and Methods
We have developed a computational model utilizing the immersed boundary method (IBM) to investigate cellular blood flow in complex vessel geometries.This approach offers the advantage of modeling flows in arbitrary geometries and accurately resolving the large deformation and dynamics of blood cells.Our model considers blood as a confined flowing suspension of red blood cells (RBCs) in the plasma.Two distinct types of boundaries are involved: deformable cellular membranes and rigid non-moving vascular walls with complex geometry.To handle these two types of interfaces, we use the continuous-forcing and direct-forcing immersed boundary methods, respectively.Specifically, the continuous-forcing IBM method couples surface stresses on the RBC membrane with the fluid flow, while the sharp-interface ghost-node immersed boundary method (GNIBM) is used to treat the rigid, non-moving vascular walls (52,53).The flow solver is based on a coupled finite-volume/spectral method, and our IBM code is parallelized using a hybrid MPI/OpenMP strategy.

Governing Equations
Assuming that blood plasma is both Newtonian and incompressible, the governing equations for its flow can be expressed as the incompressible Navier-Stokes equations.

Flow Solver
The Chorin projection method is utilized to advance the velocity field u.This method involves solving an advection-diffusion equation (ADE) to determine the intermediate velocity field u * , followed by solving a Poisson equation for pressure P to enforce the divergence-free constraints.

Re
For the ADE, we treat both the convection term and the body force term with the 2nd-order Adams-Bashforth method (AB2), and the diffusion term with the Crank-Nicholson method for numerical stability.
where N (u) = −Re∇ • uu denotes the nonlinear convection, evaluated then by the central differencing.To leverage the efficient inversion of tri-diagonal matrices, a Locally One-Dimensional (LOD) Alternating Direction Implicit (ADI) scheme is employed to solve the ADE.This scheme involves four steps: first, the explicit terms are handled; then, the x, y, and z directions are solved implicitly, one at a time.To solve the 3D pressure Poisson equation (PPE), we assume that the computational domain is periodic in two directions, namely x and y.
We begin by performing a 2D fast Fourier transform (FFT) on each x − y plane.Next, solve for the Fourier coefficients along the z direction.Finally, perform a 2D inverse fast Fourier transform (iFFT) on each x − y plane to obtain the pressure field.

Spatial Discretization
The IBM requires an Eulerian representation for the fluid flow and a Lagrangian representation for the immersed boundaries such as the RBC membrane.The governing equations are solved on the Eulerian grid.To avoid odd-even decoupling, spatial discretization is based on the staggered mesh, in which these scalar variables, such as pressure P , are located at the grid center, while vector variables, such as velocity u and force F , are located at the grid faces.All spatial derivatives are evaluated using second-order differencing.As for the Lagrangian mesh on the immersed boundaries, the membrane is discretized into piecewise flat triangular elements.An open-source software Gmsh ( 64) is used to generate the surface mesh for both vessels (Fig. S1.) and various red blood cells (Fig. S2.).

Membrane Mechanics
The RBC membrane is assumed to resist shear deformation, area dilatation, volume conservation (65), and bending resistance.The total strain energy of the RBC membrane S is given by, where K B and K B are the bending moduli, and W is the shear strain energy density; κ H , κ G are the mean and Gaussian curvature of the surface, respectively; c 0 = −2H 0 is the spontaneous curvature, where H 0 is the mean curvature of the spontaneous shape.The first two terms correspond to the Canham-Helfrich bending energy (66,67), and the third term comes from the shear strain energy stored in the RBC membrane.The strain energy density is computed based on the Skalak model (68).A finite element method (FEM) is developed to find the surface stress as a result of membrane deformation (47).

Continuous Forcing IBM
The cellular membrane is deformed by the fluid flow, while the flow is altered by the membrane deformation in turn.This fluid-structure interaction (FSI) between flow and membrane is characterized via the continuous forcing IBM.The idea of this method is to add an external force term F to the right-hand side of the Navier-Stokes equation.The external force term originates from the membrane stress f membrane .A discretized delta function is used to spread the singular force on the membrane to the surrounding fluid and interpolate the fluid velocity back to the membrane.
δ is the three-dimensional Dirac-delta function, and x and x ′ are the locations in the flow domain T and on the cell surface S, respectively.A numerical approximation of the delta function is chosen to be where ∆ is the Eulerian grid size.

Direct Forcing IBM
Direct forcing IBM is used to treat another type of immersed boundary, the rigid, non-moving but geometrically complex vessel walls.Specifically, GNIBM is used to impose no-slip velocity boundary conditions on the vessel surface (53).The idea of this method is to modify these differential operators appropriately to satisfy the no-slip boundaries condition on the geometrically-complex vascular surface while maintaining second-order accuracy.

Indicator Function for Hematocrit Analysis
To characterize the cellular spatial distribution in RBC suspension, we define an indicator function I(x, t), such that the indicator function is one inside a cell, and zero outside a cell.It can be shown (69) that the indicator function I(x, t) follows a Poisson equation as where the G(x, t) is an Eulerian variable constructed from the cell surface normals n.

Wall Shear Stress Evaluation
To compute wall shear stress, the traction vector t = τ • n at the wall is determined using the velocity field approach outlined in (70,71), where τ is the stress tensor and n is the unit normal vector.The no-slip condition on the blood vessel surface yields the expression for the local wall shear stress t s = µ∂u s /∂r.Note s denotes the local stream-wise direction.The second-order differencing method is utilized to numerically evaluate velocity derivatives.

Model Verification and Validation
Elastic Spherical Capsule in Simple Shear Flow We consider a moving deformable spherical capsule subjected to simple shear flow, see Fig. S3.(A).
The deformation of a capsule will cause stretching stress and bending stress on its surface.The

Stationary Rigid Sphere in Simple Shear Flow
We consider a rigid stationary sphere subjected to a linear shear flow, which is given by u ∞ = [ γz, 0, 0].The schematic of the model set is shown in Fig. S4.(A).The analytical solution for this problem is given by where u, v, w are the velocity component in each dimension, r = x 2 + y 2 + z 2 is the distance from a point in flow to the center of the sphere, a is the radius of the sphere.Using the direct forcing IBM, we set the velocity on the sphere surface to zero (no-slip boundary condition).
Note the solution above is for unbounded shear flow, thus we impose the Dirichlet velocity boundary conditions on both the upper and bottom walls and choose a larger computation domain to reduce the error by periodic boundaries.After the velocity field reaches its steady state, our numerical results are then compared with the analytical solutions.L 1 and L 2 error norms for different velocity components u, v, w are plotted in Fig. S4.(B).It is found that the direct forcing IBM in our model presents second-order accuracy.
Fåhraeus-Lindqvist effect: Blood Relative Apparent Viscosity Fåhraeus and Lindqvist (73) made the noteworthy observation that the apparent viscosity at shear rates ≥ 100 s −1 , determined using Poiseuille's law in a capillary viscometer, exhibited a strong dependence on the capillary tube diameter.This phenomenon is called the Fåhraeus-Lindqvist effect.For large tube diameters, a constant viscosity plateau was evident.But within the range of 10 and 1000µm, the apparent viscosity decreases substantially with decreasing tube size, before increasing sharply for tubes smaller than 10µm.The increase for very small tubes is readily explained by the relative size of red blood cells, which are approximately 8µm in diameter, but the behavior at larger tube diameters is more subtle.In this study, we investigate the behavior of healthy RBCs flowing through tubes of varying radius, ranging from R = 6µm to R = 20µm, as illustrated in in Fig. S6., it is found that our simulation results are in good agreement with the experimental findings (76).Further results on the flow of diseased blood in bifurcations and junctions will be reported elsewhere.

Impact of Membrane Spontaneous Curvature on RBC Dynamics
Although efforts to understand RBC dynamics numerically have spanned the past two decades, the majority of these works have focused on RBC dynamics in cases where the RBC shape is symmetric across the shear plane or where the dimple is centered on the shear plane.Dupont experimentally observed rolling dynamics in a dextran solution where the viscosity ratio was less than unity.The discrepancy between simulation and experiment may result from the use of a spatially uniform spontaneous curvature that corresponds to a biconcave shape.It is important to note that to model the RBC membrane correctly, an assumption of the spontaneous shape has to be made, and finding the appropriate shape has been a challenge for both theoreticians and experimentalists.Sinha et al. (47) investigate the cell dynamics' dependence on the membrane's spontaneous curvature.They found that an oblate spheroidal spontaneous curvature maintains the dimple of the RBC during tank-treading dynamics and exhibits off-shear-plane, tumbling consistent with the experimental observations of Dupire et al. (45).For a complex structure such as an RBC membrane, it is possible that the natural shape for shear elasticity may differ from that for bending elasticity so the overall natural shape of an element results from the balance of bending and shear forces.
There have been endeavors to comprehend the impact of spontaneous shape on the ultimate dynamics of RBCs.Peng et al. ( 83) conducted a study on the influence of non-biconcave spontaneous shape on RBC dynamics and concluded that in order for an RBC to maintain its biconcave shape during tank-treading, as noted by Dupire et al. (45), the spontaneous curvature must be non-biconcave.In instances where a biconcave spontaneous curvature was employed, tank-treading could not be achieved without significantly perturbing the initial shape.Additionally, Cordosco et al. ( 84) explored non-biconcave spontaneous shapes and ascertained that the spontaneous shape has a significant impact on cell dynamics, depending on the viscosity ratio.They observed that the dimple in the RBC remained intact for both biconcave and oblate spontaneous shapes.However, it should be noted that in both works, Peng et al. (83) and Cordasco et al. (84), non-biconcave spontaneous curvatures were investigated under the imposition of spatially uniform spontaneous curvature, denoted as c 0 .It is worth highlighting that RBC membranes differ from model lipid bilayers in that they possess embedded proteins with an underlying spectrin cytoskeleton and an asymmetric bilayer leaflet composition, all of which modify c 0 , with proteins in particular, having been demonstrated to preferentially bind via curvature-sensing mechanisms.Hence, it can be argued that c 0 would be spatially In our work, to investigate the effects of membrane curvature on the RBC dynamics, two types of homogeneous normal RBCs suspension are simulated: one with spontaneous bending curvature being biconcave discoid, while another being oblate spheroid.Various vascular geometries are considered, including a slit (Fig. S7.), a straight cylindrical tube (Fig. S8.), and a curved (serpentine) channel (Fig. S9.).The simulation snapshots tell that RBC rest shape has a nontrivial impact on its near-wall dynamics: RBCs with biconcave discoid rest shape tend to "rolling", while RBCs with oblate spheroid rest shape perform "tank-treading", consistent with the findings in a prior numerical investigation by Sinha et al. (47).In spite of the change in orientational dynamics when the spontaneous shape is changed, no substantial change is observed in the number density distribution for normal RBCs.

Simulation snapshots of Cell Distributions in Straight Cylindrical Tube
The

Spherocytosis RBC suspension in Straight Cylindrical Tube
Movie S5.

SCD RBC suspension in Curved Channel
Movie S6.

IDA RBC suspension in Curved Channel
Movie S7.
deformability of a capsule is characterized by the nondimensional Capillary number Ca = µ γa/G and the nondimensional bending modulus κB = K B /a 2 G, where µ is the fluid viscosity, γ is shear rate, a is the radius of the capsule, G is the shear elasticity modulus and K B is the bending modulus.Larger Ca and κB means that a capsule is more flexible, while a stiffer capsule has smaller Ca and κB .The deformation of the capsule is described by the Taylor shape parameter defined as D xz = (L − B)/(L + B), where L and B are the maximum and minimum radial distances of an ellipsoid with the same inertia tensor.In Fig.S3.(B), we show that steady state values of Taylor deformation parameter D as a function of dimensionless Capillary number Ca for different κB .Good agreement is found between our numerical results and simulation results from previous literature(72).
Fig.S5.(A).Then the computed relative apparent viscosity, as a function of tube diameter, is shown in Fig.S5.(B).Fig.S5.indicates that our numerical simulation qualitatively predicts the Fåhraeus-Lindqvist effect.For comparison, we also show the empirical relation established by Pries et al (74) based on in vitro blood flow.It is found that our model predicts apparent viscosity in good agreement with the empirical relation at the physiological length scales we study in this work.Zweifach-Fung Effect: RBC Partitioning at a BifurcationThe partitioning of blood plasma and cells near a vascular bifurcation is complex.In the microcirculation, when RBCs pass through a bifurcating region of a blood vessel, they exhibit a tendency to preferentially flow into the daughter vessel with the higher flow rate, resulting in fewer cells flowing into the vessel with the lower flow rate.This phenomenon, termed the Zweifach-Fung effect (75), plays a crucial role in shaping blood flow distribution within the microvascular network.Here we compare simulation results from our model with experimental observations from previous literature(76).The radius of the inlet vessel is taken to be 10µm, in order to be close to the value used in experiments, as shown in Fig.S6..The Zweifach-Fung effect is quantified in Fig.S6.(B).The parameter η Q is defined to characterize partition as the ratio of the volumetric flow rate at one daughter branch to the volumetric flow rate at the parent vessel.Similarly, η N is defined as the ratio for cell number flow rate.From the results et al.(77) demonstrated that an elastic capsule with a prolate spheroid rest shape, whose axis of symmetry is oriented off the shear plane, will exhibit a unique final dynamical motion for all initial orientations.Depending on the capillary number, they observed three final dynamical states: (i) rolling for lower capillary numbers, (ii) wobbling in which the capsule processes around the vorticity axis as the capillary number is increased, and (iii) a swinging-oscillating motion in which the long axis of the capsule oscillates around the shear plane with decreasing amplitude of oscillation as the capillary number increases, resulting in an in-plane swinging motion at high capillary numbers.Wang  et al. (78) investigated the off-plane motion of oblate and prolate capsules and concluded that the final dynamical state could depend on the initial inclination angle.A recent study of RBCs in shear flow (79) has demonstrated that RBCs first tumble, then roll, transit to a rolling and tumbling stomatocyte, and finally attain polylobed shapes with increasing shear rate when the viscosity contrast between cytosol and blood plasma is large enough.Minetti et al. (80) give an exhaustive description of the dynamics under a shear flow of a large number of RBCs in a dilute regime is proposed.They identify which of the characteristic parameters of motion and of the transition thresholds depend on flow stress only or also on suspending fluid viscosity.Similarly, Cordosco and Bagchi (81) studied the off-plane motion of oblate, prolate, and biconcave capsules.Unlike Dupont et al. (77) and Wang et al.(78), they included membrane bending stiffness in their formulation and considered a spatially uniform spontaneous curvature in the case of biconcave capsules.They found that rolling was the dominant mode in the physiologically relevant viscosity ratio case (i.e., 5), tank-treading or wobbling mode at λ < 1, and an intermittent regime at low capillary numbers and low viscosity ratios, where the dynamics are dependent on the initial orientation.It is noteworthy that Bitbol (82) and Dupire et al.(77) inhomogeneous.Recently, using two different simulation techniques, Mauer et al.(85) construct a state diagram of RBC shapes and dynamics in shear flow as a function of shear rate and viscosity contrast(45,86).Their studies suggest that a nearly spherical stress-free shape best reproduces experimental results for the tumbling-to-tank-treading transition at low viscos-ity contrasts.Reichel et al.(87) combined simulation and experimental investigation of RBC shapes and dynamics in microchannels to provide a consistent RBC state diagram and illustrate the complexity of RBC behavior in the microflow.The RBC model employs a stress-free shape of the elastic spring network, corresponding to a spheroidal shape with a reduced volume of 0.96.Their simulation results agree well with experimental observations, allow the characterization of RBC variability in shear elasticity, and permit us to make a significant step toward quantitative measurements of RBC mechanical properties.
simulation snapshots showing the segregation phenomenon between normal and aberrant cells in the straight cylindrical tube are given in Fig. S10.. RBC-induced Local Wall Shear Stress Fluctuation on Curved Channel Surface Cross-sectional and center-plane normalized cell number density distribution of iron deficiency RBCs, sphero-echinocytes, and spherocytes over the curved channel in each corresponding binary suspension of normal RBCs with aberrant RBCs are shown in Fig.S11.. Aberrant RBCs margination provokes local wall shear stress fluctuation on curved channel surfaces are shown in Fig. S12., S13., and S14.. Movies Movie S1.SCD RBC suspension in Straight Cylindrical Tube Movie S2.IDA RBC suspension in Straight Cylindrical Tube Movie S3.COVID-19 RBC suspension in Straight Cylindrical Tube Movie S4.

Fig. S2 .Fig
Fig. S2.The surface mesh for capsules of (A) biconcave-discoid normal RBC, (B) sickle RBC, (C) iron deficiency RBC, (D) sphere-echinocyte, and (E) spherocyte.Note the size of each cell species in this figure is corresponding to the cell used in our simulation.The radius of normal RBC is 4µm; aberrant RBCs are smaller than normal RBCs.

Fig. S6 .
Fig. S6.(A) Simulation snapshots of healthy RBCs partitioning near a bifurcation.The radius of the inlet vessel is 10µm.(B) Comparison of the Zweifach-Fung effect in our simulation with experimental results from past literature (76).

Fig. S7 .
Fig. S7.(Top) (A) Simulation diagram of homogeneous normal RBC suspension between the slit under pressure-driven flow.(Re p = 0.1, Hematocrit = 0.15)(B) Steady-state wall-normal direction cell number density profile.Note here y/a = 0 denotes the slit center and y/a = 5 close to the wall.(Bottom) The side view of RBC suspension with bending spontaneous curvature being (C) oblate spheroid and (D) biconcave discoid.

Fig. S9 .Fig. S10 .
Fig. S9.(Top) Simulation snapshots for normal RBC suspension in the curved serpentine channel with spontaneous curvature being (A) oblate spheroid and (B) biconcave discoid.(Bottom) Center-plane and cross-sectional cell number density distribution for RBC suspension with spontaneous curvature being (C) oblate spheroid and (D) biconcave discoid.