Reconciling structural and thermodynamic predictions using all-atom and coarse-grain force fields: the case of charged oligo-arginine translocation into DMPC bilayers.

Using the translocation of short, charged cationic oligo-arginine peptides (mono-, di-, and triarginine) from bulk aqueous solution into model DMPC bilayers, we explore the question of the similarity of thermodynamic and structural predictions obtained from molecular dynamics simulations using all-atom and Martini coarse-grain force fields. Specifically, we estimate potentials of mean force associated with translocation using standard all-atom (CHARMM36 lipid) and polarizable and nonpolarizable Martini force fields, as well as a series of modified Martini-based parameter sets. We find that we are able to reproduce qualitative features of potentials of mean force of single amino acid side chain analogues into model bilayers. In particular, modifications of peptide-water and peptide-membrane interactions allow prediction of free energy minima at the bilayer-water interface as obtained with all-atom force fields. In the case of oligo-arginine peptides, the modified parameter sets predict interfacial free energy minima as well as free energy barriers in almost quantitative agreement with all-atom force field based simulations. Interfacial free energy minima predicted by a modified coarse-grained parameter set are -2.51, -4.28, and -5.42 for mono-, di-, and triarginine; corresponding values from all-atom simulations are -0.83, -3.33, and -3.29, respectively, all in units of kcal/mol. We found that a stronger interaction between oligo-arginine and the membrane components and a weaker interaction between oligo-arginine and water are crucial for producing such minima in PMFs using the polarizable CG model. The difference between bulk aqueous and bilayer center states predicted by the modified coarse-grain force field are 11.71, 14.14, and 16.53 kcal/mol, and those by the all-atom model are 6.94, 8.64, and 12.80 kcal/mol; those are of almost the same order of magnitude. Our simulations also demonstrate a remarkable similarity in the structural aspects of the ensemble of configurations generated using the all-atom and coarse-grain force fields. Both resolutions show that oligo-arginine peptides adopt preferential orientations as they translocate into the bilayer. The guiding theme centers on charged groups maintaining coordination with polar and charged bilayer components as well as local water. We also observe similar behaviors related with membrane deformations.

Introduction

Translocation of a variety of chemical species across physiological cellular membranes has garnered significant attention over the last several decades.[1−5] In particular, a broad class of peptidic systems termed cell-penetrating peptides (CPPs) has generated somewhat of a controversial dialogue in the literature revolving around the exact mechanism of cell penetration.[1,4−6] Major pathways suggested include purely diffusive versus cell-mediated (energy- or ATP-associated) mechanisms.[7] The broader notion of transfer of chemical species across the membrane–water interface has been considered extensively in the recent past. For example, following seminal work on structure determination of a voltage-gated ion channel, KcSA,[8,9] tremendous effort emerged to understand the reasons for charged residues, such as arginines, in the context of what was considered the hydrophobic membrane center.[10−18] The nature of thermodynamic scales to characterize amino acid propensities toward the bilayer interior has enjoyed a rich history as well.[19,20]

In conjunction with experiment, molecular simulations have played a major role in establishing thermodynamics of partitioning and translocation of a variety of species into model bilayer membrane systems.[19,20] Molecular dynamics simulations coupled with the latest all-atom (AA) and coarse-grained (CG) force fields are routinely used to compute potentials of mean force (PMFs) for numerous peptides and small molecules, thus providing thermodynamic insights regarding stability and barriers to the translocation process.[21−33] With regard to coarse-grained force fields, the Martini coarse-grained force field for lipid bilayers[34−36] and amino acids[37,38] has become a routine element of the molecular modeler’s toolkit for studying long time scale and large length scale processes in lipid bilayer systems. The Martini force field is one of several available coarse-grain force fields as discussed in a recent review.[39] The Martini force field was recently updated to reflect the ability to treat polarizable water[40] as well as charged amino acids and ions.[37,38] In particular, the latest version of the Martini force field was reparameterized and validated against the Wimley–White peptide partitioning scale, with recent modifications to the force field providing improved agreement between computed PMFs and side-chain transfer/partitioning free energies with the Wimley–White scale.[41] In particular, the results for the transfer/partitioning free energies for charged arginine and lysine were improved dramatically (refer to Figure 7 of ref (41) and Table 1 of ref (37)). This advance in the force field is relevant to the study of charged peptidic systems that are experimentally suggested to translocate across cell membranes and model bilayers. With the ability to sample larger systems in a more efficient manner, CG models afford a way to study larger polymers in the membrane context.

Since the study of translocation thermodynamics via molecular dynamics (MD) simulation can be accomplished via computation of PMFs (i.e., the reversible work necessary for peptide transfer into the bilayer from solution under equilibrium conditions), numerous studies have appeared to address the problem.[21,22,24−29] Treatment of water in such systems is crucial, as water molecules near the lipid–water interface make a major contribution to the electrostatic potential. An improved polarizable MARTINI water model has been gained by introducing two charged sites into the existing one bead nonpolarizable MARTINI model. Recent studies suggest that the polarizable MARTINI model leads to improved results for the insertion PMF of charged oligopeptides into lipid bilayers but it does lack the ability to reproduce the interfacial dipole potential[40] predicted by all-atom force fields; factors contributing to this deficiency may be the reduced polarizability of the water model (leading to insufficient positive contribution to the dipole potential from water to offset the negative contribution from the membrane components) or lack of specific lipid electrostatic components such as glycerol dipoles which may counterbalance the net negative dipole potential. In a relatively recent study, use of big multipole water (BMW)[42] with MARTINI lipid resulted in a membrane dipole potential that is in good agreement with experimental and AA results.[43] With regard to highly positively charged peptides, however, relatively few studies have probed the issue of translocation thermodynamics using coarse-grained models.[22,26,30] Here, we consider the application of the Martini CG force field for a relevant distinction between the current Martini polarizable CG and most AA force fields (FFs) lies in the predicted PMFs describing free energetics along an order parameter (OP) connecting two states: (1) an oligo-arginine peptide in bulk aqueous solution and (2) the peptide at the center of a model lipid bilayer membrane. A significant qualitative difference is the lack of a bilayer–water interfacial free energy well in the PMFs predicted by the CG model. The interfacial well appears to be a general characteristic of all-atom based PMFs, and arises from an intricate balance of dispersion and electrostatic forces implicitly embedded in current AA force fields. Since all-atom force fields by design include a large set of partial atomic charges, the electrostatic interactions between translocating peptides and the highly charged headgroup moieties lead to a significant free energy bias for stable interfacial states. The Wimley–White scale is valuable in this sense, since it is very difficult to parametrize/calibrate force fields to bulk systems that represent the membrane–water interface. In general, force fields are calibrated by focusing on free energy differences between aqueous and nominally hydrophobic environments (to mimic the water and membrane/bilayer center end points). Because of the difficulty of isolating a pure “interfacial” phase, calibrating force fields to match partitioning free energies to this end point is not as straightforward. Interfacial binding free energies are usually not parametrized into major force fields, though the Martini force fields have included this experimental observable as an important property for force field calibration. Thus, the quantitative nature of force field based interfacial free energy wells remains somewhat ambiguous. Nevertheless, predictions of AA force fields are viewed as possible benchmarks for calibrating CG force fields, particularly in terms of free energetics and structural elements of simulations.

Presently, our intent is to explore the extent to which the Martini-based CG force field is able to recapitulate free energetic predictions based on AA PMFs of the translocation of three oligo-arginine peptides across a model DMPC bilayer. Our intent here is not to judge the validity, appropriateness, or quality of any of the current models but rather to explore whether qualitative features of predicted PMFs can be matched to those predicted by AA models (particularly for the challenging systems involving highly charged molecular species as necessary in the study of cationic cell-penetrating peptides), and whether further quantitative agreement is possible. We achieve this goal by introducing modifications to the current generation of the Martini force field embodied in modified interactions between the translocating solute and the bilayer and water components. We acknowledge that elements of the water model (i.e., representation of water electrostatic distribution and polarizability/flexibility) and the nature of electrostatic elements of the coarse-grain lipid model (glycerol dipoles are omitted in the current form of the Martini model) are also factors for consideration in attempting to reproduce the form of the PMF curves. In this work, however, we choose to consider whether maintaining the current form of the model coupled with further refinement can at the very least yield similar qualitative free energy profiles. In particular, since the Martini water and lipid force field combination yields a net negative interfacial dipole potential, a positive charge should be stabilized by the presence of the large −2 V potential;[40] since previous literature studies report no significant interfacial minimum and large barriers to charged species translocation into Martini lipid bilayers, we hypothesize that nonelectrostatic components of the force field may contribute significantly as well to the shape of the computed PMFs. We note that higher order electrostatic moments are not considered as potentially counterbalancing the effects of the dipole potential at the interface, thus leading to higher barriers for charged species. Thus, we address the nonelectrostatic parts of the force field in this study, again, acknowledging that this is a complex problem with numerous other factors contributing to observed differences within the current form of the Martini force field. This work is relevant, since models allowing access to higher time and length scales in molecular simulations of this type of system are needed.

