Structure and dynamics of Penetratin’s association and translocation to a lipid bilayer

Homeoproteins are proteins, first discovered in Drosophila, that bind to specific sites of DNA, and are functional to transcription. The domain involved in this binding (homeodomain) consists of a sequence of 60 amino acids forming three alpha helices.¹ Joliot et al.,² while studying the role of homeoproteins in the development of the nervous system, found, to their surprise, that the homeodomain of the Antennapedia protein of Drosophila was able to translocate over cell membranes. In pursuance of understanding the mechanism of its internalization, the authors altered the sequence of the homeodomain and found that the membrane translocation takes place at the third helix. This finding derived in the development of a synthetic cationic 16 aminoacid long peptide, named Penetratin (residue sequence: RQIKIWFQNRRMKWKK-NH2).³ Small cationic peptides of this kind are important not only because they can penetrate cells but, additionally, some of them have been shown to possess antibiotic,⁴ and/or antimicrobial properties.^5,6 The intriguing question since this discovery is: how can this highly cationic peptide (+7 charge) overcome the energetic barrier necessary to translocate into the lipid bilayer? Moreover, once inserted, how do the lipids in the membrane accommodate in presence of the strong electrostatic field generated by these charges? Another fundamental question is, what is the role of the peptide’s secondary structure in the insertion process?

The structure of Penetratin was first elucidated using H-NMR and CD experiments in extracelullar and in membrane mimetic environments.⁷ In a TFE/water solution, Penetratin adopts an incipient helical structure, with the 4-12 residues’ region maintaining $α −$ and 310 helices, and the terminals, β-turn structures. It was also found that the helical content lowers as the TFE concentration is decreased. Under different conditions, in a plasmon waveguide resonance experiment with different types of lipid bilayers, it was reported that Penetratin has an antiparallel β-sheet hairpin structure, with the charged amino acids pointing outwards from both sides of the bent molecule, facilitating an edge-on binding of the peptide to the charged surface of the lipid membrane.⁸ Supporting a β-sheet structure, the application of polarized infrared spectroscopy⁹ finds antiparallel β-sheet structure with a small fraction of turns, on three different bilayers. Superficial binding to the bilayers was also reported, with the charged lysine and arginine side chains forming H-bonds to the phosphate oxygens of the lipids.

Helices and β-sheets were not the only conformations reported for Penetratin; random coil structures were also found in living melanoma cells, using Raman microscopy.¹¹ 

On the computational side, structural studies reported a mix of bent conformations, 310 helices and β-turns, with both C- and N- termini highly flexible. A detailed simulation of Penetratin‘s behavior in charged and zwitterionic model bilayers shedded some light on the molecular aspects of this conformational diversity.¹² The authors adduce that the peptide adapts to the water-lipid interface by means of its specific structural plasticity. Within this model, the peptide’s conformation is determined by the side chain interactions of its aromatic and basic residues, which prepare it for insertion. Pourmousa et al.¹³ reported Penetratin’s structure as α-helical in the region between residues 2 and 6, and a variety of others, including π-helix and turns, in the rest of the peptide. The authors proposed that the peptide may adopt each structure depending on the initial configuration when it makes the first contact with the bilayer’s surface, and that all the different conformations found in the simulation are part of the penetration process, in such a way that, by switching between them, the peptide searches for the optimal shape and orientation that will allow it to penetrate the membrane. The authors formulated this model based on their findings for Transportan.¹⁴ 

In this work, we use several explicit Molecular Dynamics (MD) simulations in order to study the physical mechanisms of Penetratin association to a dioleoyl-phosphatidylcholine (DOPC) membrane with one of the latest developed force field. The conformational preferences of Penetratin close to the bilayer are fully investigated using a Markov State Model (MSM). A robust free energy analysis, with Umbrella Sampling (US), is used to study the translocation mechanisms of Penetratin, and the results are compared with findings from Yesylevskyy et al.³¹ 

