1. Introduction

Biological cell membranes are dynamic, adaptive, stimuli-responsive multi-component soft materials which separate the cytosol from the extracellular environment, and participate in vital functions, for example intracellular and extracellular traffic, sensing and cell signaling [1]. Soft materials are mesoscopic supramolecular structures assembled from the self-organization of nanoscopic building blocks such as lipid molecules whose predominant behavior occurs at room temperature. In addition to molecular-scale interactions between the nanoscopic units, the kinetics and thermodynamics of the structure, dynamics and function of cell membranes is dictated by hydrodynamic interactions which are predominant at the mesoscopic length scales. The complex interplay between these key factors across multiple scales determines the structure of the membrane and thereby the entire cell, and its responses to changes in the external environment such as temperature, pH and presence of ions, or binding or catalytic events at the membrane interface. The dynamic nature of this multi-component smart active biomaterial enables the modulation of membrane tension and thereby, its mechanical properties to facilitate various physiological processes. The molecular geometry and rigidity of the species will determine their packing and organization in the bilayer membrane [2,3], thereby influencing the mechanical and the interfacial properties of the membrane. Different types of cell membranes can be characterized by the composition of their constituent amphiphilic lipid or sterol species [4]. Hence, the innate responsive nature of cell membranes can inspire the design of robust soft materials with highly adaptive mechanical and interfacial properties. In addition, the ability to target certain cells by way of their specific membranes has lead to the development of new areas of medical research for applications in drug delivery, sensing and imaging [4,5,6,7,8]. Based on the surface chemistry of the cell membrane, nano-sized particles with active payloads can target specific cell membranes [9]. Therefore, there is a critical need to underst and the effect of the molecular geometry of the constituent species on the properties and functions of cell membranes.

There are challenges to using experimental approaches such as X-ray diffraction, NMR spectroscopy, and neutron diffraction to study the interactions of nano-sized particles with cell membranes [10], and the role of the individual molecular species on the interaction process. Computer simulations [11] can be used to elucidate the role of the membrane composition on the interaction between foreign particles and cell membranes, or the response of cell membranes to changes in external cues, through the use of suitable models and techniques. Earlier computational studies on lipid bilayers have adopted techniques such as Monte Carlo [12,13,14], and Molecular Dynamics [15,16,17,18,19,20]. Given the importance of hydrodynamic interactions on the properties of cell membranes, we adopt a Molecular Dynamics-based mesoscopic simulation technique which bridges the atomistic and continuum scales, and captures the hydrodynamic behavior of the system over extended time scales. This method is entitled Dissipative Particle Dynamics [10,20,21,22,23,24,25,26,27,28,29,30,31] and uses soft-repulsive core interaction models between the various particles which are coarse-grained representations of a group of atoms. This approach enables the study of complex dynamical and structural properties, and biological processes [32,33] which are not easily resolved by Molecular Dynamics [15,16,17]. In addition, DPD has been used to demonstrate the dependence of the phase behavior, the mechanical and interfacial properties of a lipid membrane on the membrane composition, hydrocarbon tail saturation, head group architecture, and hydrocarbon chain length [9,34,35]. Other investigations have examined the response of bilayer membranes to external stimuli [21,25,36]. Therefore, this approach can aid in the design and development of bio-inspired cells composed of phospholipid vesicles for targeted delivery of therapeutic compounds [11,37,38,39,40,41,42]. We would like to note that the DPD method has been used to provide insight into the interactions between nanoparticles and block copolymers [43,44].

