Translocation thermodynamics of linear and cyclic nonaarginine into model DPPC bilayer via coarse-grained molecular dynamics simulation: implications of pore formation and nonadditivity.

Structural mechanisms and underlying thermodynamic determinants of efficient internalization of charged cationic peptides (cell-penetrating peptides, CPPs) such as TAT, polyarginine, and their variants, into cells, cellular constructs, and model membrane/lipid bilayers (large and giant unilamellar or multilamelar vesicles) continue to garner significant attention. Two widely held views on the translocation mechanism center on endocytotic and nonendocytotic (diffusive) processes. Espousing the view of a purely diffusive internalization process (supported by recent experimental evidence, [Säälik, P.; et al. J. Controlled Release 2011, 153, 117-125]), we consider the underlying free energetics of the translocation of a nonaarginine peptide (Arg9) into a model DPPC bilayer. In the case of the Arg9 cationic peptide, recent experiments indicate a higher internalization efficiency of the cyclic structure (cyclic Arg9) relative to the linear conformer. Furthermore, recent all-atom resolution molecular dynamics simulations of cyclic Arg9 [Huang, K.; et al. Biophys. J., 2013, 104, 412-420] suggested a critical stabilizing role of water- and lipid-constituted pores that form within the bilayer as the charged Arg9 translocates deep into the bilayer center. Herein, we use umbrella sampling molecular dynamics simulations with coarse-grained Martini lipids, polarizable coarse-grained water, and peptide to explore the dependence of translocation free energetics on peptide structure and conformation via calculation of potentials of mean force along preselected reaction paths allowing and preventing membrane deformations that lead to pore formation. Within the context of the coarse-grained force fields we employ, we observe significant barriers for Arg9 translocation from bulk aqueous solution to bilayer center. Moreover, we do not find free-energy minima in the headgroup-water interfacial region, as observed in simulations using all-atom force fields. The pore-forming paths systematically predict lower free-energy barriers (ca. 90 kJ/mol lower) than the non pore-forming paths, again consistent with all-atom force field simulations. The current force field suggests no preference for the more compact or covalently cyclic structures upon entering the bilayer. Decomposition of the PMF into the system's components indicates that the dominant stabilizing contribution along the pore-forming path originates from the membrane as both layers of it deformed due to the formation of pore. Furthermore, our analysis revealed that although there is significant entropic stabilization arising from the enhanced configurational entropy exposing more states as the peptide moves through the bilayer, the enthalpic loss (as predicted by the interactions of this coarse-grained model) far outweighs any former stabilization, thus leading to significant barrier to translocation. Finally, we observe reduction in the translocation free-energy barrier for a second Arg9 entering the bilayer in the presence of an initial peptide restrained at the center, again, in qualitative agreement with all-atom force fields.

Introduction

Efficient, targeted delivery of molecular cargo (therapeutics, chemical sensors, etc.) to specific cells remains an elusive challenge.[1−5] Cell-penetrating peptides (CPPs) continue to offer an intriguing opportunity to effect cell-specific internalization of molecular cargo with minimal cytotoxicity.[1,4−6] The expansive scientific research literature in this area attests to the significance of this broad class of chemical transporters to a wide spectrum of applications ranging from cancer therapy to (bio)chemical sensors.[7,8] The human genome encodes inherent CPPs as individual elements or components of larger macromolecular entities (TAT segment).[3,9] However, numerous synthetic peptides have been proposed and studied in the literature.[5,10] Moreover, though, we are concerned here with unstructured (with respect to canonical protein secondary structure elements) cationic peptides, the class of amphipathic, helical species also presents an alternative class of internalizing vectors.[11,12] It is natural to ask if novel, more efficient CPPs can be designed de novo to exploit the power of targeted cellular delivery. One approach to achieve design success entails understanding of fundamental mechanisms and thermodynamics of the process associated with existing peptide systems. Specifically, what molecular/chemical thermodynamic and structural elements facilitate the internalization process.

Currently, two predominant mechanisms of CPP internalization are endocytosis (ATP-derived energy dependent; clathrin-mediated vesicularization) and diffusion-like process mediated by broadly defined structural perturbations of membrane/bilayer (i.e., membrane curvature induction arginine-rich CPPs, local water, and possibly other cellular membrane components; the latter mechanism is considered to be energy-independent in the CPP literature).[13] We focus in this work on issues related to CPP translocation via a diffusion mechanism, the plausibility of which was recently demonstrated in giant plasma membrane vesicles (GPMVs) with no cellular machinery for endosome formation[14,15] as well as matrix-assisted laser desorption-ionization time-of-flight mass spectrometry (MALDI-TOF MS).[15,16] Experiments using live-cell microscopy and analytical ultracentrifugation also demonstrate that the uptake efficiency, interpreted by fluorescence intensities over time within cells or vesicles after introduction of CPPs into extracellular regions, indicates that cyclic Arg9 and decaarginine translocate faster than their linear counterparts.[17] This is further taken as a demonstration of the higher uptake efficiency of the cyclic form. The cyclic form, we clarify, is the polyarginine with covalent bonds between the peptides ends. Although this effect is unambiguous, there is little understanding about the origins, causes, and molecular/chemical/thermodynamic determinants of these uptake efficiency differences. Recent advanced experimental methods such as X-ray with neutron reflection,[18] solid-state nuclear magnetic resonance spectroscopy (NMR),[19,20] optical sectioning and state-of-the-art single-molecule microscopy,[21] conductance measurements,[22] and so on explore the important interaction between cationic peptide and lipid phosphates that distorts the membrane structure and initiates pore formation. Such experiments also demonstrate that the formation of pores inside the membrane is crucial for translocation of cationic peptides. We note that the peptide length also influences the ostensible uptake efficiency,[6] although we do not address this point in the current work.

