Numerical Model for the Determination of Erythrocyte Mechanical Properties and Wall Shear Stress in vivo From Intravital Microscopy

The mechanical properties and deformability of Red Blood Cells (RBCs) are important determinants of blood rheology and microvascular hemodynamics. The objective of this study is to quantify the mechanical properties and wall shear stress experienced by the RBC membrane during capillary plug flow in vivo utilizing high speed video recording from intravital microscopy, biomechanical modeling, and computational methods. Capillaries were imaged in the rat cremaster muscle pre- and post-RBC transfusion of stored RBCs for 2-weeks. RBC membrane contours were extracted utilizing image processing and parametrized. RBC parameterizations were used to determine updated deformation gradient and Lagrangian Green strain tensors for each point along the parametrization and for each frame during plug flow. The updated Lagrangian Green strain and Displacement Gradient tensors were numerically fit to the Navier-Lame equations along the parameterized boundary to determined Lame's constants. Mechanical properties and wall shear stress were determined before and transfusion, were grouped in three populations of erythrocytes: native cells (NC) or circulating cells before transfusion, and two distinct population of cells after transfusion with stored cells (SC1 and SC2). The distinction, between the heterogeneous populations of cells present after the transfusion, SC1 and SC2, was obtained through principle component analysis (PCA) of the mechanical properties along the membrane. Cells with the first two principle components within 3 standard deviations of the mean, were labeled as SC1, and those with the first two principle components greater than 3 standard deviations from the mean were labeled as SC2. The calculated shear modulus average was 1.1±0.2, 0.90±0.15, and 12 ± 8 MPa for NC, SC1, and SC2, respectively. The calculated young's modulus average was 3.3±0.6, 2.6±0.4, and 32±20 MPa for NC, SC1, and SC2, respectively. o our knowledge, the methods presented here are the first estimation of the erythrocyte mechanical properties and shear stress in vivo during capillary plug flow. In summary, the methods introduced in this study may provide a new avenue of investigation of erythrocyte mechanics in the context of hematologic conditions that adversely affect erythrocyte mechanical properties.


