Membrane adsorption, folding, insertion and translocation of synthetic trans-membrane peptides

Abstract

Spontaneous membrane adsorption, folding and insertion of the synthetic WALP16 and KALP16 peptides was studied by computer simulations starting from completely extended conformations. The peptides were simulated using an unmodified all-atom force field in combination with an efficient Monte Carlo sampling algorithm. The membrane is represented implicitly as a hydrophobic zone inside a continuum solvent modelled using the generalized Born theory of solvation. The method was previously parameterized to match insertion energies of hydrophobic side chain analogs into cyclohexane and no parameters were optimized for the present simulations. Both peptides rapidly precipitate out of bulk solution and adsorb to the membrane surface. Interfacial folding into a helical conformation is followed by membrane insertion. Both the peptide conformations and their location in the membrane are strongly temperature dependent. The temperature dependent behaviour can be summarized by fitting to a four-state model, separating the system into folded and unfolded conformers, which are either inserted into the membrane or located at the interfaces. As the temperature is lowered the dominant peptide conformation of the system changes from unfolded surface bound configurations to folded surface bound states. Folded trans-membrane conformers represent the dominant configuration at low temperatures. The analysis allows direct estimates of the free energies of peptide folding and membrane insertion. In the case of WALP the quality of the fit is excellent and the thermodynamic behaviour is in good agreement with expected theoretical consideration. For KALP the fit is more problematic due to the large solvation energies of the charged lysine residues.

• Implicit membrane
• generalized Born
• membrane protein folding
• Monte Carlo
• WALP peptide
• KALP peptide

Introduction

Computer simulations of membrane proteins in their native lipid-bilayer environment are usually limited to the <100 ns timescale, if explicit models for solvent and membrane are used. This is due to the large number of non-bonded interactions that need to be evaluated for such large and complex systems and the short time step needed to follow covalent bond vibrations in detail [1–4]. Over the last few years several methods have been proposed to increase the timescales accessible to simulations by reducing the computational effort [5], [6]. The generalized Born implicit solvent approach [7] and its extension to biological membranes [8–11] has proved particularly useful for long timescale studies of processes which require an atomic detail solute description. This has led to a number of ab-inito folding studies successfully predicting the experimentally observed native conformations of both globular [12–15] as well as membrane associated proteins [16–18]. Combining the generalized Born based implicit solvent approach with a recently developed highly efficient Monte Carlo (MC) method [19], which is particularly good at crossing local barriers in the folding energy landscape now permits very detailed atomic resolution studies of the folding process of small peptides [20], [21]. In the present study we use our previously developed implicit membrane method to demonstrate that the timescales that can be reached by these efficient computational methods are now suitably long (> µs for the present peptides) for the thermodynamic properties of the system to be fully converged.

We study the temperature-dependent adsorption, insertion and folding of two small synthetic peptides (WALP16 and KALP16), which are part of a family of artificial membrane inserting helices that have been designed by Killian and colleagues [22], [23] to explore how lipid bilayers adapt to accommodate peptides of different hydrophobic lengths. The peptides consist of a tryptophan or lysine-flanked hydrophobic core, which can be varied in length. Both their simplicity and the large body of available experimental data provides computational studies with a valuable test system to calibrate new algorithms and validate simulation methodologies concerning membrane associated peptide insertion into lipid bilayers [24–26].

The aim and motivation of this study was to address several key scientific questions regarding the simulation of membrane associated peptides using an implicit membrane: (i) Can current simulation algorithms and molecular force fields accurately predict not just the native conformation but also the temperature dependent behaviour of the peptide? (ii) Are the timescales that can be realistically reached in such simulations sufficient for the results to be fully converged? And (iii) How do the approximations entailed by the implicit solvent approach affect the results of such studies and how do they compare to fully explicit and coarse-grain simulations?

Methods

The generalized Born membrane

The development of the present generalized Born (GB) membrane has been described in detail in a previous publication [11]. The membrane was treated as a slab in a uniform aqueous solvent with a dielectric constant ε_water=80, that becomes increasingly inaccessible to the solvent towards its centre. Both the protein interior and the membrane are assumed to have the same interior dielectric constant of ε_membrane=2. The Born radii are calculated using the fast asymptotic pair-wise summation of Qiu and Still [27]:1

where P₁-P₄ are the parameters determined by Qiu et al. [27] and ccf is a close contact function. This method has been demonstrated to yield excellent results in predicting experimental free energies of solvation as well as hydration effects on conformational equilibria [28]. By modifying the pair-wise summation to solute atoms the self-solvation terms Γ(z_i, L) as well as the atomic volumes V(z_i) were made to vary smoothly between full solvation and a limiting value for burial at the centre of the membrane. We use a Gaussian shape:2