A further aspect, emerging from molecular dynamics modeling of representative model systems, revolves around the nature of the water perturbations associated with the translocation process. Numerous simulation studies have been employed to understand the CPP internalization process from a molecular/atomic perspective.[23−30] Recent ideas invoke the formation of membrane spanning water pores (i.e., aqueous conduits across the patch of membrane under scrutiny in the simulation) as perhaps necessary structural elements for CPP translocation via diffusion-like processes.[14,22,27] Using the all-atom GROMOS87 force field for peptide, Berger force field for lipid, and SPC model for water, Huang et al. showed a 80 kJ/mol reduction in the free-energy barrier for translocation (from bulk water to bilayer center).[29] The authors conclude that a reaction coordinate (rxn. coord.) capable of accommodating a water pore spanning the membrane bilayer would be the lower free energy path relative to one in which only modest water defects (not spanning the entire membrane thickness) are possible. The necessity of a water pore is predicated on structural perturbations of the bilayer itself, and thus that study recapitulates a series of molecular dynamics studies highlighting the intimate connection between translocation of charged peptides in bilayers and some type of structural perturbation on the scales of single (or several) lipid molecules.[30−33]

Here we consider several aspects related to Arg9 translocation through a model DPPC bilayer using umbrella sampling (US) molecular dynamics simulations coupled to pairwise additive coarse-grained (CG) force fields for all component-component interactions. We first consider the models, protocols, and related issues in the Methods Section. In Results and Discussion Section, we discuss results of our potentials of mean force (PMF) for linear and cyclic Arg9 translocation using center-of-mass (COM) distance reaction coordinate along pore and pore-free paths. The effect of different conformations of peptide on PMF, decomposition of PMF into different components, and enthalpic and entropic contributions to the total PMF have also been discussed in that section. Further analysis considers the PMF for the translocation of a second Arg9 in which a single Arg9 is restrained at the center of membrane (Nonadditivity Section). The important findings and the conclusions of our study are recapitulated in the Summary.

Methods

Simulation Protocol

General Molecular Dynamics Protocol

We carried out US[34] molecular dynamics simulations to study the translocation of a positively charged peptide, Arg9, across a model membrane, 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC). We considered a fully hydrated liquid crystalline lamellar phase (Lα) of a DPPC bilayer patch of 256 lipid molecules (128 per leaflet). We included a total of 7553 water molecules for the hydration of lipid molecules. Furthermore, we added 150 mM NaCl salt to the system to mimic cellular electrolyte concentration;[35] we acknowledge that the use of low salt concentrations in all-atom simulations has been criticized for incurring sampling insufficiencies as well as creating possibly long-lived potential gradients in the simulation cell.[25,31,36] Additionally, nine Cl– ions were added for each Arg9 to the system explicitly to maintain overall charge neutrality.

We used the Martini CG model as developed by Marrink et al.[37−39] to simulate interactions between system components. The latest release version of the Martini polarizable force field for arginine and water (version 2.2P) and DPPC and ions (version 2.0) was used in combination with the polarizable water model for all simulations. Use of the current polarizable water model in conjunction with the latest model for treating charged amino acid residues in the Martini formalism is necessary because we are effectively concerned with partitioning of polar and charged groups from a high dielectric solvent into a low-dielectric bilayer medium.[37] The polarizable water and charged residue Martini force fields have been shown to, in conjunction with one another, reproduce the Wimley–White transfer free energy for charged arginine, as shown in figure 7 of ref (40) and figure 1 of ref (38). At the current time, this force-field combination appears to be the most systematically developed and validated CG force field for treating this class of molecules. We stress that the agreement of absolute energy values between CG and all-atom force field based calculations should not be expected. Here we are concerned with relative free energies of peptides in solution and in bilayer environments. For this purpose, we believe that the current CG force field (MARTINI) provides sufficiently robust and reliable free-energy differences so as to allow investigation of the influence of “pores” consisting of lipids and water molecules on the peptide translocation process. The reliability of the MARTINI force field with respect to free-energy calculations has been demonstrated in recent literature.[37,38,40] The simulation cell consists of a rectangular box of dimensions 9.0 × 9.0 × 15.0 nm, yielding about 4 nm thick slab of lipid molecules surrounded by bulk water and ions. The components of the system are depicted in Figure 1: panel a is the Martini polarizable water model, panel b is the Martini CG DPPC lipid model, panel c is the cyclic Arg9 with covalently bound terminal backbone (BB) beads, and panel d is what we refer to as the linear chain spanning the cases of no end-to-end backbone distance harmonic restraint and intermediate lengths of 0.5, 1.0, and 1.75 nm. Note that both N- and C-termini of the linear peptide are neutral in the present study; this choice removes complications arising from interactions of these charged groups with the bilayers and allows us to focus on the nature of the charged side chains. After components of system outlined in Figure 1 were initially placed in the simulation cell, the system was minimized using the steepest descent method and then equilibrated via constant particle, pressure, and temperature (NPT) ensemble molecular dynamics simulations for 1 μs at 1 atm and 323 K (above the liquid-to-gel phase transition for DPPC under experimental conditions). During the MD equilibration, the area per lipid equilibrated to a value of 63.2 Å2, in agreement with published results[37] for this force field. (See SI figures S1–S3 for details.)

Figure 1

Bead representation of CG (a) polarizable water, (b) DPPC lipid, (c) cyclic Arg9, and (d) linear Arg9 molecules. All of the water beads are red. W (type POL), WP (type D), and WM (type D) are neutral, positively, and negatively charged beads of polarizable CG water. The tails (type C1) and head groups of CG lipid molecules are yellow and green. The NC3 (type Q0), PO4 (type Qa), GL1 (type Na), and GL2 (type Na) on lipid molecule are marked for choline, phosphate, and two carbonyl beads. The backbone and nonbackbone beads of CG Arg9 peptide are blue and purple. The marked string such as BB (type P5), SC1 (type N0), SC2 (type Qd), and SCP (type D) are for backbone, two noncharged, and positively charged beads, respectively. The types of beads are adopted from Martini force field version 2.2P for water and Arg9 and version 2.0 for lipid.

All MD simulations were carried out using MPI-supported GROMACS software package (version 4.6), patched with PLUMED, version 2.0. We followed a standard scheme for simulating our system, as provided by the Marrink group.[37] The simulations were carried out at a time step of 20 fs, and we updated the neighbor list every step. In all simulations, we used 3-D periodic boundary conditions and the minimum image convention to calculate nonbonded interactions. Both the nonbonded interactions such as Lennard-Jones (LJ) and electrostatics (Coulomb) were calculated by using simple spherical cutoff at a distance of 1.2 nm with a smooth switching function of distances 0.9 and 0.0 nm, respectively. The relative dielectric constant was set to 2.5 for use with the polarizable water force field. To maintain the temperature 323 K and a pressure of 1 atm for the systems, we used the Berendsen weak coupling scheme with time constants of τ= 1.0 ps and τ = 5.0 ps, respectively.[41] We employed two temperature coupling groups: water and ions were considered as one and DPPC and Arg9 were set as the second group. To keep the bilayer in a tensionless state, we used periodic boundary conditions with a semi-isotropic pressure coupling algorithm with a 3.0 × 10–4 bar–1 compressibility. The LINCS algorithm[42] was used to apply the bond constraint present in the Martini force field.

