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]. Phospholipid molecules have been shown to be the primary component of cell membranes [1]; these amphiphilic entities typically have a hydrophilic head group and two hydrocarbon tail chains. In addition, cholesterol is an essential component of all animal cell membranes, and is required for fluidity and structural integrity of the membranes. Cholesterol molecules are composed of a small hydrophilic head group, a rigid hydrophobic ring and a single hydrocarbon tail chain. The composition of the cell membrane has been found to influence its functions [2]. The self-organization of various amphipathic molecules enables the membrane to modulate its tension and mechanical properties which facilitate various physiological processes, or promote binding or catalytic events at the membrane interface [1,3,4]. In addition, the kinetics and thermodynamics of cell membranes are highly susceptible to changes in their external environment; including but not limited to pH changes, temperature fluctuations, and binding on active sites. [5] The composition and self-organization of different molecular species in the membrane [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] have been shown to be critical for the interfacial binding events on the cell surface [21,22]. A fundamental understanding of the underlying mechanisms driving the binding and insertion of various nanoparticles and macromolecules such as proteins, peptides, oligonucleotides or synthetic counterparts to the membranes will enhance the development of bionanomaterials via the use of interfacial binding for functional integration.

Experimental methods for studying the interactions of membranes with peptides and nanoparticles pose many challenges [23]. These challenges include limitations of experimental methods in capturing dynamics of peptide-membrane and nanoparticle-membrane interactions as these biological processes take place in a relatively short time scale. Most peptide-membrane interactions involve both electrostatic and hydrophobic contributions, yet it is still not clear which one plays a more dominant role in driving the initial binding and subsequent insertion of peptides during their interactions with model membranes. Computational simulation techniques can be used to model these membranes, along with their interactions with peptides and nanoparticles. All-atom simulations of lipid bilayers that resolve the dynamics of the various molecules are computationally expensive and limit the investigation to small spatio-temporal scales [24,25,26,27,28,29,30]. These tools are not suitable for addressing phenomena occurring on the mesoscale, such as membrane fusion and rupture, domain formation in the multicomponent membranes as well as the structural and dynamical effects of nanoparticle adsorption and self-organization onto the membrane [31]. Dynamics spanning large length and time scales can be resolved via coarse-grained [32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52] implicit solvent models which are used along with Monte Carlo [24,25,33,34,35,53,54], Molecular Dynamics [31], Brownian dynamics simulation methods [55,56] or mean field theoretical approaches [57,58,59,60,61,62,63,64,65].

The interactions between nanoparticles and biomembranes have been extensively studied previously by using various computational techniques. [37,66,67,68,69,70,71]. We have adopted two different Molecular Dynamics (MD)-based mesoscopic simulation techniques for the investigations presented in this paper. The first simulation technique is entitled Dissipative Particle Dynamics (DPD) [7,37,42,72,73] which simultaneously resolves both the molecular and continuum scales as well as the hydrodynamics of the system. The DPD method has been used to investigate the dynamics and morphology of self-assembly, phase separation and phase transition in lipid systems [44,46] block co-polymers [74,75], dense colloidal suspensions [76], polymers in dilute solution or in a melt [77] and chains in microfluidic channels [78,79]. The technique is previously used to study the translocation process of the nanoparticles across lipid vesicles [80] and membranes [81,82]. The method has been also used to investigate the different pathways of spontaneous nanoparticle penetration into a vesicle [83] and examine spontaneous insertion of amphiphilic nanotubes into a lipid vesicle and membrane, and their organization in the bilayer [36,38]. DPD is limited in its ability to capture the chemical specificity in peptide-membrane interactions because the highly coarse-grained approach makes it difficult to differentiate various types of amino acids in peptides. In addition, DPD model does not capture the electrostatic interactions between charged groups explicitly. As a comparison, the second simulation technique used for the investigations presented in this work is the coarse-grained (CG) MD simulation technique using the implicit solvent Dry Martini model [84]. Martini is a coarse-grained force field which has been used to investigate various biomolecular processes [52]. The Martini model [85,86] has been used to investigate the insertion and organization of nanoparticles (such as peptides) in lipid membranes. To our knowledge, there is no previous study on the membrane-nanoparticle interactions by using Dry Martini model. It is also worthy to note that since no solvent is involved in the system, the hydrodynamics of the peptide in aqueous environment is not captured using Dry Martini technique. Thus, we are not able to investigate the diffusion of peptide molecules in solvent prior to their binding to the membrane.