where g_bulk is the limiting value of Γ at a large distance from the membrane (i.e., z≫L) corresponding to the self-solvation term of the unmodified generalized Born method , while g_centre is the value of Γ at the membrane centre. We used a Gaussian with γ = − 3.0 and membrane half width of L = 15 Å, which corresponds roughly to the hydrophobic profile of a DPPC lipid bilayer. A value of g_centre= − 7.67 kcal/mol was used, as reported previously [9], [29].

The non-polar part of the solvation free energy is modelled using an effective surface tension associated with the solvent accessible surface area (SA) [27]. As it is moved towards the centre of the membrane the surface energy contribution of each atom is scaled down by an exponential switch at the membrane interface (z = ±15 Å). Thus for distances far from the membrane (i.e., z≫L) the nonpolar contribution is included with the positive surface tension of solvation in water, while in the centre of the membrane the surface tension is negative, (i.e., energy is gained by moving into the membrane from the gas phase) as determined experimentally [30]. The surface tension contribution of each atom was varied between the limiting values of 12 cal/molŲ in bulk solvent and −19 cal/molŲ at the membrane centre, determined by fitting to the water→cyclohexane transfer free energies of hydrophobic side chain analogs [30]. No parameters were changed or optimized for the present simulations.

Simulation set-up

The peptides have the sequence AXX-(LA)₅-XXA, with X = W for WALP16 and X = K for KALP16. The C- and N-termini were acetylated, and amidated, respectively. Twenty replicas of each peptide were simulated for 2×10⁹ Monte Carlo steps with extended starting configurations. Two replica exchange simulations (REMC) were carried out each for WALP and KALP, first one with the peptide initially spanning the membrane (TM) and second with the peptide fully solvated and oriented parallel to the plane of the membrane 20 Å away from the interface (SB) (see Figures 1 and 2). For the KALP peptide the SB starting simulation has both N-terminal lysine residues in their neutral form to allow membrane insertion, while in the TM starting simulation all four lysine residues are charged. In addition a standard simulation (SMC) of the WALP peptide was performed at 353 K for 10×10⁹ MC steps. A recent study comparing the sampling of Monte Carlo and molecular mechanics protein folding simulations shows that this corresponds to ∼1.25 µs of molecular dynamics [21].

Figure 1.  Folding of the WALP16 peptide at 300 K. Panel A: Starting from an extended conformation in solution (35 Å from the centre of the membrane) the peptide rapidly adsorbs to the membrane surface retaining its extended configuration. Interfacial folding is followed by membrane insertion. Panel B: Starting from an extended trans-membrane configuration the peptide first exits the membrane to fold at the interface, it subsequently reinserts to form the native state. The 3 chief states of the system are highlighted. (This Figure is reproduced in colour in Molecular Membrane Biology online).

Figure 2.  Representative conformers in the folding process of the KALP16 peptide in the 300 K replicas. Panel A: Surface-bound folding and repeated insertion of KALP with neutral N-terminal lysine residues. Panel B: Trans-membrane folding of KALP with all four terminal lysine residues in their charged states. The trans-membrane simulation displays little to no beta-sheet content, while it is more prominent in the surface-bound simulation. The native state conformers are identical (highlighted). (This Figure is reproduced in colour in Molecular Membrane Biology online).

The simulations were run with a Monte Carlo program developed by the authors specifically for the simulation of proteins in a GB/SA continuum solvent. An efficient concerted rotation sampling technique [19] is used to move the protein backbone; in addition there are single rapid side-chain moves, with a ratio of 3 side-chain moves per backbone move. The potential energy was evaluated with the OPLS all-atom force field [31], and the Monte Carlo simulations used Metropolis sampling. All non-bonded interactions as well as the generalized Born energy were fully evaluated (i.e., no cut-off), and the Born radii were recomputed for every configuration. The implicit solvent and membrane stretch to infinity, thus no periodic boundary conditions are applied. The described set-up has been shown to perform well in sampling DNA [32] and protein folding simulations [20], where the native state of several polypeptides was rapidly determined starting from extended conformations. We recently demonstrated that this method is equivalent to molecular dynamics sampling, with both methods able to find the native state of several polypeptides with comparable computational effort [21].

Parallel tempering

The replica exchange method has recently been reviewed in detail [33]. Twenty replicas of each system were used with identical initial configuration and exponentially spaced temperatures in the range 300–500 K. Every 10⁴ Monte Carlo moves a replica swap with transition probability3

