Modelling of Red Blood Cell Morphological and Deformability Changes during In-Vitro Storage

Featured Application: Red blood cell (RBC) storage lesion is a critical issue facing transfusion treatments, and signiﬁcant changes in RBC morphology and deformability are observed due to the storage lesion. RBC membrane biomechanics, and cytoskeletal structural integrity are explored numerically in this study. Abstract: Storage lesion is a critical issue facing transfusion treatments, and it adversely a ﬀ ects the quality and viability of stored red blood cells (RBCs). RBC deformability is a key indicator of cell health. Deformability measurements of each RBC unit are a key challenge in transfusion medicine research and clinical haematology. In this paper, a numerical study, inspired from the previous research for RBC deformability and morphology predictions, is conducted for the ﬁrst time, to investigate the deformability and morphology characteristics of RBCs undergoing storage lesion. This study investigates the evolution of the cell shape factor, elongation index and membrane spicule details, where applicable, of discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies during 42 days of in-vitro storage at 4 ◦ C in saline-adenine-glucose-mannitol (SAGM). Computer simulations were performed to investigate the inﬂuence of storage lesion-induced membrane structural defects on cell deformability and its recoverability during optical tweezers stretching deformations. The predicted morphology and deformability indicate decreasing quality and viability of stored RBCs undergoing storage lesion. The loss of membrane structural integrity due to the storage lesion further degrades the cell deformability and recoverability during mechanical deformations. This numerical approach provides a potential framework to study the RBC deformation characteristics under varying pathophysiological conditions for better diagnostics and treatments. such as number of spicules, average spicule height, average spicule surface area and average curvature, were observed between the storage lesion-produced and echinocytogenic shape-transforming conditions-produced echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies. Furthermore, the results predict the adverse inﬂuence of membrane structural defects due to storage lesion on cell deformability and recoverability during mechanical deformations, and indicate further degradation of cell viability.