Via the Dissipative Particle Dynamics (DPD) approach, we will investigate the role of the architecture of amphiphilic lipids on the interfacial, structural and dynamical properties of lipid bilayer membranes [21,25,26,27,28,29]. We focus on four lipid architectures with differing molecular geometries. The first lipid model represents DPPC (Dipalmitoylphosphatidylcholine) which is a cylindrical-shaped molecule with a large head group and two hydrocarbon tails [20,25,26,27,28,29,31]. The other three lipid models represent DMPC (Dimyristoyl-sn-glycero-3-phosphocholine) which is an inverted wedge-shaped molecule with a small head group and two-hydrocarbon tails [45]. The three models of DMPC have different hydrocarbon tail lengths. In addition, we have investigated a coarse-grained model for Cholesterol which has a small head group and a single rigid hydrocarbon tail with a steroid ring [46]. We have examined the interfacial properties of the Cholesterol molecule as a function of the inter-molecular spacing. Finally, we demonstrate the design of a stable mixed bilayer membrane composed of DPPC and cholesterol molecules using their average area per molecule corresponding to a tensionless bilayer. Our approach can be adopted for the conception and design of multi-component bilayer membranes with tunable mechanical and interfacial properties.

2. Materials and Method

2.1. Dissipative Particle Dynamics

DPD is a mesoscopic MD-based simulation technique that uses soft-sphere coarse-grained (CG) models to capture both the molecular details of the system components and their supramolecular organization while simultaneously resolving the hydrodynamics of the system over extended time scales [9,10,11,21,47,48]. In order to capture the dynamics of the soft spheres, the DPD technique integrates Newton’s equation of motion via the use of similar numerical integrators used in other deterministic particle-based simulation methods [21,22]. The force acting on a soft sphere i due to its interactions with a neighboring soft sphere j (j ¹i) has three components: a conservative force, a dissipative force and a r and om force, which operate within a certain cut-off distance r_c from the reference particle i. These forces are pairwise additive and yield the total force acting of particle i, which is given by ${{\bf{F}}_i} = \sum\limits_{j \ne i} {{{\bf{F}}_{c,ij}} + {{\bf{F}}_{d,ij}} + {{\bf{F}}_{r,ij}}}$. The soft spheres interact via a soft-repulsive force (${{\bf{F}}_{c,ij}}={{a}_{ij}}\left( 1-{}^{{{r}_{ij}}}\!\!\diagup\!\!{}_{{{r}_{c}}}\; \right){{\bf{r}}_{ij}}$, for r_ij < r_c and and F_c,ij=0 , for r_ij ≥ r_c), a dissipative force (${{\bf{F}}_{d,ij}}=-\gamma {{\omega }^{d}}\left( {{r}_{ij}} \right)\left( {{{\bf{\hat{r}}}}_{ij}}\bullet {{\bf{v}}_{ij}} \right){{{\bf{\hat{r}}}}_{ij}}$) and a random force (${{\bf{F}}_{r,ij}}=-\sigma {{\omega }^{r}}\left( {{r}_{ij}} \right)\theta {{{\bf{\hat{r}}}}_{ij}}$), where w^d (r) =[w^r (r) ]²= (1-r) ² (for r<1) , and ${\sigma ^2} = 2\gamma {k_B}T$. a_ij is the maximum repulsion between spheres i and j, v_ij = v_i - v_j is the relative velocity of the two spheres, r_ij = r_i - r_j, r_ij = |r_i - r_j|, ${{\bf{\hat r}}_{ij}}$= r_ij/r_ij, r = r_ij/r_c, γ is viscosity related parameter used in the simulations, q_ij (t) is a r and omly fluctuating variable from Gaussian statistics, w^d and w^r are the separation dependent weight functions which become zero at distances greater than or equal to the cutoff distance r _c. Each force conserves linear and angular momentum. Since the local momentum is conserved by all of these three forces, even the small systems exhibit hydrodynamic behavior [21]. The constraints imposed on the r and om and dissipative forces by certain relations ensure that the statistical mechanics of the system conforms to the canonical ensemble [21,22]. The relation between the pair repulsion parameter a_ij and the Flory interaction parameter c for a bead number density r = 3 is given by [48].

As shown in Figure 1 (a) - (e), the individual lipid and cholesterol molecules are represented by bead-spring models. Two consecutive beads in a chain are connected via a bond that is described by the harmonic spring potential E_bond=E_bond ((r-b)/r_c)², where K_bond is the bond constant and b is the equilibrium bond length. The constants, K_bond and b are assigned to the values of 64 and 0.5, respectively [20,24,25,26,27,28,29,85,49,50]. For the lipid molecules, a weaker bond is inserted (K’_bond = 16) between the first beads on the two tails to ensure that the tails are positioned in the same direction. The three-body stiffness potential along the lipid tails has the form E_angle=K_angle (1+cosθ) where % $\theta$is the angle formed by three adjacent beads. The coefficient K_angle is set to be 20 in our simulations. This stiffness term increases the stability and bending rigidity of the bilayers [47].

