Nonlinear modeling for predicting red blood cell morphological transformations

Developing theoretical models to describe cell behavior is a crucial aspect in cell mechanics for gaining a deep understanding of the physical principles governing cellular processes. These models can provide a framework for predicting and analyzing cellular responses to various mechanical stimuli and environmental conditions, which is essential for various fields such as bioengineering, disease diagnostics and treatments, and drug development. Taking the deformation of red blood cells (RBCs), for example, a theoretical model can offer invaluable insights into the underlying mechanisms that enable RBCs to navigate through the complex vascular network and perform their pivotal role in oxygen transport. In addition, it can aid in predicting how changes in RBC mechanics may manifest in various diseases, such as anemia, sickle cell disease, and diabetes.

Currently, there are two main models for describing cell deformation, continuum, and discrete models.^1–4 Continuum models treat the cell as a continuous material, applying principles of mechanics to understand how forces induce deformations.^5–10 They often employ constitutive equations to describe the stress–strain relationships within the cell membrane, allowing for the analysis of large-scale deformations and overall cell behavior under various mechanical conditions. Among these equations, the neo-Hookean (NH) and Skalak (SK) laws are particularly prominent. The NK law assumes the membrane to be a three-dimensional isotropic and incompressible material, often used for vesicles,⁹ whereas the SK law assumes the membrane to be isotropic in its plane, initially proposed for RBCs.¹¹ Discrete models, such as particle-based approaches, represent the cell membrane as a collection of interconnected particles or elements.^12–17 Each particle follows specific physical laws, providing a detailed representation of the cell’s mechanical properties at a microscopic level. They often incorporate the microstructure of the membrane, offering deeper insights into the local deformations and structural changes within the cell. A cell membrane is comprised of three main components: the phospholipid bilayer, transmembrane proteins, and the skeleton.^18,19 It is generally accepted that the lipid-bilayer contributes to the out-of-plane bending resistance and surface area incompressibility, while the skeletal spectrin network contributes to the in-plane shear deformation, and the cytoplasm contributes to the volume incompressibility of the cell.¹⁵ The discrete models allow to describe a more realistic and accurate RBC representation by incorporating the properties of these membrane components. While continuum models are favored for their computational efficiency and simplicity for large-scale analysis, discrete models provide a more detailed understanding of the micro-mechanical behavior of cells, making them widely used for studying specific cellular components or localized phenomena. Both continuum and discrete models have been extensively employed in analyzing cell deformation across various applications, such as blood-related diseases^20–22 and cellular responses in microfluidic chips.^23,24

In the present work, we aim to develop a nonlinear deformation model based on a discrete model, and investigate the red blood cell morphological transition. It is expected that the nonlinear model is more concise and more robust, offering a streamlined framework that not only accurately captures the essential characteristics of the cell deformation but also facilitates easier application in both research and practical scenarios.

The HK law assumes that the membrane tensions depend linearly on the surface strain, while the NH law shows a nonlinear dependence but without an explicit parameter to control the area dilatation. The SK law not only provides a nonlinear dependence on shear deformation but also has an explicit parameter $C$ to control the area dilation. Several studies have shown that both the NH and SK laws are close to the HK law in the small-deformation regime, while in the large-deformation regime, the NH law is strain-softening but the SK law is strain-hardening.^7,25 Hence, the SK law may be more accurate to model the RBC deformation because available experiments suggest that the RBC membrane is strain hardening.

Compared with the continuum model in Eq. (4), the discrete model in Eqs. (6)–(11) includes both the bending deformation and volume constraint, except the shear deformation and surface area constraint. The more constraints a deformation model integrates, the stronger its ability to distinguish cell shapes. Based on the theory of differential geometry, a surface is uniquely determined, up to isometry by its first and second fundamental forms. This means that if two surfaces have the same first and second fundamental forms, they are essentially identical. The shear energy $U s$ is intended to constrain the arc (spring) length on the surface, and the area- and volume-constraint energies $U a$ and $U v$ are designed to restrain the surface area and volume, respectively. These geometric properties are associated with the first fundamental form of a surface. The bending energy $U b$ is considered to constrain the surface curvature, which is related to the second fundamental form. Although Eqs. (6)–(11) are associated with the first and second fundamental forms of surfaces, they are still not sufficient to uniquely determine a surface. Hence, more constraints are required to be integrated into the continuum or discrete models.