where4

is attempted. E₁ and E₂ are the total energies of two conformers at temperatures T₁ and T_2, respectively. High temperature replicas facilitate the crossing of energy barriers, while low temperature replicas extensively sample low energy conformations. This enables the efficient and increased sampling of the entire system by frequent crossing of high energy barriers.

Thermodynamic analysis

After a 300 million MC step equilibration period all states are found to be either bound to the membrane surface or inserted into the membrane. The system was partitioned into four thermodynamic states: S = surface adsorbed unfolded, SF = adsorbed folded, M = membrane inserted unfolded and MF = inserted folded. Each sampled configuration was assigned to one of these states. Inserted configurations are defined as having a centre of mass position along the membrane normal (z-axis) of −8 Å<z<8 Å and conformers are classified as folded if at least 60% of their sequence is helical.

The populations p_i of the four states and the free energy differences5

are calculated for each replica temperature over the equilibrium phase of the simulation (i.e., excluding the 300 million MC step equilibration period). Linear fits of the ΔG_ij with respect to T are then used to determine the energetic component ΔE_ij and the solute entropic component ΔS_ij for each transition ij:6

Note that ΔE_ij is not the internal solute energy, but contains the generalized Born solvation free energy, including entropic terms due to the solvent ordering near the solute surface. ΔE_ij is also determined by subtracting the bin averages of E for each temperature replica: . Finally, the populations of the individual states at each temperature as predicted by the 4-state fit are obtained by:7

Results

Folding pathway

Figures 1 and 2 show the folding pathways for the trans-membrane and solvent starting configurations of the WALP and KALP peptides, respectively. In the solvent simulations both peptides are initially located 35 Å from the membrane centre and oriented parallel to the membrane plane. All replicas rapidly adsorb to the membrane surface within 10⁶ MC steps. Low temperature replicas form continuous helices with the long axis parallel the membrane surface, while high temperature replicas remain unfolded or trapped in beta-sheet conformations. Spontaneous insertion is observed only for helical conformations, which orient to span the membrane. The helices frequently exit and reinsert into the membrane in a temperature dependent manner, with the lowest temperature replica having the largest fraction of inserted states. Insertion and exit always occurs to the same interface, i.e., no translocation takes place for helical conformers (see below). For both peptides insertion only occurs from the N-terminal side, which in the case of KALP has neutralized lysine residues.

In the trans-membrane starting configuration the peptide spans the membrane perpendicularly in an extended conformation. For WALP all replicas rapidly exit the membrane (<10⁶ MC steps) in an unfolded state to adsorb at the interfaces. Both interfaces are almost equally populated. The subsequent behaviour is exactly identical to the solvent simulations. KALP remains firmly integrated into the membrane due to the large penalty associated with translocating two charged lysine residues across the hydrophobic membrane core.

Therefore, the overall folding pathway observed can be summarized as: (i) Peptide adsorption, (ii) interfacial folding, and (iii) folded insertion. KALP follows this pathway only if the lysine residues of at least one terminus are neutralized (as in the solvent starting system, see Simulation set-up), to allow it to translocate the hydrophobic membrane core.

Membrane location and conformation

Figure 3 shows the states sampled by the WALP REMC systems at 300 K. The Figure shows a very strong correlation between tilt angle and insertion depth. Two distinct clusters are observed corresponding to a surface bound configuration at z = 11.2±1.1 Å, with an average tilt angle of 81±7° and an inserted configuration at z = 2.3±1.2 Å, with low tilt angles of 16±9°. The REMC simulations starting from fully solvated peptides sample one side of the membrane only (the insertion side), while the simulations starting from a trans-membrane conformation sample states on both interfaces since different replicas exit the membrane to different interfaces.

Figure 3.  Scatter plot showing the distribution of tilt angles and centre of mass position of the WALP 16 peptide along the membrane normal for the 300 K replicas. Panel A: The surface bound starting configuration only translocates very late in the simulation and hence only explores one interface (see text). Panel B: The trans-membrane starting configuration has an almost equal probability of exiting the membrane on each leaflet exploring both interfaces. Inserted and surface bound configurations show clear tilt preferences matching the energy landscape of a rigid body scan (see Figure 4A below). The roman numerals correspond to locations/orientations with respect to the membrane, depicted in Figure 4 (I & IV = M & MF, II & III = S & SF). (This Figure is reproduced in colour in Molecular Membrane Biology online).