INTRODUCTION
Red Blood Cells (RBCs) are highly deformable biconcave particles with diameter 8 µm that constitute ∼45% of total blood volume. The mechanical properties and deformability of RBCs are some of the most important determinants of blood rheology and microvascular hemodynamics, and decreased RBC deformability is characteristic of several hematologic disorders (Rand and Burton, 1964;Nash et al., 1984;Fedosov et al., 2014). RBC mechanical properties are derived from their membrane structure, which consists of a spectrin-actin cytoskeleton linked to the lipid bilayer and transmembrane glycoporins via ankyrin proteins (Skalak and Chien, 1982;Stokke et al., 1986). Membrane embedded proteins, including Band 3 proteins, and lipid microdomains have demonstrated to be important determinants erythrocyte membrane mechanics (Markin, 1981;Gimsa and Ried, 1995;Kralj-Iglič et al., 1996;Seifert, 1997;Huttner and Zimmerberg, 2001;McMahon and Gallop, 2005;Hägerstrand et al., 2006;Baumgart et al., 2011;Walani et al., 2014;Fosnaric et al., 2019). Additionally, the actin components are themselves atypical, consisting of short, uniform length filaments arranged in a highly irregular hexagonal lattice, with subpopulations of dynamic filaments highly dependent on the shear modulus and modulus of area compression (Mohandas and Evans, 1994;Gokhin et al., 2015). The spectrin-actin network itself is important for the stabilization of many erythrocyte shapes, including spiculated RBCs (Iglič, 1997;Iglič et al., 1998;Mukhopadhyay et al., 2002). Additionally, the non-linear viscoelastic mechanical properties of the RBC membrane results in context-dependent elastodynamic behavior that stems from changes in shear modulus (Stokke et al., 1986;Secomb and Hsu, 1995). In large vessels, RBCs aggregate together along the central streamline in Rouleaux-stack formations. As vessel diameter decreases and shear rate increases, these formations reversibly break down into flowing discocytes, which successively roll, tumble, and deform in areas proximal to the vessel cell free layer (Lanotte et al., 2016). Furthermore, increasing shear rate results in deformations of RBCs along streamlines, contributing to the non-Newtonian shear-thinning behavior of blood (Secomb, 1991;Baskurt and Meiselman, 2003). Within capillaries, RBCs undergo deformations unlike those in larger microvessels, adopting a nonaxisymmetric parachute-like conformation, which minimizes the shear modulus (Secomb et al., 1983;Secomb and Hsu, 1995). The RBC membrane itself is not static and undergoes tank-treadmilling to compensate for the increase in fluid viscosity. While tank-treadmilling motion increases membrane energy dissipation, cytoplasmic energy dissipation is minimized, preventing hemolysis and capillary plugging (Lanotte et al., 2016). Together, these aspects of RBCs impart the ability for these cells to withstand large shear stresses in the larger vessels while being able to deform through capillaries with ease.
Changes in RBC deformability are characteristics of mechanical insult to RBC membranes in various hemoglobinopathies, vascular disease, and during storage lesion. In Sickle Cell Disease (SCD), HbS polymerization forces erythrocyte dependent shape changes, increasing membrane viscosity and membrane rigidity (Nash et al., 1984). These mechanical changes contribute to the pathogenesis of vasoocclusive crises in patients with SCD, which remains a clinical problem (Messmann et al., 1990). Similar changes in mechanical properties and the mechanical stability of RBCs have been observed in symptomatic thalassemia (Schrier et al., 1989). Changes in RBC deformability have also been found to be significantly associated with the pathogeneis of microvascular stenosis and retinal vascular disorders (Agrawal et al., 2016;Vahidkhah et al., 2016). During storage lesion, irreversible changes to the RBC membrane result in decreased deformability, contributing to the increased hemolysis, microvascular injury, and decreased perfusion observed post transfusion. The pathogenesis of these donor-dependent mechanical changes during storage lesion results from oxidative damage to lipid membranes, which are the primary component of RBC membrane elasticity (Jani et al., 2017a,b). These pathologic changes in RBC mechanical properties reflect the need for the development of in vivo methods for quantification of RBC mechanical properties in the microvasculature.
The viscoelastic mechanical properties of the human RBC membrane, namely deformation, volume loss, and rupture, have been traditionally assessed via the method of micropipette aspiration (Hochmuth, 2000;Baskurt and Meiselman, 2003). Micropipette aspiration involves application of an external pressure on a local area of the membrane, resulting in membrane deformation. Quantification of this deformation allows for the determination of the membrane elastic modulus via the Law of Laplace. However, micropipette aspiration fails to account for the context-dependent mechanics of the RBC and the parachute-like morphology adopted by these cells in capillaries. More sophistical techniques, namely ektacyotmetry, allow for the determination of the mechanical response to shear stress (Parrow et al., 2018). Molecular techniques, namely atomic force microscopy and optical tweezers, allow for increased spatial resolution by determining spatial maps of stress along the RBC membrane (Fedosov et al., 2014). Recently, microfluidic devices have been developed for probing the mechanical properties of RBCs and the effects of mechanical stress on RBCs in hematologic conditions, namely sickle cell disease (Iragorri et al., 2018;Ye et al., 2019). The mechanical properties and the strain energy function for the RBC membrane have been defined for a two-dimensional membrane with variable membrane thickness, physiologically dependent on the density of glycoproteins and transporters at any one given section of the RBC membrane (Skalak and Chien, 1982). In these models, the bending moments are assumed to be isotropic and proportional to the gradient of curvature along the surface, neglecting internal stresses (Secomb, 1991;Secomb and Hsu, 1995). The aforementioned model, known as Skalak's Law, determines the path integral of stress along the membrane axis, or membrane tension, and predicts mechanics dependent on changes in the shear modulus (Mohandas and Evans, 1994). More sophisticated models have incorporated the gel like behavior of the RBC by quantifying cytoplasmic viscosity from hemoglobin (Stokke et al., 1986). Recent studies have incorporated the molecular basis for membrane deformation by creating finite element models for actin-spectrin interactions. These models have been used to support measurements of RBC mechanical properties ex vivo and simulate RBC dynamics during capillary plug flow (Stokke et al., 1986). Particle tracking and dissipative particle dynamics have recently been applied for diagnostic purposes by utilizing semi-continuum models using immersive boundary or front tracking techniques to exploit pathologic changes in RBC membrane rigidity that result from parasitic infection (Fedosov et al., 2014). To our knowledge, however, these models have been all applied in vitro and no methods have been designed to estimate RBC membrane mechanics and wall shear stress in vivo during capillary plug flow.
The objective of this study is to quantify the mechanical properties and wall shear stress experienced by RBC membranes during capillary plug flow in vivo utilizing intravital imaging and computational methods. Capillaries were imaged in the rat cremaster muscle pre-and post-RBC transfusion with 2week stored cells, and mechanical properties, shear stress, and shear strain were quantified. Briefly, RBC membrane shapes were extracted utilizing image processing and a custom tracing and parametrization algorithm. RBC parameterizations were used to determine updated deformation gradient and Lagrangian-Green strain tensors for each point along the parametrization and for each frame during plug flow. The updated Lagrangian Green strain and Displacement Gradient tensors were numerically fit to the Navier-Lame equations along the parameterized boundary to determined Lame's constants. The RBC membrane itself is modeled as a two-dimensional Hookean Isotropic Solid, and the Cauchy stress tensor was determined utilized the Lagrangian Green strain tensor and the numerically determined Lame's constants. This study is the first to provide a comprehensive methodology for quantifying the deformation of RBCs in vivo without the need for physical tracking mechanisms such as fluorescent beads or tags.