Here, we consider free energetics and kinetics aspects related to force field parameters for the translocation of a short cationic oligo-arginine peptide across a model DMPC membrane using free energy simulation techniques. In particular, we compared the results obtained from all atom (AA), Martini nonpolarizable coarse grained (nonpol-CG), and Martini polarizable coarse grained (pol-CG) models in this study. The rest of the article is organized as follows. In section 2, we describe in brief the setup of the systems and the simulation methods employed. The potentials of mean force (PMFs) that are obtained from three different models (AA, nonpolarizable CG, and polarizable CG) for three oligo-arginines are presented and discussed in the following section (section 3). Decomposition of the total PMF into the system’s components is also discussed in this section. Further, the time scale of such a translocation process has also been explored and presented by estimating the rate constants in that section. A summary of important findings and conclusions are highlighted in section 4.

Methods

All Atom (AA) MD Simulation

We consider translocation of mono-, di-, and triarginine peptides into a model DMPC bilayer, with the aim of computing atomistic force field-based PMFs along an OP spanning configurational states with peptide fully solvated in aqueous solution to those with the peptide residing in the bilayer center (OP to be defined below). The membrane system was constructed in a rectangular box of initial dimensions 93.5 Å × 93.5 Å × 120 Å with 288 DMPC molecules (144 molecules per leaflet), 24 635 water molecules, one cationic oligo-arginine peptide (charge = +1, +2, +3), and sufficient numbers of chloride ions to neutralize the overall system. The peptides are patched at ends with an NH2 (CT2) group at the C-terminus and an acetyl group (ACE) at the N-terminus; we chose to only consider systems where the charge is constrained to reside only on the side chains. The peptides in each system were initially placed about 15 Å away from the membrane/water interface. The membrane center of mass was centered at (x = 0.0, y = 0.0, z = 0.0). All simulations were performed using NAMD 2.93b[44] with CMAP corrected CHARMM22 all atom force fields for peptides, CHARMM36 for lipids,[45−49] and the rigid TIP3P model for water.[50] Following standard procedures for equilibration, the systems were first minimized using the conjugate gradient energy minimization method and then equilibrated at constant temperature (T = 303 K) and at constant pressure (P = 1 atm) ensemble (NPT) for about 10 ns.

For all systems, the cutoff for van der Waals (VDW) interactions was set to 12 Å with smoothing functions activated from 10 to 12 Å. The pairlist distance was set at 14 Å. The conditionally convergent long-range electrostatic interaction was modeled by using the particle mesh Ewald (PME) method.[51] The simulations utilized SHAKE for constraining the hydrogen–hydrogen and hydrogen–oxygen distances of water.[52] Nosé–Hoover barostat[53,54] and Langevin dynamics methods were used to maintain constant pressure and temperature for each simulation. A time step of 1 fs was used to integrate the equation of motion.

The adaptive biasing force (ABF) module of the NAMD program along a chosen OP[55,56] was used to compute the PMFs of the translocation of oligo-arginine peptide into the model DMPC membrane. Here, we chose the z component of the distance vector formed between the center of mass of the peptide and a dummy atom located at the Cartesian position (x = 0.0 Å, y = 0.0 Å, z = 0.0 Å) close to the center of mass of the DMPC bilayer as an OP for this study. The use of the dummy atom is made essential due to the algorithmic constraints of the NAMD software. We note that there is little difference between the use of a membrane center of mass to peptide center of mass (com–com) OP for the AA system and an OP defined by the distance from a dummy atom placed at the absolute Cartesian point (0, 0, 0). The z-component of the bilayer to peptide center of mass distance correlates very well with the distance from the point z = 0 to the peptide center of mass, as shown in Figure S1 of the Supporting Information. For each system, a total of 16 windows ranging from OP values of −1.5 to 40 Å at a spacing of 3 Å were constructed along the OP. Additionally, we performed six replicas of each ABF simulation by changing the initial velocities of each atom in the system. Each window was simulated for about 100 ns by using the above protocol. The details of the system setup and ABF simulations are outlined in our previous work.[29]

To verify whether the current combination of solvent, ion, peptide, and lipid force fields can facilitate unbiased translocation of the peptides, we performed MD simulations of each peptide in solution. The simulations were performed using a harmonic boundary potential wall on one edge of the box (in the z-dimension) in order to keep the peptide positioned in one-half of the central simulation cell (this effectively reduces the sampling volume without incurring loss of generality of the observations related to spontaneous translocation of the peptide). Specifically, we set an upper wall at 30.0 Å for Arg1, 32.0 Å for Arg2, and 35.0 Å for Arg3 along the positive z direction with a force constant of 25.0 kcal/mol/Å2, and the width of the upper boundary wall was extended by 0.5 Å; these specifications are per the NAMD wall potential implementation. After a 10 ns equilibration period, the systems were continued for 50 ns production. The final production trajectory was used to analyze the spontaneous binding between oligo-arginine and DMPC lipid as well. The average minimum distance between the peptides and bilayer interface was about 2.5 Å. The values are pretty small and may be assumed as the upper limit for the binding between peptide and bilayer. We also noticed that such minimum distances fluctuate very little around their average values over the simulation period (see Figure S2 in the Supporting Information). Thus, our results suggest that, once the oligo-arginine binds with the bilayer interface, the dissociation event rarely occurs. Although the oligo-arginine spontaneously binds with the bilayer interface, we did not observe a direct translocation event in this simulation.

Coarse Grained (CG) Simulation

General Molecular Dynamics Protocol

The Martini CG model as developed by Marrink et al.[37,38,40] was used to simulate the lipid–peptide–counterion–water interaction. We used both the polarizable (version 2.2P[37,40]) and nonpolarizable (version 2.2[37,38]) Martini force fields for the oligo-arginine and water models. A widely used standard P5 bead was chosen to map the neutral terminal residue and backbone of CG oligopeptides as in the recent work by Marrink and Tieleman.[41,57,58] We used the nonpolarizable Martini DMPC and ion model (version 2.0[34−36]) for this simulation. The CG simulations were performed using the MPI supported GROMACS software package (version 4.6)[59,60]

The initial structures of three CG systems were constructed by converting the above well equilibrated all atom systems using the VMD tool, CGbuilder plugin, version 0.1. After this conversion, the numbers of particles of the systems are almost reduced by a factor of 12 for nonpol CG and a factor of 4 for pol-CG. The simulation cells consist of a rectangular box of dimensions around 9.3 × 9.3 × 12.0 nm3, yielding about a 3.6 nm thick slab of lipid molecules surrounded by bulk water and ions. The components of the system are depicted in Figure 1. Following the conversion of AA systems to coarse-grained representations, the systems with peptide in bulk solution were minimized by using the steepest descent method, followed by an equilibration run at 1 atm pressure and 303 K for 500 ns in the NPT ensemble. During the MD equilibration, the area per lipid equilibrated to the values of 61.11 and 60.04 Å2 in agreement with published results[40] for the polarizable and nonpolarizable force fields (see Figure S3 in the Supporting Information for details).

Figure 1

Molecular structures of water, DMPC lipid, Arg1, Arg2, and Arg3 in the all atomic model, coarse-grained nonpolarizable model, and polarizable model. In the all atomic model, oxygen, hydrogen, carbon, nitrogen, and phosphate atoms are colored in red, gray, green, blue, and purple, respectively. In the coarse-grained nonpolarizable model, one bead W represents four atomic water molecules in the all-atom model. Choline (NC3), phosphate (PO4), carbonyl (GL1, GL2), and tail beads are colored in blue, purple, red, and green, respectively. Arginine amino acids are represented by three beads, one backbone bead BB (red), and two side chain beads SC1 and SC2 (green). In the coarse-grained polarizable model, W, WP, and WM are neutral, positively charged, and negatively charged beads. The off-center charge model is used for the charged CG arginine. Backbone (BB), noncharged side chain beads (SC1, SC2), and charged (SCP) beads of CG arginine peptide are colored in red, green, and blue, respectively.

We used a time step of 20 fs and updated the neighbor list every 10th step. The Lennard-Jones (LJ) and electrostatic (Coulomb) interactions were calculated by using a 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 conditionally convergent long-range electrostatic interactions were modeled by using the PME method with a fourth-order spline and a 0.12 nm grid spacing. The relative dielectric constants were set to 2.5 and 15.0 for use with the polarizable and nonpolarizable water force fields. To maintain a temperature of 303 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.[61] We employed two temperature coupling groups: water and ions were considered as one, and DMPC and peptide were set as the second group. To keep the bilayer in a tensionless state, periodic boundary conditions with a semi-isotropic pressure coupling algorithm with a 3.0 × 10–4 bar–1 compressibility were used. The LINCS algorithm[62] was used to apply the bond constraint present in Martini force fields.

