Computational Biorheology of Human Blood Flow in Health and Disease

Hematologic disorders arising from infectious diseases, hereditary factors and environmental influences can lead to, and can be influenced by, significant changes in the shape, mechanical and physical properties of red blood cells (RBCs), and the biorheology of blood flow. Hence, modeling of hematologic disorders should take into account the multiphase nature of blood flow, especially in arterioles and capillaries. We present here an overview of a general computational framework based on dissipative particle dynamics (DPD) which has broad applicability in cell biophysics with implications for diagnostics, therapeutics and drug efficacy assessments for a wide variety of human diseases. This computational approach, validated by independent experimental results, is capable of modeling the biorheology of whole blood and its individual components during blood flow so as to investigate cell mechanistic processes in health and disease. DPD is a Lagrangian method that can be derived from systematic coarse-graining of molecular dynamics but can scale efficiently up to arterioles and can also be used to model RBCs down to the spectrin level. We start from experimental measurements of a single RBC to extract the relevant biophysical parameters, using single-cell measurements involving such methods as optical tweezers, atomic force microscopy and micropipette aspiration, and cell-population experiments involving microfluidic devices. We then use these validated RBC models to predict the biorheological behavior of whole blood in healthy or pathological states, and compare the simulations with experimental results involving apparent viscosity and other relevant parameters. While the approach discussed here is sufficiently general to address a broad spectrum of hematologic disorders including certain types of cancer, this paper specifically deals with results obtained using this computational framework for blood flow in malaria and sickle cell anemia.


INTRODUCTION
Parasitic infections or genetic factors can drastically change the viscoelastic properties and the biconcave (discocyte) shape of red blood cells (RBCs). 36 For example, the parasite Plasmodium falciparum that invades the RBCs (Pf-RBCs) of most malaria patients markedly affects the RBC membrane properties resulting in up to a ten-fold increase of its shear modulus and a spherical shape in the later stages of the intra-cell parasite development. 36 Sickle cell anemia is another blood disorder caused by the polymerization of the hemoglobin inside the RBCs, which, in turn, leads to dramatic changes in their shape and deformability. These changes combined with the increased internal viscosity affects the flow of sickled RBCs through the post-capillary venules leading to flow occlusion. 36,77 Other hereditary diseases with similar effects are spherocytosis and elliptocytosis. 14 In the former, RBCs become spherical with reduced diameter and carry much more hemoglobin than their healthy counterparts. In the latter, RBCs are elliptical or oval in shape and exhibit reduced deformability.
These hematologic disorders, despite their differing origins as infectious diseases arising from external vectors or as hereditary abnormalities ascribed to genetic defects, also reveal some common characteristics in terms of the remodeling of cytoskeleton. Such molecular remodeling of the spectrin cytoskeleton induces a change in the structure and viscoelastic properties of individual RBCs. Therefore, studying the rheological properties of blood and its components such as the RBC can aid greatly in the understanding of many major diseases. To this end, new advanced experimental tools are very valuable in obtaining the biophysical properties of single RBCs in health and disease, which are required in formulating multiscale methods for modeling blood flow in vitro and in vivo. Specifically, advances in experimental techniques now allow measurements down to the nanometer scale, and include micropipette aspiration, 39,159 RBC deformation by optical tweezers, 21,76,147 optical magnetic twisting cytometry, 131 three-dimensional measurement of membrane thermal fluctuations, 118,124 and observations of RBCs immersed in both shear and in pressuredriven flows. 2,66,135,148,154 Micropipette aspiration and optical tweezers techniques tend to deform the entire RBC membrane directly, while optical magnetic twisting cytometry and measurements of membrane thermal fluctuations probe the membrane response locally.
Several numerical models have been developed recently including a continuum description 46,52,69,125 and a discrete approximation on the spectrin molecular level 38,96 as well as on the mesoscopic scale. 42,44,115,122 Some of the models suffer from the assumption of a purely elastic membrane, and are able to capture only the RBC mechanical response, but cannot quantitatively represent realistic RBC rheology and dynamics. Fully continuum (fluid and solid) models often suffer from non-trivial coupling between nonlinear solid deformations and fluid flow with consequential computational expense. Semi-continuum models 46,125 of deformable particles, which use immersed boundary or front-tracking techniques, can be quite efficient. In these models, a membrane is represented by a set of points which are tracked in Lagrangian description and are coupled to an Eulerian discretization of fluid domain. However, these models employ the same external and internal fluids and do not take into account the existing viscosity contrast between them. In addition, continuum models omit some mesoscopic and microscopic scale phenomena such as membrane thermal fluctuations, which affect RBC rheology and dynamics. 114 On the microscopic scale, detailed spectrin molecular models of RBCs are limited by the demanding computational expense. Therefore, we will focus here on an accurate mesoscopic modeling of RBCs.
There exist a few mesoscopic methods 42,44,115,122 for modeling deformable particles such as RBCs. Dzwinel et al. 44 model RBCs as a volume of elastic material having an inner skeleton. This model does not take into account the main structural features of the RBC, namely a membrane filled with a fluid, and therefore it cannot accurately capture the dynamics of RBCs, such as the observed tumbling and tank-treading behavior in shear flow. 2,142 Three other methods 42,115,122 employ a conceptually similar approach to the method we will present here, where the RBC is represented by a network of nonlinear viscoelastic springs in combination with bending rigidity and constraints for surface-area and volume conservation. Dupin et al. 42 coupled the discrete RBC to a fluid described by the Lattice Boltzmann method 145 obtaining promising results. However, this model does not consider external and internal fluids separation, membrane viscosity, and thermal fluctuations. Noguchi and Gompper 115 employed Multiparticle Collision Dynamics 102 and presented encouraging results on vesicles and RBCs; however, they do not use realistic RBC properties and probe only very limited aspects of RBC dynamics. Pivkin and Karniadakis 122 used dissipative particle dynamics (DPD) 79 for a multiscale RBC model which will be the basis of the general multiscale RBC (MS-RBC) model we will present here, following the earlier versions developed in Fedosov et al. 56 The MS-RBC model described here is able to successfully capture RBC mechanics, rheology, and dynamics. Potential membrane hardening or softening as well as the effects of metabolic activity can also be incorporated into the model leading to predictive capabilities of the progression of diseases such as malaria. Theoretical analysis of the hexagonal network yields its mechanical properties, and eliminates the need for ad hoc adjustments of the model parameters. Such models can be used to represent seamlessly the RBC membrane, cytoskeleton, cytosol, the surrounding plasma and even the parasite. This paper is organized as follows: In ''Materials and Methods'' section, we review the basic DPD theory and the MS-RBC models. In ''Healthy Blood Flow'' section, we present rheology results of healthy blood flow in capillaries and arterioles, and comparisons with available experimental observations. In ''Diseased Blood Flow'' section, we review recent results on modeling blood flow in malaria and in sickle cell anemia. We conclude in ''Discussion'' section with a brief summary and a discussion on the potential of multiscale modeling in predicting the onset and progression of other hematologic disorders.