The soft repulsive pair potential parameters for the lipid and cholesterol hydrophilic and hydrophobic beads were selected to capture their amphiphilic nature. The interaction parameters between the like components, ${a_{ij}}$, are based on the property of water [21]. The repulsion parameter between two beads of the same type is set at a_ii = 25 (measured in units of ${k_B}T$) which is based upon the compressibility of water at room temperature [21] for a bead density of r = 3. The soft repulsive interaction parameter a_ij between hydrophobic and hydrophilic beads is set at ${a_{ij}} = 100$, and is determined by using the Flory-Huggins interaction parameters, c, as ${a_{ij}} = {a_{ii}} + 3.496\chi$, for r = 3.

The soft repulsive interaction parameters between the head (h), tail (t) beads of lipid molecules, and the solvent (s) beads are assigned the following values: a_ss = 25, a_hh = 25, a_tt = 25, a_hs = 25, a_ts = 100 and a_ht = 100. The values of the inter-lipid species head-head a_h1h2 and tail-tail a_t1t2 soft repulsive interaction parameters will mimic mixtures of lipid species with different head or tail groups. The physical properties of the lipid molecules are summarized in Table 1. These parameters are selected to model the effective distinct chemistry of the molecular species, thereby capturing the differences in the melting temperature of the individual species [13,51,52]. This approach enables us to develop a simple representation of mixtures composed of lipids with two hydrocarbon tails.

In our simulations, the respective characteristic length and energy scales are r_c and k_BT. We run our simulations at a reduced temperature T = 1 which corresponds to the temperature at which the DPPC bilayer is in the fluid state. The transition temperatures for pure DPPC and DMPC bilayers are respectively, 41°C and 23°C [53,54], as shown in Table 1. As a result, our characteristic time scale can be described as $t = \sqrt {mr_c^2} /\mathop k\nolimits_b T$. Finally, σ = 3 and Δt = 0.02τ are used in the simulations along with the total bead number density of ρ = 3 and a dimensionless value of r_c = 1 [24].

To develop a correspondence between the reduced and physical units, we relate the experimental measurements of the area per lipid [54] for a tensionless membrane and the diffusion coefficient of a lipid molecule in a membrane in the fluid state [55,56]. We obtain the characteristic length scale for our model through the comparison of experimental measurements of the interfacial area per lipid of the corresponding bilayer with similar measurements from our simulations. To compute the average area per lipid, the vesicle is divided into many small rectangular patches so that each patch can be treated as a bilayer membrane. By summing the areas of all the patches and then averaging them over the systems total number of lipid molecules in the whole system, the average area per lipid for the vesicle bilayer is computed. The diffusion coefficient of the lipid molecule in the simulations can be found by tracking the mean squared displacements of 10 lipid molecules in a vesicle bilayer. We use the re lation $\frac{{\partial \left\langle {{r^2} (t) } \right\rangle }}{{\partial t}} = 2dD$to relate the diffusion coefficient D to the mean square displacement of a particle in a time intervalt [23]. The variable d is the dimensionality of the system that is given to be 3 for our system. We calculate the diffusion coefficient D using the slope of the time evolution of the mean square displacement. Table 1 lists the physical correspondence of the reduced units for the four lipidmodels. We note that the lipid model A is based upon the DPPC molecule and lipid model B represents the DMPC molecule.

2.2. Lipid Models and Architecture