Initial coordinates of Penetratin were taken from the Protein Data Bank, reference code 1OMQ.¹⁵ This positively charged peptide was initially placed in a TIP3P²⁹ water box, containing 4720 molecules, and the system was neutralized adding 7 Cl^- counter ions. The Amber force field ff14SB was used to describe the peptide.¹⁶ The energy of the system was minimized using the steepest descent method first, followed by a conjugate gradient method. The system was then heated to 300K and pressure was equilibrated at 1 atm. After equilibration, the following protocol was used to select the initial configuration for Penetratin in membrane: An accelerated molecular dynamics (aMD) simulation was performed¹⁷ for 200 ns, with The Amber15 Package,¹⁸ generating an NVT ensemble. Conformations from this preliminary simulation were grouped in clusters, based on root mean square deviations of the structure (rmsd) using Cpptraj.¹⁹ Penetratin coordinates corresponding to the node with greater population (67%), a pure $α −$helix, were selected and used as initial coordinates for the following simulations in membrane.

Next, a DOPC lipid bilayer was built using the Charmm-Gui Lipid Builder,²⁰ and the Lipid14 force field.²¹ We should note that Lipid14 was tested and validated for several models of lipid membranes, including the one used here, the DOPC bilayer. Four copies of the previously equilibrated Penetratin were taken, randomly rotated, and placed in a TIP3P water box, maximally away from each other and from the membrane. Also, a KCl salt was added, resulting in a neutral solution, with salt concentration 0.02 M. All the molecular manipulation (except the creation of the bilayer) was done with tleap.¹⁸ The size of the final system was 150 ×148×118 ų (large enough for allowing the peptides to fit in a vertical position outside the membrane), with a bilayer containing 1536 residues (512 phospholipids) in a rectangular shape, and a total of 46775 water molecules.

This system was minimized and equilibrated with constant NPT, during 100 ns, using a Langevin thermostat with collision frequency of 2 ps ^-1, and a Berendsen barostat with semi-isotropic pressure scaling and constant surface tension. Then, a production run of 400 ns was ran. In parallel, three other systems were constructed, similar to the previous one, but having only one peptide each. The size of these systems was 67×67×94 ų, and they were also minimized and equilibrated with constant NPT, and used to perform production runs of about 800 to 950 ns each. Data was saved every 1 ps for later analysis.

The resulting trajectories were analyzed using MSM, as implemented in MSMbuilder 2.10.²² Within this model, each trajectory is initially partitioned into a discrete number of structurally similar states, called microstates. Microstates were computed using a hybrid cluster k-center algorithm,²² considering only Cα atoms as reference, with a fixed distance of 2 Å. The Spectral Robert Perron Cluster Analysis (PCCA⁺) method²³ was used to find the minimal number of states, or macrostates, that represent the dynamics of the system.

Starting with a representative snapshot from one of the previously described systems with one peptide, we set up a simulation with US, where Penetratin was pushed, in order to force its translocation. The protocol consisted in using 50 sampling windows, separated 1 Å, where the distance between the center of mass (COM) of peptide and opposite layer of membrane were forced together, with the help of a harmonic constraint with an elastic constant of 20 kcal/(mol Ų). A MD simulation was performed for each window, for 5 ns, using the same parameters as in the production MDs. This procedure was ran four times, resulting in four independent US simulations, with an overall time of 1 μs. Two of them had the upper layer of the membrane constrained, by means of harmonic restraints on 33 residues, randomly taken throughout the layer. And the other two, had the same kind of restraints, but on the bottom layer. The results were analyzed with Grossfield lab’s implementation of the Weighted Histogram Analysis Method, WHAM (Ref. 28). Two similar but longer US runs were also performed, with windows of 20 ns (1 μs total simulation time each), in order to confirm the findings that we will describe in the next section.