In this study, we are interested in investigating the principle mechanisms driving the capture, binding and insertion of peptides and peptide-mimetic nanoparticles into bio-inspired membranes. We will also examine the effect of explicit and implicit solvent approaches on the interaction of natural and synthetic nanoparticles with membranes.

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 [36,37,38,39,73]. 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 [73,87]. 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 random 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 ${{\mathbf{F}}_{i}}=\sum\limits_{j\ne i}{{{\mathbf{F}}_{c, ij}}+{{\mathbf{F}}_{d, ij}}+{{\mathbf{F}}_{r, ij}}}$. The soft spheres interact via a soft-repulsive force (${{\mathbf{F}}_{c,ij}}={{a}_{ij}}(1-{{{r}_{ij}}}/{{{r}_{c}}}\;){{\mathbf{\hat{r}}}_{ij}}$, for r_ij < r_c and ${{\mathbf{F}}_{c, ij}}=0$, for r_ij ≥ r_c), a dissipative force (${{\mathbf{F}}_{d, ij}}=-\gamma {{\omega }^{d}}({{r}_{ij}})({{\mathbf{\hat{r}}}_{ij}}\bullet {{\mathbf{v}}_{ij}}){{\mathbf{\hat{r}}}_{ij}}$) and a random force (${{\mathbf{F}}_{r, ij}}=-\sigma {{\omega }^{r}}({{r}_{ij}}){{\theta }_{ij}}{{\mathbf{\hat{r}}}_{ij}}$), where ${{\omega }^{d}}(r)={{[{{w}^{r}}(r)]}^{2}}={{(1-r)}^{2}}$ $\text{(for }r<1)$, ${{\omega }^{d}}(r)={{[{{w}^{r}}(r)]}^{2}}=0$ $\text{(for }r\ge 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|, ${{\mathbf{\hat{r}}}_{ij}}$= r_ij/r_ij, r = r_ij/r_c, γ is viscosity related parameter used in the simulations, σ is the noise amplitude, θ_ij(t) is a randomly fluctuating variable from Gaussian statistics, ω^d and ω^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 [73]. The constraints imposed on the random and dissipative forces by certain relations ensure that the statistical mechanics of the system conforms to the canonical ensemble [73,87].

The relation between the pair repulsion parameter a_ij and the Flory interaction parameter χ for a bead number density ρ = 3r_c⁻³ is given by $\chi =(0.286\pm 0.002)({{a}_{ij}}-{{a}_{ii}})$ [73].

We model a ternary membrane encompassing two distinct amphiphilic phospholipid species, 1, 2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC), 1, 2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) and a sterol species, cholesterol, as respectively shown in Figure 1(a), (b) and (c). Individual lipid molecules are represented by bead-spring models, and are modeled by a head group comprised of three hydrophilic beads and two hydrocarbon tails represented by three hydrophobic beads each. DPPC is a cylindrical-shaped molecule with a large head group and two hydrocarbon tails. DMPC is an inverted wedge-shaped molecule with a smaller head group and two-hydrocarbon tails [88]. The amphiphilic sterol molecule is modeled by a head group comprising of two hydrophilic beads and, a rigid ring and a single hydrocarbon tail encompassing seven hydrophobic beads. Two consecutive beads in a chain are connected via a bond that is described by the harmonic spring potential ${{E}_{bond}}={{K}_{bond}}{{((r-b)/{{r}_{c}})}^{2}}$, 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.5r_c, respectively [36,37,38,39,42]. The three-body stiffness potential along the lipid tails has the form ${{E}_{angle}}={{K}_{angle}}(1+\cos \theta )$ where θ 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 [39].