We explore the effect of the lipid molecular architecture in the bilayer properties. We study four models of double-tail lipids with variations in the molecular packing parameter and the hydrocarbon chain length [2]. The molecular packing parameter for each lipid architecture depends upon its area per lipid, hydrocarbon chain volume, and the chain length [2]. Lipid model A has a cylindrical shape with a packing parameter around 1, and is comprised of three head group beads (which are hydrophilic in nature) and two hydrocarbon tails with three hydrophobic beads each (see Figure 1). Lipid models B, C and D have an inverted wedge shape and a molecular packing parameter greater than 1, with successively longer hydrocarbon chain lengths (see Figure 1) The head group architectures for models B, C and D are identical and are comprised of three hydrophilic beads. Each hydrocarbon tail is comprised of three, four and five hydrophobic beads respectively, for lipid models B, C and D. In summary, lipid models A and B differ in their head architecture; lipid models A, C and D differ in their head and tail architectures, and lipid models B, C and D differ in their hydrocarbon tail length.

2.3. Cholesterol Model and Architecture

The cholesterol molecule is modeled by two hydrophilic head beads and seven hydrophobic beads that are organized in a ring and a short tail (see Figure 1) To model the sterol ring, we introduce bonds between diametrically opposite beads encompassing the ring, in addition to the bonds between successive beads at the periphery of the ring. A schematic of the bond topology of the molecule is overlaid on the coarse-grained representation (see Figure 1.) For simplicity, we have maintained the harmonic bond potential K_bond = 64 and equilibrium bond length r₀ = 0.5. The angle potential energy between three consecutively bonded beads is set at K_angle = 20.

3. Results

3.1. Characterization of Lipid Bilayer Membranes

We use a cubic simulation box of dimension 20r_c with periodic boundary conditions along all three directions, with 24, 000 beads. We begin with a preassembled lipid bilayer membrane with a predetermined average area per lipid [32,57], equilibrate the membrane for a time interval of 102, 000t, and measure interfacial properties over a duration of 20, 000t. We repeat this process for each of the molecular models. Following the equilibration phase, we measured the interfacial properties for the lipid bilayer membranes as a function of the area per lipid. We measure contributions to the interfacial properties by computing the pair, bond and angle energies of the lipid bilayer membranes composed of the different models. The enthalpic component of the free energy of a bilayer in solution encompasses contributions from the inter-head, head-tail, inter-tail, head-solvent, and the tail-solvent pair interaction energies. We would like to note that the pair, bond and angle energies presented in the paper are averaged over the total number of beads in the system. An increase in the inter-lipid spacing or area per lipid is accompanied by a growth of voids in the hydrophobic regions in the membrane. The response of the bilayer to this increasing number of voids in the hydrophobic region of the membrane is determined by the shape of the lipids. Our measurements of the interfacial tension of the bilayers composed of lipids B, C and D demonstrate lower sensitivity to the increased inter-lipid spacing or area per lipid. The lipid molecules are able to increase the splaying of their hydrocarbon tails with the area per lipid, thereby filling in the voids in the hydrophobic region with minimal decrease in the inter-head, head-tail, and inter-tail enthalpic interactions, as shown in Figure 2. We attribute this response to the smaller head groups of the lipids and the splayed hydrocarbon tails. We would like to note that longer hydrocarbon tail lengths might contribute to tighter packing of the lipid molecules, and thereby increase the inter-head interactions. Furthermore, the differences in the inter-tail energies are attributed to the architecture and the tail length of the lipid molecule, and increases with the packing of the hydrocarbon tails and their length. Lipid A has a cylindrical shape with straight hydrocarbon tails and larger-sized head groups which span the lateral area (with respect to the bilayer normal) of the hydrocarbon tails. Hence, the bilayer composed of lipid A has higher inter-head and head-tail interfacial energies, and smaller head-solvent and tail-solvent interfacial energies (see Figure 3) Our measurements of the head-solvent and tail-solvent interactions energies also highlight the differences in the molecular architectures of the different lipid models. The protruding head groups of lipids B, C and D are responsible for the higher head-solvent and tail-solvent interaction energies, in comparison to similar measurements for lipid A. However, both the head-solvent and tail-solvent interactions decrease with chain length due to the higher inter-head interactions resulting from the tighter packing between the lipid molecules. The bond and angle energies highlight the role of the molecule shape and hydrocarbon tail length, as shown in Figure 4. We observe the cylindrical shaped lipid architecture to have higher bond and angle energies. In addition, there is an increase in these energies with the hydrocarbon tail length.