Animal Model and Tissue Preparation
Animal handling and care procedures were in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals, 2011. The University of California, San Diego Animal Subjects Committee approved the experimental protocol. Twelve Sprague-Dawley rats weighing 140-190 g were used for intravital microscopy examination; additional animals were used as plasma and RBC donors. Anesthesia was induced using 60 mg/kg i.p of pentobarbital sodium (Nembutal, Lundbeck Inc, Deerfield, IL) and supplemented as needed by i.v. infusion of 20% of the initial dose. Animals were placed on a heating pad to maintain body temperature at 37 • C. The jugular vein and femoral artery were catheterized, and a tracheal tube was inserted to support ventilation. Catheters were pre-filled with heparinized saline solution (30 IU/ml).
The rat cremaster muscle was exteriorized as previously described by Baez (1973). The rat was positioned on a Plexiglas plate, and the cremaster muscle was secured to a heated platform for viewing. The platform enabled muscle temperature maintenance at 35 • C during surgery and the experiment. The muscle was moistened with warm Plasma-Lyte solution (≈ 37 • C) and following surgical preparation, exposure to atmospheric oxygen was prevented by a covering of polyvinyl film (Saran, SC Johnson & Son, Racine, WI).

Intravital Microscopy
The experimental setup consisted of an intravital microscope (Olympus-BX51WI) equipped with a matching long working distance condenser (NA = 0.8, Thorlabs, Newton, NJ). Two magnifications of 2X and 1.5X were installed between the objective (60X, LUMPFL-WIR, NA = 0.8; Olympus) and the high-speed camera, providing a total magnification of 1800X. A mercury arc lamp (100 W, Walker Instruments, Scottsdale, AZ) was used to illuminate the tissue. A 400-nm interference filter (Spectra Physics, no. 59820) was placed above the condenser in the light path to maximize contrast between blood and the surrounding tissue. A high-speed video camera (Fastcam 1024 PCI, Photron USA), equipped with a one-megapixel chip was used for video recording. Additionally, all videos were recorded between 1,000 and 2,000 frames per second and the camera shutter speed was optimized to obtain the highest quality image possible. Videos were saved as 8-bit graysclae image stacks. A representative intravital microscopy time-lapse video is shown in

Blood and Plasma Collection
Whole blood was obtained from rat donors, adult male Sprague-Dawley rats (280-300 g, n = 24). Briefly, rats were anesthetized, left carotid artery catheter was implanted, and blood was allowed to flow into heparinized tubes (sodium heparin 15 IU/ml). RBCs and plasma were separated by centrifugation (2,700 rpm, 7 min). The buffy coat was discarded. RBCs morphology was checked prior to infusion and discarded if cells were crenated. Blood was stored for 2 weeks prior to transfusion.

Isovolumic Exchange Transfusion Protocol
To adjust the Hct syst to the target Hct level the volume of the exchange-transfusion was calculated from the formula, exchange volume = 0.024 x Hct, where V is the equal to the rat's blood volume in ml (estimated as 6% of the weight in grams) and Hct is the difference between the initial and the target Hct (Cabrales et al., 2008). Packed RBCs or plasma were infused through the jugular vein catheter and an equal volume of blood was withdrawn from the femoral artery, simultaneously, at a rate of 0.3 ml/min to prevent rapid shifts in blood volume.

Systemic Parameters and Blood Measurements
Mean arterial pressure (MAP) and heart rate (HR) were recorded using a pressure transducer connected to the femoral artery catheter using a MP-150 (BIOPAC Systems, Goleta, CA). Blood samples were collected from the femoral artery catheter in heparinized capillary tubes (Fisher Scientific, Pittsburgh, PA) to determine Hctsyst by centrifugation (5 min at 8,000 rpm, Readacrit; Clay Adams, Division of Becton Dickinson, Parsipanny, NJ). Additionally, blood mixtures were prepared using packed RBCs and fresh plasma to determine blood and plasma viscosities (at 37 • C and shear rates of between 50 and 450 s −1 ) using a cone-plate viscometer (DV II + Pro, Brookfield, Middleboro, MA, USA).

Experiment Protocol
Imaging for intravital microscopy of isolated capillaries with good image quality, and image contrast were selected. Blood sample and microcirculation videos were taken 10 min after the end of the 40% exchange transfusion. Video recording of these capillaries were taken over 1,632 frames (1,000-2,000 fps).

Cell Tracking and Image Processing
RBCs were manually identified and tracked by the user utilizing a user interface developed in MATLAB (Mathworks, Natick, MA). Following cell selection and tracking, a region of interest containing the cell was generated, and image processing was performed to obtain the parametrization. Briefly, image processing to obtain binary images involved the following: Smoothing (kernel: 5, σ = 7), (3) Thresholding (T = 68), and (4) morphological opening (kernel: 3x3). All image processing was performed in MATLAB (Mathworks, Natick, MA). Contour extraction involved edge detection of the processed binary images, obtained utilizing a MATLAB implementation of the Moore-Neighbor algorithm for binary image edge detection (Reddy et al., 2012). Following edge detection, parameterizations of edge pixel positions were obtained utilizing MATLAB (Mathworks, Natick, MA). Image transformations are expressed in detail in Figure 2.