Fluid Flow Modeling
Fluid flow modeling is referred here to the modeling of the Newtonian solvent flow, which mimics blood plasma. In particle-based methods a fluid is represented by a collection of interacting particles, which recovers hydrodynamics on the length scales several times larger than the particle size. Examples include molecular dynamics, 6 DPD, 51,75,79 multi-particle collision dynamics, 74,102 and smoothed particle hydrodynamics (SPH). 99,110 The DPD system consists of N point particles, which interact through three pairwise forces-conservative (C), dissipative (D), and random (R)-such that the force on particle i is given by where the sum runs over all neighboring particles j within a certain cutoff radius. The conservative force mainly controls fluid compressibility, while the dissipative force governs fluid viscosity. In addition, the pair of dissipative and random forces forms a thermostat such that an equilibrium temperature can be imposed. 51 Recently, the smoothed dissipative particle dynamics (SDPD) method 50 has been proposed, which has several advantages over the conventional DPD. The SDPD method is derived from SPH in combination with a thermostat similar to that used in DPD. In comparison to DPD, the SDPD method allows one to directly input fluid transport coefficients (e.g., viscosity) and to select an arbitrary equation of state, and therefore to have full control over fluid compressibility. Other methods for fluid flow modeling include continuum approaches such as the Navier-Stokes equation or its modifications. 161 In addition, the other two methods which are quite popular for fluid flow modeling are the lattice Boltzmann method (LBM) method 145 and Brownian dynamics. 49

Modeling Blood Cells
Blood cells are modeled by a flexible membrane with constant area and volume. As an example, the RBC membrane is represented by a network of springs, which corresponds to a triangulation on the membrane surface. 34,38,43,56,57,96,113,122 The free energy for each cell is given by where V s is the spring's potential energy, V b is the bending energy, and V a+v corresponds to the area and volume conservation constraints. The V s contribution supplies proper cell membrane elasticity. A ''dashpot'' is attached to each spring, and therefore, the spring forces are a combination of conservative elastic forces and dissipative forces, which provide membrane viscous response similar to a membrane viscosity. The bending energy mimics the bending resistance of the cell membrane, while the area and volume conservation constraints mimic area-incompressibility of the lipid bilayer and incompressibility of a cytosol, respectively. More details on the cell model can be found in Fedosov et al. 56,57 Linear analysis for a regular hexagonal network allows us to relate the model parameters uniquely with the network macroscopic properties, see Fedosov et al. 56,57 for details. Thus, in practice, the given macroscopic cell properties serve as an input to be used to calculate the necessary mesoscopic model parameters without any manual adjustment. We also employ a ''stress-free'' model, 56,57 which eliminates existing artifacts of irregular triangulation. It is obtained by computational annealing such that each spring assumes its own equilibrium spring length adjusted to be the edge length after triangulation. Both internal and external fluids are simulated by a collection of free DPD particles and are separated by the cell membrane through bounce-back reflections at the membrane surface. Moreover, a dissipative force between fluid particles and membrane vertices is set properly to account for the no-slip boundary conditions at the membrane surface. More details on boundary conditions can be found in Fedosov et al. 56,57 The other class of methods for cell membrane modeling corresponds to a continuum representation. Thus, cells are modeled using some type of a constitutive equation from solid mechanics. Examples here include boundary-integral formulation 125,165 and finite element discretization. 100 More details on RBC and blood flow modeling can be found in recent reviews. 62,85,97 HEALTHY BLOOD FLOW Blood flow in microcirculation 123 affects a number of vital processes in the human body including oxygen delivery to tissues, removal of waste products, immune response, haemostasis, etc. Detailed experimental measurements of blood flow in vivo down to a cell level are very difficult and often not feasible. This has motivated a large worldwide scientific effort to derive reliable numerical models of blood flow in order to study blood hydrodynamics in vessel networks, its rheology, and the effect on various blood-affected processes. Here, we will mainly focus on mesoscopic modeling of blood flow down to the single blood cell level; see also a recent review. 62 We also present typical average plasma and RBC properties in Table 1.

Blood Flow in Tubes
Blood flow in tubes or rectangular channels is a relatively simple problem which mimics the main characteristics of blood flow in a in vivo vessel. It allows the quantification of blood flow resistance for various conditions by estimating such factors as the hematocrit, flow rate, and RBC aggregation. This model also permits quantitative simulation of such processes as the migration, deformation, and dynamics of single cells within the flow. As an example, Fig. 1 shows a snapshot of flow of blood in a tube with diameter D = 20 lm.