Umbrella Sampling

US methods require biased sampling of a chosen reaction coordinate within appropriately narrow bounds along the domain of possible values of the reaction coordinate. Presently, we choose the distance between the COM of the nine charged beads (bead name SCP in Martini nomenclature) of the peptide (SCPCOM) and the COM of the bilayer as the reaction coordinate (SCPCOM–COM). Furthermore, the sampling of the multiple regions along the reaction coordinate requires initial configurations in each region (in each window corresponding to a particular value of the reaction coordinate). Before discussing the method for generating initial configurations, we discuss the various systems we wish to consider. We performed five sets of US simulations defined by how we treated the end-to-end distance of the peptide backbone. We refer to the cyclic peptide as one in which the terminal backbone beads are covalently bound as in the experimental work of Lättig-Tünnemann.[17] We also study four cases for linear Arg9, which is not covalently bound at the termini. One case is a fully unrestrained linear peptide (no restraint), and the other three cases have varying end-to-end distances of 0.5, 1.0, and 1.75 nm. (For details, see SI Figure S4.) In summary, the simulated systems are: (A) linear, (B) cyclic, and linear configurations with three restrained end-to-end distances at a value of (C) 0.5 nm, (D) 1.0 nm, and (E) 1.75 nm in this study. In addition, we address the free energy of translocation of a second cyclic Arg9 across the membrane in the presence of a cyclic Arg9 placed at the center of bilayer patch for (F) cyclic Arg9. (See Table 1 for details.)

Table 1

System Setup for All the Umbrella Sampling Systems

US	peptide	DPPC	water	Na+	Cl–	time (ns)a	details
A	1 linear Arg9	256	7553	82	91	500	pore (0.0 to 0.5 nm)
B	1 cyclic Arg9	256	7553	82	91	500	pore (0.0 to 0.5 nm)
C	1 linear Arg9	256	7553	82	91	500	0.5 nm,b pore (0.0 to 0.5 nm)
D	1 linear Arg9	256	7553	82	91	500	1.0 nm,b pore (0.0 to 0.5 nm)
E	1 linear Arg9	256	7553	82	91	500	1.75 nm,b pore (0.0 to 0.5 nm)
F	2 cyclic Arg9	256	7544	82	100	500	pore (all windows)
G	1 cyclic Arg9	256	7553	82	91	500	no pore (except 0.0 nm)
H	1 linear Arg9	256	7553	82	91	500	no pore (all windows)

Plain simulation time for each US window.

Restrained backbone end-to-end distance.

We first considered generating initial configurations in the windows along the reaction coordinate SCPCOM–COM by growing in a Arg9 at the center of bilayer patch of the previously equilibrated systems and further equilibrating the peptide–bilayer–water–ion system for 100–200 ns after the growing in bilayer center phase. In a separate simulation, we grew in a cyclic Arg9 (designated as G) and a linear Arg9 (designated as H) in the bulk water phase and followed the same equilibration protocol. To prevent unnecessary drift of membrane, we applied a position restraint along the z dimension with a force constant of 1000 kJ/mol/nm2 on the charged groups of lipid molecule (NC3, PO4) during the growing-in and equilibration phase for simulations A–H. The growing of peptide inside the system (either in bulk water phase or at the bilayer center) was done in two steps. We first slowly raised the Lennard-Jones interactions up to normal strength over the course of a 10 ns simulation period using the method of thermodynamic integration as implemented in GROMACS, where step length (dλ) is set to 2 × 10–6 per time step, and soft-core potential was used to prevent bead overlap. In the following step, we slowly grew in the Coulomb interactions using the same protocol. To calculate the equilibrium PMF for Arg9, we constructed a total of 61 US windows (ranging from 0.0 to 6.0 nm) at a spacing of 0.1 nm along a reaction coordinate.

So far, we have only addressed how we generate initial configurations (starting configurations) for the windows where the peptide resides in bilayer center and in bulk water phase. We next discuss the protocols for generating the remaining windows. We considered generating the remaining windows using two protocols. First, starting with the peptide at the center window, we performed MD simulations harmonically biasing the COM of the nine charged beads of Arg9 to remain at bilayer center; from this simulation, we selected a configuration in which the peptide position had fluctuated to the next adjacent window moving outward to the bulk water phase. We repeated this protocol for the remaining windows, thus generating initial configurations for each window moving from bilayer center to bulk water; this we refer to as Protocol 1. The same protocol was followed in simulations G and H, except that we started with the peptide in bulk water and generated initial configurations for windows moving inward to the bilayer center; this is referred to as Protocol 2. Once the initial configuration of each window was generated, a harmonic potential with a force constant of 3000 kJ/mol/nm2 was applied to restrain the peptide along the reaction coordinate. Each window was sampled for 500 ns. Note that we considered the first 50 ns data as an equilibration period for each window and excluded those data during the computation of final PMF of our system. A detailed description of nine US systems is presented in Table 1. The additional harmonic restraint between first and last backbone bead with a force constant of 500 kJ/mol/nm2 for C–E simulations was applied by PLUMED2.0 package.

Simulations A–F form stable pore (water+lipids) inside the bilayer, which was validated by careful experimentation and visualization. The pore is considered to be a generic term referring to locally complexed lipids and water as the Arg9 enters the bilayer; the combined presence of lipid and water, as will be shown later, assists in stabilizing the highly charged peptide relative to the situation where no pore is attainable. However, extensive testing to generate stable pores inside the bilayer by trial and error demonstrated a viable region of the reaction coordinate where such configurations could persist indefinitely. Moreover, it is important to note that the translocation of peptide from center to bulk water was considered for simulations A–F, whereas the reverse process was considered for simulations G and H.

