Erythrocyte Aging, Protection via Vesiculation: An Analysis Methodology via Oscillatory Flow

We demonstrate that erythrocyte deformations, specifically of a type as occur in splenic flow (Zhu et al., 2017), and of the type that promote vesiculation can be caused by simple, yet tailored, oscillatory shear flow. We show that such oscillatory shear flow provides an ideal environment to explore a wide variety of metabolic and biochemical effects that promote erythrocyte vesiculation. Deformation details, typical of splenic flow, such as in-folding and implications for membrane/skeleton interaction are demonstrated and quantitatively analyzed. We introduce a theoretical, essentially analytical, vesiculation model that directly couples to our more complex numerical, multilevel, model that clearly delineates various fundamental elements, i.e., sub-processes, that are involved and mediate the vesiculation process. This analytical model highlights particulary important vesiculation precursors such as areas of membrane/skeleton disruptions that trigger the vesiculation process. We demonstrate, using flow cytometry, that the deformations we experimentally induce on cells, and numerically simulate, do not induce lethal forms of cell damage but do induce vesiculation as theoretically forecasted. This, we demonstrate, provides a direct link to cell membrane/skeletal damage such as is associated with metabolic and aging damage. An additional noteworthy feature of this approach is the avoidance of artificial devices, e.g., micro-fluidic chambers, in which deformations and their time scales are often unrepresentative of physiological processes such as splenic flow.


INTRODUCTION AND BACKGROUND
Without a nucleus, a mature erythrocyte (or RBC) contains a cytosol enclosed within a highly flexible cell membrane. This composite membrane, consisting of a lipid bilayer supported by a membrane skeleton, is essential to its structural integrity and stability. The basic picture is that of a skeleton created from junctional complexes (JCs) bound to each other via head-tohead associations of spectrin (Sp), which is anchored to the fluidic lipid bilayer at linkage sites (Mohandas and Evans, 1994;Mohandas and Gallager, 2008;Lux, 2015). RBCs possess one of the best characterized molecular architectures among all cell types and are thus often chosen as a model system to study various functions of cells, especially those that involve mechanical behavior and response. A particularly intriguing process that erythrocytes exhibit is vesiculation, during which part of the phospholipid bilayer separates from the skeleton and the rest of the membrane to create a vesicle (Peng et al., 2010;Zhu et al., 2017). The occurrence of vesiculation depends on the magnitude of the skeleton-bilayer dissociation stressthat develops during deformation -as well as the strength of the skeleton-bilayer connectivity-which may be compromised during the aging process as discussed below. To date there is incomplete knowledge about the detailed process or physiological significance of vesiculation although evidence has mounted as to its role in erythrocyte aging (Fox et al., 1991;Schwarz-Ben Meir et al., 1991;Glaser et al., 1994;Willekens et al., 2003b;Hattangadi and Lodish, 2007;Bosman et al., 2012;Zhu et al., 2017). Herein we present a novel methodology for a coupled experimentaltheoretical study of human erythrocyte vesiculation.
Our recent studies have demonstrated the strong prospects for human erythrocyte vesiculation during flow through the spleen (Zhu et al., 2017). Our numerical simulations, along with our underlying analysis of skeleton-membrane interactions provided vital insight into the coupling of the RBC aging process with metabolic processes that occur during the aging process. These aging processes include, among others, denaturing of hemoglobin, including such occurring via oxidative damage, and binding of denatured Hb to the skeleton-membrane attachment sites causing disruption of the skeleton-membrane attachment (Willekens et al., 2003a,b;Bosman et al., 2012). Erythrocyte aging, involving vesiculation, in turn involves processes such as oxidative stress (ROS) and other reactive stresses such as involves, inter alia, nitrates (NOS), phenyldrazine, Ca 2+ uptake, or Cd. Figure 1 illustrates examples of various key steps of oxidative stress and damage via ROS (Low et al., 1985;Fox et al., 1991;Nagababu and Rifkind, 2000;Cao et al., 2009;Rifkind and Nagababu, 2013;Mohanty et al., 2014). Such effects can be readily explored via our developed methodology as they contribute to vesiculation (Willekens et al., 2003a;Cao et al., 2009;Rifkind and Nagababu, 2013;Mohanty et al., 2014;Zhu et al., 2017).
In the context of splenic flow, our simulations led to a description of what is most likely required to induce vesiculation during splenic flow and provided a paradigm that, in fact, supported the view that erythrocyte vesiculation, including vesiculation during splenic flow, may be a selfprotective mechanism (Willekens et al., 2003b;Bosman et al., 2012;Zhu et al., 2017). In this context, self-protection involves the elimination of removal molecules such as denatured Hb as well as phosphatidylserine (PS) and IgG that are known to be associated with cell removal (Willekens et al., 2003b;Bevers and Williamsonl, 2010;Wieschhaus et al., 2012;Kostova et al., 2015;Bevers and Williamson, 2016;Bevers et al., 2017). Our results, in addition, revealed that as vesiculation occurs, presumably in younger deformable cells, and hemoglobin concentration increases and membrane area decreases, the prospects for vesiculation decreases; hence the self-protective mechanism may be effectively shut off with aging. This is closely linked with a loss in cell deformability that is often related to a loss in cell viability. In addition, our methods should be expected to shed new light on the effects of oxidative damage, caused by reactive oxidative species (ROS), on the vesiculation process (Hattangadi and Lodish, 2007;Marinkovic et al., 2007). Thus, the continued study of the vesiculation process is warranted as it appears so closely tied to cell aging, to cell viability, and cell death. Particularly important is to directly link vesiculation to the vital factors of aging, such as those associated with oxidative damage and a methodology to confirm the various hypotheses of the mechanisms involved.

