Effects of membrane reference state on shape memory of a red blood cell

Abstract By using a three-dimensional continuum model, we simulate the shape memory of a red blood cell after the remove of external forces. The purpose of this study is to illustrate the effect of membrane reference state on cell behavior during the recovery process. The reference state of an elastic element is the geometry with zero stress. Since the cell membrane is composed of cytoskeleton and lipid bilayer, both the reference states of cytoskeleton (RSC) and lipid bilayer (RSL) are considered. Results show that a non-spherical RSC can result in shape memory. The energy barrier due to non-spherical RSC is determined by the ratio of the equator length to the meridian length of the RSC. Thus different RSCs can have similar energy barrier and leading to identical recovery response. A series of simulations of more intermediate RSCs show that the recovery time scale is inversely proportional to the energy barrier. Comparing to spherical RSL, a spheroid RSL contributes to the energy barrier and recovery time. Furthermore, we observe a folding recovery due to the biconcave RSL which is different from the tank treading recovery. These results may motivate novel numerical and experimental studies to determine the exact RSC and RSL.


Introduction
The red blood cells (RBCs) are the primary constituent of blood. A mature RBC is like a droplet of cytosol enclosed within the cell membrane. A lot of work has been made to examine the mechanical responses of RBCs. For example, two basic cell motion in shear flow are observed depending on the shear rate: (a) a solidlike tumbling motion in which the cell rotates around an axis perpendicular to the shear plane; and (b) a droplike tank treading motion in which the shape and orientation of the cell remain almost steady when its membrane circulates around the cytosol (Fischer et al. 1978;Keller and Skalak 1982;Fischer 2007;Fischer and Korzeniewski 2013). As the shear rate increases, the cell motion changes from tumbling to tank treading. Together with the steady tank treading motion in moderate shear flow, Abkarian et al. (2007) revealed the swinging motion, that the inclination angle of the cell oscillates. In these tank treading motions, the cell maintains its biconcave shape (Abkarian et al. 2007;Dupire et al. 2012). These studies indicate the existence of the energy barrier between the dimple and the rim elements of the membrane.
The membrane energy barrier may be determined by the reference state. Generally, the reference state of an elastic element is defined as the geometry with zero stress, which is also known as 'reference configuration' or 'stress-free state'. The cell membrane is composed of cytoskeleton and lipid bilayer. Both components have their reference states and can be different from each other.
In order to bring energy barrier to the membrane, the reference state must be non-spherical. To date, the main idea for the energy barrier comes from the non-spherical reference state of cytoskeleton (RSC). Two limited RSCs are proposed to reproduce the cell dynamics: (a) the biconcave resting shape (Skotheim and Secomb 2007;Dodson and Dimitrakopoulos 2010;Fedosov et al. 2010;Yazdani and Bagchi 2011;Peng and Zhu 2013), and (b) the spheroidal shape with the same surface area (Lim et al. 2002;Tsubota et al. 2014;Peng et al. 2014;Cordasco et al. 2014;Peng et al. 2015). Nevertheless, these studies do not exclude the possibility of a non-spherical reference state of lipid bilayer (RSL). Sui et al. (2010) showed that a non-uniform RSL can also cause the tumblingto-tank-treading transition. Systematic analysis on the effect of RSL conducted by Sinha and Graham (2015) showed that the spheroid RSL leads to better agreement with experiments than a spherical RSL. However, in the above numerical studies the differences due to RSCs or RSLs are either subtle or hard to be observed experimentally.
In a classic experiments, Fischer (2004) observed the shape memory of RBCs in shear flow. After the flow cessation the cell recovered its biconcave resting shape, and the rim was always formed by the same part of the membrane. Similar recovery motion was also observed in experiments by Braunmuller et al. (2012). Since the effect of the reference state appears at a relatively small deformation, further studies on the shape memory effect may help determine the reference states of both components.
Using a 3D continuum model, Cordasco and Bagchi (2017) showed that the shape memory is related to non-spherical RSC. They also found that as the RSC becomes more non-spherical, the recovery time decreases significantly. However, only two limit RSCs were investigated. There are other intermediate shapes and their effects on shape memory are still unclear. The roles of RSL were not considered either. A non-spherical RSL might change the recovery time or even the recovery mode.
Towards this end, we perform 3D simulations of shape memory focusing on the effect of both RSC and RSL. In the present work, two limited RSCs (spheroidal and biconcave) and other intermediate shapes are investigated. The roles of spherical and non-spherical RSLs are also taken into account. The rest of the paper is organized as follows. We first present a brief description of the membrane model and the numerical method. This is followed by a validation through the optical tweezers stretching test. Numerical results, including the effect of both RSC and RSL on shape memory are illustrated together with discussions.

