1 Introduction

Biomembrane is a vital component of living cells since it helps preserve the cell integrity^[1-2]. It plays a key role on various functions of cells^[3]. The lipid vesicles formed by budding off from a cell membrane can transport the biological materials between cells. These transport vesicles drop off their loads through fusion with the membrane of the target cells. Cell membrane is also essential in multiple intercellular processes, e.g., cell differentiation, cell-cell adhesion, and cell migration^[1]. Biomembrane consists of various lipid molecules and proteins. A typical lipid molecule is constituted by two hydrophobic hydrocarbon chains and a hydrophilic polar head, such that the lipid molecules can spontaneously aggregate into bilayer and then form vesicles in aqueous environment because of the hydrophobic effect. Biomembrane behaves similar to two-dimensional (2D) fluid. Despite that the thickness of membrane is only about 5 nm, the size of lipid vesicles can be up to micrometers^[2].

Red blood cells (RBCs) are unique amongst eukaryotic cells since they bear no nucleus, cytoplasmic structures, or organelles^[4]. Therefore, the structural properties are linked to the cell membrane. In addition to the lipid bilayer and integral proteins, RBC membrane possesses a 2D cytoskeleton tethered to the lipid bilayer. The RBC membrane cytoskeleton has a 2D six-fold structure, and is made of the spectrin tetramers connected at the actin junctional complexes, forming a 2D sixfold structure. The cytoskeleton is connected to the lipid bilayer via the "immobile" band-3 proteins at the spectrin-ankyrin binding sites and the glycophorin protein at the actin junctional complexes^[4]. In RBC membranes, the stiffness and elasticity of the RBCs arise primarily from the cytoskeleton. Given this particular membrane structure, RBCs have remarkable deformability, which allows them to undergo repeated severe deformation when they are traversing through small blood vessels and organs^[5].

In spite of the significant advances made in computing power in the past few decades, it is still computationally prohibitive or impractical to perform such atomistic simulations on cell membrane at time scales and length scales that the obtained results can be used to directly compare with typical laboratory experimental studies. However, atomistic simulation techniques are limited by the number of the involved atoms/molecules, typically 10⁴~ 10⁸ corresponding to a length scale on the order of tens of nanometers. Moreover, the maximum time step in atomistic simulations is limited by the smallest oscillation period of the fastest atomistic motion in a molecule, which is typically several femtoseconds (10^-15 s). Most atomistic membrane models are limited to the study of only a few hundred lipids or single proteins for a period of a few nanoseconds because of the excessive computational cost^[6-11].

At the opposite end of the length scale and time scale spectra, continuum models^[12-14] were developed to simulate the morphologies of vesicle and to measure the thermal fluctuations of fluid membranes at much larger length scales^[15-16]. Boundary integral method^[17-20] and immersed boundary method^[21-24] are two popular ways to simulate the RBCs based on the assumptions that the RBC membranes and embedding fluids are homogeneous materials. These continuum-based RBC models facilitate the study of single cell dynamics, e.g., tank-threading^[25], RBC passing splenic slit^[26], and blood flow on macroscopic length and time scales^[27-28]. Although the continuum-based RBC models provide an accurate description of the RBC deformation at cellular level, they are not able to describe mesoscopic-and microscopic-scale phenomena, e.g., representing the structural defects in the RBC cytoskeleton, and capture the topological changes of membranes or RBCs, e.g., self-assembling, fusion of lipid vesicles, and membrane budding.

The limitations of atomistic and continuum methods have motivated a continuous effort on developing the coarse-grained (CG) particle models that bridge atomistic and continuum models^{[10, 29-51]}. CG particle models can substantially simplify the atomistic dynamics by eliminating the fast degrees of freedom while preserving the mesoscopic structures and properties of the simulated system^[52]. It has been applied broadly in polymer dynamics^[53-61] and biological fibers^[62-68]. In the following sections, we will review the previously developed biomembrane models, RBC membrane models, and RBC whole cell models, and demonstrate their power and versatility in simulating lipid self-assembly vesicle fusion and membrane budding.

2 CG particle models 2.1 Biomembrane models

