Molecular dynamics studies of interactions between Arg9(nona-arginine) and a DOPC/DOPG(4:1) membrane
aa r X i v : . [ q - b i o . B M ] J un Molecular dynamics studies of interactions between Arg (nona-arginine) and aDOPC/DOPG(4:1) membrane Seungho Choe ∗ School of Undergraduate Studies, College of Transdisciplinary Studies,Daegu Gyeongbuk Institute of Science & Technology, 333 Techno Jungang-daero,Hyeonpung-eup, Dalseong-gun, Daegu 42988, South Korea (Dated: June 4, 2020)It has been known that the uptake mechanisms of cell-penetrating peptides(CPPs) depend onthe experimental conditions such as concentration of peptides, lipid composition, temperature,etc. In this study we investigate the temperature dependence of the penetration of Arg s intoa DOPC/DOPG(4:1) membrane using molecular dynamics(MD) simulations at two different tem-peratures, T = 310 K and T = 288 K. Although it is difficult to identify the temperature dependencebecause of having only one single simulation at each temperature and no evidence of translocationof Arg s across the membrane at both temperatures, our simulations suggest that followings arestrongly correlated with the penetration of Arg s: a number of water molecules coordinated byArg s, electrostatic energy between Arg s and the lipids molecules. We also present how Arg schange a bending rigidity of the membrane and how a collective behavior between Arg s enhancesthe penetration and the membrane bending. Our analyses can be applicable to any cell-penetratingpeptides(CPPs) to investigate their interactions with various membranes using MD simulations. I. INTRODUCTION
Cell penetrating peptides(CPPs), typically compris-ing 530 amino acids, can translocate across cell mem-branes. They have been extensively studied for a coupleof decades since they are capable of transporting vari-ous cargoes(e.g., proteins, peptides, DNAs, and drugs)into cells [1]. However, the uptake mechanisms of CPPsare still controversial. It is known that the translo-cation mechanisms of different families of CPPs arenot the same, and most CPPs can have more than asingle pathway depending on the experimental condi-tions [2]. There are two major uptake mechanisms,i.e., the endocytotic(energy-dependent) pathway and thedirect(energy-independent) translocation. The detaileduptake mechanisms can be found in Refs. [1–3].The factors that control the uptake pathways of CPPswere categorized into two groups [3]: First, the physicaland chemical properties, concentration of CPP and itscargo. Second, the properties of the membrane, its lipidand protein composition. It was also mentioned that theconcentration of CPPs is one of key factors, because itcan result in different uptake pathways. It is generallyaccepted that endocytosis occurs at low concentrationand it can switch to the direct translocation when theconcentration becomes higher. However, the effect of theconcentration on the uptake mechanisms can not be con-clusive. In the case of penetratin, it has been shown thatthe direct translocation occurs at low concentration whileendocytosis occurs at high concentration [4]. There arevarious models to explain the direct translocation, suchas the carpet-like model [5], pore formation [6], invertedmicelle formation [7], the membrane thinning model [8], ∗ [email protected] etc. Detailed mechanisms of those models are also givenin Refs. [2, 3].Among various CPPs, arginine-rich peptides seem veryinteresting. As mentioned above the switch between dif-ferent uptake mechanisms might depend on the concen-tration of the peptides. It has been shown in experi-ments that endocytosis is only dominant to a specificpeptide concentration and the direct translocation occursat levels higher than this threshold [9]. One of possiblescenarios for the insertion of arginine-rich peptides intothe lipid bilayer is the ability of the peptides to recruitnegatively charged phospholipid heads and thus decreasethe surface pressure of the membrane [3]. The chargeneutralization results in the membrane deformations andthis can lead to a decrease in a bending rigidity of themembrane. The rapid cytoplasmic entry has been ex-plained by the accumulation of the peptides at certainmembrane areas, called nucleation zones [10]. It was alsoreported that a rapid temperature decrease from 37 ◦ Cto 15 ◦ C induces an efficient entry of Arg (called R9 ornona-arginine in different literature) even at a low pep-tide concentration [11]. This is why two temperatures(T= 310 K and 288 K) were chosen in our simulations. Al-though the mechanism of temperature dependence shownin Ref. [11] is closely related to the calcium signalling,we focus primarily on the temperature dependence of in-teractions between Arg and the lipids molecules.It was shown by Eir´ıksd´ottir et al. [12] that severalCPPs can adopt a specific secondary structure in stabi-lizing their interactions with membranes. Among themArg exhibited a random coil in both water and mem-brane environments, and it means that Arg was unableto adopt any specific secondary structure in both envi-ronments. Eir´ıksd´ottir et al. [12] also carried out leakagetests to assess the consequence of interactions betweenArg and membranes on membrane stability, and foundthat Arg did not induce any leakage, suggesting no mem-brane perturbation and probably no insertion into thebilayer.It was suggested that the entry of Arg intoDOPG/DOPC-GUV(GUV: the giant unilamella vesicle)and DLPG/DTPC(2/8)-GUV occurs through transienthydrophilic prepores [13]. A prepore is an area of de-creased density of lipids due to the fluctuations of lo-cal lateral density. Prepores immediately close becauseof the large line tension. The interaction of CPPs withthese prepores may stabilize the prepores by decreasingtheir line tension. More detailed physical views on pre-pores and comparison of the line tension of several lipidmixtures are given in Ref. [13] and references therein.Herce et al. [14] found that the ionic permeability ofphospholipid bilayers and cell membranes is increased bythe presence of the Arg peptides. The authors showedthat Arg follows the same mechanism previously de-scribed for Tat [15] using molecular dynamics(MD) sim-ulations. The Arg peptides are attracted by the phos-phate groups of the phospholipids and these interactionsproduce large local distortions to the bilayer. This leadsto the formation of a toroidal pore. Note that they usedthe system containing a single phospholipid(DOPC) intheir MD simulations. However, the authors also showedin experiments that anionic phospholipids such as DOPGlipids are not strictly required because the peptides arestill able to interact and permeabilize purely zwitterionicDOPC phospholipid bilayers [14].MD simulations have been used to investigate func-tional properties of Arg and its interactions with manydifferent lipids, however, only a few simulations usedmixtures of lipids(e.g., a mixture of DOPC/DOPG orDOPC/DOPE) so far and the mechanism of transloca-tion of Arg and interactions with lipids are still not con-clusive. It was found that a DOPC/DOPG membraneappeared the better host for the translocation of KR9Cin MD simulations [16].In our study we investigate the penetration of Arg peptides into a model membrane(e.g., a mixture ofDOPC/DOPG(4:1)) using MD simulations at two dif-ferent temperatures(T = 310 K and T = 288 K). Ourfirst goal is to see if Arg can translocate across theDOPC/DOPG mixture by making a pore or a prepore,and the second goal is to investigate the temperaturedependence of interactions between Arg and the lipidsmixture. II. SIMULATION METHODS
All simulations were performed using the NAMD pack-age [17] and CHARMM36 force field [18]. CHARMM-GUI [19] was used to setup a DOPC/DOPG(4:1) mix-ture and TIP3P water molecules. The lipids mixtureconsists of 76 DOPC and 19 DOPG lipids in each layer.K and Cl ions were added to make a concentration of150 mM. One Arg peptide was made by the molefac-ture plugin in VMD [20] and it was duplicated to make four Arg peptides. The total number of Arg peptideswas the same as in Herce et al.’s [14] to make a similarconcentration of Arg s. All Arg peptides were initiallylocated in the upper solution and they were bound to theupper layer during the equilibration. The total numberof water molecules after making the system neutralizedwas 9643. The NPT simulations were performed at twodifferent temperature T = 288K and 310K in order toinvestigate the temperature dependence of penetrationof Arg s. Temperature and pressure were kept constantusing Langevin dynamics. An external electric field(0.05V/nm) was applied in the negative z direction(from theArg s peptides to the membrane) as suggested in previ-ous work [14, 21] to account for the transmembrane po-tential [22]. The particle-mesh Ewald(PME) algorithmwas used to compute the electric forces and the SHAKEalgorithm was used to allow a 2 fs time step during thewhole simulations. At the first stage, the system wasequilibrated at T = 310 K for 500 ns long. At the sec-ond stage, two systems were prepared: One was a systemwith the same temperature(T = 310 K) as in the equili-bration at the first stage. The other was a system withT = 288 K. When the temperature was down from T =310 K to T = 288 K in the second system, the systemwas quickly equilibrated at T = 288 K. Both systems(T= 310 K and 288 K) were ran for another 500 ns. Thefirst 200 ns data in both simulations were discarded forthe equilibration and only the last 300 ns data were usedfor the analyses. III. SIMULATION RESULTSA number of water molecules coordinated by Arg sis strongly correlated with the penetration depth We don’t see any translocation of Arg across theDOPC/DOPG bilayer during the simulations. However,it can penetrate into the lipids mixture along with aquite number of water molecules. We define a pene-tration depth as a distance in the z direction betweenthe center of mass of each Arg (taking only C α s for thesake of simplicity) and that of the upper leaflet. Here,a leaflet means a surface which consists of 95 phospho-rus(P) atoms in the upper or the lower layer.Fig. 1(a) shows the penetration depth of each Arg atT = 310 K. The solid lines denote the penetration depthof each Arg , and the dotted line represents an averageover those four values. The positive value means thatthe center of mass of Arg is located above the center ofmasses of phosphorus(P) atoms, and the negative meansbelow the phosphorus(P) atoms. One of Arg s(Arg3) islocated below the upper leaflet during the most of simu-lation time.Fig. 1(b) shows the same quantity as in Fig. 1(a),but at T = 288 K. It shows that Arg can penetrateslightly more into the lipids mixture at T = 288 K. One ofArg s(Arg3) is below the upper leaflet during the wholesimulation time, and it reaches even -5 ˚A in this figure.Note that each point in the figures corresponds to a blockaverage over 1 ns. We find that the largest penetrationdepth is -6.3 ˚A during the simulation at T = 288 K,and Fig. 2 shows a snapshot at that moment. C α s ofArg3 are shown as the same color code in Fig. 1(b).The green dots by the vdW representation of VMD cor-respond to phosphorus(P) atoms in the DOPC/DOPGmixture. The lipids molecules are shown as gray linesand water molecules are represented by vdW of VMD.Fig. 3 represents a deformation of the upper leaflet dueto the presence of Arg s and a relative distance betweenthe upper leaflet and each Arg using the snapshot inFig. 2. The lower leaflet is not shown in this figure forclarity. The blue color in the leaflet denotes a trough andthe yellow a crest. Some of C α s located below the leafletare not seen in this figure.A quite number of water molecules also penetrate intothe DOPC/DOPG mixture and this results in a mem-brane thinning and bending. Fig. 4 shows the totalnumber of water molecules located 15 ˚A apart from thecenter of mass of Arg3 at T = 288 K. We count onlywater molecules below the center of masses of phospho-rus(P) atoms. Fig. 4 shows a strong correlation betweenthe penetration depth of Arg and the number of watermolecules surrounding it. The number of water moleculesincreases as the penetration of Arg3 increases. This con-firms the previous results from MD simulations [23–25].At the current stage we can not conclude that the low-ering temperature will really help Arg s to make morepenetration into the lipids mixture. First of all, the dif-ference in average penetration depths at two tempera-tures is not that significant(about 1 or 2 ˚A in average).Second, we have only one single simulation at each tem-perature, and the total simulation time(300 ns) may notbe long enough to see the temperature dependence. How-ever, it turns out that the penetration depths are affectedby electrostatic energy between Arg s and the lipids mix-ture and we want to focus on how Arg s interact with thelipids mixture in the following analyses. Electrostatic energy between Arg s and theDOPC/DOPG mixture is important for thepenetration It is difficult to imagine how Arg s can pass throughthe DOPC/DOPG mixture because DOPG is a nega-tively charged lipid and therefore one can expect strongattractive interactions between Arg s and DOPG lipids.It was suggested that the electrostatic interaction be-tween positive charges of Arg s and negative charges oflipid molecules is crucial at the early stage of translo-cation [21]. We want to investigate this feature in oursimulations.We calculate electrostatic energy between each Arg and the DOPC/DOPG mixture at both temperaturesusing the NAMD Energy plugin in VMD. We want to Time (ns) -6-4-20246810 D i s t an c e ( Å ) Arg1Arg2Arg3Arg4Avg (a)
Time (ns) -6-4-20246810 D i s t an c e ( Å ) Arg1Arg2Arg3Arg4Avg (b)
FIG. 1. The distance between the center of mass of each Arg and the center of masses of phosphorus(P) atoms in the upperleaflet(in the z direction) at (a) T = 310 K (b) T = 288 K,respectively. The color code represents each Arg . Avg(greendotted line) means an average value over four Arg s. see if there is any correlation between the strength ofelectrostatic energy and the penetration depth. Fig.5(a) shows electrostatic energy between each Arg andthe DOPC/DOPG mixture at T = 310 K. Each colorcode(the same color scheme as in the previous figures)represents the interaction between one of Arg s and thelipids mixture. Avg denotes an average energy over allfour Arg s. Fig. 5(b) shows the same energy, but at T= 288 K. The average electrostatic energy at T = 288K is slightly lower(more stable) than that at T = 310 Kby 100 kcal/mol. This can be explained by comparingthe average penetration depths in Figs. 1(a) and 1(b).Arg can penetrate slightly more at T = 288 K becauseof a lower electrostatic energy. Note that in Fig. 5(b)one of the lowest energy regions(170-200 ns) correspondsto a region which shows the largest penetration depth inFig. 1(b). One can see a similar behavior in the caseof Arg1(blue) at about 200 ns in Figs. 1(b) and 5(b). FIG. 2. A snapshot, which shows the largest penetrationdepth of Arg3(yellow), at T = 288 K. Only Arg3(yellow), thephosphorus(P) atoms(green dots), the lipids molecules(greylines), and water molecules in vdW of VMD are shownFIG. 3. Membrane curvature and a relative position of eachArg (showing C α s only) from the snapshot in Fig. 2. Repre-sentation in 3D(left) and in the x-y plane(right), respectively. In both Figs. 5(a) and 5(b) electrostatic energy betweenArg3(yellow) and the lipids mixture increases(less stable)at the end of simulation because the penetration depth ofArg3 decreases as shown in Figs. 1(a) and 1(b), respec-tively. Our simulations clearly suggest that electrostaticenergy plays a key role to control the penetration depthof Arg .Then, the next question will be whether Arg can makea pore or a prepore in the presence of strong electrostaticenergy between Arg s and the lipids mixture. At thisstage it is not clear how this can happen, but we thinkthat Fig. 2 can give us a clue. In Fig. 2 one can observea quite number of water molecules coordinated by Arg3when it makes the largest penetration through the lipidsmixture at T = 288 K. It was known by molecular dynam-ics simulations that water molecules can enhance a mem-brane bending and disorder the lipid molecules [23–25].
175 180 185 190 195
Time(ns) -7-6-5-4-3-2-10 D i s t an c e ( Å ) N u m be r o f W a t e r M o l e c u l e s FIG. 4. The distance between the center of mass of Arg3 andthe center of masses of 95 phosphorus(P) atoms in the upperleaflet(solid blue line, the left axis) and a number of watermolecules within 15 ˚A of Arg3(dotted red line, the right axis)at T = 288 K. We count only water molecules below the centerof masses of P atoms.
We conjecture that the concentration of Arg s will be oneof key factors to coordinate more water molecules. Morewater molecules will accelerate the membrane bendingand thinning, and thus there may be more interactionsbetween Arg s and the lipids molecules in the lower layer.However, the concentration of Arg s in our simulationsseems not large enough to see this behavior. In the Ap-pendix we present the membrane thinning due to thepenetration of Arg3. It is obvious that the increase ofpenetration of Arg s results in the membrane thinning.The detailed method to obtain the distance between thelayers is explained in the Appendix. More on the mem-brane bending is given in the below. Bending rigidity of the membrane is decreased dueto the presence of Arg s It was shown that there is a correlation between CPPhydropathy and the penetration of water molecules in thelipid bilayer and the presence of water molecules locallydestabilizes the lipid order and reduces a bending mod-ulus of the membrane [26]. The reduction in membranebending rigidity has also been reported in other literature[3, 16].We compare bending moduli of the DOPC/DOPGmixture at two different temperatures(T = 288 K and 310K). We compute the bending moduli using the methoddescribed in Ref. [26–28]. A similar approach to calcu-late the bending moduli of various lipid molecules wassuggested by Levine et al. [29] The authors also triedKhelashvili et al.’s method [28] to analyze their data andfound that their bending moduli are consistent with Khe-lashvili et al.’s. We use Khelashvili et al.’s method forobtaining the bending modulus of the membrane and this
Time (ns) -1800-1600-1400-1200-1000-800-600-400-200 E ( kc a l / m o l ) Arg1Arg2Arg3Arg4Avg (a)
Time (ns) -1800-1600-1400-1200-1000-800-600-400-200 E ( kc a l / m o l ) Arg1Arg2Arg3Arg4Avg (b)
FIG. 5. Electrostatic energy between Arg s and theDOPC/DOPG lipids at (a) T = 310 K (b) T = 288 K. Thesame color code in Fig. 1 is used. Avg(green dotted line)means an average value over four Arg s. will be enough to see the temperature dependence of thebending moduli in our simulations.It was shown that one can calculate the splay modulus χ by calculating splay angles( α ) between lipids molecules.We use the same definition of a splay angle in Ref. [26–28](see more in the Appendix). There are two criteriato analyze the splay angles as explained in the previouswork [28, 29]: 1) Only lipid pairs separated by less than10 ˚A are considered. 2) At least one director shouldbe oriented less than or equal to 10 degrees from thebilayer normal. In this approach a normalized probabilitydistribution is used to calculate the potential of meanforce(PMF). PMF is defined byPMF( α ) = − k B T ln P ( α ) P ( α ) (1)where P ( α ) = sin( α ) is the probability distribution ofa hypothetical non-interacting particle system, k B theBoltzmann factor and T the temperature. The splay modulus χ is obtained by a quadratic fit of the PMFdata, and it corresponds to the bending modulus of amonolayer K m . The presence of Arg , e.g, in the up-per layer, also affects the bending modulus of the lowerlayer. Therefore, we compute the bending modulus K c of a bilayer by averaging two moduli from each layer.Fig. 6(a) shows normalized probability distributionsof α , P( α ), at T = 310 K. The blue solid line is a dis-tribution from the upper layer, while the red dotted linecomes from the lower layer. Two lines are almost over-lapped with each other and the same behavior is seen inprobability distributions at T = 288 K(see Fig. 6(b)).From the second order polynomial fitting of Eq.(1) asdescribed in Ref. [26, 28, 29], we obtain K c (310 K) = 7.6 × − J and K c (288 K) = 6.8 × − J as bendingmoduli of the DOPC/DOPG mixture(see Figs. 11 and12 for the quadratic fits). The bending modulus at T =288 K is slightly smaller, but the difference is minimal.Our values( 17 ∼ k B T) are slightly smaller than anexperimental value, 20 k B T ± k B T [30]. This is dueto the penetration of water molecules along with Arg sand water molecules play a role to make softening of thelipids mixture. A collective behavior between Arg s enhances thepenetration and the membrane bending Li et al. [31] found that it is difficult for a single pep-tide to penetrate through the membrane while multiplepeptides can translocate across membranes by making apore. Hu et al. [32] also showed that the water pore for-mation is facilitated by the aggregated charge peptides,and the authors claimed that the pore may be a result ofthe cooperative effect from global and local disturbanceof membrane by random and aggregated peptides.Although we don’t see any translocation of Arg s dur-ing our simulations, our simulation results are consistentwith above findings that the collective behavior betweenArg s significantly changes a local membrane curvatureand results in more penetration into the lipids mixture.In Fig. 3 one can see that two Arg s(red and yellow)come closer and make a significant change in the mem-brane bending and thinning(see also the Appendix andFig. 10 for the membrane thinning).First, we investigate a relation between the penetra-tion depth of Arg3(yellow) and the distance between thecenter of mass of Arg2(red) and that of Arg3(yellow).In Fig. 7(a) the solid blue line denotes the penetrationdepth of Arg3 and the dotted red line shows the distancebetween Arg2 and Arg3. We show only data from 175ns to 195 ns, which include the moment of the largestpenetration depth of Arg3. Whenever Arg3 is trying toreach one of the largest penetration depths(e.g., about-6˚A at 184 ns or at 193 ns), the distance between Arg2and Arg3 is going to increase rapidly.We have another plot which shows a comparison be-tween the penetration depth of Arg3 and electrostatic (degree) P () Upper LayerLower Layer (a) (degree) P () Upper LayerLower Layer (b)
FIG. 6. Normalized probability distributions of α at (a) T =310 K (b) T = 288 K. The blue solid line denotes a distributionfrom the upper layer while the red dotted line from the lowerlayer. energy between Arg2 and Arg3. We expect that thefinding in Fig. 7(a) regarding the distance between Arg2and Arg3 could be related to electrostatic energy betweenArg2 and Arg3, and this analysis will give an idea whythe distance between Arg2 and Arg3 is slightly reducedwhen Arg3 makes the maximum penetration depth at184.68 ns in Fig. 7(a). In Fig. 7(b) the solid blue lineis the same as in Fig. 7(a)(the penetration depth) andthe dotted red line is electrostatic energy between Arg2and Arg3. Whenever Arg3 reaches the maximum pene-tration, electrostatic energy between Arg2 and Arg3 be-comes the minimum(or less repulsive).The induction of curvature by association of proteinsto the membrane has been well known for the so-calledcurvature sensing proteins [33] For example, it was knownthat if one protein induces a given membrane curvature,it attracts more proteins that favor a similar curvature[34]. It was also shown that curvature-inducing pro-teins experience attractive short-range interactions due to a membrane curvature, and this curvature-mediatedinteraction leads to protein aggregation [35]. In orderto estimate membrane-mediated interactions(attractiveor repulsive potential) between Arg s we need to con-sider long-range interactions such as entropic membranefluctuations,[36, 37]however, we want to focus primarilyon a cooperative effect between Arg s which affects thepenetration depth and the membrane bending. Our elec-trostatic energy calculations seem to be enough to showthis feature.We want to briefly mention the usefulness of mem-brane curvatures such as the mean curvature( H ) or theGaussian curvature( K ) in our simulations [38]. Each cur-vature is defined by H = 12 ( c + c ) (2) K = c c where c and c are the principal curvatures at a pointon the membrane. Fig. 8 represents the mean curvatureand the Gaussian curvature of the membrane in the pres-ence of Arg s when Arg shows the largest penetrationat T = 288 K(see also Fig. 3). The blue color in the fig-ure means a negative curvature while the yellow denotesa positive curvature. As we expect the mean curvaturenear Arg2(red) and Arg3(yellow) is negative which meansa concave deformation. We have checked that the Gaus-sian curvature near Arg3 is always positive during thesimulation at T = 288 K, which means no evidence ofthe pore formation [38].We expect that both the concentration of Arg s andthe collective behavior are important in translocation.More simulation studies with different concentration ofArg s may be necessary to gain more insights into thecollective behavior. IV. DISCUSSION
Our simulations show that Arg can penetrate into theDOPC/DOPG mixture, but it can not pass through themixture both at T = 310 K and at T = 288 K. It doesn’tshow any pore or prepore as suggested in Refs. [13, 14].Although we don’t see any translocation of Arg , our sim-ulations suggest that electrostatic energy between Arg and the lipids mixture plays a role in determining thepenetration depth. A few remarks regarding our resultsare given below.The role of the Arg side chain was identified in the pre-vious MD simulations [14]. The Arg side chains initiallybind to the phospholipid phosphate groups, producingstrong distortions of the bilayer which lead to the for-mation of a pore. Our simulation shows that more than70 water molecules also penetrate into the lipids mixturewhen Arg reaches the maximum penetration depth(seeFig. 4), however, Arg s in our simulations don’t makesuch strong distortions in the DOPC/DOPG mixture
175 180 185 190 195
Time(ns) -7-6-5-4-3-2-10 D i s t an c e A r g3 - P ( Å ) D i s t an c e A r g2 - A r g3 ( Å ) (a)
175 180 185 190 195
Time(ns) -7-6-5-4-3-2-10 D i s t an c e A r g3 - P ( Å ) E ( kc a l / m o l ) (b) FIG. 7. The penetration depth of Arg3(solid blue line, leftaxis) at T = 288 K and (a) the distance between Arg2 andArg3(dotted red line, right axis) (b) electrostatic energy be-tween Arg2 and Arg3(dotted red line, right axis), respectively. which can lead to a pore. We expect that the transloca-tion may be difficult at low concentration of Arg s dueto the attractive potential energy between Arg s(positivecharge) and DOPG lipids(negative charge). We find thatall Arg s were strongly bound to the DOPG lipids dur-ing the simulations. We conjecture that the ability ofthe translocation can be highly dependent on the con-centration of Arg s. To our knowledge, it has not yetknown how positively charged Arg s pass through theDOPC/DOPG mixture within all-atom MD simulations.One of crucial reasons will be a simulation timescale. Itwas claimed that spontaneous translocation is not ex-pected within timescales accessible to MD simulations,unless nonphysical restraints are applied to acceleratethe underlying process [39]. In that case, a path sam-pling approach(e.g., the weighted ensemble method [40])might be one of solutions to reveal the underlying mech-anism of the translocation of Arg s.It was known that charged surfaces are known to be FIG. 8. The mean curvature(left) and the Gaussian curva-ture(right) of the membrane in the presence of Arg s(showingonly C α s) when Arg3(yellow) makes the largest penetration atT = 288 K. The blue color in each curvature means a negativevalue while the yellow a positive. more rigid than neutral ones, therefore, the decrease inmembrane rigidity could be explained by considering thatthe arginine-rich peptide partially neutralizes these mem-branes [16]. Crosio et al. [16] claimed based on their ex-periments in the DOPC/DOPG(1:1) environment thatthe recruiting of the anionic lipid in the regions closeto the peptide would lead to an enrichment of DOPCin other regions, subsequently decreasing the observedbending rigidity. The incorporation of water coupledwith peptide translocation and local membrane thinningcan be another explanation for the reduced bending rigid-ity of the membrane [26]. Our simulations suggest thatboth the partial neutralization and the membrane thin-ning affect a softening of the DOPC/DOPG mixture. Aswe mentioned above, Arg s are strongly bound to DOPGlipids during the simulations and this results in the par-tial neutralization of the membrane. We also showedthe membrane thinning due to the penetration of watermolecules in the previous section. It will be interesting toinvestigate which factor(the neutralization or the mem-brane thinning) can contribute more on the softening ofthe membrane. This can be done using different concen-tration of Arg s as well as using different DOPC/DOPGratios.Although the total number of Arg s and the systemsize may not be large enough to clearly show the col-lective behavior in our simulations, our results suggestthat it can affect the membrane bending and thinningand thus the penetration of Arg s. Movie S1 showshow each Arg makes membrane deformations duringthe time frame between 175 ns and 195 ns which in-cludes the moment that Arg3 shows the largest pene-tration depth at about 185 ns. One can see in thismovie that there are two main local minima made byArg1(blue) and Arg3(yellow). Arg2(red) is already gotcloser to Arg3(yellow) and it is trapped in a local mini-mum during this time frame. One interesting finding isthat sometimes two local minima are connected to eachother, and this results in a drastic change in the mem-brane bending(see, e.g., at 188.6 ns, 190.4 ns, and 191.6ns in the movie). One can find an increased penetra-tion depth and a drastic decrease of electrostatic energyof Arg1(blue) at about 190 ∼
200 ns in Figs. 1(b) and5(b), respectively. We compute electrostatic energy be-tween Arg1(blue), Arg2(red), and Arg3(yellow) as we didpreviously for the interaction between Arg2 and Arg3.Fig. 9 represents these calculations with the penetra-tion depth of Arg1 for comparison. The blue solid lineis the penetration depth of Arg1 and the red dotted lineis the sum of two electrostatic energy: One is electro-static energy between Arg1 and Arg2, and the other oneis that between Arg1 and Arg3. There is an interest-ing energy drop at about 185 ns, and this time framecorresponds to the moment that Arg3 makes the maxi-mum penetration depth. One can see that the penetra-tion increases at the time frames mentioned above, i.e.,at 188.6 ns, 190.4 ns, and 191.6 ns, and electrostatic en-ergy becomes minimum (or less repulsive) at these timeframes. It shows that there exits the collective behav-ior between Arg s which enhances the disorder of themembrane. Further simulations with more detailed cur-vature analyses [33, 38] will be helpful to identify the fac-tors(e.g., concentration of Arg s, a distribution of DOPGlipids in the DOPC/DOPG mixture, etc.) which controlthe collective behavior between Arg s and their effectson the membrane bending.
175 180 185 190 195
Time(ns) -3-2-101234 D i s t an c e A r g1 - P ( Å ) -20020406080100120 E ( kc a l / m o l ) FIG. 9. The penetration depth of Arg1(solid blue line,left axis) at T = 288 K and the sum of electrostatic en-ergy between Arg1 and Arg2, and that between Arg1 andArg3(dotted red line, right axis)
V. CONCLUSION
The uptake mechanisms of CPPs have been extensivelystudied in both experiments and computer simulations. However, they are not still conclusive. The difficultycomes from the fact that the uptake pathway is not thesame for different families of CPPs and depends on theexperimental conditions [2]. There is no doubt that anMD simulation is still one of valuable tools to reveal func-tional and mechanical properties of different families ofCPPs and their interactions with various types of lipidmembranes.Our simulations suggest that the electrostatic interac-tion between Arg s and the lipids molecules is importantat the early stage of the penetration, and water moleculescoordinated by Arg s make a membrane thinning whichresults in a softening of the bending rigidity of the mem-brane. The collective behavior between Arg s also playsa role to help the penetration of Arg s into the mem-brane. Although we don’t see any direct translocation ofArg across the DOPC/DOPG bilayer in our simulations,our analyses will be helpful for understanding functionalproperties of other CPPs and their interactions with var-ious membranes using MD simulations. Appendix A: Membrane thinning due to thepenetration of Arg We investigate a relation between the penetrationdepth of Arg3(one of Arg9 peptides in our simulations)and the membrane thickness at T = 288 K(see Fig.10). When the penetration depth of Arg3 increases,the membrane thickness near Arg3 is expected to be de-creased. The first panel represents the penetration depthof Arg3(the distance between the center of mass of Arg3and the center of masses of phosphorus(P) atoms in theupper leaflet, where a leaflet is a surface which consistsof 95 phosphorus(P) atoms in each layer) between 150 nsand 300 ns. The negative value means that the center ofmass of Arg3 is below than the center of masses of phos-phorus(P) atoms, and the positive means the reverse.The second panel shows time variations of the membranethickness due to the penetration of Arg3. The smallestthickness is about 27 ˚A, which is almost 10 ˚A smallerthan that of a membrane without any interactions withArg9s(a rest state). We compute the membrane thick-ness as follows: First, we find the three nearest phospho-rus(P) atoms in the upper leaflet from the center of massof Arg3 and calculate an average z position of these threeatoms. Next, we do the same thing in the lower leaflet,namely, we obtain another average z position of threenearest phosphorus(P) atoms in the lower leaflet. Then,we define a membrane thickness as a distance betweenthose two z positions. One can see that two figures showa very similar trend. This is what we exactly expected.We choose Arg3 because it shows the largest penetrationdepth at T = 288 K. The other Arg9s can have differentpatterns in the membrane thickness.
150 200 250 300
Time(ns) -8-6-4-2024 D i s t an c e A r g3 - P ( Å )
150 200 250 300
Time(ns) M e m b r ane T h i ck ne ss ( Å ) FIG. 10. The penetration depth of Arg3(the upper panel)and the membrane thickness(the lower panel) at T = 288 K.
Appendix B: Bending modulus of the DOPC/DOPGmixture
One of mechanical properties of a membrane is itsbending modulus. It was known that the penetrationof a quite amount of water molecules plays a role to de-crease a bending rigidity of a membrane. We investigatethe bending modulus of the DOPC/DOPG mixture inthe presence of Arg9s both at T = 310 K and at T = 288K. We use Khelashvili et al.s approach [28] to get thebending modulus from MD simulations. One can com-pute a bending modulus from a probability distributionof splay angles. Here is a general method to obtain splayangles [28]. First, local lipid director vectors ~v is defined,where ~v is the vector that connects the midpoint betweenthe phosphate and backbone C2 atoms to the geometric center of mass of the three terminal carbons on the twolipid chains. The splay modulus for a particular pairof lipid molecules is calculated by defining a splay angle( α ) between ~v vectors and constructing probability dis-tributions P( α ) (see the main text and Ref. [28]). Fig.11 represents the potential of mean force (PMF) pro-file PMF( α )= - k B T ln (P( α )/sin( α )) and a quadraticfit at [10, 30] degrees at T = 310 K, where sin ( α ) is theprobability distribution of a hypothetical non-interactingparticle system. A figure in the left represents the PMFprofile(the blue solid line) from the upper layer and aquadratic fit(the blue dotted line). A figure in the rightshows the same plots from the lower layer.In order to obtain a quadratic fit of the PMF profilewe use the Curve Fitting Tool of MATLAB. We obtainthe following bending modulus at T = 310 K: 17.408 k B T(the upper layer) and 17.990 k B T(the lower layer). Wetake an average of these two values, 17.699 k B T , as thebending modulus of the DOPC/DOPG mixture at T =310 K. Fig. 12 is the same as in Fig. 11, but at T =288 K. As for the probability distributions of the splayangles, there are no significant changes at T = 288 K.Using the same approach, we get the following bendingmodulus: 17.428 k B T(the upper layer) and 16.600 k B T(the lower layer). Again, we take an average of thesetwo values, 17.014 k B T, as the bending modulus of theDOPC/DOPG mixture at T = 288 K.
ACKNOWLEDGMENTS
The author thanks Hanna Salman and Xiao-lun Wu forvaluable comments on the early stage of this work, andhe also thanks the University of Pittsburgh for supportduring his stay. This work was supported in part bythe start-up fund from the Daegu Gyeongbuk Instituteof Science & Technology(DGIST). [1] G. Guidotti, L. Brambilla, and D. Rossi, Trends in Phar-macological Sciences , 406 (2017).[2] F. Madani, S. Lindberg, U. Langel, S. Futaki, andA. Gr¨aslund, J. Biophys. , 1 (2011).[3] I. Ruseska and A. Zimmer, Beilstein J. of Nanotechnol. , 101 (2020).[4] J. Pae, P. S¨a¨alik, L. Liivamagi, D. Lubebets, P. Aruku-usk, U. Langel, and M. Pooga, J. Controlled Release ,103 (2014).[5] Y. Pouny, D. Rapaport, A. Mor, P. Nicolas, and Y. Shai,Biochemstry , 12416 (1992).[6] K. Matsuzaki, S. Yoneyama, O. Murase, and K. Miya-jima, Biochemistry , 8450 (1996).[7] D. Derossi, S. Calvet, A. Trembleau, A. Brunissen,G. Chassaing, and A. Prochiantz, J. Biolo. Chem. ,18188 (1996).[8] M. T. Lee, W. C. Hung, F. Y. Chen, and H. W. Huang,Biophys. J. , 4006 (2005). [9] M. M. Fretz, N. A. Penning, S. Al-Taei, S. Futaki,T. Takeuchi, I. Nakase, G. Storm, and A. T. Jones,Biochem. J. , 335 (2007).[10] R. Wallbrecher, T. Ackels, R. A. Olea, M. J. Klein,L. Caillon, J. Schiller, P. H. Bov´ee-Geurts, T. H. vanKuppevelt, A. S. Ulrich, M. Spehr, M. J. W. Adjobo-Hermans, and R. Brock, J. Controlled Release , 68(2017).[11] K. Melikov, A. Hara, K. Yamoah, E. Zaitseva, E. Zaitsev,and L. V. Chernomordik, Biochem. J. , 221 (2015).[12] E. Eir´ıksd´ottir, K. Konate, U. Langel, G. Divita, andS. Deshayes, Biochimica Biophysica Acta , 1119(2010).[13] M. Z. Islam, S. Sharmin, M. Moniruzzaman, and M. Ya-mazaki, Appl. Microbio. Biotech. , 3879 (2018).[14] H. D. Herce, A. E. Garcia, J. Litt, R. S. Kane, P. Martin,N. Enrique, A. Rebolledo, and V. Milesi, Biophys. J. ,1917 (2009). (radian) P M F () ( k B T ) upper layerquadratic fit (radian) P M F () ( k B T ) lower layerquadratic fit FIG. 11. PMF profiles and quadratic fits for both the upper and the lower layers at T = 310 K. The blue solid line and theblue dotted lines in the left figure are a profile from the upper layer and a quadratic fit at [10, 30] degrees, respectively. Theright figure is the same, but from the lower layer. (radian) P M F () ( k B T ) upper layerquadratic fit (radian) P M F () ( k B T ) lower layerquadratic fit FIG. 12. The same figure as in Fig. 11, but at T = 288 K.[15] H. D. Herce and A. E. Garcia, Proc. Natl. Acad. Sci.USA , 20805 (2007).[16] M. A. Crosio, M. A. Via, C. I. C´amara, A. Mangiarotti,M. G. D. P´opolo, and N. Wilke, Biomolecules , 625(2019).[17] J. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhor-shid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, andK. Schulten, J. Comput. Chem. , 1781 (2005).[18] B. R. Brooks, C. L. B. III, A. D. M. Jr., L. Nilsson, R. J.Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels,A. C. S. Boresch, L. Caves, Q. Cui, A. R. Dinner, M. Feig,S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera,T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pas-tor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M.Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York,and M. Karplus, J. Comput. Chem. , 1545 (2009).[19] S. Jo, T. Kim, V. G. Iyer, and W. Im, J. Comput. Chem. , 1859 (2008).[20] W. Humphrey, A. Dalke, and K. Schulten, J. Molec.Graphics , 33 (1996).[21] A. Walrant, A. Vogel, I. Correia, O. Lequin, B. E. S. Olausson, B. Desbat, S. Sagan, and I. D. Alves, Biochim-ica Biophysica Acta , 1755 (2012).[22] B. Roux, Biophys. J. , 4205 (2008).[23] J. A. Freites, D. J. Tobias, G. von Heijne, and S. H.White, Proc. Natl. Acad. Sci. USA , 15059 (2005).[24] S. Dorairaj and T. W. Allen, Proc. Natl. Acad. Sci. USA , 4943 (2007).[25] J. L. MacCallum, W. F. D. Bennett, and D. P. Tieleman,J. Gen. Physiol. , 371 (2007).[26] G. Grasso, S. Muscat, M. Rebella, U. Morbiducci, A. Au-denino, A. Danani, and M. A. Deriu, J. Biomech. , 137(2018).[27] G. Khelashvili, G. Pabst, and D. Harries, J. Phys. Chem.B , 7524 (2010).[28] G. Khelashvili, B. Kollmitzer, P. Heftberger, G. Pabst,and D. Harries, J. Chem. Theory Comput. , 3866 (2013).[29] Z. A. Levine, R. M. Venable, M. C. Watson, M. G.Lerner, J.-E. Shea, R. W. Pastor, and F. L. H. Brown, J.Am. Chem. Soc. , 13582 (2014).[30] J. Steink¨uhler, P. D. Tillieux, R. L. Knorr, R. Lipowsky,and R. Dimova, Sci. Rep. , 11838 (2018). [31] Z. Li, H. Ding, and Y. Ma, Soft Matter , 1281 (2013).[32] J. Hu, W. Tian, and Y. Ma, Macromol. Theory Simul. , 399 (2015).[33] O. Maniti, H. Piao, and J. Ayala-Sanmartin, Int. J.Biochem & Cell Biolo. , 73 (2014).[34] P. Sens, L. Johannes, and P. Bassereau, Current Opinionin Cell Biology , 476 (2008).[35] B. J. Reynwar, G. Illya, V. A. Harmandaris, M. M.M¨uller, K. Kremer, and M. Deserno, Nature , 461(2007).[36] M. Goulian, R. Bruinsma, and P. Pincus, Europhys. Lett. , 145150 (1993). [37] H. Agrawal, L. Liu, and P. Sharma, Soft Matter , 8907(2016).[38] N. Schmidt, A. Mishra, G. H. Lai, and G. C. L. Wong,FEBS Letters , 1806 (2010).[39] S. Yesylevskyy, S.-J. Marrink, and A. Mark, Biophys. J. , 40 (2009).[40] M. C. Zwier, J. L. Adelman, J. W. Kaus, A. J. Pratt,K. F. Wong, N. B. Rego, E. Suarez, S. Lettieri, D. W.Wang, M. Grabe, D. M. Zuckerman, and L. T. Chong, J.Chem. Theory Comput.11