Membrane mechanics
In the present work, the RBC membrane is considered as a continuum with zero thickness. The membrane force results from a combination of shear strain, area conservation, and bending energies. The in-plane strain energy density is modeled using the Skalak law (Skalak et al. 1973): I 1 and I 2 refer to the two invariants as where k 1 and k 2 are the principal stretch ratios. RSC corresponds to the shape that W SK = 0. The parameters G ¼ 2:5 Â 10 À6 N m À1 and C ¼ 0.5 are the shear modulus and the energy penalty for area dilation, respectively (Dodson and Dimitrakopoulos 2010). The corresponding force f SK is calculated following the method of Charrier et al. (1989). Apart from the Skalak law which includes local area constraints, another global area conservation force is also considered: where k S ¼ 1000G is the area modulus, S is the surface area of the membrane, S 0 ¼ 134lm 2 is the reference surface area, j is the mean curvature and n is the outward normal vector to the surface. The bending resistance follows Helfrich's formulation (Helfrich 1973) for bending energy as where E B ¼ 6 Â 10 À19 J is the bending modulus (Evans et al. 2008) associated with the mean curvature j, and E G is the bending modulus associated with the Gaussian curvature j G . The value of c 0 corresponds to the curvature of RSL for which the bending energy is zero. A non-spherical RSL has spatially varied c 0 .
The bending force f B is calculated as (Zhong-can and Helfrich 1989; Sinha and Graham 2015) where D LB is the Laplace-Beltrami operator (Sinha and Graham 2015). The calculation of j and j g follows the method proposed by Petitjean (2002) and Boedec et al. (2011) and is not shown for simplicity. The internal volume of the cell is conserved due to the incompressibility of the cytosol. The volume constraint force is described as where k V =250 N m À2 is the penalty coefficient, V denotes the cell volume, and V 0 ¼ 94 mm 3 the reference volume (Peng et al. 2010).
The total elastic force of the RBC membrane can be described as The values of k S and k V ensure the variation of area and volume within the range of 1%.

Reference states
The RSC is characterized by the volume ratio V I =V S , where V I is the volume of the RSC and V S is the volume of the sphere with surface area S 0 . The biconcave RSC is also investigated using the expression of Evans and Fung (1972): where g 2 þ f 2 ¼ r 2 . R refers to the cell radius, which is 3.91 lm for a human RBC. The coefficients C 0 , C 2 , and C 4 are 0.207, 2.003, and À1.123, respectively. The volume ratio for this biconcave shape is 0.65. Spherical RSL corresponds to uniform dimensionless spontaneous curvature c Ã 0 ¼ c 0 r 0 , where r 0 is the radius of the sphere with volume V 0 (in the case we consider, r 0 ¼ 2:82 lm). The volume ratio V II =V S is used to characterize the non-spherical RSL, where V II is the volume of the RSL. The spatially non-uniform c 0 are then calculated. Table 1 summarizes the parameter sets studied in this article. Detailed explanation on how to generate the mesh of reference states is given in supplementary materials.

Simulation method
The finite element method is used to solve the mechanical response of the cell. The membrane is triangulated into 1280 surface elements with zero mass. The viscosity of fluid is represented with a damping coefficient l. After discretization of finite elements, the assembled governing equation can be expressed in matrix form as where C is the assembled damping matrix, v g is the assembled velocity matrix, and f int g and f ext g are the global internal and external force matrices. The damping coefficient l is embedded in C, which is related to the shape function of finite element (see the supplementary materials for more information). It should be noticed that the value of l only affects the time scale of the simulation. In the following simulations, the value of l is set as 1 Â 10 5 N m À1 , and the time t refers to a dimensionless time. To model the cell motion under a steady shear flow in the x direction, the external force is applied to the cell as where f x is the force gradient with respect to the y coordinate and e x is the unit vector in the y direction.