This behaviour corresponds remarkably well to the insertion energy landscape (tilt versus membrane position with optimized rotation angle) for a helical conformer (see Figure 4), since at 300 K the peptide is helical for nearly 80% of the simulation time (c.f. Figure 5). The insertion energy landscape gives barrier estimates, which are difficult or impossible to extract from the folding simulations as they correspond to regions of very low sampling. The surface has four distinctive minima, two corresponding to trans-membrane configurations (I and IV) and two representing surface adsorbed helices (II and III). These correspond exactly to the four states sampled in the folding simulations (c.f. Figure 4).

Figure 4.  Panel A: Insertion energy landscape of a helical WALP16 conformer. The distinctive 4 minima are indicated. There are two trans-membrane configurations (I and IV), which are separated by a barrier, which is much higher than the barrier between the trans-membrane and its associated surface bound configurations (i.e., I↔II and III↔IV). The insertion energy minima match those of the 300 K replicas (c.f. Figure 3). Panel B: Orientations of the four minima in the membrane, the minimum energy pathway on the tilt/z surface connecting the four states shows the high translocation barrier, separating the states in each membrane leaflet (see text). (This Figure is reproduced in colour in Molecular Membrane Biology online).

Figure 5.  Thermodynamic analysis of the trans-membrane WALP peptide simulation, the surface bound simulation shows identical behaviour. Panel A: Average peptide helicity and beta-sheet content as a function of replica temperature (excluding equilibration). Panel B: Four-state fitting of all sampled states (S = surface unfolded, SF = surface folded, M = membrane unfolded, MF = membrane folded). The M state is hardly populated, reducing the system to three main states: S, SF and MF. The dots represent the direct populations measured for each temperature replica, while the thick line is the prediction from the thermodynamic four state model fit (see Methods). (This Figure is reproduced in colour in Molecular Membrane Biology online).

The two inserted configurations are not centred on the membrane but one terminus is always at an interface. Thus the helix does not fully span the membrane but is still oriented along the membrane normal. Each inserted configuration is separated from its associated surface bound state by a ∼5 kcal/mol barrier (e.g., I↔II), which is much lower than the ∼17 kcal/mol barrier to its opposite surface bound configurations (i.e., I↔III). The helix therefore inserts and exits the membrane remaining in one bilayer leaflet only. The profile of the barrier on the tilt/z surface is given in Figure 4B.

The height of the barrier is also evident from the very low number of translocation events, which restrict sampling to one leaflet (c.f. Figure 3A).

Spontaneous translocation was observed for two high energy replicas of the WALP peptides (see Supplementary data online).

The KALP peptide shows similar behaviour, but the energy landscapes are dominated by the large penalty associated with inserting lysine residues into the hydrophobic core. Realistic treatment of this system therefore requires dynamic protonation (see Discussion).

Thermodynamic behaviour

The replica exchange method significantly decreases the simulation time required for the system to reach equilibrium. In addition it provides thermodynamic information over a wide temperature range. Analysis of the REMC and SMC WALP trajectories showed that at equilibrium the system could be subdivided into four distinct states: (i) A surface adsorbed unfolded conformation (S), (ii) A surface adsorbed alpha helix (SF), (iii) An unfolded membrane inserted configuration (M), and (iv) A folded trans-membrane helix (MF) (see Methods for the definition of ‘unfolded’ and ‘membrane inserted’ and Figure 6 for a graphical visualization of each state). Since membrane adsorption is rapid and completed for all replicas after 10⁶ MC steps, no state that is significantly distanced from the membrane is observed in the equilibrium phase. Included in the S state are many beta-structures. Figure 5A shows the average helicity and beta-sheet content of the WALP peptide as a function of the system temperature. Low temperatures strongly favour helical conformations, while beta-sheets and unfolded states are the dominant secondary structure at elevated temperatures. Figure 5B shows the temperature dependence of the populations of the four states. The trans-membrane helix is the dominant conformation for low temperatures. Slight elevation of temperatures increases the fraction of surface bound helices, while at high temperatures unfolded surface bound conformers dominate. The population of unfolded membrane inserted configurations increases with temperature but is generally very low.

Figure 6.  Temperature dependence of the transition free-energy between the four states (M, S, SF, MF). The four states and main associated transitions are shown at the bottom of the figure. Surface folding (S→SF) can be seen to be favourable below ∼350 K and the insertion of folded helices (SF→MF) is favourable below ∼370 K. Unfolded insertion (S→M) is strongly unfavourable at all temperatures as is intra-membrane unfolding (MF→M). Dots and thin lines represent the directly measured values, while the thick straight lines are the predicted thermodynamic behaviour from the four-state fitting (see Methods). (This Figure is reproduced in colour in Molecular Membrane Biology online).