Fahraeus Effect
Fahraeus performed in vitro experiments on blood flowing through glass tubes 53 and found an elevated value of discharge hematocrit (H d ), defined as the RBC volume fraction exiting a tube per unit time, in comparison with the RBC concentration within a tube called tube hematocrit (H t ). This finding, which is called the Fahraeus effect, is a consequence of the cross-stream migration of RBCs in tube flow. RBCs in blood flow migrate away from the vessel walls towards the centerline 28,72 due to hydrodynamic interactions with the walls and their deformability 4,19,40,108 ; such interactions are also referred to as lift force. Crowding of RBCs around the tube center results in their faster motion along the tube with respect to the average blood flow velocity, thus leading to an increased H d value measured experimentally.
The Fahraeus effect is one of the standard tests for any blood flow model, and has been modeled in simulations of blood flow in rectangular channels 43 and in tubes. 59,93 In contrast to experiments, 53,128 where H d is the main controlled parameter and H t has to be measured, in simulations H t is an input and corresponds to the volume fraction of RBCs placed within the channel, while H d is computed. In order to compute the ratio H t /H d in simulations a simple mass balance can be employed in flow, which yields where " v is the average flow velocity and v RBC is the average velocity of red cells. Figure 2a shows the ratio H t /H d for two hematocrits and various tube diameters D. The simulation results 93 are compared with the empirical fits to experimental data on blood flow in glass tubes. 128 A strong effect of RBC migration to the tube center on the discharge hematocrit is observed for tubes with diameters of up to about 200 lm. For larger tube diameters, the layer next to the wall void of RBCs, called the cell-free layer (CFL), becomes rather small in comparison with D, and therefore the Fahraeus effect is negligible.

Fahraeus-Lindqvist Effect
Another phenomenon directly related to RBC migration towards a vessel centerline and the formation of a CFL next to the wall is the Fahraeus-Lindqvist effect, 55 describing a decrease in blood flow resistance with decreasing tube diameter. 128 The flow resistance is affected by the CFL, which serves as an effective lubrication layer for relatively viscous flow at the core consisting primarily of RBCs. For large tube diameters the CFL thickness is small with respect to D, and therefore the bulk blood viscosity is essentially measured. As the tube diameter is reduced, the CFL thickness becomes more significant with respect to D leading to an effective decrease in flow resistance. The blood flow resistance in glass tubes is quantified by an effective viscosity, termed the apparent blood viscosity, using the Poiseuille law in tubes where DP=L is the pressure drop along the tube length L and Q is the volumetric flow rate. For convenience, experimental results 128 are presented in terms of the relative apparent viscosity defined as where g 0 is the blood plasma viscosity. Thus, g rel characterizes the effective viscosity of blood with respect to that of plasma. Similar to the Fahraeus effect, the Fahraeus-Lindqvist effect also serves as a standard validation test for any blood flow model. Figure 2b shows the relative apparent viscosity for different H t values and tube diameters, where the solid curves are the fits to experimental measurements 128 ; the symbols correspond  to simulation results. 93 The minimum value of g rel is found for the tube diameter of D m = 7-8 lm. In tubes with diameters smaller than D m , RBCs have to strongly deform and squeeze through the tubes attaining bullet-like shapes, 3,140,155 which leads to a steep increase in the relative apparent viscosity. On the other hand, when D > D m , the thickness of the CFL in comparison with the tube diameter becomes smaller, leading to an effective viscosity increase. The Fahraeus-Lindqvist effect has been quantitatively captured in 3D blood flow models in cylindrical vessels 5,59,68,93 as well as in rectangular channels. 41,43 Furthermore, this effect has also been quantitatively predicted by 2D blood flow models 9 and by the simplified blood models, 84,105,116 which do not properly capture RBC membrane deformability. This clearly indicates that the deformability and dynamics of individual cells do not significantly affect the flow resistance. Thus, it appears computationally advantageous to divide blood flowing in tubes into two regions: (i) a relatively viscous flow core consisting mainly of cells, and (ii) the CFL region devoid of RBCs. This assumption leads to the so-called twophase model 28,71,93,138 for blood flow in tubes or vessels, where the flow is divided into core and near-wall regions having different viscosities. Although such models describe well the flow resistance in straight tubes, they might not provide realistic simulations in complex geometries (e.g., vessel networks) due to spatial variability of the CFL 88 and non-trivial partitioning of RBCs within the network. 130 Finally, the blood flow resistance depends on flow rate and aggregation interactions between RBCs. Experimental data compiled in Pries et al. 128 and simulations of blood flow in small tubes 68 suggest that the dependence of flow resistance on flow rate disappears at sufficiently high flow rates exceeding the pseudo-shear rate _ c ¼ " v=D of approximately 50 s 21 . In fact, fits to experimental data presented in Fig. 2 are based solely on experimental studies, where the pseudo-shear rates are higher than 50 s 21 . 128 As the flow rate is decreased such that _ c < 50 s À1 ; the blood flow resistance increases for RBC suspensions 68,134 without aggregation interactions between cells. RBC aggregation in whole blood occurs in equilibrium or at sufficiently low shear rates. Viscometric measurements of whole blood viscosity 25,106,141 have shown a steep viscosity increase for shear rates below 10 s 21 due to RBC aggregation in comparison with analogous nonaggregating RBC suspension. Thus, the flow resistance of blood in tubes is not affected by RBC aggregation for sufficiently high flow rates, and this has been confirmed experimentally 28,134 and in simulations. 93 At sufficiently low flow rates, experimental studies 28,54,111,134 report a decrease in the relative apparent viscosity due to RBC aggregation, but this has not been yet explored in simulations.

Cell-Free Layer
The CFL is a layer of blood plasma near the wall that is devoid of RBCs due to their migration to the vessel center. The viscosity in the CFL region is equal to that of blood plasma, while the viscosity in the tubecore region populated with RBCs is several times larger. Hence, the formation of CFL leads to a better efficiency for the RBC core to flow resulting in the aforementioned Fahraeus and Fahraeus-Lindqvist effects. The thickness of CFL d may serve as an alternative indicator for blood flow resistance. In in vitro 17,134 and in vivo 88,101,163 experiments, the outer edge of the RBC core has been tracked to calculate its average position and deduce the CFL thickness. A similar approach has also been employed in The variability in CFL thicknesses also exists across various CFL experimental measurements, 88,101,163 so the agreement of the simulation results with in vivo data is expected to be qualitative. Several factors may contribute to the variability of in vivo measurements of CFL and the discrepancies with simulations including the glycocalyx layer at vessel walls, variations in vessel diameter and length, close proximity of the site of CFL measurements to vessel bifurcations, vessel elasticity, and spatial resolution of the recorded data. 47,59,88,127 An alternative means to compute the thickness of CFL in simulations is to analyze the local hematocrit distribution (see Fig. 4a), which drops quickly as we exit the RBC core region and enter the CFL moving along the tube radial direction. 9,41,68 The simulation of non-aggregating RBC suspension in tube flow 68 showed that d increases as the pseudo-shear rate is increased from 0 to approximately 50 s 21 . This result indicates that the flow resistance is decreasing in agreement with the experimental measurements in Reinke et al. 134 For flow rates above 50 s 21 , the CFL thickness remains nearly constant, 68 which supports the argument that the flow resistance becomes independent of the flow rate as discussed in ''Fahraeus-Lindqvist Effect'' section. In contrast to a constant d for high flow rates, in vivo experiments 88 and 3D simulations 59 have shown a mild decrease in the CFL thickness with increasing flow rate. Finally, the effect of RBC aggregation, which is important only at low shear rates, leads to an increase in the CFL thickness as the flow rate is reduced due to the possibility to form a more compact RBC core of the flow. 89,134 Consequently, RBC aggregation results in a lower flow resistance as the flow rate is decreased. Note that RBC aggregability reverses the trend of the flow resistance for low flow rates in comparison with a non-aggregating suspension; the RBC aggregation effect has not been yet investigated by simulations.