CG membrane models treat a group of atoms in a single phospholipid molecule into CG particles since phospholipids are the most abundant membrane lipids, which have a polar head group and two hydrophobic hydrocarbon tails. An atomistic model for a phospholipid molecule is shown in Fig. 1(a). As a result, each lipid molecule is modeled as a short chain of CG particles. The levels of coarse-graining in different CG membrane models are determined by the targeted problems of the simulation. Markvoort et al.^[45] introduced a CG molecular dynamics (CGMD) model for lipid molecules, where two chains of four CG particles represented the two hydrophobic tails, and one chain of four CG particles represented the lipid headgroup (see Fig. 1(b)). The authors used two types of lipids to introduce the asymmetry in the lipid bilayers. At a high level of coarse-graining, Wang and Frenkel^[43] proposed a membrane model, where each lipid molecule was represented by only three CG particles with one particle being the hydrophilic head and the other two being the hydrophobic tails (see Fig. 1(c)). A number of CGMD membrane models were developed, following the similar strategy^{[30-31, 33-35, 38-39, 44, 69]}. To achieve a high level of coarse-graining, Yuan et al.^[70] coarse-grained a group of lipid molecules to one single CG particle, and built a one-particle thick lipid bilayer model (see Fig. 1(d)). The hydrophilic and hydrophobic properties of the lipid molecules were represented by a directional vector.

CG particle models for biomembranes can also be grouped into explicit solvent models and implicit solvent models^{[29, 44, 48]}. In explicit solvent models, the hydrophobic interactions are implemented between the lipids and the solvent particles to maintain the stability of the lipid membrane^{[34, 42]}. Explicit solvent models frequently employ dissipative particle dynamics (DPD), which is a very efficient method representing the large volume of the solvents with a soft bead and thus significantly accelerate the computations^[71-73]. Three types of forces are applied in traditional DPD models, i.e., conserved soft repulsion force, pairwise dissipation forces, and pairwise random forces. The temperature of the DPD system is preserved through balancing the dissipation and random forces. The momentum of DPD particles is also preserved, and thus provides accurate hydrodynamics. Besides increasing the length scale, DPD simulations can apply a larger timestep because of the soft repulsion conservative forces. However, explicit solvent methods can dramatically boost the computational expense since the solvent particles are required to fill a three-dimensional (3D) simulation box. The number of solvent particles can grow quickly when the size of the simulation box is increased. Implicit solvent models can overcome these disadvantages, and thus can achieve larger length and time scales^{[44, 74]}. In implicit solvent models, instead of modeling the solvent particles, the effect is represented by effective multi-body interaction potentials based on either the local particle density^{[30, 33, 43]} or the application of different pair-potentials between the particles representing the hydrophobic tails and the particles representing the hydrophilic head of the lipids^{[35, 38-39]}. The development of solvent-free CGMD lipid membrane models requires the representations of hydrophobic interactions by interacting the potentials so that the hydrophobic behavior of the lipid molecules can be captured. One group of models applied a pair potentials, which are softer than the Lennard-Jones (LJ)-type potentials, between the tail beads to preserve the stability of the fluid membrane, and then implemented LJ-type potentials for all the other inter-bead interactions^[39]. Along this line, Drouffe et al.^[30] and Noguchi et al.^[46] developed one-particle thick CGMD models to simulate biomembranes. In Drouffe's model, an LJ-type pair potential is used to simulate the interactions between the CG particles. This potential is a function of the distance between the particles and their directionality. Noguchi et al.^[46] introduced a multibody potential that did not require the directionality of particles. Yuan et al.^[70] introduced a membrane model with the application of a soft-core potential in order to illustrate the particle self-diffusion.

2.2 Lipid self-assembling and membrane fusion

The self-assembly processes of lipids from a randomly dispersed state to an ordered bilayer structure have been captured by CG particle models^{[36, 45, 75-77]}. Starting from a random dispersion, lipids rapidly aggregate into micelles and small bilayers referred to bicelles. Subsequently, these small aggregates merge into a large disklike bilayer. In order to minimize the line tension arising at the bilayer edge, this bilayer tends to seal and form a vesicle by gradually encapsulating water (see the self-assembly process of lipids in Fig. 2). In fact, such bent lamellar structure is essential to cell membrane related processes, especially for exo-endo-cytosis. For entropic reasons, a bilayer does not stay at a simply planar state, it exhibits thermal fluctuations (undulations), and its bending rigidity is responsible for manipulating the fluctuations in the average position of the bilayer surface^[78]. Although a vesicle can persist for a long time, it is merely in a metastable state instead of a thermodynamic equilibrium state. This is attributed to the implicitly higher pressure and chemical potential of water in the interior domain than those in the exterior domain of a vesicle^[76]. In addition, the smaller the vesicle is, the larger the pressure difference exists, and thus, small vesicles are thermodynamically less stable against fusion than large vesicles.