Figure 1. (a-d) Coarse-grain model (for DPD) of (a) DPPC, (b) DMPC, (c) cholesterol, (d) nanopin. (e-f) Dry martini coarse-grain model (for MD) of (e) DPPC, (f) peptide PMLKE with residue gly ( blue ), glu ( yellow ), lys ( magenta ), leu ( black ), met ( red ) and pro ( lime ).

DownLoad: Full-Size ImgPowerPoint

We draw a correspondence between our model and physical systems via the experimental properties of biological lipid bilayers. We obtain the characteristic length scale (r_c = 0.76 nm) for our model through the comparison of experimental measurements of the interfacial area per lipid of a DPPC bilayer at 50 °C with similar measurements from our simulations [90,91]. Similarly, the time scale τ is found to be 6.0 ns by comparing the experimental measurement of the diffusion coefficient of dipalmitoylphosphatidylcholine (DPPC) bilayer at 50 °C with that obtained from the simulations [90,91]. Using a temperature of 50 °C, the energy scale is calculated to be ε = k_BT = 4.5 × 10⁻²¹ J.

We model the amphiphilic nanoparticles after short peptides that have the capability to bind or penetrate the membrane of a cell. Most of these peptides are amphipathic and can be categorized into two different families: (1) antimicrobial peptide (AMP) and (2) cell-penetrating peptide (CPP). In this paper, we specifically model one type of short cell-penetrating peptide, called PMLKE [89]. We construct the model such that the helical moieties of the peptide are modeled by a rigid hydrophobic rod composed of 6 hydrophobic beads with a spacing of 0.5r_c between the center of mass of two consecutive beads, the hydrophobic rod is attached with a hydrophilic segment that serves to regulate the overall hydrophobicity of the nanopin. As PMLKE has 3 polar amino acid residues, we model the hydrophilic segment with 3 flexible beads and 1 bead that is joined to the rigid hydrophobic rod, as shown in Figure 1(d). Hence, the PMLKE is effectively modeled as an amphiphilic nanopin. We focus our investigations on short peptides to capture their spontaneous interaction with cell membranes and the mechanisms underlying their insertion into the membrane. We would like to note that the nanopin can also represent peptide-mimetic nanoparticles. Our results can be extrapolated to provide insight into experimental observations on the interactions of short peptides or synthetic counterparts with cell membranes [85,86].

The soft repulsive pair potential parameters for the lipid molecule head and tail beads were selected to capture its amphiphilic nature. The interaction parameters between the like components, ${{a}_{ij}}$, are based on the property of water [73]. The repulsion parameter between two beads of the same type is set at ${{a}_{ii}}=25$ (measured in units of ${{{k}_{B}}T}/{{{r}_{c}}}\;$) which is based upon the compressibility of water at room temperature [73] for a bead density of $\rho =3r_{c}^{-3}$. The soft repulsive interaction parameter a_ij between hydrophobic and hydrophilic beads is set at ${{a}_{ij}}=100 {{{k}_{b}}T}/{{{r}_{c}}}\;$, and is determined by using the Flory-Huggins interaction parameters, χ, as ${{a}_{ij}}={{a}_{ii}}+3.496\chi$ [73], for $\rho =3r_{c}^{-3}$.

The soft repulsive interaction parameters between the head (h), tail (t) beads of DPPC, DMPC, cholesterol, nanopins and the solvent (s) beads are assigned such thata_hh = 25, a_ht = 100, a_hs = 25, a_tt = 25, a_ts = 100 and a_ss = 25. The interaction parameter between cholesterol head beads and solvent beads is set to be a_hs = 15, due to the stronger affinity of the cholesterol hydroxyl head group with water.