The corresponding probabilities obtained from a four state thermodynamic fit are shown as thick lines in the same figure, and the fitted free energies are shown in Figure 6A. The quality of the four state fit is excellent, with r²>0.98 and rms errors of ∼0.15 kcal/mol. The results of the fit are shown in Table I. The MF state is stabilized by ΔE = 4.9 kcal/mol with respect to the SF state, and also has the lowest solute entropy, as expected. Note that ΔE contains both the intra-solute energy as well as the solute-solvent free energy. Both the S and M states are much higher in energy. The S state has the highest solute entropy, as expected. The values of ΔE were also obtained by averaging E for each of the four states at all 20 temperatures. The match to the estimate of ΔE from the thermodynamic fit is good, except for the poorly sampled M state, which exhibits considerable fluctuation. In the same table, the free energy differences at 300 K are given, both directly calculated by Equation 5 and as obtained from the fit. At this temperature, the folded states clearly dominate, with a stabilization of 0.75–0.91 kcal/mol for the MF state over the SF state. This is in good agreement with the values from the insertion energy landscape (Figure 4).

Table I.  Entropic and energetic contributions to the free energies of the four states (M = membrane inserted unfolded, S = surface adsorbed unfolded, SF = surface adsorbed folded, and MF = membrane inserted folded) of the WALP peptide, as determined from a thermodynamic four-state fit. Shown are the relative energies ΔE (solute internal energy plus generalized Born solvation free energy, including solvent entropy), and solute entropies ΔS. The fitted free energy ΔG is given at 300 K. Also shown is ΔG as determined directly from Equation 5, and ΔE as determined from direct bin averages. The folded membrane inserted state (MF) has the lowest free energy, energy and entropy and is defined as the reference state zero.

State	ΔG300K [kcal/mol]	[kcal/mol]	ΔE [kcal/mol]	ΔEfit [kcal/mol]	ΔSfit [cal/mol·K]
M	3.84	4.01	27±5.3	10.9	22.9
S	1.64	1.83	14±4.0	11.6	32.6
SF	0.75	0.91	4.7±3.5	4.9	13.4
MF	0.00	0.00	0.0±3.7	0.0	0.0

The transition probabilities between the states are highly temperature dependent (see Figure 6): Unfolded insertion (S→M) is unfavourable for all temperatures, while for folded helices insertion is favourable below ∼375 K (SF→MF). Helical folding of surface bound configurations is favourable below ∼350 K (S→SF). The lowest curve shows that the penalty for intra-membrane unfolding is extremely high for all but the highest temperatures (M→MF). This thermodynamic description matches the quantitative pathway description presented above. It should be noted that the M state is very sparsely populated resulting in deviations from the thermodynamic straight line fit. All data presented here was from the trans-membrane simulation. Results for the surface bound simulation are identical, demonstrating that the system has been thoroughly sampled in both cases.

Figure 7 compares the behaviour of the 353 K replica with a standard simulation at the same temperature. The initial non-equilibrium structural collapse to compact conformations takes 1–2 billion MC steps and is very similar in both simulations. After reaching equilibrium the 353 K simulation swaps between folded surface bound (SF) and membrane inserted (MF) configurations. Table II shows that the probability for the dominant MF state is identical in both simulations. However, the 353 K replica trajectory contains a significant amount of unfolded surface bound peptide (S) due to the occasional swapping with only partially folded conformations from other replicas. This also lowers the average helicity compared to the 353 K simulation. Nevertheless the overall helicity is reduced when the peptide is located at the surface as is evident from the correlations seen when plotting helicity together with the position along the membrane normal in Figure 8. The standard simulation displays long-timescale oscillation (τ≈2.5 billion MC steps) between surface bound and trans-membrane configurations, not observed in the replica trajectory.

Figure 7.  Panel A: WALP peptide helicity (α) and position along the membrane normal (Z) at 353 K for the first 2 billion MC steps of both MC and REMC simulations. The REMC simulation exhibits smaller overall helicity due to the occasional swapping with only partially folded conformations from other replicas. Panel B shows the time-evolution of probabilities of the four states for the 353 K replica, while Panel C shows the equivalent behaviour for the 353 K simulation. The initial non-equilibrium folding process takes 1–2 billion MC steps and is similar in both simulations. The REMC simulation has a non-zero probability of the S state, while it decays to zero in the 353 K simulation. (This Figure is reproduced in colour in Molecular Membrane Biology online).