It should be noted that the computational study of this kind of systems, highly charged molecules moving through a membrane, is not free of strong technical difficulties, particularly convergence and membrane deformation (see Ref. 30 for discussion about translocation of long amphiphilic chains moving through a membrane, a system with many similarities to ours). Specifically, in our study we found that the lipid bilayer was too easily deformed while pushing the peptide inside of it, in some cases to the point of reaching the boundary of the simulation box, such that it could affect the assumptions of periodic boundary conditions (PBC) (one imagines a smooth transition from one periodic cell to the other, not transitions where the membrane forms sharp angles). Ideally, only very long simulations should be performed, so that the peptide is able to find the best way to move through the lipidic medium, without excessively deforming it. We found this was not possible for our systems, at least with total simulations lengths in the order of a μs per US run. As a consequence, we decided to introduce positional restraints on the outer layers of the membrane. We found that harmonically restraining the COM of $≈ 30$ residues in one of the layers (we tried both the upper and lower layers, with similar results) was enough to achieve our goal of conserving the overall shape, thus not significantly affecting the assumptions of PBC.

While this procedure seems reasonable, inasmuch as it forces the system to remain in the assumptions of PBC, we have to also consider the possibility of that deformation being real. For example, Yesylevskyy et al³¹ observed a strong deformation of the membrane in simulations with multiple peptides. After comparing their simulations to another one with a reaction field and no PBC, they concluded the disruption of the membrane is real, and proposed it actually is the mechanism–micropinocytosis–that allows the peptide to travel through the membrane. But one must note that the PBC artifacts are real and present, since a high deformation makes the upper bilayer interact with the lower one, appearing from the contiguous periodic cell, and viceversa. The final answer to this problem would be to do simulations without PBC, or with PBC but in a very large box, such that even when highly deformed, the interactions between a membrane in one periodic cell and one in the next cell is negligible. Unfortunately, such a system needs to be very large, due to the charged nature of the membrane (phosphate groups), making it computationally too expensive to be used for free energy calculations. In any case, Yesylevskyy didn’t observe any significant deformation when studying only one peptide, so our system appears compatible with theirs.

With this understanding, we decided to use restraints on the bilayer; but the result of such procedure is that the free energy calculated is not the one for the translocation alone, but for translocation of a peptide while restraining the membrane. This results in much larger values of the free energy, but it shouldn’t significantly affect the relative values in different runs, since they have the same set of extra restraints. Thus, the conclusions of the US study should be taken qualitatively, and only relative–not absolute–values should be considered.

Analysis of the different trajectories generated, shows that Penetratin wanders freely in the liquid medium, with no apparent specific interaction with the membrane; but once it comes in contact with it, it associates and tends to stay in a membrane-bound state. More specifically, in the MD total simulation time of more than 4 μs, 70% of the time the peptide was found in the bound state. And in all of those cases, R1 is seen as the responsible of keeping the membrane-bound state of the peptide: Penetratin contacts the membrane with this aminoacid first, and this contact remains during most of the peptide-membrane association. R10 appears to be somewhat important, in terms of interaction with the bilayer, as the peptide can sometimes be observed changing from a perpendicular (with R1 making the contact) to a (approximately) parallel conformation, where R10 also reaches the membrane. This will be further analyzed in the coming sections.

The structure of Penetratin near the DOPC lipid bilayer was analyzed using the MSM methodology. All MD trajectories were grouped into a single trajectory with coordinates saved every ps. From this combined single trajectory, microstates were calculated using a hybrid cluster k-center algorithm, as explained in the methods section. A variety of structurally different conformations describing the peptide as it moves towards the lipid bilayer was found. These states consist of pure α, π and 310 helices and mixed conformations (half of the peptide in an α-helix conformation and half in turns or β sheets, and conformations where the first twelve residues close to the N terminus (1-12) are helical and residues at the C terminus remain in a turn conformation). Using PCCA⁺, microstates were re-grouped in four macrostates that reasonably describe all of the physical processes occurring during the sampled dynamics. Fig 1 shows the two representative conformations that constitute the most populated macrostates. The one in the left panel, with a population of 50 %, is a pure α-helix, with the C-terminus in a coil structure and the N-terminus slightly disordered; we call this cluster ’macrostate S0’. The other state, with a population of 48% (right panel of the figure), is the second most populated macrostate and it has an α-helical structure, except for the N-terminus, again, slightly disordered. We call this state ’macrostate S1’. The difference in population values between S0 and S1 is almost negligible, so we can assert that Penetratin, near the bilayer, is mostly helical. The only effect of the bilayer on the structure of the peptide is the disorder of the C-terminus, where the majority of lysines are localized. The structures observed in these simulations are in good agreement with previous computational findings,¹³ and with NMR experiments;⁷ there is a small difference in terms of the both termini, where some studies find turns instead of coils. Unlike other experimental results^8,24 no $β −$sheet conformations were found here.