The pair interaction, bond and angle energies between the different components of the system will determine the interfacial tension of the bilayer. We utilize the principle pressure components P_xx, P_yy, and P_xx to compute the interfacial tensionusing the following relation $\sigma {{r}^{2}}/{{k}_{b}}T=\left( {{L}_{z}}\times \left( {{P}_{xx}}-\frac{1}{2}{{P}_{yy}}-\frac{1}{2}{{P}_{zz}} \right) \right)$[14] where is the dimension of the cubic simulation box [40]. The average value of the interfacial tension was computed via the box averaging. In this techniq ue, five consecutive measurements of the interfacial tension are averaged. This process is repeated on the subsequent average values until there is one final value of tension, as opposed to computing one singular average. This technique is useful because it demonstrates a reduction in st and ard deviation (in comparison to experimental results) while maintaining fixed the total number of measurements of the interfacial tension. A similar approach was adopted by [58,59]. The errors bars for the interfacial tension measurements are computed using the st and ard deviation. From Figure 5, we can determine the areas per lipid corresponding to the tensionless bilayers for models A, B, C and D that are respectively, 1.205 r_c² (664 lipids), 1.303 r_c² (582 lipids), 1.328 r_c² (560 lipids) and 1.287 r_c² (584 lipids.) These results are also shown in Table 1. We surmise that the cylindrical shape of lipid A prevents the lipid tails from splaying with increase in the spacing between the lipid molecules, thereby resulting in a lower area per lipid corresponding to a tensionless membrane. Using the molecular shape argument, one would expect the area per lipid corresponding to a tensionless membrane to decrease with the hydrocarbon tail length. However, we note that the head- and tail-solvent interaction energies decrease with hydrocarbon chain length on account of the tighter packing of the lipids in the bilayers. Therefore, we hypothesize that the interplay between inter-tail, head- and tail-solvent interactions will also determine the area per lipid corresponding to a tensionless bilayer membrane.

We measured the thickness of the bilayer membranes associated with each lipid model. The average thickness of each tensionless lipid model was calculated by dividing the bilayer membrane into many patches, measuring the width of the patch at multiple sites (using the lipid head bead positions), and averaging over all the measurements. These measurements were done for a range of average area per lipid ranging from 1.17 to 1.34, for each lipid model. Figure 5 shows the bilayer thickness measurements for the different lipid architectures, and shows theaverage thickness for the tensionless membranes for models A, B, C and D to respectively be, 3.439 r_c, 4.181r_c , 5.018r_c and 6.080r_c. We expect that the combination of a large head group and the cylindrical shape of lipid model A to favor tight hydrocarbon chain packing in the hydrophobic region of the membrane, and not exhibit interdigitation of the lipid tails. Interdigitation is known to occur in bilayers composed of the inverted wedge-shaped lipids where the hydrocarbon tails of lipids in opposing monolayers interlock with each other leading to a lower bilayer thickness [2]. Our calculation of the membrane thickness does not demonstrate significant variation with the area per lipid. However, we do observe the membrane thickness to increase with the hydrocarbon tail length on a proportional scale. This proportional increase is quite surprising since it does not show any signs of interdigitation. One plausible explanation for this could be the fact that these simulations were run above the lipid’s transition temperature. In agreement with earlier studies, if the temperature of the system is above the lipid’s transition temperature the tails will not remain straight and are therefore unable to interdigitate [60].

We have also investigated the role of molecular architecture on the lateral and transverse diffusion of lipids in a bilayer membrane. Earlier theoretical studies [61,62,63,64,65,66] suggest that the lateral diffusion of lipids in a bilayer membrane depends upon the mean free volume as well as area per lipid, chain length, and lipid packing. Recent results from atomistic simulations and quasi-elastic neutron scattering experiments have shown that the theory of mean free volume is not sufficient for predicting the lateral diffusion of lipid molecules [61,62,63,64,67,68,69]. Experimental studies [3,62,70,71] have shown the area per lipid and hydrocarbon chain length to affect the lateral diffusion coefficient. Changes in the area per lipid and hydrocarbon chain length can influence the packing of the lipid molecules in the bilayer, and hence determine the lateral diffusion coefficient [3,62,70,71].