Contrast Stretching
Contrast enhancement was applied by mapping all the intensity values within a range [a, b], of an image f , to the full range of grayscale intensity values [0, 255] in a new image g. Defining the stretching function as S(f ij ), the mapped range of the intensity values is determined as follows: The histogram of image f is assumed to be normally distributed with mean µ f and standard deviation σ f . Optimal contrast enhancement was determined empirically as [µ f − σ f , µ f ] is mapped to [0,255]. Results for contrast stretching are shown in Figure 2.

Gaussian Smoothing
Gaussian smoothing was applied to obtain well-defined boundaries for the erythrocyte. A Gaussian Kernel K of size [MxM] was utilized, where Gaussian smoothing was performed by convolving the image, f , with the kernel K, g = K • f .

Binary Thresholding and Morphological Opening
Binary thresholding for edge detection involved application of a thresholding function T(f ij ), as follows: A threshold intensity of β = 68 was found to empirically yield the best results. Morphological opening was applied to thresholded binary images to eliminate overlapping boundaries between neighboring erythrocytes. Briefly, a structural element, S, is defined as a binary [MxN] matrix. The kernel is then moved over an image f and if f ij = s ij = 1 for all i ∈ [1, M] and j ∈ [1, N] then, then S is said to "fit" f . If f ij = S ij = 1 for some i ∈ [1, M] and j ∈ [1, N], then S is said to "hit" f . A morphological erosion of an image f , defined as g = f ⊖ s, creates a new image g defined by: where P is defined as: (2.8) A morphological dilation of an image f , defined as g = f ⊕ s, creates a new image g defined by: where P is defined as: (2.10) The process of morphological opening to eliminate shared boundaries between neighboring cells is thus defined as: Morphological opening was applied to obtain isolate object boundaries for edge detection.

Binary Edge Detection
Contour extraction involved edge detection of the processed binary images, obtained utilizing a MATLAB implementation of the Moore-Neighbor algorithm for binary image edge detection (Reddy et al., 2012). Erythrocytes were manually tracked utilizing MATLAB, and contours were extracted. Figure 3 demonstrates the erythrocyte tracking and contour extraction methods utilized.

Determination of the Deformation Gradient and Lagrangian Green Strain Tensor
A discretized deformation gradient tensor, F, was determined from the discrete parametrization by the method of Gullett et al. (2007) and Zimmerman et al. (2009). Briefly, the discrete deformation gradient tensor was defined by minimizing the error between clusters of L points weighted by a cubic spline function. In the case of this analysis, there are N deformed states (different time points) and M points that parametrize the edge, which is being analyzed. Therefore, there are a total of M(N − 1) deformation gradient tensors F iR . As a result, the F iR is a function of the initial position X R = (X 1 , X 2 ) from a Lagrangian perspective.
FIGURE 3 | Cell tracking and contour extraction. A custom program for manual cell tracking was developed in MATLAB. Briefly, an image processing transformation was applied to an intravital microscopy time-lapse video of capillary plug flow, followed by binary edge detection and contour smoothing. The identified contour of the smoothed cell is observed in red.
The deformation gradient tensor is therefore numerically determined between deformed states. Instead of referencing the deformation relative to the initial position, an updated Lagrangian methodology will be used, such that the deformation gradient tensor for all the points in a single frame, F iR = F iR (X 1 , X 2 ), is a function of the previous frame. Mathematically, for the Nth frame, this is written as: Let the coordinate system of the (N − 1)th frame be denoted as: ). The deformation gradient tensor of the Nth frame is therefore a function of the coordinate system of the (N − 1)th frame: (2.14) Using the formulation defined by Gullett et al., a discrete deformation gradient tensor F iR can be defined. The discrete deformation gradient tensor can be defined by minimizing the error in the mapping between clusters of L points. Let X il represent the change in the vector X in the ith point in the jth frame relative to the lth point in the cluster of L points of the undeformed state, and x il represent the change in the vector x in the ith point in the jth frame relative to the lth point in the cluster of L points of the deformed state. The discretized deformation gradient tensorF iR can be calculated as: In matrix notation, this is: Therefore, the deformation gradient tensor,F is: where w M is a weighting function for all M points, giving more weight to the ith point in the cluster of L points. A cubic spline was used for the weighting function w(r): where r is the adjusted radius of the cluster defined as: where r g1 and r gl are the distances between the current point and the first and lth points respectively, and r cut is the cutoff radius, which is the maximum distance between the reference point and any point in the cluster. A large cutoff radius was utilized to account for membrane treadmilling during erythrocyte deformation. Furthermore, to account for discontinuities from image processing at horizontal and vertical edges, posts near this vicinity were ignored. The Cauchy strain tensor was determined assuming the infinitesimal strain assumption from the deformation gradients: To adjust for membrane viscoelasticity, the membrane was treated as a Maxwell solid, such that the elastic strain was an empiric fraction β of the total strain, where β = 0.1. Due to the short time scale, effects of shape memory are neglected.
As consistent with the Maxwell model, we assume only shape memory during the time of imaging. To adjust for membrane geometry, the strain tensor was multiplied by a ratio of the membrane area, with constant thickness w n = 5 nm, and total cell area, determined from image processing: where A is total cell area and l is total cell perimeter.