Umbrella Sampling in CG MD Simulation

To obtain a PMF for the transfer of oligo-arginine in each system, we use 41 umbrella sampling (US) windows that range from 0.0 to 4.0 nm at a spacing of 0.1 nm along our chosen OP. The same OP is the z-dimension distance between the center of mass of the peptides and bilayer. The use of this ostensibly different OP for the coarse-grained systems warrants caution, particularly in the context of comparing the AA and CG results further below. We address this by noting that there is little difference between the use of a com–com OP for the AA system and an OP defined by the distance from a dummy atom placed at the point (0, 0, 0). The z-com of the bilayer correlates very well with the point z = 0, as shown in Figure S1 of the Supporting Information.

We first considered generating initial configurations in the windows along the specified OP by growing an oligo-arginine in the bulk water of the above equilibrated systems, and further equilibrating the peptide–bilayer–water–ion system for about 200 ns after growing in the aqueous phase. We performed a 200 ns production run for each system to study the binding between oligo-arginine and DMPC. In order to prevent unnecessary drift of membrane, a position restraint along the z dimension with a force constant of 1000 kJ/mol/nm2 was applied on the charged groups of the lipid molecule (NC3, PO4) during the growing-in phase for all the simulations. The growing of peptide inside the system 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 a soft-core potential was used to prevent bead overlap. In the following step, we slowly grew in the Coulomb interactions using the same protocol. Each window was simulated for about 200 ns, and the total simulation time period is about 8.2 μs. For US MD simulations, we applied a harmonic potential with a force constant of 1500 kJ/mol/nm2 to restrain the peptide at each window. The details of the window setup and US method have been described in detail in our recent work.[30]

The weighted histogram analysis method (WHAM) was used for postsimulation unbiasing of umbrella sampling data.[63] We use the Gromacs tool “g_wham” to generate the final PMF. Order parameter histograms in adjacent US windows exhibited sufficient overlap.

Potentials of Mean Force and System Component Contributions

Taking our order parameter, η = rc.o.m.protein – rc.o.m.bilayer, as the z-component of the vector between the center of mass of the peptide and the center of mass of the lipid bilayer, the potential of mean force can be written as (see the Appendix for details):where W(η) is the potential of mean force along the chosen OP, η, η0 is the value of the OP in the reference state (taken to be the peptide in bulk aqueous solution), η′ is the dummy variable of integration, and ⟨F⟩η′ is the average z-component of the total system force on the peptide center of mass, with the average evaluated at a particular value of the OP (z-component of the total system force on the peptide center of mass averaged over the ensemble positions of all other particle positions). This average force is explicitly dependent on the OP.

As shown in the Appendix, we can make use of the fact that the interaction potential energy, U(r), is pairwise additive to write the expression for the contribution to the total PMF from an individual system component, α (i.e., α = water molecules, α = lipids, α = ions), aswhere ⟨Fα⟩η′ is the average z-component of the total force on the peptide center of mass arising from interactions with system component α. The total PMF is then obtained as a sum over the system component contributions:The instantaneous force on a peptide from system component α, Fα, was computed postsimulation by processing the trajectories of each US window using the Gromacs “mdrun” module. We excluded the interactions between the peptide and system components other than α. We follow an approach by Zhang and van der Spoel.[64] The LJ and real space part of the PME interactions were excluded by using the energy group exclusions parameter in the Gromacs input file, whereas the reciprocal space electrostatic interactions were excluded by setting the charges on each of the particles of the other components to zero.[64] In a second approach, we extracted the coordinates of the pair of system components of interest (peptide and component α) from the US trajectories and computed the nonbonded interactions between them. Results obtained with both approaches matched numerically.

The final PMF and its standard error (uncertainty) were estimated by block averaging consecutive 50 ns time periods from the production run of each US window.[28] We ensured that the block size was significantly larger than the correlation time in each umbrella window. The Visual Molecular Dynamics (VMD) package[65] was used to monitor the simulation, visualization, and graphics preparation for this work.

Tuning of CG Force Field

As motivated in the Introduction, a relevant distinction between the current Martini polarizable CG and most AA force fields lies in the predicted PMFs describing oligo-arginine (defined as a polymer of arginine residues along with backbone atoms) translocation free energetics along an OP connecting two states: (1) an oligo-arginine peptide in bulk aqueous solution and (2) the peptide at the center of a model lipid bilayer membrane. A significant qualitative difference is the lack of a bilayer–water interfacial free energy well in the PMFs predicted by the CG model. The interfacial well is a general characteristic of AA based PMFs, and arises from an intricate balance of dispersion and electrostatic forces implicitly embedded in current AA force fields. In general, force fields are calibrated by focusing on free energy differences between aqueous and nominally hydrophobic environments (to mimic the water and membrane/bilayer center end points). Because of the difficulty of isolating a pure “interfacial” phase, calibrating force fields to match partitioning free energies to this end point is not as straightforward. Interfacial binding free energies are usually not parametrized into major force fields. Thus, the quantitative nature of force field based interfacial free energy wells remains somewhat ambiguous. Nevertheless, predictions of AA force fields are viewed as possible benchmarks for calibrating CG force fields, particularly in terms of free energetics and structural elements of simulations. Presently, our intent is to explore the extent to which the Martini-based CG force field is able to recapitulate free energetic predictions based on AA PMFs of the translocation of three oligo-arginine peptides across a model DMPC bilayer. Our intent here is not to judge the validity, appropriateness, or quality of any of the current models but rather to explore whether qualitative features of predicted PMFs can be matched (particularly for the challenging systems involving highly charged molecular species as necessary in the study of cationic cell-penetrating peptides), and whether further quantitative agreement is possible. In the following, we present a modification to the existing Martini polarizable CG force field. We strongly assert that we are not claiming our tuning to be a final force field, though we will demonstrate the applicability of the model for modeling a cationic cyclic arginine nonamer during the latter part of this work. This work is relevant, since models allowing access to higher time and length scales in molecular simulations are needed.

We thus consider PMFs predicted by all atom force fields to have a bilayer–water interfacial free energy minimum as a reference for tuning the CG force field in this study. We systematically tuned the parameters of the Martini polarizable CG force field for the peptide. We did this because the original combination of the polarizable Martini force field for water and an arginine residue (including its backbone bead, referred to as CG-pol-std-P5) failed to predict an interfacial minimum in the PMF. However, using the original polarizable Martini force field for water and arginine side chain beads, Tieleman and co-workers predicted a PMF minimum at the DOPC/water interface. We repeated their study with polarizable Martini CG water and arginine side chain beads and obtained similar results (see Figure 2). Within the context of the Martini force field, this suggests that the combination of the backbone bead (taken to be neutral) along with the side chain representation is not capturing a proper balance between the interactions of the total residue and the rest of the system components. For instance, the observed PMF may be a result of weak interaction between the backbone bead of the arginine residue and the headgroup of the lipid but a strong interaction between the backbone and water. Thus, we focused on the backbone bead and side chain beads of the CG arginine residue and modified relevant interaction parameters. We strongly caution that our choice to manipulate these particular beads is in some sense defined by the complex interplay of all molecular interactions of the force field, thus providing the context dependence of the choice of free interaction parameters one has at their disposal to modify. That is to say, in this particular case, the lack of an interfacial PMF minimum was traceable (at least in this very specific case) to a small, finite set of interactions, the backbone (name BB, type P5) bead and the arginine side chain beads (including bead names SC1 type N0, SC2 type Qd, and SCP type D) in the current version of the polarizable Martini CG force field.

Figure 2

(left) PMFs for side chain analogues across a DOPC bilayer interface. The peak of phosphate density on the top leaflet (which is 21.6 Å relative to the center of the bilayer) is set to 0. The blue line (triangle symbols) denotes the result from the standard Martini 2.2P force field, the black line (no symbols) is obtained from the OPLS force field (the plot is recovered from the Tieleman et al. work by using a plot digitizer tool), and the red line (circle symbols) represents the result by using the CG-Pol-Z5 force field. (right) PMFs of cyclic Arg9 translocate across the DPPC bilayer with the standard polarizable CG force field type P5 and modified z5 CG force field. The error bars represent the standard error of PMFs computed from each 50 ns data block.