Resting shapes
Different parameter sets lead to different resting shapes. The deflation process proposed by Lim et al. (2002) is used to determine these resting shapes.
Starting from an RSC, the cell volume gradually decreases to V 0 ¼ 94 mm 3 , while the surface area remains unchanged as S 0 ¼ 134 mm 2 . During this process, the cell deforms and finally reaches its resting shape. Detailed description of the deflation process is given in the supplementary materials. Figure 1a demonstrates the resting shapes for different parameter sets. The resting shapes for case 2-5 are not shown for clarity, since these shapes are between those for case 1 and case 6. It can be seen that the resting shape for case 1 almost coincides with the parametric cell shape based on Equation 8. As the RSC becomes more non-spherical, the dimple of the resting shape becomes shallower. The cell with nonspherical RSL shows contrary trends. The dimple becomes deeper as the RSL becomes more non-spherical.
The distributions of membrane deformations are shown in Figure 1b and c. The in-plane deformations are quantified through two independent parameters: a 1 ¼ k 1 k 2 and a 2 ¼ k 1 =k 2 . Here k 1 and k 2 satisfy k 1 ! k 2 . According to Section 2.1, a 1 is the area change of membrane elements (a 1 >1 for area expansion and a 1 <1 for compression) and a 2 is the shear deformation. The maximum area and shear deformations occur near the rim, similar to the numerical results obtained by Peng et al. (2015).

Validation
The optical tweezers experiments conducted by Mills et al. (2004) are used to validate the membrane V I =V S , volume ratio of RSC; c Ã 0 , dimensionless spontaneous curvature of spherical RSL; V II =V S , volume ratio of non-spherical RSL. The value 0.65 refers to the biconcave shape based on Evans and Fung (1972).
model. In this set up, two silica beads are attached to the cell at opposite points. One bead is anchored to the surface of a glass slide, and the other is trapped by a laser beam. The cell is stretched by moving the slide while keep the trapped bead at rest. The axial diameter D A (in the direction of stretching) and the transverse diameter D T are measured.
The computational setup follows the methodology proposed by Siguenza et al. (2017). The resting shapes obtained in Section 2.4 are used as the initial shapes. Two contact areas of diameter 2 lm in the opposite direction are chosen, and the stretching force is applied only to the nodes on the edges of the contact areas. As the cell is stretched to its extremity, the axial diameter D A is evaluated by calculating the mean position of each loaded edge, rather than the longest distance between the stretched cell.
The simulation results for different parameter sets are shown in Figure 2. Excellent agreement with experiments is obtained for both D A and D T . It is seen that spherical RSLs barely affect the simulation results, which has been observed by Sinha and Graham (2015) numerically. As the RSC becomes more spherical, both the values of D A and D T decrease, in consistent with the results of Siguenza et al. (2017). They compared simulation results for biconcave RSC and quasi-spherical RSC, and the biconcave RSC led to larger values of D A and D T . More information about the optical tweezer simulation is given in the supplementary materials.

Existence of shape memory
We follow the go-and-stop method promoted by Fischer (2004) to simulate the shape memory, as shown in Figure 3a. The initial shape of a cell is the same as its resting shape. The major axis refers to the longest distance between two points of the membrane on the xy plane. The incline angle u is defined as the angle between the major axis and the positive x direction. The phase angle b is the angle between the   Evans and Fung (1972). All resting shapes are biconcave including cases 2-5. Thus only parts of the shape profiles are shown for clarity. (b) The area deformation a 1 and (c) the shear deformation a 2 of the membrane at different resting shapes. major axis and the vector of the material point and the barycenter. The initial position of the material point is b ¼ 0. At the 'go' process, a force gradient f x ¼ 2 Â 10 À7 N m À1 is exerted on the cell, in order to ensure tank treading. The force gradient is removed when the material point moves to the dimple of the cell, corresponding to the 'stop' process. The starting time of the 'stop' process is labeled as t ¼ 0. u 0 and b 0 represent the corresponding angles at t ¼ 0. Figure 3b demonstrates the recovery process for case 1 when the force gradient is removed at b 0 ¼ p=2. The cell first recover to a biconcave shape and then rotates to its resting shape. The material point can return to the rim of the cell, in agreement with experimental and numerical results (Fischer 2004;Cordasco and Bagchi 2017). Thus our method can reproduce shape memory. It should be noted that the recovery motion is an intrinsic property of the cell. Thus the magnitude of the force gradient has no influence on it (see the supplementary materials for more details about the effect of the force gradient).