RBC Distribution and Deformation in Flow
Migration of RBCs in blood flow leads to a variation in local RBC density and hematocrit. Departures from normal blood flow characteristics are potential indicators of disease pathology. For tube flow, the local hematocrit H(r) is defined as the volume fraction of RBCs within an annular region with the radius r and thickness Dr: Figure 4a presents the local hematocrit distributions for various H t values and tube diameters. For small tube diameters, H(r) in the tube center is significantly larger than H t in agreement with the simulation results in Freund and Orescanin 68 Thus, RBCs are mainly located at the tube center, however their dynamics and structure depend on H t . 93 At low H t , RBCs attain a parachute shape and move in a trainlike arrangement, which has been also reported in McWhirter et al. 104 In contrast, at high H t , RBCs arrange into the so-called zig-zag mode found experimentally in Gaehtgens et al. 70 and in simulations. 104 As the tube diameter is increased, the local hematocrit converges towards the H t value. However, for low H t a peak in H(r) may still exist around the tube center due to close packing of RBCs within the centerline region where local shear rates are small. For tube diameters DJ200 lm a nearly homogeneous distribution of RBCs is expected, and therefore a continuum approximation for blood flow in large tube is appropriate. 93 Quantitative analysis of local RBC deformation in tube flow can be conveniently performed using the RBC gyration tensor. 103 Eigenvalues of the gyration tensor characterize cell shape along the directions defined by the tensor eigenvectors. The distributions of the largest eigenvalue k max for the RBC across the tube are shown in Fig. 4b indicating that RBCs are slightly compressed in the tube center and get stretched more and more by the flow as we approach the CFL region. The k max distributions can be qualitatively divided into the three regions: i) a centerline region, where k max is nearly constant; ii) the region between the tube center and CFL, where k max slowly increases; iii) a region next to the CFL, where k max steeply increases. Similar conclusions for RBC deformation in tube flow have also been reported in Alizadehrad et al. 5 Finally, for sufficiently large D the function of k max (r) converges to a common curve, and therefore the dependence on tube diameter disappears.

Blood Rheology
A major aim of computational cell biorheology is to predict the macroscopic flow properties of a suspension (e.g., shear viscosity, yield stress) from the mesoscopic properties of the constituent particles (e.g., size, deformability, inter-particle interactions). Bulk blood properties under shear have been measured in a number of laboratories 25,106,121,137,141,150 predicting a non-Newtonian shear thinning behavior. Here, we draw attention to two different sets of viscometric measurements. The first set corresponds to whole blood, where no significant changes are introduced to freshly drawn blood except for a necessary initial stabilization with an anti-coagulant. Such procedure for sample preparation leaves the aggregation properties among RBCs virtually unaffected, and therefore cell aggregates, known as rouleaux, can be observed. 23,24,106,107 The second set of viscometric measurements corresponds to the so-called non-aggregating RBC suspensions, where RBCs are washed and re-suspended into a neutral solvent. RBC aggregation interactions strongly depend on the nature and concentration of available proteins or polymers 24,107 and can be triggered by adding them to non-aggregating RBC suspensions. 107,126,144 Whole blood is also known to exhibit a yield stress 27,29,106 defined as a threshold stress that must be exceeded for flow to begin.

Reversible Rouleaux Formation
Aggregation interactions among RBCs lead to the formation of rouleaux structures in equilibrium or in relatively weak flows. A rouleaux resemble stacks of coins as illustrated in Fig. 5 showing a simulation snapshot in equilibrium. There exist two different hypotheses for the RBC aggregation in polymer or protein solutions: (i) the bridging model, which suggests that cross-linking between cells may be achieved by polymer adsorption or adhesion on a cell surface, 16,22 and (ii) the depletion model, which suggests that RBC aggregation arises from depletion interactions. 13,112 In equilibrium, rouleaux formation consisting of a few RBCs is first observed, followed by further coalescence into branched rouleaux networks. 106,136,137 In flow, the large rouleaux break up into smaller ones. However, at sufficiently high shear rates, rouleaux structures no longer exist. 168 The aggregation process is reversible once the shear rate is decreased.
RBC aggregation has been simulated in 2D, 11,158,164 where aggregation and disaggregation of several cells in shear flow has been studied. The effect of RBC aggregation on blood rheology has also been investigated in 3D simulations. 98 However, the marked viscosity increase at low shear rates has not been captured due to the very small size of a simulated system with only six RBCs. Finally, recent 3D simulations 63  quantitative predictions for the strength of RBC aggregation and its effect on blood viscosity.