The nonlinear model in Eq. (16) is employed to reproduce the biconcave equilibrium shape of a RBC from an initial stress-free state (SFS). It is reported that the SFS shape 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,³² where the reduced volume is defined as $ν= V 0/(4/3π R 0 3)$ and $R 0$ is an equivalent radius. We here choose $v=0.94$⁠,¹⁵ and the force strength coefficients are listed in Table II. The velocity-Verlet algorithm is used to numerically solve Eq. (5) with taking the initial conditions of $x i 0$ as the SFS shape and $v i 0= 0$⁠. In order to non-dimensionalize the simulation system, we choose the characteristic length, Boltzmann’s temperature, and mass as $l ′=1μ$m, $ε ′=4.142× 10 − 21$ J, and $m ′= 10 − 9$ kg, respectively. The particle mass and viscosity are set as $m= 10 − 9$ kg and $μ= 10 − 7$ Pa $$s, and these values do not affect the equilibrium shape as they only dictate the convergence speed from the initial to the equilibrium state.¹⁰ 

To calibrate the nonlinear model, we conduct comprehensive parametric studies, listed in Table III. Figure 2 shows the effect of shear modulus $k s$⁠. Within the range $k s∈[6.0× 10 2,2.4× 10 3]$⁠, the RBC equilibrium shapes are nearly identical, exhibiting similar values of $ΔA$ and $Δ A ′$⁠. For $k s≤1.2× 10 2$⁠, the biconcave shape becomes asymmetric, with one side less concave and the other side more concave. As $k s$ exceeds $6.0× 10 3$ and increases further, the RBC shrinks with smaller radius. This is consistent with the fact that the spring becomes less extensible as the shear modulus is increased. Note that the linear shear modulus is measured in experiments about $k s=6.0× 10 2$ (i.e., $2.5× 10 − 6$ N/m) for a healthy RBC.³² 

Figure 3 shows the effect of bending modulus $k b$⁠. Within the range $k b∈[4.8,96]$⁠, the RBC equilibrium shapes are nearly symmetrical and have the similar $ΔA$ and $Δ A ′$⁠. As $k b$ exceeds 384 and continues to rise, the equilibrium shapes become less symmetric and more spherical. It is found that the bending modulus is a crucial factor in determining the local morphology of the equilibrium configuration. When it is large enough, the cell becomes so rigid that it rotates several rounds in the simulation domain. Consequently, the area differences $ΔA$ and $Δ A ′$ exhibit deviations from the reference values. Experimental measurements have determined the bending modulus to be $k b=48$ (⁠ $2.0× 10 − 19$ J).³² 

Figure 5 shows the effect of DDAD coefficient $k a d$⁠. The RBC equilibrium shape is not significantly changed as $k a d$ increases from 18 to $1.8× 10 5$⁠. Previous studies^31,32,40,41 have suggested that $k a d$ should be of the same order as $k b$⁠, so that the bending deformation results from the competition between membrane curvature and DDAD. However, Geekiyanage et al.¹⁵ have indicated that using $k a d$ and $k b$ on the same scale results in a significant difference between the reference DDAD and the instantaneous DDAD at equilibrium. Therefore, they recommended using a large value of $k a d=1.8× 10 4$ (⁠ $7.5× 10 − 17$ J), which is considerably higher than $k b$⁠. In such a condition, the bending energy can be considered negligible in comparison to the DDAD energy, which does not align with the effect of $k b$⁠. Consequently, we recommend a middle value with $k a d≈3.75 k b$ in the present work.

Figure 6 shows the effect of DIAD coefficient $k a d ′$⁠. Within the range $k a d ′∈[0.18,1800]$⁠, the RBC equilibrium shapes are nearly identical, with the similar $ΔA$⁠, $Δ A ′$⁠, and $κ$⁠. For $k a d ′<0.18$⁠, the equilibrium shape is no longer symmetrical, indicating that DIAD plays a non-negligible role in achieving the equilibrium shape. Moreover, Geekiyanage et al.¹⁵ have pointed out that the presence of DIAD can more accurately predict RBC shapes in agreement with experimental observations. We set $k a d ′= k a d$ in the present work.

Figure 7 illustrates the effects of the length limits $α L$ and $α U$⁠. Despite considerable adjustments to these limits, their effect is negligible, meaning that our nonlinear model is relatively insensitive to such variations for achieving a biconcave equilibrium shape. Likewise, the DDAD and DIAD limits $β L$ and $β U$ exert a negligible influence on the RBC equilibrium configuration, given that $β L∈[0,0.45]$ and $β U∈[5,40]$⁠. Moreover, the power $n ^$ does not significantly affect the biconcave equilibrium shape, if $n ^∈[1,4]$⁠. Hence, we set $α L=0.2$⁠, $α U=2.5$⁠, $β L=0.05$⁠, $β U=15$⁠, and $n ^=2$ as default values in the present work.