Thus, the goal is to attempt to reproduce the interfacial PMF minimum observed from AA simulations. We propose five modification scenarios. We increased van der Waals (vdw) interactions between the backbone bead of arginine (type: P5) and all the DMPC lipid beads (each DMPC has 10 beads, see Figure 1) by 20%, and decreased the vdw interactions between the backbone bead of arginine (P5) and the central bead of polarizable water (type: POL) by 10% (we call this modified force field CG-pol-B5). We increased vdw interactions between the backbone bead and all the lipid beads by 20%, and decreased the vdw interactions between all the arginine beads (both backbone bead and side chain beads) and the central bead of polarizable water by 10% (we call this model CG-pol-Z5). The effective size of the particles in the Martini force field is governed by the LJ parameter σ. We changed the c6, c12 value in the LJ 12-6 potential by the same percentage. Because σ = (C6/C12)1/6, the vdw interactions were tuned without changing the size of each bead (see Table 1). The ϵ = (C6)2/4C12 for the interaction pairs are also included in the table. Moreover, the parameters for electrostatic and bonded interactions were not changed during the tuning process. The PMF for the Z5 side chain analogue across a DOPC bilayer interface is also shown in Figure 2. Compared to standard arginine side chain beads, we found a deeper minimum in PMF at the interfacial region for this new parameter. This initial observation motivates us that the Z5 side chain parameter might capture better peptide–lipid interactions at the membrane–water interfacial region as compared to its standard one. Recently, Gao et al.[66] have shown that hydrophilicity of the peptide backbone in the coil state is overestimated if the standard P5 backbone bead is used to recover the binding event of the oligo-arginine with the membrane surface; we again systematically tuned the standard Martini polarizable parameter of the backbone bead of the arginine residue. Specifically, we changed the parameter of the backbone bead from P5 to P4, P3, and P2 and refer to these modified force fields as CG-pol-P4, CG-pol-P3, and CG-pol-P2 throughout our study. Since arginine is a polar residue, we decided not to change the parameter (type) of the backbone bead from polar to apolar or neutral beads. Detailed descriptions of all the different parameter sets that are used in our study are presented in Table 1.

Table 1

Overview of the VDW Interaction Parameters Assessed in the CG Polarizable Force Fieldsa

Down arrow and up arrow are used to indicate the parameter increased and decreased relative to the default CG-pol-P5 force field, respectively. The units of C6, C12, and ϵ are kJ/mol·nm6, kJ/mol·nm12, and kJ/mol, respectively.

Results and Discussion

Density Profile

We first briefly address structural similarities between the AA and CG simulations. Snapshots of AA, nonpolarizable, and polarizable CG bilayer configurations obtained from corresponding free simulations (NPT bilayer simulations in the absence of peptides) are shown in the top panel of Figure 3. The structures of bilayers do not change within our simulation period, and water does not penetrate inside the core of the lipid. In the bottom panel of Figure 3, we show the partial mass density profiles for the three model bilayers. The partial mass density profiles for the three models are qualitatively similar, showing that water penetrates into the headgroup region and smaller amounts into the hydrocarbon tail group region of lipid bilayers. The densities of water and lipid intersect at around 19 Å, which is common for the three models. However, compared to the AA model, the mass density of the whole system, carbonyl and tail groups of lipid molecules are relatively higher for both nonpolarizable and polarizable CG models. This is because the masses of a DMPC lipid in atomistic and CG models are different (678 amu vs 720 amu). The density profiles we report for both the AA and CG models are in agreement with numerous published literature results, thus providing some validation of the models.

Figure 3

Snapshots and partial mass density profiles for an AA and CG nonpolarizable and polarizable DMPC bilayer systems at equilibrium. Water is shown as red small beads, headgroups as large beads, and lipid tails as thin gray lines. The headgroup includes the choline, phosphate (blue), and carbonyl (brown) groups, and tails include all the acyl chains. Distributions of the system particle mass density for different AA and CG groups of the DMPC bilayer, water, and the whole system, with respect to the bilayer center (Z = 0). A schematic representation of the structure of DMPC mapping is shown in Figure S4 (Supporting Information).

Total Potential of Mean Force

Figure 4 shows the PMFs for the set of original and modified force fields as well as the all-atom result. Panel A is the all-atom result, panel B shows the PMF for the standard nonpolarizable Martini model, and panel E shows the PMFs for the standard polarizable Martini force field. Panels C, D, F, G, and H show PMF results using various modified force field parameters as outlined in Table 1. The PMFs indicate that the energy scales for the AA model and both polarizable and nonpolarizable CG models are comparable which is somewhat nonintuitive, though satisfying. We found the PMFs for the AA and nonpolarizable CG models reflect an interfacial minimum, whereas no such stable state is predicted using the standard polarizable Martini force field. The presence of such a minimum in the PMF at the interface region indicates that the oligo-arginines bind with lipid, which may be crucial for the translocation process,[4,5,67−69] though we admit that the presence of anionic lipids is more likely a stronger determinant to oligo-arginine binding strength. Translocation barrier heights (bulk to bilayer center) predicted by each force field are shown in Table 2. For the AA model, we obtain an interfacial well-depth for mono-, di-, and triarginine cases to be about 0.83, 3.33, and 3.29 kcal/mol, respectively. The corresponding values for the nonpolarizable CG model of three oligo-arginine systems are 1.96, 3.57, and 4.39 kcal/mol. Thus, both the AA and CG models are able to capture the trend of increasing interfacial stability with increasing chain length. Furthermore, both AA and CG models indicate a shift of the location of the interfacial minimum to larger values of the OP, reflecting the influence of solute size on relative separation.

Figure 4

PMFs of translocating single oligo-arginine Arg1, Arg2, and Arg3 into the center of the model DMPC bilayer by using different FF parameters. For clarity, all the PMF curves are offset by 5 kcal/mol, and dashed lines show the free energy reference value in the bulk. The uncertainties are depicted as shadows. They represent the standard error among PMFs computed from a data block size of 30 ns for the all atom model and 50 ns for the CG model.

Table 2

Free Energetic Barriers from Bulk Water to the Bilayer Center and Interfacial Minima in Oligo-Arginine Translocation across the DMPC Membrane (kcal/mol)

	barrier ΔGtotal			interfacial minima
force field	Arg1	Arg2	Arg3	Arg1	Arg2	Arg3
AA	6.94	8.64	12.80	–0.83	–3.33	–3.29
CG-nonpol-std-P5	13.73	17.49	20.92	–1.96	–3.57	–4.39
CG-pol-B5	13.18	17.27	20.81	–0.37	–2.14	–2.28
CG-pol-Z5	11.71	14.14	16.53	–2.51	–4.28	–5.42
CG-pol-std-P5	15.78	21.81	27.15
CG-pol-P4	14.75	20.42	25.52
CG-pol-P3	13.93	19.28	24.03
CG-pol-P2	12.76	17.14	21.34