Blood Viscosity and Yield Stress
As noted earlier, the viscosity of whole blood is strikingly different from that of non-aggregating RBC suspension due to aggregation interactions between RBCs. Figure 6 presents the relative viscosity of blood, defined as the ratio of RBC suspension viscosity to the solvent viscosity, vs. shear rate at H t = 0.45. The blood viscosity was calculated in shear flow 63 using the MS-RBC model. 56 The model results compare very well with the experimental viscometric measurements. 25,106,141 The correct predictions of viscosity by the model facilitates calculation of the strength of aggregation forces between two RBCs, which appears to be in the range of 3.0 pN to 7 pN. Note that these forces are much smaller than those used in single RBC stretching tests with optical tweezers. 147 Whole blood is considered a fluid with a yield stress (a threshold stress for flow to begin), 27,29,106,121 which is usually calculated by extrapolation of available viscometric data to the zero shear rate. For blood the Casson's equation 20 is often used where s xy is the shear stress, s y is a yield stress and g is the suspension viscosity at large _ c: The Casson's relation has been used to compute yield stress for pigment-oil suspensions, 20 Chinese ovary hamster cell suspensions, 81 and blood. 107 Figure 7 shows the simulation results from Fedosov et al. 63 where computed shear stresses are extrapolated to the zero shear rate using a polynomial fit in Casson coordinates. The extrapolated s y values are virtually zero for non-aggregating RBC suspensions, while the presence of RBC aggregation leads to a nonzero yield stress. The simulation results support the hypothesis that whole blood has a non-zero yield stress due to aggregation interactions between RBCs. 27,29,106 Moreover, the simulated values of s y lie within the range of 0.0015 Pa and 0.005 Pa for yield stress measured experimentally. 106 Another aspect which has been discussed in Fedosov et al., 63 is a link between RBC suspension macroscopic properties and its microscopic characteristics such as cell structure, deformation, and dynamics. The analysis of the suspension's structure has confirmed that a steep rise in viscosity for the aggregating RBC suspension at low shear rates is due to the existence of small aggregates consisting of 2-4 RBCs. As the shear rate is increased, the suspension's structure completely disappears, and therefore the viscosity at high shear rates becomes independent of RBC aggregation. Furthermore, the shear-thinning property of a non-aggregating RBC suspension is strongly associated with the transition of RBCs from tumbling to the tank-treading dynamics 2,142 and the alignment of tank-treading RBCs with the shear flow. Simulations of a single elastic oblate-shape capsule in shear flow, 10 which mimic a dilute solution, have also shown a shear thinning behavior which mainly occurs within the capsule tank-treading regime. Another simulation study of dense near-spherical capsules 26 have also shown a slight shear-thinning behavior due to the deformation of capsules in shear flow. In addition, normal stress differences have been computed for the capsule suspensions. Numerical studies of the rheology of dense suspensions of deformable particles and cells are still rather limited. Such simulations can be used to investigate a wide class of complex fluids (e.g., cell, vesicle, and capsule suspensions) and to tune their properties to modulate cell behavior. In addition, such 3D highfidelity simulations will allow one to study in detail the connection between macroscopic and microscopic properties of such suspensions.

Margination of White Blood Cells and Platelets
Margination of different solutes (e.g., white blood cells, platelets, drug carriers) in blood flow is the migration process of the solutes towards vessel walls. The margination process is essential for many blood solutes in order to perform their function, as they come into contact with vessel walls when necessary. For example, white blood cells (WBC) have to marginate towards the walls 12,65,73 in order to be able to efficiently adhere to vascular endothelium and eventually transmigrate into the surrounding tissues. 143 The margination mechanism is mediated by RBCs, which migrate towards the vessel center 28,72 due to the hydrodynamic lift force, 4,19,40,108 and effectively push out various solutes from the central region into the CFL. More precisely, the margination mechanism is a consequence of the competition between lift forces on RBCs and blood solutes, and their interactions in flow. Recently, a margination theory for binary suspensions 90,91 has been proposed through two mechanisms: (i) wall-induced particle migration due to lift force, and (ii) particle displacement due to pair collisions of different solutes or their interactions in flow; the latter mechanism is also often referred to as shear-induced diffusion.