Determination of Body Forces
Body forces were determined with the following assumptions: (1) Each RBC is a thin, 2D membrane with constant thickness w n = 5 nm, (2) Internal stresses are replaced by a single constant internal pressure P int = 85 dynes/cm 2 acting perpendicular to the membrane (Tran- Son-Tay et al., 1987), (3) The external pressure gradient is replaced by a single horizontal acting pressure acting left of the center of area, P ext . All calculations were performed in Lagrangian coordinates such that the cell center of area was defined as the origin. The free body diagram for the system is shown in Figure 4. Expressions for body forces from the free body diagram were determined to be the following: where ρ is membrane density, approximated as 1,300 kg/m 3 (Tran- Son-Tay et al., 1987). The relative external pressure for an RBC flowing in a capillary was adopted from the method developed by Fung et al., assuming that capillary diameter is approximately equal to that of RBC diameter (Fung and Zweifach, 1971). In this model, arterial pressure is replaced by the arterial-venous pressure gradient, allowing for a spatial external pressure function.
where µ is blood viscosity (2.63 cP), V c is capillary velocity, determined by pixels traveled by a cell in t frame, L t is the length of the capillary, D t is the diameter of the capillary, and α is the RBC displacement shift relative to plasma.

Determination of Mechanical Constants
From the Navier-Cauchy Equations of Elastodynamics: where u = u(X 1 , X 2 ), v(X 1 , X 2 ) is the vector of displacements in the X 1 and X 2 directions, ρ is the average membrane density, approximated as 1,300 kg/m 2 , b are the body forces acting on the membrane, and λ and µ are Lamé's First and Second parameters, respectively (Tran- Son-Tay et al., 1987). Multiplying the equation by dA = w n dl, which is the area of the infinitesimal piece of membrane for which these equations must hold yields: When this occurs, the body forces simplify to the following expressions: ρb X 1 dA = −(P ext + P int cosθ )dl (2.28) ρb X 2 dA = −w n P int sinθ dl (2.29) Dividing by w n yields the following simplification of the Navier-Cauchy equations: where dl can be approximated as: where M is the number of points in a frame, which is constant between frames, and l is the perimeter of the RBC in a single frame. Written out for a problem in two dimensions, these equations simplify to: When this is rewritten as a function of the components of the displacement gradient tensor, as well as the body forces b, the equations simplify to: Using Cauchy boundary conditions, as the values of u = u(X 1 , X 2 ), v(X 1 , X 2 ) are known along the boundary of the membrane at all time points (i.e., frames), as well as the fact that for any of the j of N frames, the boundary condition X j 2 = f j (X j 1 ) must hold along the membrane, the boundary value problem numerically simplifies to a series of 2(N − 1)M equations, where the index i corresponds with one of the M points, and the index j corresponds with one of the N frames. Let G j RS (X 1,i , X 2,i ) denote the displacement gradient tensor G RS at point i in frame j. The boundary value problem can thus be solved at all M points and N frames under Cauchy boundary conditions via second order numerical difference approximations.
In matrix form, these equations are written as A x = b, where A is the matrix of convective forces x is the vector of Lamé's Constants, λ µ , and b is the vector of the sum of the inertial and body forces. Utilizing Equations (2.34) and (2.35), the matrix form can be written as: Taking the Moore-Penrose pseudoinverse of the convective forces yields a numerical approximation for Lame's constants: where A + is the Moore-Penrose pseudoinverse of the matrix A, which is solved at every point, and averaged over all points to obtain a solution for a mean value for λ and µ.

Determination of Wall Shear Stress and Membrane Tension
The updated Cauchy stress tensor was determined utilizing Hooke's Law, assuming the RBC membrane to be a 2D Hookean Isotropic Solid.
where ǫ ij k is the updated Lagrangian Green strain tensor for the kth frame. Wall shear stress, τ 1 is therefore geometrically determined as: τ 1 = T 12 cos θ − (T 11 − P ext ) sin θ (2.39) To account for the assumed constant membrane thickness, the shear component of membrane tension, t 1 , was approximated by integrating wall shear stress over the thickness of the membrane:

Statistical Analysis
Unless otherwise noted, variables are represented as mean ± standard error of the mean (SEM). All statistical calculations were performed with Prism 8 (GraphPad Software, San Diego, CA) and/or MATLAB (Mathworks, Natick, MA). SEM for mechanical properties, shear stress, and shear strain were calculated from the total number of erythrocytes used for analysis. Comparisons between mechanical constants, wall shear stress, and wall shear strain were performed using a student's t-test. Following principle component analysis, comparisons between mechanical constants, wall shear stress, and wall shear strain were performed using a one-way analysis of variance (ANOVA). Post-hoc analyses were performed using the Newman Keuls multiple comparison test. For all statistical tests, P < 0.01 was accepted as statistically significant. All animals passed the Grubbs test, ensuring that all the measured values at baseline were within a similar population (P < 0.01).