Since we do not observe interfacial minima in the PMFs predicted using the standard Martini polarizable model (panel E), we consider the relevant force field parameters. As discussed earlier, an interfacial minimum in the translocation PMF for a single arginine side chain (without a backbone bead) is predicted if the arginine residue (excluding the backbone bead) is used as the translocating solute. Once the backbone bead is added, the interfacial minimum vanishes (see Figure 2). Initially, we focused on modifying the interaction parameters of the backbone bead of the arginine residue so as to decrease its overall interactions with the entire system. This approach exploits the fact that the arginine backbone bead is considered an extremely polar particle, and as such, its interactions with polar and charged components are parametrized to be strong. Thus, the standard backbone bead, given the type P5 in the Martini atom-typing scheme, is supra-attractive in its interactions with water, and has a varied mix of interaction strengths with components of the DMPC lipid bilayer (as shown in Table 1). Our consideration initially was to decrease the interaction strength of the P5 type by simply changing the bead type in the Martini.itp file of arginine to other bead types already present in the standard polarizable Martini model. These types are the P2, P3, P4, and P5 beads. By computing the PMFs of peptide translocation using these modified force fields (panels E, F, G, and H for bead types P5, P4, P3, and P2, respectively), we found that this modification did not result in the prediction of an interfacial minimum. In order to reproduce the minimum in the PMF using the polarizable CG model, we further argued that a stronger interaction between oligo-arginine and the membrane components and a weaker interaction between oligo-arginine and water are necessary. Thus, we followed another protocol where we increased the interaction between the backbone bead of oligo-arginine and the membrane and decreased the interaction between water and the backbone bead. Specifically, we increased the backbone–membrane van der Waals and decreased the backbone–water van der Waals interactions by 20 and 10%, respectively (leading to modified force field type B5). Since we did not change the interaction between the membrane and water, the time average area per lipid in this new force field is unchanged. The B5 force field predicted an interfacial minimum for all oligo-arginines (panel C) . The well depths for the three oligo-arginines range from 0.37 kcal/mol to about 2.28 kcal/mol, which is consistent with the AA and nonpolarizable CG models. However, we observed the barrier height of the bulk water to bilayer center barrier is higher for the B5 model compared to the AA model result (see Table 2). To attempt an improvement, we used the B5 model and tuned the arginine residue side chain parameters to reduce their interaction with water. Overall, we increased van der Waals interactions between the backbone bead and all the lipid beads by 20%, and decreased the van der Waals interactions between all the arginine beads and the central bead of polarizable water by 10% (this model is called type Z5). In Figure 2, we have presented the PMF profiles for the transfer of arginine side chain beads across the DMPC membrane by using the AA, polarizable CG P5, and CG Z5 force fields. The CG Z5 model is respectably successful in reproducing the AA PMF. The side chain parameter of Z5 produces a deeper minimum and reduces the water to the bilayer center barrier. Moreover, we comment on how well the Z5 model recapitulates Wimley–White interface partitioning free energies as discussed by Singh et al.[37,41] The relative partitioning free energy (ΔΔGWW) of Wimley–White (WW) peptides WLXLL with respect WLALL were considered as the hydrophobic scale to gauge the relative interface affinity of different amino acids in experiment. We adopted the protocol established by Singh and de Jong[37,41] to calculate the ΔΔGWW for our new parameters. The partitioning free energy for arginine was calculated at the POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine)/water interface using the thermodynamic cycle (see Figure 2 in ref (41)). The thermodynamic integration (TI) method was use to compute the free energy. The derivative of the Hamiltonian is computed and then integrated numerically to obtain the free energy difference along a predefined order parameter describing the transformation. Here, we built a system with 72 POPCs (36 lipids per leaflet), 1200 polarizable water, and equilibrated it at 1 atm and 300 K for 500 ns in the NPT ensemble using periodic boundary conditions. Then, we added the Wimley–White peptide and equilibrated the peptide for 1 ns at bulk water and membrane interface position. The Z distance of the center of mass of peptide and center of mass of membrane was restrained at 4.0 nm for the bulk water window and 1.9 nm for the interface window with a force constant of 1000 kJ mol–1, respectively. The interface is defined as the intersection of the mass density of water and POPC carbonyl groups (GL1,GL2).[41] A total of 21 windows with 0.05 spacing of the couping parameter λ were used to switch off the vdw and Coulomb interactions of the perturbed side chain from 1 to 0. The side-chain beads of arginine in WLRLL have been converted to dummy beads, which contain no vdw or Coulombic interactions. Soft-core interactions were used to avoid particle overlap. We used the Z5 bead type for all the backbone parameters. For alanine, we made a new bead type Z4 by using the same protocol as that for the backbone bead type P5. The Z4 bead type was used for alanine in the full interaction state, and it is converted to the Z5 type in the dummy bead state, as suggested in Singh’s work.[41] Each λ window was run for 500 ns, and the first 10 ns was considered as equilibration. The Gromacs tool “g_bar” was used to estimate the free energy and uncertainty. The free energy of Z5 we found is −2.69 ± 0.16 kJ/mol; the P5 model and experimental value in ref (37) are −2.5 ± 0.3 and −3.4 ± 0.7 kcal/mol. To provide further evidence for the suitability of the Z5 parameter, we calculated the arginine side chain solvation free energy for standard (P5) and modified MARTINI beads (Z5) using the TI method.[70] The details of the simulation protocol are presented in the Supporting Information. The estimated solvation free energies for the standard MARTINI polarizable force field and modified Z5 force field are −62.5 ± 0.01 and −44.56 ± 0.02 kJ/mol, respectively. Clearly, the modified Z5 force field reduced the arginine water interaction which lowers the solvation free energy of arginine. However, compared to the experimental value −44.8 kJ/mol,[71,72] the modified Z5 force field predicts a consistent solvation free energy. Since we changed the interaction parameters of arginine beads in our simulation, there is always the open question about whether the new parameters drastically perturb the existing arginine–arginine interactions. To explore this, we computed the arginine–arginine pairing free energy for standard (P5) and modified MARTINI beads (Z5) using umbrella sampling methods. We investigated three side chain orientations as follows: collinear, parallel, and orthogonal (see Figures S6, S7, S8, and S9 in the Supporting Information).[73] We do not observe a big difference of the Arg–Arg interaction for the original model, and Z5. Such a result convinces us that the inherent properties of the arginine residue are not affected by our parameter modification. Thus, we expect that the combination of both the backbone and side chain parameter of Z5 will produce a PMF which is close to the results obtained from the AA model (see panel D of Figure 4).

In the above PMF analysis, we took the results obtained from AA simulation as a reference and systematically varied the parameters ranging from the AA to modified polarizable CG model. We found a strong interaction between oligo-arginine and the membrane headgroup and a weak interaction between oligo-arginine and the water are necessary to reproduce the AA model. Further, we found that the increase in van der Waals interaction between backbone and lipid beads and decrease in van der Waals interaction between all the peptide beads and the central bead of polarizable water (type: Z5) provide a PMF which is close to the AA model. We further computed the PMF for cyclic nona-arginine (Arg9) translocation using the Z5 parameter set, Figure 2. The Z5 parameter of Arg9 is capable of reproducing a deep minimum at the water–membrane interface in the all-atom simulation,[28] whereas its standard version (P5) is unable to predict a minimum which we explored in a recent study.[30] Further, the height of the barrier is reduced significantly for the Z5 parameter. It was shown that coarse-grained and all-atom force fields are actually able to reflect very similar thermodynamics of arginine peptide translocation into a model DPPC bilayer. We next decompose the total PMFs into contributions from system components (water, lipids, protein, ions) in the next section.

Decomposition of Total PMF

We follow an approach by Zhang and van der Spoel[64] to decompose the total PMFs of the previous section into system component contributions. Figure 5 shows the contributions from the membrane, counterions, and water to the total PMFs. We computed the components of PMF for all of the force fields described in the previous section, but here we have present results from four models: AA, nonpolarizable CG, P5, and Z5 polarizable CG parameters. Results of other models are presented in Figure S10 in the Supporting Information. The summation of contributions from individual components to the total PMFs match with the total PMF obtained directly from the ABF and US methods (see Figure S11 in the Supporting Information). A similar justification of the decomposition approach is given by Zhang and van der Spoel.[64] All force fields predict destabilization of the translocating solute arising from ion and water contributions. This is not surprising, as the negatively charged counterions and water act to effectively solvate the arginines in aqueous solution. The membrane as a whole acts to stabilize the peptides at the bilayer center. Stabilization is effected by membrane deformations, as will be discussed below. Again, we observe a striking similarity in the profiles of the contributions of the individual components along the OP.

Figure 5

Free energetic contributions arising from system components (a–d) ion, (f–i) water, and (j–m) membrane are shown. Each contribution is depicted with four models, such as the AA, nonpolarizable CG, polarizable CG and CG Z5 models from left to right.

Briefly, we address the nature of stabilization and destabilization of the membrane, ion, and water depicted in Figure 5. The membrane undergoes deformation, the head groups of lipid molecules reorient, and the stabilization from the membrane primarily arises from the interaction between negatively charged phosphate and positively charged peptide. This idea is in part supported by counting the average number of phosphate groups present within a distance of 0.67 nm of all the oligopeptide beads or atoms (considered as the first solvation shell[30]), as a function of OP, shown in Figure S12 in the Supporting Information. Destabilization from water is related to the dehydration of highly solvated charged peptide inside the bilayer, illustrated by the average number of water molecules present within a distance of 0.67 nm of all the peptide beads or atoms, as a function of OP, shown in Figure S13 in the Supporting Information. To address the destabilization from chloride anions, we evaluate radial distribution functions (RDFs) between the amino acid residue and the chloride anions shown in Figure S14 of the Supporting Information. We have provided RDFs for the case when the residue is in the bulk water position as well as when it is located at the bilayer center. Our analysis is performed for the well-equilibrated later portions of the MD trajectory for the relevant umbrella sampling windows. We observe that, in the bulk solution state, the residue has a high probability of having first solvation shell anions in its vicinity. This is intuitive, since the positively charged residue has strong electrostatic interactions with the anion. When the residue is at the bilayer center, the residue becomes “desolvated” from the anion. The electrostatic free energy penalty is quite large because of the close interaction of the ions in solutions. If the ions did not have a strong interaction in bulk solution, then there would be a smaller free energy penalty for the residue losing its counterion associations with the chloride anions. Furthermore, we have presented the analysis for the original P5 and the new Z5 parameter sets. Since the new Z5 set has weakened water–residue interactions, the result is that the residue–anion interactions are effectively enhanced (effectively by virtue of the water not being attracted to the residue as strongly); the RDF for the Z5 set shows an enhanced probability of the anions in the residue first solvation shell compared to that for the original P5 set. This is consistent with our observations that the energy penalty for separating the residue (positive charge) from the anions is nontrivial.