We pause here to note that the formation (or dissolution) of a pore in the bilayer is an example of a slowly evolving, orthogonal degree of freedom (orthogonal to the chosen reaction coordinate). We thus clarify that the calculations we performed and discuss here are free-energy profiles along local paths that do not completely sample the slowly evolving, orthogonal, pore-forming degree of freedom (perhaps as well as others that are unknown). Thus, the nature of the paths we have chosen are well-described by Figure 1 in the work of Huang et al.[29] On the basis of that work, we essentially are addressing two distinct local paths in this work. First, a path with a water/lipid pore formed when the peptide nears the bilayer center (path B in figure 1 of Huang et al.[29]) and second, a local path without a water/lipid pore when the lipid nears the bilayer center (i.e., path A or C in figure 1 of Huang et al.[29]). The spirit of our study is the same as that of Huang et al.,[29] although we are attempting to use CG force fields that allow us to probe a few other characteristics of the peptide system such as multiple end-to-end restraints, decomposition of the PMF, and the role of the solvation shell. We caution that the PMF along the pore-forming path we have chosen may not be the global, lowest free energy path; however, our aim is not to calculate that path but rather to simply comment on the effect of the presence of a pore. Furthermore, we acknowledge the nonuniqueness of our chosen reaction coordinate and the undefined (unknown) influences of slowly evolving degrees of freedom along this coordinate.

Postprocessing of Umbrella Sampling Simulations