White Blood Cell Margination
WBC margination has been studied in a number of experiments 1,12,65,73,83,119 and simulations. 61,67,146 Several blood and flow properties contribute to WBC margination including H t , flow rate, RBC aggregation, and vessel geometry. In vivo experiments on mesenteric venules of rat 65 showed an increase of WBC adhesion rate (and consequently margination) for high H t J0:45 as the flow rate was decreased. A recent microfluidics study 83 concluded that the most efficient WBC margination occurs within intermediate H t . 0.2-0.3, while at lower or higher H t values WBC margination is attenuated. In contrast, in vitro experiments in glass tubes 1 and 2D simulations 67 have reported no significant sensitivity of WBC margination on blood hematocrit. The effect of flow rate on WBC margination is consistent across various studies, 61,65,67,83,119 which show strong WBC margination at low flow rates, characteristic for venular blood flow. Furthermore, RBC aggregation leads to an enhancement of WBC margination, 1,61,65,83,119 while the effect of complex geometries (e.g., vessel contraction, expansion, and bifurcation) on WBC margination has a very limited interpretation so far. The 2D simulation study 146 and microfluidic experiments 83 mimicking a vessel expansion suggest that the WBC margination is enhanced in post-capillary venular expansions. Finally, a combined in vitro and in vivo study of blood flow around bifurcations 156 have shown preferential adhesion of WBCs near bifurcations.
Recent 2D simulations 61 provided a systematic study of WBC margination for various conditions including flow rate, H t , WBC deformability, and RBC aggregation, and attempted to reconcile existing contradicting observations. Using blood flow simulations, WBC center-of-mass distributions were computed for various flow rates and H t values with an example shown in Fig. 8. Clearly, the strongest WBC margination is found within a range of intermediate H t = 0.25-0.35, while at lower and higher H t values WBC margination is attenuated. The weak WBC margination at low H t has been expected due to low concentration of RBCs. Surprisingly, however, no significant WBC margination at high H t was observed. The mechanism of WBC margination attenuation at high H t appeared to be also related to the presence of RBCs in flow. Thus, at low enough H t the region in front of a marginated WBC in blood flow remains virtually void of RBCs as they pass above the slowly moving WBC. As we increase H t to a certain value, RBCs may often enter that region due to high cell crowding, and effectively lift off the WBC away from the walls. 61 This lift-off mechanism is different from the lift force due to cell-wall hydrodynamic interactions and is governed by the particulate nature of blood.
Having calculated WBC distributions, one can define a WBC margination probability by integrating the probability distributions, for instance, up to a distance of 1.1R WBC away from the walls, where R WBC is the WBC radius. Figure 9 presents a WBC margination diagram for a wide range of flow rates and hematocrits. It shows that efficient WBC margination occurs only within intermediate ranges of flow rates and H t values. These results are consistent with experimental observations, which report WBC adhesion mainly in venular (not arteriolar) part of microcirculation, since the characteristic values of _ c Ã in venules of a comparable diameter are in the range _ c Ã ¼ 1À25; while in arterioles _ c Ã J30. 61,123,129 The simulation results also agree with in vitro experiments on WBC margination, 83 which identified optimal WBC margination in the range H t = 0.2-0.3. The discrepancies with previous simulations 67 and experiments, 1 which found WBC margination and adhesion to be independent of H t , can also be rationalized by noting that the studied flow rates and H t values in Abbitt and Nash 1 and Freund 67 fell almost entirely into the region of strong WBC margination.
The effect of WBC deformability has also been explored in Fedosov et al., 61 where the region of strong WBC margination shrinks substantially as a WBC becomes more deformable. This is due to an increase of the lift force on a deformed WBC, since it may significantly depart from a spherical shape, which leads to a reduction in its margination. RBC aggregation has also been shown to enhance WBC margination. 61 In particular, no significant differences in WBC margination between aggregating and non-aggregating RBC suspensions have been found for high flow rates similar to the effect of RBC aggregation on blood viscosity discussed in ''Blood Rheology'' section. However, WBC margination becomes more pronounced at high H t values due to RBC aggregation. Here, RBC aggregation interactions lead to a less dispersed flow core, and therefore the shedding of RBCs into the region in front of a marginated RBC is much reduced. Future studies of WBC margination will likely focus on 3D realistic models, and WBC margination and adhesion in complex microvascular networks.
Similar margination mechanisms are of importance in a number of diseases. For instance, margination of circulating tumor cells, which affects their adhesion to vascular endothelium, is expected to resemble that of WBCs due to the similarities in cell shapes and sizes. The WBC margination results are also relevant in spherocytosis disorder and in malaria, where spherically-shaped cells are expected to be mainly present near vessel walls resulting in an increased blood flow resistance and an enhanced cell adhesion when possible. Furthermore, the margination effect directly applies to segregation of various cells and particles in flow due to the differences in their deformability and shape. Recent numerical and theoretical investigations 90,91 have shown that less deformable cells or capsules are preferably segregated toward channel walls. Such an effect is expected to play a role in malaria and sickle-cell anemia, where margination of stiff cells (e.g., Pf-RBCs) contributes to their adhesive potential. On the other hand, margination and segregation properties of different cells in blood flow can be exploited in diagnostics and treatment which might require rare cell detection and separation. 15,78,80 Finally, solutes smaller than WBCs (e.g., platelets, large globular proteins such as von Willebrand factor) are also subject to margination in blood flow, which is relevant for margination (and therefore adhesion) of micro-and nano-carriers used for drug delivery. Next, the margination effect for platelets is discussed.

Margination of Platelets
Preferential migration of platelets towards the vessel walls in blood flow is critical for haemostasis and thrombosis, since marginated platelets can better respond to injuries at the vessel wall. Margination of platelets as well as of WBCs is strongly affected by local hematocrit, flow rate, vessel geometry, and RBC aggregation properties. These characteristics might be significantly altered in different diseases in comparison with a healthy state leading to irreversible changes in cell margination properties which might affect their proper function. An increased concentration of platelets near the walls has been confirmed in in vitro 45,151,157 and in vivo 149,162 experiments. Experiments in glass channels with a suspension of RBCs and latex beads mimicking platelets 45,151 for different H t values have shown strong bead margination at high H t J0:3 values. The experiments for various flow rates 45,151 demonstrated a non-trivial bead margination dependence. At low shear rates platelet margination was weak and it increased with increasing flow rates, while the margination decreased again at very high shear rates; this behavior is qualitatively similar to that of WBCs. In vivo experiments 162 have found that platelet margination is different in arterioles and venules confirming the effect of shear rate on margination. In view of these observations, it is plausible to expect that the mechanism of platelet margination is similar to that for WBCs described above, however its detailed understanding is still an open question.
There exist a number of numerical studies of platelet margination and its mechanisms. 7,32,33,133,152,153,166,167 The 2D simulations of blood flow using ellipsoid-like RBCs and circular platelets 7 have shown an increased platelet margination with increasing H t as well as with an increase of shear rate in agreement with the experimental observations. The other set of 2D simulations 32,33 has led to similar conclusions for the effects of H t and shear rates on platelet margination. The local drift and diffusion of platelets have also been measured 32,33 in an attempt to describe the platelet margination process using the continuum drift-diffusion equation. The drift-diffusion equation results in a good qualitative description of platelet margination, where the drift can be hypothesized to arise from wallplatelet hydrodynamic interactions (i.e., lift force), while the shear-induced diffusivity of platelets is due to their collisions (or interactions) with RBCs in blood flow. Recent 3D simulations on platelet margination 166,167 have found the margination process to be diffusional, which fits well into the drift-diffusion mechanism. Another 3D simulation study 133 have confirmed the effect of H t on platelet margination, and, in addition, reported that spherical platelets marginate better than those having an ellipsoidal shape. Although the drift-diffusion mechanism for platelet margination seems to provide plausible explanation, it is still not clear whether it will lead to a universal description of the margination process. For instance, a recent theoretical model 152 for platelet adhesion in blood flow suggests a rebounding collision mechanism for platelets due to their interactions with the flow core consisting primarily of RBCs.