RESULTS
A total of seven (n = 7) animals were entered into the study; two or three capillaries were selected for each animal after application of systemic and vessel reactivity inclusion criteria. Data from sixty-seven (n = 67; 32 in the before transfusion and 35 after transfusion were included in the results. There were no significant differences in hematocrit, systemic, or microhemodynamic parameters at baseline among animals.

Principle Component Analysis
Principle component analysis was utilized to identify and separate heterogeneous populations of erythro-cytes introduced in the experimental group (40% exchange transfusion). The two principle components with the largest variance were isolated from the population (n = 32) of native RBCs at baseline to construct a model to identify RBCs native to the circulation post transfusion (Figure 5A). In the population of erythrocytes analyzed from animals transfused with stored cells, RBCs with both principle components greater than 3 standard deviations from the mean principle component of native RBCs were identified and analyzed as a separate population for all calculations (Figure 5B). 68.6% (24/35) of RBCs analyzed from animals subjected to exchange transfusion were found to have both components within 3 standard deviations of the native RBC principle components (SC1). Additionally, 31.4% (11/35) of RBCs analyzed from animals subjected to exchange transfusion were found to have components greater than 3 standard deviations of the native RBC principle components (SC2). These results reflect the relative proportions of erythrocytes in these two populations expected after a 40% exchange transfusion.

Erythrocyte Mechanical Properties
Erythrocyte mechanical properties were determined via numerical methods utilizing the discretized displacement gradient tensor across the erythrocyte membrane over time and the Navier-

Erythrocyte Wall Shear Stress and Shear Strain
The Cauchy stress tensor was determined for every pixel along the erythrocyte membrane over time. Shear stress was estimated from the stress tensor and averaged over time. Spatial maps of temporal average wall shear stress along the erythrocyte membrane are shown in Figure 7. A representative time lapse of a single erythrocyte is shown in Figure 8. Stress concentrators were observed at contact sites between the RBC and the capillary endothelium despite not being specified in the model. Average, maximum, and minimum wall shear stress and membrane tension were quantified and are shown in Figure 9. In these results, maximum shear stress at contact sites between the RBC and the capillary endothelium is suggestive of the von Mises stress criterion of the erythrocyte during capillary plug flow. No statistically significant differences were observed for the average, maximum, and minimum wall shear stress and membrane tension across the three populations of erythrocytes investigated. Furthermore, no statistically significant differences were observed in average shear strain and maximum shear strain in the three populations of erythrocytes investigated (Figure 10). As expected, minimum shear strain was observed to be 0 at points of only tensile stresses for all erythrocytes investigated. Together, these results suggest that each of the cells were subjected to a

DISCUSSION
This study proposes a methodology that allows for the determination of the mechanical properties of RBCs in an in vivo setting, using a purely image processing-based approach. While an attempt was successfully made in order to get absolute quantities capable of characterizing the mechanical properties of RBCs, the main benefit of this approach is its ability to determine differences between distinct groups of RBCs present in the circulation. Furthermore, this approach is also novel in the sense that it does not require the use of physical markers, such as fluorescent beads or tags, in order to determine the deformation of the RBC membrane. This allows to analyze the deformations in the RBC membrane in a highly relevant physiological environment.
This study demonstrates a novel method for the quantification of erythrocyte mechanical properties and shear stress in vivo utilizing image processing and computational methods.
Furthermore, application of PCA to the computationally determined shear stress, shear strain, and mechanical properties, allowed for the identification of heterogeneous populations of RBCs with statistically significant changes in mechanical properties after transfusion. These results indicate that both native erythrocytes and heterogeneous populations of erythrocytes identified by PCA were subjected to similar mechanical environment and forces during capillary plug flow despite the observed decrease in membrane deformability. Consequently, these results indicate that changes to perfusion and oxygen delivery observed by altering erythrocyte mechanical properties may result solely from changes in membrane deformability and physiologic processes important for its regulation, including the reticuloendothelial system. Furthermore, to our knowledge, this is the first study to quantify erythrocyte mechanics and membranes stresses in vivo during capillary plug flow utilizing measured strains. The mechanical properties and membrane stresses quantified from this study may therefore more accurately reflect shear dependent erythrocyte mechanics observed in vivo. The methods introduced in this study may provide a new avenue of investigation of erythrocyte FIGURE 7 | Wall shear stress distributions across erythrocyte membranes. Temporal averages of wall shear stress were visualized across the cell membrane for 9 different RBCs. Top Row-native RBCs (NC). Middle Row-RBCs from the SC1 population. Bottom Row-RBCs from the SC2 population. Stress concentrators were observed at contact sites between the RBC and the capillary endothelium despite not being specified in the model. mechanics in the context of hematologic conditions that adversely affect erythrocyte mechanical properties.
Principle component analysis allowed for the separation and identification of a heterogeneous population of RBCs in the experimental group. Further quantification of the mechanical properties of these heterogeneous RBCs demonstrated decreased compliance. We hypothesize that the erythrocytes with decreased compliance are those that were transfused post storage. It is well-known that RBCs undergo time dependent changes that impair deformability and various biochemical properties, namely pH, ATP, and 2,3 DPG concentration, known as storage lesion. While decay in biochemical properties is uniform across multiple populations of RBCs, our group has recently demonstrated that changes in deformability are variable and donor dependent (Jani et al., 2017a). As lipids are the most compliant of biomolecules, an increase of lipid oxidation injury and subsequent increase in membrane cell density during storage result in decreased compliance. Studies have also suggested that oxidative damage during storage is a pH dependent process and dependent on extracellular chloride concentration during storage, as the chloride-bicarbonate transporter is critical for the maintenance of intracellular pH (Donadee et al., 2011;Kim-Shapiro et al., 2011). Together, changes in pH and membrane composition seem to be the primary mediators of the decrease in deformability observed during storage lesion. To our knowledge, our study is the first to demonstrate, in vivo a heterogeneous population of erythrocytes post transfusion, suggesting that the observed impaired oxygen delivery and perfusion may be dependent on these mechanical changes, further demonstrating the wide applicability of this method.
Various methods have been utilized to quantify erythrocyte mechanical properties and membrane stresses due to the increasing importance of erythrocyte mechanics in the context of hematologic conditions. Of these techniques, micropipette aspiration is the most established and widely accepted (Hochmuth, 2000). Rand et al. reported an elastic modulus in the range of 0.7 to 3 MPa for erythrocytes measured via micropipette aspiration (Rand and Burton, 1964). These results, however, did not account for changes in elastic modulus in response to shear flow observed in vivo. Attempts by Hocmuth et al. to measure flow dependent elastic modulus utilizing a model cell tethered to the surface of a Couette flow channel reported an elastic modulus in the same range, ∼1 MPa (Hochmuth et al., 1973). In the development of a strain energy function for the RBC membrane, Skalak et al. utilized a constant elastic modulus of 1 MPa, which was found to fit experimentally determined strain energy dissipation and membrane tension (Skalak et al., 1973). More recent techniques, namely dual beam optical stretcher and atomic force microscopy (AFM) have attempted to quantify the erythrocyte membrane elastic modulus (Bremmell et al., 2006;Bareil et al., 2007). While optical methods have largely agreed with older measurements, the elastic modulus determined from AFM was found to be 0.35-0.45 MPa and underestimates older techniques by approximately one order of magnitude (Bremmell et al., 2006). Our method approximated a Young's Modulus and Shear Modulus of approximately 1 MPa and 3 MPa, respectively, which agrees with those estimated from micropipette aspiration, rheologic methods, and optical stretching. We hypothesize that the discrepancy in elastic modulus relative to AFM is a consequence of the increase in elastic modulus in response to increase shear stress, as observed in capillary plug flow in vivo. Other mechanical properties traditionally used to asses erythrocyte membrane stiffness include the modulus of area compression, which is reported as the bulk modulus in our experiments. It is well-established that the modulus of area compression and elastic area compressibility are approximately one order of magnitude larger than both Young's Modulus and Shear Modulus (Evans et al., 1976). Our results demonstrate similar findings. In summary, erythrocyte membrane properties determined from our methods closely align with those from other methods despite the incorporation of experimentally measured dynamics observed during capillary plug flow in our method.
Determination of erythrocyte membrane wall shear stress is difficult due to spatial variability in membrane thickness. Membrane tension, defined as stress integrated over the membrane thickness, is therefore a more widely accepted metric. While membrane tension is generally assessed via micropipette aspiration, few studies have been able to obtain the spatial precision required to assess membrane tension during flow. Recently, Omori et al. reported a membrane tensions in the range of 10 −6 and 10 −4 N/m, with the upper threshold correlating with the minimum von Mises stress required for hemolysis (Omori et al., 2012). The membrane tension reported in this study support these findings. We attribute the increased range in our findings to errors in temporal averaging and discontinuities inappropriately corrected for in our numerical methods. As our findings for membrane tension agree with those presented in the literature, we hypothesize that the estimated wall shear stress is representative of shear stress observed in vivo. Of note in our study is the lack of statistically significant differences in wall shear stress, membrane tension, and shear strain. We hypothesize that during capillary plug flow, erythrocyte dynamics are constrained by the mechanical environment within the capillary due to the inability for active capillary dilation or constriction, a consequence of the lack of smooth muscle. Therefore, for RBCs in plug flow, the strain tensor must be conserved across all cells. The increase in mechanical properties by approximately one order of magnitude observed in the secondary heterogeneous population of RBCs identified by PCA is insufficient to result in significant changes in shear stress, thus accounting for the conserved stress among the populations of cells investigated.
The reticuloendothelial system (RES) is responsible for clearance of RBCs with membrane changes, known as eryptosis. Briefly, phagocytic macrophages in the RES mediate eryptosis and red cell clearance via complement and IgG dependent mechanisms (Halpern et al., 1957;Kelton, 1987). Consequently, the RES is responsible for clearance of schistocytes and hemolytic debris. However, recent studies have suggested that oscillatory flow conditions during splenic flow may induce vesiculation to increase RBC lifespan as a protective mechanism (Asaro et al., 2018). While it is well-established that this mechanism is responsible for the pathogenesis of membranopathies, namely hereditary spherocytosis, the protective mechanism of vesiculation for erythrocyte aging explains the increased flux through the RES during transfusion (Perrotta et al., 2008). Heterogeneous populations are identified by splenic macrophages and targeted for clearance. Consequently, the increased lifespan of RBCs exogenously introduced is likely a result. Vesiculation of transfused RBCs may be responsible for their increased lifespan in vivo. The results from this study demonstrate modest increase in membrane stiffness without changes in membrane wall shear stress and shear strain in erythrocyte mechanics analyzed in the experimental group. Our results suggest that protective vesiculation of transfused cells in the RES is responsible for altering erythrocyte membrane mechanics and dampening adverse mechanical injury from heterogeneous RBCs.
Strain was calculated by numerically determining discretized deformation gradient tensors for each pixel identified on the RBC membrane via image processing. To account for treadmilling, the weighting function employed to determine discretized deformation gradients spanned a substantial area of the erythrocyte membrane. Consequently, one of the main advantages of the methods presented here involved approximation of mechanical properties from experimentally estimated strains in vivo. Comparable methods to calculate shear strain and wall shear stress involve in vitro isolation of the erythrocyte or extensive modeling of erythrocyte mechanics. However, the methods presented here account for the unique environment to which the cell is subjected in vivo. As erythrocyte mechanics are shear dependent, we hypothesize that the mechanical properties estimated in these experiments are more representative of RBC mechanics in the circulation. The methodology developed in this study additionally has much utility of understanding erythrocyte mechanics in the context of various hematologic pathologies, namely sickle cell disease (SCD), membranopathies, and hemoglobinopathies, utilizing animal models in vivo. Clinically, the only treatment for membranopathies, like hereditary spherocytosis and elliptocytosis, is splenectomy (Perrotta et al., 2008). Understanding the mechanics of the RES and changes in mechanics at the level of the microcirculation may allow for the development of new therapeutics. Furthermore, in SCD, clinical management for vaso-occlusive crisis due to RBC agglutination and hypoxia involves opioid therapy, resulting in increased rates of addiction among these patients (Shapiro et al., 1997). Application of the methods presented in this study have the potential to discover new methods and therapies for acute vasoocclusive crisis in SCD by understanding changes in erythrocyte mechanics and perfusion.
Limitations of this study include assumptions to account for viscoelasticity and shape memory, due to the short time scales observed, and of Hookean isotropy for the determination of the stress tensor. The strain energy function for erythrocyte membranes is non-linear as assessed by Skalak et al. (1973). In this model, we utilized Lagrangian coordinates to estimate the required mechanical properties given a pressure gradient. Rather than assuming that mechanical shape is coupled with flow, we assumed that the external pressure gradients required for flow mechanically constrain the shape of the erythrocyte membrane. Consequently, given an observed shape, we can estimate the mechanical properties. As a result of these assumptions as well as the 2D analytic methodology utilized, the quantified mechanical constants should be interpreted cautiously. Future directions should utilize the proper adjustments of the stressstrain relationships employed in this study and account for the viscoelastic properties of the membrane. Furthermore, the study assessed RBC mechanics by determining two-dimensional strain tensors at the imaging plane. Future directions should employ more sophisticated modeling and image processing techniques to extrapolate the three-dimensional deformation of the erythrocyte in the capillary. Due to the non-equivalence of 2D and 3D mechanical properties of erythrocytes, this is a significant limitation of this study. While our results demonstrate the presence of a secondary heterogeneous population of cells with different mechanical properties identified by PCA, we are unable to conclude that the rigid cells were indeed the transfused cells. Future directions should aim to identify the etiology and biological mechanisms responsible for the heterogeneity quantified.

CONCLUSION
The study presented here quantified the mechanical properties and wall shear stress experienced by RBC membranes during capillary plug flow in vivo utilizing intravital imaging and computational methods. Erythrocyte mechanical properties were determined via numerical methods utilizing the discretized displacement gradient tensor across the erythrocyte membrane over time and the Navier-Lame equations of elastodynamics. To our knowledge, the methods presented here are the first erythrocyte mechanical properties and shear stress in vivo during capillary plug flow. The methodology developed in this study additionally has much utility of understanding erythrocyte mechanics in the context of various hematologic pathologies, namely sickle cell disease (SCD), membranopathies, and hemoglobinopathies. Future directions should aim to confirm the findings presented here and apply them to investigate pathologies with impaired erythrocyte mechanics.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.