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.

1. 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 non-axisymmetric 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 vaso-occlusive 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 2-week 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.

2. Materials and Methods

2.1. 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).

2.2. 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 Figure 1.

FIGURE 1

Figure 1. Capillary intravitral microscopy. Intravital microscopy of capillary plug flow in vivo in the rat cremaster muscle.

2.3. 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.

2.4. 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.

2.5. 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⁻¹) using a cone-plate viscometer (DV II + Pro, Brookfield, Middleboro, MA, USA).

2.6. 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).

2.7. 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: (1) contrast thresholding ([μ_i − σ_i, μ_i] → [0, 255]), (2) Gaussian 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.

FIGURE 2

Figure 2. Image processing and contour extraction. Novel image processing were applied to obtain the erythrocyte contour for analysis. Briefly, image processing to obtain binary images involved the following: (1) contrast thresholding ([μ_i − σ_i, μ_i] → [0, 255], (2) Gaussian Smoothing (kernel: 5, σ = 7), (3) Thresholding (T = 68), (4) morphological opening (kernel: 3x3), and (5) binary edge detection for contour extraction.

2.7.1. 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:

$\begin{array}{ll}S(f_{ij})=(f_{ij}-a)\frac{(b-a)}{255}& (2.1)\end{array}$

$\begin{array}{ll}{g}_{ij}=\{\begin{array}{cc}S(f_{ij})\hfill & 0\le S(f_{ij})\le 255\hfill \\ 0\hfill & S(f_{ij})<0\hfill \\ 255\hfill & S(f_{ij})>255\hfill \end{array}\hfill & (2.2)\end{array}$

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.

2.7.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

$\begin{array}{ll}{L}_{ij}=G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x^2+y^2)}{2{\sigma }^2}}& (2.3)\end{array}$

$\begin{array}{ll}{K}_{kl}=\frac{1}{{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}(L_{ij})}{L}_{kl}& (2.4)\end{array}$

Gaussian smoothing was performed by convolving the image, f, with the kernel K, g = K ∘ f.

2.7.3. Binary Thresholding and Morphological Opening

Binary thresholding for edge detection involved application of a thresholding function T(f_ij), as follows:

