A coarse-grained red blood cell membrane model to study stomatocyte-discocyte-echinocyte morphologies

An improved red blood cell (RBC) membrane model is developed based on the bilayer coupling model (BCM) to accurately predict the complete sequence of stomatocyte-discocyte-echinocyte (SDE) transformation of a RBC. The coarse-grained (CG)–RBC membrane model is proposed to predict the minimum energy configuration of the RBC from the competition between lipid-bilayer bending resistance and cytoskeletal shear resistance under given reference constraints. In addition to the conventional membrane surface area, cell volume and bilayer-leaflet-area-difference constraints, a new constraint: total-membrane-curvature is proposed in the model to better predict RBC shapes in agreement with experimental observations. A quantitative evaluation of several cellular measurements including length, thickness and shape factor, is performed for the first time, between CG-RBC model predicted and three-dimensional (3D) confocal microscopy imaging generated RBC shapes at equivalent reference constraints. The validated CG-RBC membrane model is then employed to investigate the effect of reduced cell volume and elastic length scale on SDE transformation, to evaluate the RBC deformability during SDE transformation, and to identify the most probable RBC cytoskeletal reference state. The CG-RBC membrane model can predict the SDE shape behaviour under diverse shape-transforming scenarios, in-vitro RBC storage, microvascular circulation and flow through microfluidic devices.


Introduction
Red blood cell (RBC) is a unique cell without any nucleus or mitochondria [1][2][3] and is remarkably simple in its structure.The most important function of RBC is the transfer of oxygen to body tissues and RBC is adapted with many features to maximize its performance as a gas carrier.It is extremely deformable and elastic to sustain its passage through narrow capillaries of the microvasculature [3,4].RBC holds ~40% excess surface area as compared to a sphere with the same volume [2,5], and the higher surface area allows an increased gas transfer across its surface [6].The structural integrity and the stability of the RBC is due to its membrane which is comprised of three main components; the phospholipid bilayer, transmembrane proteins and the cytoskeleton [4,7].The cytoskeleton is mainly composed of spectrin tetramers connected at actin junctional complexes in a triangular network form [7,8], and is the major contributor to the highly elastic nature of the RBC, whereas the lipid bilayer contributes to the viscosity and area preserving nature of the RBC [9,10].
It has been found that the normal biconcave shape of RBC can be modified into stomatocytes and echinocytes through variety of agents, some of them such as amphipaths, extracellular ionic strength and pH, causing reversible transformations at constant membrane surface area and cell volume [11][12][13][14][15][16].Stomatocytes are cup-like concave shapes whereas echinocytes are crenated shapes with a spiculated cell surface [9,[12][13][14][15]17].Stomatocytes and echinocytes shapes are further categorized into subclasses I, II, III and IV considering different stages of concavity and crenation.The Bessis' nomenclature is commonly used to identify different RBC shapes [18], and the different stages of stomatocyte-echinocyte can be described as follows.Stomatocyte I is a cup shape with a shallow circular invagination; Stomatocyte II is a cup shape with a deeper invagination, still at least approximately circular; Stomatocyte III is a cup shape with a deep invagination, often elongated into a mouth-like slit and sometimes accompanied by other pit-like invaginations; Sphero-stomatocyte (Stomatocyte IV) is a spherical shape with small interior buds still attached to the membrane; Echinocyte I is a disc with several undulations around its rim; Echinocyte II is a flattened elliptical body with rounded spicules distributed more or less uniformly over its surface; Echinocyte III is an ovoid or spherical body with sharper and more numerous  spicules distributed evenly over its surface; and Sphero-echinocyte (Echinocyte IV) is a sphere with small sharp projections still attached to its surface [19].Refer to Fig 1 for a representation of RBC morphology at different stages of stomatocytosis and echinocytosis.
Among the identified many physical and chemical conditions, cationic amphipaths, low salt, low pH and cholesterol depletion can transform normal biconcave RBC into stomatocytes whereas anionic amphipaths, high salt, high pH, adenosine triphosphate (ATP) depletion, cholesterol enrichment, and proximity to a glass surface can transform normal biconcave RBC into echinocytes [12,20].These shape changes are reversible to the stage where sphero-stomatocyte and sphero-echinocyte become smooth spheres when small interior buds and vesicles are pinched off from the RBC surface [17].These reversible RBC shape changes from biconcave shape to-and-from stomatocyte and echinocyte shapes are widely known as stomatocyte-discocyte-echinocyte (SDE) shape transformations.SDE shape transformations have a universal nature for the resultant shape and do not depend on the engaged shape-transforming agent [12].
Sheetz and Singer [21] proposed the widely known bilayer couple hypothesis to explain SDE transformation, and according to the bilayer-couple hypothesis, the two leaflets of the membrane bilayer responds differently to various perturbations and remain coupled to each other.Therefore, the expansion of the outer bilayer-leaflet relative to the inner bilayer-leaflet due to any shape-transforming agent results in convex forms on the membrane surface and induces echinocytes.Similarly, the expansion of the inner bilayer-leaflet relative to the outer bilayer-leaflet due to any shape-transforming agent results in concave forms on the membrane surface and induces stomatocytes.In this manner, the universal nature of SDE transformation is explained through the bilayer-couple hypothesis, and any change in relaxed area-difference between bilayer-leaflets is the necessary and sufficient factor for shape changes.