The weighted histogram analysis method (WHAM) was used for postsimulation unbiasing of US data.[43] The free-energy parameters, f, for each umbrella window can be estimated by the following equation in an iterative manner[44]where N is the total number of umbrella windows, the variable ξ is the reaction coordinate, and the values n, W(ξ), and P( (ξ) are the number of sampling, biased potential, and biased probability distribution functions of the ith window. Once the f values are estimated, the unbiased probability of each configuration of the US simulation iswhere P0(R) is the unbiased probability of configuration R obtained from umbrella simulations. Finally, the unbiased probabilities can be used to project along any reaction coordinate (η) by using

We use the WHAM utility of Grossfield[45] to generate the final PMF. We ascertained that sufficient overlap of reaction coordinate values between adjacent windows was maintained. Using a modified in-house WHAM code implementing eq 3, we projected our original reaction coordinate, SCPCOM–COM to a new reaction coordinate, the distance between COM of peptide and bilayer (COM–COM). All results will be presented in terms of this alternate reaction coordinate (COM–COM). It allows us to compare the PMF contributions on the whole peptide with the PMF projected on COM–COM because it is complicated and not meaningful to calculate contribution between noncharged beads and charged beads of the same peptide on SCPCOM–COM rxn. coord.

To quantify free-energy contributions of different system components to the PMF, we use the relation[36]where Fα(η) is the instantaneous force of the component α acting along the chosen reaction coordinate, η0 is the lower limit, and η1 is the upper limit of this reaction coordinate. Because our chosen reaction coordinate is the z component distance between COM of peptide and bilayer, the net force along the reaction coordinate should be the difference between the instantaneous force acting between component α and the protein and the instantaneous force acting between component α and lipid bilayer. Because our reaction coordinate is just the relative distance between COM of peptide and bilayer, we can estimate the PMF along that reaction coordinate from the relative force between them. Therefore, PMF of component α can be computed directly from the relative instantaneous force acting between peptide and bilayer for that component, that is

The instantaneous relative force between the protein and bilayer for component α was computed by rerunning the trajectory of each window using Gromacs “mdrun”.

The final PMF and its standard error was estimated by using the block averaging method obtained from each consecutive 50 ns time block in the production run of each US window.[29] We ensured that the block size was significantly larger than the correlation time in each umbrella window. Details of the convergence tests have been presented in SI Figures S5–S12. The visual molecular dynamics (VMD) package[46] was used to monitor the simulation, visualization, and graphics preparation for this work.

Results and Discussion

Figure 2 shows snapshots of configurations from simulations A, B, G, and H at three different windows corresponding to bulk water, the lipid–water interface, and the bilayer center. There is a significant difference in the bilayer center region. The formation of a pore consisting of water and lipids is observed when the peptide resides in the interior of the bilayer for simulations A and B. However, no pore materializes in the interior region of bilayer for simulations G and H. Moreover, we have confirmed via analyses of simulation trajectories that once the pore is formed, both the upper and lower leaflets contribute individual lipid molecules that participate in forming the local pore around the translocating peptide. We thus consider the former two simulations as transferring Arg9 along a pore-forming path, whereas we consider the latter two as along a pore-free path. We adopt such language from the work of Huang et al.,[29] in which the authors demonstrate that the PMF of cyclic Arg9 translocation through a PC bilayer using all-atom molecular dynamics US simulations is reduced when the reaction coordinate of focus is able to accommodate a pore-formation pathway for translocation.

Figure 2

Representative snapshots of simulation (A) transfer of linear Arg9 along pore-forming path, (B) transfer of cyclic Arg9 along pore-forming path, (G) transfer of cyclic Arg9 along pore-free path, and (H) transfer of linear Arg9 along pore-free path. For each simulation, peptides are restrained at (a) bulk water (6.0 nm), (b) interface (2.0 nm), and (c) center (reaction coordinate 0.1 nm). The phosphate groups of lipid are green, water is red, and the peptide is blue. For clarity, other beads of lipid and ions are not shown in the Figure.

Figure 3a shows the PMFs for cyclic and linear (fully free peptide, no end-to-end restraint) Arg9. The largest values of the rxn. coord. distance correspond to the peptide in the bulk solution; peptide at bilayer center corresponds to a value of zero. Two characteristic sets of curves are evident. First, the two curves with a flattened region near the center of bilayer correspond to the linear and cyclic Arg9, with US windows generated using Protocol 1, and the other two curves without flattened region correspond to the linear and cyclic Arg9, with US windows generated using Protocol 2. For both Protocols, the current force field and methodological combination predict substantial barriers ranging from 240 to 330 kJ/mol. Clearly, all PMFs show that the free energy of transfer of Arg9 monotonically increases from bulk water to the interior of bilayer. The presence of a large barrier in all PMFs suggests that the transfer process occurs rarely, and in fact, our results would suggest that translocation of highly charged and hydrophilic Arg9 (cyclic or linear) peptide through a PC bilayer is not purely diffusion-based but may involve anionic lipid composition or some degree of cellular internalization machinery in real systems.[21,47,48]

Figure 3

(a) PMFs with standard errors for the transferring of linear Arg9 and cyclic Arg9 along the pore and pore-free paths. (b) Corresponding curves for the transferring of linear Arg9 spanning the cases of three end-to-end harmonic restraint of lengths 0.5, 1.0, and 1.75 nm. The PMFs at the pore region are highlighted in the inset. For the sake of clarity, a vertical offset of 25 kJ/mol is added for all PMFs. No shift has been made for the inset Figures.

The inset shows the region of the reaction coordinate between 0 and 0.6 nm. In this region, we observe striking differences in the PMFs predicted using the two protocols. Protocol 2 selects a pore-free translocation pathway (although a pore is formed only in the window where the peptide resides at bilayer center for cyclic Arg9). Protocol 1 selects a pathway that includes a pore (previously defined). The consequences of the pore-forming pathway are clear; there is a nontrivial stabilization of ∼90 kJ/mol once a pore is formed within the bilayer. This behavior qualitatively agrees with the conclusions of Huang et al. A reaction coordinate capable of allowing a pore stabilizes the highly charged peptide at bilayer center. We will turn to contributions to the PMF further later. Second, we observe that the differences between the cyclic and linear peptides as they translocate the bilayer, via either pathway, are rather small in relation to the total barrier heights; along the pore-free path, the PMFs of linear and cyclic Arg9 from bulk water to the center of bilayer are about 330 and 329 kJ/mol; along the pore-forming path, the PMFs of linear and cyclic Arg9 from bulk water to the center of bilayer are about 237 and 242 kJ/mol, respectively. This appears to be in contradiction with recent results, indicating the higher uptake efficiency of cyclic polyarginine in C2C12 mouse myoblasts cell systems.[17] The difference may arise from the different physical systems studied in the two cases as well as lack of the atomistic details of the coarse-grained force field we use. However, considering that the desolvation of nine positive charges is a fairly unfavorable process, it is self-consistent that differences in desolvating a linear or cyclic Arg9 would be small in the context of the overall desolvation free-energy penalties. Furthermore, the barrier heights we observe are large, again consistent with the prediction of the all-atom calculations of Huang et al.[29]

Moreover, the PMFs for all cases show a plateau at the bulk water region, that is, around the 5 to 6 nm region. We did not observe any local minima in PMF profiles at the lipid–water interface. The absence of such local minima in PMF at the lipid–water interface region for both cyclic and linear Arg9 indicates that the positively charged peptide does not prefer to bind with the head groups of neutral PC lipid molecule. However, such lack of binding of Arg9 peptide with the headgroup of DPPC is in agreement with experimental observation.[49] We also observe in free simulation of Arg9 in a solution bathing a bilayer patch that very little binding of the peptide to the bilayer surface occurs. (See Figure S13 in the SI.) This has been previously demonstrated for this force field when compared with the use of the BMW water model.[50] It appears that the Martini CG model used here agrees qualitatively with the suggested lack of binding to pure PC bilayers from experiments.

To explore the effect of peptide compactness on such free-energy profile, we further calculated the PMF for linear Arg9 restrained at three different end-to-end distances (simulations C–E). The corresponding PMFs are presented in Figure 3b. The overall features of the PMF profiles are quite consistent with the above discussion. However, it is important to note that the PMF profile for cyclic Arg9 is well-matched with the PMF profile of linear Arg9 restrained at the smallest end-to-end distance.

PMF Decomposition

We consider specific contributions from system components (i.e., water, lipid, ion, and peptide) to the total PMF; we decompose contributions based on the approach of Allen et al.[36] Specifically, we decomposed the total PMF for linear and cyclic Arg9 translocating from the bilayer center to bulk solution along the pore-forming path. The results of the decomposition are presented in Figure 4. To validate the decomposition procedure (the individual component contributions to the PMF were obtained by integration of the average force from that component along the reaction coordinate), we added the contributions from the four different components of the system. The sum of the component contributions to the total PMF obtained from force integration matches the calculated PMF obtained from WHAM analysis. (See Figure S14 in the SI.) For comparison, we also decomposed the total PMF profiles obtained from simulations of pore-free path and included them in the same Figure. The corresponding average force profiles acting along the reaction coordinate for each component are presented in SI Figure S15.

Figure 4

Components of PMF arising from water, ion, membrane, and self-contribution of peptide are shown. The inset figure shows the contribution of sodium and chloride ions. A vertical offset of 50 kJ/mol is added for clarity.

From Figure 4, water and ion have large destabilizing contributions, whereas lipid confers stabilization. However, such stabilizing contribution arising from lipid is not sufficiently adequate to balance the large destabilizing contributions arising from water and ions. We observe that the slope of PMF for cyclic Arg9 is nearly opposite for pore and pore-free paths at the central region of bilayer (i.e., 0.0 to 0.5 nm). However, an identical PMF profile for each component is observed in the remaining corresponding regions. It is apparent from our calculation that the Arg9 experiences major penalties from the water as it approaches the bilayer center. The slope of each profile is indicative of the sign of z-component force acting on Arg9 at different distances from the center of bilayer. The contribution of water to the total PMF is continuously destabilizing, moving from bulk into the bilayer interface and reaching a plateau at the headgroup region, and becomes nearly flattened toward the center of bilayer. We found the barrier heights for linear and cyclic Arg9 as obtained along pore-forming paths are ∼530 kJ/mol. However, the corresponding barrier height for linear and cyclic Arg9 as obtained from pore free path are ∼510 kJ/mol. Therefore, the difference in free energy for this component between pore and pore-free path is ∼20 kJ/mol. This large free-energy penalty is related to the dehydration of highly solvated Arg9 inside the bilayer. To illustrate this, we compute the average number of water molecules (central bead, W, of water was considered for the calculation) present within a distance 0.67 nm of all peptide beads as a function of reaction coordinate, shown in Figure 5a. Such width of shell has been chosen by the calculation of radial distribution function (rdf), which has been calculated between the peptide and water beads (see Figure S16 in the Supporting Information), the first minima of which was found at 0.67 nm.

Figure 5

Average number of species present within 0.67 nm of all peptide beads. (a) Average number of water, (b) average number of chloride ions (Cl–) and sodium ions (Na+), and (c) average number of phosphate (PO4) and choline (NC3) beads. Vertical offsets of 20, 0.5, and 5 are added to panels a–c for clarity. Symbol coding scheme adopted in the Figure is the same as in Figure 4.

The average number of water molecules present in the first hydration shell of Arg9 falls significantly from bulk water to the interior region of bilayer in all the cases. However, a substantial difference in average water number of a value of ∼10 has been found between the pore and pore-free paths for both Arg9 at the central region of bilayer. Such difference in average water number is in agreement with the larger destabilization contribution arising from water for the pore-forming path. Contributions of membrane to the total PMF are increasingly stabilizing as the Arg9 approaches the bilayer region from bulk water. The change in values of PMF for this component from bulk water to the bilayer center as obtained along pore-forming paths are about −465 kJ/mol and along the pore-free path are −355 kJ/mol. Therefore, the difference in free energy for this component between pore and pore-free path is ∼110 kJ/mol. This large negative contribution of lipid to the total PMF arises from the interaction between negatively charged phosphate and positively charged peptide. This is corroborated by analysis of the number of phosphate groups around the shell of peptide shown in Figure 5c. Once again, we used the distance between phosphate groups and peptide as 0.67 nm for calculating the average number of phosphate groups around the peptide. Counter to the trend in solvation by water, the peptide interaction with negatively charged PO4 groups increases along the path into the bilayer, thus offsetting the loss of hydration water. We also found that the average number of choline groups around the peptide inside the membrane is fairly small and has little influence on the total PMF. We found the ion contribution to the total PMF is ∼180 kJ/mol in all four cases. Decomposition of this component into positively and negatively charged ions revealed that the major contribution is arising from chloride ions. (See the inset of Figure 4.) As the positively charged nonaarginine travels from bulk to interior of membrane, it loses chloride ions and adds destabilizing contribution to the total PMF.

Thus, the previous results indicate that the translocation of positively charged Arg9 from bulk water to the interior region of membrane is not a barrierless process. However, the barrier is lower along a reaction coordinate that includes a pore consisting of water and lipids whose negative charges can stabilize the highly positively charged Arg9. Recently, García and coworkers proposed a transient water pore model that explains the energy-independent transfer of positively charged Arg9 across a lipid bilayer.[29] (Note that in their water pore they also include local lipid molecules.) They also mentioned that pore-forming degrees of freedom (DOF) might belong to the slow-relaxing DOFs and there is a requirement of activation energy for that process. From PMF decomposition analysis, it is observed that the stabilizing effect arising from membrane is not sufficient to overcome the destabilizing effect that originates from water and ion. Finally, our analysis indicates that the conformation of Arg9 peptide has little influence on the PMF profile.

Enthalpy and Entropy Decomposition

We next consider enthalpy/entropy decomposition of the PMF. Enthalpy of each configuration was calculated by taking summation of total energy and mechanical energy (i.e., pressure (P) – volume (V) energy) of the system. We computed average enthalpy of system as a function of COM distance between peptide and membrane. For each distance, we took the average of all production configurations obtained from US windows. We used the “g_energy” utility program implemented in Gromacs to compute the enthalpy. We remove the biasing potential energy for our analysis. For consistency with the PMF, the reference zero enthalpy is taken to be the peptide in bulk water state; thus, we subtract the enthalpy of the bulk water windows (reaction coordinate at z = 6 nm) from components of the membrane system. System entropy relative to the bulk water state is the difference between the PMF and estimated enthalpy along the reaction coordinate.

In Figure 6, we have presented the change in enthalpy, ΔH, and entropy in terms of −TΔS contributions to the total PMF as a function of reaction coordinate for the transferring of linear and cyclic Arg9 from bulk water to interior of membrane as obtained from pore and pore-free paths. The change in enthalpies for transferring linear and cyclic Arg9 along the pore-forming path are ∼600 kJ/mol. Similarly, the obtained values for linear and cyclic Arg9 along the pore-free path are about 810 and 656 kJ/mol. Thus, the major contribution to the barrier of PMF profile arises from the enthalpy component for all cases. Additionally, we did not find any significant difference in enthalpy profile between linear and cyclic Arg9 for the transferring along pore-forming path. Furthermore, although the change in enthalpy profile for cyclic Arg9 along pore and pore-free paths shares similar features at the bulk water and interface region, there are significant differences at the pore forming region. We observed that the change in enthalpy profile at the pore-forming region for pore-forming path is nearly flattened, which is arising due to the formation of a pore. In the previous section, we argued strong interaction of Arg9 with water and negatively charged phosphate groups of lipid molecules are the origin of flattened of PMF profile. This is once again validated from the present enthalpy calculation. The enthalpy of linear Arg9 is monotonically increasing along the pore-free path. However, enthalpy of this transferring process decreases as the peptide moves toward the center of the bilayer. Changes in entropies in terms of −TΔS show that it has a favorable contribution to the PMF for all cases. In particular, the values of −TΔS for transferring linear and cyclic Arg9 along the pore-forming path are about −367 and −387 kJ/mol, respectively. Similarly, the obtained values for linear and cyclic Arg9 along pore-free path are about −480 and −327 kJ/mol, respectively. Decrease in enthalpy at that region occurs due to the formation of pore at the central window. Because the sign of −TΔS is negative, hence the net change in entropy along the reaction coordinate is positive. The increase in entropy change along the reaction coordinate is rationalized by consideration of the increase in microstates upon membrane deformation and rearrangement of water molecules and ions between bulk water and membrane phases of the system. However, the gain in −TΔS is not sufficient to offset the increase in enthalpies for all cases.

Figure 6

Enthalpic ΔH and entropic −TΔS contributions to the total PMF obtained from the four simulations: pore-forming path of linear and cyclic Arg9 and pore-free path of linear and cyclic Arg9. A vertical offset of 100 kJ/mol is added for clarity.

We further decompose the change in total enthalpy (ΔH) into different components of the system. Because the mechanical energy (i.e., PΔV work) and the change in bonded energies of each component for the transferring process are negligible and the total kinetic energy of the system is constant (i.e., the simulations are performed at constant temperature), it is considered that the enthalpic contribution of each pair of component is arising from the nonbonded interactions between them. Therefore, this is done by the splitting of total interaction energy into four self and six cross-pair interaction energies. We have displayed the change in enthalpic contribution of the different components of the system for pore and pore-free paths in Figure 7. The Figure shows that membrane–membrane, water–water, and water–peptide interaction energies increase dramatically, while membrane–water and membrane–peptide interaction energies decrease for all the cases as the Arg9 travels from bulk water to interior of membrane. Other interaction energies have minor contribution to the total enthalpy, which is shown in Figure 7. The details of each contribution can be found in SI Figure S17. Once again, we observed that the change in enthalpy profile at the pore-forming region for pore-forming path is nearly flattened. Moreover, we noticed that peptide–peptide interaction energy does not contribute anything to the total enthalpy change in the course of translocation. Thus, the entire free energy of destabilization mainly arises from the membrane–membrane, water–water, and water–peptide repulsive interactions. We found that such destabilization interaction energies are around 1700, 800, and 1100 kJ/mol, respectively, for all cases. The result is quite consistent with our expectation. As the Arg9 moves toward the center of the bilayer, the repulsive interaction among the lipid molecules increases due to the severe deformation of two leaflets of DPPC bilayer from its perfect form. Moreover, as the peptide travels toward the center of bilayer, a few water molecules move inside membrane along with the peptide during the transferring process. Therefore, these water molecules experience a penalty of water–water interaction. The water–water repulsive interaction also increases because configuration of such water molecules present inside the membrane is far from equilibrium.

Figure 7

Decomposition of the change in total enthalpy into 10 different components of the system: interaction energies of (a) membrane-membrane, (b) water–water, (c) water–peptide, (d) membrane–water, (e) membrane–peptide, and (f) the remainder, which is the summation of ion–membrane, ion–water, ion–peptide, ion–ion, and peptide–peptide interaction energies. Details of each contribution in the remainder are shown in SI Figure S17. A vertical offsets of 500 kJ/mol is added for clarity. The symbol coding scheme adopted in the Figure is the same as in Figure 4

Increase in peptide–water interaction along the reaction coordinate can be explained by loss of peptide solvation energy. This is evident from the analysis of average water number along that reaction coordinate (Figure 5a). Membrane–water interactions contribute stabilization enthalpy to total PMF. This is because as the peptide travels toward the center of membrane, it deforms the membrane and consequently negatively charged phosphate groups of lipid molecule become exposed to water. Similarly, the membrane–peptide interaction is becoming favorable as the peptide moves toward the center of bilayer. We obtained a value of about −1000 kJ/mol. This stabilization appears to be due to the strong interaction between negatively charged phosphates groups of lipid molecules of deformed bilayer and positively charged Arg9.

The above analysis suggests that the barrier in PMF is dominated by the enthalpic contribution of the system. Moreover, our analysis revealed that water–peptide interaction is almost compensated by the membrane–peptide interaction, and net enthalpy change for the translocation process arises mainly from the uncompensated self-interaction energy of water and lipid molecules and the cross interaction between water and peptide. Furthermore, we learned that entropic contribution of the system is favorable to the total PMF. We argued that such favorable entropy arises due to the membrane deformation and rearrangement of ions and water molecules between bulk water and membrane phases. Moreover, the enthalpic and entropic contributions to the total PMF are nearly flattened at the pore-forming region which suggests that the formation of a pore inside the bilayer changes the enthalpy and entropy barriers abruptly.

Nonadditivity

In this section, we consider thermodynamics of translocation of a second Arg9 peptide into a pore created by the prior insertion of Arg9 into the bilayer center. In particular, we have studied this issue only for cyclic Arg9. It has been previously observed in all-atom molecular dynamics simulations that the transfer of an additional monoarginine molecule into an existing bilayer-peptide system occurs via sharing the water defect created by first monoarginine.[25]

Figure 8 shows snapshots from US MD simulations with one Arg9 restrained at bilayer center and the second Arg9 at bulk water, interface, and membrane center regions in the presence of a pore created by the first one. Although the structure of the pore is not much affected by the second Arg9 at the bulk water and interface regions, the size of pore increases significantly at the interior region of membrane (shown in Figure 9). We roughly estimated the radii of the pore for the single and double Arg9 systems, and the calculated values are about 1.2 and 2.5 nm.

Figure 8

Snapshot of configurations for the transferring of second peptide with the first one restrained at the center of the bilayer: (a) bulk water (6.0 nm), (b) interface (2.0 nm), and (c) center (0.1 nm). Color coding scheme adopted in the Figure is the same as in Figure 2

Figure 9

Two-dimensional normalized density maps of the lipid beads for the single (left) and double (right) Arg9 systems are shown in the Figure. The Figures were obtained by considering the central windows of the US sampling simulations.

We further estimated PMF for the transferring of second cyclic Arg9 from bulk water to the center of membrane in the presence of the first one placed at the center of bilayer. The results are presented in Figure 10a. For comparison, we also included the PMF profile for transferring single Arg9 from bulk water to the center of bilayer. The free-energy cost for the transfer of a second Arg9 from bulk water to the center of membrane is much lower than the first one. The height of barrier has been reduced by a factor of 2. As the second peptide moves toward the center of the bilayer, it shares the pore formed by the first one. Decrease in free-energy barrier for the transferring of second Arg9 occurs due to the increase in stabilizing interactions of peptide with water and negatively charged phosphate groups of lipid molecules present inside the pore. This stabilization energy is further enhanced due to the increase in pore size when the peptide is approaching toward the center of bilayer. The result recapitulates the observations of MacCallum and coworkers using atomistic simulation.[25]

Figure 10

(a) Red solid line (diamond symbol) shows the PMF of transferring a second cyclic Arg9 with the first one restrained at the center of the bilayer, black solid line (circle symbol) shows the PMF of single Arg9 along pore-forming path, and blue dashed line (no symbol) shows the PMF of single Arg9 along the pore-free path. (b) Total PMF of two Arg9 system is decomposed into the different contributions arising from water, ions, membrane, and first and second peptides. The inset shows the contribution of sodium and chloride ions.

The PMF for transferring second Arg9 was again decomposed by the same method as we used for the single Arg9. The results are displayed in Figure 10b. Similar to single Arg9 (systems A and B), the Figure once again shows that water and ions keep destabilizing contributions (519 and 37 kJ/mol), whereas membrane has stabilizing contribution (−536 kJ/mol) to the PMF for the transferring of the second Arg9 across the membrane and they nearly compensate each other. However, the contributed values are relatively smaller as compared with PMF of first Arg9. Such destabilizing contributions mainly arise from the desolvation and deionization of peptide, which is in agreement with the results obtained from the transferring of single Arg9 to the membrane. The stabilizing contribution arising from membrane was further enhanced for second Arg9 because the existing pore reduces the cost of deforming membrane. Interestingly, we found that the first Arg9 has destabilizing contribution to the PMF of second Arg9 (104 kJ/mol). This is not surprising, and it mainly arises due to the repulsion between two peptides as the second Arg9 moves toward the center of membrane. Thus, our decomposition analysis suggests that repulsion between two peptides has an important role on the PMF barrier of a second peptide.

Summary

We have studied translocation thermodynamics of a positively charged Arg9 peptide across a DPPC bilayer using US molecular dynamics simulations coupled to pairwise additive CG Martini force fields for all component–component interactions. We have computed PMF in the presence and absence of water/lipid pores. We note that the formation (or dissolution) of a pore in the center of a bilayer is an example of a slowly evolving orthogonal degree of freedom. We thus stress that the calculations we perform are free-energy profiles along local paths that do not sample both pore and nonpore states. In particular, we explored the associated changes in free energy, enthalpy, and entropy. Focus has been made on the formation of transient pore inside the membrane that changes the landscape of free energy profile. The potential of mean force shows a significant barrier to translocation, on the order of 240 kJ/mol, in the presence of a lipid/water pore that bathes the translocating peptide as it crosses the bilayer. There is no PMF minimum at the water–bilayer interface observed. The tight binding with membrane may require anionic lipid composition.[51] The absence of a minimum at the water–bilayer interface seems to contradict all-atom studies that invariably predict some degree of stability of polar or charged solutes at this interface. We stress that based on our current work as well as in the context of current knowledge and state of the art in modeling membranes at various resolutions, we are not in a position to argue whether the presence or absence of the interfacial minimum is capable of discerning the accuracy or reliability of current force fields. For the particular case of charged CPPs, there is experimental literature indicating that the binding of such highly charged species to bilayers is facilitated by some degree of anionic character;[52−54] binding to model bilayers/membranes that are zwitterionic is much weaker (or nonexistent) in experiment. Thus, taken in the context of this information, the current force field does seem to capture the behavior qualitatively in that we do observe (results not shown) that the addition of anionic lipids to the bilayer leads to emergence of an interfacially stable state. Decomposition of the PMF indicates that although there is significant entropic stabilization arising from the enhanced configurational entropy exposing more states as the peptide moves through the bilayer, the enthalpic loss (as predicted by the interactions of this CG model) far outweighs any former stabilization, thus leading to significant barrier to translocation. The underlying physicality of this result is the severe net desolvation of a highly charged solute. The present results qualitatively recapitulate observations from all-atom MD simulations of similar systems. Furthermore, in comparing the translocation process in the presence and absence of a pathway that accommodates the formation of a long-lived, stable pore, we find, as in previous studies, that the pore-forming pathway is the lower free-energy pathway (lowered by an amount of 90 kJ/mol). In this sense, the CG model further reinforces the conclusions of the all-atom simulations.[29] Nevertheless, the high barriers predicted contradict the efficient internalization of this peptide into cells and model GPMV.[14] This suggests several areas of concern regarding force-field calculations and the interpretation of such calculations. First, because of the lack of atomistic detail in CG model, it may not reflect the actual energy scales found in nature. Furthermore, system size effects may contribute to free energetics of internalization. As indicated by Hu et al.,[30] system size effects are significant in all-atom force field based calculations of PMF for arginine translocation in model PC bilayers. Of course, the effects of system size are difficult to assess, particularly with the CG models studied here as the reaction coordinate chosen a priori become degenerate when used in systems in conjunction with larger lateral system dimensions. Another important factor such as the composition of lipid is missing in our model. Recently Ciobanasu et al.[21] have shown that including anionic phosphatidylserine, which is common in eukaryotic cells,[55] in phosphatidylcholine GUV facilitates the CPP uptakes. Thus, we reasoned that incorporating anionic lipid in our model may further reduce the free-energy barrier for the CPP translocation. Currently, this research is ongoing in our laboratory. Moreover, our CPP is composed entirely of arginine amino acids, is charged and highly hydrophilic, and may not undergo translocation through PC bilayer. However, experiments show that the introduction of hydrophobic groups such as fluorophores, tryptophan, or hydrophobic sequence HAtag (YPYDVPDYA) in polyarginine induces the translocation.[56,57] This could be another possible reason why we predict a high translocation barrier.

Additionally, we observed that not only water but also phosphate groups of lipid molecules and Cl– ions are present inside the pore, lowering the PMF in the pore-containing region along the reaction coordinate. Furthermore, we found that once the pore is formed, the free energy, enthalpy, and entropy of the system are nearly constant, implying that the system is in a metastable state. This result also indicates that the formation of transient pore and diffusion are the key steps for the energy-independent translocation of CPP across a membrane. Because the formation of water pore involves only one of the many slow-relaxing DOFs, exploring other slow-relaxing degrees of freedom should further lower the free energy. Decomposition of PMF analysis revealed that the stabilizing effect arising from membrane is insufficient to overcome the destabilizing effect originating from water and negatively charged chloride ions. We note that chloride ions do show a condensation effect because they are more probable near the membrane and peptide compared with when there is no peptide in the system. (See Figure S20 of the Supporting Information.) Moreover, we found that the barrier in PMF is dominated by the enthalpic contribution of the system. Such destabilizing enthalpic contribution mainly arises from the self-interaction energy of water and lipid molecule and cross-interaction between water and peptide, which was evident from our analysis. We found that the process is favored by the entropic contribution of the system. The contribution of such favorable entropy mainly arises due to the membrane deformation and rearrangement of ions and water molecules between bulk water and membrane phase. Additionally, our analysis revealed that the cost of transferring of an additional Arg9 in the presence of pore formed by first one is minimal. This indicates that the cost of forming a pore for the first Arg9 is the dominant energetic factor. These findings complement previous studies on the transferring of arginine-rich peptide across the model membrane and provide further support to a complex picture of the translocation process.