In our simulations, the respective characteristic length and energy scale are r_c and k_BT. As a result, our characteristic time scale can be described as $\tau =\sqrt{m{{r}_{c}}^{2}/{{k}_{B}}T}$. Finally, σ = 3 and Δt = 0.02τ are used in the simulations along with the total bead number density of ρ = 3r_c^-3 and a dimensionless value of r_c = 1 [42]. The mass of all the beads is set to unity [37,42,44,72,73,90].

2.2. Molecular Dynamics with Dry Martini Coarse-Grained Model

We use a Martini coarse-graining method which groups approximately four non-H atoms as beads of similar radii [52]. Groups are categorized by their overall charge and polarity which leads to an established table of pairwise interaction parameters [52]. This form of coarse-graining has been used in other simulations involving membranes [52] and proteins [92]. To simulate longer length and time scales without increasing computational cost, an implicit solvent can be utilized with appropriate adjustments to the model. Dry Martini coarse-graining is one such model and has been used to simulate vesicles composed of ~22, 000 DOPC lipids [84] and KALP peptides in DPPC bilayers [84].

We model DPPC by coarse-graining the head groups with 4 beads and tail groups with 8 beads, as shown in Figure 1(e). Among the 4 head beads, the top two beads respectively possess a negative and a positive charge, the other two beads are neutral. We develop a CG model for the PMLKE peptide using the Dry Martini scheme. PMLKE is a cell-penetrating penta-peptide (CPP5) that has been studied for protein transduction [89]. It is a sequence of proline, methionine, leucine, lysine, and glutamic acid residues that makes it amphiphilic. In physiological pH conditions, lysine and glutamic acid side chains are positively and negatively charged, respectively, and the remaining side chains are only slightly polar. The coarse-grained structure is composed of a glycine-like nonpolar backbone with unique groups branching off the side. Figure 1(f) illustrates this coarse-grained representation.

3. Results and Discussion

3.1. Creating a Tensionless Membrane

We begin with a stable pre-assembled 1:1:2 DPPC, DMPC and cholesterol ternary component membrane encompassing of 15, 998 molecules. The molecular species are organized such that the hydrophobic tails point to the mid-plane of the bilayer to form the hydrophobic core with the hydrophilic head groups exposed to the solvent. The simulation is run in the canonical ensemble with three dimensional periodic boundary conditions. Following the equilibration phase, the simulation box dimensions are adjusted to change the area per lipid of the lipid molecules, without changing the volume of the simulation box [38] to obtain a tensionless membrane. The corresponding simulation box dimensions is determined to 98 × 98 × 41.6494 r_c³. These dimensions are found through stretching the x and y dimensions of the box while keeping the volume constant. We calculate the surface tension using the following equation $\gamma \text{=}{{L}_{\text{z}}}({{P}_{zz}}-\frac{1}{2}{{P}_{XX}}-\frac{1}{2}{{P}_{yy}})$ [93].

3.2. Nanopin Spontaneous Insertion using the Dissipative Particle Dynamics Simulation Technique

We introduce a total of 80 nanopins into the simulation box at random positions outside the interaction range from each other as well as the ternary membrane, as shown in Figure 2(a). Nanopins are introduced in increments of 10 in order to prevent significant aggregation between the nanopins prior to their insertion into the lipid bilayer. At early times, some of the nanopins are observed to interact with each other and form aggregates encompassing 2-6 individual nanopins. The aggregation of the nanopins is promoted by the unfavorable enthalpic interactions between the hydrophilic solvent and the hydrophobic segments of the nanopins. Following the aggregation, we observe spontaneous insertion of nanopin aggregates into the membrane bilayer, as shown in Figure 2(b-d). We would like to note that other nanopins that are not part of an aggregate are observed to interact directly with the membrane.

