article: Healthy red blood cells (RBCs) have remarkable deformability, squeezing through narrow capillaries as small as 3 microns in diameter without any damage. However, in many hematological disorders the spectrin network and lipid bilayer of diseased RBCs may be significantly altered, leading to impaired functionality including loss of deformability. We employ a two-component whole-cell multiscale model to quantify the biomechanical characteristics of the healthy and diseased RBCs, including Plasmodium falciparum-infected RBCs (Pf-RBCs) and defective RBCs in hereditary disorders, such as spherocytosis and elliptocytosis. In particular, we develop a two-step multiscale framework based on coarse-grained molecular dynamics (CGMD) and dissipative particle dynamics (DPD) to predict the static and dynamic responses of RBCs subject to tensile forcing, using experimental information only on the structural defects in the lipid bilayer, cytoskeleton, and their interaction. We first employ CGMD on a small RBC patch to compute the shear modulus, bending stiffness, and network parameters, which are subsequently used as input to a whole-cell DPD model to predict the RBC shape and corresponding stress field. For Pf-RBCs at trophozoite and schizont stages, the presence of cytoadherent knobs elevates the shear response in the lipid bilayer and stiffens the RBC membrane. For RBCs in spherocytosis and elliptocytosis, the bilayer-cytoskeleton interaction is weakened, resulting in substantial increase of the tensile stress in the lipid bilayer. Furthermore, we investigate the transient behavior of stretching deformation and shape relaxation of the normal and defective RBCs. Different from the normal RBCs possessing high elasticity, our simulations reveal that the defective RBCs respond irreversibly, i. e. , they lose their ability to recover the normal biconcave shape in successive loading cycles of stretching and relaxation. Our findings provide fundamental insights into the microstructure and biomechanics of RBCs, and demonstrate that the two-step multiscale framework presented here can be used effectively for in silico studies of hematological disorders based on first principles and patient-specific experimental input at the protein level.
  Blood is a biological fluid that delivers nutrients and oxygen to living cells and removes their waste products. Due to its particulate nature, blood is considered as a complex non-Newtonian fluid exhibiting intriguing dynamic and rheological behavior depending mainly on flow rate, vessel geometry, and volume fraction of suspending particles, especially red blood cells (RBCs). The deformability of a RBC is determined by the geometry, elasticity, and viscosity of its membrane [1,2]. A healthy RBC (H-RBC) has a biconcave shape when not subject to any external stress and is approximately 8. 0 μm in diameter and 2. 0 μm in thickness. The membrane of a RBC consists of a lipid bilayer contributing to the bending resistance, an attached spectrin network (cytoskeleton) responsible for the shear stiffness, and transmembrane proteins such as band-3 and glycophorin C, bridging the connections between lipid and spectrin domains. Experimental and numerical observations of RBC behavior in flow mimicking the microcirculation reveal dramatic deformations and rich dynamics. The extreme deformability allows RBCs to squeeze through narrow capillaries in microcirculation without any damage. However, this feature of RBCs can be critically affected by parasitic conditions such as malaria [2,3], or genetic factors, e. g. , in sickle cell disease (SCD) [4,5], hereditary spherocytosis (HS) [6], and hereditary elliptocytosis (HE) [7]. For many diseases involving RBCs, it is known that the membrane damage and interactions associated with the lipid bilayer and the cytoskeleton including elastic strength, viscous friction, and integral protein linkages strongly influence the biomechanics of RBCs [8]. For example, RBCs infected with Plasmodium falciparum (Pf-RBCs) become progressively less deformable and more spherical during the intraerythrocytic cycle. Remarkable nanoscale protrusions (knobs) have recently been identified causing significant stiffening effects on the cell membrane [9,10]. Knobs mainly composed of two exported parasite proteins, knob-associated histidine-rich protein (KAHRP) and parasite-derived erythrocyte membrane protein (PfEMP1), deposit at the cytoplasmic face of RBC membrane and form vertical links to spectrin network. Consequently, both the lipid bilayer and the spectrin network are stiffened by the knobs as it is implied by the notable enhancement in shear resistance. In addition to the presence of knobs, actin mining introduced by the parasite invasion may result in an enhanced spectrin network at the trophozoite stage or in a deficient spectrin network at the schizont stage [10,11]. SCD is an inherited blood disorder exhibiting heterogeneous cell morphology and abnormal rheology, due to the polymerization of sickle hemoglobin at high enough concentrations, forming long fibers that distort the RBC shape and dramatically alter their biomechanical properties [4,12]. HS RBC is usually caused by the defects in anchoring proteins involved in vertical interactions between lipid bilayer and spectrin network, whereas HE RBC is a result of defects in spectrin filaments related to lateral interactions in the spectrin network [13]. Protein mutations associated with membrane defects subsequently lead to aberrant cell shape and impaired deformability. Thus, quantifying the deformability of RBCs could play a key role in understanding RBC related diseases. Recent advances in computational modeling and simulation enables us to tackle a broad range of dynamics and rheology problems associated with RBCs [14–16]. Several computational models at the whole-cell level, including spectrin-level and multiscale RBC models [17–23], have been developed and employed to quantify the biomechanical properties and dynamic behavior of RBCs in malaria and other hematological diseases [24–29]. Examples include dynamic cell deformability for various stages of the Pf-RBCs [24–26,30–33] and cell morphological sickling [34–36] and vaso-occlusion phenomena in SCD [27,37]. In these existing models, the membrane is usually considered as a single shell with effective properties that represent the combined effects of lipid bilayer and spectrin network. Under normal conditions, the cytoskeleton is attached to the lipid bilayer from the cytoplasmic side. However, under certain conditions, such as RBCs in SCD and other hereditary disorders, the cytoskeleton may become dissociated from the lipid bilayer [12,13]. Also, in malaria disease, the Pf-RBC undergoes irreversible structural modifications with the deposition of knob structures on the membrane surface and the actin remodeling in the skeletal network [38,39]. The biomechanical properties associated with the bilayer-cytoskeleton interactions strongly influence biorheology, cell function and the onset and progression of RBC diseases. However, the one-component whole-cell models cannot facilitate detailed whole-cell exploration of diverse biophysical and biomechanical problems involving RBCs, such as the bilayer loss in HS, the bilayer-cytoskeleton uncoupling in SCD and the aforementioned multiple stiffening effects of knob structures in Pf-RBCs. For these reasons, there is a compelling need to develop a more realistic RBC representation, e. g. , to endow the spectrin-based RBC models with more accurate structure, hence considering separately the lipid bilayer and cytoskeleton but also include the transmembrane proteins. Recent efforts have been directed towards this approach. For example, a two-component composite model of RBC membrane with explicit descriptions of lipid bilayer, cytoskeleton, and transmembrane proteins has been developed and implemented using coarse-grained molecular dynamics (CGMD) [10,41,42]. This CGMD membrane model has been successfully applied to study the membrane-related problems in RBCs such as protein diffusion and vesiculation in defective RBC membrane [43,44], and the stiffening effects of knobs on Pf-RBCs [10]. Although changes on the biomechanics of RBC membrane, including bending rigidity and shear modulus, in certain diseases can be evaluated by modeling a small piece of cell membrane with the two-component composite model, the whole-cell characteristics strongly related to RBC biomechanics and biorheology are not efficiently depicted by modeling only a portion of the RBC membrane. Recently, a two-component whole-cell model has been developed and implemented using dissipative particle dynamics (DPD) [45,46]. The DPD RBC model also accounts separately for the lipid bilayer and cytoskeleton, but it includes implicitly the transmembrane proteins, Thus, the DPD model is computationally more efficient than CGMD model for RBC modeling at the whole-cell level, which has been applied to investigate RBC response and dynamics in flow. However, the lack of the molecular details in this two-component whole-cell model may limit its predictive capacity in identifying the key factors that cause the reorganization of the RBC membrane. An effective way to address this problem is to incorporate only the necessary molecular information from a molecular-detailed composite membrane model into a more coarse-grained whole-cell model. In this work, we develop a two-step multiscale framework by employing the two-component models; see Fig 1. The only experimental input required is information about the structural defects of the lipid bilayer, the cytoskeleton, and their coupling via the transmembrane proteins. This information can be obtained by Scanning Electron Microscopy (SEM) [47] or Transmission Electron Microscopy (TEM) [40]. We then apply this framework to study the biomechanical characteristics of the healthy and diseased RBCs under static and dynamic tensile forcing. Specifically, we probe the alterations in cell deformability from the development of knob structures and remodeling of the spectrin network of Pf-RBCs. We also assess the implication of dynamic deformation and recovery response of defective RBCs from the bilayer-cytoskeleton dissociation in hereditary disorders such as HS and HE. The rest of the paper is organized as follows. In Section 2, we briefly describe the RBC model and simulation method. In Section 3, we present and discuss our numerical results. Finally, in Section 4, we summarize the findings and present the conclusion. We also provide supplementary material in order to better explain the simulation results.
  In the two-component CGMD model of a normal RBC membrane, the lipid bilayer and cytoskeleton as well as the transmembrane proteins are explicitly represented by CGMD particles [41]. Specifically, three types of CG particles are introduced to represent the lipid bilayer of the RBC membrane (Fig 2A and 2C). The red particles represent clusters of lipid molecules with a diameter of 5 nm, which is approximately equal to the thickness of the lipid bilayer; the black-color particles signify band-3 complexes; the light blue particles immersed in the lipid bilayer are glycophorin C. The volume of one black particle is similar to the excluded volume of the membrane domain of a band-3 protein. However, when band-3 proteins interact with the cytoskeleton, the effect of the cytoplasmic domain has to be taken into account and thus the effective radius is ≈ 12. 5 nm. One third of band-3 particles, which are connected to the spectrin network, are depicted as yellow particles (Fig 2C). The spectrin tetramer is modeled by a chain of 39 beads (grey particles) connected by spring bonds [48]. The corresponding potential has the form, V c y s - s = 1 2 k 0 (r - r e q s - s) 2 (1) where k0 and r e q s - s are the spring constant and equilibrium distance between two spectrin particles, respectively. Spectrin particles that are not connected by the spring potential interact with each other via the repulsive term of the Lennard-Jones potential as follows V r e p (r i j) = 4 ϵ σ i j r i j 12 - σ i j r i j 6 + ϵ r i j < r e q s - s 0 r i j ≥ r e q s - s (2) where ϵ is the energy unit and σij is the length unit, and rij is the distance between spectrin particles. To couple the lipid bilayer and spectrin network, actin junctional complexes (blue particles) are connected to the glycophorin C and the middle beads of the spectrin network are bonded to the band-3 complexes (black particles) which are specifically rendered in yellow particles, as shown in Fig 2C. These bonds are modeled as harmonic springs, given by V c y a - s = 1 2 k 0 (r - r e q a - s) 2 (3) where r e q a - s = 10 nm is the equilibrium distance between an actin and a spectrin particle. For a detailed description of the configuration of the cell membrane and the employed potentials, we refer to Ref. [42]. In the two-component whole-cell model, the cell membrane is modeled by two distinct components, i. e. , the lipid bilayer and the spectrin network (Fig 2B). Specifically, through the DPD approach, each component is constructed by a 2D triangulated network on a membrane surface that is characterized by a set of points with Cartesian coordinates xi, i ∈ 1 ⋯ Nv which are vertices of the 2D triangulated network. Different from the one-component whole-cell model [20,21], the lipid bilayer of the two-component whole-cell model has no shear stiffness at healthy state but only bending stiffness and a very large local area stiffness, whereas the cytoskeleton has no bending stiffness but possesses a finite shear stiffness. The whole-cell DPD model takes into account the elastic energy, bending energy, bilayer-cytoskeleton interaction energy, and constraints of fixed surface area and enclosed volume, hence V (x i) = V s + V b + V a + V v + V i n t (5) where Vs is the elastic energy that mimics the elastic spectrin network, given by V s = ∑ j ∈ 1. . . N s k B T l m (3 x j 2 - 2 x j 3) 4 p (1 - x j) + k p (n - 1) l j n - 1, (6) where lj is the length of the spring j, lm is the maximum spring extension, xj = lj/lm, p is the persistence length, kBT is the energy unit, kp is the spring constant, and n is a specified exponent. The shear modulus of the RBC membrane, μ0, is determined by μ 0 = 3 k B T 4 p l m x 0 x 0 2 (1 - x 0) 3 - 1 4 (1 - x 0) 2 + 1 4 + 3 k p (n + 1) 4 l 0 n + 1, (7) where l0 is the equilibrium spring length and x0 = l0/lm. The bending resistance of the RBC membrane is modeled by V b = ∑ j ∈ 1. . . N s k b 1 - c o s (θ j - θ 0), (8) where kb is the bending constant, θj is the instantaneous angle between two adjacent triangles having the common edge j, and θ0 is the spontaneous angle. Constraints on the area and volume conservation of RBC are imposed to mimic the area-preserving lipid bilayer and the incompressible interior fluid. The corresponding energy is given by V a + v = ∑ j ∈ 1. . . N t k d (A j - A 0) 2 2 A 0 + k a (A cell - A cell, 0 tot) 2 2 A cell, 0 tot + k v (V cell - V cell, 0 tot) 2 2 V cell, 0 tot, (9) where Nt is the number of triangles in the membrane network, A0 is the triangle area, and kd, ka and kv are the local area, global area and volume constraint coefficients, respectively. The terms A cell, 0 tot and V cell, 0 tot represent the specified total area and volume, respectively. The bilayer-cytoskeleton interaction potential, Vint, is expressed as a summation of harmonic potentials given by V i n t = ∑ j, j ′ ∈ 1. . . N b s k b s (d j j ′ - d j j ′, 0) 2 2, (10) where kbs and Nbs are the spring constant and the number of bond connections between the lipid bilayer and the cytoskeleton, respectively. djj′ is the distance between the vertex j of the cytoskeleton and the corresponding projection point j′ on the lipid bilayer, with the corresponding unit vector njj′; djj′, 0 is the initial distance between the vertex j and the point j′, which is set to zero in the current simulations. Physical view of the local bilayer-cytoskeleton interactions include the major connections via band-3 complex and ankyrin, as well as the secondary connections via glycophorin C and actin junctions (Fig 2A), here we consider them together as an effective bilayer-cytoskeleton interaction and model it as a normal elastic force, f j j ′ E, and a tangential friction force, f j j ′ F (Fig 2B). The corresponding elastic force on the vertex j of the cytoskeleton is given by f j j ′ E = k b s (d j j ′ - d j j ′, 0) n j j ′ d j j ′ < d c 0 d j j ′ ≥ d c (11) where dc ≈ 0. 2 μm is the cutoff distance. The tangential friction force between the two components is represented by f j j ′ F = - f b s [ v j j ′ - (v j j ′ · n j j ′) n j j ′ ], (12) where fbs is the tangential friction coefficient and vjj′ is the difference between the two velocities. The RBC membrane interacts with fluid particles through DPD conservative forces, which are defined as F i j C = a i j (1 - r i j / r c) s n i j, (13) where rij = ri − rj, rij = |rij|, nij = rij/rij, aij is the maximum repulsion between particles i and j, rc is the cutoff distance, s = 1 is the most widely adopted for the classical DPD method. However, other choices (e. g. , s = 0. 25) for the envelopes have also been used. Detailed description of these interations can be found in Ref. [21]. In addition, in combination with the total energy V (xi) expressed in Eq 5, we are able to derive the total force acting on particle i of the RBC membrane, F i = - ∂ V (x i) ∂ x i + ∑ j ≠ i F i j C. (14) The RBC model is multiscale, as the RBC can be represented on the spectrin level (NvDPD, S = 23,867), where each spring in the network corresponds to a single spectrin tetramer with the equilibrium distance between two neighboring actin connections of ∼ 75 nm [20,59], equal to the average length of one edge of triangular mesh of RBC membrane in the CGMD model. In such case, both CGMD and DPD simulate the RBC membrane at the spectrin level [21,44]. Thus, the size and density of nanoscale knobs of Pf-RBCs applied in both CGMD and DPD models can be directly derived from measured data in experiments. In addition, the RBC membrane properties, such as shear modulus μ0 and bending stiffness kc can be passed directly from CGMD to DPD, l 0 DPD, S = l 0 CGMD, S, μ 0 DPD, S = μ 0 CGMD, S, k c DPD, S = k c CGMD, S. (15) On the other hand, for more efficient computation, the RBC network can also be coarse-grained by using a smaller number of vertices. The equilibrium spring length of the coarse-grained RBC model is then estimated as: l 0 DPD, C = l 0 DPD, S (N v DPD, S - 2 N v DPD, C - 2). (16) Using a similar geometric argument, the spontaneous angle is adjusted as, θ 0 DPD, C = θ 0 DPD, S (l 0 DPD, S / l 0 DPD, C). (17) In addition, the property parameters of the RBC membrane can be estimated as, μ 0 DPD, C = μ 0 DPD, S, k c DPD, C = k c DPD, S. (18) In this study, we adopt a whole-cell DPD model with NvDPD, C = 9128 in order to simulate the stretching behavior of the whole RBC more effectively. We model an H-RBC using the whole-cell DPD model with the following properties: A cell, 0 tot = 134 μm2, V cell, 0 tot = 94 μm3. Based on the CGMD and previous one-component whole-cell simulations of [21], we choose μ0 = 4. 73 μN/m for H-RBC. Regarding the elastic contribution to the interaction energy, we use the value of the bending modulus derived directly from CGMD simulations, i. e. , kc = 31. 9 kBT for H-RBC, which is approximately 1. 3 ×10−19 J. The corresponding bending constant is set to k b = 2 3 k c = 36. 8 kBT. It has been shown that the bending resistance contributes little to the cell deformation in the stretching test in previous studies [23,26,33], so we keep kb constant in all simulation cases. The simulations are performed using a modified version of the atomistic code LAMMPS. The time integration of the motion equations is computed through a modified velocity-Verlet algorithm [60] with λ = 0. 50 and time step Δt = 1. 0 × 10−5 τ ≈ 0. 23 μs. It takes 2. 0 ×106 time steps for a typical simulation performed in this work.
  Here we investigate the elastic properties of the modeled RBC in health but also at different stages of malaria. To probe the RBC mechanical response and the change of its mechanical properties at different malaria stages, we subject the cell to stretching deformation analogously to that in optical tweezer experiments. As shown in Fig 2D, the T-RBC has a lower knob density (ρknob, DPD = 7 knobs/μm2) and enhanced spectrin network, whereas the S-RBC bears a higher knob density (ρknob, DPD = 12 knobs/μm2) and deficient spectrin network. The total stretching force, Fs, is applied in opposite direction to ϵNv (ϵ = 0. 05) vertices of the cell membrane at diametrically opposite directions. First, we examine the response of H-RBCs and Pf-RBCs to large deformation in comparison with previous experimental measurements [2] and computational simulations based on the one-component whole-cell model [24]. We analyze the changes in DA and DT. Our simulation results show that both DA and DT are in agreement with previous experimental and computational results, see Fig 4A. Since the model has separate components for the cytoskeletion and lipid bilayer, we can get the DA and DT values for each component. The consistency of the DA (or DT) values calculated from the lipid bilayer and the cytoskeleton indicate the strong association between the two components. To investigate the influence of the knob density on RBC deformability, we change the knob density from 4 to 9 knobs/μm2 for T-RBC and 10 to 18 knobs/μm2 for S-RBC. For an H-RBC under tensile force Fs = 110 pN, we obtain that the EI value is about 0. 51. For a less deformable Pf-RBC, we find a considerable decrease in EI value, i. e. , EI ≈ 0. 20 for T-RBC at ρknob, DPD ≈ 4 knobs/μm2 and EI ≈ 0. 15 for S-RBC at ρknob, DPD ≈ 10 knobs/μm2. The decrease in EI values from H-RBC to T-RBC then to S-RBC is associated with a reduction of RBC deformability as the progression of the parasite maturation in Pf-RBC. With the increase of knob density, a further decrease in EI values for both T-RBC and S-RBC is obtained, see Fig 4B, which indicates a further reduction in cell deformability of the Pf-RBC. Our simulation results demonstrate that the knobs, being rigid, contribute to cell membrane stiffness. In other words, we model the Pf-RBC with normal spectrin network by removing the influence of actin remodeling. We find that the EI values increase slightly for T-RBC but decrease slightly for S-RBC (Fig 4B). Nevertheless, in comparison with the change of EI values resulting from the stiffening effects of knobs, the difference of EI values between the modified (enhanced or deficient) and normal spectrin network is relatively small. These simulation results further demonstrate that the presence of the knobs in Pf-RBCs is the primary stiffening factor in the loss of cell deformability, which is also consistent with the recent CGMD simulation study by Zhang et al. [10]. More importantly, we efficiently scale up the particle-based RBC model from a portion of cell membrane to a whole-cell structure leading to more realistic simulations of RBC dynamics and blood rheology. Using the virial theorem, we can obtain the average virial stress tensor over a volume Ω, Π = 1 Ω ∑ i ∈ Ω - m i (v i - v ¯) (v i - v ¯) + 1 2 ∑ j ∈ Ω r i j · F i, (20) where mi and vi are the mass and velocity of particle i, respectively, and v ¯ is the average velocity of particles in the volume Ω. The stress contours of H-RBC, T-RBC, and S-RBC at Fs = 110 pN are shown in Fig 4C. In general, large stress response is observed around the longitudinal axis of the stretched cell, by contrast, small stress response emerges at the two sides of the transverse axis [65]. For H-RBC, the lipid bilayer has little shear resistance but large bending stiffness and helps to maintain the membrane surface area. The cytoskeleton is primarily responsible for the shear elastic properties of the RBC. Consequently, the principal stress of the stretched RBC is essentially reflected on the cytoskeleton, while only a small stress response is observed at the ends of lipid bilayer where the tensile force is imposed. As shown in Fig 4C, more than 90% of RBC’s total tensile stress comes from the cytoskeleton component, while the rest (less than 10%) from the lipid bilayer component. For T-RBC and S-RBC, the presence of the knobs immersed in the lipid domain stiffens the cell membrane, which is reflected in the stress contour of the lipid bilayer of the stretched RBC, see Fig 4C. The lipid bilayer with stiff knobs contributes ∼ 33% of the total tensile stress of T-RBC membrane. Compared with the results from T-RBC, the knobs with a higher density in S-RBC withstand an elevated stress of the stretched RBC. The stiffened lipid bilayer bears ∼ 45% of the total tensile stress of the S-RBC membrane. The associations between the lipid bilayer and the cytoskeleton mediated by specific molecular interactions are essential for the mechanical stability of RBCs [66]. The primary interaction is the connection between the transmembrane protein band-3 and the spectrin via the ankyrin binding sites. The secondary interaction is supported by another transmembrane protein glycophorin C and the actin junctional complexes. To involve these bilayer-cytoskeleton interactions into the two-component whole-cell model, we have simply considered them together as normal (elastic) and tangential (friction) interactions. Two parameters, the elastic interaction coefficient, kbs = 46 pN/μm, and the tangential friction coefficient, fbs = 0. 194 pN⋅μm−1s−1, are accordingly introduced for the normal RBC membrane [45]. We examine cell biomechanics and deformation of H-RBCs subject to stretching tests with the default values of kbs = 46 pN/μm and fbs = 0. 194 pN⋅μm−1s−1 of the bilayer-cytoskeleton interactions. First, we probe the stretching deformation of both the lipid bilayer and the cytoskeleton under different stretching force (Fs), see S3 Fig. From this figure, we find that the values of DA and DT obtained from the two-component whole-cell model are in agreement with those from experimental measurements [2], and the one-component whole-cell model [21]. In addition, we find that there is a small difference in the detachment length, which is defined as the distance along the longitudinal axis from the rightmost part of the lipid bilayer to that of the cytoskeleton. However, the deviations are sufficiently small or purely due to statistical fluctuations under normal stretching forces (Fs ≤ 200 pN); hence, there is no significant bilayer-cytoskeleton detachment in these cases. A visible detachment of the lipid bilayer from the cytoskeleton appears in extreme cases (for example, Fs > 300 pN) due to abnormal stretching forces exerted on the cell membrane. Nevertheless, under normal conditions, the large deformability and the strength of the bilayer-cytoskeleton interactions of the H-RBC can be computationally described by the two-component whole-cell model. To investigate the effect of weakened bilayer-cytoskeleton interactions in defective cell membrane on the deformability of RBCs with stretching test, we present an extensive simulation study on the biomechanical behavior of defective RBCs by varying the values of bilayer-cytoskeleton elastic interaction coefficient, kbs, and tangential friction coefficient, fbs, see Figs 5 and 6. First, we consider the cell membrane with the default value of kbs but different fbs (Fig 5A). Specifically, we chose the value Fs = 140 pN, for which the RBCs undergo large biomechanical deformation without bilayer-cytoskeleton detachment under normal physiological condition. We find that when using the default value of fbs = 0. 194 pN⋅μm−1s−1 or other values of similar magnitude, there is a strong coupling between the lipid bilayer and the cytoskeleton. This is indicated by the overlapped black solid lines (with black squares) and red dashed lines (with red spheres) in Fig 5A. However, assuming a pathological RBC state where fbs is decreased by one or two orders of magnitude, an apparent uncoupling between bilayer and cytoskeleton occurs. Compared to the DA value of cytoskeleton, the DA value of the lipid bilayer is relatively larger when fbs is smaller than 0. 01 pN⋅μm−1s−1, and the detachment length between the lipid bilayer and the cytoskeleton increases as fbs reduces. Moreover, the DT values obtained from the lipid bilayer and cytoskeleton remain almost the same regardless of the changes in fbs, which can also be observed in Fig 5A. It seems that the DT value is insensitive to the variation of fbs. Next, we study the effect of the bilayer-cytoskeleton elastic interaction coefficient as shown in Fig 5B. Similarly, there is no obvious distinction in the deformation of RBC for different kbs at the default value of fbs (0. 194 pN⋅μm−1s−1). However, the bilayer-cytoskeleton uncoupling along the stretching axis occurs when we decrease the values of both fbs and kbs by one or two orders of magnitude. Consistent with the previous simulations on RBC traversing across microfluidic channels and tank-treading dynamics [45,46], the strength of the bilayer-cytoskeleton interaction responsible for the association between lipid bilayer and cytoskeleton is again verified by the RBC stretching tests. Next we compare our simulation results with experimental data for normal and defective RBCs. The results on the functional dependence of the stretching response of normal and defective RBCs on the stretching force, Fs, are shown in Fig 6A. In general, DA grows while DT decreases with increasing Fs. For H-RBC, no bilayer-cytoskeleton detachment is observed due to the strong bilayer-cytoskeleton interactions. The corresponding stress contours of deformed RBCs at Fs = 0,80 and 160 pN are presented in Fig 6B. The tensile stress of H-RBC inherently coming from the cytoskeleton is enhanced as Fs rises, which is in reasonable agreement with other computational results [65]. On the contrary, the lipid bilayer is observed to separate from the underlying spectrin network when a defective RBC undergoes large deformation due to the weakened bilayer-cytoskeleton interactions (Fig 6A). The stress contours of the defective RBC at Fs = 160 pN in Fig 6C show significant shear response at two ends of the stretching sites in the lipid domains coinciding with the location for bilayer-cytoskeleton detachment. Upon external tensile forces, a normal RBC undergoes large mechanical deformation, and it restores to its original state when the external tensile force is no longer applied. However, the weakened bilayer-cytoskeleton interactions due to the vertical or lateral defects disturb the biomechanical stability of the RBC membrane. Therefore, when a defective RBC experiences large biomechanical deformation, it may be unable to recover its original shape even though the stretching force causing the shape change is removed. Here, we examine the elastic relaxation of H-RBC and defective RBC when Fs = 110 pN is turned off. The simulation results are shown in Fig 7 and more details are available in the video clips in the Supporting Information. In Fig 7A, we find that both the lipid bilayer and the cytoskeleton are deformed simultaneously under stretching-relaxation cycles and their DA and DT values approach the initial values at equilibrium state, which demonstrates that a normal RBC with high elasticity can recover its original shape. The recovery process of the RBC can be characterized by the dynamic recovery expression, R (t), which is described by an exponential decay [67], R (t) = (λ - λ ∞) (λ 0 + λ ∞) (λ + λ ∞) (λ 0 - λ ∞) = exp - (t - t 0) t c, (21) where λ = D A D T, λ0 and λ∞ correspond to the ratios at release and recovery states; t0 is the time when Fs is turned off and tc is the characteristic time. Using the relaxation results from Fig 7A into R (t), we obtain the best fitting decay as illustrated in S4 Fig and estimate the characteristic time, tc, for both lipid bilayer and cytoskeleton, i. e. , tc ≈ 0. 12 s. This value falls within the range of 0. 1–0. 3 s measured from the shape relaxation of H-RBCs in experiments [65,67] and also computed in simulations [33]. However, for the defective RBC with weakened bilayer-cytoskeleton interactions, once the detachment of the lipid bilayer from the cytoskeleton occurs, the DA values of both lipid bilayer and cytoskeleton are always larger than those at equilibrium state, even for a sufficient relaxation time, which indicates that the defective RBC loses its ability to recover its original shape. We note that not only the lipid bilayer but also the cytoskeleton has a permanent deformation because of the slow response distortion as shown in Fig 7B and the S2 video. In summary, our simulation results show that the elastic deformations of Pf-RBCs match those obtained in optical tweezer experiments for different stages of intraerythrocytic parasite development. We found that the developed knobs on the cell membrane can effectively stiffen the lipid bilayer and the spectrin network, leading to a decrease in RBC deformability. We also investigated the effects of bilayer-cytoskeleton interactions and our simulation results reveal that these interactions play a key role in the determination of cell membrane biomechanical properties. We analyzed the biomechanical response of the RBCs during loading and upon release of the tensile force, in order to mimic the stretching and compression that RBCs experience as they pass through small capillaries in the microcirculation. Our results indicate that H-RBCs with normal bilayer-cytoskeleton interactions recover their original discocyte shape when the stretching force is turned off. However, defective RBCs with weakened bilayer-cytoskeleton interactions lose their ability to recover their original shape due to the irreversible membrane damage. Overall, our findings demonstrate that the two-step multiscale framework we have developed, combining coarse-grained molecular dynamics (CGMD) and dissipative particle dynamics (DPD), can be used to predict the altered biomechanical properties of RBCs associated with their pathophysiological states.
summary: Red blood cells (RBCs) and their mechanical properties play a crucial role in the dynamic and rheological behavior of blood in normal and disease states. However, the precise determination of RBC membrane properties is hard to be achieved experimentally. To this end, accurate numerical modeling can be used to provide valuable information for quantifying the biomechanical properties of RBCs. In this paper, we have developed and validated a two-step multiscale framework for RBC modeling, by performing molecular dynamics simulations to compute the shear modulus, bending stiffness and network parameters of a small RBC patch, which we then use as input to dissipative particle dynamics simulations to predict the stress field and morphology of defective RBCs, including Plasmodium falciparum-infected RBCs as well as RBCs in hereditary spherocytosis and elliptocytosis.
keywords: classical mechanics
  spectrins
  membrane proteins
  damage mechanics
  cellular structures and organelles
  cytoskeleton
  lipids
  proteins
  deformation
  membrane-associated cytoskeletal proteins
  cell membranes
  physics
  biochemistry
  cytoskeletal proteins
  biochemical simulations
  cell biology
  integral membrane proteins
  lipid bilayer
  biology and life sciences
  physical sciences
  computational biology
title: MD/DPD Multiscale Framework for Predicting Morphology and Stresses of Red Blood Cells in Health and Disease