We address differences in component contributions to the PMF among different force field models. We notice that the PMF contributions of water and membrane almost totally compensate each other for the AA and polarizable CG models but not for the nonpolarizable CG model. An unambiguous minimum in the membrane contribution is observed at the membrane–water interface for the nonpolarizable CG model, but such a minimum is replaced by an almost plateau region for the AA and polarizable CG models. However, the destabilizing contributions arising from water and ion are not of sufficient magnitude to offset the stability gained from membrane contributions at the interface region for all models. Thus, the origin of interfacial minima for all models is from the stabilizing contribution of the membrane. The PMF contributions of the membrane, counterions, and water for all polarizable CG models are almost identical. We summarize the PMF decomposition data for the three components at the center of the bilayer for the three different models in Table S1 in the Supporting Information. The data show that introduction of parameter Z5 in the peptide reduces the destabilization contribution of water significantly (overall peptide interactions with water are reduced), which ultimately lowers the height of the barrier in the total PMF.

We further decompose contributions to the total PMF from components into van der Waals and electrostatic. The electrostatic and vdw contributions from ions are shown in Figure 6 and Figure S15 in the Supporting Information. Attractive electrostatic interaction between chloride ion and oligo-arginine strongly disfavors the peptide translocation universally. Although the energy scales are somewhat different, the form of the ion contribution along the OP is similar across the models tested. We find, not surprisingly, that the electrostatic contribution of AA and nonpolarizable CG models contribute the most and least destabilizing effect to the total PMF. Since the force depends on the distance between chloride ion and oligo-arginine, we calculated that for each peptide. We find that the minimum distance between the ion and peptide along the OP is very large (average value around 50 Å, see Figure S16 in the Supporting Information) for the central window. Such large distances lead to the destabilization effect of the ions.

Figure 6

(a–d) Electrostatic (ELEC) and (e–h) van der Waals (VDW) contributions of ions with the AA, nonpolarizable CG, polarizable CG, and CG Z5 models from left to right.

The water contribution is also decomposed into electrostatic and vdw contributions, shown in Figure 7 and Figure S17 in the Supporting Information. Similar to the ion contribution, the electrostatic interactions between oligo-arginine and water for the AA and polarizable CG models also disfavor the translocation process. Due to the absence of charge, the nonpolarizable CG water has zero electrostatic contribution to the total PMF. However, the contribution from vdw interactions between water and oligo-arginine favors the translocation process for the AA and polarizable CG models. Using the standard polarizable CG model (type: P5), we obtained that the vdw contribution of water stabilizes the mono-, di-, and triarginine at the center of the membrane by roughly −53.0, −113.0, and −179.0 kcal/mol. The corresponding values are −81.0, −164.0, and −253.0 kcal/mol for the Z5 parameter. The result is consistent with our expectation as we reduced the vdw interaction between peptide and water by 10% going from the P5 to Z5 model. However, the vdw contribution of nonpolarizable CG water disfavors the bilayer center state. This is because the interaction of water with the peptide and membrane for the nonpolarizable CG parameter is explicitly modeled only by vdw interaction. Thus, the vdw interaction between the water and peptide, which is attractive, turns out to be destabilizing along the OP. Further, we notice that the AA and polarizable CG models reproduce qualitatively similar features for both of these interactions.

Figure 7

(a–d) Electrostatic (ELEC) and (e–h) van der Waals (VDW) contributions of water with the AA, nonpolarizable CG, polarizable CG, and CG Z5 models from left to right.

The membrane contribution is decomposed into electrostatic and vdw contributions in Figure 8 and Figure S18 in the Supporting Information. For all models, electrostatic contributions are stabilizing for some part of the domain of the OP. The AA model shows the largest electrostatic stabilization contribution from the membrane. The polarizable and nonpolarizable CG models predict minima in the electrostatic contribution within the interfacial region of the bilayer. This is because the highest concentration of charged species resides at the bilayer–water interface, thus naturally giving rise to this free energy minimum. Once the peptide passes through the interfacial region, the locally coordinated membrane-based charged species decreases and the contribution from the higher charge density in the interfacial region attracts the peptide away from the bilayer center, thus giving rise to the steep increase in the PMF profile; this changes the net membrane contribution to being destabilizing as the force now works against the peptide being within the bilayer. The increase is also observed in the AA model, though there is sufficient electrostatics (charged species) locally coordinated (the all-atom lipids have glycerol groups that also have partial charges which can interact with the positively charged peptide) to maintain the net stability further into the bilayer (the electrostatic minimum occurs around 10 Å in the profile of the AA model). Passing closer to the bilayer center than 10 Å sees the AA model also give rise to a repulsive contribution from membrane components. This contribution has been linked to membrane deformation as well as direct interaction sources.[23] We find the greatest disparity in membrane contribution arising from the van der Waals interactions along the OP. The AA model for Arg1 and Arg2 displays vdw contributions that are stabilizing at the bilayer center, while, for the trimer, the vdw contribution from the membrane is highly destabilizing. For the two polarizable CG models, we observe very similar vdw contributions, both of which are different in nature from the AA model. The two polarizable CG models predict fully stabilizing contributions beginning near the interfacial region and extending the rest of the domain of the OP into the bilayer. This is not surprising, as the membrane–peptide interactions are increased by 20% relative to the original. This is a large driving force for stability of peptides in the bilayer. The nonpolarizable CG model displays yet another unique profile, first conferring stabilization until the peptide reaches the interfacial region and then becoming steeply repulsive (destabilizing) as the peptide resides at the bilayer center.

Figure 8

(a–d) Electrostatic (ELEC) and (e–h) van der Waals (VDW) contributions of membrane with the AA, nonpolarizable CG, polarizable CG, and CG Z5 models from left to right.

PMF decomposition shows that the energy scales of each component for three models are different, but the corresponding scale of the total PMF for three models are almost identical. Such results imply that the forces are well balanced among the different components of the system for all of these models. In general, we found the membrane has a stabilizing contribution, whereas water and ions have a destabilizing contribution to the total PMF for all the model parameters. However, a noticeable difference in membrane and water contributions to the total PMF has been found in the nonpolarizable CG model. Further, our analysis revealed that the Z5 parameter of the peptide reduces the destabilization contribution of water significantly, which ultimately lowers the height of the barrier in the total PMF. Decomposition of each component into vdw and electrostatic contributions provides further insights of those force field parameters which we observed from the above analysis.

We consider that, apart from our choice to test the P4, P5, P1, P2, and P3 parameter values from the original Martini polarizable force field, we have the option to test other parameters associated with alternative atom types defined in the Martini force field. These are the Na, N0, C5, Nd, and Nda atom types. We do not explicitly consider PMFs associated with using these parameters for the backbone. We offer a straightforward analysis suggesting that these other atom types and their parameters are not sufficient for the purpose of this paper. The major thesis of the present argument (and supported by our results so far) is that the residue–water interactions must be weakened and that residue–lipid interactions should be strengthened in order to lower the overall translocation PMFs as well as introduce interfacial minima. To estimate the appropriateness of the various Lennard-Jones type interactions between the residue backbone atom/bead and the other system components (water, lipid headgroup, lipid tail, etc.), we computed the bead–bead interaction curves as a function of separation (Figures S19–S23, Supporting Information). These curves show that the parameter sets B5 and Z5 we find to provide good agreement with the all-atom force field results represent limits of suitability. The alternative parameter combinations for residue–lipid headgroup attractive interactions are weaker than our final parameter set, Z5. Since the goal was to have a strengthened lipid headgroup–residue interaction, we see that the alternative parameters as present in the current Martini force field are not appropriate. Furthermore, we see from these curves that the water–residue and lipid–residue interactions are essentially equivalent in magnitude. This will result in little preferential free energetic stability of the residue in bulk solution versus the bilayer interior (or center). Thus, any effect introduced by the use of the alternative parameters existing in the Martini polarizable CG force field will be minimal, if at all.

Core Contribution to Total PMF