Effect of RSC
In this subsection we consider the effect of RSC (case 1-6). The shear forces are removed at b 0 ¼ p=2. The angular displacements jDuj ¼ juÀu 0 j and jbj are shown in Figure 4a and b. The intermediate phase and rotational phase observed by Cordasco and Bagchi (2017) are also denoted. As the RSC becomes more spherical, the recovery time t 1 becomes longer, but the rotational angle of the cell is not affected. Particularly, case 5 and case 6 have almost identical rotational phase, while the decay of jbj for case 5 is slightly faster than that for case 6 during the intermediate phase. The cell length in vorticity direction L z is also plotted (see Figure 4c). The compression phase after the remove of shear forces is shown. After that, the value of L z first shows a plateau, and then a fast increase with slightly overshoot.
The time histories of the total membrane energy W total , the total in-plane strain energy W T , and the bending energy W B are shown in Figure 4d-f. W T is obtained by integrating equation 1 over the entire membrane, and W total ¼ W T þ W B . In the compression phase, both W total and W T decrease while W B increases within a short time. This is followed by a slower decay in the intermediate phase and then a fast decay in the rotational phase for all three energy values. The decrements for W T and W B during the intermediate phase and the rotational phase are almost at the same level. Again, case 5 and case 6 have similar decrements of energy.
It is intuitive to relate the anisotropy of RSC to its volume ratio V I =V S . However, case 5 and case 6 show similar recovery process with different volume ratio. As the RSC changes with fixed surface area S 0 , its equator C e and meridian C m also changes (see Figure  5a, inset). By defining the ratio of the equator to the meridian as d ¼ C e =C m , the variation of d with different RSCs is shown in Figure 5a. As the RSC approaches the sphere, the value of d decreases towards 1, and the degree of anisotropy also decreases. When the cell tank treads in shear flow, its equator and meridian also rotate about the z axis. During the recovery process, the equator and the meridian tend to move to their original states. For case 5 and 6, although their RSCs and corresponding volume ratios are different, the values of d are almost the same. Thus the recovery processes for case 5 and 6 are similar.
According to Cordasco and Bagchi (2017) the rotational phase corresponds to a slower time scale t r . They suggested a reciprocal relationship between t r and the in-plane energy barrier DE 0 of RSC. The calculation of DE 0 involves two 'ground states': the resting shape with the material point at the rim, and the same shape but with the material point at the dimple. However, the latter state can not be obtained due to shape memory. In order to estimate DE 0 , we use the shape after the compression phase instead. The time scale of the compression phase is independent of RSC, and we choose t ¼ 0.5 as the end of the compression phase for all cases. The estimated energy barriers are DE 0 % 2$9 Â 10 À18 J, within the same range suggested by Cordasco and Bagchi (2017). As plotted in Figure 5b, the variation of t r is almost proportional to E B =DE 0 , in great agreement with the hypothesis of Cordasco and Bagchi (2017). For case 5 and 6, the values of d are almost the same, thus these two cases have similar in-plane energy barrier DE 0 . The deviations for case 4-6 may come from the difference  between the shape after the compression phase and the resting shape with marker point at the dimple.

Effect of RSL
Now we turn to the effect of RSL. The comparisons between case 6 and case 7 are shown in Figure 6. On the one hand, both cases show similar trends, indicating that the recovery motion is dominated by the biconcave RSC. On the other hand, the spheroid RSL leads to slower decay of energy in the intermediate phase, and more cell rotation. It is interesting to see that for case 7, the decrement of W T is smaller while the decrement of W B is larger. Figure 7a shows the recovery response of case 8 after the remove of the shear forces at b 0 ¼ p=2. The cell shows the folding mode assumed by Fischer (2004) rather than tank treading mode as case 1. After the remove of the shear forces, the dimples first appear in the middle of the cell. Then they become shallower and move outward. At the same time, another pair of dimples appear in the perpendicular direction to the former pair of dimples. Finally the newly generated dimples stay while the former  dimples become the rim. Due to the folding motion, the recovery time is far less than all the above cases.
If the shear forces are removed at b 0 ¼ 6p=4, the cell would deform in tank treading mode, as shown in Figure 7b. Under these circumstances, the recovery time is longer than the case with b 0 ¼ p=2. However, this tank treading recovery is different from those for case 1-6. As shown in Figure 3b, the cell maintains biconcave during the recovery process. While in Figure 7b, the dimples first appear in the vicinity of the material point. Then they become shallower and disappear, leading to a flat disk shape. Finally another pair of dimples appear away from the material point. Thus the recovery response exhibits some degree of folding mode.
The comparison of L z is shown in Figure 8a. It can be seen that there is no plateau as Figure 4, and the values of L z just increase and then decrease. When b 0 ¼ p=2, the value of L z even show some oscillations before reaching stable. Another significant difference between case 1 and 8 is the trace of the material point as shown in Figure 8b. The circles represent the starting points, and the dots the end points. The material points first slant to the left and up direction, indicating the compression phase. Then the material point for case 1 moves in a reverse curve to the right at b 0 ¼ p=2 or p=4. However, the material point for case 8 slants down in a straight line at b 0 ¼ p=2, and arcs towards the left at b 0 ¼ p=4.
The variations of membrane energies are shown in Figure 9, for which the shear forces are removed at b 0 ¼ p=2. Through the comparisons we can see that for case 8 the recovery response is dominated by the biconcave RSL. The decay of membrane energy is faster for biconcave RSL, and the decrement of W B is about twice as much as that of W T . Just like the comparison of case 6 and 7, the biconcave RSL leads to lower W total , W T , and W B when comparing to spherical RSL.