It has been previously reported that Penetratin has the ability of interchanging its conformation in the membrane-bound state, as a necessary condition for membrane translocation.²⁵ Additionally, it has been shown for another cell penetrating peptide, Transportan, that the peptide’s secondary structure that minimizes the free energy is a function of the vertical distance to the membrane.²⁶ This observation might indicate the existence of an optimal structure for insertion. In such a scenario, it is reasonable to expect changes in Penetratin’s secondary structure over the course of the simulation, as it approaches the bilayer surface. To investigate this issue, we computed the temporal evolution of Penetratin’ secondary structure as it approaches the bilayer. This is represented in Fig 2. The figure shows that, as the peptide approaches the membrane, the region close to the C-terminus (R10-K16) changes its secondary structure, switching between helices (α, π and 310), turn, bend and coil structures. A well defined helical structure for the region I3-R10 is clearly observed for the membrane-bound state (marked in the picture by the white vertical line). The first two residues in the N-terminus form a coil, and the region close to the other end, R11-K16, contains both turn and coil conformations. The same behavior repeats in all other trajectories: mostly a helix-like conformation with the majority of structural variations after R10, particularly turns.

Penetratin’s orientation relative to the bilayer, while on the membrane-bound state, was also investigated from the simulations. The left panel of Fig 3 shows the coordinate, relative to an axis perpendicular to the membrane (z axis), of the peptide’s COM, the first arginine (R1), and the last lysine (K16), as a function of time, for the time interval in which the peptide is in the membrane-bound state. COM coordinates of the phosphate groups are also shown, for reference. As previously mentioned, we see that, in average, R1 is closer to PO₄ groups than K16; i.e., most of the conformations associated with the membrane-bound state, have R1 interacting directly with the phosphate groups at the bilayer, while lysines at the C-terminus tend to adopt a coil conformation pointing towards the water phase and interacting with ions present in the solution. Conformations with lysines closely interacting with the membrane are also observed, but they are more transient. This difference in interaction strength and penetration depths between arginine and lysine residues, with respect to the membrane, can also be illustrated using radial distribution functions (rdf).

Fig 4 shows the rdf of the heavy atoms at the end of all the arginine and lysine side-chains of Penetratin, with respect to the phosphorous atoms in the phospholipids, at the membrane’s surface. Both positive residues have a first peak around 3.8 Å, but the one corresponding to arginine is significantly higher, evidencing a stronger contact between itself and the phospholipids at the bilayer. The interaction between each positive side chain and the core of the membrane, represented by the acyl chain residues, is also shown (Fig 4, right panel). As expected, R1 and R10 are the residues showing deepest insertion, while residues at the C terminus (K15, K16) are the ones that locate the furthest away from the membrane’s core.