The transverse diffusion, or flip-flop, of lipids is the movement of a lipid molecule from one monolayer to another, across the hydrophobic region of the bilayer. Studies on the transverse diffusion have related the interaction energies of the lipids with the requirement of an activation energy barrier that has to be overcome for a flip-flop event to occur [72]. Furthermore, lipid molecules constituting bilayer membranes with lower inter-head interactions energy are more suitable to participate in transverse diffusion. We measure the lateral and transverse diffusion coefficient of the lipid molecules, for the different lipid models [73], and provide them in Table 1.

3.2. Interfacial Tension Measurements of Cholesterol-based Bilayers

En route to building a mixed phospholipid-cholesterol membrane, we investigated the interfacial properties of a pure cholesterol bilayer membrane. Where as pure cholesterol bilayers are not found in nature, we use the bilayer to determine the area per cholesterol molecule corresponding to a tensionless membrane. We repeat the process detailed earlier to build a pure cholesterol bilayer membrane with an average area per molecule varying from 1.00 to 1.15 r_c². Our measurements of the interfacial tension (as shown in Table 2) demonstrate the area per molecule corresponding to a tensionless cholesterol bilayer membrane to be 1.045 r_c². We note that the average area per molecule corresponding to a tensionless bilayer is smaller for cholesterol than for the different phospholipid models. We surmise that the relative higher rigidity of the cholesterol molecule, on account of its sterol ring, prevents conformational changes to fill in the voids created in the hydrophobic region of the bilayer while under positive tension. Hence, the inter-molecular spacing between the cholesterol molecules at which the membrane attains a tensionless state is smaller than for its phospholipid counterparts.

3.3. Mixed Bilayer: DPPC-Cholesterol

Human cell membranes have been found to have a 50 % concentration of cholesterol that plays an important role in regulating the fluidity in a membrane, as well as several other key functions [1]. We build a mixed bilayer membrane comprised of DPPC and cholesterol en route to designing a bio-inspired membrane. We use the area per molecule corresponding to tensionless bilayers for both molecules, and vary the relative concentration of cholesterol from 10 % to 50 %. Given the total lateral area of the simulation box we can determine the number of Cholesterol and DPPC molecules, for a given concentration of each of the species. Using an identical simulation box to the earlier investigations, we place the DPPC and cholesterol molecules in the pre-assembled bilayer and equilibrate the system for a duration of 20, 000t. Given that the tensionless area per molecule of cholesterol is smaller than that for the DPPC molecule, we observe the mixed bilayers to have negative tension, resulting in excess area of the membrane. To obtain a tensionless mixed membrane, we gradually stretch the membrane by 1% of its original area (while preserving the simulation box volume) until we get a tensionless mixed bilayer. Figure 6 shows images of a stable mixed 1:1 DPPC:Cholesterol membrane (a) before and (b) after stretching. Our approach has demonstrated the formation of a stable tensionless mixed lipid bilayer, and can be used for generating multi -component bilayers composed of other species of surfactant molecules.

4. Conclusion

In summary, we explore the role of molecular geometry of phospholipid molecules on the properties of lipid bilayer membranes as a function of inter-lipid spacing through the characterization of the interfacial tension, the interaction energies, diffusivity of the molecules and the membrane thickness. We have shown how both the molecular shape and the hydrocarbon tail length can determine properties of the lipid bilayer membranes. In addition, we have examined the properties of a pure cholesterol bilayer and found the molecular stiffness to be responsible for a relatively smaller average area per molecule. We have further demonstrated the design of stable mixed bilayers composed of DPPC and cholesterol using their average area per lipid corresponding to tensionless bilayer membranes. Our results can be potentially used to design bio-inspired membranes for probing the interactions between therapeutic peptides or delivery vehicles and the biological cell membranes, for applications in biomedicine, sensing and imaging.