Computational simulations of RBC
Computational modelling and experimental studies have been carried out to investigate the biomechanical and rheological behaviour of RBC [9,[35][36][37][38][39][40].Due to the simplicity of the RBC structure, they can be numerically approximated as a bag of concentrated haemoglobin solution surrounded by a thin macroscopically homogeneous membrane [13].The thickness of the lipid bilayer being ~4 nm [28] and the offset between lipid bilayer-cytoskeleton being ~30-40 nm [28], RBC membrane can be treated as a two-dimensional (2D) viscoelastic surface in a threedimensional (3D) space in the cellular length scale [9,13].The heterogeneous nature of the RBC membrane can be reasonably approximated as homogeneous in its properties for length scales above 100 nm [13,41], and is suitable for mesoscopic and macroscopic scale investigations.However, a more realistic and accurate RBC representation requires detailed description of its structure, and therefore the incorporation of the properties of lipid bilayer, cytoskeleton and transmembrane proteins.There are several types of computational modelling and simulation approaches for RBC studies, and the prime approaches are continuum, particle-based and hybrid continuum-particle [41] based techniques.RBC membrane and associated fluid components are considered as homogeneous material under continuum approach, whereas the particle-based modelling represents these components via a network of springs and particulate assembly [42], respectively.It is difficult for the continuum approaches to incorporate the detailed structure of RBC membrane, whereas particle-based methods are capable of detailed representation of the RBC membrane but require high computational cost to represent it at the molecular level through molecular dynamics approaches.Therefore, particle-based modelling approach in macroscale and mesoscale are widely used to study RBC behaviour.Besides, hybrid continuum-particle technique is the merging between continuum and particle-based approaches and can lead to more accurate and computationally efficient numerical simulations to study a broad range of biochemical and rheological behaviour of RBC.
For computational simulation purposes, the equilibrium RBC shape was derived by minimizing the in-plane shear energy and the out-of-plane bending energy under the reference constraints of cell membrane area and cell volume [9,42].There are several forms of Helfrich energy; spontaneous curvature model (SCM), bilayer coupling model (BCM), and area-difference-elasticity (ADE), which have been used to describe the RBC membrane bending energy [9,23].The SCM describes the membrane bending energy (E SCM Bending ) for a membrane with surface area, A and local bending modulus, k, such that [23,24]; where, C 1 ðrÞ and C 2 ðrÞ are the principal curvatures at the point r on the membrane surface, and � C 0 is the spontaneous curvature.� C 0 indicates any asymmetry between the two bilayer-leaflets.The SCM model derived, vesicle shape is obtained at the minimum of E SCM Bending for given A and vesicle volume, V.
The BCM is based on the bilayer-couple hypothesis, and assumes fixed area for a membrane lipid molecule and no molecular exchange between the two bilayer-leaflets.Therefore, the area of each bilayer-leaflet remains constant, and the area-difference between the two bilayer-leaflets ðDAÞ can be determined from the integrated mean curvature over the membrane surface, such that [23]; where, D is the distance between bilayer leaflets.The membrane bending energy E BCM Bending in this instance is determined for a defined reference surface with fixed DA as another constraint, such that [23]; The BCM model derived vesicle shape is obtained at minimum E BCM Bending for given A; V and DA.It has been proven elsewhere [24,27,43] that both SCM and BCM models lead to same shape equations, and the vesicle shape behaviour is an extensively studied aspect [13,15,19,23,24,27,28,31,33,34,[43][44][45][46][47][48][49].The ADE model is a combined representation of SCM and BCM, and the ADE model determined membrane energy (E ADE ) for a vesicle having A; V and DA is as follows [23], where, � k is the non-local bending modulus and DA 0 is the reference area-difference between bilayer-leaflets.The ADE model converges to SCM model at � k =k !0 , and into BCM at � k=k ! 1 [23,46,48].The discocyte RBC shape is usually derived from SCM approach at � C 0 e 0 [50], and the following approaches have been followed to predict the SDE transformation.

Computational simulations of SDE transformation
The shape, behaviour and deformability of a RBC have been extensively studied through computational simulation methods [9,35,36] for various conditions and scenarios.However, despite the numerous experimental and theoretical investigations on RBC shape transformations, the numerical predictions on full SDE transformations of RBC are limited [15,47].The existing limited number of numerical shape predictions use ADE membrane bending model to represent the bending resistance of the lipid-bilayer, and no information of the use of SCM and BCM in numerical implementations of SDE shape predictions are available.Nevertheless, the coarse-grained (CG)-RBC membrane model is developed on BCM approach and demonstrates the potential of BCM approach to accurately predict SDE transformation.
The studies by Lim et al. [19,28,49], and Khairy and Howard [31,33] have used a continuum numerical approach based on ADE model to predict different stages of stomatocytes and echinocytes of SDE transformations.These studies have used a single parameter; effectivespontaneous-curvature as the control parameter to predict the different stages of SDE transformation under the constraints of membrane surface area and cell volume.The effective-spontaneous-curvature is a combination of both spontaneous membrane curvature and the reference area-difference of bilayer-leaflets, and at different values of effective-spontaneous-curvature, the corresponding RBC shape has been determined.It has also been identified that the effects of both spontaneous membrane curvature and bilayer-leaflet-area-difference are present in RBC membrane, and their effect is comparable in magnitude [19,31].Furthermore, any lipid asymmetry can be cancelled out by manipulating the bilayer-leaflet-area-difference, while any bilayer-leaflet-area-difference can be cancelled out by manipulating the lipid asymmetry of the lipid bilayer.Therefore, the spontaneous membrane curvature and the area-difference of bilayer-leaflets can be combined, and since, the spontaneous membrane curvature and the bilayer-leaflet-area-difference cannot yet be measured directly, it is also convenient to employ their combined effect through a single control parameter.However, the use of k and � k in the same scale could lead to a significant difference [11] between the reference spontaneous membrane curvature and bilayer-leaflet-area-difference with the instantaneous membrane curvature and bilayer-leaflet-area-difference at the equilibrium.Furthermore, a detailed RBC membrane representation is required to study RBC shape transformations in the presence of membrane inhomogeneity, defective cytoskeleton and blood disorders that affect the membrane structure [8] such as hereditary-spherocytosis (HS) and hereditary-elliptocytosis (HE) [51,52].Therefore, the studies by Lim et al. [19,28,49], and Khairy and Howard [31,33] are difficult to be extended to study the RBC morphology at above diseased conditions.
Chen and Boyle [15] have also developed a CG based spring-particle (SP) numerical approach to predict the SDE transformations, and have suggested a surface curvature approximation that is independent of the network topology to estimate the membrane bending energy.Their ADE model is comprised of the two distinct parameters; spontaneous membrane curvature and reference area-difference between bilayer-leaflets.However, the spontaneous membrane curvature and the reference area-difference between bilayer-leaflets, both being mathematical relationships of the membrane curvature, these control parameters for SDE transformation can be combined together as in the studies by Lim et al. [19,28,49], and Khairy and Howard [31,33].In addition, Chen and Boyle [15] have proposed shape boundaries for a wide range of RBC cell volume and membrane-area-difference.For all these SDE shape predictions; Lim et al. [19,28,49], Khairy and Howard [31,33] and Chen and Boyle [15], the minimum energy configuration of the RBC membrane has been determined by the in-plane elastic energy and out-of-plane bending energy, under the constraints on cell surface area and cell volume.Lim et al. [19,28,49], and Khiary and Howard [31,33] have both employed the non-linear Skalak law to estimate the in-plane shear deformation of the RBC membrane.Chen and Boyle [15], in their CG approach, have adopted the worm-like-chain (WLC) type springs to estimate the resulting elastic energy between these CG particles.The membrane bending energy is estimated through discrete approximations of Helfrich's energy for membrane triangulation.
These computational simulation approaches qualitatively agree well with the experimental observations of RBC SDE transformation.However, all these studies have used the ADE model in which � k =k � 1:0 , and the actual bilayer-leaflet-area-difference at the equilibrium differs from the reference conditions.Therefore, it is difficult to quantitatively compare numerical predictions against the experimental observations as the experimentally extracted bilayer-leaflet-area-difference is at the equilibrium.In addition, it is difficult to find any implementation of SCM method to successfully predict discocyte-echinocyte transformation.Therefore, it is reasonable to expect the presence of any bilayer-leaflet-area-difference with � k=k 6 ¼ 1:0 in order to predict discocyte-echinocyte transformation.Furthermore, if � k=k � 1:0 as in the case of BCM, then the deviation of the equilibrium bilayer-leaflet-area-difference form the reference conditions would be minimal.By this means, a framework to quantitatively validate the developed computational simulations against experimentally observed SDE transformation can be achieved, and will enhance the applicability of numerical predictions to investigate exact SDE transformation conditions.The universal nature of SDE transformation in the presence of echinocytic and stomatocytic shape-transforming agents can be adapted to suit distinctive scenarios if independent measurements of control parameters are available experimentally or computationally.Therefore, a CG-RBC membrane model inspired from the previous research by Lim et al. [19,28,49], is proposed with enhanced applicability to address some of the above concerns, namely, the significant difference between effective-spontaneouscurvature and instantaneous membrane curvature and bilayer-leaflet-area-difference at the equilibrium for ADE model, and the lack of framework to quantitatively validate the developed computational simulations against experimentally observed SDE transformation.The proposed improved CG-RBC membrane model predicts the SDE transformation and is first validated against experimentally observed RBC shapes.It is then applied to study the effect of reduced cell volume and elastic length scale on SDE transformation.Furthermore, the deformability of different stages of SDE transformation and the most probable cytoskeletal reference state of the RBC are investigated using the present CG-RBC membrane model.

Methods
The RBC membrane is the primary contributor for its mechanical nature since a RBC does not possess any internal structure.The lipid-bilayer contributes to its out-of-plane bending resistance and surface area incompressibility while the cytoskeletal spectrin network contributes to the in-plane shear deformation [53], and its cytoplasm contributes to the volume incompressibility of the RBC.The Helmholtz free energy of the RBC membrane is the collective contribution of out-of-plane membrane bending energy, in-plane shear energy and the energy penalty due to cell surface area and volume constraints relative to the specified reference membrane configuration.The equilibrium RBC shape is determined at the minimum free energy state of the RBC membrane.

Free energy of the CG-RBC membrane model
The developed CG-RBC membrane model employs N V particles to represent the actin junctional complexes of the RBC membrane and form a 2D triangulated surface of N t triangles.The N S adjacent particle-particle connections of the triangulated surface represent the cytoskeletal spectrin links.The in-plane shear energy of the RBC membrane ðE Stretching Þ is estimated based on the coarse-graining approach implemented by Fedesov et al. [50], and is composed of an attractive potential ðE WLC Þ in the form of the WLC potential and a repulsive potential ðE POW Þ in the form of a power function.E Stretching can be expressed as follow [50]; where, l j is the length of j th link, k B is the Boltzmann constant, T is the absolute temperature, l max is the maximum link extension, p is the persistence length, k p is the power function coefficient, and m is an exponent such that m > 0: x j is defined as x j ¼ l j = l max .The experimentally estimated RBC membrane shear modulus ðm 0 Þ lies between 4-12 μNm -1 [9,50], and can be expressed as following for the CG-RBC membrane model [50]: where, l 0 is the equilibrium length of spectrin link and defined as x 0 ¼ l 0 = l max .The parameters k p and p can be estimated for a given m 0 and x 0 using Eqs ( 5) and ( 8) at the equilibrium state of specified cytoskeletal reference state.The out-of-plane bending energy of the RBC membrane ðE Bending Þ is estimated based on a discrete approximation of the Helfrich energy model [54] for a zero spontaneous membrane curvature, such that: M j is the membrane curvature at the triangle-pair that shares the j th link, and DA j represents the surface area associated with the j th link.M j and DA j can be given as: where y j is the angle between outward normal vectors to the triangles sharing j th link, and A T1 and A T2 are the planer area of T1 and T2 triangles respectively that share the j th link.The concave arrangement of a triangle-pair results in positive M j whereas the convex arrangement of a triangle-pair results in negative M j .An illustration of the triangle-pair made of T1 and T2 triangles that shares the j th link is provided in Fig 2.
The energy components due to RBC surface area constraint ðE Surface Area Þ and the cell volume constraint ðE Volume Þ are estimated as follows [50,53]: 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, A k is the instantaneous area of k th triangle, V 0 is the reference cell volume and V is the instantaneous cell volume.k A ; k a and k v represent the total surface area, local surface area and volume constraint coefficients, respectively.The resistance of the lipid-bilayer for surface area change is considered for both the whole RBC surface and for individual triangular element, as the lipid-bilayer is anchored to the cytoskeleton through transmembrane proteins and therefore, the movement of lipid molecules over the membrane is restricted [55].The first part of Eq (12) represents the total surface area constraint whereas the second part of Eq (12) represents the triangular element surface area constraint.
A BCM-based membrane bending energy approach is implemented in the CG-RBC membrane model to predict SDE transformations.Therefore, in addition to above energy components, another two constraints are included to the system to maintain a reference areadifference between bilayer-leaflets and to maintain a reference total-membrane-curvature.The area-difference between bilayer-leaflets with D 0 monolayer thickness can be estimated based on a discrete approximation as follows [19,49]: The energy penalty of the membrane due to the constraint to maintain reference bilayerleaflet-area-difference ðE AreaÀ difference Þ is considered to be proportional to the square of the difference between DA and the reference bilayer-leaflet-area-difference ðDA 0 Þ, and can be given as follows: where, k ad is the bilayer-leaflet-area-difference constraint coefficient.Following-up on the discussion in section; 'Computational simulations of RBC', ADE model converges to BCM at � k=k ! 1 , and therefore, it is reasonable for k ad to adopt a significantly higher value than k.
Furthermore, the comparison between Eqs (4) and (15) shows that � k � 4:0 k ad.In the ADE model, k and � k are in the same order [19,49], and therefore the competition between spontaneous membrane curvature and bilayer-leaflet-area-difference, leads to consistent RBC shapes as minimum energy configuration for specified reference conditions.However, for the BCM model, inconsistent RBC shapes can arise at similar DA 0 as � k is of higher order than k, and therefore is the dominant component.DA is the integrated mean curvature over the membrane surface, and therefore DA 0 can be achieved at multitude combinations of directiondependent mean curvature at any point on the membrane surface.Equivalently, in the case of present discrete approximation of DA, the reference DA 0 can be achieved at multitude convex and concave combinations of triangular-pair arrangement over the membrane surface.Therefore, to restrict the inconsistency of resultant RBC shape, the total-membrane-curvature constraint is introduced to the system.The total-membrane-curvature can be presented as the integrated direction-independent mean curvature over the membrane surface.Consequently, the minimum energy configuration of a RBC at specified DA 0 can be obtained at a consistent convex and concave combination of triangular-pair arrangement over the membrane surface, which is also in agreement with the specified reference total-membrane-curvature ðC 0 Þ.The instantaneous total-membrane-curvature ðCÞ for any RBC shape is higher than DA since it considers the absolute value of mean curvature over the membrane surface.C can be estimated as follows: Therefore, the energy penalty of the membrane arising due to the constraint to maintain specified C 0 ðE TotalÀ curvature Þ is defined as follows: where k C is the total-membrane-curvature constraint coefficient, and one can expect k C to be within the same order of � k .The total free energy of the membrane ðEÞ, can therefore be determined by the summation of all these energy components, such that: It is assumed that the membrane particles move over the RBC membrane to achieve the minimum free energy state.The force ðF i Þ acting on the i th membrane vertex point at the point r i on the surface can be derived from the principle of virtual work, such that, The motion of i th point can be described from the Newton's second law of motion as follows: where, m i is the mass of i th particle, dot ð:Þ is the time derivative and c is the viscosity of RBC membrane.

Model parameters for RBC shape prediction
The initial spherical geometry is built up on an icosahedron of equilateral triangular faces that is inscribed within a sphere [19,49,56] of radius R RBC : R RBC is determined for a sphere having equivalent surface area to a RBC.The physiological RBC surface area being ~140 μm 2 [2, 5], A 0 is selected as 140 μm 2 and therefore the estimated R RBC = 3.34 μm.The refinement of triangulation is obtained by generating new vertices at midpoints of each triangular face edge, and combining the vertices together such that the preceding triangular surface is divided into four smaller triangles.The new vertices are then projected radially onto the spherical surface with the radius R RBC .The minimum required particle resolution is determined such that the percentage error ðεÞ for estimating E Stretching ; E Bending ; A; V and DA for the initial spherical geometry falls below 1.0% against the analytical estimations.The present CG-RBC membrane model consists of N V = 2,562 membrane vertex points, N t = 5,120 triangular elements and N S = 7,680 adjacent particle-particle connections.The quality of the triangulation is determined by the distribution of adjacent particle-particle connection length and the distribution of adjacent particle-particle connection for each vertex point [50].The former is attributed to the ratio of standard deviation of the adjacent particle-particle connection length to the average adjacent particle-particle connection length.A smaller value is preferred for this ratio, and for the present consideration the value of this ratio is only 0.065.It is also preferred to have a higher percentage of vertex points having 6 adjacent particle-particle connections.In the present CG-RBC membrane, only the 12 vertices (= 0.47%) on the initial icosahedron have 5 particleparticle connections, whereas all the remaining membrane vertex points (= 99.53%) have 6 particle-particle connections.Therefore, the CG-RBC membrane has high quality of triangulation, and the effects of inhomogeneity can be considered minimal during SDE transformation.The reference state of the cytoskeleton is assumed as an ellipsoid at 0.94 of volume of a sphere having A 0 surface area, and the equilibrium length of the spectrin link l 0 is determined at this reference state.The membrane shear modulus, m 0 is set at 4.0 μNm -1 and agrees with the physiological shear modulus [9,50], whereas the parameters k p ; p and l max are estimated at the room temperature, T = 296.15K, and at m = 2 and x 0 = 0.45 [50].
The experimentally estimated RBC membrane bending modulus lies in the range of 1 x 10 −19 -7 x 10 −19 Nm [9,57], and therefore k is selected to be 2.5 x 10 −19 Nm.The constraint coefficients k A ; k a ; k v ; k ad and k C are set to 1 x 10 −3 Nm -1 , 5 x 10 −5 Nm -1 , 100 Nm -2 , 7.5 x 10 −17 Nm and 2.5 x 10 −17 Nm, respectively to achieve reference constraint conditions.The reference triangular element surface area ðA k;0 Þ is set at the corresponding triangular element area at cytoskeletal reference state.The reference cell volume ðV 0 Þ is considered as 93.48 μm 3 in agreement with the physiological RBC volume [2,5], and is 0.6 of volume of the sphere of radius R RBC .The lipid monolayer thickness ðD 0 Þ is considered as 2 nm.DA 0 =A 0 and C 0 =A 0 are set within the ranges of 0.05% -0.30% and 0.30% -0.90% respectively to predict different stages of SDE transformation.
The motion of RBC membrane particles to reach the equilibrium state, is estimated at c = 1 x 10 −7 Nsm -1 and m i = 1 x 10 −9 kg.The parameters c and m i do not affect the equilibrium RBC shape as they only control the speed of achieving the equilibrium state [54].The updated velocity (_ r i ) and the position (r i ) of the i th vertex at time (t þ Dt) from time (t) is given as follows; where, Dt is the time-step for succeeding iteration.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 of CG-RBC membrane model, the equilibrium cell state is acknowledged when the change between each analogous energy component (E Stretching ; E Bending ; E Surface Area ; E Volume ; E AreaÀ difference ; E TotalÀ curvature ) at two successive iterations, is less than 1 x 10 −7 in the order of energy component in consideration.

Validation of CG-RBC membrane model
Qualitative validation of CG-RBC membrane model predicted SDE transformation.RBC SDE shape transformation is predicted at reference reduced cell volume ðnÞ and fitting combinations of DA 0 and C 0 conditions.n specifies the cell volume relative to the sphere having equivalent cell surface area, and is defined as follows: Several shape classes can be observed as the minimum energy configuration at combinations of different n; DAo and C 0 .The SDE shape transformation presented in Fig 3 is predicted at n = 0.6, 0.05% � DA 0 =A 0 � 0.30% and 0.30% � C 0 =A 0 � 0.90%, and is compared with both experimental observations and different numerical predictions.Experimental observations from SEM imaging (refer to S1 Appendix for detailed information.The laboratory protocol is available at http://dx.doi.org/10.17504/protocols.io.yvhfw36), and the numerical predictions by Lim et al. [19,28,49] and Chen and Boyle [15] are used for the comparison in Fig 3 .It can be observed that the predicted SDE transformation through the implemented CG-RBC membrane model agrees well with the experimentally observed and numerically predicted results.The implemented CG-RBC model is also capable of predicting the sphero-stomatocyte (stomatocyte IV) and sphero-echinocyte (echinocyte IV) stages as stable minimum energy configurations, which have not been predicted by other numerical approaches.The predicted sphero-stomatocyte shape consists with small interior buds that are attached to the membrane surface whereas the sphero-echinocyte shape consists of small sharp projections as expected.Therefore, the CG-RBC membrane model can predict the complete SDE transformation at constant cell surface area and cell volume constraints.
Quantitative validation of CG-RBC membrane model predicted SDE transformation.In addition to the qualitative validation presented in Fig 3, a quantitative analysis is performed on the predicted RBC shapes obtained using the CG-RBC membrane model against experimental data extracted from 3D confocal microscopy imaging of RBC.Different shapes were obtained by changing buffer composition before labelling and fixing the cells (refer to S2 Appendix for detailed information.The laboratory protocol is available at http://dx.doi.org/10.17504/protocols.io.yjyfupw).Imaging data were then analysed using numerical methods to recreate a triangulated mesh over the surface of the cells.Experimental measurements were calibrated using fluorescent beads with a known diameter and imaged at the same time as the RBC.Calibration enabled accurate measurements of RBC surface area and volume.The cell surface area ðA ex Þ, cell volume ðV ex Þ, bilayer-leaflet-area-difference ðDA ex Þ and total-membrane-curvature ðC ex Þ are estimated from these generated triangulated surfaces for four randomly selected RBC having discocyte, echinocyte I, echinocyte II and echinocyte III shapes, and are summarized in Fig 4 .The estimated DA ex =A ex and C ex =A ex for these experimentally obtained RBC are introduced to the CG-RBC membrane model as reference DA 0 =A 0 and C 0 =A 0 parameters, and the resultant shape is determined at equivalent reduced cell volume ðn ex Þ constraint.n ex for the randomly selected RBC each having discocyte, echinocyte I, echinocyte II and echinocyte III shapes, is defined according to Eq (24).

Volume of a sphere having equivalent A ex ð24Þ
Assuming rigid body conditions, the centre of mass of the whole cell and its three-principle axes of inertia are determined for each experimentally observed and numerically predicted cell shapes.H 1 ; H 2 and H 3 are defined as the distance between the furthest vertex points on RBC membrane surface along these three-principle axis of inertias of the cell such that H 1 � H 2 � H 3 , and used to estimate the cellular measurements: the normalized cell length predicted RBC shape.Representative images for the shapes predicted by Lim et al. [19,28,49] through continuum approach (G.H. W. Lim, personal communication, March 19, 2019), and by Chen and Boyle [15] through SP-RBC model ((M.Chen, personal communication, September 13, 2018)) are presented.These images are similar but not identical to the original images from Lim et al. [19,28,49] and Chen and Boyle [15], and therefore for illustrative purposes only.ðH x Þ, normalized cell thickness ðH z Þ and shape factor ðSFÞ. H x is defined as the ratio between H 3 and the equivalent spherical radius ðR � Þ, where R � is the radius of the sphere having equivalent cell surface area.Similarly, H z is defined as the ratio between H 1 and R � : SF ¼ H 1 = ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi H 2 � H 3 p [58], and indicates the sphericity of the cell.The cell becomes more and more spherical as SF reaches the value 1, and becomes more and more flattened disc as SF reaches 0. These three cellular measurements are used to quantitatively compare the corresponding experimentally observed and numerically predicted RBC shapes.The estimated percentage error values ðεÞ for H x ; H z and SF between corresponding experimentally observed and numerically predicted RBC shapes are presented in Fig 5 .It can be observed that the resulting shapes obtained through CG-RBC membrane model agree well with the shapes derived from 3D confocal microscopy imaging.The values of ε for H x ; H z and SF are reasonable and the maximum ε (= 21.65%) is for the H z in the case of echinocyte III RBC shape.However, any experiment error during 3D confocal microscopy imaging and image analysis can lead to ε while triangulated surface generation can also be another contributing factor.Therefore, the resulting ε is the resultant effect of any experimental errors, any triangulated surface generating errors and any error in CG-RBC shape predictions.It can be observed that ε is generally higher for H z than H x , which contributes to ε for SF as well.The experimentally observed cell being not completely lying on the bottom surface of the imaging chamber leads to over prediction of H z .Therefore, better fixation of the cell during 3D confocal imaging experiments can improve the accuracy of H z measurements.In addition, the accuracy of DA ex =A ex and C ex =A ex estimations can be improved by generating a more homogeneous triangular network having similar triangular edge lengths.By this means, it is possible to improve the accuracy of experimentally derived measurements, which can improve the comparability between analogous experimental observations and numerical predictions.However, the maximum ε being 21.65%, is acceptable based on all these uncertainties, and therefore the CG-RBC model is capable of quantitatively represent different SDE morphologies.

Effect of RBC reduced volume on SDE transformation
The CG-RBC membrane model is used to investigate the effect of reduced cell volume ðnÞ on SDE transformation.In an isotonic environment, n is equivalent to 0.6 whereas n increases or decreases in hypotonic and hypertonic environments, respectively.The insight on SDE behaviour at different tonicity conditions facilitate the better understanding of SDE behaviour at specific shape-transforming conditions.Therefore, RBC shapes are predicted at similar reference conditions of DA 0 and C 0 as presented in It can be observed that the predicted RBC shapes become more spherical as expected with the increase in n: H x for the SDE transformation gradually reduces whereas H z increases, and therefore the SF also increases gradually with the increase in n.The membrane curvature corresponding to i th membrane vertex point ðC V;i Þ is defined as follows: where, N S;i is the number of associated spectrin links with the i th membrane vertex point, and M jj and DA jj represent the membrane curvature at the triangle-pair that shares the jj th spectrin link and the surface area associated with it, respectively.It can also be observed (refer to Fig 6) that the maximum C V;i for echinocyte II, echinocyte III and echinocyte IV shapes increase rapidly whereas the minimum C V;i for stomatocyte II, stomatocyte III and stomatocyte IV decrease with increase in n.Therefore, it can be summarized that the absolute C V;i for stomatocyte II, III and IV stages and echinocyte II, III and IV stages increases rapidly with the increase in n whereas no significant change can be observed for stomatocyte I, discocyte and echinocyte I shapes.In general, for the stomatocyte shape, the oval shaped concavity mouth becomes more circular and its minimum cell thickness increases with n.Similarly, the predicted discocyte shape shows asymmetric.biconcavity at low n, which becomes more symmetric with increasing n.The minimum cell thickness of discocyte also increases with n.For the echinocyte shapes, the spicules become narrower and sharper with uneven distribution over the membrane surface with the increase of n.Therefore, the outcomes of the study on the effect of n on SDE transformation suggest that though the predicted SDE shape class has not been affected by the changes in n, the features of the resultant RBC shape is affected.Therefore, the SDE transformations with equivalent reference conditions of DA 0 and C 0 behave differently at different n conditions, and one should consider all the relevant reference conditions to accurately predict the resultant RBC morphology.

Effect of RBC elastic length scale on SDE transformation
It has been found that the RBC membrane bending modulus differs for different SDE morphology [59], whereas the RBC membrane stiffness increases with in-vitro cell ageing and several diseased conditions such as malaria infection, SCA, HS and HE [9].Therefore, the CG-RBC membrane model is used to investigate the effect of elastic length scale ðL el Þ on SDE transformation.L el is the relation of RBC membrane bending modulus ðkÞ to its shear modulus ðm 0 Þ, and is defined as follows [13]: Li et al. [60] have suggested that the L el should be above a critical value (~0.37 μm) at the spherical cytoskeletal reference state in order to obtain the discocyte shape as the possible minimum free energy configuration.Otherwise, stomatocyte shape is resulted as the minimum free energy configuration.However, Mukhopadhyay et al. [13] have used L el = 0.28 μm at zero prestressed cytoskeletal reference state, whereas Lim et al. [19,49] have used L el = 0.28 μm at the ellipsoidal cytoskeletal reference state, in order to investigate RBC shapes and shape transformations.Moreover, Tsubota [54] has investigated the critical value of L el at the shape transition between stomatocyte and discocyte for several forms of out-of-plane bending energy estimations.In the present study, SDE transformation presented in Fig 3 is obtained at L el = 0.25 μm at an ellipsoidal cytoskeletal reference state, and agrees with the critical limits of L el suggested by Tsubota [54] for the corresponding out-of-plane bending energy form.
Several L el ½L el (μm) = 0.11, 0.18, 0.25, 0.35 and 0.56] values are achieved at 5.0, 2.0, 1.0, 0.5 and 0.2 multiplications of m 0 ðm 0 = 4.0 μNm -1 ) at constant k ðk = 2.5 × 10 −19 Nm), where m 0 and k are in the physiological range for cytoskeletal shear modulus and membrane bending modulus.RBC shapes are predicted at similar reference conditions of DA 0 ; C 0 , and n as presented in Fig 3 but at varying L el .The behaviour of SDE transformation at varying L el is illustrated in Fig 7 for stomatocyte II, discocyte and echinocyte II stages, which each represents the stomatocyte, discocyte and echinocyte stages, respectively.The cellular measurements; H x ; H z and SF do not show significant variation with the change in L el .However, the maximum C V;i for echinocyte II, echinocyte III and echinocyte IV shapes decreases rapidly whereas the minimum C V;i for stomatocyte II, stomatocyte III and stomatocyte IV increases with the increase in L el .Therefore, it can be summarized that the absolute C V;i for stomatocyte II, III and IV stages and echinocyte II, III and IV stages decreases rapidly with the increase in L el whereas no significant change can be observed for stomatocyte I, discocyte and echinocyte I shapes.In general, for the stomatocyte shape, the oval shaped concavity mouth becomes more circular and its minimum cell thickness decreases with L el .Similarly, the predicted discocyte shape shows asymmetric biconcavity at high L el while the more symmetric biconcavity becomes less symmetric with increasing L el .However, the minimum cell thickness of discocyte does not indicate significant change with L el .For the echinocyte shapes, the spicules become broader and the spicule size become irregular with the increase of L el .At lower L el values, the contribution from membrane bending energy component becomes less significant and therefore result in higher values of absolute C V;i for equilibrium RBC shapes.On the contrary, the membrane bending energy becomes more dominant at higher L el values and restricts C V;i for equilibrium RBC shapes.Therefore, smoother RBC membrane can be observed and it becomes difficult for equilibrium RBC shapes to reach the exact reference area-difference between bilayer-leaflets and total-membrane-surface-curvature conditions at higher L el values.These results reveal how the features of SDE transformation have been affected by the changes in L el , and therefore provide valuable insight to SDE behaviour at different shape-transforming scenarios.

Deformability of SDE transformation
The RBC membrane elasticity is a critical physiological index [61], and a RBC needs high elasticity to carry out its function as a gas carrier in the microcirculation.The main contributors of RBC deformability are cell geometry, intra-cellular and extra-cellular fluid viscosity and elastic properties of the cell membrane [62,63].Therefore, any change in cell deformability

Cytoskeletal reference state of RBC
The stress-free reference state of the cytoskeleton is a controversial subject [4,42].Several studies [4,60] have suggested that a significantly lower cytoskeletal elasticity is needed to achieve a biconcave shape for a spherical cytoskeletal reference state.Otherwise, a non-spherical reference state; biconcave or ellipsoidal cytoskeletal reference states, must be assumed to achieve stable biconcave shape at physiological cytoskeletal elasticity.Lim et al. [19,28,49] numerically proposed that the cytoskeletal reference state of the RBC is most likely to be an ellipsoid with a reduced volume in the range of 0.925-0.976 of a sphere having equivalent surface area, and Khairy and Howard [31] have also suggested that the cytoskeletal reference state is most likely to be an ellipsoidal shape.
The investigation for the cytoskeletal reference state is carried out also with the implemented CG-RBC membrane model.The equilibrium spectrin length, l 0 is extracted at several cytoskeletal reference states, and used in CG-RBC membrane model to predict discocyte, echinocyte I, echinocyte II and echinocyte III shapes that are analogous to 3D confocal microscopy imaging observations (refer to section 'Quantitative validation of CG-RBC membrane model predicted SDE transformation').A similar approach as in Lim et al. [19,49] is performed to generate cytoskeletal reference states and the CG-RBC membrane model is adopted to represent only the cytoskeletal spectrin network such that the stable minimum energy state is determined at set reference cytoskeletal surface area ðA 0; Cyto Þ, cytoskeletal volume ðV 0;Cyto Þ and cytoskeletal reduced volume n Cyto � � : A 0; Cyto is assumed to be equivalent to that of the RBC ðA 0 Þ whereas the reference triangular element surface area of the cytoskeleton, ðA k;0;Cyto Þ is set at the corresponding triangular element area at the initial spherical geometry having the radius R RBC .Several V 0;Cyto conditions are studied such that n Cyto = 1.00, 0.99, 0.94, 0.89, 0.83, 0.78, 0.73, 0.68 and 0.63.The presence of E Bending , in the form of a stronger bending modulus weakens the contribution from shear modulus and leads to an unstressed cytoskeletal state.In addition, the presence of cytoskeletal shear modulus though in weaker form avoids any numerical inconsistency.Therefore, a significantly higher bending modulus ðk = 5.0 × 10 −18 Nm) is used at the physiological cytoskeletal shear modulus ðm 0 = 4.0 μNm -1 ) in order to predict the resultant cytoskeletal equilibrium state.The constraint coefficients; k A ; k a and k V are set to 1 x 10 −3 Nm -1 , 5 x 10 −5 Nm -1 and 100 Nm -2 , respectively.The influence of any bilayer-leaflet-area-difference and total-membrane-curvature is neglected for cytoskeletal reference state generation.The stress-free equilibrium cytoskeletal state is acknowledged at the minimum free energy state, and corresponding l 0 is extracted at set n Cyto .Next, discocyte, echinocyte I, echinocyte II and echinocyte III predictions described in the section; 'Quantitative validation of CG-RBC membrane model predicted SDE transformation', are repeated at estimated l 0 values.The error ðεÞ is once again estimated for the cellular measurements of H x ; H z and SF between corresponding experimentally observed and CG-RBC predicted discocyte, echinocyte I, echinocyte II and echinocyte III shapes at stated cytoskeletal reference states.It can be observed (refer to Fig 9) that ε is very low for all H x ; H z and SF cellular measurements at n Cyto = 0.83,0.89and 0.94 cytoskeletal reference states, whereas ε for H x and H z is significantly higher at n Cyto = 1.00 and 0.99 cytoskeletal reference state conditions.The near spherical cytoskeletal reference state at n Cyto = 1.00 and 0.99 results in higher cell thickness, which leads to the significantly higher ε at these cytoskeletal reference states.
The qualitative comparison (refer to Fig 10) of CG-RBC predicted discocyte, echinocyte I, echinocyte II and echinocyte III shapes at n Cyto = 0.83, 0.89 and 0.94 cytoskeletal reference states, indicates that results at n Cyto = 0.94 cytoskeletal reference state agree well with the corresponding experimental observations.There is no significant difference between predicted discocyte shape at n Cyto = 0.83, 0.89 and 0.94.Although the echinocyte I shape at n Cyto = 0.83 agrees best with the experimental observation, echinocyte II and echinocyte III shapes at n Cyto = 0.83 do not agree well with the corresponding experimental observations.There is no significant difference between predicted echinocyte II shape at n Cyto = 0.89 and 0.94, though the echinocyte III shape at n Cyto = 0.94 has evenly distributed spicules over the membrane surface which agrees better with the experimental observations.Therefore, the most likely cytoskeletal reference state obtained from the CG-RBC membrane model is an ellipsoidal shape at around n Cyto = 0.94.