The stronger attraction observed between arginines and phosphate groups can be explained by the fact that the number of hydrogen bonds (HB) formed between these species is about one order of magnitude larger than the corresponding to lysine-phosphate, with the former being about 5 per peptide (Fig 5). This observation, somewhat expected as arginine has more available hydrogens than lysine (5 vs 3) in the charged end of the side-chain, serves as an explanation of the intriguing result from Paramagnetic Relaxation Enhancement (PRE) experiments,¹⁰ where the reported arginine-phosphate distance is 4.0 Å for (⁠$R 10 , C γ ) − P O 4$⁠, and 4.2 Å for (⁠$K 13 , C ϵ − P O 4$⁠), a small but significant difference. Related, there is strong evidence that lysine to arginine mutations result in improved translocation efficiency.²⁷ This improvement could be explained by the larger number of HB formed by the mutant arginines, suggesting a penetration mechanisms led by arginines and their direct strong interactions with the membrane.

As discussed in previous sections, analysis of the different MD trajectories shows that the peptide’s positive residues, in particular arginines, are the ones that make the first stable contacts with the phosphate heads of the bilayer; they insert themselves as hooks in the membrane (see Fig 3, right panel), and start penetrating. In order to better understand the details of this association and later crossing of the peptide through the lipid bilayer, we performed six US simulations of the translocation of Penetratin, with a total time of 3 μs.

Visual inspection of the generated US trajectories shows that the charged residues are the ones that lead the way throughout the whole translocation. There are differences between the six US runs performed but, in general, after association the peptide starts entering the membrane in an approximately parallel (to the membrane) orientation. The region around R10-K13 is the one that goes ahead, pulling the peptide (see left panel of Fig 6). As it traverses the membrane, the orientation slowly changes, becoming perpendicular to the layers. Lastly, one of the charged residues in the 2^nd half of the peptide (R10, R11, K13, K15 or K16) breaks out of the final layer, with the rest of the peptide following.

The passage of the peptide through the bilayer is not accompanied by water molecules. Contrary to the findings of Yesylevskyy et al,³¹ where they observed the formation of a water-filled pore (same was found by Huang and Garcia,³² using a cyclic Arg₉ peptide) through which the peptide was translocated, we didn’t observe the passage of water molecules. Instead, the peptide traversed dehydrated. Since both the formation of a pore and the dehydration of the peptide are energetically unfavorable processes, this is the result of the competition between those two effects: a small transient pore with dehydration vs a larger and stable pore where water can pass too. Our study shows, thus, that under the present conditions, the first option is energetically favorable.

The secondary structure of Penetratin in one of the US simulations is shown in the right panel of Fig 6. It has a mainly helical form, when associated with the bilayer, and it retains that shape between the 3^rd or 4^th residues from the N-terminus, up to the 10^th or 11^th aminoacid, while it immerses in the membrane. In five of the US simulations that region is alpha-helical, while in the other it appears more disordered, with a combination of turns and $α −$⁠, $π −$ and 310- helices. This picture is qualitatively reproduced in all six US simulations. The left panel of the same figure shows a representative snapshot of the peptide while it is fully immersed in the membrane.

Fig 7 displays the free energy profile of translocation (with a restrained layer, as previously mentioned). More precisely, it is such a profile starting at a time when the peptide is completely out, but close to the membrane, to the time when part of the peptide, not all of it, breaks out of the opposite layer. The abscissa measures the distance of the peptide’s COM from the layer initially at the other side of the membrane. That is, Penetratin starts around 49 Å away from the furthest layer, from which it is being pulled, and a few Å from the closest one. Then the pull starts, and the free energy barrier starts increasing with decreasing acceleration, until the peptide’s COM is around 20 Å, in the middle of the membrane. Then, the increase continues at a constant rate. At approximately 18 Å, four of the simulations sharply reach a plateau (red tonality lines in Fig 7), while the other two keep on growing to a value $≈ 2.7$ times larger than the other (blue tonalities). These two last cases represent, thus, highly unfavorable translocation paths. But, what is it that the paths with lower barriers have in common? Inspection of the corresponding trajectories shows that the peptide’s C-terminus fully extends in the last part of the penetration, right before the break-out occurs. On the contrary, this does not happen in high barrier cases, where the extension is prevented by the formation of one or two HBs between residues R10/R11 and K16: they stay connected, with antiparallel alignment, not allowing the C-terminus to break out from the membrane. This is true also in one of the low barrier cases, where the C- and N-termini get together, with a strong R1-K16 interaction (HB mediated), the difference here being that the $α −$helical structure of the peptide is mostly lost to a varying combination of different helices, turns, bends and coils, highly disordered (Fig 8, right panel). Thus, it appears as plausible that the formation of HBs between the charged aminoacids of the peptide prevents the final opening it needs in order to break the outer layer of the membrane; the only way Penetratin can still succeed to emerge is by loosing its $α −$helical shape, and becoming more disordered. Examining the respective trajectory, it is observed that R10 (Fig 8, left panel) then points outward and takes the lead, carrying the rest of the system behind it.