Rheology in Malaria
Red blood cells parasitized by the malaria-inducing Plasmodium falciparum (Pf) parasites, referred to here as Pf-RBCs, are subject to changes in their mechanical and rheological properties as well as in their morphology 30,139,147 during intra-erythrocytic parasite development, which includes three intra-erythrocytic stages that span the 48-h asexual cycle: from ring fi trophozoite fi schizont stages. Gradual progression through these stages leads to considerable stiffening of Pf-RBCs as found in optical tweezers stretching experiments 35,109,147 and in diffraction phase microscopy by monitoring the membrane fluctuations. 37,118 Pf development also results in vacuoles formed inside of RBCs possibly changing the cell volume. Thus, Pf-RBCs at the final stage (schizont) often show a ''near spherical'' shape, while in the preceding stages maintain their biconcavity. These changes greatly affect the rheological properties and the dynamics of Pf-RBCs, and may lead to obstruction of small capillaries 139 impairing the ability of RBCs to circulate.
In this section, we present results of the application of the computational framework we discussed for healthy RBCs to Pf-RBCs. In particular, we first consider single RBCs for validation purposes and subsequently we present simulations for whole infected blood as suspension of a mixture of healthy and Pf-RBCs in order to investigate its rheological behavior.
In Bow et al. 15 a new microfluidic device with periodic obstacles to red blood cell flow was employed to perform experiments for the late ring-stage P. falciparum-infected RBCs that are infected with a gene encoding green fluorescent protein (GFP). This device consisted of triangular obstacles (in converging and diverging form). For both the converging and diverging geometries infected RBCs exhibit lower average velocities than healthy RBCs (see Fig. 10a). In the DPD simulations, the infected cells are modeled with increased shear modulus and membrane viscosity values obtained from optical tweezers experiments. The parasite was modeled as a rigid sphere of two microns in diameter 48 placed inside the cell. Time lapse images from the simulations showing passage of an infected RBC through the periodic obstacles with converging and diverging opening geometries are shown in Fig. 10b.
The DPD model is able to capture the effect of changes of RBC properties arising from parasitization on the movement of healthy RBCs and Pf-RBCs quite accurately. A quantitative comparison of the simulation results with experimental data for the average velocity of uninfected RBCs and Pf-RBCs as a function of applied pressure gradient is shown in Fig. 11.
In order for the DPD simulation to provide insights into the design of microfluidic devices that are potentially capable of diagnosing disease states, it is essential to develop computational capabilities for the biorheology of whole blood (containing multiple components) where only a small volume fraction of the total cell population contains diseased cells. Next we present simulation results of blood flow in malaria as a suspension of healthy and Pf-RBCs at the trophozoite stage and hematocrit H t = 0.45. Several realistic (low parasitemia levels) and hypothetical parasitemial levels (percentage of Pf-RBCs with respect to the total number of cells in a unit volume) from 5 to 100% are considered in vessels with diameters 10 and 20 lm. The inset of Fig. 12a shows a snapshot of RBCs flowing in a tube of diameter 20 lm at a parasitemia level of 25%.
The main result in Fig. 12a is given by the plot of the relative apparent viscosity in malaria obtained at different parasitemia levels. The effect of parasitemia level appears to be more prominent for small diameters and high H t values. Thus, at H t = 0.45 blood flow resistance in malaria may increase up to 50% in vessels of diameters around 10 lm and up to 43% for vessel diameters around 20 lm. These increases do not include any contributions from the interaction of Pf-RBCs with the glycocalyx 123,160 ; such important interactions are complex as they may include cytoadhesion.
In Fig. 12b we also present the bulk viscosity of infected blood (schizont stage) simulated in a Couette type device at shear rate _ c ¼ 230 s À1 . The DPD simulations compare favorably with the experimental data obtained with a corresponding rheometer in Raventos-Suarez et al. 132 These validated predictions were obtained without an explicit adhesion model between Pf-RBCs. It seems that such cell-cell interactions are not important at this high shear rate value.
With regards to the cytoadherence of Pf-RBCs, microfluidic experiments have been conducted in Antia et al. 8 to investigate the enhanced cytoadherence of Pf-RBCs in flow chambers. These experiments revealed that the adhesive dynamics of Pf-RBCs can be very different than the well-established adhesive dynamics of leukocytes. For example, the adhesive dynamics of Pf-RBCs on purified ICAM-1 is characterized by stable and persistent flipping (rolling) behavior for a wide range of wall shear stresses (WSS) 8 but also by intermittent pause and sudden flipping due to the parasite mass inertia. This interesting adhesive dynamics was simulated in Fedosov et al. 60 where a systematic study was conducted to estimate the magnitude of the adhesion force. It was found to be in agreement with the adhesion force measured in related experiments by Cravalho et al. 31 using an atomic force microscope (AFM). A Pf-RBC may exhibit firm adhesion at a WSS lower than a certain threshold and can completely detach from the wall at higher WSS. At low WSS, adhesion forces are strong enough to counteract the stress exerted on the cell by the flow resulting in its firm sticking to the wall. On the contrary, at high WSS existing bonds do not provide sufficiently strong adhesive interactions which yields RBC detachment from the wall. RBC visualizations showed that its detachment at high WSS occurs during the relatively fast motion of RBC flipping, since the contact area at that step corresponds to its minimum. 58 However, in experiments 8 Pf-RBCs which moved on a surface coated with the purified ICAM-1 showed persistent and stable rolling over long observation times and for a wide range of WSS values. This suggests that there must be a mechanism which stabilizes rolling of infected RBCs at high WSS. This fact is not surprising since, for example, leukocyte adhesion can be actively regulated depending on flow conditions and biochemical constituents present. 64