Figure 8.  Secondary structure and membrane position of the WALP 353 K simulations. Panel A: Peptide helicity (blue line, left axis) and position along the membrane normal (z-position, red line, right axis) versus simulation time (MC step). Panel B: Simulation averages of secondary structure for each residue. (This Figure is reproduced in colour in Molecular Membrane Biology online).

Table II.  Helicity, position along the membrane normal (Z) and probability of the 4 states of the system for the 353 K replica (REMC) and a standard 353 K simulation (SMC).

	Helicity [%]	Z-pos [Å]	S	SF	M	MF
REMC	71.4±3.8	5.6±0.6	0.19±0.04	0.21±0.04	0.00±0.00	0.60±0.06
SMC	88.9±4.8	5.7±3.1	0.00±0.00	0.40±0.35	0.00±0.00	0.60±0.35

Figure 9 shows that the KALP peptide also exhibits temperature dependent behaviour. However, the fluctuation in peptide helicity are much larger than those of the WALP peptide (c.f. Figure 5), and the loss of helicity at higher temperatures is not compensated for by an increase in beta-structure. The temperature dependence of each of the four states probability roughly follows that of the WALP peptide. Nevertheless, the corresponding free-energy analysis reveals that of the six possible transitions only two (S→M and SF→MF) show thermodynamic behaviour, suggesting that the simulations are not fully converged. In particular the S→SF transition, which represented the major folding pathway for the WALP peptide, exhibits a strongly non-linear temperature variation. Grouping the system into inserted (SA = S + SF) and extra-membrane states (MA = M + MF) regardless of the fold shows that unlike the WALP peptide the membrane inserted configuration does not correspond to the free-energy minimum of the system. This is due to the large penalty associated with partial burial of the charged lysine residues at the membrane interfaces, suggesting that a dynamic treatment of ionizable residues in combination with a more realistic model of the complex lipid head-group region is required to accurately simulate this system (see Discussion).

Figure 9.  Thermodynamic analysis of the surface bound KALP peptide simulation. Panel A: Average peptide helicity and beta-sheet content as a function of replica temperature (excluding equilibration). Panel B: Temperature dependence of the four chief states of the system. Panel C: Temperature dependence of the transition free-energy between states. SA = surface bound, MA = membrane inserted (irrespective of the fold). (This Figure is reproduced in colour in Molecular Membrane Biology online).

Discussion

Peptide adsorption, insertion and folding

The prevailing current theory of spontaneous peptide insertion into membranes [34–37] can be summarized as a three-stage process. First the peptide adsorbs to the membrane interface in an unfolded state. Interfacial helix formation lowers the solvation energy by largely burying the polar peptide backbone. Once folded the helix inserts spanning the hydrocarbon core. This model is based on theoretical and experimental considerations derived from bulk hydrophobic solvents, which show that unfolded insertion would entail enormous energy penalties of the order of 6–8 kcal/mol/residue, associated with de-solvating the peptide backbone [38–40]. The helical conformation largely buries the backbone, thus removing the insertion penalty.

The present simulations follow this model closely (see Results). A similar pathway was also found by a recent study of WALP16. Using a different generalized Born-based implicit membrane algorithm [10] in combination with the param22 all-atom force field [41], [42] Brooks and Im were able to simulate spontaneous insertion and folding of WALP16 by simulating 32 replicas in temperature range 300–800 K for 14 ns [25]. They found that membrane insertion required the peptide to have at least 80% helicity. Similarly a recent long timescale (600 ns) coarse-grain simulation of helically restrained WALP16 conformers in a 256 DPPC lipid bilayer found the system to oscillate between a trans-membrane (population = 41%) and surface bound state (59%), with average lifetimes of ∼10 ns each [43].

Woolf and co-workers recently reported folding of WALP16 in a fully explicit membrane set-up. 64 replicas spaced between 350–800 K each consisting of the peptide, a 36 DPPC lipid bilayer solvated with 1048 TIP3P water molecules [24]. NVT simulations of the systems for a total of 2.6 ns using the param22 force field resulted in a different insertion pathway than described above. Initial interfacial adsorption is followed by an unfolded insertion with subsequent intra-membrane helix formation. No folded insertion was observed. Even though unfolded insertion entailed a large 5–10 kcal/mol energy penalty per residue, peptide insertion was made possible by entropic compensation. Subsequent intra-membrane folding was driven by a very sharp energy decrease, consistent with our findings.