Through the parametric analysis, it is demonstrated that the shear modulus $k s$ affects the cell’s extensibility, and the bending modulus $k b$ is pivotal in dictating local membrane bending. The local area coefficient $k a$⁠, despite its minimal impact on the biconcave shape, is crucial for numerical stability, especially in cells undergoing significant deformation. The DDAD coefficient $k a d$ and DIAD coefficient $k a d ′$ are crucial for controlling the global concavity or convexity of the equilibrium shape. The parameters $α L$⁠, $α U$⁠, $β L$⁠, and $β U$⁠, which set the lower and upper bounds for geometric properties, exhibit a wide acceptable range, indicating the robustness of the nonlinear model.

It has been found that the normal biconcave shape of RBC can be modified into stomatocytes and echinocytes through a variety of agents, such as anionic amphipaths, high salt, high pH, and ATP depletion,^15,36,42 known as stomatocyte–discocyte–echinocyte (SDE) transformation. Stomatocytes are cup-like concave shapes while echinocytes are crenated shapes with a spiculated cell surface. They are further categorized into subclasses I–IV, according to different levels of concavity and crenation.

We apply the nonlinear model to reproduce the SDE transformation, as shown in Fig. 8. The model’s predictions closely match the experimental observations and the numerical results of Geekiyanage et al.¹⁵ The slight differences with the experiments are caused by the parameter setting of $ΔA$ and $Δ A ′$⁠. Figure 9 shows quantitative comparisons of geometrical properties. During the transformation from the SFS ellipsoid to the discocyte, all the geometrical properties, except for $V$⁠, undergo minor variations that are indicative of linear deformation. Consequently, the parameters determining the nonlinear behavior, such as $δ φ L$⁠, $δ φ U$ and $n ^$ in Eq. (17), exert only minor effects on the biconcave equilibrium shape. In contrast, the parameters directly associated with linear behavior, $k ^ φ$⁠, have obvious effects on this shape. Likewise, for the equilibrium shapes of S I–III and E I, the geometrical properties of $l ¯$⁠, $a ¯$⁠, $m ¯$⁠, $A$⁠, $ΔA$⁠, and $Δ A ′$ also experience minor variations. As a result, the linear parameters $k ^ φ$ have a pronounced effect, whereas the nonlinear parameters have a less noticeable impact. For the equilibrium shapes of S IV and E II–IV, the geometrical properties of $m ¯$⁠, $ΔA$⁠, and $Δ A ′$ are largely varied, meaning that their deformation is in the range of nonlinear region. In such cases, the nonlinear parameters also have strong effects on their equilibrium shapes, especially for S IV and E IV. This is exactly the reason why most of previous model^32,36 cannot predict the equilibrium shapes of S IV and E IV.

During the transformation from the discocyte to the stomatocyte, the value of $ΔA$ gradually decreases, indicating an increasing concavity of the RBC. In contrast, $Δ A ′$ is almost unchanged for S I–III, predominantly influencing the shape to have one concave side and one convex side. As the stomatocyte stage develops, the cup radius decreases, and several buds are formed inside the RBC. At the stage of S IV, the cell becomes more round, the cup is almost closed, and the interior buds are attached to the membrane. From the discocyte to the echinocyte, the value of $ΔA$ gradually increases, indicating an increasing convexity of the RBC, as well as the $Δ A ′$⁠. As the echinocyte stage develops, the spicules are formed over the membrane. At the stage of E IV, the cell also becomes round and a lot of sharp projections still attached on its surface. Therefore, both S IV and E IV have the nearly spherical shape, with the small deformation index $κ$⁠, as shown in Fig. 9(h). All the equilibrium shapes have the same surface area $A$ and volume $V$⁠, yet they exhibit different curvature $m ¯$⁠, which results in distinct area differences $ΔA$ and $Δ A ′$⁠.