Background on Extracellular Vesicles: viz. Microvesicles (MV's)
The extracellular space of multicellular organisms contains a variety of species including, inter alia, metabolites, ions, proteins, polysaccharides, etc. and also a large number of mobile membrane bound vesicles that have been collectively referred to as extracellular vesicles (EV's) (Morel et al., 2010;György et al., 2011;Raposo and Stoorvogel, 2012). There are ongoing attempts at classification of EV's where distinctions are based on, for example, size, constituency, and mechanisms of formation (Morel et al., 2010;György et al., 2011;Raposo and Stoorvogel, 2012;Alaarg et al., 2013). For example, exosomes are generally placed in a size range with diameters <100 nm whereas microvesicles (MV's) are generally placed in the diameter range of 100-1,000 nm. Moreover, exosomes are generally created intracellularly and excreted, whereas MV's are formed through budding from the bilipid membrane. Exceptions, however, may exist as we note, for example, the report by Booth et al. (2006) of exosomes being in the size range 50-100 nm budding from T cells. Herein we focus on what we call (following e.g., Morel et al., 2010;György et al., 2011;Raposo and Stoorvogel, 2012;Alaarg et al., 2013), MV's budded from erythrocyte membranes, and generally expected to be in the size range 100-250 nm. Our analysis, however, does not preclude budded vesicles in a size range <100 nm, yet probably not smaller than 40-50nm as discussed below.
leads to membrane loss via vesiculation. These two studies were particularly concerned with aging, not just of the red blood cell, but of people over 70 years of age in whom the effects were accelerated. Still again, Kostova et al. (2015) have recently discussed how the Ca 2+ ionophore ionomycin provides an effective pathway, due to membrane disruption, to vesiculation. Wieschhaus et al. (2012) report of the effects of increased Ca 2+ uptake in degrading attachment proteins such as band 3, 4.1, and ankyrin in the presence of calpain-1 and consequently reducing cell deformability. These, among other observations, provide a basis for our model of budding that begins with localized, i.e., small 20-40 nm, sections of the membrane whose attachments to the skeleton have been disrupted. Specifically, we study herein how vesiculation can be explored using a remarkably simple imposed oscillatory shear flow that allows for readily tailored modes and intensities of deformation.
Finally, we note that microvesicles budded from erythrocytes are typically characterized by exposure of PS on their outer bilayer leaflet (Dasgupta et al., 2009;Bevers and Williamsonl, 2010;Bevers and Williamson, 2016;Bevers et al., 2017). Hence translocation of PS from the inner-to-outer leaflets may be viewed as a part of the vesiculation process or, in fact, a contributor to vesiculation. Two possible pathways for a PS translocation role are: (i) it induces a shape change involving initial curvature of the membrane out from the skeleton (Sheetz and Singer, 1974), and/or (ii) that it acts as an energy source as does skeletal deformation (Sheetz and Singer, 1974;Bevers and Williamson, 2016). Lipid asymmetry represents stored free energy, given the work required to create it against a steep concentration gradient, that can be made available to induce deformation and aid the vesiculation process as described below in section 4.2. Also as we show in section 4.3, precurvature provides a "boost" in the membranes bulging process by which a vesicle forms. This view is supported by Bevers et al. (2017), who demonstrated that in RBCs affected by Scott's syndrome microvesiculation was diminished in frequency (see also Lhermusier et al., 2011). Dasgupta et al. (2009) studied microvesiculation in mice treated with lactadherin and showed a lactadherin deficiency led to increased numbers of microvesicles at steady state.

Background on Recent Theoretical Innovations
Our simulations (Zhu et al., 2017) of erythrocyte flow through simulated venous-like slits of the spleen were performed with a hierarchical model that followed the cell's deformation process in detail. Noteworthy is that this model considers the cell membrane and skeleton separately and explicitly accounts for the lateral mobility of the skeleton's attachment trans-membrane proteins through the membrane. Hence during deformation, the areal density of skeletal attachments may change, viz. decrease due to large areal deformation of the skeleton, and thereby weaken the overall skeleton-membrane attachment as described in detail elsewhere (Peng et al., 2010;Zhu et al., 2017). However, as the trans-membrane attachment proteins have finite (but documented) mobility, the process of changing skeletal density has its own time scale that operates within the deformation's time scale; this too was predicted (see Zhu et al., 2017 and its discussion of the results of Peng et al., 2010Peng et al., , 2011Peng and Zhu, 2013). Now, in this we find that deformation processes in splenic flow, such as what we have dubbed in-folding, produce a "tension" (dubbed negative pressure in Zhu et al., 2017) between the skeleton and membrane that can promote separation that leads to vesiculation (Zhu et al., 2017). However, a key feature of this is that the time scales of splenic flow are such that large changes (decreases) in the areal density of attachment points are unlikely (Zhu et al., 2017). Hence the vital role of attachment disruption via, e.g., binding of denatured Hb, is required. That is, since areal density reduction of attachments does not occur in time during splenic flow deformation restructuring, it must occur by metabolically induced disruption of anchorage points.
Our simulations were quite detailed as to what was required with respect to such disruptions to cause separation. But to study vesiculation by inducing erythrocyte flow through artificial slits (as can be made by pores or slits in commercial filters can be) however, problematic. For example, flow through artificial slits in the splenic range of 1 µm generally leads to cell damage and fragmentation under very different time scales. This brings us to the following novel discovery, prospect and proposal.
We have found, that under certain ranges of frequency and shear amplitudes, simple oscillatory shear flow causes cell deformations whose characteristics (e.g., magnitude of shear deformation and even in-folding) and time scales are comparable to those in splenic flow, and thereby may induce vesiculation. This provides the prospect of following vesiculation without the use of artificial slits or resorting to deformation methods such as micropipette aspiration -that, as it happens, we have already shown can cause vesiculation but whose time scales are quite unlike those of splenic flow. Hence we have discovered the prospect of studying vesiculation in an open environment, devoid of artificial structure, e.g., the slits or pores of man-made filters, in which control of chemistry is readily achievable. This is now illustrated by way of past and new results specific to our modeling in sections 3 and 4.

OVERVIEW OF THE SIMULATION MODEL AND RESULTS
One of the most severe physiological deformations a RBC sustains occurs inside the spleen, where it "squeezes" through slits as narrow as 0.6 micron (Chen and Weiss, 1973;Mebius and Kraal, 2005;Lux, 2015;Zhu et al., 2017). Under this extreme condition the cell undergoes dramatic shear variations, including a novel deformation mode called infolding discovered in recent computational studies (Freund, 2013;Salehyar and Zhu, 2016;Zhu et al., 2017). The mechanical loads on the cell membrane include the external loads (traction and normal stress) from the external/internal fluids, and the internal stresses within the membrane itself. Hereby the internal stresses are categorized into two parts, viz. in-plane and out-of-plane. The in-plane stresses contain an isotropic component, a shear component, and friction between the skeleton and the bilayer. The out-of-plane stress is the normal interaction stress between the skeleton and the bilayer that either pulls them together (association) or pushes them apart (dissociation). No prestress in the skeleton is considered so that the stress-free state of the cell coincides with its natural biconcave state. The membrane thus possesses shape memory.
The model we use is based on a multiscale multi-physics framework as described in detail elsewhere (Zhu et al., 2007(Zhu et al., , 2017Zhu and Asaro, 2008;Peng et al., 2010Peng et al., , 2011Peng and Zhu, 2013). As illustrated in Figure 2, this approach includes three models at different scales: in the complete cell level (Level III) the membrane is modeled as two layers of continuum shells using the finite element method; the constitutive properties of the inner layer (the membrane skeleton) are obtained from a 3D molecular-detailed JC model (Level II); the mechanical properties of Sp, including its folding/unfolding reactions (Bustamante et al., 1994;Lee and Discher, 2001;Zhu and Asaro, 2008), are obtained with a stress-strain model based on Arrhenius equation (Level I). The fluid-cell interaction is mathematically formulated within a low-Reynolds number Stokes/Oseen flow framework and solved using a boundary-element method (for details see Appendix).
As demonstrated in Figure 3A (Salehyar and Zhu, 2016;Zhu et al., 2017), under certain conditions right after the cell passes through the slit a part of its membrane bends inward to form a concave region on its rear surface. Our model also illustrates that due to the increased surface curvature and skeletal energy density, infolding usually coincides with significantly increased dissociation stress between the skeleton and the lipid bilayer. The values of stress we found are orders of magnitude higher than those in normal conditions such as in a steady shear flow field (Zhu et al., 2017).
Further examination shows that in addition to splenic flow, infolding may happen in other scenarios, e.g., in an oscillatory shear flow field. As shown in Figure 3B, when a RBC is put in a shear flow field with sinusoidally varying   (Salehyar and Zhu, 2016) and (B) oscillatory shear flow. In (B) the cytosolic/medium viscosity ratio was 1; the average shear stress applied to the cell was in the rangeτ ∼ 6Pa. The peak shear rate is 1,000 s −1 , and the frequency is 67 Hz. In both figures the color contour corresponds to shear deformation χ (see section 3.2 for the definition of χ) of the skeleton with the value ranging from 1 (in blue) to 3 (in red). χ = λ 1 /λ 2 with λ i , i = 1, 2 being the principal skeletal stretches (see section 3.2).
FIGURE 4 | Contours of skeleton energy density, χ, in an oscillatory shear flow field. The shear rate was σ =2,000s −1 at a frequency of 100 Hz. The cytosolic/medium viscosity ratio was equal to 1; the average shear stress applied to the cell was in the rangeτ ∼ 12Pa. Areas in red color are judged to be those with a high probability for vesiculation (see section 4.4).
shear strength and direction, it undergoes similar infolding behavior, resulting in increased skeleton-bilayer dissociation stress inside the membrane. Indeed, our simulations suggest that within the range of parameters achievable in laboratory conditions, this dissociation stress in oscillatory shear flows attains levels comparable to splenic flow and may thus contribute to vesiculation (see also Figure 7). Meanwhile, the deformations of the cell remain moderate so that the risk of cell bursting is minimal. Figure 4 shows a snapshot of the deformed cell subject to a shear flow field characterized by frequency 100 Hz with a shear rate σ = 2,000s −1 . The contours shown are for skeleton deformation energy density as so labeled with the insert. We note that in this case the maximum skeleton energy density is nearly 5×10 −4 Jm −2 and is located on Figure 10 as point "a"; this is used in later discussion of vesiculation in section 4.4.

Examples of Oscillatory Flow
Still another case is shown is shown in Figure 5 where in this case the shearing rate is well above physiological rates and is set FIGURE 5 | Contours of skeleton energy density, χ, in an oscillatory shear flow field. The shear rate was σ = 5,000s −1 at a frequency of 67 Hz. The cytosolic/medium viscosity ratio is 5; the average shear stress applied to the cell was in the rangeτ ∼ 6Pa. The energy densities range from about 1 × 10 −4 Jm −2 (blue) − 5 × 10 −4 Jm −2 (red).
at σ = 5,000s −1 ; the frequency of flow was 67 Hz. This case differs from those also shown, i.e., in Figures 3, 4, 6 in that a higher viscosity is prescribed for the cell interior. In this case, the contours are of skeleton shear deformation and like those for Figures 3, 6 the area deformations are just above unity. That means that the maximum skeletal energy density is near to point "b, " close to point "a, " of Figure 10.

Additional Perspective on Oscillatory Shear Flow
We believe this ability to tailor deformation via oscillatory flow presents a novel methodology for probing the metabolic effects FIGURE 6 | Deformation of a red blood cell in a strong high-frequency oscillatory shear flow field. The amplitude of the shear flow is 2,000 s −1 and the frequency is 100 Hz. The cytosolic/medium viscosity ratio was equal to 1; the average shear stress applied to the cell was in the rangeτ ∼ 12Pa. The color contour shows the shear deformation χ of the skeleton. Areas in red color are judged to be those with a high probability for vesiculation (see section 4.4). χ = 1 (blue) to 3 (red).
of membrane/skeleton degradation that defines aging. For this reason we have run additional simulations to investigate the response of red blood cells to strong oscillatory shear flows. This has uncovered additional interesting phenomena. For example, in Figure 6 we plot the deformation of a cell to a sinusoidally varying shear flow of 100 Hz and peak shear rate of 2,000 s −1 . It is seen that the cell undergoes complicated deformations with multiple infolded regions on its membrane. These regions are characterized by large curvature as well as significant dissociation stress between the cytoskeleton and the lipid bilayer. Indeed, our simulations indicate that in cases like these the dissociation stress reaches levels of O(100)Pa, comparable to the level that can be reached as the cell passes through inter-endothelial slits in a spleen (Zhu et al., 2017). This provides evidence that, in terms of skeleton-bilayer dissociation, an oscillatory flow field is able to create conditions similar to the physiological conditions of splenic flow.
Furthermore, we carefully monitored maximal areal deformations of the lipid bilayer during these processes and note that even under extreme conditions, such as the case of Figure 6, they were only ∼ 1 − 2%; this is significantly smaller than the approximate areal strain of 10% a cell can sustain under dynamic (i.e., impact) conditions (Li et al., 2013). It is also at, or below, the somewhat lower values in the range 2-3% found under different deformation conditions (Leverett et al., 1972;Sandza et al., 1974;Sutera and Mehrjardi, 1975;Evans et al., 1976;Sutera et al., 1977;Daily et al., 1984;Watanabe et al., 2006;Wantanbe et al., 2007;Meram et al., 2013;Hashimoto, 2014;McNamee et al., 2016;Horobib et al., 2017). We thus conclude that pressure-induced cell bursting due to large area deformation of the bilayer is not likely to occur in these conditions as discussed further in section 4.5 along with additional detail and perspective.

Vesiculation Forecasts With Oscillatory Flow
Based on the analysis of our next section, based in turn on our simulation results of cell deformations presented in this section, we suggest that areas shaded in red and zones of infolding as shown, for example, in Figure 7 are most probable vesiculation sites. The areas in red generally experience skeletal energy densities in the range ǫ 0 ∼ 3 − 5 × 10 −4 Jm −2 and shear stresses in the rangeτ ∼ 5 − 10Pa that as discussed in sections 3, 4.4, and 4.5 promote vesiculation.

PROSPECTS FOR VESICULATION
The model developed herein focuses on the mechanisms of vesiculation as a natural, essential, mechanistic component of the erythrocyte aging process. However, in a more general context the focus is on cellular aging per se, and on the role of oxidative stress and other biochemical stress on aging (Iuchi et al., 2007;Rizvi, 2010, 2011;Mohanty et al., 2014). Such aging effects are manifest in at least two ways, viz. a weakening of the skeletal/bilayer connection as measured by γ and in the size of the initial skeletal/bilayer separation. The particular contribution here is that we provide a detailed, quite specific, yet quite simple, mechanistic pathway by which erythrocyte aging occurs and by which erythrocytes are removed. This process involves a general metabolic aging via ROS and other biochemical activity and in this way its analysis provides still another clue as to the effects of biochemical activity on the aging process (Dröge, 2002;Rattan, 2006;Harman, 2009;Tsuda, 2010). Our unique approach follows vesiculation as promoted by biochemical aging and the imposition of splenic-like deformations using our novel tool of oscillatory shear flow as illustrated in Figure 7.

Vesiculation Model
Imagine that initially a quite localized circular patch of membrane, of radius ℓ 0 , separates from the skeleton, as in the bottom most figure, Figure 8Ai. We can imagine that this occurs due to the binding of Hb (Willekens et al., 2003a,b;Bosman et al., 2012), or by oxidative damage (Rifkind and Nagababu, 2013;Mohanty et al., 2014), or other forms of biochemically induced damage (Allan and Mitchell, 1977;Low et al., 1985;Nagababu and Rifkind, 2000;Cao et al., 2009), as noted above. This induces a bulging out and some curvature -also note that if there is a negative (i.e., separating) pressure this will also contribute to an initial curvature as sketched in the next upward figure -call this curvature C 0 as indicated in Figure 8Aii. We address this in more detail below.
The initial bleb, or "blister, " may propagate due to (i) a negative pressure induced by membrane/skeleton deformation as appears in our simulations (Zhu et al., 2017), (ii) additional binding of Hb or translocation of anionic lipids such as PS (Sheetz and Singer, 1974) leading to additional curvature caused fueled by biochemical energy, and (iii) the release of skeletal energy attached to the membrane following more intense skeletal deformation. Later we suggest a potentially important contribution arising from stored chemical energy associated with the asymmetry in lipid composition between the leaflets of the bilayer. Let us define ǫ 0 as the elastic energy stored in the skeleton per unit area -this energy is released as the membrane separation spreads. This energy indeed depends on the deformation state as discussed in section 3.2. In a simple sense, the skeleton exerts a kind of "pinching action" on the membrane causing membrane/skeleton separation. The effect of skeletal deformation has been considered in several studies, e.g., those concerned with erythrocyte shape (Waugh, 1996;Lim et al., 2002;Mukhopadhyay et al., 2002) and vesiculation (Li and Lykotrafitis, 2015). In these citepd studies of minimum energy erythrocyte shapes Mukhopadhyay et al., 2002) the same skeletal deformation that we consider here actually acts to inhibit budding, in their case to form echinocytes. In these studies skeletal energy was not estimated as we do by simulating actual cell deformations. This view is entirely consistent with ours, except that here we explicitly account for release of skeletal energy as membrane/skeleton separation occurs.
Note the geometric factors ℓ(t) and h(t), where ℓ is the current radius of the membrane separation and h is the height of the bulging membrane; here t may be viewed as a "time like parameter" marking the progression of the process (see Figure 8Aii). Vesiculation occurs soon after h → ℓ as discussed below. Note also, that there is the "work of adhesion" that must be done to de-bond the area of the membrane still attached to the skeleton; this called γ (per unit area); this has been discussed at length in the context of splenic vesiculation (Zhu et al., 2017). It is important to note that the membrane considered in Figure 8A is not open, i.e., finite, but is considered part of the overall cell's closed membrane. As budding occurs, additional membrane is drawn in as the vesicle's membrane area increases; this process requires the additional work of membrane/skeleton separation, γ , described below. Now, for perspective, note that the bending energy of a spherical vesicle is simply 8πκ b , independent of its radius; κ b is the membrane's bending modulus. On the other hand, the elastic energy stored in a circular membrane/skeleton patch of radius ℓ is FIGURE 8 | Schematic of a budding vesicle. In (Aii) we illustrate a pre-curvature as discussed in the text. From (Aiii to Av) we have progressing membrane bulging leading to vesiculation. (B,C) depict a model expansion analyzed below. Vesiculation begins within a damaged membrane/skeleton area of diameter 2ℓ 0 . Driven by release of stored skeletal energy a bulge forms and possibly grows as in (B,C) where the energy pathway is given in Figure 11. The final stage of pinching-off and vesicle release in illustrated in Figure 12. πℓ 2 ǫ 0 . This would seem to set a lower limit to the size of a vesicle that may form, i.e., we need We next discuss the skeletal energy, i.e., ǫ 0 , and then the debonding energy γ .

Skeleton Stored Energy
As shown in Figure 9, we consider a single JC consisting of six spectrin dimers (Sp), numbered 1 to 6. In its undeformed state the length of each Sp is Ł 0 so that the area covered by this hexagon is We now consider a case in which the JC is stretched in x direction by a factor of λ 1 and in y direction by a factor of λ 2 (λ 1 ≥ λ 2 ), the corresponding area deformation is ψ = λ 1 λ 2 and the shear deformation is χ = λ 1 /λ 2 . In the deformed state the lengths of the Sp become L i (i = 1, · · · , 6).
Invoking the Worm Like Chain (WLC) model (Weiner, 1983;Zhu et al., 2007) the strain energy stored in the i-th Sp is Zhu et al. (2007) where k B is the Boltzmann constant. T is the temperature. L c and p are the contour length and persistence length of Sp, respectively. In the following example we choose L c = 144 nm, L 0 = 40 nm, and p = 0.8 nm. The strain energy density is then calculated as 6 i=1 φ i /(ψA 0 ). We note additionally that the erythrocyte skeleton is in a deformed state while attached to the cell's membrane "at rest, " i.e., without imposed deformations. This is revealed clearly in the observations of Vertessy and Steck (1989) and Shen et al. (1986) FIGURE 9 | Simplified sketch of a junctional complex consisting of six Sps. The initial length of the JC spokes are taken as L 0 = 40 nm here. who observed skeletal shrinkage in Triton X-100 membranes. Such collapse of released membrane skeletons was reported by Levin and Korenstein (1991) and Truvia et al. (1998) in their studies of membrane fluctuations. In such studies skeletal contractions were at least on the order of a factor of 2.

Skeletal Energy Density vs. Deformation Modes, ǫ 0
In Figure 10 we present contours of ǫ 0 as it depends of the deformation parameters χ = λ 1 /λ 2 (representing shear deformation) and ψ = λ 1 λ 2 (representing area deformation), where λ 1 and λ 2 are principal in-plane stretches and by definition λ 1 > λ 2 . We note that, within the range of deformations explored, we find We note this range would lead using Equation (1) to minimal separated membrane patch and hence vescile sizes in the range This is indeed a relevant, albeit lower, size range as noted, for example, by measured size ranges for vesciles formed during splenic flow (Bosman et al., 2012). It is also noteworthy that this areal size scale is comparable to that covered by a JC unit as depicted in Figure 9. This provides a basis for the comment made above in section 1.1 that vesicles with diameters less than, say, 40-50 µm are unlikely. For example, a membrane patch with a area of πr 2 m would produce a vesicle of diameter r v ∼ 0.5r m as follows from the membrane's area incompressibility, i.e. πr 2 m = 4πr 2 v for a spherical vesicle.
FIGURE 10 | Stored skeletal energy density in units of Jm −2 vs. deformation mode. The heavy line indicates a possible path of reduced skeletal area and shear deformation. The points labeled "a" and "b" correspond to maximum energy densities achieved in simulations that are described in section 4. The broken line arrow suggests an energy path aiding promoting vesiculation.

Work of Membrane/Skeleton Adhesion, γ
Estimates for the work required to separate the membrane from the skeleton have been reviewed and analyzed recently by Zhu et al. (2017). Their analysis was based on experimental data measured for experiments on the formation of tethers. The results are somewhat tentative, yet would estimate γ to be of order γ ∼ O(1 × 10 −4 )Jm −2 . This order is noteworthy when compared to our simulated range of the stored elastic energy in the skeleton, viz. ǫ 0 ∼ O(3 − 8 × 10 −4 )Jm −2 . This implies that if (ǫ 0 − γ ) < 0 this combination represents a driving force for vesiculation, although this assumes that as an area element of membrane/skeleton separates all of the stored elastic energy is released.
The heavy line in Figure 10 drawn at slope β represents a possible deformation path taken by the skeleton during its release from the membrane during a vesiculation process as discussed below.

Vesiculation Simulations
Here we consider a simple scenario of spherical cap-like bud emerging from an initially flat membrane/skeleton section. The cap initiates from a circular patch, or raft, of radius l 0 as shown in Figures 8(Aii,iii) and expands as a growing spherical cap with height h(t) and a fixed cord 2ℓ 0 . As the cap expands in area it draws in lipid bilayer and requires additional area of skeleton/membrane interface to separate.
Let the cap be of a sphere of radius r(t). The emerging cap's area is A cap (t) = π(h 2 + ℓ 2 0 ) and its radius is such that 2r(t) = (h 2 +ℓ 2 0 )/h; of course this means that as h → ∞, 2r → h. Also let ℓ be defined as the radius of a circular raft of membrane/skeleton with area of the spherical cap bud of Figures 8B,C. Hence the change in the area of skeleton/membrane area separated upon the cap's growth is δA ske = δ(πℓ 2 ). But the change in the cap's area is δA cap = δ{π(h 2 + ℓ 2 0 )}. Since the membrane is incompressible, A ske = A cap , and thus ℓδℓ = hδh and ∂ℓ/∂h = h/ℓ. Now the change in energy stored in the skeleton that is currently in the cap is and this leads, upon differentiation with respect to h noting ∂ℓ/∂h = h/ℓ, to δǫ ǫ = −2πhǫ 0 δh. Next we note that as δA ske = 2πhδh from above that the work expended in membrane/skeleton separation is δǫ γ = 2πγ hδh.
The energy of bending in the cap is given as with κ taken as the mean curvature given by With Equation (6) we have (7) Note that as h → ∞, ǫ κ = 8πκ b as expected. This yields upon differentiation Taken together we have with, as above, ζ = ǫ 0 − γ . To resolve this we may write h = αℓ 0 and obtain Upon reintegrating Equation (10) we obtain
Here expansion continues as h increases along a downward slope of a ǫ vs. α(h) curve. At some stage the "bridge, " ℓ 0 , may itself retreat as the bud becomes more compliant with increasing radius, r = (h 2 + ℓ 2 0 )/2h. As the angle θ defined in Figure 12 decreases, the bending stresses in at the junction of the bud and bilayer main body increase and eventually rupture the bilayer. Spectrin may be retained within the vesicle, e.g., bound to Hb (Snyder et al., 1985), and/or released into the medium so as to remove skeletal proteins from the cell in rough proportion to lost membrane area in the vesicle as suggested by Ciana et al. (2017b). Resealing on the membrane is favored, of course, as is reattachment of the skeleton to the membrane.
As for estimates of vesicle size, we may proceed as follows. Consider as a first example Figure 11A; here the activation barrier is crossed at α ∼ 1.2(ζ = 4) or α ∼ 2(ζ = 1) with ℓ = 2 andκ = 1.5. Now analysis of Figure 12 shows that and hence if θ reduces so that θ → θ crit , 2α crit Let us take, to explore the numerology, θ crit = π/8; this yields α crit ≈ 5. At this stage we would have h v ≈ 5 × 2 × 10 −8 m = 100nm. On the other hand ifl = 4, as in Figure 11C, we would find h v ∼ 200nm. As a purely geometric criterion, Equation (14) would effectively define the length scale ℓ 0 as the mediator of vesicle formation (at least vesicle size). Yet the vanishing of G act that depends also on κ b , ǫ 0 and γ as well as on ℓ 0 all play vital roles. However, it may be judged, that as long as G act → 0, ℓ 0 may be judged as a (the) prime determinant of vesiculation. Of course, the very factors that influence γ as well as ǫ 0 , and possibly even κ b are involved in causing an ℓ 0 in the first place. With the numerology used here we find that predominantly vesicles are expected in the size range 100nm h v 200nm; however, vesicles as small as h v ∼ 40 − 50nm are possible if sufficient driving energy, ǫ 0 is available.
Finally we discuss the release of a severed bud, i.e., its closing to a spherical vesicle. The question is: once free, does the bud completely fold or revert to a flat disc with an exposed edge? This question is depicted in Figure 13. The severed bud ("B") is depicted in Figure 13b and it may revert to a planar disc of area "A" as in Figure 13a or close to a vesicle ("V") as in Figure 13c. The bud and the flat circular patch, however, possess an "edge, " Ŵ, with energy γ s measured in units Jm −1 and which provides a driving force to close the bud or, if insufficient, fold a flat patch as in "A".  The issue of the folding of an initially flat patch was considered by Helfrich (1974), and later by Fromherz (1983), and here we list a key result as phrased by Hu et al. (2012) who developed an approach based on these results for estimating values for the key parameters. In Equation 15 E is the energy along the path given by curve A → V in Figure 13.
Here ξ is a parameter that measures the progress along path A → B → V where ξ = 0 represents the initial flat patch of membrane of area A and ξ = 1 a full spherical vesicle of radius R (whose area is equal to A). The modulus κ b has already been defined and κ g is the so-called Gaussian modulus (Helfrich, 1974;Fromherz, 1983;Hu et al., 2012) (κ g ≈ −κ b ). Now as noted by Hu et al. (2012) -and is readily verified -the path A → B → V displays an activation barrier if η < 2 at ξ act = 1 − (η/2) 2 ; this amounts to Ẽ act = (1 − η/2) 2 . With this, our question is: as per Figure 13b, does the B revert to A or proceed to vesicle V after the driving energy supplied by, inter alia, ǫ 0 and G mix has been expended?
Now if at the point of bud severing α were in the range, say, α ≥ 2 andl ≥ 2 we would judge that η 2.6. In such cases no activation barrier would exist (from a flat disc) to form a closed vesicle. If, on the other hand,l were modest in size, sayl ∼ 1 and α 2 at budding then η 1.22 and a barrier would exist to form a closed vesicle. In such cases one might suspect that severed bud would revert back to a flat disc, i.e., a "membrane fragment" as are often found.

DISCUSSION OF VESICULATION SIMULATION RESULTS OF SECTION 3
The results shown in Figures 11A-F provide an initial assessment of the key parameters, κ b , ζ = ǫ 0 − γ , and ℓ 0 . These are, in turn, affected by metabolic and biochemical influences as discussed above. Now given the magnitudes of the activation barrier, G act as indicated in Figure 11D, we expect vesiculation to essentially require G act → 0; we thereby take this as a possible criterion for vesicle release. For example, Figures 11A,D clearly indicate the effects of high values of κ b and/or small initial membrane/skeleton separation ℓ 0 in inhibiting vesiculation; the effect of ℓ 0 is noteworthy. This perhaps obvious result is nontrivial as it suggests that the initial "triggering event" may have a mediating effect on the eventual size of the vesicle released, as opposed to the details of progression of skeleton/membrane separation and ongoing membrane bulging. Figures 11B,C,F support this view as they show that, even at modest-to-high values of κ b and low values of ζ , G act may effectively vanish if ℓ 0 is at least as large as a typical JC complex, viz. if ℓ 0 ∼ 40 − 50nm. Next, we note that the magnitude of κ b (Harmandaris and Deserno, 2000;Hu et al., 2012), the bending modulus, clearly plays a role in setting the height of G act , and even determining if G act ≥ 0. In general, κ b depends on membrane thickness and also on the bilayer membrane's lipid composition, including the lipid composition asymmetry. Here we introduce Figure 14, taken from Zwaal et al. (1997) that helps explain these issues. Yet, the prospect for vesiculation is also determined by the relative magnitudes of the terms ǫ 0 and γ given in Equation (11) and (12) via ζ = ǫ 0 − γ . We continue to discuss the prospects for G act → 0 below. For now we note that Equation (12) reveals an activation barrier at α 2 = √ 80/µ − 1 with µ =l 2ζ /κ. This shows, after substitution into Equation (12), that if µ/80 < 1 an activation barrier exists with G act > 0; conversely, we have that FIGURE 14 | Regulation and physiology of membrane phospholipid asymmetry and vesiculation, taken from Zwaal et al. (1997). The effects of Ca 2+ uptake on phospholipid randomization and on calpain activation that facilitates blebbing and release of PS exposing vesicles are addressed. We use the figure to also illustrate the geometry of a budding vesicle and to highlight how large curvatures develop and lead to pinching off of the vesicle bud. Note also that the large tensile strains in the membrane may increase the kinetics of randomization that initiates before vesiculation (Zwaal et al., 1997).
implies that, once initiated, a bud will spontaneously grow. Now the roles of calpain and Ca 2+ uptake in vesiculation have been already discussed and are referred to in Figure 14. We take the depiction of budding as shown there to represent a stage in any of our Figure 11 but along a downward slope of a ǫ vs. α(h) curve; note the curvatures at such stages, most particularly at the sites where the bud meets the membrane body. As such stages we envision that phospholipid asymmetry has been lost (or partially lost) and that the membrane is under considerable tension at what is to be the "pinching off " sites located as just referred to above. This is addressed in more detail below as we suggest how vesiculation is completed, but after a relevant diversion concerning the spectrin skeleton.

Disposition of Skeletal Proteins During Vesiculation
A question arises suggested by the Zwaal-Schroit depiction of Figure 14 as well as by Figure 12 below, viz. what happens to the spectrin contained in the released skeleton? It has been reported that shed vesicles (MV's) are deficient in skeletal proteins, viz. spectrin (Lutz et al., 1977;Dumaswala and Greenwalt, 1984;Knowles et al., 1997;Willekens et al., 2003a,b;Bosman et al., 2012) and this has suggested that spectrin may remain in the cells. Ciana et al. (2017a,b) have, however, recently shown that the cells that shed vesicles, in vivo, lose spectrin in rough proportion to their membrane area; no data was presented on the spectrin content of the shed vesicles, however. The fact is that the number of vesicles found in venous blood is orders of magnitude too low to account for the cell membrane lost. This means that vesicle membrane is recycled without entering general circulation, i.e., presumably within the organ where the vesicle forms. Yet Snyder (Snyder et al., 1985) had already shown that vesicles shed from "young cells" (i.e., cells so classified by their low density) contained "Hb-spectrin complexes" suggesting that spectrin was also shed into vesicles in at least some amounts. Hence we speculate that spectrin may indeed be lost to either the medium and/or retained within shed vesicles along with Hb as is repeatedly reported. In fact, examining either Figures 14, 12 or Figure 2 of Alaarg et al. (2013) suggests that, since the pinching off process produces very large local deformations, the skeleton would bed severed and retained within the vesicle or the medium of the organ producing vesicle. As many reports of spectrin deficient vesicles are those produced in-vitro we discuss this further and will suggest, in the Discussion, that our oscillatory shear methods be utilized to study the character of vesicles produced via our splenic-like deformations. This logically leads us to more detailed consideration of γ and ǫ 0 .

Biochemical Energetics
As noted in section 3 most prominently, and elsewhere, attachment between the skeleton and membrane can be reduced by a number of biochemical and metabolic causes. Hence we realize that although imposed deformations that increase ǫ 0 will stimulate vesiculation, e.g., during splenic flow or under our conditions of oscillatory shear flow, even modest states of deformation are sufficient to induce the process. An intriguing possibility is that the free energy released by losing phospholipid concentration asymmetry may provide some of the required energy associated with ǫ 0 . This remains to be explored, but here we make a simple estimate of such a contribution to test the idea's viability. The fact is that significant metabolic free energy is expended in creating lipid concentration asymmetry (Zwaal et al., 1997;Fadeel and Xue, 2009;Bevers and Williamson, 2016) and we wonder what portion of it may be available to perform the work of vesiculation.
We take a simple view and consider the primary lipids that "mix" in a final state of symmetry from an initial state of asymmetry between the outer and inner bilayer leaflets to be grouped into 2 groups. Let group #1 consist of sphingomyelin and phosphstidylcholine (initially located on the outer leaflet) and group #2 contain phosphatidylethanolamine, phosphatidylserine, and phosphatidylnositol (initially on the inner leaflet). We will call the molecular fraction of either group "x" since we assume a perfect solution within the leaflets.
We use the well known relation for the entropy of mixing for a perfect solution, viz.
where, in context, n is the total number of lipids in a leaflet, k is Boltzman's constant, and x is either group #1's or #2's molecular fraction. We use Equation (18) to compute and note that the free energy change, just associated with perfect mixing, is G a→s = −T S mix ; we will take T = 300 K. We take for n, the number of lipids per leaflet, 1.4 × 10 18 m −2 ≤ n ≤ 1.5 × 10 18 m −2 (White and King, 1985;Giang and Schick, 2014) (this being computed from the average area per lipid, a ∼ 0.68 − 0.70nm 2 ) and kT ∼ 4 × 10 −21 J. For × we take x ∼ 0.74 (Zwaal et al., 1997;Fadeel and Xue, 2009). These numbers yield G mix ∼ −(6.7 − 7.2) × 10 −4 Jm −2 , which is in the range of ǫ 0 , a primary driving force for vesiculation (see Figure 10). Even if a portion of this stored free energy were made available it could drive vesiculation. For perspective, we recall that the standard free energy of hydrolysis of ATP is in the range G • ATP ∼ −(11 − 13)kcal/mole ∼ −(7.5 − 9) × 10 −20 J/ATP (Alberts et al., 2002). Hence one ATP per 10 2 − 10 3 lipids could fuel vesiculation. We note Beleznay et al. (1994) and Zachowski (1993) report that transport of lipids such as PS consumes about one ATP per lipid, i.e., an order of magnitude larger than we estimate for G mix per lipid -as such lipid transport must have a substantial activation barrier, this seems quite reasonable. Interactions and association of PS with spectrin have been documented (An et al., 2004a,b;Grzybek et al., 2006) and those may be disrupted during membrane/skeletal separation, promoting translocation.
Thus for the above reasons we expect reductions in γ and increases in ℓ 0 are prime contributors to, and even mediators of, vesiculation. A simple view would be that the skeletal energy density, ǫ 0 , would be augmented by − G a→s as a prime driver of vesiculation. As imposed deformations, i.e., ǫ 0 , have the prospect of promoting the process as we explore next regarding imposed oscillatory shear flows in the specific context of the above vesiculation model prospects. Finally we address the "pinching off " process that releases the vesicle.

Pre-curvature and Biochemically Induced Curvature
Returning to Figure 11, we recognize that when G act > 0 we are typically in scenarios where α 1 − 1.5 yet if we began the process at, say α 1.5 in most cases G act would vanish. We thereby ask about the effects of pre-curvature, or spontaneous curvature, as may arise, for example, by the expected trends to randomize the lipid concentration between the inner and outer leaflets (Zwaal et al., 1997). As noted Sheetz and Singer (1974) anionic molecules such as PS that are known to partition to the outer leaflet with vesiculation will indeed tend to cause crenation, i.e., outward budding (see also Deuticke, 1968). Just above we have shown that the thermodynamic driving force for this can be, indeed, substantial. A simple analysis of such induced curvature might go as follows. Simply consider Equation (7), the bending energy, and equate it to a portion of the free energy of mixing G mix ; call this portion G mix . This yields or (21) For some numerology, take α = 1, κ b = 10 −19 J and recall the definition ofl to obtain 5 × 10 −4 ,l = 2 2.2 × 10 −4 ,l = 3 1.25 × 10 −4 ,l = 4, etc. (22) Hence we may envision that in all the cases shown in Figure 11, taken as examples, G act may be driven to near zero. Now the reader will notice that as we take ℓ 0 larger and larger, the required G mix reduces, eventually as G mix ∼ 1/α 2 . These arguments also demonstrate that the effects of the cell membrane's natural curvature are quite modest due to the much larger curvatures involved in creating a budding vesicle; that is, the cell's surface appears to be relatively flat compared to the emerging vesicle.

Vesiculation in Oscillatory Flow: Final Analysis
Taken together, the results herein and those of sections 2 and 3 would support a scenario such as this: imagine a deformation state such as represented by points "a" or "b" of Figure 10 that can arise from the simulations we described. Now imagine an augmentation of the associated ǫ 0 ∼ 5 × 10 −4 Jm −2 of order, say − G a→s ∼ O(2 − 3 × 10 −4 Jm −2 ). This places the total prime vesiculation driving force at ∼ 7 − 8 × 10 −2 Jm −2 , that according to the result scenario's of Figure 11 could very probably produce vesiculation.
A further assessment of vesiculation prospects with the examples of Figures 4-6 could go as follows. Consider the vesiculation simulations as, say, in Figure 11. From Figure 11C, for example, looking at the contours of energy density we would conclude that if there were regions-call them aged regionswhere skeletal/membrane disruptions had occurred over areas of, say,l ≥ 3 the probability of vesiculation at regions whereζ ≥ 2 would be quite high; at regions whereζ ≥ 4 vesiculation would be judged to be nearly assured. This follows since in this case G act → 0. Hence, in Figure 6 we would conclude that at any region displaying a red color we expect vesiculation to occur in cells containing aged areas wherel 3 − 4. In Figure 4 which also displays deformation states quite similar to those developed in splenic flow (Zhu et al., 2017) we would likewise focus on those regions in the lower-right foreground. Note, here we take γ = 1 × 10 −4 Jm −2 (Zhu et al., 2017) and that the red areas displayζ ≥ 3. On the other hand, the case shown in Figure 3B displays somewhat reduced levels of stored energy, i.e.,ζ 3 and hence we judge that vesiculation may require larger aged areas involving, sayl 4 − 5.
Similar case studies can be made and from such we could well conclude that during splenic flow -as closely replicated by our oscillatory flow -vesiculation is highly probable providing aged regions exist covering nearly the area of a single JC unit (i.e., corresponding tol ∼ 3 − 4). This would surely classify the spleen as an effective filter as expected -it might also support the notion that the majority of such self-protection vesicles are produced in vivo in such flow conditions, e.g., in the spleen and elsewhere where such deformations are imposed.

Shear Induced Erythrocyte Membrane Trauma
Given the flexibility and inherent time and rate dependence of RBC deformation, it is unsurprising that cell damage induced by imposed shear displays a complex phenomenology. For perspective vis-à-vis oscillatory shear, as we have presented it, we briefly review some background. We begin by noting that the physiological range of shear stress imposed upon the RBC is often put in the range ofτ ∼ 5 − 10Pa (Meram et al., 2013); it is hence noteworthy that our forecasted average shear stresses presented herein for oscillatory flow are in that range. However, as explained early on by Leverett et al. (1972), the question of membrane rupture depends not only on stress level but also on the time of exposure suggesting perhaps "creep-like" damage mechanisms are operative. Moreover, stressing rate is important (Li et al., 2013) and in the case of cyclically imposed stresses, both frequency and amplitude are important (Watanabe et al., 2006;Wantanbe et al., 2007;Hashimoto, 2014;McNamee et al., 2016); this suggests that even wave form for periodic stressing would likely be still another factor. In fact, patterns using stress (flow) bursts coupled to pauses (i.e., hold times) have been employed (McNamee et al., 2016;Horobib et al., 2017). Nevertheless, when stresses are imposed for long times ( 60 min) a critical shear stress to induce cell lysis in the range τ crit ∼ 150Pa has been quoted (Leverett et al., 1972); other reports place this higher at τ crit ∼ 250Pa (Sutera and Mehrjardi, 1975). Yet-and this is important!-in a study involving exposure to shear stresses τ ∼ 10Pa, but for 2 h, it was found that such exposed cells were recognized and sequestered by rabbit spleens (Sandza et al., 1974). That is, such cells were sequestered and removed during splenic flow. Here we may speculate that although such exposure was sublethal (i.e., did not cause bursting, cell distortion, or other residual damage), it did induce cell removal processes, perhaps vesiculation. Note that the time to induce possible vesiculation in Sandza et al. (1974) was unknown, but the stress level at τ ∼ 10Pa is just at the levels we forecast for oscillatory flow as in the case of Figure 5. Finally, we note that human splenic flow passage times are < 1s (MacDonald et al., 1987;Zhu et al., 2017). Given that multiple passages may be involved, we take the stress protocol to involve multiple passes with durations of, say 0.02-1 s; these to employ shear stress levels of τ ∼ 5 − 10Pa. In oscillatory flow we would suggest frequencies in the range ν ∼ 1 − 50Hz.
We also note that our simulations methods can be used to analyze the deformations experienced by RBC's subject to all the flow patterns imposed in the studies cited above, and can thereby be used to assess potential damage, lethal or sublethal, that may be caused. Moreover, it is possible to include interactions between cells and/or cells and structural features such a channel walls so that more quantitative assessments can be made of possible influences they may have on cell response.

SUPPORTING EXPERIMENTAL EVIDENCE
To confirm the effects of oscillatory flow as we have described them and to assure that our proposed imposed deformations do not cause lethal damage to a significant population of cells we used flow cytometry to explore the process of vesiculation and possible fragmentation following typical programs of shear. We imposed both modest and somewhat severe deformations to more completely explore the extent of possible cell damage via fragmentation and vesiculation. We suggest below, however, that in exploring aging induced vesiculation only the more modest deformations we have described be imposed. For example, an oscillatory frequency of 10 Hz appears to be reasonable since it coincides with the typical duration (0.1 sec) of slit passage in splenic flows (MacDonald et al., 1987). The shear stress on the cell depends on both the applied shear rate and the viscosity of the surrounding fluid. To match the typical shear deformation of the skeleton in splenic flow (the area deformation is negligibly small) this stress should be 5-10 Pa as described in section 4.5 and by the various plots of deformation in oscillatory flow we have included.

Blood Samples and Reagents
Adsol (Fenwal Laboratories, Deerfield, IL) preserved and leukoreduced red cell units were purchased from the San Diego Blood Bank (San Diego, CAI). Glycophorin A-PE, Annexin V FITC, and annexin V-binding buffer concentrate were purchased from BD Pharmingen (San Jose, CA). Flow Cytometry Absolute Counting Standard microbeads (6-7 µm) were purchased from Bangs Laboratories, Inc. (Fishers, IN). Megamix, a blend of monodisperse fluorescent beads of three diameters (0.5, 0.9, and 3 µm) was purchased from BioCytex (American Diagnostic, Hauppauge, NY, United States).

Oscillatory Shear
The steady and oscillatory shear conditions were applied on a stress-controlled shear rheometer (T.A. Instruments, model AR-G2), with a cone-plate of 60 mm diameter and 1 • of angle, and with a plate-plate geometry of 60 mm diameter and a gap of 100 µm. The shear flows are imposed with a duration of 1-2 s, and repeated for multiple cycles, and carried out at 37 • C.

Isolation and Purification of Red Cell Vesicles From Shear RBC Suspension
Isolation of vesicles from the RBCs suspensions were completed by high-speed centrifugation as previously described (Kriebardis et al., 2008). Briefly, 1-mL aliquots were removed and centrifuged at 2,000 g at 4 • C. The supernatant was centrifuged once again to ensure the absence of any RBCs and immediately filtered through sterile 0.8 mm pore size syringe-driven nitrocellulose filter units (Millipore, CA, United States). The supernatant was ultra-centrifuged at 37,000 g at 4 • C for 1 hour, and the pellet of vesicles was resuspended in PBS and ultra-centrifuged twice under the same conditions.

Quantitation of Micro-Vesicles Using Flow Cytometry
Red cell vesicle analysis and enumeration was obtained using a FACS Aria flow cytometer (BD Biosciences, San Jose, CA, United States). To confirm the formation of vesicles after imposed shear, the side scatter and forward scatter events from size calibration beads (0.5, 0.9, 3, and 7.6µm,) was compared to the side scatter and forward scatter events of a small volume of supernatant and erythrocytes of blood centrifugated at 2,000 g at 4 • C ( Figure 15A). These results were used to determine the resolution of the instrument, and to confirm the presence of vesicles and smaller cell fragments as the amount of imposed shear deformation increased ( Figure 15A). The size range was confirmed and corrected by spiking a blood sample with calibration beads to establish the effect of the calibration beads' higher index of refraction on size relative to membrane vesicles (lower size threshold than membrane vesicles) (Yuana et al., 2011). Purified vesicles from RBC suspensions were labeled with glycophorin A-PE and annexin V FITC at ice temperature and protected from light. Red cell vesicles were discriminated by size and further defined as CD235a-PE and Annexin V positive event. Vesicle counts were calculated by adding a predetermined number of calibration beads to each blood sample, and from the nominal number of beads added per volume of sample a total number of vesicles was calculated.

Results
Two-color flow cytometric analysis, employing a combination of antiglycophorin A and Annexin-V-FITC, has demonstrated that shear deformations induce vesiculation as theoretically forecasted herein. Flow cytometry results are shown in Figure 15, where size is one of the key factors defining these events. Flow cytometry allowed us to obtain information on the morphology (size and granularity) of vesicles formed due to imposed shear, as evidenced by the forward scatter and side scatter of vesicles compared to calibration beads of appropriate diameters ( Figure 15A). Forward scatter properties were used as an estimate of size, and side scatter properties were used as an estimate of particle shape. Forward and side scatter events from size calibration beads were used to resolve the instrument sensitivity and detection range. Figure 15A shows that the subpopulation of vesicles increased with shear deformation as forecasted. Examination of the top-left scatter plot in Figure 15A indicates that 0.2 µm is the lower limit of detection for the beads. Vesicles are defined as membranous vesicles of arbitrary size from 0.2 to 1.0 µm. In contrast to RBCs, vesicles derived from RBCs may expose negatively charged phospholipids toward their surface. Living cells, however, employ energy-dependent mechanisms to actively shift negatively charged phospholipids to the inner membrane leaflet as discussed in section 4.2. Flow cytometric observation allows for the determination of the fraction of vesicles that bind annexin V, which suggests exposure of negatively charged phospholipids. Figure 15B shows that the fractions of glycophorin A positive and exposed PS vesicles (observations in the upper right quadrant) increased with shear deformation. Fluorescence events from antiglycophorin A and Annexin-V-FITC showed that the erythrocyte vesicles that expose PS, also had a high expression of glycophorin A.
Moreover, these results prove that oscillatory shear flows do not damage a significant portion of the RBC's. As an example, the results in Figure 15 indicate that after a protocol of 100 cycles, of shear rate σ = 200s −1 , at 100 Hz we found the number of vesicles to be O(10 4 ). As this protocol corresponds, based on the deformations experienced by the cells, to approximately 10,000 in-vivo splenic slit passages for each cell we would expect that given our 10 7 cells this would have resulted in a number of vesicles of O(4.8 − 6.5 × 10 9 ) in-vivo. In this we assumed that in-vivo a vesicle is produced at a rate of approximately 2-2.71 per day per cell (Willekens et al., 2003b;Ciana et al., 2017b). Hence as this expected in-vivo vesicle count is clearly much larger than our result of O(10 4 ) we conclude our methods do not induce spurious cell damage or unexpected vesiculation. In terms of in-vivo splenic passages the deformation histories experienced by cells in this protocol would have occurred over about 250 days (if cell deformations outside of the spleen are not considered), or more than twice of a cell's average lifetime. Thus, a reason for the difference is that cells would not be expected to metabolically age during as they would in-vivo in the quite short durations of our tests. Additionally, we imposed somewhat more severe deformations (shear rate amplitude 2,000 s −1 at 50 Hz and with 100 cycles with durations of 2 s) to further explore the extent of possible cell damage. The basic hemologic parameters demonstrate that the majority of the cells (> 95%) survive the treatment.
Vesicle counts found per microliter of plasma were as follows: 217 for our control with no-shear and 986, 23,034, and 46,026 Top left panel shows light scatter parameters of the mixture of standard calibration beads (7.6, 3, 0.9, and 0.5 µm). Other panels show light scatter parameters of the blood samples before mechanical shear (no-shear), and after being subjected to different frequencies, cycles, and degrees of oscillatory shear. Inspection of the scatter from the calibrated beads and the blood subjected to mechanical shear indicates the increasing presence of vesicles below 1 µm in size. Direct size comparison cannot be made between beads and vesicles, since beads have a higher index of refraction, and therefore lower size threshold, than vesicles (Lacroix et al., 2010;Yuana et al., 2011). (B) Flow cytometry analysis of RBC vesicles extracted from blood before oscillatory shear (no-shear), and after being subjected to different frequencies and number of cycles of oscillatory shear. Only events positive to anti-CD235a and Annexin V were defined as RBC vesicles (Kriebardis et al., 2008). The number of calibration beads was counted to determine the absolute number of RBC vesicles. The number of cells was approximately 10 7 .
for the cases of 1 cycle at 10 Hz, 10 cycles at 100 Hz and 100 cycles at 100 Hz, all with a shear rate of σ = 200s −1 , respectively. The number of 217 vesicles per µl of plasma is consistent with and helps confirm previous such findings, e.g., those of Willekens et al. (2003b) who reported the range 61-308 with an average of 169 vesicles per µl of plasma.
Measurements have shown that the volume rate of blood flowing through 100gms of spleen is approximately 170 ml/min/100 gms (Oguro et al., 1993). Assuming the spleen weighs, on average, 150gms and the body contains 4.5l of blood we estimate that the splenic passage rate of a typical erythrocyte is about 80 per day. However, not all these passages involve passages through the venous slits of the red pulp; in fact depending on mammalian species and degree of health versus diseased states, up to 90% of the blood flow may bypass the filtration beds of the red pulp (Schmidt et al., 1993;Cesta, 2006). Thus, the precise rate of erythrocyte passage through the venous slits is difficult to specify and here we make the assumption of 50%. Hence we arrive at an estimate for, discussion sake, of ∼ 40 slit transits per cell per day. We thereby envision they undergo a pulsedlike deformation with a frequency of say ν f ∼ O(0.0005)Hz. Hence we suggest that our lower frequency test protocols are more appropriate for studies of age induced vesiculation. Now the most severe deformation occurs within a time scale of 0.1 s. A realistic deformation protocol would thus be to subject cells to one oscillation cycle, pause for 2160s, start a second cycle, and repeat that for ∼ 40 × 120 ∼ 4, 800 times (as a cell does in typical circulation). As this may not be practical, we believe the approach we use, i.e., subject a population of cells to high frequency oscillatory flow for 1-2 s and repeat for O(10 − 100) cycles is an attractive, viable alternative; this is consistent with physiological deformation time scales as discussed in section 4.5. Indeed, we also suggest that lower shear rates such as σ = 200 − 500s −1 be explored as we have done.

CONCLUDING DISCUSSION
To illustrate the versatility of our proposed methods, we revisit the discussion of section 4.1 concerning the disposition of skeletal proteins, viz. spectrin, during vesiculation. As noted there, Ciana et al. (2017b,a) have shown that cells that shed vesicles in-vivo lose spectrin in proportion to their loss of membrane; this makes sense given the cell's need to maintain a functional mechanical membrane/skeleton structure. Now consider the in-vitro experimental findings of e.g., Knowles et al. (1997) who induced vesiculation via micro-pipette aspiration. Indeed, (i) skeleton/membrane separation, (ii) membrane vesicle ejection, and (iii) skeleton retraction back into the cell, were all observed. Hence, these vesicles were indeed deficient in spectrin and we speculate also in anchoring proteins. This process was successfully simulated and explained by simulation with an earlier version of our computational models (Peng et al., 2010). In short, what we found in Peng et al. (2010) was that the deformations and their time scales cause very significant skeletal restructuring at regions of eventual vesiculation. This was due to the viscous drag of anchoring proteins (viz. band 3 and glycophorin C/JC) and results in very large (∼ 80-90%) reductions in areal densities of skeleton anchoring sites. This, combined with the large disassociation stresses at the aspiration's tip, causes the membrane to "tear away, " or rather "pull away, " from the cell and skeleton where the latter retracts back into the cell. But we later showed that such a scenario does not occur during splenic flow due to the very different time scales and hence reductions in anchoring density must occur otherwise (Zhu et al., 2017) . That is, during in-vivo splenic flow, there is not time for mechanically driven reductions in anchoring density. This explains the important role our analysis parameter ℓ 0 plays in deciding the progress of vesiculation. Indeed, if ℓ 0 were very large, say ℓ 0 40 − 80nm, as it may well be in in-vitro experiments, vesiculation would readily occur via a different pathway as suggested by Ciana et al. (2017b,a). Since our methods are designed to reproduce both the deformation histories, with the appropriate time scales, as occur in splenic flow, we suggest that the character of vesicles produced in oscillatory flow be examined with various degrees of membrane/skeleton disruption.
We add that as far as the effects of viscosity and viscosity ratios is concerned, that the hydrodynamic load on the cell depends on the frequency of oscillatory shear flow as well as the amplitude of the shear stress, η 1 σ ; here η 1 is the viscosity of the outside medium and σ is the amplitude of the shear flow. Numerical algorithms work more efficiently when the ratio η 2 /η 1 = 1 with η 2 being the cell viscosity. For that reason simulations were often performed with η 2 /η 1 = 1 with the understanding that has little effect as long as the amplitude of the shear stress, η 1 σ , is unchanged. This correlation has been confirmed by simulations with viscosity ratios in the range 1 ≤ η 2 /η 1 ≤ 7. This observation is most useful when exploring a large parameter space as mentioned above. Our previous work provides much further detail as to the broader range of numerical methods (see e.g., Peng et al., 2010Peng et al., , 2011.
The focus herein has been on oscillatory shear flow states that produce splenic-like conditions of deformation, stress distributions around the cell, and their respective time scales. It would be of interest, however, to more systematically explore a wider range of scenarios involving various conditions of frequency, shear rate amplitude, and cell/medium viscosity. In this manner we may explore the potential effects of time scales that may allow, for example, skeletal restructuring to occur. This is a large, yet interesting, topic of future research.
The scenarios presented just above, along with all the results and analysis of section 4 support the view that vesiculation can be driven by either metabolically, or mechanically (i.e., imposed deformation), induced energetics. In either case a reduction in skeletal/bilayer binding, (i.e., aging induced reductions in γ ), although not strictly required, are strongly causative influences that promote vesiculation. This view is consistent with our earlier conclusions in Zhu et al. (2017). Purely mechanically driven vesiculation is unlikely under physiological conditions since such large separating stresses are unlikely to develop. Exceptions are, for example, flow in aspiration as demonstrated by Peng et al. (2010), where we have noted the time scales as well as the severe deformations were sufficient to effectively reduce γ by dynamically restructuring the skeleton that reduces the areal density of attachment points. However, such restructuring does not occur during splenic flow or our oscillatory flow introduced herein (Zhu et al., 2017) due to the inherently shorter time scales in these types of flow. Indeed, our modeling described in section 3 defined a key parameter, ℓ 0 , that represented a membrane/skeleton area damaged by metabolic and/or biochemical means that played a strong role, if not mediated, the vesiculation process. The importance of such damaged membrane zones has been recently emphasized (Leal et al., 2018).
This view is consistent with the fact that although a given RBC transits the spleen on the order of 40 times per day, such cells only produce about 2.71 vesicles per day (Ciana et al., 2017b). Hence vesiculation, even under conditions of large deformations, should be viewed as a rather improbable event! An efficient way to generate vesicles systematically in laboratory conditions is to process vast number of cells simultaneously, which happens to be the advantage of our oscillatory flow device.
We have noted that the prime driving force for vesiculation in our model is ǫ 0 , the skeleton's elastic energy, may be augmented by − G a→s , or at least by some part of it. This free energy is released upon a transition to a more symmetric lipid leaflet composition. But on the one hand, the kinetics for lipid flipping is not normally large and thus times scales are a question. Yet, as inner leaflet lipids like PS are exposed on vesicles, it appears that flipping does indeed occur during vesiculation; it is possible that the normally slow kinetics of lipid flipping is increased by membrane deformation. But this will alter the membrane energetics due to, e.g., non-local bending energy as well as through the pre-curvature energy term (Miao et al., 1994;Waugh, 1996;Lim et al., 2002;Mukhopadhyay et al., 2002;Bozic, et al.); specifically, it relates to, for example, how one assesses the terms involving A 0 in the cited works. Hence, this question awaits future detailed study.
It was noted that RBCs of those affected by Scott's syndrome display a decreased rate of vesiculation (Bevers and Williamson, 2016;Bevers et al., 2017). Indeed, Rosing et al. (1985) reported that the syndrome originates from a defective scramblase activity that suppresses PS exposure on the outer membrane leaflet. Hence these observations may support the idea that a loss in lipid asymmetry may play a direct role in promoting vesiculation as suggested above and in section 4.2.
We add that our approach allows for a thorough study of the effects of not only cytosol viscosity but also medium, or plasma, viscosity that will effect the forces and deformations experienced by the cell during a given imposed flow (Williams and Morris, 1980;Chien, 1987;Tuvia et al., 1997;Pozrikidis, 2005;Freund et al., 2012;Freund, 2013;Zhu et al., 2017). Effects of altered viscosity, including plasma viscosity, have been associated with a variety of dysfunctions and disease (e.g., Késmárky et al., 2008;Sapmaz et al., 2011;Toprak et al., 2012). In general, with increased plasma viscosity deformations, and the associated ǫ 0 , will be augmented and hence so will the prospect for vesiculation. Hence our approach offers novel prospects for pursuing such lines of future inquiry.
We herein introduced a novel theoretical/simulation approach to subject cells to tailored shear deformations that, mimic, and expand on the types of deformations that promote vesiculation (Zhu et al., 2017). This clearly suggests an experimental plan, based on these findings, that would subject erythrocytes to biochemical stress including, e.g., oxidative, nitrite (NOS), and Ca 2+ uptake stress, inter alia, so as to induce "aging damage, " e.g., disruptions to the cell's skeleton/bilayer membrane and to follow the resulting course of promoted vesiculation as influenced by carefully tailored modes of deformation with varying shapes and intensity. The findings presented in sections 3.1 and 4 demonstrate that such an approach is viable and would clearly provide novel insights into the aging/vesiculation process.
Finally we add that, aside from the existence of nonbiological walls that typically are part of flow chambers, our methods involve no artificial structures and naturally subject large numbers of cells to our tailored deformations. We demonstrate how splenic flow induced deformations can be produced and thereby provide conditions to study events such as age induced vesiculation. Our methods may also provide a valuable compliment observations made within micro fluidic devices (e.g., Deplaine et al., 2011;Picot et al., 2015) provided these can be fabricated with splenic-like slits with dimensions ≤ 1.5µm and through which red cells can actually flow through and not simply become "jammed"; such behavior involves deformations and time scales quite unlike splenic flow and is hence nonrepresentative. Moreover, although of use in visualizing individual passage through even sub-micron artificial slits (Gambhire et al., 2017), micro-fluidic chambers do not naturally deal with the very large number of cells required for vital statistical studies; recall that, in-vivo, vesiculation is a quite rare event, occurring roughly 1 in every 15 splenic passages.

AUTHOR CONTRIBUTIONS
RA and QZ wrote the manuscript and performed the analysis and simulations. PC contributed to the experimental portions of the research.