The lipids in bilayer membrane can be either at a solid-ordered state or a liquid-disordered state. It is known that lipid bilayers exhibit a phase transition from a gel phase to a liquid phase when the temperature increases^[79-80]. For some lipids, e.g., DPPC and dimyristoylphosphatidylcholine (DMPC), there could exist an intermediate state (ripple phase) between the gel and liquid phases^[81] (see Fig. 3), and the ripple phase can be either symmetric or asymmetric, depending on the cooling rate from the liquid state down to a crystalline structure^[82]. The phase behavior of a lipid bilayer affects its membrane structure and mechanical properties, and therefore the functioning of biological membranes is also temperature-dependent. In addition to the change of molecular volume (or area per lipid) and heat capacity with temperature as shown in Fig. 3, several DPD simulations^[83-85] have demonstrated that the membrane thickness of the lipid bilayer (h) and the orientation of lipid tail (S) can be functions of the temperature (T), and the main transition temperature (T_m) is located right at the inflection point of h(T) and S(T). CG membrane models can facilitate the estimation of membrane physical properties as well. For example, the water permeability of the modeled membrane can be predicted by the release rate of the lumenal substance. Wu et al.^[84] showed that the water permeability of bilayer membrane grew exponentially with increasing T for T > T_m, and the obtained simulation results agreed well with the experimental measurements^[86]. Besides, the fluctuation spectra of the lipid membrane can provide the information of the bending modulus. The intensity I = 〈u²〉 of the undulations follows the q⁴ behavior in the long wavelength regime, i.e., where u is the fluctuation amplitude, A is the membrane area, q is the wave number, and k_c is the bending modulus. By fitting the equation to the simulation data, Marrink et al.^[87] obtained k_c= 4 ± 2 × 10^-20 J, which is consistent with the experimentally measured 5.6 ± 0.6×10^-20 J for DMPC lipids^[88].

Vesicle fusion is a principal function of cells to communicate and execute a series of bio-physical activities from vesicular trafficking to cell-cell fusion^[91]. Besides, lipid fusion has been utilized to several technological applications such as drug delivery^[92] and the surface modifications of supported lipid bilayer for biosensors^[93-94]. Many in vitro fusion experiments have been performed with lipid vesicles and the dynamic process is in milliseconds and micrometer length scale^[95-96]. However, the molecular structures and transitions during fusion cannot be captured and visualized by currently experimental techniques. Numerical framework can allow itself to elucidate the vesicle fusion in more detail. Several pathways have been proposed to describe the fusion processes of lipid vesicles with the help of CG modeling, and the stalk-pore hypothesis has been widely used^{[37, 97-101]} (see Fig. 4 (left)). First, the emergence of a neck-like structure (usually regarded as a stalk) is initiated by splayed and tilted lipids that connected the outer monolayers of two approaching vesicles. Then, the stalk grows quickly to form a hemifusion diaphragm (HD), which has the thickness equal to a bilayer thickness. The lumenal contents enclosed by each vesicle still remain separated in the HD stage. Finally, the pores are generated to end up the fusion process. By using a CG lipid model, Marrink and Mark^[37] found that the alternations of lipid composition can result in distinct speeds for the stalk formation and the opening of the fusion pore. Another pathway for lipid fusion is called the direct stalk-pore model^[100-104] (see Fig. 4 (right)). Similar to the original stalk model, the modified one is defined when there is no evident HD stage during the fusion. It occurs usually when vesicles are under high membrane tension and the pores are initiated from the radially expanding stalk. Both the original and modified stalk models are applicable to the fusion of vesicles formed by phospholipids and lipid-like copolymers.

2.3 RBC membrane models

The aforementioned CG particle membrane models only simulate the lipid bilayer without considering the membrane cytoskeleton. Thus, their applications are limited in the study of the biological problems regarding RBC membrane, e.g., interactions between the lipid bilayer and cytoskeleton and between the cytoskeleton and diffusing membrane proteins. To model the membranes for RBCs, Li and Lykotrafitis^[105] introduced a two-component RBC membrane model by combining a cytoskeleton model with the lipid bilayer model. The cytoskeleton model was comprised of particles representing actin junctional complexes, which were connected by the weighted linear combination (WLC) potentials to represent the spectrin filaments. Although this two-component model exhibits membrane shear modulus comparable with the experimental measurements, the implicit representation of the spectrin filaments by the WLC potential did not allow them to consider the interactions between the lipid bilayer and cytoskeleton. Therefore, Li and Lykotrafitis^[106] extended this two-component RBC membrane model by representing the spectrin filaments explicitly by CG particles^[107]. This modified membrane model describes the RBC membrane as a two-component system with cytoskeleton and lipid bilayer. The lipid bilayer is represented by three types of CG particles, and the cytoskeleton is represented by two types of CG particles (see Fig. 5). The cytoskeleton consists of the hexagonal spectrin network and actin junctions. The actin junctions represented by the black particles (see Fig. 5) are connected to the lipid bilayer via glycophorin (yellow particles). 39 spectrin particles (grey particles) connected by unbreakable springs are used to represent the spectrin filament. The spring potential is represented by