It is known from the force magnitude in Eq. (17) that the nonlinear model exhibits linear behavior at $δφ≈0$ but demonstrates nonlinear behavior when $δφ$ approaches both the lower and upper limits $δ φ L$ and $δ φ U$⁠. This aligns with the biomechanical response of a cell susceptible to small deformation but significantly resists large deformation. For the linear deformation, such as S I, D, and E II, the force strength parameters $k ^ φ$ have an obvious effect on the RBC equilibrium shape. In the case of nonlinear deformation, such as S IV and E IV, the nonlinear parameters, including lower limit $δ φ L$⁠, upper limit $δ φ U$⁠, and power index $n ^$⁠, are also crucial in determining an equilibrium shape.

Considering the force direction in Eq. (16), the shear force of spring is directed along $∂ l j/∂ x i$⁠. The forces resulting from local and global area variations are directed along $∂ a j/∂ x i$⁠. The force due to volume variation is mainly directed along $∂ n j/∂ x j$⁠, where $n j$ is the outer normal to triangle $j$⁠. Additionally, the forces arising from variations in dihedral curvature, DDAD and DIAD are directed along $∂ m j/∂ x i$⁠. For an equilibrium configuration, the dihedral curvature, DDAD and DIAD are in a balanced state, meaning that bending deformation results from the competition among the local curvature, global DDAD, and global DIAD. Hence, assuming $k a d≈ k a d ′≫ k b$ is not reasonable,¹⁵ as it overemphasizes the contribution of global DDAD and DIAD. It is widely recognized that the local curvature, global DDAD, and global DIAD are linked to the lipid bilayer, and consequently, the lipid bilayer is responsible for the bending deformation.^15,32 The shear force is self-balancing among all the springs and does not directly compete with other forces. Since it is linked to the membrane skeleton, a SFS configuration characterized by a reference length $l 0$ is required to maintain its self-balance. Hence, the membrane skeleton is responsible for the shear deformation.

The forces arising from local and global area variations are in equilibrium, leading to the area dilation as a consequence of the competition between these areas. Lim et al.³² suggested that the local area is associated with the membrane skeleton, with $k a≈2 k s$⁠, contrasting with the majority of previous studies that recommended $k a≈ k A≫ k s$⁠, implying a possible connection to the lipid bilayer. As shown in Fig. 4, the effect of $k a$ on the discocyte’s equilibrium shape is negligible, given that the local area is almost unchanged. However, in situations where the local area experiences large fluctuations, such as in cases S IV and E IV, $k a$ significantly influences the equilibrium shape, as shown in Fig. 10. At a small $k a$⁠, the triangle area undergoes large fluctuations, and the force due to the local area variation is smaller than bending force. As a result, the interior buds and exterior spicules on the membrane is rounder at the small $k a$⁠. It is observed that the results with $k a≫ k s$ are more consistent with experimental data than those with $k a≈2 k s$⁠. Therefore, we agree with the prevailing viewpoint that the local area is also linked to the lipid bilayer as the global area, such that the lipid layer contributes to the area dilation.

The force due to volume variation does not compete with other forces. In fact, it is also different with other forces. Unlike the shear force, it is a global force and not requires the balance among all the local elements. Unlike the forces due to the global variations, such as global area, DDAD and DIAD, it does not compete with any other local forces, but an independent force, guiding the cell towards a reference volume $V 0$⁠, either through contraction or expansion.

We have developed a nonlinear model to describe the cell deformation behavior, including the shear, bending, dilation, budding, and spiculation. The model exhibits linear behavior under small deformation but transits to nonlinearity under large deformations, precisely mirroring the biomechanical response of actual cells. Hence, the nonlinear model is accurate to capture stomatocyte–discocyte–echinocyte transformation with reasonable physical parameters. Additionally, the model incorporates seven geometrical parameters: spring length, triangle area, dihedral curvature, surface area, volume, direction-dependent area difference, and direction-independent area difference. These parameters are expressed in a unified mathematical framework, which endows the nonlinear model with conciseness and user-friendliness. Apart from those, we have observed that the model is insensitive to parameter variations when a cell experiences small deformation, exhibiting the low reliance on exact parameter specification.

In the present work, the nonlinear model has been applied to only reproduce the stomatocyte–discocyte–echinocyte transformation of RBCs. In near future, we anticipate expanding its application to additional scenarios, including cell deformation under shear flow, tube flow, and other relevant conditions. This extension will not only demonstrate the model’s versatility but also its potential to predict shape changes in other types of cells by recalibrating the geometrical parameters based on experimental data or theoretical considerations specific to the cell type of interest. We are confident that this model will facilitate a deeper understanding of cell mechanics and deformation responses.

The authors have no conflicts to disclose.

Sisi Tan: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.