Discussion
It has been shown by Cordasco and Bagchi (2017) that the recovery time is longer as the RSC becomes more spherical. They also suggested that the slower time scale t r is inversely proportional to the in-plane energy barrier DE 0 . However, they only simulated two RSCs without intermediate states between these two limits. The present results fill the blank and we propose an estimation method for energy barrier. We  show that the energy barrier due to non-spherical RSC is related to the value of dthe ratio of the equator length to the meridian length of the RSC. Thus different RSCs may have similar energy barrier, which leads to identical recovery response. The estimated t r and DE 0 also coincide with the inversely proportional relationship.
In numerical simulations, a spatially non-uniform RSL can also contribute to the energy barrier of membrane, thus lead to tumbling-to-tank-treading transition (Sui et al. 2010) and shape memory (Gounley and Peng 2014;Niu et al. 2015;Gou et al. 2018). The contribution of RSL seems qualitatively equivalent to RSC. We demonstrate that a non-spherical RSL results in the increase of rotational angle and recovery time comparing to spherical RSL. Thus there might be a time scale t b related to the bending energy barrier DE b0 . As the RSL becomes the spherical shape, DE b0 ! 0. Since a spherical RSL has no shape memory, t b ! 1. Therefore t b could also be inversely proportional to DE b0 .
The present study uncovers a folding mode of shape recovery related to biconcave RSL. Unlike the previously observed tank treading mode, an important feature of folding mode is the appearance and disappearance of dimples, and the recovery time is significantly shortened. Cordasco and Bagchi (2017) also reported a folding-like recovery when the bending modulus decreases significantly and the RSC becomes biconcave. With these results, we hypothesize that the folding mode is determined by the difference between two energy barriers DE 0 and DE b0 . Considering a dimensionless parameter defined as DE Ã ¼ jDE 0 ÀDE b0 j E B . On the one hand, in case 8 the biconcave RSL results in a bending energy barrier as DE b0 %10 À17 J, greater than the in-plane energy barrier due to the spheroidal RSC as DE 0 %10 À18 J. The calculated DE Ã is in the order of 10 2 , higher than some threshold value. Thus case 8 recovers in folding mode. On the other hand, in the study of Cordasco and Bagchi (2017) the biconcave RSC results in the in-plane energy barrier as DE 0 %10 À18 J, while the spherical RSL with zero bending energy barrier. The corresponding DE Ã value is in the order of 10, less than the threshold value, thus the recovery response shows tank treading mode. With about 10 times lower bending modulus, the DE Ã value becomes 10 times larger and beyond the threshold. Then the folding-like recovery becomes prominent. Additional simulations with other non-spherical RSLs and elastic modulus may help examine the transition from tank treading mode to folding mode.
The present study mainly focuses on the recovery response started from an elongated discocyte. In shear flow, the cell is observed to show other shapes like stomatocyte, trilobe, hexalobe, etc. (Lanotte et al. 2016). After leaving a constriction the cell adopts a parachute shape (Braunmuller et al. 2012). The recovery response from these shapes may demonstrate other properties of the membrane that help determine both the RSC and the RSL. For these purposes, further numerical and experimental studies on the recovery response of RBCs are needed in future works.