Due to technical reasons–the peptide was being pulled from the membrane’s layer (the one located at 0 Å)–a new set-up was needed in order to continue the simulations further; but this introduced a discontinuity in the profile that was not easy to remove. As a consequence, we stopped the simulations at this point, when Penetratin was starting to exit the membrane. But there is no reason not to expect that the highest point in the profile was already reached in the low barrier cases (which happen to be the relevant ones), as the peptide is already partially out, and that free energy will start decreasing from this point on, as the translocation process completes.

As detailed in the Introduction, the secondary structure reported in the literature for Penetratin, in the vicinity of a lipid bilayer, is varied and controversial. In this work, this issue was investigated by performing several MD simulations of the peptide in a box of water containing a DOPC bilayer. The net MD simulation time was 4 μs. The trajectories were analyzed using an exhaustive MSM analysis. It was found that the preference for helical structures observed in water, remains near a DOPC bilayer, except for both termini which show a tendency for turn or coil conformations: the C-terminus and residues close to it, have the highest mobility, and they loose their structure when Penetratin approaches the membrane.

In the membrane-bound state, the peptide orients with its principal axis forming a small angle with the interface. R1 and R10 closely interact with the phosphate groups at the bilayer, as evidenced by the constructed radial distribution functions in Fig 4. This interaction between the arginine side-chains and the phosphate groups is enhanced by the greater number of HBs between them, as compared to the lysine-phosphate interaction. The C-terminus, on the other hand, through K15 and K16, seems to control Penetratin’s direction of entrance. The evaluation of the number of HBs between the charged residues and the phosphate groups, R-PO4 and K-PO4, with the largest value being on the side of the first one, shows the stronger interaction of arginines and suggests a mechanism for penetration that is initiated by them. Thus, we proceeded to the study of Penetratin’s translocation, using Umbrella Sampling for the evaluation of free energy. Six simulations were performed, with a net computation time of 3 μs. The US results don’t show that conclusion (that arginines are more important in the translocation than lysines) as correct; in fact, analysis of the trajectories shows that both positive residues, arginines and lysines, are responsible for pulling the peptide across the membrane. More specifically, the residues close to the N-terminus (but not the last ones) are the ones leading the way: R10, R11 and K13. Once Penetratin gets a few Å away from the second layer, the C-terminus overtakes those pulling residues and takes the lead, finally breaking out of the membrane. This final break-out of the C-terminus was repeated in four out of the six US simulations.

Thus, the arginine residues appear to be most important in the association stage, where their stronger interaction via formation of HB stabilizes the membrane-bound state. During the translocation process, though, both charged residue types appear equally important, as they both pull the peptide. And the they are equally important, too, in building HBs while immersed in the membrane, fact that can cause the failure of the translocation process. In the final stage, the situation appear opposite to the initial one, as it is the C-terminus containing three lysines the one that eventually pierces the membrane’s last layer, and let’s Penetratin out. We could summarize these findings by saying that, although the two types of residues closely interact in this process, the effect of arginines is, to some extent, to bring the peptide into the membrane, while the effect of lysines is, once the peptide is inside, to take it out, thus completing the translocation process.

This work was supported by Agencia Nacional de Promoción Científica y Tecnológica, under awards PICT-2015-1706 and PICT-2015-3832.