$\begin{array}{ll}T(f_{ij})=\{\begin{array}{cc}1& f_{ij}\ge \beta \\ 0& f_{ij}<\beta \end{array}& (2.5)\end{array}$

$\begin{array}{ll}{g}_{ij}=T(f_{ij})& (2.6)\end{array}$

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:

$\begin{array}{ll}{g}_{ij}=\{\begin{array}{cc}1& S“\text{fits}”P\\ 0& otherwise\end{array}& (2.7)\end{array}$

where P is defined as:

$\begin{array}{ll}{P}_{kl}=f_{mn}\{\begin{array}{c}k\in [1,M]\hfill \\ l\in [1,N]\hfill \\ m\in [i-\frac{(M-1)}{2},i+\frac{(M-1)}{2}]\hfill \\ n\in [j-\frac{(N-1)}{2},j+\frac{(N-1)}{2}]\hfill \end{array}\hfill & (2.8)\end{array}$

A morphological dilation of an image f, defined as g = f ⊕ s, creates a new image g defined by:

$\begin{array}{ll}{g}_{ij}=\{\begin{array}{cc}1& S“\text{hits}”P\\ 0& otherwise\end{array}& (2.9)\end{array}$

where P is defined as:

$\begin{array}{ll}{P}_{kl}=f_{mn}\{\begin{array}{c}k\in [1,M]\hfill \\ l\in [1,N]\hfill \\ m\in [i-\frac{(M-1)}{2},i+\frac{(M-1)}{2}]\hfill \\ n\in [j-\frac{(N-1)}{2},j+\frac{(N-1)}{2}]\hfill \end{array}\hfill & (2.10)\end{array}$

The process of morphological opening to eliminate shared boundaries between neighboring cells is thus defined as:

$\begin{array}{ll}g=S\oplus (S\ominus f)& (2.11)\end{array}$

Morphological opening was applied to obtain isolate object boundaries for edge detection.

2.7.4. 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.

FIGURE 3

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.

2.8. 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₁, X₂) from a Lagrangian perspective.

$\begin{array}{ll}{F}_{iR}=F_{iR}(X_1,X_2)& (2.12)\end{array}$

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₁, X₂), is a function of the previous frame. Mathematically, for the Nth frame, this is written as:

$\begin{array}{ll}{F}_{iR}^N=f(F_{iR}^{N-1})& (2.13)\end{array}$

Let the coordinate system of the (N − 1)th frame be denoted as: $X_R^{N-1}=(X_1^{N-1},X_2^{N-1})$. The deformation gradient tensor of the Nth frame is therefore a function of the coordinate system of the (N − 1)th frame:

$\begin{array}{ll}{F}_{iR}^N=F_{iR}^N(X_1^{N-1},X_2^{N-1})& (2.14)\end{array}$

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 $\overrightarrow{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 $\overrightarrow{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 tensor $\widehat{F_{iR}}$ can be calculated as:

$\begin{array}{ll}\widehat{F_{ab}}{\displaystyle \sum _M}\Delta X_b^{il}\Delta X_j^{il}{w}_M={\displaystyle \sum _M}\Delta x_a^{il}\Delta X_j^{il}{w}_M& (2.15)\end{array}$

In matrix notation, this is:

$\begin{array}{ll}\widehat{\text{F}}{\displaystyle \sum _M}\Delta X^{il}\Delta X^{i{l}^T}{w}_M={\displaystyle \sum _M}\Delta x^{il}\Delta X^{i{l}^T}{w}_M& (2.16)\end{array}$

Therefore, the deformation gradient tensor, $\widehat{F}$ is:

$\begin{array}{ll}\widehat{\text{F}}=A{D}^{-1}& (2.17)\end{array}$

where

$\begin{array}{ll}A={\displaystyle \sum _M}\Delta x^{il}\Delta X^{i{l}^T}{w}_M& (2.18)\end{array}$

and

$\begin{array}{ll}D={\displaystyle \sum _M}\Delta X^{il}\Delta X^{i{l}^T}{w}_M& (2.19)\end{array}$

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):

$\begin{array}{ll}w(r)=\{\begin{array}{cc}1-6r^2+6r^3\hfill & r\le \frac{1}{2}\hfill \\ 2-6r+6r^2-2r^3\hfill & \frac{1}{2}<r<1.0\hfill \\ 0\hfill & r\ge 1.0\hfill \end{array}\hfill & (2.20)\end{array}$

where r is the adjusted radius of the cluster defined as:

$\begin{array}{ll}r=\frac{r_{gl}-r_{g1}}{r_{cut}}& (2.21)\end{array}$

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:

$\begin{array}{ll}ϵ=\frac{1}{2}(\text{F}+{\text{F}}^T)-\text{I}& (2.22)\end{array}$

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:

$\begin{array}{ll}ϵ=\frac{\beta w_{n}l}{A}(\frac{1}{2}(\text{F}+{\text{F}}^T)-\text{I})& (2.23)\end{array}$

where A is total cell area and l is total cell perimeter.

2.9. 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² 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:

$\begin{array}{ll}\rho {\text{b}}_X=\langle -\frac{1}{w_n}(P_{ext}+P_{int}cos\theta ),-\frac{1}{w_n}{P}_{int}sin\theta \rangle & (2.24)\end{array}$

where ρ is membrane density, approximated as 1,300 kg/m³ (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.

$\begin{array}{l}{P}_{ext}(t)\hfill \\ =\{\begin{array}{cc}\frac{23\mu (V_C(t)+\alpha (t))}{D_T}\frac{L_T}{D_T}(1-\frac{D_T}{L_T}\frac{0.8V_C(t)}{V_C(t)+\alpha (t)})\hfill & X_1\le 0\hfill \\ 0\hfill & X_1>0\hfill \end{array}\hfill & (2.25)\end{array}$

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.

FIGURE 4

Figure 4. Free body diagram for determination of body forces. Body forces were calculated assuming constant internal pressure and external pressure determined as a function of time.

2.10. Determination of Mechanical Constants

From the Navier-Cauchy Equations of Elastodynamics:

$\begin{array}{ll}(\lambda +\mu )\nabla (\nabla ·\overrightarrow{u})+\mu {\nabla }^2\overrightarrow{u}+\rho \overrightarrow{b}=\rho \frac{{\partial }^2\overrightarrow{u}}{\partial t^2}& (2.26)\end{array}$

where $\overrightarrow{u}=\langle u(X_1,X_2),v(X_1,X_2)\rangle$ is the vector of displacements in the X₁ and X₂ directions, ρ is the average membrane density, approximated as 1,300 kg/m², $\overrightarrow{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_ndl, which is the area of the infinitesimal piece of membrane for which these equations must hold yields:

$\begin{array}{ll}(\lambda +\mu )\nabla (\nabla ·\overrightarrow{u})dA+\mu {\nabla }^2\overrightarrow{u}dA+\rho \overrightarrow{b}dA=\rho \frac{{\partial }^2\overrightarrow{u}}{\partial t^2}dA& (2.27)\end{array}$

When this occurs, the body forces simplify to the following expressions:

$\begin{array}{ll}\rho b_{X_1}dA=-(P_{ext}+P_{int}cos\theta )dl& (2.28)\end{array}$

$\begin{array}{ll}\rho b_{X_2}dA=-w_n{P}_{int}sin\theta dl& (2.29)\end{array}$

Dividing by w_n yields the following simplification of the Navier-Cauchy equations:

$\begin{array}{l}(\lambda +\mu )\nabla (\nabla ·\overrightarrow{u})dl+\mu {\nabla }^2\overrightarrow{u}dl\\ -\frac{dl}{w_n}\langle P_{ext}+P_{int}cos\theta ,P_{int}sin\theta \rangle =\rho \frac{{\partial }^2\overrightarrow{u}}{\partial t^2}dl& (2.30)\end{array}$

where dl can be approximated as:

$\begin{array}{ll}dl\approx \frac{l}{M}& (2.31)\end{array}$

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:

$\begin{array}{l}\left(\lambda +\mu \right)(\frac{{\partial }^{2}u}{\partial X_1^2}+\frac{{\partial }^{2}v}{\partial X_1\partial X_2})dl+\mu (\frac{{\partial }^{2}u}{\partial X_1^2}+\frac{{\partial }^{2}u}{\partial X_2^2})dl\\ -\frac{dl}{w_n}\left(P_{ext}+P_{int}cos\theta \right)=\rho \frac{{\partial }^{2}u}{\partial t^2}dl& (2.32)\end{array}$

$\begin{array}{l}\left(\lambda +\mu \right)(\frac{{\partial }^{2}u}{\partial X_2\partial X_1}+\frac{{\partial }^{2}v}{\partial X_2^2})dl+\mu (\frac{{\partial }^{2}v}{\partial X_1^2}+\frac{{\partial }^{2}v}{\partial X_2^2})dl\\ -\frac{dl}{w_n}{P}_{int}sin\theta =\rho \frac{{\partial }^{2}v}{\partial t^2}dl& (2.33)\end{array}$

When this is rewritten as a function of the components of the displacement gradient tensor, as well as the body forces $\overrightarrow{b}$, the equations simplify to:

$\begin{array}{l}(\lambda +2\mu )\frac{\partial G_{11}}{\partial X_1}dl+(\lambda +\mu )\frac{\partial G_{22}}{\partial X_1}dl+\mu \frac{\partial G_{12}}{\partial X_2}dl\\ -\frac{dl}{w_n}(P_{ext}+P_{int}cos\theta )=\rho \frac{{\partial }^{2}u}{\partial t^2}dl& (2.34)\end{array}$

$\begin{array}{l}(\lambda +2\mu )\frac{\partial G_{22}}{\partial X_2}dl+(\lambda +\mu )\frac{\partial G_{11}}{\partial X_2}dl+\mu \frac{\partial G_{21}}{\partial X_1}dl\\ -\frac{dl}{w_n}{P}_{int}sin\theta =\rho \frac{{\partial }^{2}v}{\partial t^2}dl& (2.35)\end{array}$

Using Cauchy boundary conditions, as the values of $\overrightarrow{u}=\langle u(X_1,X_2),v(X_1,X_2)\rangle$ 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_2^j=f_j(X_1^j)$ 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_{RS}^j(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\overrightarrow{x}=\overrightarrow{b}$, where A is the matrix of convective forces $\overrightarrow{x}$ is the vector of Lamé's Constants, $\left[\begin{array}{l}\lambda \\ \mu \end{array}\right]$, and $\overrightarrow{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:

$\begin{array}{l}dl\left[\begin{array}{cc}\frac{\partial G_{11}}{\partial X_1}+\frac{\partial G_{22}}{\partial X_1}& 2\frac{\partial G_{11}}{\partial X_1}+\frac{\partial G_{22}}{\partial X_1}+\frac{\partial G_{12}}{\partial X_2}\\ \frac{\partial G_{22}}{\partial X_2}+\frac{\partial G_{11}}{\partial X_2}& 2\frac{\partial G_{22}}{\partial X_2}+\frac{\partial G_{11}}{\partial X_2}+\frac{\partial G_{21}}{\partial X_1}\end{array}\right]\left[\begin{array}{c}\lambda \\ \mu \end{array}\right]\\ =\left[\begin{array}{c}\rho \frac{{\partial }^{2}u}{\partial t^2}dl+\frac{dl}{w_n}(P_{ext}+P_{int}cos\theta )\\ \rho \frac{{\partial }^{2}v}{\partial t^2}dl+\frac{dl}{w_n}{P}_{int}sin\theta \end{array}\right]& (2.36)\end{array}$

Taking the Moore-Penrose pseudoinverse of the convective forces yields a numerical approximation for Lame's constants:

$\begin{array}{ll}\left[\begin{array}{c}\lambda \\ \mu \end{array}\right]=\frac{1}{dl}\left[\begin{array}{c}{A}^{+}\end{array}\right]\left[\begin{array}{c}\rho \frac{{\partial }^{2}u}{\partial t^2}dl+\frac{dl}{w_n}(P_{ext}+P_{int}cos\theta )\\ \rho \frac{{\partial }^{2}v}{\partial t^2}dl+\frac{dl}{w_n}{P}_{int}sin\theta \end{array}\right]& (2.37)\end{array}$

where [A⁺] isthe 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 μ.

2.11. 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.

$\begin{array}{ll}{\text{T}}_{\text{ij}}=\lambda {{ϵ}_{\text{aa}}}^k{\delta }_{ij}+2\mu {{ϵ}_{\text{ij}}}^k& (2.38)\end{array}$

where ϵ_ij^k is the updated Lagrangian Green strain tensor for the kth frame. Wall shear stress, τ₁ is therefore geometrically determined as:

$\begin{array}{ll}{\tau }_1=T_{12}cos\theta -(T_{11}-P_{\text{ext}})sin\theta & (2.39)\end{array}$

To account for the assumed constant membrane thickness, the shear component of membrane tension, t₁, was approximated by integrating wall shear stress over the thickness of the membrane:

$\begin{array}{ll}{t}_1={\displaystyle {\int }_0^{w_n}}{\tau }_{1}dw\approx {\tau }_1{w}_n& (2.40)\end{array}$

2.12. Principle Component Analysis

Upon a 40% exchange transfusion with 2-week stored RBCs results, a heterogeneous population of RBCs was present in the circulation. Principle component analysis (PCA) was performed on vectors containing the 2 Lame constants, the Young's modulus, Bulk modulus, shear stress and shear strain, for each cell. The two principle components with the largest variance were separated, for each cell, in transfusion videos, and compared with a PCA model developed from fresh cell videos, used as a baseline. Principle components greater than 3 standard deviations from the mean baseline principle components were treated as a separate population for all calculations. All PCA was performed utilizing MATLAB (Mathworks, Natick, MA).

2.13. 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).

3. 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.

3.1. 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.

FIGURE 5

Figure 5. Principle component analysis of native and stored cell mechanical properties post transfusion. Principle component analysis of 6 numerically determined mechanical property variables was utilized to distinguish populations of cells post transfusion. The two principle components with largest variance were considered for analysis. (A) A two component PCA model was constructed for the population of native RBCs (n = 32). Two (red) and three (blue) standard deviations from the mean principle component vector are shown. (B) Two component PCA was performed on the population of erythrocytes analyzed in animals transfused with stored cells. Two (red) and three (blue) standard deviations from the mean principle component vector of native RBCs are plotted for reference. RBCs with both principle components greater than 3 standard deviations from the mean principle component of native RBCs were considered a separate population and are shown in red. Approximately 11/35 cells were found to satisfy this criteria.

3.2. 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-Lame equations of elastodynamics. Results for Lame's first constant, Shear Modulus, Bulk Modulus, and Young's Modulus are summarized in Figure 6. Average Lame's first constant was calculated as 35 ± 8, 22 ± 4, and 170 ± 74 MPa for native RBCs (NC), SC1, and SC2, respectively. Similarly, average shear modulus was calculated as 1.1 ± 0.2, 0.90 ± 0.15, and 12 ± 8 MPa for NC, SC1, and SC2, respectively. Average young's modulus was calculated as 3.3 ± 0.6, 2.6 ± 0.4, and 32 ± 20 MPa for FC, SC1, and SC2, respectively. Average bulk modulus was calculated as 36 ± 8, 23 ± 4, and 180 ± 73 MPa for NC, SC1, and SC2, respectively. Lame's first constant, Shear Modulus, Bulk Modulus, and Young's Modulus were found to be statistically significantly higher (P < 0.01) for erythrocytes in population SC2 compared to NC. No statistically significant differences were observed for erythrocytes in population SC1 compared to NC, as expected. Together, these results suggest erythrocytes in population SC2 are less compliant and less deformable compared to NC.

FIGURE 6

Figure 6. Erythrocyte mechanical properties. Erythrocyte mechanical properties were determined for Native Cells (NC), RBCs analyzed from animals with 40% exchange transfusion with principle components within 3 standard deviations of the control RBC principle components (SC1), and RBCs analyzed from animals with 40% exchange transfusion with principle components greater than 3 standard deviations of the control RBC principle components (SC2). The following mechanical constants were determined assuming Istotropic Hookean dynamics: (A) Lame's first constant (λ), (B) Shear Modulus (μ), (C) Young's Modulus (E), (D) Bulk Modulus (K).

3.3. 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.

FIGURE 7

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.

FIGURE 8

Figure 8. Time-lapsed wall shear stress distribution for a single erythrocyte.

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 similar mechanical environment and forces during capillary plug flow despite changes in membrane deformability.

FIGURE 9

Figure 9. Wall shear stress and membrane tension. Temporal averages of WSS distributions were quantified and integrated over membrane thickness to determine membrane tension for native (NC), RBCs analyzed from animals with 40% exchange transfusion with principle components within 3 standard deviations of the control RBC principle components (SC1), and RBCs analyzed from animals with 40% exchange transfusion with principle components greater than 3 standard deviations of the control RBC principle components (SC2). (A) Average, (B) Maximum, and (C) Minimum Wall shear stress were not found be statistically different across the three populations of erythrocytes investigated. Furthermore, (D) Average, (E) Maximum, and (F) Minimum membrane tension were not observed to statistically different across the three populations of erythrocytes investigated.

FIGURE 10

Figure 10. Erythrocyte wall shear strain. The updated Lagrangian Green strain tensor was quantified and temporally averaged for native (NC), RBCs analyzed from animals with 40% exchange transfusion with principle components within 3 standard deviations of the control RBC principle components (SC1), and RBCs analyzed from animals with 40% exchange transfusion with principle components greater than 3 standard deviations of the control RBC principle components (SC2). No significant differences were observed in (A) Average and (B) Maximum shear strain across the three populations of eyrthrocytes investigated.

4. 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 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⁻⁶ and 10⁻⁴ 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 vaso-occlusive 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 stress-strain 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.

5. 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.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This work was supported by NIH grants from the Heart Lung and Blood Institute, R01HL126945, and R01-HL138116.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank Cynthia Walser (UC San Diego) for surgical preparation of the animals.

References

Agrawal, R., Sherwood, J., Chhablani, J., Ricchariya, A., Kim, S., Jones, P. H., et al. (2016). Red blood cells in retinal vascular disorders. Blood Cells Mol. Dis. 56, 53–61. doi: 10.1016/j.bcmd.2015.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Asaro, R., Zhu, Q., and Cabrales, P. (2018). Erythrocyte aging, protection via vesiculation: an analysis methodology via oscillatory flow. Front. Physiol. 9:1607. doi: 10.3389/fphys.2018.01607

PubMed Abstract | CrossRef Full Text | Google Scholar

Baez, S. (1973). An open cremaster muscle preparation for the study of blood vessels by in vivo microscopy. Microvasc. Res. 5, 384–394. doi: 10.1016/0026-2862(73)90054-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Bareil, P. B., Sheng, Y., Chen, Y.-Q., and Chiou, A. (2007). Calculation of spherical red blood cell deformation in a dual-beam optical stretcher. Opt. Express 15, 16029–16034. doi: 10.1364/OE.15.016029

PubMed Abstract | CrossRef Full Text | Google Scholar

Baskurt, O. K., and Meiselman, H. J. (2003). “Blood rheology and hemodynamics,” in Seminars in Thrombosis and Hemostasis, Vol. 29, ed O. K. Baskurt (New York, NY: Thieme Medical Publishers, Inc.), 435–450.

Google Scholar

Baumgart, T., Capraro, B. R., Zhu, C., and Das, S. L. (2011). Thermodynamics and mechanics of membrane curvature generation and sensing by proteins and lipids. Annu. Rev. Phys. Chem. 62, 483–506. doi: 10.1146/annurev.physchem.012809.103450

PubMed Abstract | CrossRef Full Text | Google Scholar

Bremmell, K. E., Evans, A., and Prestidge, C. A. (2006). Deformation and nano-rheology of red blood cells: an afm investigation. Colloids Surf. B Biointerfaces 50, 43–48. doi: 10.1016/j.colsurfb.2006.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabrales, P., Tsai, A. G., and Intaglietta, M. (2008). Isovolemic exchange transfusion with increasing concentrations of low oxygen affinity hemoglobin solution limits oxygen delivery due to vasoconstriction. Am. J. Physiol. Heart Circ. Physiol. 295, H2212–H2218. doi: 10.1152/ajpheart.00751.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Donadee, C., Raat, N. J., Kanias, T., Tejero, J., Lee, J. S., Kelley, E. E., et al. (2011). Nitric oxide scavenging by red blood cell microparticles and cell-free hemoglobin as a mechanism for the red cell storage lesion. Circulation 124, 465–476. doi: 10.1161/CIRCULATIONAHA.110.008698

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, E. A., Waugh, R., and Melnik, L. (1976). Elastic area compressibility modulus of red cell membrane. Biophys. J. 16, 585–595. doi: 10.1016/S0006-3495(76)85713-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedosov, D. A., Dao, M., Karniadakis, G. E., and Suresh, S. (2014). Computational biorheology of human blood flow in health and disease. Ann. Biomed. Eng. 42, 368–387. doi: 10.1007/s10439-013-0922-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Fosnaric, M., Penic, S., Iglič, A., Kralj-Iglic, V., Drab, M., and Gov, N. (2019). Theoretical study of vesicle shapes driven by coupling curved proteins and active cytoskeletal forces. Soft Matter 15, 5319–5330. doi: 10.1039/c8sm02356e

PubMed Abstract | CrossRef Full Text | Google Scholar

Fung, Y. C., and Zweifach, B. W. (1971). Microcirculation: mechanics of blood flow in capillaries. Annu. Rev. Fluid Mech. 3, 189–210. doi: 10.1146/annurev.fl.03.010171.001201

CrossRef Full Text | Google Scholar

Gimsa, J., and Ried, C. (1995). Do band 3 protein conformational changes mediate shape changes of human erythrocytes? Mol. Membr. Biol. 12, 247–254. doi: 10.3109/09687689509072424

PubMed Abstract | CrossRef Full Text | Google Scholar

Gokhin, D. S., Nowak, R. B., Khoory, J. A., de la Piedra, A., Ghiran, I. C., and Fowler, V. M. (2015). Dynamic actin filaments control the mechanical behavior of the human red blood cell membrane. Mol. Biol. Cell 26, 1699–1710. doi: 10.1091/mbc.E14-12-1583

PubMed Abstract | CrossRef Full Text | Google Scholar

Gullett, P. M., Horstemeyer, M. F., Baskes, M. I., and Fang, H. (2007). A deformation gradient tensor and strain tensors for atomistic simulations. Model. Simul. Mater. Sci. Eng. 16:15001. doi: 10.1088/0965-0393/16/1/015001

CrossRef Full Text | Google Scholar

Hägerstrand, H., Mrówczyńska, L., Salzer, U., Prohaska, R., Michelsen, K. A., Kralj-Iglič, V., et al. (2006). Curvature-dependent lateral distribution of raft markers in the human erythrocyte membrane. Mol. Membr. Biol. 23, 277–288. doi: 10.1080/09687860600682536

PubMed Abstract | CrossRef Full Text | Google Scholar

Halpern, B. N., Biozzi, G., Benacerraf, B., and Stiffel, C. (1957). Phagocytosis of foreign red blood cells by the reticulo-endothelial system. Am. J. Physiol. Legacy Content 189, 520–526. doi: 10.1152/ajplegacy.1957.189.3.520

PubMed Abstract | CrossRef Full Text | Google Scholar

Hochmuth, R. M. (2000). Micropipette aspiration of living cells. J. Biomech. 33, 15–22. doi: 10.1016/S0006-3495(73)86021-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hochmuth, R. M., Mohandas, N., and Blackshear, P. Jr. (1973). Measurement of the elastic modulus for red cell membrane using a fluid mechanical technique. Biophys. J. 13, 747–762.

PubMed Abstract | Google Scholar

Huttner, W. B., and Zimmerberg, J. (2001). Implications of lipid microdomains for membrane curvature, budding and fission: commentary. Curr. Opin. Cell Biol. 13, 478–484. doi: 10.1016/S0955-0674(00)00239-8

CrossRef Full Text | Google Scholar

Iglič, A. (1997). A possible mechanism determining the stability of spiculated red blood cells. J. Biomech. 30, 35–40. doi: 10.1016/S0021-9290(96)00100-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Iglič, A., Kralj-Iglič, V., and Hägerstrand, H. (1998). Amphiphile induced echinocyte-spheroechinocyte transformation of red blood cell shape. Eur. Biophys. J. 27, 335–339. doi: 10.1007/s002490050140

PubMed Abstract | CrossRef Full Text | Google Scholar

Iragorri, M. A. L., El Hoss, S., Brousse, V., Lefevre, S. D., Dussiot, M., Xu, T., et al. (2018). A microfluidic approach to study the effect of mechanical stress on erythrocytes in sickle cell disease. Lab Chip 18, 2975–2984. doi: 10.1039/C8LC00637G

CrossRef Full Text | Google Scholar

Jani, V. P., Mailo, S., Athar, A., Lucas, A., Williams, A. T., and Cabrales, P. (2017a). Blood quality diagnostic device detects storage differences between donors. IEEE Trans. Biomed. Circuits Syst. 11, 1400–1405. doi: 10.1109/TBCAS.2017.2749304

PubMed Abstract | CrossRef Full Text | Google Scholar

Jani, V. P., Yalcin, O., Williams, A. T., Popovsky, M. A., and Cabrales, P. (2017b). Rat red blood cell storage lesions in various additive solutions. Clin. Hemorheol. Microcirc. 67, 45–57. doi: 10.3233/CH-170248

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelton, J. G. (1987). Platelet and red cell clearance is determined by the interaction of the igg and complement on the cells and the activity of the reticuloendothelial system. Transfusion Med. Rev. 1, 75–84. doi: 10.1016/S0887-7963(87)70008-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim-Shapiro, D. B., Lee, J., and Gladwin, M. T. (2011). Storage lesion: role of red blood cell breakdown. Transfusion 51, 844–851. doi: 10.1111/j.1537-2995.2011.03100.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralj-Iglič, V., Svetina, S., and Žekž, B. (1996). Shapes of bilayer vesicles with membrane embedded molecules. Eur. Biophys. J. 24, 311–321. doi: 10.1007/BF00180372

PubMed Abstract | CrossRef Full Text | Google Scholar

Lanotte, L., Mauer, J., Mendez, S., Fedosov, D. A., Fromental, J.-M., Claveria, V., et al. (2016). Red cells' dynamic morphologies govern blood shear thinning under microcirculatory flow conditions. Proc. Natl. Acad. Sci. U.S.A. 113, 13289–13294. doi: 10.1073/pnas.1608074113

PubMed Abstract | CrossRef Full Text | Google Scholar

Markin, V. (1981). Lateral organization of membranes and cell shapes. Biophys. J. 36, 1–19. doi: 10.1016/S0006-3495(81)84713-3

PubMed Abstract | CrossRef Full Text | Google Scholar

McMahon, H. T., and Gallop, J. L. (2005). Membrane curvature and mechanisms of dynamic cell membrane remodelling. Nature 438:590. doi: 10.1038/nature04396

PubMed Abstract | CrossRef Full Text | Google Scholar

Messmann, R., Gannon, S., Sarnaik, S., and Johnson, R. M. (1990). Mechanical properties of sickle cell membranes. Blood 75, 1711–1717. doi: 10.1182/blood.V75.8.1711.bloodjournal7581711

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohandas, N., and Evans, E. (1994). Mechanical properties of the red cell membrane in relation to molecular structure and genetic defects. Annu. Rev. Biophys. Biomol. Struct. 23, 787–818. doi: 10.1146/annurev.bb.23.060194.004035

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukhopadhyay, R., Lim, H. G., and Wortis, M. (2002). Echinocyte shapes: bending, stretching, and shear determine spicule shape and spacing. Biophys. J. 82, 1756–1772. doi: 10.1016/S0006-3495(02)75527-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Nash, G. B., Johnson, C. S., and Meiselman, H. J. (1984). Mechanical properties of oxygenated red blood cells in sickle cell (HbSS) disease. Blood 63, 73–82. doi: 10.1182/blood.V63.1.73.bloodjournal63173

PubMed Abstract | CrossRef Full Text | Google Scholar

Omori, T., Ishikawa, T., Barthès-Biesel, D., Salsac, A.-V., Imai, Y., and Yamaguchi, T. (2012). Tension of red blood cell membrane in simple shear flow. Phys. Rev. E 86:056321. doi: 10.1103/PhysRevE.86.056321

PubMed Abstract | CrossRef Full Text | Google Scholar

Parrow, N. L., Violet, P.-C., Tu, H., Nichols, J., Pittman, C. A., Fitzhugh, C., et al. (2018). Measuring Deformability and Red Cell Heterogeneity in Blood by Ektacytometry. J. Vis. Exp. 131:e56910. doi: 10.3791/56910

CrossRef Full Text | Google Scholar

Perrotta, S., Gallagher, P. G., and Mohandas, N. (2008). Hereditary spherocytosis. Lancet 372, 1411–1426. doi: 10.1016/S0140-6736(08)61588-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rand, R. P., and Burton, A. C. (1964). Mechanical properties of the red cell membrane: I. Membrane stiffness and intracellular pressure. Biophys. J. 4, 115–135. doi: 10.1016/S0006-3495(64)86773-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Reddy, P. R., Amarnadh, V., and Bhaskar, M. (2012). Evaluation of stopping criterion in contour tracing algorithms. Int. J. Comput. Sci. Inform. Technol. 3, 3888–3894. doi: 10.4108/eai.21-12-2018.159334

CrossRef Full Text | Google Scholar

Schrier, S. L., Rachmilewitz, E., and Mohandas, N. (1989). Cellular and membrane properties of alpha and beta thalassemic erythrocytes are different: implication for differences in clinical manifestations. Blood 74, 2194–2202. doi: 10.1182/blood.V74.6.2194.bloodjournal7462194

PubMed Abstract | CrossRef Full Text | Google Scholar

Secomb, T. W. (1991). Red blood cell mechanics and capillary blood rheology. Cell Biophys. 18:231. doi: 10.1007/BF02989816

PubMed Abstract | CrossRef Full Text | Google Scholar

Secomb, T. W., Fischer, T. M., and Skalak, R. (1983). The motion of close-packed red blood cells in shear flow. Biorheology 20, 283–294. doi: 10.3233/BIR-1983-20303

PubMed Abstract | CrossRef Full Text | Google Scholar

Secomb, T. W., and Hsu, R. (1995). Red blood cell mechanics and functional capillary density. Int. J. Microcirc. 15, 250–254. doi: 10.1159/000179026

PubMed Abstract | CrossRef Full Text | Google Scholar

Seifert, U. (1997). Configurations of fluid membranes and vesicles. Adv. Phys. 46, 13–137. doi: 10.1080/00018739700101488

CrossRef Full Text | Google Scholar

Shapiro, B. S., Benjamin, L. J., Payne, R., and Heidrich, G. (1997). Sickle cell-related pain: perceptions of medical practitioners. J. Pain Symptom Manage. 14, 168–174. doi: 10.1016/S0885-3924(97)00019-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Skalak, R., and Chien, S. (1982). Rheology of blood cells as soft tissues. Biorheology 19, 453–461. doi: 10.3233/BIR-1982-19306

PubMed Abstract | CrossRef Full Text | Google Scholar

Skalak, R., Tozeren, A., Zarda, R., and Chien, S. (1973). Strain energy function of red blood cell membranes. Biophys. J. 13, 245–264. doi: 10.1016/S0006-3495(73)85983-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Stokke, B. T., Mikkelsen, A., and Elgsaeter, A. (1986). Spectrin, human erythrocyte shapes, and mechanochemical properties. Biophys. J. 49, 319–327. doi: 10.1016/S0006-3495(86)83644-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Tran-Son-Tay, R., Sutera, S. P., Zahalak, G. I., and Rao, P. R. (1987). Membrane stress and internal pressure in a red blood cell freely suspended in a shear flow. Biophys. J. 51, 915–924. doi: 10.1016/S0006-3495(87)83419-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Vahidkhah, K., Balogh, P., and Bagchi, P. (2016). Flow of red blood cells in stenosed microvessels. Sci. Rep. 6:28194. doi: 10.1038/srep28194

PubMed Abstract | CrossRef Full Text | Google Scholar

Walani, N., Torres, J., and Agrawal, A. (2014). Anisotropic spontaneous curvatures in lipid membranes. Phys. Rev. E 89:062715. doi: 10.1103/PhysRevE.89.062715

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, T., Shi, H., Phan-Thien, N., Lim, C. T., and Li, Y. (2019). Numerical design of a microfluidic chip for probing mechanical properties of cells. J. Biomech. 84, 103–112. doi: 10.1016/j.jbiomech.2018.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimmerman, J. A., Bammann, D. J., and Gao, H. (2009). Deformation gradients for continuum mechanical analysis of atomistic simulations. Int. J. Solids Struct. 46, 238–253. doi: 10.1016/j.ijsolstr.2008.08.036

CrossRef Full Text | Google Scholar