The difference in folding pathways is intriguing and highlights a major shortcoming of the implicit membrane model, which cannot currently account for entropic effects due to lipid rearrangements or hydrophobic mismatch. At high temperatures unfolded insertion is observed frequently in the implicit model. However, since the insertion barriers are significantly lowered due to the continuum representation these conformations rapidly exit the membrane upon cooling (i.e., at lower replica temperatures). This means that the unfolded insertion pathway is almost certainly under-represented in the implicit model. Quantifying the importance of each pathway would require very long timescale explicit lipid-bilayer simulations observing a large number of insertion events. Unfortunately this remains utterly unfeasible with present simulation methodologies, and no experimental data is currently available to ultimately clarify this point.

Effects due to the lipid bilayer

Experimentally WALP16 is know to form stable membrane spanning alpha-helices with low tilt angles [44], [45]<15–20° irrespective of hydrophobic mismatch [46], [47]. Several computational approaches have been used to determine the tilt angle and location in the membrane. Poisson-Boltzmann rigid body scans of helical WALP19 and WALP20 peptides, which have slightly longer hydrophobic segments in a 25 Å width low dielectric slab (ε_membrane=2) gave tilt angles of ∼0° [48]. Short timescale (1.5 ns) molecular dynamics simulations of helical WALP16 conformers in explicit DPPC and DMPC lipid bilayers found mean tilt angles of ∼9±3° [26]. Brooks and co-workers report an average tilt of 12±7° for the 300 K WALP16 replica of their implicit membrane simulation similar to the 16±9° found in the present study [25]. A recent coarse-grain simulation of helically restrained WALP16 conformers in a DPPC lipid bilayer found mean tilt angles of 23±13° [43].

The trans-membrane configurations encountered in the present study do not fully span the membrane since one pair of tryptophan residues firmly anchors one helix terminus in one interface, restricting the motion of the peptide to one leaflet (see Results). This is a result of ‘negative’ mismatch, where the peptide length is shorter than the hydrophobic membrane core. Similar behaviour was also observed in coarse-grain simulations of WALP16 in DPPC lipid bilayers [43], which have a hydrophobic core width of ∼34 Å, exceeding the ∼25 Å length of helically constrained WALP16. Interestingly the study found no translocation for WALP16, while simulations of WALP19, which has a longer hydrophobic core, showed several such events, indicating that ‘negative’ mismatch may cause a barrier for membrane crossing events.

Insertion for WALP was found to always occur from the N-terminal side, with the C-terminus remaining at the interface. This behaviour was also observed by Brooks and co-workers [25]. Since the peptide sequence is completely symmetric they suggest that the larger dipole moment of the carbonyl groups at the C terminus might explain this result.

Experimentally, the bilayer width in the vicinity of WALP does not change significantly with hydrophobic mismatch [46], [47]. Instead, peptide localization in the membrane seems to be determined chiefly by interfacial tryptophan anchoring [49]. Both the coarse-grain and the implicit membrane simulations seem to reproduce this result. However, the helical restraint used in the coarse-grain simulation may mask other insertion and translocation pathways involving unfolded or partially folded conformers. The present implicit membrane method on the other hand samples a large number of structurally diverse protein conformations, but does not model entropy changes of the bilayer itself. The simulations thus cannot account for local variations in hydrophobic width as well as changes in lipid ordering in the vicinity of the peptide, which might be important factors in peptide insertion and translocation [24].

Simulation accuracy and temperature dependence

We report a strong temperature dependence of secondary structure and preferred peptide location in accordance with expected thermodynamic behaviour (see Figure 5). Higher replica temperatures exhibit a rapidly decreasing helicity, which drives the peptide out of the membrane due to the high cost of exposing the peptide bond [37]. This interplay between secondary structure and membrane localization causes the overall temperature dependence of the transition free energies between the four chief states of the system. The fact that for the WALP peptide both surface and trans-membrane starting configurations converge to yield both identical and fully converged thermodynamic behaviour (i.e., the free energy of each transition varies linearly with temperature) is remarkable, and demonstrates that the system has been sufficiently sampled (see Figure 6). At equilibrium the OPLS all-atom force field captures a large range of conformations ranging from unfolded, partial helical folds, misfolded beta-sheets to fully folded alpha-helices, with non-native folds dominating at higher temperatures as expected. This behaviour differs drastically from previous simulations with a related implicit membrane model employing the CHARMM force field, which have consistently reported that all replicas rapidly fold continuous helices up to temperatures of 1000 K [10], [18], [25]. The authors suggest this is most likely due to a helical bias of the CHARMM force field [50]. Such biases may veil important thermodynamic behaviour as well as insertion pathways involving unfolded states.