We consider contributions (total electrostatic and vdw) to the PMF from water and lipid molecules that are considered to be in the “core” region of the membrane. Core region lipids and water are defined in such a way that the lipid P or water O atoms for the AA model and PO4 or water beads for the CG model reside in a region ±13 Å from the center of bilayer.[23,29] The calculated membrane and water core contribution to the total PMF for three oligo-arginines as obtained from the AA and CG models are shown in Figure 9 and Figure S24 in the Supporting Information. The corresponding data are tabulated in Table 3. For the AA and polarizable CG models, we find that the core membrane contributions somewhat parallel the total membrane contribution. As the peptide moves into the bilayer, the core lipid contribution is stabilizing with all models, though the energy scales are not equivalent. In all cases, a minimum is predicted, after which the core lipid contribution becomes destabilizing. The position, incidentally, is strikingly similar for all peptides across the AA and polarizable force fields. In all three models, Arg3 experiences a substantially larger core lipid generated stability moving into the bilayer compared to Arg1 and Arg2. This is a consequence of the larger size of the peptide, allowing it to maintain interactions with the polar groups of the lipid, while not forcing the lipid molecules to be perturbed away from the canonical spatial distributions. In the case of the smaller Arg1 and Arg2, for core lipids to influence the peptide stability while in the bilayer, larger deformations are required (i.e., the lipid molecule must move further away to directly interact with the smaller peptides). The total core water contribution to the stability of the peptide is again very similar for the AA and two polarizable models. The AA result recapitulates several studies in the literature.[23,25,74,75] Core water is stabilizing, and core water stabilization increases with increasing peptide length (increasing charge). The AA model exhibits a barrier in the core water contribution profile at 10 Å; this barrier is not present in the CG models. The nonpolarizable CG model stands out uniquely, as the profiles are dramatically different from those of the other models. In this sense, one might consider that, at a high resolution, the polarizable CG models are a closer representation of the AA model results. The thermodynamic agreement between AA and CG models in the present study is again shown to be strikingly good.

Figure 9

PMF contributions from core-located (a–d) membrane and (e–h) water in the AA, nonpolarizable CG, polarizable CG, and CG Z5 models from left to right. Water or DMPC lipid molecules are considered as core water or core membrane, when the z distance of water central beads or lipid phosphate beads enters the region within ±13 Å of the center of the bilayers.

Table 3

PMF Contributions of Oligo-Arginine Translocation from Core-Located Water and Membrane (kcal/mol)

	core watera			core membranea
force field	Arg1	Arg2	Arg3	Arg1	Arg2	Arg3
AA	–52.60	–75.06	–81.44	24.39	–6.19	–123.45
CG-nonpol-std-P5	0.06	–0.28	–1.31	6.88	9.38	11.23
CG-pol-B5	–10.63	–16.07	–21.20	9.88	1.64	–15.23
CG-pol-Z5	–11.29	–17.02	–22.34	9.95	3.02	–12.34
CG-pol-std-P5	–12.11	–21.56	–30.78	9.26	2.10	–11.67
CG-pol-P4	–12.26	–21.69	–29.59	8.23	1.50	–14.00
CG-pol-P3	–11.90	–19.74	–29.05	8.21	1.32	–12.86
CG-pol-P2	–11.95	–18.96	–27.26	7.91	1.06	–12.73

Water or DMPC lipid molecules are considered as core water or core membrane, when the z distance of water central beads or lipid phosphate beads enters the region within ±13 Å of the center of the bilayers.

Further Implications: Structural Features of Membrane and Peptide along OP

Since the membrane undergoes structural changes to accommodate the oligopeptides, we next consider aspects related to such changes. Specifically, we compare the results for membrane deformation, orientational order parameter of lipid tail groups, and side chain orientation of oligopeptides among the AA, CG, and modified CG model systems.

Membrane Deformation

We consider the structural perturbations of the bilayer as solutes translocate. Structural perturbations of the bilayer leaflet at the local region around the translocating species have been widely reported by Li et al.,[76] Vorobyov et al.,[77] MacCallum et al.,[24] and Hu et al.[29,30] Such deformation helps to spatially modulate the interface between low and high electrostatic potential regions, thus stabilizing the charged oligo-arginine peptides in a favorable electrostatic environment.[77] Formation of water pores/defects in the membrane is facilitated by dramatic membrane deformations during charged oligo-arginine translocation.[24,30]

To visualize the deformation of the membrane surface, the z position of the heavy atoms of headgroups (including cholines, phosphates, and carbonyl groups for the AA model and choline (bead name NC3), phosphate (bead name PO4), and carbonyls (beads names GL1 and GL2) for the CG model) of the membrane were considered for all the snapshots of the window where the peptide resides at the bilayer center. At first, the lateral coordinates of all the atoms in the systems were shifted relative to the center of oligo-arginine and then the z-coordinates of the relevant atoms were averaged in bins spaced at a resolution of 3.1 Å × 3.1 Å for both the AA and CG systems in the x–y plane. Average deformation of membrane surfaces at the local region of translocating peptides was observed for both the AA and CG systems, shown in Figure 10. Detailed deformation surfaces of the top and bottom layers are shown in Figure S25 in the Supporting Information. Qualitatively, the figure shows that the membrane surfaces for the AA, CG nonpolarizable, CG polarizable, and modified CG polarizable model parameters have very similar surface features. Such results are important because the average surface behavior of the membrane is insensitive to the model. Moreover, the figure shows that the surfaces obtained from all CG models are smoother than those obtained from the AA model. This is reasonable because the resolution of the CG model is much lower than that of the AA model. The shape of the deformed membrane surface is characterized as a funnel with an average radius of around 20 Å and height of 12 Å of all oligopeptide systems. However, a compact structure of oligopeptide is expected to produce a narrow funnel shape deformed structure of membrane. We have characterized the compactness of oligopeptide by calculating the radius of gyration (Rg) for the AA, CG, and modified CG parameters, and the results are shown in Figure S26 in the Supporting Information. Consistent with expectation, the values of Rg are slightly increased from Arg1 to Arg3 for both the AA and CG models, which further affects the radii of the funnels for three oligopeptide systems.

Figure 10

3D surface of the average deformation of membrane. The surface is drawn from the z-coordinates of the headgroup heavy atoms averaged in bins spaced at a resolution of 3.1 Å × 3.1 Å for both AA and CG systems in the x–y plane.

The deformations of membrane surfaces are comparable for the AA, CG nonpolarizable, CG polarizable, and modified CG polarizable models. Despite this structural similarity, we do not discount the possibility of different deformation free energies contribute by each model to the overall translocation PMF. Such a deformation contribution has been estimated by using the membrane elasticity theory.[78] Note that the nature of the membrane deformation obtained from our study is similar to numerous earlier studies.[78] Using elasticity theory (Helfrich functional for a deformable sheet), Choe et al.[78] estimated the deformation free energy of the membrane which can be up to 5 kcal/mol. The deformation free energy of the membrane may have an important contribution to the total PMF of oligopeptide translocation.

Side Chain and Lipid Tail Groups Orientation

The orientation of lipid molecules is likely to deform the membrane surface to a greater extent. To explore the extent of such a deformation, we have calculated the local orientational order parameter for both the all atomistic and CG DMPC bilayer. Since the radius of the funnel-shape membrane deformation surface is up to 20 Å in the xy plane, we defined the local membrane as any lipid whose distance of the phosphate bead to the center of mass of the oligo-arginine is within 20 Å. The order parameter, Stail, is obtained by the following equation.The second-rank Legendre polynomial, P2(cos θ) = (1/2)(3 cos2 θ – 1), was computed for consecutive bonds, with θ being the angle between the direction of the bond and the bilayer normal. The values of P2 = 1, −0.5, and 0.0 correspond to perfect alignment, perfect antialignment, and a random orientation, respectively. Finally, we have calculated the average value of Stail over all the bonds present in lipid tails. Figure S27 in the Supporting Information shows the magnitude of calculated average order parameters, ⟨Stail⟩ as a function of OP for the all-atom, CG, and modified CG model systems. Figures S28, S29, and S30 in the Supporting Information show the average local order parameter value of each bond present in each lipid tail sn-1 and sn-2 in each simulation window of oligo-arginines. All of these models show very similar characteristic features and have a good agreement between atomistic and CG models. Further, our calculation shows that the bilayer becomes disordered when the oligopeptides approach the bilayer center. Moreover, we observed that the extent of disorder is increased form Arg1 to Arg3 for all the systems in the three models. As discussed earlier, such a difference in disorder enhances the possibility of membrane deformation and accelerates water penetration into the bilayer core as well.

The strong interactions between the guanidinium group of arginine and the membrane are expected to influence the peptide side chain orientation relative to the membrane headgroup region. To explore this, we define an angle between the membrane normal and the vector formed between the backbone and charged side chain charged for the AA model and beads for the CG model. We found that the change in angle distributions occurs as the oligo-arginines translocate into the bilayer. The change is independent of the model parameters. The results are shown in Figure 11 for each arginine residue separately. Further, we noticed that, when arginine is located at the bulk water region, the distance between oligo-arginine and the membrane is too far to affect the orientation of the side chain, so the angles are found to be distributed randomly. The values are near 40° at the center of the bilayer. Thus, as the oligo-arginine peptide approaches the membrane–water interface, the interaction of the side chain with polar/charged headgroup moieties helps to orient the side chain relative to the membrane and the angle is found to be larger than 90°. However, when oligo-arginines pass the headgroup region and move into the center of the membrane, the interactions reverse the orientation and direct the charged components of the side chain toward the positive z direction, and the angle becomes smaller than 90°.