Summary and outlook
A CG-RBC membrane model was developed in this study to improve the understanding of the stomatocyte-discocyte-echinocyte (SDE) transformation of RBC.In contrast to the existing ADE-based numerical simulation methodologies, a bilayer coupling model (BCM)-based membrane bending energy was employed in combination with a new additional constraint: the total-membrane-surface-curvature, in order to predict each stage of stomatocyte IV-echinocyte IV RBC shapes.
The CG-RBC membrane model produced similar RBC shapes with comparison to the existing numerically predicted SDE transformation that were based on either continuum or particle-based approaches.In addition to the commonly used qualitative validation of SDE transformation against experimental observations, a quantitative analysis is conducted in details to evaluate the accuracy and reliability of the shape predictions by the CG-RBC membrane model.Both qualitative and quantitative analyses of the CG-RBC membrane model predicted shapes against the experimentally observed shapes have proved that the newly developed CG-RBC membrane model can predict equivalent RBC shapes with reasonable accuracy.The membrane model was then employed to investigate the effect of cell reduced volume ðnÞ and elastic length scale ðL el Þ.Furthermore, the CG-RBC membrane model was successfully employed to identify the suitability of shape factor ðSFÞ to estimate the deformability of each stage of SDE transformation.The predicted RBC shapes confirmed that the cytoskeletal reference state of the RBC is most likely an ellipsoidal shape, and is agreeable with the numerical predictions by Lim et al. [19,28,49] and Khairy and Howard [31] as well.
The CG-RBC membrane model is a general framework to predict the SDE transformation of RBC, and employs a coarse-grained approach to reduce the excessive computational cost for detailed description of the membrane structure.However, the model can be easily refined to represent the detailed RBC membrane structure and any molecular level changes in the presence of shape-transforming agents.The CG-RBC membrane model can be extended to predict the membrane budding and their detachment from RBC to form spherocyte; the RBC shape at higher concentrations of shape-transforming agents or at latter stages of in-vitro RBC storage.Similarly, it is possible to extend the current model to predict the RBC behaviour at defective cytoskeletal conditions; i.e.HE conditions.The equilibrium length of the spectrin links can be set individually, whereas the disruption of vertical connections between lipidbilayer and cytoskeletal actin junctions can be incorporated into the model by cancelling E Stretching for the spectrin links attached to these actin junctions.The loss of membrane surface under HS conditions is possible to be introduced to the model by appropriate adjustments to n.However, the cytoskeletal is under compression as it is attached to a lipid-bilayer having a lower surface area than a healthy cell.Therefore, the cytoskeletal reference state would need adjustments to discuss HS cell behaviour though current CG-RBC membrane model.
In addition, extending the present CG-RBC membrane model to discuss distinct RBC shape-transforming conditions, is a potentially interesting future research avenue.Each shapetransforming agent affects the bilayer-leaflets-area-difference in a unique manner and lead to different stages of SDE transformation.However, the process of changing the area-difference between bilayer-leaflets is unique for each shape-transforming agent.The present CG-RBC membrane model can be extended to incorporate equivalent model parameters at distinct shape-transforming conditions; disease conditions such as HS and HE; occurrence of echinocytic and stomatocytic shape-transforming agents; and RBC aging during in-vitro storage.Furthermore, the proposed CG-RBC membrane model can be easily adapted to consider other hypotheses such as Band 3 conformation [17] and cytoskeleton-bilayer interactions [9] induced RBC shape transformations.
The introduced quantitative analysis framework between experimentally observed and CG-RBC membrane model predicted RBC shapes, forms a feasible approach to better incorporate the biochemical effects into the model.Therefore, the presented CG-RBC membrane model facilitate better understanding of RBC behaviour at different shape-transforming conditions.