Introduction
A red blood cell (RBC) transfusion is a critical and lifesaving medical procedure for anaemia caused by disease leading to inadequate bone marrow production or increased RBC breakdown, or by blood loss due to trauma or surgery [1][2][3][4]. Therefore, in-vitro storage is important to fulfil both constant and sudden demands of RBCs for transfusion. According to the current blood bank procedures, RBCs can be stored up to 42 days at 2-4 • C in a preservative solution [1,5], and studies on alternative preservative solutions, rejuvenation procedures, anaerobic storage conditions and temperature effects [2,6,7] are being conducted to ensure the maximum storage time for a healthy and safe RBC supply. The stored RBCs are suitable for transfusion if the mean haemolysis is below 0.8% and more than 75% RBCs survive the first 24 h of post-transfusion period [5,[8][9][10][11]. However, concerns have been raised about the quality of stored RBCs as several transfusion complications [12] are observed with respect to the shelf life of the stored RBCs. During the storage, RBCs undergo several biochemical, structural and functional changes, usually known as storage lesion [2,4,5,[12][13][14][15][16][17][18][19][20]. These changes include slowed metabolism with a decrease of adenosine triphosphate (ATP) concentration, oxidative damage with structural changes to band 3 protein and lipid peroxidation, functional loss of cation pumps and consequent loss of intracellular potassium and accumulation of sodium, and apoptotic changes of membrane phospholipids and membrane vesiculation [3,19,[21][22][23]. Most of these changes that occur due to decreased metabolism such as depletion of ATP and 2,3-diphosphoglycerate (DPG) and cation imbalance are reversible under in vivo conditions. However, protein and lipid modifications, morphology changes, and membrane vesiculation are not reversible under in vivo conditions [15,22]. These irreversible modifications can affect the in vivo cell survival leading to adverse post-transfusion outcomes. Therefore, it is important to measure the quality of each RBC unit prior to use, and the RBC deformability is a key indicator of RBC quality and post-transfusion viability [8,[24][25][26]. RBC deformability is primarily influenced by mechanical and geometrical factors of the cell such as cell surface area and volume, elasticity and viscosity of the cell membrane, and volume and viscosity of the cytosol [27][28][29][30][31][32][33]. Therefore, the changes in the cell membrane structure and its mechanical properties adversely affect the cell deformability, and less deformable RBCs can obstruct capillaries and require significantly higher transit time to navigate through the microvasculature leading to decreased levels of oxygen delivery to organs [34][35][36]. In addition, less deformable RBCs are promptly removed from the circulation at the spleen [14,37].
These investigations improve the basic understanding of RBC behaviour under healthy and pathophysiological conditions, in order to be integrated into clinical applications. However, a more detailed understanding on RBC ageing during storage, ageing-influenced changes in RBC metabolism, structural and functional aspects, is required to improve the quality, efficacy and safety of stored RBCs and thereby to prevent or reduce the adverse post-transfusion outcomes. Numerical modelling is a desirable approach to investigate the fundamental mechanics of cell deformability over the experimental investigations, which can be challenging and costly [67,68,73]. Accurate numerical modelling combined with experimental measurements can provide valuable insight on RBC quality and viability. This paper proposes a numerical approach inspired from the previous research for RBC morphology predictions [74][75][76] and for analysis of RBC deformations [39,58] to investigate, for the first time, the deformability and morphology changes of RBCs subjected to storage lesion. A three-dimensional coarse-grained (CG)-RBC membrane model based on the bilayer-couple model (BCM) [76] was employed to predict RBC echinocytosis during 42 days of in-vitro storage at 4 • C in saline-adenine-glucose-mannitol (SAGM) preservative solution. The evolution of RBC morphology, its membrane shear modulus and bending modulus, and cell surface area and volume based on analogous experimental measurements [8,47], was incorporated into the model to better predict the deformability and morphological changes of stored RBCs. The deformability and morphology characteristics: cell shape factor, elongation index and membrane spicule details where applicable, of stored RBCs having discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies, were studied during the 42 days of in-vitro storage, and compared against analogous morphologies obtained via echinocytogenic shape-transforming conditions. The shape-transforming conditions have the potential to recover the morphological and deformability changes of stored RBCs, and to lower the influence of storage lesion. Therefore, the present study investigates the morphology and deformability characteristics of RBCs during in-vitro storage and compares them against echinocytic conditions. In addition, the influence of storage lesion-induced membrane structural defects on cell deformability and its recoverability of sphero-echinocyte morphology under optical tweezers stretching deformation was investigated. The following subsections briefly outline the evolution of RBC morphology and associated CG-RBC membrane model parameters for the prediction of echinocytosis, and discuss the deformability and morphology changes of RBCs during in-vitro storage. This study provides valuable predictions of the influence of storage lesion on RBC quality and viability, and is currently being executed to investigate the recoverability of deformability and morphology changes of RBCs undergoing storage lesion.

Methods
The storage lesion-induced changes in RBC morphology, membrane mechanical properties and cell geometry were implemented on a CG-RBC membrane model [76] to predict the analogous echinocytosis process during in-vitro storage. In the first step, the correlations of storage time-dependent evolution of cell morphology index, shear modulus and bending modulus of the membrane, and cell surface area and volume were determined based on analogous experimental observations (i.e., the data presented in Figures 3, 4 and 5 of Kozlova et al. [8], Figure 4a of Park et al. [47] and based on the data extracted from 3D confocal microscopy imaging of RBCs during in-vitro storage by D. Santangelo (D. Santangelo, personal communication, 10 October 2018)). In the second step, the corresponding RBC morphologies were numerically predicted through the CG-RBC membrane model to investigate the deformability and morphology characteristics of RBCs during in-vitro storage.

Evolution of RBC Morphology during In-Vitro Storage
In this section, the storage time-dependent evolution of RBC morphology index (MI) is determined for the numerical modelling of RBC echinocytosis during in-vitro storage. MI is based on the Bessis's classification [75,77], and is 0 for the biconcave discocyte; +1 for the echinocyte I (a flattened or elongated cell with several undulations around its rim); +2 for the echinocyte II (a flattened elliptical cell with rounded spicules distributed more or less uniformly over its surface); +3 for the echinocyte III (an ovoid or spherical cell with sharper and more numerous spicules distributed evenly over its surface); and +4 for the sphero-echinocyte (a spherical cell with small sharp projections still attached to its surface). Figure 1 summarizes the evolution of RBC morphology during in-vitro storage and was determined based on the experimental observations by Kozlova et al. [8]. It can be observed (refer to Figures 4 and 5 of Kozlova et al. [8]) that five main forms of RBC morphologies arise during in-vitro storage, which are discocyte, echinocyte, sphero-echinocyte, dacrocyte and swelled cells, which are in the form of spheres and ovals. However, only the evolution of discocyte, echinocyte and sphero-echinocyte cells were considered for the present study as these are the main forms of RBC morphologies observed during storage [3,20,21]. The lack of information on the evolution of dacrocytes and swelled cells during in-vitro storage makes it difficult to identify their contribution to the RBC morphology and deformability changes during storage and their recoverability. It can be observed that the major fraction of RBCs up to 20 days of in-vitro storage is composed of discocytes, whereas the major fraction of RBCs from 26 days of storage is composed of sphero-echinocytes. Discocyte RBCs follow the morphology transformation of discocyte → echinocyte I → echinocyte II → echinocyte III → sphero-echinocyte during days 20 to 25, and this period can be considered as the transition period of RBC morphology. The periods of echinocyte I, echinocyte II and echinocyte III morphologies are estimated based on the assumption that the transformation of discocyte to sphero-echinocyte takes place at a steady rate during the transition period. The corresponding RBC morphology in terms of MI as a function of in-vitro storage duration (t), is presented in Figure 1.

Evolution of Coarse-Grained (CG)-RBC Membrane Model Parameters during In-Vitro Storage
The CG-RBC membrane model is composed of a triangular surface of N t triangles representing N V cytoskeletal actin junctions and N S spectrin links. The model predicts the minimum free energy state of the RBC membrane, which is the equilibrium RBC morphology at the given reference conditions (e.g., cell surface area (A 0 ), cell volume (V 0 ), reference bilayer-leaflet-area-difference (∆A 0 ) and reference total-membrane-curvature (C 0 )). The CG-RBC membrane model is briefly reviewed in Appendix A: CG-RBC Membrane Model, whereas detailed description and discussion on the model development and shape predictions of the discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte, are available in Geekiyanage et al. [76]. Geekiyanage et al. [76] presents the numerical predictions of stomatocyte-discocyte-echinocyte morphology transformations of a healthy RBC under stomatocytogenic and echinocytogenic shape-transforming conditions, and the corresponding model parameters and constraint coefficients for discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies are summarised in Appendix B: Prediction of Discocyte, Echinocyte I, Echinocyte II, Echinocyte III and Sphero-Echinocyte RBC Morphologies.
The CG-RBC membrane model requires appropriate adjustments to compensate for the evolution of RBCs during in-vitro storage. Therefore, the numerical evolution of the membrane shear modulus (µ(t)), bending modulus (κ(t)), cell surface area (A 0 (t)), cell volume (V 0 (t)), reference bilayer-leaflet-area-difference (∆A 0 (t)) and reference total-membrane-curvature (C 0 (t)), is determined as a function of in-vitro storage time. The CG-RBC membrane model parameters mentioned in Appendix B: Prediction of Discocyte, Echinocyte I, Echinocyte II, Echinocyte III and Sphero-Echinocyte RBC Morphologies, were adopted as the analogous RBC properties at day 0 of in-vitro storage, and the corresponding membrane shear modulus (µ(t 0 )), bending modulus (κ(t 0 )), cell surface area (A 0 (t 0 )), and cell volume (V 0 (t 0 )) were established accordingly.
The evolution of the CG-RBC membrane µ(t) and κ(t) during in-vitro storage is presented in the Figure 2a and 2b, respectively. The estimation of µ(t) is based on the experimentally observed evolution of RBC membrane Young's modulus by Kozlova et al. [8], and it is assumed that the evolution of Young's modulus during in-vitro storage is comparable to that of the membrane shear modulus. Young's modulus is one of the key indicators of cell membrane mechanical integrity, and thereby the cell stiffness [8]. However, the storage lesion affects the Young's modulus. It can be observed (refer to Figure 3 of Kozlova et al. [8]) that the Young's modulus of RBC membranes is approximately constant up to 20 days of in-vitro storage, and increases rapidly during days 20 to 25, where it remains approximately constant during the period of days 26 to 42. This behaviour is comparable with the evolution of RBC morphology during in-vitro storage. The evolution of µ(t) was determined in agreement with the outline of the evolution of RBC membrane Young's modulus, and is presented in Figure 2a, where the evolution of µ(t) during storage was assumed to be directly proportional to that of the Young's modulus (i.e., Young's modulus ∝ µ). It can be observed that µ(t) increases rapidly during the transition period, and reaches a value which is nearly 2.6 times higher than µ(t 0 ). These observations agree with available experimental results [18,38]. The µ(t) corresponding to echinocyte I, echinocyte II and echinocyte III during the transition period was determined through linear interpolation of µ(t). Similarly, the evolution of κ(t) was determined based on the measurements of RBC membrane mechanics by Park et al. [47]. It can be observed (refer to Figure 4a of Park et al. [47]) that the spherocyte morphology has the highest membrane bending modulus, whereas the echinocyte has an intermediate membrane bending modulus. The evolution of κ(t) was determined in agreement with these observations such that the measured membrane bending modulus of the discocyte [47] is comparable with a κ(t) of up to 20 days of in-vitro storage (i.e., bending modulus of the discocyte ∝ κ(t); 0 ≤ t <20) and the measured membrane bending modulus of the spherocyte [47] is comparable with a κ(t) from day 26 of in-vitro storage (i.e., bending modulus of the spherocyte ∝ κ(t); 26 ≤ t ≤ 42). The κ(t) corresponding to echinocyte I, echinocyte II and echinocyte III during the transition period was determined through linear interpolation of κ(t). The evolution of κ(t) during in-vitro storage is presented in Figure 2b.
The cell surface area is an important property of RBCs as the exchange of oxygen and carbon dioxide takes place at the cell surface. A higher cell surface area is required for improved gas exchange in the lungs and body tissues [78]. However, RBCs lose their cell membrane during cell aging, and the cell becomes smaller and denser during storage [38]. The evolution of A 0 (t) and V 0 (t) was determined based on the data extracted from 3D confocal microscopy imaging of RBCs during in-vitro storage (D. Santangelo, personal communication, 10 October 2018). The RBC surface area and volume measurements for the discocyte, echinocyte I, echinocyte II and echinocyte III cell morphologies corresponding to three donors during 42 days of in-vitro storage in SAGM additive solution are presented in Appendix C: Measurements of Cell Surface Area and Volume of RBC during In-Vitro Storage-Experimental Observations. Figure 3 represents the general evolution of cell surface area and volume during in-vitro storage. The loss of membrane surface area and volume results in an increased concentration of micro-particles in the extracellular additive solution, and a notable increase of micro-particle concentration is observed from day 28 onwards [9,20,79]. In agreement with these observations, the evolution of cell surface area and cell volume was considered for two in-vitro storage durations such that 0 ≤ t < 28 and 28 ≤ t ≤ 42, and the expressions as in Equations (1) and (2) respectively, were established to reflect this evolution. The coefficients of Equations (1) and (2) are estimated based on the linear regression relationships between storage duration (i.e., dependent variable) and cell surface area and volume (i.e., independent variable) for the two in-vitro storage durations (i.e., 0 ≤ t < 28 and 28 ≤ t ≤ 42). The local surface area of the RBC membrane was assumed to follow the evolution of A 0 (t), and the reference area of the k th triangle of the CG-RBC membrane model (A k,0 (t)) was estimated according to Equation (3). A k,0 (t 0 ) is the reference area of the k th triangle at day 0, and is adopted from Appendix B: Prediction of Discocyte, Echinocyte I, Echinocyte II, Echinocyte III and Sphero-Echinocyte RBC Morphologies.
The evolution of ∆A 0 (t) and C 0 (t) during in-vitro storage was assumed to be equivalent to that of the CG-RBC membrane model-predicted discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies (refer to Appendix B: Prediction of Discocyte, Echinocyte I, Echinocyte II, Echinocyte III and Sphero-Echinocyte RBC Morphologies), and the corresponding ∆A 0 (t) and C 0 (t) are summarised in Figure 4. It can be observed that ∆A 0 (t) and C 0 (t) increase quickly during days 20 to 25 in agreement with the discocyte → echinocyte I → echinocyte II → echinocyte III → sphero-echinocyte morphology transformation.  These storage-time-dependent model parameters were developed based on the limited available experimental data from RBC concentrates that are collected, processed and stored under varied conditions. However, considering these effects as well as any influence of the donor characteristics and heterogeneity of the RBC population can strengthen the computational approach and improve the numerical predictions. In addition, the availability of a fundamental computational approach facilitates more specific experimental data to be obtained at a larger scale to carry out more statistical and stochastic modelling.

CG-RBC Membrane Model-Predicted RBC Morphologies during In-Vitro Storage
The in-vitro storage effects were implemented on a CG-RBC membrane model to predict the analogous RBC time-dependent echinocytosis process. Identical model parameters and constraint coefficients as in Appendix B: Prediction of Discocyte, Echinocyte I, Echinocyte II, Echinocyte III and Sphero-Echinocyte RBC Morphologies, except for µ(t), κ(t), A 0 (t), A k,0 (t), V 0 (t), ∆A 0 (t) and C 0 (t), were adopted to predict the RBC morphology at t (Days) = 0, 7, 14, 21, 23, 25, 28, 35 and 42. The RBC morphology at the above time points during in-vitro storage was predicted as the stable minimum energy configuration at equivalent reference conditions, and is presented in Table 1. The numerically predicted RBC morphologies at day t qualitatively agree well with both the corresponding morphology characteristics given by Bessis's classification [75,77] and with the corresponding experimental observations [54,[74][75][76] for discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies. It can be observed that the thickness of the biconcavity of the discocyte morphology increases with t, and the cell diameter reduces in response to the increased reduced cell volume at t.
In addition, the discocyte cell shape at t = 7 (Days) and t = 14 (Days) indicate uneven cell thickness on its rim, and the echinocyte I at t = 21 (Days) has uneven undulations around its rim. The loss of RBC membrane surface area during in-vitro storage induces increased levels of stress on cytoskeletal spectrin links, and this may affect the equilibrium cell state and morphology at t. Moreover, the cytoskeletal reference state is predicted at an equivalent membrane shear modulus (µ(t 0 )), membrane bending modulus (κ(t 0 )), membrane surface area (A 0 (t 0 )) and cell volume (V 0 (t 0 )) of a healthy RBC at t = 0. Therefore, the cytoskeleton is subjected to increased lateral compression due to the loss of membrane surface area during storage. These changes during in-vitro storage may have influenced the predicted RBC morphology at t. Table 1. The CG-RBC membrane model-predicted morphology of in-vitro stored RBC. Shear modulus (µ(t)), bending modulus (κ(t)), cell surface area (A 0 (t)) and cell volume (V 0 (t)) during in-vitro storage are presented in the form of multiplications of µ(t 0 ), κ(t 0 ), A 0 (t 0 ), V 0 (t 0 ), respectively, which represent the corresponding parameters at t = 0 (Days) (refer to Figures 2 and 3). ∆A 0 (t) and C 0 (t) are the reference bilayer-leaflet-area-difference and reference total-membrane-curvature, respectively (refer to Geekiyanage et al. [76]). Cross-sectional view of discocyte morphology is presented.

Evolution of RBC Morphology Characteristics during In-Vitro Storage
The functionality of a RBC is linked with its cell geometry. The evolution of cell dimensions (i.e., cell length, cell diameter and cell thickness) with respect to in-vitro storage effects can be described through the cell shape factor (SF) [76]. The cell is assumed to be a rigid body having uniform density, and the centre of mass of the whole cell and its three-principal axes of inertia are determined for each numerically predicted cell shape at t (Days) = 0, 7, 14, 21, 23, 25, 28, 35 and 42. H 1 , H 2 and H 3 are defined as the distance between the furthest vertex points on the RBC membrane surface along the three-principal axes of inertia of the cell such that H 1 ≤ H 2 ≤ H 3 . SF is given by [76,80], and indicates the sphericity of the cell. The cell becomes more and more spherical as SF reaches the value 1, and becomes a more and more flattened disc as SF reaches 0. SF at t (Days) = 0, 7, 14, 21, 23, 25, 28, 35 and 42 was estimated and is presented in Figure 5. The SF of stored RBCs gradually increases at 0 ≤ t < 20 and at 26 ≤ t ≤ 42, whereas a sudden increment of SF was observed at 20 ≤ t < 26. The evolution of SF is comparable to the evolution of cell morphology, µ(t) and κ(t) during storage (refer to Figures 1 and 2), and can be attributed primarily to the evolution of cell morphology as µ(t) and κ(t) are influenced by the cell morphology as well [8,47]. It is important to distinguish the morphology characteristics of stored RBCs against other pathophysiological conditions where similar morphologies are observed, for better diagnostic purposes. Each shape-transformation scenario might affect the cell details differently. Therefore, Figure 5 also compares the morphology characteristics of storage lesion-produced discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte RBC morphologies against that of echinocytogenic shape-transforming conditions (e.g., amphipaths, extracellular ionic strength and pH) produced RBC morphologies. The CG-RBC membrane model [76] is capable of predicting RBC morphologies under varying shape-transformation scenarios, and Geekiyanage et al. [76] discussed the morphology characteristics of discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies of a healthy RBC at t = 0 (Days), which were obtained via echinocytogenic shape-transforming conditions. The echinocytogenic shape-transforming conditions transform the biconcave discocyte shape into echinocyte forms at constant cell surface area and volume [48][49][50][51][52][53][54]. SF values of discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte RBC morphologies obtained via the echinocytic shape-transforming condition (from Geekiyanage et al. [76]) are presented in Figure 5 at t where analogous morphology was observed during storage.
The dependency of cell SF on its morphology was confirmed by the similar behaviour of SF values for the two considered shape-transforming scenarios; in-vitro storage and echinocytogenic condition. Overall, SF of stored RBCs is higher than that of echinocytogenic conditions-produced RBCs, which shows the influence of the loss of membrane surface area and cell volume during storage. The cell surface area and volume of the stored RBC gradually decreases during storage, whereas the morphological changes in the presence of echinocytogenic conditions take place at constant cell surface area and volume. However, the SF of echinocyte I morphology at t = 21 (Days) is comparable to that of echinocyte I morphology obtained by echinocytogenic conditions. This discrepancy can be attributed to the cell details of the echinocyte I morphology obtained at these two shape-transforming scenarios. The echinocyte I at t = 21 (Days) has uneven undulation around its rim, whereas the echinocytogenic conditions produced echinocyte I has even undulations around its rim. This can influence SF estimation. In general, the evolution of SF of RBCs under both these shape-transforming scenarios indicate the rise of cell sphericity leading to lowered cell deformability [76]. These results agree with the experimentally observed deformability changes of RBCs [26,78], and indicate the potential of SF as a quality measure of stored RBCs.
At latter stages of in-vitro storage, the spicules on sphero-echinocyte RBC morphology form membrane budding and detach from the cell such that the spherocyte is formed. The progression of sphero-echinocyte to spherocyte is characterized by shortened spicules such that complete spheres are formed [77,81]. The SF of a spherocyte is 1.0 and indicates the null deformability of the RBC, and it is important to predict and understand the progression of spherocyte formation during storage. The evolution of membrane spicules of echinocyte II, echinocyte III and sphero-echinocyte was therefore investigated to study the progression of a RBC towards a spherocyte during storage. Number of spicules (N Spicules ), average spicule height (H Spicules ), average spicule surface area (A Spicules ) and average spicule curvature (C Spicules ) of echinocyte II, echinocyte III and sphero-echinocyte at 22 ≤ t ≤ 42 were studied and compared against that of the echinocytogenic conditions-produced echinocyte forms (refer to Figure 6). H Spicules of the spiculated cells was estimated as follows; where r kk,max is the distance to the membrane vertex on the spicule that has the maximum C V,i from the centroid of the cell, whereas r kk,min is the distance to the membrane vertex on the spicule that has the minimum C V,i from the centroid of the cell. C V,i is the membrane curvature at i th membrane vertex (refer to Equation (5)), and the maximum C V,i is given by the membrane vertex at the peak of the spicule, whereas the minimum C V,i is given by the membrane vertex at the foot of the spicule. A Spicules of the spiculated cells was estimated according to Equation (6); where N V is the number of membrane vertices, N V−T,i is the number of neighbouring triangles to which the i th vertex is attached to, and M jj and ∆A jj are the membrane curvature and the membrane surface area at the triangle-pair that shares the j j th triangle edge, respectively. In general, C V,i > 0 occurs on the spicule surface, whereas C V,i ≤ 0 occurs at the valleys between the membrane spicules. Similarly, C Spicules of the spiculated cells was estimated according to Equation (7).
It can be observed (refer to Figure 6) that the number of spicules on echinocyte II and echinocyte III are similar for both shape-transforming scenarios. However, the number of spicules on the sphero-echinocyte at t (Days) = 28, 35 and 42 is less than the number of spicules on the sphero-echinocyte achieved under echinocytogenic conditions, and increases with the storage duration. This indicates the higher sensitivity of sphero-echinocyte morphology on µ(t), κ(t), A 0 (t), A k,0 (t), and V 0 (t) than the echinocyte II and echinocyte III morphologies. In particular, the loss of membrane surface area and cell volume during storage constraints the formation of spicules on the cell surface. In general, H Spicules and A Spicules of stored RBCs are lower and decrease with the storage duration. This demonstrates the smaller size of membrane spicules of stored RBCs than that of echinocytogenic shape-transforming conditions-produced RBC shapes, and the spicule size decreases with the storage duration. Similarly, C Spicules of stored RBCs is higher and increases with the storage duration. The gradual decrement of spicule size and increment of spicule curvature of stored RBCs together indicate the progression of membrane vesiculation during in-vitro storage, which leads the cell towards spherocyte formation. In particular, the rapid increment of C Spicules of sphero-echinocytes during 26 ≤ t ≤ 42 confirms the increased potential for spherocyte formation during this period. ATP depletion during storage adversely affects the structural integrity of the cytoskeleton, and the membrane vesiculation contributes to removal of defective membrane patches from the cell [26,82] such that the spherocyte is formed. The present CG-RBC membrane model treats the RBC membrane as a single constituent, and the probable detachment of the lipid-bilayer from the underlying cytoskeleton for membrane vesiculation cannot be observed. However, the model can be extended to consider the remodelling of the cell membrane under membrane vesiculation.

Evolution of Deformation Behaviour of RBCs during In-Vitro Storage
The evolution of deformation behaviour of stored RBCs was investigated for CG-RBC membrane model-predicted RBC echinocytosis. RBC deformability is a measure of cell viability, and gradually decreases with storage duration [29,83,84]. The optical tweezers stretching forces were applied on the discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies at t (Days) = 0, 7, 14, 21, 23, 25, 28, 35 and 42, and the evolution of cell deformability in terms of elongation index (EI) during storage was determined. EI was estimated at the equilibrium cell stretched state based on axial (D A ) and transverse (D T ) diameters of the stretched cell [30,62], such that EI = (D A − D T )/(D A + D T ). The methods of implementation of optical tweezers stretching forces on the RBC and the estimation of EI are discussed in Appendix D: Implementation of Optical Tweezers Stretching. The effect on EI of the discocyte and echinocyte cells at echinocytogenic conditions undergoing optical tweezers stretching at varying membrane shear modulus, spectrin link extensibility, membrane bending modulus and d RBC membrane-bead contact diameter, is discussed in Geekiyanage et al. [85]. The evolution of EI of stored RBCs at the total stretching forces of F ext = 60, 100 and 200 pN are presented in Figure 7. The EI of each RBC morphology increases with the increase of F ext , and reduces with the storage duration. The EI of day 42 sphero-echinocyte is significantly lower than that of the day 0 discocyte; 54.44% at F ext = 60 pN, 43.39% at F ext = 100 pN, and 29.19% at F ext = 200 pN. Therefore, the stored RBCs demonstrate significantly lower deformability especially at low F ext . A rapid decrement in EI can be observed during the transition period (20 ≤ t < 26), and demonstrates the reducing cell deformability with respect to its morphology transformation. These results indicate that the RBC deformability remained unaltered during 0 ≤ t < 20, and RBCs become difficult to elongate when stored for longer periods and when morphology transformations occur. These observations are consistent with available analogous observations as well [9,17,26,38]. The deformation behaviour of RBCs undergoing optical tweezers stretching is primarily due to the cytoskeletal stretching characteristics [86,87], and therefore the evolution of EI is primarily lead by cytoskeletal stretching characteristics (i.e., membrane shear modulus (µ) and maximum extensibility of spectrin links (x 0 )). The gap between EI values at t decreases with increasing F ext , and can be attributed to the influence of x 0 such that the stretchability of the cell is inhibited at higher F ext . The EI of echinocyte I morphology at t = 21 (Days) indicates a lower EI in comparison to the predicted general EI evolution during storage. The definition of EI becomes less suitable for irregular shapes [26], and therefore this observation may result from the uneven undulations on the rim of echinocyte I cells at t = 21 (Days). A slight increment in EI was observed between similar morphologies with respect to the storage duration; for the discocyte morphology at t (Days) = 0, 7 and 14, and also for the sphero-echinocyte morphology at t (Days) = 28, 35 and 42.
The evolution of µ(t) and κ(t) during storage is comparable with the morphology transition (refer to Section 2). Therefore, the discocyte morphology at 0 ≤ t < 20 has a similar µ(t) and κ(t), and the sphero-echinocyte morphology at 26 ≤ t < 42 has a similar µ(t) and κ(t) as well. However, the cell surface area and volume gradually decrease throughout the storage. Therefore, the observed increment of EI for the similar morphology can be attributed to the loss of cell surface area and volume. The results indicate that the stored RBCs have acceptable deformability during 0 ≤ t < 24 up to echinocyte II morphology. However, echinocyte III and sphero-echinocyte morphologies during 24 ≤ t ≤ 42 have degraded deformability, and contribute to poor quality and viability of RBCs at prolonged storage.

Influence of Cytoskeletal Actin Junctional Defects during In-Vitro Storage on RBC Morphology and Deformability
The ATP depletion during in-vitro RBC storage can lead to cytoskeletal protein defects such that the connections between the cytoskeletal spectrin network and the lipid-bilayer are disrupted [32,68,82]. RBCs with defective vertical interactions between the lipid-bilayer and the cytoskeleton can lose the unsupported membrane area contributing to membrane vesiculation. In this section, the influence on the equilibrium cell shape and its deformation behaviour is investigated for stored RBCs with defects in the membrane proteins that cause disruption in the vertical interactions between the cytoskeleton and the lipid-bilayer. The protein defects in the vertical interaction are presented in the CG-RBC membrane model by broken connections of spectrin links at randomly selected actin junctional complexes. Loss of vertical connectivity at the actin junctional complexes (A Vertical ) from 0% to 25% was introduced to the CG-RBC membrane model such that the spectrin links attached to these defective actin junctional complexes do not contribute to the cytoskeletal stretching response. Therefore, the elastic energy of the cytoskeleton (E Stretching ) for a cell membrane with loss of vertical connectivity between cytoskeleton and lipid-bilayer, is given by where S A Vertical indicates the spectrin links attached to defective actin junctional complexes. k B is the Boltzmann constant, T is the absolute temperature, l max is the maximum extension of j th link, p is the persistence length and k p is the power function coefficient. x j is defined as x j = l j /l max for the j th link having l j length, and m is an exponent such that m = 2. The cytoskeletal spectrin network of A Vertical = 0, 1, 10 and 25 (%) of defective vertical interactions at randomly selected actin junctional complexes is presented in Figure 8. The equilibrium cell shapes and membrane spicules become irregular with the severity of membrane defects (i.e., increase of A Vertical ). The equilibrium cell shapes of day 42 sphero-echinocytes with A Vertical = 0% and 25% subjected to optical tweezers stretching forces of F ext = 0 and 200 pN are presented in Figure 9. Spicules with sharp projections can be observed (refer to Figure 9) in the presence of A Vertical , especially if the actin junctional defects and spicule tips coincide together. At higher A Vertical , the unsupported lipid-bilayer forms sharp projections as it can easily contribute to ∆A 0 (t) and C 0 (t) in the absence of a cytoskeletal stretching response. The potential to form sharp projections increases with the increase of surface area of the lipid-bilayer segment that is unsupported by the cytoskeleton. Therefore, the loss of cytoskeleton and lipid-bilayer vertical interactions at adjacent actin junctional complexes promotes the formation of sharp projections on the cell surface, and can promote vesiculation at the tips of these sharp projections. This observation is in agreement with the experimental and computational predictions, which show that vesiculation does not require fragmentation of the cytoskeleton, although it can be enhanced by the uncoupling of the lipid-bilayer and the cytoskeleton [32,82]. Therefore, the sphero-echinocyte cells at the latter stages of storage have increased potential for vesiculation due to the higher number of membrane spicules and depleted ATP levels.
In addition, the deformation behaviour of stored RBCs at t (Days) = 0, 7, 14, 21, 23, 25, 28, 35 and 42 was also investigated for RBCs with a defective cytoskeleton of A Vertical = 0, 1, 10 and 25 (%). The evolution of change in the axial (D A ) and transverse (D T ) diameters of stored RBCs with varying A Vertical with respect to its analogous recovered cell state (i.e., recovered RBC shape at t and A Vertical = 0%) at F ext = 60.0 and 200.0 pN is presented in Figure 10.
It can be observed (refer to Figure 10) that the change in D A and D T cell diameters in comparison to the recovered cell state increases with A Vertical , and indicates a similar form to that of the evolution of cell EI during storage. As expected, the change in D A and D T at equivalent A Vertical increases with F ext as a response to the increased external stretching forces. The increasing change in D A and D T with A Vertical demonstrates the increased cell fluidity in the presence of decreased cytoskeleton and lipid-bilayer interactions. The cells showed a higher D A than D T to accommodate the applied stretching forces while maintaining cell reference conditions, and stretch along the axial direction of stretching. The observed lower change in D A of echinocyte I morphology at t = 21 days can be attributed to the slight deviations of CG-RBC membrane model-predicted echinocyte I morphology at t = 21 days from an ideal echinocyte I morphology. However, it is difficult to distinguish any significant difference in the evolution of D A of stored RBCs in the presence of high F ext (= 200.0 pN) and A Vertical (= 25%). This demonstrates the insignificant dependency of changes in cell morphology, µ(t), κ(t), A 0 (t), A k,0 (t), and V 0 (t) on RBC stretching response at very high F ext and A Vertical , and the probability of cell rupture is very high for these RBCs. The present CG-RBC membrane model does not consider cell rupture, but can be extended to investigate the mechanisms of cell rupture of defective cells subjected to mechanical deformation.

Influence of Cytoskeletal Actin Junctional Defects during In-Vitro Storage on Cell Recovery following Deformations
It is also important to investigate the influence of storage-induced membrane defects on cell recovery subjected to mechanical deformation. Among those factors that influence the cell stretching response (e.g., spectrin network, intracellular viscosity and cholesterol content), the cytoskeletal spectrin network is the primary contributor for cell stretching response, and differences in cytoskeletal spectrin densities can be observed for different RBC morphologies [28]. During the cell recovery process, the cell can achieve equilibrium and corresponding reference conditions (i.e., ∆A 0 (t), C 0 (t), A 0 (t), A k,0 (t), and V 0 (t)) at different cell features to that of the original state. However, a minimal change in cell features indicates a higher integrity of the cell. Therefore, the recoverability of stretched cells at t storage-period with a defective cytoskeleton of A Vertical = 0, 1, 10 and 25 (%) to its original cell state with A Vertical = 0% was investigated in terms of cell dimensions. H 1 , H 2 and H 3 cell dimensions were estimated and used to estimate the average percentage error (ε) for H 1 , H 2 and H 3 for the analogous original and recovered cell. The evolution of ε of recovered stretched RBCs with varying A Vertical with respect to its analogous recovered cell state (i.e., recovered RBC shape at t and A Vertical = 0%) at F ext = 60.0 and 200.0 pN is presented in Figure 11. It can be observed (refer to Figure 11) that the trend of evolution of ε increases with storage-period for similar morphologies and decreases during the transition period where irregular RBC morphologies occur. The discocyte morphology during 0 ≤ t < 20 and sphero-echinocyte morphology during 26 ≤ t ≤ 42 each has equivalent µ(t), κ(t), ∆A 0 (t) and C 0 (t), but varying A 0 (t), A k,0 (t), and V 0 (t). The cytoskeletal reference state of the CG-RBC membrane model is considered to be an ellipsoidal shape at around 0.94 of the volume of a healthy RBC. Therefore, the cell recovery process after direct axial cell stretching is affected by the increasing cell sphericity and cytoskeletal reference state such that a cell with increased cell thickness results. In addition, the stretching response of cells with similar morphology indicates increasing change in cell diameters with respect to the applied stretching forces. Therefore, the discocyte at day 19 and the sphero-echinocyte at day 42 need to recover a higher change in cell diameters than that of the day 0 discocyte and day 26 sphero-echinocyte, which can be difficult. Similarly, the decreasing cell stretchability with the storage time with respect to increasing µ(t) reduces the cell elongation to be recovered, and therefore indicates better cell recoverability during the transition period. In addition, the cell sphericity increases with echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies, and therefore the cytoskeletal reference state and the evolution of A 0 (t), A k,0 (t), and V 0 (t) supports better cell recovery during the transition period. The sphero-echinocyte during 26 ≤ t ≤ 42 with increasing A Vertical at F ext = 60.0 pN indicates better recoverability than that of a healthy RBC. Similarly, it is difficult to achieve a trend to represent the evolution of ε for sphero-echinocytes during 26 ≤ t ≤ 42 at A Vertical = 0%, 1% and 10% and at F ext = 200.0 pN, such that all the numerically predicted ε values are contained within. These observations for sphero-echinocytes during 26 ≤ t ≤ 42 with increasing A Vertical may be due to the increasing spicule irregularity that can influence the estimation of cell dimensions, and therefore ε (refer to Figure 12). The recovered stretched cell shapes of day 42 sphero-echinocyte RBCs with A Vertical = 0% and 25% at F ext = 60.0 pN, and that of RBCs having A Vertical = 0% and 25% at F ext = 200.0 pN, are presented in Figure 12. In general, the numerical predictions indicate that the recoverability of a defective RBC membrane with A Vertical is adversely affected due to its low cytoskeletal integrity. As expected, the day 0 discocyte with low or zero A Vertical indicates the best recoverability, whereas the day 42 sphero-echinocyte with higher A Vertical indicates the worst recoverability at high stretching forces. The RBCs with prolonged storage have limited recoverability following mechanical deformations. RBCs are required to deform significantly and repeatedly during their passage through microvasculature, and RBCs with lower recoverability may not survive in-vivo circulation leading to adverse post-transfusion outcomes.

Conclusions
A numerical investigation was performed to systematically investigate the deformability and morphology changes of RBCs undergoing storage lesion. It is hypothesized that discocyte, echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies and associated changes in cell membrane mechanics and cell geometry during in-vitro storage, can be attributed to a CG-RBC membrane approach integrated with the BCM-based membrane free energy model. The CG-RBC membrane model [76] captures the essential features of RBC echinocytosis during in-vitro storage, and enables valuable predictions of the influence of storage lesion on RBC quality and viability. The predicted morphology and deformability characteristics of RBC echinocytosis confirms the decreasing cell deformability, and thereby the decreasing quality and viability of stored RBCs with increasing storage duration. Differences in deformability and morphology characteristics: cell shape factor, elongation index and membrane spicule details such as number of spicules, average spicule height, average spicule surface area and average spicule curvature, were observed between the storage lesion-produced and echinocytogenic shape-transforming conditions-produced echinocyte I, echinocyte II, echinocyte III and sphero-echinocyte morphologies. Furthermore, the results predict the adverse influence of membrane structural defects due to storage lesion on cell deformability and recoverability during mechanical deformations, and indicate further degradation of cell viability.
RBC storage lesion is a significant issue facing transfusion treatments, and despite the improved in-vitro storage procedures, storage lesion still occurs, affecting the maximum storage time for a healthy and safe RBC supply. Therefore, improved understanding of the interrelation between quality and viability of stored RBCs and storage lesion is essential to prevent or reduce transfusion related adverse outcomes. The present numerical study provides a strengthened approach for such understanding, and can be extended for a detailed discussion of the relationship between storage lesion and RBC membrane mechanics, where any associated changes to the lipid-bilayer and/or cytoskeletal spectrin network should be considered. In addition, the availability of any physicochemical models that discuss the influence of storage environment on red blood cells along with these numerical investigations on RBC membrane mechanics would facilitate better predictions of the storage lesion. The present study assumes a homogeneous population of healthy young RBCs at day 0 of storage. However, RBCs are a heterogenous population with different ages, and they are also affected by donor characteristics such as donor age, gender and health. Therefore, this study can be extended further to investigate the influence of these donor characteristics on RBCs during storage. Insights gained from such investigations will facilitate better storage procedure and transfusion practices. Moreover, this numerical approach provides a potential framework to study the membrane vesiculation and RBC deformation characteristics under varying pathophysiological conditions, and the influence of altered membrane proteins on RBC function such as hereditary haemolytic disorders (i.e., spherocytosis, elliptocytosis and ovalocytosis) for better diagnostics and treatments. Finally, this study is currently being extended to investigate the recoverability of deformability and morphology changes of RBCs undergoing storage lesion.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. CG-RBC Membrane Model
The CG-RBC membrane model is composed of N V vertices to represent the cytoskeletal actin junctional complexes, which form a 2D triangulated surface of N t triangles. The N S adjacent vertex-vertex connections represent the cytoskeletal spectrin links. The free-energy of the CG-RBC membrane model (E) takes into account the elastic energy of the cytoskeletal spectrin network (E Stretching ), bending resistance of the lipid-bilayer (E Bending ), and constraints of fixed total surface area (E Sur f ace−area ) and enclosed volume (E Volume ). A BCM-based membrane bending energy approach [68,[88][89][90][91] was implemented via bilayer-leaflet-area-difference constraint (E Area−di f f erence ) to achieve the stomatocyte, discocyte and echinocyte morphologies at corresponding reference bilayer-leaflet-area-difference (∆A 0 ). In addition, the total-membrane-curvature constraint (E Total−curvature ) was implemented to achieve numerically consistent RBC morphologies at corresponding reference total-membrane-curvature (C 0 ). Therefore, E is given by, E Stretching was estimated based on the coarse-graining approach by Fedosov et al. [58], and is given by, where k B is the Boltzmann constant, T is the absolute temperature, l max is the maximum extension of j th link, p is the persistence length and k p is the power function coefficient. x j is defined as x j = l j /l max for the j th link having l j length, and m is an exponent such that m > 1. The experimentally estimated membrane shear modulus (µ 0 ) of a healthy RBC lies between 4-12 µNm −1 [58,68], and can be determined as below for the CG-RBC membrane model [58].
where l 0 is the equilibrium length of the j th link, and defined as x 0 = l 0 /l max . E Bending of the lipid-bilayer for a zero spontaneous curvature is estimated by where κ is the membrane bending modulus, M j is the membrane curvature at the j th link, and ∆A j is the membrane surface area associated with the j th link. M j and ∆A j , corresponding to the triangle-pair composed of T1 and T2 triangles that share the j th link, are estimated as follows; where θ j is the angle between outward normal vectors to the triangles T1 and T2, and A T1 and A T2 are the planer area associated with T1 and T2, respectively. θ j is defined such that the concave arrangement of T1 and T2 corresponds to a positive θ j , whereas the convex arrangement corresponds to a negative θ j , and results in positive or negative M j , respectively. The energy components due to constraints: E Sur f ace−area , E Volume , E Area−di f f erence and E Total−curvature are estimated by [58,76,92], where A 0 is the reference membrane surface area, A is the instantaneous membrane surface area, A k,0 is the reference area of k th triangle, and A k is the instantaneous area of k th triangle. V 0 is the reference cell volume and V is the instantaneous cell volume. D 0 is the monolayer thickness of the lipid-bilayer, ∆A is the instantaneous bilayer-leaflet-area-difference, and C is the instantaneous total-membrane-curvature. k A , k a , k V , k AD and k C represent the total surface area, local surface area, volume, bilayer-leaflet-area-difference and total-membrane-curvature constraint coefficients, respectively. It is assumed that the membrane vertices move over the RBC membrane surface to achieve the minimum free energy state, which is the equilibrium RBC shape. The force (F i ) acting on the i th vertex at point r i on the surface is derived from the principle of virtual work, such that; The resulting motion of the i th vertex is then estimated from the Newton's second law of motion as follows; where f ext i is the contribution from any external forces on the i th vertex point, m i is the mass of i th vertex point, dot (.) is the time derivative and c is the viscosity of the RBC membrane. The updated velocity (ṙ i ) and the position (r i ) of the i th vertex at time (t + ∆t) from the time (t) is given as; In these simulations, m i and c in Equations (A12) and (A13) do not affect the equilibrium shape of the RBC as these parameters only control the convergence speed of the simulation from the initial to the equilibrium state, and are purely numerical inputs for the present steady-state numerical predictions. Therefore, c in Equation (A13) acts as a damping coefficient to achieve the numerical consistency during the simulation. The iterations are continued until the RBC membrane reaches the equilibrium state, which is the minimum free energy state of the RBC membrane at given reference conditions. In the present computational implementation, the equilibrium cell state is acknowledged and the derivation is terminated when the change between each analogous energy component at two successive iterations is less than 1 × 10 −7 in the order of energy component in consideration.

Appendix D. Implementation of Optical Tweezers Stretching
The optical tweezers stretching forces were implemented on the CG-RBC membrane model-predicted equilibrium RBC morphologies. Assuming rigid body conditions, the centre of mass of the whole cell and its three-principal axes of inertia were determined for each RBC morphology. The total stretching force (F ext ) was applied on N + = a N V vertices, whereas −F ext is applied on N − = a N V vertices along the principal axis of inertia 3 [39,58] (refer to Figure A2). N + and N − are the vertices that locate within the circular region of radius d C /2 on the initial spherical geometry, and from the 2 vertices (i X max and i X min ) on the furthest ends of the equilibrium cell shape along the principal axis of inertia 3 . d C is the contact diameter between the cell membrane and the attached silica beads. The vertex fraction a corresponds to d C and is given by a = π d 2 C /(4 A 0 ). Therefore, f ext i is applied on the i th vertex, such that; Figure A2. Optical tweezers stretching implementation on the discocyte cell; (a) identification of i X max and i X min vertices, (b) identification of N + vertices, and (c) contact region between the RBC membrane and binding silica beads.
F ext was gradually applied on the cell via regular force increments of ∆F ext , where the cell is provided with sufficient time to converge to the equilibrium stretched state after each force increment. The equilibrium stretched cell state was determined at equivalent A 0 , A k,0 , V 0 , ∆A 0 and C 0 for stomatocyte, discocyte and echinocyte cell shapes and at corresponding F ext . The resultant force and motion of i th membrane vertex was determined according to Equations (A11) and (A12) for F ext external force, and the equilibrium cell state was determined. The axial diameter (D A ) was measured along the principal axis of inertia 3 , whereas the transverse diameter (D T ) was measured along the principal axis of inertia 2 of the equilibrium stretched cell state at F ext . Refer to Figure A3 for a graphical representation of D A and D T .