Rheology in Sickle Cell Anemia
The abnormal rheological properties of sickle cell RBCs (SS-RBCs) are correlated with the changes in their shapes and with the stiffened cell membrane. This was measured by micropipette experiments in Itoh et al. 82 at different deoxygenated stages. The shear modulus of the full deoxygenated sickle cell falls within a wide range of values depending on the intracellular sickle hemoglobin HbS polymerization and it is more than 100 times the value of healthy cells. In the simulation studies in Lei and Karniadakis, 94 the shear modulus of the full deoxygenated sickle cell was set at very high value, e.g., 2000 times the value of the healthy cells to ensure a fully rigid SS-RBC. The bending rigidity of sickle cells under different deoxygenations was also set to be 200 times the value of the healthy RBC as no detailed values are available in the literature. In Lei and Karniadakis 94 sickle blood was modeled by a suspension of SS-RBCs in a solvent, which is represented by a collection of coarse-grained particles with DPD interactions. Different sickle cell morphologies were obtained based on a simulated annealing procedure and medical image analysis of clinical data. 94 To investigate the relationship between the effect of the rate of deoxygenation and the rheology of sickle blood, Kaul et al. 87 examined the shear viscosity of sickle blood subjected to both fast and gradual deoxygenation procedures. Sickle blood subjected to gradual deoxygenation procedure showed monotonic elevation of shear viscosity and the formation of the sickle shape of blood cells over a period of 30 minutes until the full deoxygenated state was achieved. On the contrary, sickle blood subjected to the fast deoxygenation procedure exhibits two distinct phases. The shear viscosity of the sickle blood showed fast elevation within the first 7 minutes of deoxygenation accompanied with the cell morphology transition to granular shape. However, the shear viscosity decreased gradually during further deoxygenation. A large portion of cells appears extremely elongated with the intracellular HbS fibers aligned in one direction. To study the morphology effect on the rheological behavior of the sickle blood, in Lei and Karniadakis 94 simulation results were presented for shear flow of the sickle blood with the three distinct types of sickle cell reported in the experiment (Hct = 40%). Figure 13 plots the shear viscosity of the deoxygenated sickle blood; the shear modulus of the cell membrane is the same for all the three types. The sickle blood shows elevated and shearindependent viscosity values for all three types. Moreover, the sickle blood with granular and sickle morphology shows larger viscosity compared with the elongated shape, which explains the progressive decrease of the viscosity value with further deoxygenation, since a large portion of granular cells transforms into the elongated shape during the procedure. This result is probably due to the different effective volume for each type of the sickle blood in the shear flow system, 87 which affects the momentum transport ability between the cells.
The hemodynamics of sickle blood was studied in an isolated vasculature in Kaul et al. 86 While the oxygenated sickle blood exhibits hemodynamics similar to healthy blood, the deoxygenated sickle blood shows distinctive dynamic properties for different RBC types. In the simulation, sickle blood in a tube flow system with Hct = 30% was considered similar to the experiment although precise morphology details are difficult to extract from the experiments. The deoxygenated sickle blood flow was represented by a suspension of RBCs with sickle and granular shapes. Blood plasma and cytosol are explicitly represented by DPD, and they are separated by the cell membrane through the bounceback reflection on membrane surface. The viscosity of the cytosol is set to 4g 0 and 50g 0 for the healthy and deoxygenated blood flow, where g 0 is the viscosity of the blood plasma. Figure 14 plots the increase of the flow resistance with different oxygen tension for the sickle and granular shapes. While both types of blood flow show further increase in flow resistance at deoxygenated state, the granular type of blood flow shows a more pronounced elevation compared with the sickle shape. sickle and granular cells in the tube flow. The cells of sickle shape tend to flow along the axis of the tube as observed in experimental studies in LaCelle 92 ; specific distributions of SS-RBC orientation were presented in Lei and Karniadakis. 94 In this section we presented DPD simulations of sickle blood flow and compared them with existing experimental results primarily on flow resistance. New microfluidic experiments, like the ones presented in the previous section on malaria, will be very useful in further validation of the DPD models, especially as a function of deoxygenation rate. Such experiments on a microfluidic cytometer are currently under way by the MIT coauthor of this paper. In addition, quantifying-both experimentally and through DPD simulations-the SS-RBC membrane fluctuations will characterize better the precise state of the SS-RBCs. In a recent study in Byun et al., 18 it was found that the membrane fluctuations are significantly reduced for SS-RBCs taken from a sickle cell disease patient at ambient oxygen concentration. Careful analysis identified that high cytoplasmic viscosity is primarily responsible for the decreased membrane fluctuations, while irreversibly sickled cells were found to have higher membrane stiffness than that from other types of RBCs in sickle cell disease. Similarly, new microfluidic experiments are needed, in addition to the published results, 77 to elucidate the mechanism of cytoadhesion and occlusion in sickle cell blood flow, and to validate the recent findings of Lei and Karniadakis 95 on the precise mechanisms of vaso-occlusion in postcapillaries.

DISCUSSION
We have presented a new computational framework for simulating hematologic disorders based on DPD, derived from the molecular dynamics method. This is a Lagrangian mesoscopic approach that can model seamlessly blood cells, proteins, plasma, and arterial walls, and can predict accurately blood flow in health and disease. To make this approach patient-specific we rely on single-RBC experiments, e.g., using optical tweezers, atomic force microscopy, etc., and novel microfluidic devices in order to obtain data from which we can extract the macroscopic parameters of the model, which, in turn, can be related to the microscopic parameters required in DPD modeling. Moreover, these single-RBC data sets can serve as a validation testbed over a wide range of operating conditions. The utility of the DPD models is then to predict whole blood behavior in health or disease without any further ''tuning'' of the models' parameters. We demonstrated that this is indeed the case for healthy and malaria-infected as well as sickle whole blood. We can also employ the DPD models for sensitivity studies, hence interacting more meaningfully with the experiments, and suggest specific measurements or parameter variations or even probe certain biophysical interactions, e.g., the interaction of the lipid bilayer with the cytoskeleton in a RBC, 120 not accessible with experimental techniques alone.
The two diseases we studied have some common biophysical characteristics, however, their biorheology is distinctly different. Hence, combining computationalexperimental studies, as we have been pursuing in our research groups, can help greatly in quantifying fundamental biophysical mechanisms and possibly identifying unique biomarkers for these diseases. Our computational framework is general and employs patient-specific model input, hence it can be also used for other hematologic diseases, such as diabetes, HIV, and even in cancer metastasis. In the flow resistance results we presented in the current paper, we modeled FIGURE 14. Increase of the flow resistance induced by the sickle blood flow for both sickle (left) and granular (right) shapes. The inset plot shows a snapshot of the sickle cells in the tube flow. The figure is reproduced with permission from Lei and Karniadakis. 94 whole blood as suspension of RBCs in plasma, hence ignoring the effect of white cells (about 0.7%) and platelets (less than 0.5%), also the effect of other proteins in the plasma, or the interaction with the endothelial surfaces of the blood vessel walls. From the numerical modeling standpoint, there is no particular difficulty in also modeling these other cells as we have done for margination studies in the subsection ''Margination of White Blood Cells and Platelets'', or modeling the endothelium, which may contribute significantly in specific biomedical studies, e.g., in thrombosis, immune response. From the biorheological view point, however, as we demonstrated for both malaria and sickle cell anemia, the presence of the other cells is not important unlike the endothelium that can modify significantly the blood flow resistance.