Figure 11

Distribution of side chain orientational preference in (a–f) the AA model, (g–l) the nonpolarizable CG model, (m–r) the polarizable CG model, and (s–x) the CG Z5 model. Side chain orientation is represented by the average angle between the membrane normal and the guanidinium-backbone vector for each single side chain of the oligo-arginines.

The above discussion suggests that the average orientation behavior of the peptide side chain and lipid tail order for the AA, CG, and modified CG models follow the similar trends.

Peptide Translocation Kinetics

So far, we have addressed the free energy and structural features of oligo-arginine translocation into the membrane. In this section, we integrate free energetic information to estimate the translocation rate of a single oligo-arginine into the membrane using different models. To assess the rate of peptide translocation, we follow a diffusion based model as proposed by Hummer et al.[79,80] According to the model, the rate of translocation of a single peptide from umbrella sampling MD simulations can be obtained by combining the 1D free energy surface for peptide translocation with the Smoluchowski diffusion model. The resulting diffusion model is parametrized in terms of a 1D PMF and 1D position-dependent diffusion coefficient along a specified OP. The unidirectional rate of peptide translocation can be obtained from the following equationwhere ρ is the number density of peptide in water, S is the cross-sectional area in x and y dimensions, and G(z) and D(z) are the 1D PMF and position-dependent diffusion coefficients along z. To calculate the diffusion coefficient from the US window, we again follow their protocol.[80] According to the protocol, the local diffusion coefficient of the peptide under narrow harmonic potential can be estimated fromwhere the relaxation time, τ, is obtained by the following equationwhere var(z) and var(z̅) are the variance of z and variance of average z coordinate and Δt is the time interval between two data points. We used the Gromacs tool “g_analyze” to compute variances. Further details of the theory and method for the computation of the rate of such a translocation process can be found in Hummer’s work.[79]

The estimated rate constants and the time required to translocate a single oligo-arginine across the membrane obtained from three different models are tabulated in Table 4. The position-dependent diffusion constants D() calculated from the AA model and three CG umbrella sampling simulations (CG-nonpol, CG-pol-P5, CG-pol-Z5) are shown in the Supporting Information (Figure S31). Note that the z-position-dependent diffusion constants for the AA simulation were obtained from 41 MD simulations (4 ns trajectory each) by using a harmonic potential with a force constant 5 kcal/mol/Å2 range from 0.0 to 4.0 nm at a spacing of 0.1 nm along our chosen OP.[81] The kinetic data clearly shows that the estimated rate constants obtained from the AA model and both CG models are within the experimentally accessible range. The data further reveals that the results obtained for the AA, standard CG polarizable (bead type, P5), and nonpolarizable (bead type, P5) models are not of the same order of magnitude. We notice that such a difference is more prominent for Arg3. However, with the change of bead type from P5 to B5 for the CG polarizable model, the rate constants obtained from the AA and polarizable CG models are not also in the same order. However, both rate constants are in the same order of magnitude if the bead type of the polarizable CG model was changed from P5 to Z5. Thus, our current results explore that a short cationic peptide might translocate across a pure bilayer. However, such a translocation is practically feasible by changing the type of lipid and peptide type.[82−85] In fact, it has been shown recently using fluorescence experiment that the rate of translocation of the cationic amphiphilic peptide across the membranes of pure phospholipid giant vesicles is in the order of minutes scale.[86]

Table 4

Estimated Oligo-Arginine Translocation Rate Constant and Time Scale

	rate constant			translocation time
force field	Arg1	Arg2	Arg3	Arg1	Arg2	Arg3
AA	2.61 × 10–10/ps	1.89 × 10–11/ps	1.83 × 10–14/ps	3.8 ms	52.8 ms	54.7 s
CG-nonpol-std-P5	1.74 × 10–13/ps	5.55 × 10–16/ps	3.73 × 10–18/ps	5.7 s	30.1 min	74.4 h
CG-pol-std-P5	2.95 × 10–15/ps	1.85 × 10–19/ps	3.76 × 10–23/ps	339.5 s	1506.4 h	855 d
CG-pol-B5	2.19 × 10–13/ps	4.07 × 10–16/ps	2.56 × 10–18/ps	4.6 s	40.9 min	108.4 h
CG-pol-Z5	3.42 × 10–12/ps	5.78 × 10–14/ps	2.57 × 10–15/ps	292.0 ms	17.3 s	389.4 s

Summary and Conclusions

In this work, we have explored the free energetics and structural and kinetics features of single oligo-arginine translocation into the DMPC bilayer with the AA and polarizable and nonpolarizable CG Martini force fields. We used the ABF method for sampling our AA system, whereas the US method was used to sample the CG system. Moreover, the standard backbone bead parameter of peptide (P5) of the CG polarizable model was systematically changed to P4, P3, P2, B5, and Z5. We obtained the PMFs for the reversible transfer of a single oligo-arginine across the bilayer. We further decomposed each PMF, obtained from different force fields, into contributions from various system components including lipid, water, and ion. The effects of membrane deformation on the PMF have also been explored by calculating the PMF of the membrane and water core. Furthermore, combining the obtained 1D PMF for peptide translocation with the Smoluchowski diffusion model, we computed the corresponding kinetics of the process.

In the PMF analysis, we acknowledge the results obtained from AA simulation as a reference and then systematically varied the force field parameters ranging from the atomistic to modified polarizable CG model. The total PMF shows a significant barrier, on the order of 10 kcal/mol, for translocation of a single oligo-arginine across the DMPC bilayer. The free energy cost for transferring the peptides across the membrane increases from Arg1 to Arg3, and it shows a qualitatively similar trend for all of the force fields. A minimum in PMF has been found at the water–bilayer interface for the AA model, but the standard polarizable CG models failed to reproduce such a minimum. The presence of such a minimum indicates that the oligo-arginine binds with the bilayer at the interface region. Absence of an interfacial minimum in the PMF for the polarizable CG models may be related to the overly large dipole potential predicted by the polarizable model at the membrane–water interfacial region. Use of a different CG water model such as the BMW,[42] which produces the correct dipole potential at the interface, might lead to an interfacial minimum in computed PMFs; this represents important future work. Our analysis further revealed that the strong interaction between oligo-arginine and the membrane headgroup and a weak interaction between oligo-arginine and the water are sufficient to qualitatively reproduce such minima. However, 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 AA force fields. Some recent experimental studies on charged cell-penetrating peptides show that the binding of such highly charged species to bilayers is facilitated by some degree of anionic character.[83−85] Moreover, we found that the increase in vdw interaction between the backbone and lipid beads and the decrease in vdw interaction between all the peptide beads and the central bead of polarizable water provide a PMF which is close to the AA model (type: Z5). We further tested the Z5 backbone bead parameter on a longer peptide, Arg9, for validating it. The obtained PMF for Arg9 shows a minimum at the interface, and the barrier is reduced significantly. Such a result is in accordance with previous AA simulation studies which suggest that choosing the Z5 parameter for the backbone bead is more reasonable than the standard one (P5) to model the backbone bead of the peptide in the polarizable CG model. Decomposition of the total PMF into system components revealed that the membrane has a stabilizing contribution whereas ion and water have a destabilizing contribution. We found that the energy scales of each component for the three different models are different but the total PMFs are almost identical, which suggests that forces are well balanced among the three different components of the system for all of these models. Implications of the peptide size and force field parameters on sustainable structural perturbations of the bilayer have also been explored in our study. We found that the nature of deformations of arginine associated lipids is ostensibly quite different for the three oligo-arginines. Further, our analysis revealed that membrane deformation is larger for large peptides. The trend of membrane deformation free energy as obtained from membrane and water core contribution for the AA and polarizable CG models is qualitatively similar, but the nonpolarizable CG model shows different behavior. Further, estimation of the rate constant for the translocation of three oligo-arginines across the pure bilayer explores that a short cationic peptide might translocate through the membrane. Experimental studies have shown that the above translocation is practically feasible with change of lipid and peptide type.[83−85] In fact, with use of a fluorescence experiment, it has been shown that the rate of translocation of the cationic amphiphilic peptide across the membranes of pure phospholipid giant vesicles is in the order of minutes scale.[86,87] However, the free energetics of the translocation process is influenced by the structural and dynamical properties of the system’s components. Therefore, it would be interesting to explore how those properties are influenced by the different force fields. We are currently investigating this aspect in great detail.