Fig 2 .
Fig 2. Illustration of M j at j th link.The shaded planer area is the area associated with the j th link.n T1 and n T2 are the normal vectors to the triangles T1 and T2 respectively.https://doi.org/10.1371/journal.pone.0215447.g002

Fig 3 .
Fig 3. Comparison of SDE shape transformation predicted using CG-RBC membrane model with scanning electron microscopy (SEM) experimental observations and different numerical prediction techniques.The corresponding CG-RBC membrane model parameters of reference reduced cell volume ðnÞ, bilayer-leaflet-area-difference ðDA 0 Þ and total-membrane-curvature ðC 0 Þ are indicated next to the

Fig 5 .
Fig 5. Summary of estimated percentage error ðεÞ for normalized cell length ðH x Þ, normalized cell thickness ðH z Þ and shape factor ðSFÞ between corresponding experimentally observed and CG-RBC membrane model predicted RBC shapes.https://doi.org/10.1371/journal.pone.0215447.g005 Fig 3 and at different n condition ðn = 0.55, 0.60, 0.65, 0.70 and 0.75).The stable minimum energy configurations at specified DA 0 ; C 0 and n are illustrated in Fig 6 for stomatocyte II, discocyte and echinocyte II stages.The BCM model can predict stable vesicle shapes if the ratio of � k=k is above the critical value at n [23, 48], and for the present consideration, � k =k = 150, a value well above the critical ratio of � k =k = (0, 20) for n = [0.55,0.75].