The KALP peptide simulations exhibit only partially converged thermodynamic behaviour (see Figure 9), highlighting one of the major problems of current membrane protein simulations, namely: the correct treatment of ionizable residues. Currently the ionization state is set at the start of the simulation and remains fixed throughout. This prevents charged residues entering the membrane due to the large desolvation energy associated with burial in a hydrophobic environment. Experimentally, however, lysine residues can be buried deeply in membrane spanning helices with surprisingly small energy penalties [51], indicating that their integration involves deprotonation or accompaniment by a water shell. This implies that for membrane proteins containing charged residues realistic simulation of the folding and insertion process requires a dynamic treatment of ionization states.

In the present simulations the N-terminal lysine residues were neutralized to allow membrane insertion in the first place. However, the overall system is not fully sampled since the energy landscape is dominated by the two C-terminal lysine residues, which remained charged. Even small movements of these residues cause large energy changes. Indeed the simulations show that the peptide core is not hydrophobic enough to make the trans-membrane helix the dominant conformation of the system (see Figure 9). The problem is aggravated by the absence of a realistic treatment of the highly polar lipid bilayer interfaces. Explicit lipid head-groups could lower the overall energy by formation of hydrogen bonds. Together with membrane deformation this allows some degree of burial through snorkelling even when protonation states are fixed. Unfortunately such effects are not accounted for by the implicit membrane, suggesting that further developments are needed.

In summary: An implicit membrane model based on the generalized Born theory of solvation was used in combination with an all-atom force field and an efficient Monte Carlo sampling scheme to explore the spontaneous inserting and folding of two small membrane associated peptides, containing a hydrophobic core and aromatic (WALP) or charged termini (KALP) respectively. For each peptide, two simulations with 20 replicas were performed for 2×10⁹ Monte Carlo steps, giving an accumulated simulation time of 16×10¹⁰ Monte Carlo steps, corresponding to ∼20 µs of molecular dynamics [21].

Fully solvated extended starting configurations rapidly adsorb to the membrane surface, while extended trans-membrane starting structures quickly exit the membrane to locate themselves at the interfaces. Interfacial folding is strongly temperature dependent, with low temperatures driving helix formation, while elevated temperatures retain a large amount of extended and beta-structure. Spontaneous membrane insertion is observed for helical conformers at low temperatures. The very long timescales of the simulations together with the exponentially spaced 20 replicas allowed a thorough analysis of temperature dependent properties. Fitting showed that the system can be described by four chief states: Surface bound unfolded, surface bound folded, membrane inserted unfolded and membrane inserted folded.

In the case of the WALP peptide the thermodynamic behaviour was found to be fully converged. At low temperatures the membrane inserted folded state dominates, while intermediate temperatures favour the surface bound folded state. At high temperatures unfolded surface bound states prevail. The relative transition free energies between the states are in good agreement with theoretical expectations and experimental estimates of similar systems. For the KALP peptide the strong energy penalties associated with partial membrane burial of charged residues prevent a convergence of the thermodynamic averages over the timescales investigated, highlighting a general problem in simulations with fixed ionization states, which needs further attention.

Nevertheless, the improved timescales and sampling provided by implicit membrane representations in conjunction with accurate force fields and efficient replica exchange methods is providing a useful tool for previously unfeasible studies of the thermodynamic behaviour of membrane associated peptides at an atomic level.

Supplementory material

Spontaneous translocation

Two spontaneous translocation events were observed in the REMC solvent starting stimulations of the WALP peptide. Translocation occurred at high temperatures (487–500 K), taking less than 10⁶ MC steps (see Figure 10). In each case the intermediary membrane inserted conformers contained significant beta-structure and no helical content. Since both translocations occurred very late in the simulation (at ∼1.9×10⁹ steps), the remaining simulation time was insufficient to allow the cooling required for sampling the low temperature regions of the second interface (c.f. Figure 3A). However, the events demonstrate that in principle all states of the system are accessible starting from a bulk solvated peptide.

Figure Supplementory.  Two spontaneous translocation events occurred for the WALP peptide during the surface bound simulations. In both cases translocation occurred at high temperatures (487–500 K), taking less than 10⁶ MC steps. Panel A: The 362 K starting simulation. Panel B: The 461 K starting simulations. Both intermediary conformers contained large amounts of beta-structure and no helical content. Panel C: Simulation temperature (green = running average over 100 frames) of the 362 K starting simulation, the translocation event is indicated by a circle (at MC step 1.89×10⁹).

Acknowledgements

This work was partially funded by an International Fellowship of the Wellcome Trust (MBU) and an Emmy-Noether Fellowship of the Deutsche Forschungs-gemeinschaft (JPU). We thank Professors Antoinette Killian, Charles Brooks and Thomas Woolf for interesting discussions. Declaration of interest: The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.