with the equilibrium distance defined by

where L_max is the spectrin contour length, and The spectrin chains are connected to the band-3 particles (white particles) in the middle, and are linked to the actin junction via the spring potential, i.e.,

at two ends. The equilibrium distance is 10 nm. The spring constant k₀ can be calculated by

An LJ potential is applied between the spectrin particles, i.e.,

(1)

where r_ij is the interparticle distance. The energy unit and the length unit are ε and σ, respectively. The cutoff distance R_{cut, LJ} is selected to be the equilibrium distance between the spectrin particles.

The lipid bilayer and transmembrane proteins are represented by three types of CG particles (see Fig. 5). A cluster of lipid molecules is grouped into a red particle, which has a diameter of 5 nm. The yellow particles denote the glycophorin proteins. One third of the band 3 particles (white particles) denoting the band 3 dimmers are linked to the cytoskeleton. Two thirds of the band 3 particles (green particles) simulate the mobile band 3 proteins. Both the translational and the rotational degrees of freedom (x_i, n_i) are assigned to the CG particles, which form the lipid bilayer and transmembrane proteins. n_i and x_i are the direction and position vectors of the particle i, respectively. x_ij= x_j-x_i is defined as the interparticle distance vector between the particle i and the particle j. Correspondingly, r_ij ≡ |x_ij| is the distance, and is a unit vector. The CG particles, forming the lipid membrane and membrane proteins, interact through a pair-wise soft potential as follows:

(2)

(3)

where u_A(r_ij) and u_R(r_ij) are the attractive and repulsive parts of the pair potential, respectively. α is a parameter that can tune the attractive part of the potential. Figure 6(a) illustrates that when A(α, a) = +1, the energy well of the applied potential is wider than the LJ6-12 potential. Therefore, facilitating the CG particles pasting each other can enhance the fluidity of the membrane. The function

can adjust the energy well of the potential to regulate the fluidity of the membrane (see Fig. 6(b)). In the simulations, α is chosen to be 1.55, and the cutoff distance R_{cut, mem} is 2.6σ. The parameters α and R_{cut, mem} are selected to maintain the fluid phase of the lipid bilayer. The directional vector carried by each CG particle is important for the self-assembly of the membrane in the absence of solvent particles and lipid chains. The detailed information about the applied potentials and the selection of the potential parameters can be found in Refs. [105] and [106].

This modified two-component membrane model can simulate the molecular structures of the normal or diseased RBC membranes, and thus has been used to study the membrane stiffness of the protein diffusion and membrane vesiculation (see Fig. 7) in healthy and pathological conditions^[108-111].

However, it is computationally expensive to simulate large membrane domains such as an entire RBC or the whole blood involving large numbers of RBCs. To increase the length scale of the modified two-component membrane model, Tang et al.^[112] recently developed OpenRBC, i.e., a multi-thread CGMD code, which was capable of simulating an entire RBC with explicit representations of lipid bilayer and cytoskeleton by multiple millions of CG particles with a single shared memory commodity workstation (see Fig. 8).

2.4 RBC whole cell models