Figure 2. Images of the system with 30 nanopins and a three-component lipid membrane. Only portion of the membrane and nanopins are shown. The lipid head and tail beads are respectively in blue and red. All the images track a nanopin which is identified with a red arrow. (a) At time t = 0 τ, the nanopins are initially introduced outside the interaction range from each other and the membrane. (b) At time t = 2000τ, a number of nanopins aggregate with each other. (c) At time t = 5500τ, a nanopin (indicated by an arrow) is captured by the membrane. (d) At time t = 6500τ, the nanopin (indicated by an arrow) is inserted in the membrane and remains in the membrane thereafter.

DownLoad: Full-Size ImgPowerPoint

We investigate the dynamics of the nanopin insertion into the lipid bilayer by tracking the distance between the center of mass of a single nanopin and the membrane midplane in the vicinity of the nanopin, during an insertion event (as shown in Figure 3). The membrane region in the vicinity of the nanopin is defined by a cylindrical section of the membrane with a diameter of 4 r_c with approximately 80-100 amphiphilic molecules. The location of this region is dynamic and is updated according to the position of the specific nanopin. Our results show that initially the nanopins diffuse in the solution, outside the interaction range from the membrane. Once the nanopin is within a distance of 6 r_c from the midplane of the membrane, a rapid decrease of the nanopin-membrane mid-plane distance is observed as the nanopin inserts itself into the membrane within a time interval of 330τ. The steady state value of the distance of the nanopin to the membrane midplane is defined as the “insertion depth”. The penetration of the membrane by the nanopin is reflected by a small value of the insertion depth. Our results demonstrate the average insertion depth for 80 nanopins to be given by 1.5 ± 0.4 r_c (1.1 ± 0.3 nm).

Figure 3. Time evolution of (a) distance from center of mass of nanopin (identified in Figure 2) to the midplane of local membrane bilayer; (b) interaction count of hydrophobic and hydrophilic components of the nanopin with different components of membrane during an insertion event in a time interval of 10000τ.

DownLoad: Full-Size ImgPowerPoint

We observe the binding events to be initiated by interactions between the hydrophobic components of the nanopins and the bilayer. Our results demonstrate that the favorable interactions between the hydrophobic components of the nanopins and the lipid drive the spontaneous capture of the nanopins by the membrane bilayer and their subsequent transmembrane insertion. These observations are supported by the time evolution of the interaction count between hydrophilic and hydrophobic components of the phospholipids, cholesterol and nanopin, as shown in Figure 3(b). A pair of beads is considered to be interacting if their center-to-center distance is smaller than the interaction cut-off distance. The initial binding of a single nanopin to the membrane is characterized with a rapid increase in the number of interactions between the hydrophobic components of the bilayer and the nanopins in comparison to that for the hydrophilic components. The complete insertion of the nanopin is demonstrated by the steady value of the interaction count between the hydrophilic and hydrophobic components of the nanopin and the membrane. Our results indicate that the interactions between the hydrophobic components of the membrane and nanopin are responsible for the insertion of the nanopin into the membrane and its retention by the membrane. In addition, we find the favorable enthalpic interactions between the hydrophobic components of the nanopin and the membrane is responsible for enabling the capture of the nanopin by the membrane.

3.3. Peptide Insertion into Membrane Bilayer using the Molecular Dynamics Simulation Technique with Dry Martini Coarse-Grained Model

We introduce a single peptide molecule, PMLKE, into a simulation box encompassing a single-component (DPPC) tensionless lipid membrane that is derived by using the same protocol as detailed in the earlier section. The initial distance between the center of mass of the peptide molecule and hydrophobic core of the membrane is varied from 4 to 5 nm. The distance is selected to ensure that there is no initial interaction between the peptide and membrane. The simulation is run for a time interval of 9 × 10⁻¹¹ s. and is repeated 10 times by introducing the peptide at different initial locations above the membrane.

We observe that the peptide molecule spontaneously inserts itself into the membrane bilayer, as shown in Figure 4. At earlier stages of the simulation, we observe the peptide approaches to the bilayer until it is captured by the surface of the bilayer. During that process, the peptide adopts an orientation with its backbone parallel to the surface of bilayer in its vicinity. After initially binding to the bilayer, the peptide is further accommodated into the hydrophobic region of the membrane. We measure the distance between the center of mass of peptide and the mid-plane of the membrane (in its vicinity) as a function of time, as shown in Figure 5(a). Our results indicate the insertion dynamics of the peptide into the membrane is very similar to that for the nanopin using the coarse-grained explicit solvent DPD approach. The entire insertion process begins with the gradual approach of the peptide to the membrane followed by the initial interaction between the peptide and membrane. The capture of the peptide by the membrane is indicated by the sudden increase in the interaction count between the peptide and the membrane. The insertion of the peptide into the bilayer and its subsequent retention by the membrane is characterized by the steady state values of the interaction count. The average insertion depth is calculated to be 1.8 ± 0.3 nm, which is in agreement with corresponding measurements using the DPD approach.

Figure 4. Images of the system encompassing a single peptide PMLKE and one-component (DPPC) lipid membrane. Only portion of the membrane is shown. The lipid head and tail beads are respectively in green and orange. (a) At time t = 0 τ, the peptide is initially introduced outside the interaction range from the membrane. (b) At time t = 140 ps, the peptide has bound with membrane surface. (c) At time t = 200 ps, part of the peptide has inserted into the membrane. (d) At time t = 300 ps, the insertion process is complete and the peptide remains inserted in that orientation.

DownLoad: Full-Size ImgPowerPoint

Figure 5. Time evolution of (a) the distance of the center of mass of the peptide from the mid-plane of the local bilayer region and (b) the interaction count between the hydrophilic and hydrophobic components of the peptide and membrane. The inset shows time evolution of the interaction count. The Dry Martini coarse-grained model is applied for the system with a single peptide molecule and DPPC membrane.

DownLoad: Full-Size ImgPowerPoint

We measure the number of interactions between the different components of the peptide and the membrane in order to elucidate the driving force of the peptide insertion, as shown in Figure 5(b). The initial approach of the peptide to the membrane results from the occasional interactions between hydrophobic components of peptide and lipid head groups of the membrane, as shown in the inset of Figure 5(b). The subsequent capture of the peptide and its insertion into the membrane is primarily characterized by a rapid increase in the interaction count between the hydrophobic components of the peptide and the membrane. Due to proximity of the peptide hydrophobic groups to the lipid head groups, the corresponding interaction count also reflects the capture and insertion of the peptide. Our results demonstrate that the peptide insertion into the membrane is initiated, and sustained by the favorable enthalpic interactions between the hydrophobic groups of the peptide and membrane.

4. Conclusion

In summary, we have modeled the interactions between an amphiphilic cell-penetrating peptide PMLKE and its synthetic counterpart with a lipidmembrane, using the explicit and implicit solvent coarse-grained models. We demonstrate spontaneous binding and subsequent insertion of the peptide (and nanopin) into the membrane. Our results show that the insertion of nanopins and the peptide are driven primarily by the favorable interactions between their hydrophobic components. In addition, insertion depths of the nanopins and the peptide using respectively, explicit and implicit solvent models are found to be in a good agreement with each other. These results can be used to design materials that functionally integrate cell membranes with natural or synthetic nanoparticles for applications in medicine, sensing and energy.

Acknowledgments

The authors would like to thank Jia Li for assistance in providing data for running portion of the simulations. We would also like to acknowledge the use of high performance computational resources at the Rutgers Discovery Informatics Institute (http://rdi2.rutgers.edu/).

Conflict of Interest

All authors declare no conflicts of interest in this paper.