Although OpenRBC allows us to simulate an entire RBC at protein resolution, it is still computationally prohibitive to simulate RBC suspensions and study the blood rheology. Thus, a more efficient RBC model with higher level of coarse-graining is required. Discher and his co-workers^[113] developed a CG particle model for the spectrin cytoskeleton of RBC membranes to study the elasticity of RBC in micropipette aspiration. In this spectrin-level RBC model, a RBC is treated as a 2D canonical hexagonal network of CG particles (see Fig. 9(a)). The neighboring CG particles are connected via a WLC potential to represent a spectrin filament. The bending rigidity induced by the lipid bilayer is represented by a bending potential applied between two neighboring triangles. Moreover, the constrains on the cell volume and surface area are implemented to model the incompressibility of the cytosol and lipid bilayer. Li et al.^[114] extended this model by introducing a dynamic spectrin network spontaneous curvature for the lipid membrane and in-plane shear energy relaxation. Pivkin and Karniadakis^[115] introduced the hydrodynamic effect into this elastic RBC model by combining it with DPD. More importantly, this RBC model is further coarse-grained from 23 867 CG particles to 500 CG particles while the elastic properties are still preserved. This DPD-based whole cell model is computationally very efficient and it can be used not only in the simulations of single cell dynamics^[116] and morphological changes of diseased RBC (see Fig. 9(b))^[117] but also the rheology of blood flow^[118-120] and blood flowing through micro-devices (see Fig. 9(c))^[121]. However, it is not capable of modeling the interactions between the cytoskeleton and the lipid bilayer, which is frequently essential. To overcome this issue, a RBC model made of two layers of 2D triangulated networks is developed to study the lipid bilayer and cytoskeleton interactions in healthy and diseased RBCs^[122-125]. Both one-component and two-component whole-cell models consider the bending energy (V_b), the elastic energy (V_s), the constraints of the enclosed volume (V_v), and the fixed surface area (V_a) of a RBC. In the two-component whole-cell model, the bilayer-cytoskeleton interaction energy (V_int) is additionally considered. V_s represents the elastic spectrin network, and is given by

(4)

where l_m is the maximum spring extension, while l_j is the length of the spring j. x_j = l_j/l_m. p represents the persistence length. k_BT designates the energy unit. k_p is the spring constant of the potential. n is a specified exponent. V_b results from the bending stiffness of the lipid bilayer, and it is expressed as follows:

(5)

where k_b is the bending constant. θ_j is the instantaneous angle between two neighboring triangles. θ₀ is the spontaneous angle between the two neighboring triangles. V_a and V_v are implemented to represent the incompressible lipid bilayer and interior fluid. The potentials for constraining the surface area and volume are given by

(6)

(7)

where A₀ is the triangle area. N_t is the number of triangles in the membrane network. k_d, k_a , and k_v are the constraint coefficients for the local area, global area, and volume, respectively. The specified total area and volume are and respectively. V_int represents the interaction between the lipid bilayer, and it is expressed as a summation of the harmonic spring potentials as follows:

(8)

where N_bs is the number of the bond connections between the lipid bilayer and the cytoskeleton. k_bs is the spring constant. d_jj' measures the distance between the vertex j on the cytoskeleton and its projection point j' on the lipid bilayer. d_{jj', 0} is the initial distance between the vertex j and the point j'.

3 Summary

Although atomistic simulations provide the detailed structural properties of lipid membranes, they are limited to studying small membrane patches due to the high computational expense. Many interesting biological phenomena cannot be accessed by these simulations when individual atoms or lipid molecules are resolved. Continuum-based methods can model big membranes or entire RBCs. However, they cannot address the problems concerning the detailed structural changes in RBC membranes and the specific protein defects in diseased RBCs. Moreover, typical continuum models are not suitable to simulate the dramatic topology changes such as membrane budding, fusion, and self-assembly of lipids to vesicles. CG particle models are bridging in certain aspects of the MD and continuum methods. First, the CG particle-based method is similar to the MD method. However, when each CG particle represents a lump of molecules, the atomistic simulations will be simplified by eliminating the fast degrees of freedom while preserving the mesoscopic structures and correcting the dynamics properties of the membrane or RBCs at large spatial and temporal scales. Second, the CG particles can be treated as a mesh-free representation of continuum. Particularly, the CG particle model can be used to simulate the dynamical processes involving significant topological changes, e.g., self-assembly lipid particles, vesicle fusion, and nanoparticles transport through the erythrocyte membrane. It is worth noting that in CGMD simulations, the employed length scale and time scale do not have an immediate correlation with a physical system. Such a correspondence can be established only through the comparison with a natural procedure. An emerging particle-based method, i.e., the smoothed dissipative particle dynamics (SDPD)^[126-128], has also been used to model the rheology of red and white blood cell behaviors so as to elaborate the interaction of the surrounding fluid and membrane^[129-130]. The SDPD method has advantages over conventional DPD, including well-defined physical scale of particles, direct inputs of transport properties, and arbitrary equation of state. As an alternative of the DPD, the SDPD can be used to study the transport phenomenon involved in blood flows, e.g., drug delivering and cell loading in saline solution. Computational simulations have become increasingly important to improve our understanding in the rheology and dynamics of RBCs, especially in blood diseases. The CG membrane and RBC models developed in the last decade have been successfully implemented to uncover various membrane functions and characteristics. For the complete application picture of the CG particle models for membranes and RBCs under physiological and pathological conditions, one can refer to Refs. [131], [132], and [133] for the details.