The interaction between membrane proteins and the surrounding membrane is becoming increasingly appreciated for its role in regulating protein function, protein localization, and membrane morphology. In particular, recent studies have suggested that membrane deformation is needed to stably accommodate proteins harboring charged amino acids in their transmembrane (TM) region, as it is energetically prohibitive to bury charge in the hydrophobic core of the bilayer. Unfortunately, current computational methods are poorly equipped for describing such deformations, as atomistic simulations are often too short to observe large-scale membrane reorganization and most continuum approaches assume a flat membrane. Previously, we developed a method that overcomes these shortcomings by using elasticity theory to characterize equilibrium membrane distortions in the presence of a TM protein, while using traditional continuum electrostatic and nonpolar energy models to determine the energy of the protein in the membrane. Here, we linked the elastostatics, electrostatics, and nonpolar numeric solvers to permit the calculation of energies for nontrivial membrane deformations. We then coupled this procedure to a robust search algorithm that identifies optimal membrane shapes for a TM protein of arbitrary chemical composition. This advance now permits us to explore a host of biological phenomena that were beyond the scope of our original method. We show that the energy required to embed charged residues in the membrane can be highly nonadditive, and our model provides a simple mechanical explanation for this nonadditivity. Our results also predict that isolated voltage sensor segments do not insert into rigid membranes, but membrane bending dramatically stabilizes these proteins in the bilayer despite their high charge content. Additionally, we use the model to explore hydrophobic mismatch with regard to nonpolar peptides and mechanosensitive channels. Our method is in quantitative agreement with molecular dynamics simulations at a tiny fraction of the computational cost.
The membrane serves as a barrier between the interior of a cell and its environment, allowing the cell to retain essential components and control its internal chemistry. The primary constituents of the membrane are amphipathic lipids composed of a polar headgroup that faces solution and a hydrophobic tail that forms a low-dielectric barrier preventing ions and molecules from penetrating the membrane (Parsegian, 1969). Much of the communication and exchange of material between the inside and outside of the cell is mediated by embedded membrane proteins that enable a variety of biological phenomena, from small molecule and ion transport to cell signaling. The transmembrane (TM)-spanning regions of membrane proteins are characteristically made up of hydrophobic amino acids, which are energetically compatible with the hydrophobic environment of the membrane core. However, charged residues within the TM domains of a variety of proteins present a situation that is electrostatically challenging to stable inclusion in the membrane. For example, voltage sensor domains of voltage-dependent ion channels contain charged arginines and lysines in the membrane-spanning S4 segment that make it possible for the protein to respond to electric fields (Papazian et al., 1987; Aggarwal and MacKinnon, 1996; Seoh et al., 1996; Long et al., 2005). A seminal study of the thermodynamic properties of amino acid analogues suggests that it requires 65–80 kcal/mol of energy for charged residues to enter organic phases from water (Cabani et al., 1981), and numerical calculations support these high energies (Sitkoff et al., 1994). Various models have been proposed to explain how voltage-dependent proteins are able to stably accommodate charged residues in the hydrophobic span of the bilayer. For instance, the S4 segment may be positioned within a canaliculus of the protein, where it avoids interaction with the bilayer core and makes charge pairs with other regions of the protein (Armstrong, 1981; Catterall, 1986a; Tiwari-Woodruff et al., 1997; Durell et al., 1998; Lecar et al., 2003; Pathak et al., 2007). Yet electron paramagnetic resonance experiments on KvAP suggest that some portions of the S4 segment directly interact with lipid, whereas other portions are protected from the lipid (Cuello et al., 2004). Additionally, biotin-avidin accessibility experiments on the KvAP S4 segment suggest even more extensive interactions of S4 with the lipid (Ruta et al., 2005). Studies have also attempted to explain the stability of these charged segments in terms of lipid composition (Schmidt et al., 2006; Johansson and Lindahl, 2009), but it is still unclear to what extent charged voltage sensor segments are exposed to lipid and whether they are truly stable in any lipidic environment.
Recent translocon- (Hessa et al., 2005a, 2007) and porin folding–based (Moon and Fleming, 2011) experiments have presented data in which charged residues only destabilize a membrane protein by a few kilocaries per mole. In light of these results, it is not surprising that voltage-dependent membrane proteins harbor charged residues in the membrane-spanning region, but how do we reconcile these low insertion energies with the large values derived from bulk partitioning experiments (Cabani et al., 1981)? Molecular dynamics (MD) simulations suggest that membrane bending helps stabilize charges in the membrane by allowing water access to the buried amino acid (Freites et al., 2005; Dorairaj and Allen, 2007; MacCallum et al., 2007). This is in accord with experiments and MD simulation on the voltage sensor from KvAP (segments S1–S4) showing that the bilayer thins in the presence of the voltage sensor and that the voltage sensor is significantly hydrated in the membrane (Krepkiy et al., 2009). Based on such observations, we developed a solvation model for protein insertion into the membrane that treats the membrane as a deformable continuum (Choe et al., 2008), similar to classical studies (Helfrich, 1973; Kim et al., 1998; Nielsen et al., 1998), but we couple membrane bending to protein electrostatics and hydrophobic forces in an analogous manner to theoretical treatments of small molecule solvation (Luo and Tucker, 1995; Dzubiella et al., 2006). Interestingly, both MD simulations (Dorairaj and Allen, 2007; MacCallum et al., 2007) and our continuum-based molecular calculations (Choe et al., 2008) predict significantly larger destabilization energies, on the order of 10–18 kcal/mol, than those predicted by the translocon and porin folding scales. Nonetheless, these values are much smaller than those based on bulk experiments (Cabani et al., 1981; Sitkoff et al., 1994) as a result of bilayer deformations in the presence of buried charged amino acids that expose charged groups to water and polar lipid headgroups (Freites et al., 2005; Dorairaj and Allen, 2007; MacCallum et al., 2007). Previously, we predicted the energetics of charged TM segments by using linear elasticity theory to allow for membrane bending, and our results clearly showed that these distortions facilitate favorable electrostatic interactions between charged protein moieties and solvent and polar headgroups at a minimal cost to the membrane bending energy (Choe et al., 2008). However, our method was severely limited by the ability to only allow distortions in the upper leaflet and the need to posit a priori the boundary curve of the protein–membrane interface using predetermined geometries. For simple protein sequences harboring a single charged residue, we show here that we made good guesses as to the contact curve in our original study, but for more complex peptide sequences a systematic approach must be adopted. We have therefore expanded our previous work by introducing a search on the inner boundary curve of the protein–membrane interface to find optimal distortions in the membrane. Additionally, we now allow both the upper and lower leaflets to deform.
These extensions to our original model have allowed us to probe several phenomena central to membrane biophysics. Our method successfully identifies the TM segment of a membrane protein, and it predicts the vertical position of the protein within the membrane that minimizes the total insertion energy. Unlike other continuum membrane models, our method allows for large membrane deformations that cover hydrophobic stretches of the TM protein, and it dramatically reduces the system energy. We explore the influence of varying the membrane thickness on the stability of the mechanosensitive channel MscL, and we predict a degree of water accessibility that is in good agreement with experiment. We show that insertion energy values are not additive, particularly when considering the insertion of multiple charged residues into the TM region. Finally, our method predicts that the S4 segments of voltage-dependent potassium channels are stable in the membrane, and we calculate stability values and membrane distortions that are in semi-quantitative agreement with fully atomistic and coarse-grained MD simulations; however, we show that stability values can vary greatly depending on the material properties of the bilayer. In the Discussion, we suggest additional areas of biology to which we could apply our method, as well as current limitations and future changes to the algorithm.
MATERIALS AND METHODS
As in our previous work (Choe et al., 2008), we assume that membrane protein stability is dominated by three energetic components: membrane bending, Gmem; electrostatics, Gelec; and nonpolar interactions between the protein and water, Gnp. All free energy changes are calculated by comparing the energy of the protein in a reference state completely in solution, far from an unstressed membrane with the energy of the protein embedded in the membrane. Hence, Gmem is zero for the reference state. We assume that the protein structure is the same in the reference state and the membrane-embedded state. Membrane protein stability is then given by the difference in the total energy of the protein in solution compared with the energy in the environment of the membrane:
All three energy terms on the right hand side of Eq. 1 are coupled in a complex manner through the shape of the membrane. As the membrane shape is changed, this influences the electrostatic energy of the system and the nonpolar energy by altering the amount of protein surface exposed to water. These second two energy terms drive changes in the shape of the membrane as our search algorithm minimizes ΔGtotal.
The shape and energy of the membrane are determined using linear elasticity theory, in which each leaflet is described by a thin surface, u(x,y), illustrated in Fig. 1, with material properties that can be tuned to the membrane of interest. The total membrane bending energy is expressed as follows:
where u+ is the shape of the upper leaflet and u− is the shape of the lower leaflet, L0 is the equilibrium length of the membrane, Kc is the membrane bending energy, α is the surface tension, and Ka is the compression modulus. The energy of the upper and lower leaflets are coupled through the compression modulus. The functional derivative of Eq. 2 with respect to variations in u+ and u− gives the partial differential equation (PDE) that determines the shapes of each leaflet. Assuming the opposite leaflet is flat, each surface obeys the following PDE:
where γ = α/Kc and β = 2Ka/(L02Kc). We solve this equation separately for u+(x,y) and u−(x,y) applying height and contact angle boundary conditions at the membrane–protein interface as described below and shown in Fig. 1 B, and far from the protein, we assume that the membrane is asymptotically flat at its equilibrium length, L0. After solving for u+(x, y) and u−(x, y), the total membrane energy is determined by carrying out the integral in Eq. 2. We have investigated the error associated with assuming the upper leaflet is uncoupled from the lower leaflet during the solution step, and for the class of problems examined here, the absolute error in the total energy is <0.5 kcal/mol (data not shown).
The electrostatic energy component is calculated with Poisson-Boltzmann electrostatics, and the nonpolar contribution is calculated from the solvent accessible surface area (SASA). Both of these methods take into account all of the atomic detail of the protein, and therefore, a PDB file, or equivalent file, is required to carry out these calculations. To compute the SASA, the rough surface of the protein is computed by running a water probe of 1.4 Å over the van der Waals surface created by the union of all atomic van der Waals surfaces for the individual atoms in the protein using the Shrake-Rupley algorithm (Shrake and Rupley, 1973). In the presence of the membrane, we carefully keep track of which protein atoms are embedded in the membrane, and modify the SASA accordingly. After standard convention, we ignore hydrogen atoms when calculating the SASA. To compute the electrostatic energy, we use the software APBS (Baker et al., 2001). The protein is treated in atomic detail, and the atomic partial charges are set using the PARSE parameter set, which was parameterized to reproduce solvation transfer free energies for small molecules (Sitkoff et al., 1996). The protein–water interface is determined using the protein molecular surface, which is calculated by running a 1.4-Å water probe over the van der Waals radii of the protein atoms as discussed above. The dielectric value of points inside the protein were set to 2.0, in accord with the PARSE parameterization (Sitkoff et al., 1996), and values in water were set to 80.0. In the presence of the membrane, the membrane shape determined by u+(x, y) and u−(x, y) is used to modify the local dielectric environment of the protein, and points in the membrane are given a low dielectric value of 2.0 if in the core and 80.0 if in the headgroup region. The electrostatic energy difference for insertion, ΔGelec, is then computed by subtracting the total electrostatic energy of the protein in pure solution from the value computed in the presence of the membrane. These energy terms are described in further detail in the Supplemental text, as well as in our previous work (Choe et al., 2008).
Construction of transmembrane segments
α-Helical peptides were constructed with the VMD plug-in Molefacture version 1.1 (Humphrey et al., 1996). We used SCWRL 4 (Krivov et al., 2009) to optimize side-chain rotamer conformations, and MODELLER 9v8 (Sali and Blundell, 1993) to orient the principal axes of the helix to the z axis, perpendicular to the plane of the membrane. The structures were converted to PQR format using PDB2PQR 1.4 (Dolinsky et al., 2004, 2007) with the PARSE radii parameter set (Sitkoff et al., 1994). To be consistent with most biochemical experiments, the WALP23 peptide was constructed with neutral N and C termini in PDB2PQR by turning on the neutraln and neutralc flags. All other helical segments contained charged N and C termini.
Search algorithm for identifying optimal boundary conditions
We implemented a workflow (illustrated in Fig. 2) to identify the shape of the membrane that minimizes the total energy in Eq. 1. For each cycle, we started with a given discretized contact boundary curve for the upper leaflet and the lower leaflet as shown in Fig. 1 C. We solved Eq. 3 once for each leaflet using the height boundary conditions imposed by the posited boundary curve. The imposed contact angle that the membrane makes at the protein surface was treated differently for each of the problems investigated below. For the WALP peptide, we assumed a contact angle of zero, and for the MscL channel, we assumed that the contact angle was linearly proportional to the displacement from equilibrium with a coefficient of 1. For the voltage sensor helices, we performed contact angle searches to minimize the total energy. The solution surfaces were then provided to Eq. 2 to calculate the elastostatic energy.
The membrane shape was used to create dielectric and ion-accessibility maps that were subsequently read into the Adaptive Poisson-Boltzmann Solver (APBS) version 1.2 software package to calculate the total electrostatic energy (Baker et al., 2001). Parameters for the system geometry and electrostatics calculations are shown in Table S1.
After calculating the total energy for a given membrane shape, it was used as a cost function in a Powell’s-based search (Powell, 1964) to generate a new boundary curve for the next cycle of the flowchart. Powell’s is a conjugate gradient-based method and is particularly suited to this problem because it operates well on many dimensions and does not require the calculation of derivatives. For neutral proteins, the search began from a flat membrane; however, for proteins bearing charged residues, we started the search from a bent configuration that exposes the charged residues to water. This initial guess at the boundary curve overcomes the nonpolar energy barrier associated with exposing hydrophobic residues to solvent, as shown in Fig. S1. We ended a search when the relative error was less than a small tolerance value of 5 × 10−3, following the typical stop condition for Powell’s method. In our original implementation, each boundary curve was discretized in θ, the angle around the cylindrical TM protein, and the search was performed by vertically moving nodes on the boundary curve. However, because the displacement of adjacent points was uncorrelated, the search produced high-curvature kinks in the membrane that hindered the ability to identify global minima. Therefore, we used a Fourier representation of the contact boundary:
where a0 is a constant offset, an is the amplitude of the nth mode, r0 is the radius of the TM protein, and θi is the angular position discretized with 10 points along the membrane–protein boundary. We determined that including terms above n = 4 did not improve the minimum energy configuration; however, it did increase the search space and number of iterations required to find the minimum. For n = 4, the search space is 8 + 1 for the upper leaflet and 8 + 1 for the lower leaflet, for a total dimension of 18.
Online supplemental material
In the supplemental text, we provide a more in depth discussion of how the electrostatic and nonpolar energies are calculated, as well as a table with all parameter values (Table S1). Details of three different approaches to the search algorithm are discussed along with a figure demonstrating the success of each method (Fig. S1). Fig. S2 is an updated amino acid insertion energy scale showing that the search algorithm provides very similar, but smaller, insertion energy values than those deduced in our previous study (Choe et al., 2008). This figure also demonstrates that insertion energy scales are context dependent. Finally, we provide additional information concerning the Generalized Born calculations.
Our model captures large-scale membrane rearrangements
Integral membrane proteins are characterized by stretches of hydrophobic residues that are well suited for incorporation into the low-dielectric core of the membrane. The chemical nature of these regions is important for the initial targeting of the chain to the membrane from the translocon and for the ultimate stability of the protein in the lipidic environment. However, the boundaries of TM segments are often poorly delineated because these stretches usually contain a few polar or charged residues. Therefore, it is important when modeling these systems to consider that the membrane–protein boundary is likely to have a complex shape that covers hydrophobic residues while exposing polar and charged residues to water. Moreover, the hydrophobic nature of the TM protein can induce large-scale conformational rearrangements in the bilayer when anchored to cytoskeletal elements or proteins attached to apposing membranes. For instance, SNARE, BAR, ESCRT-III, and coat proteins all cause significant deformations in the membrane (Snapp et al., 2003; Itoh et al., 2005; Lee et al., 2005; Wollert et al., 2009; Vrljic et al., 2010).
Most methods for identifying optimal protein placement and stability in the membrane assume that the membrane is static, ignoring its dynamic and flexible nature (Im et al., 2003a; Lomize et al., 2006; Ulmschneider et al., 2007). In principle, our method not only captures small deformations of the membrane that occur around buried charged residues, but it also allows for large-scale, low-energy conformations that may occur when an embedded protein is under load because of attachment to the cytoskeleton, for instance. To test our ability to identify such large distortions, we solved for the membrane shape that maximizes the stability of a WALP23 peptide (sequence GWWLALALALALALALALALWWA). WALPs are ideal for exploring membrane–protein interactions since their midsection has a strong hydrophobic signature and flanking tryptophan residues anchor the protein in the membrane by partitioning into the headgroup–water interface (Yau et al., 1998; Weiss et al., 2003).
We translated the peptide from one side of the membrane to the other, moving the center of mass (COM) from −48 to +48 Å in 3-Å steps. At each position, we used the search algorithm outlined in Fig. 2 to determine the optimal membrane configuration and corresponding energy as shown in Fig. 3. All energy differences are computed with respect to the reference state in which the membrane protein is free in solution far from an unperturbed membrane. The system energy takes on a minimum value of approximately −60 kcal/mol when the peptide COM coincides with the middle of the bilayer at z = 0. In this configuration, Fig. 3 (panel 1) shows that there is no hydrophobic mismatch because the membrane is flat, yet the hydrophobic residues (white) are maximally embedded in the membrane and the tryptophan residues (cyan) are in the headgroup region. When the peptide is moved to +24 Å (Fig. 3, panel 2), the total energy increases because the membrane bends to cover the hydrophobic residues. At +30 Å, the membrane bending energy becomes larger than the hydrophobic energy required to uncover the TM, and the bilayer exhibits a snap-through instability that extracts the TM (Fig. 3, panel 3). Our analysis shows that our search algorithm can successfully identify large-scale deflections in the membrane that bring about drastic reductions in the total energy of the system.
Predicting optimal membrane thickness for a mechanosensitive channel
Even when the hydrophobic domain of a membrane protein is clearly delineated, mismatch between the length of the hydrophobic stretch and the equilibrium length of the membrane can lead to bilayer distortions and an increase in the energy of the system, and conversely, mismatch can lead to distortions in the protein (Kleinschmidt and Tamm, 2002; Soubias et al., 2008). X-ray lamellar diffraction studies show that the membrane expands or compresses at the edge of the protein to accommodate for hydrophobic mismatch (Harroun et al., 1999), and mismatch has been shown to cause proteins to segregate to specific locations in the cell (Ravazzola and Orci, 1980). Changes in membrane thickness have also been shown to influence the functional state of membrane proteins such as mechanosensitive channels and voltage-gated ion channels (Morris et al., 2006). The mechanosensitive channel of large conductance (MscL) is a homomeric pentamer that is thought to open and close like the aperture of a camera (Sukharev et al., 2001). The application of membrane tension biases the channel from a closed state, in which the TM α-helices are primarily perpendicular to the plane of the membrane, to an expanded, open state, in which the helices are significantly tilted. This tilt decreases the hydrophobic length of the membrane-spanning region in the open state with respect to the closed state. Perozo et al. (2002) hypothesized that this structural change in the channel should lead to greater stabilization of the open state in thinner membranes compared with thicker membranes, and they verified their claim by showing that the channel open probability increased as the bilayer thickness decreased at a fixed pressure.
We used our model to determine the membrane thickness that optimally stabilizes the closed state of the MscL structure from Mycobacterium tuberculosis (Chang et al., 1998). We removed the cytoplasmic helices and embedded the TM domain in the bilayer. We then used our search algorithm to determine the optimal membrane contact curve and total insertion energy for a given membrane thickness as shown in Fig. 4 B. Next, we varied the equilibrium membrane length over the range of values suggested by the experiments performed by Perozo et al. (2002) and plotted the energy values with respect to the minimum value (Fig. 4 A). The channel is stable in the membrane over the entire range, with the optimal thickness being 38 Å. For values larger than 38 Å, the membrane thins as it approaches the channel surface. The OPM method predicts an optimal thickness that is several Ångstroms larger than our method (Lomize et al., 2006); however, our value is in better agreement with fluorescence spectroscopy studies showing that the hydrophobic thickness of the bilayer is 25 Å (Powl et al., 2005). Ignoring the headgroup region, we predict a thickness of 24 Å. Previous low-resolution models have been used to calculate the influence of membrane thickness on the open probability of the channel (Wiggins and Phillips, 2004), but this is not possible with our method because the structure of MscL in the open state is unknown.
Amino acid insertion energies are not additive
Hydrophobicity scales provide a straightforward means for assessing the stability of transmembrane proteins by adding up the individual energetic contribution from each amino acid in the transmembrane domain to determine the overall stability. However, insertion energy scales based on in vitro translation and insertion via the Sec61 translocon suggest that the apparent transfer-free energy for an amino acid depends on the amino acid sequence of the transmembrane segment (Hessa et al., 2007). Thus, it is possible that amino acid insertion energy values are not additive, which would severely limit the value of any hydrophobicity scale.
Based on molecular simulations (Dorairaj and Allen, 2007; MacCallum et al., 2011), nonadditivity has been suggested to arise because there is no additional membrane bending energy to insert a second or third charged residue after the membrane has bent to expose the first charged residue. To further examine this hypothesis, we constructed α-helical peptides containing zero to three charged residues and probed their stability in the membrane. We find that the energetic cost of inserting each additional charged amino acid is significantly less than that of inserting the first (Table 1). Whereas a peptide with a single central arginine is 10 kcal/mol less stable than one with alanine, the addition of a second arginine requires only 1 additional kcal/mol, and a third requires only an additional 2 kcal/mol. These low energies support results from recent MD simulations which found that there is essentially no additional energetic cost required to insert an arginine once the first has already formed a water defect (MacCallum et al., 2011). Experimental support for nonadditivity comes from membrane protein folding experiments, which show that the cost to insert two charged residues is less than twice the sum of inserting a single charged residue (Moon and Fleming, 2011). Nonetheless, we predict cooperativity values of ∼9 kcal/mol for inserting two arginines, whereas the predicted cooperativity based on the porin folding studies is only 1.6 kcal/mol (Moon and Fleming, 2011). These energy values are highly dependent on system geometry, and it may be difficult to compare values obtained for a single pass α-helix with those obtained from a β barrel. Our results highlight the nonadditivity inherent in this system and suggest that a simple hydrophobicity scale may lead to incorrect conclusions, especially when considering highly charged proteins or peptides.
Although our method captures the nonadditivity inherent in these systems, we wanted to compare our calculations to Generalized Born models for electrostatics in the presence of the membrane, which do not account for changes in membrane geometry (Im et al., 2003a; Lomize et al., 2006; Ulmschneider et al., 2007). To make a proper comparison with these models, we extracted the electrostatic component of the insertion energy for helices containing 0–3 arginines from Table 1 and reported these values in Table 2. Energy values are reported as ΔΔGn = ΔGn − ΔGn−1, such that the value represents the energy required to insert one more charged residue into the TM helix. Please see the Supplemental materials and methods for details concerning these calculations. Even the electrostatic component of the energy alone from our model is highly nonadditive; however, the energy values calculated using the model of Im et al. (2003a,b) is quite linear, indicating that subsequent arginines are as costly to insert as the first. Additionally, the electrostatic component of the energy for inserting even a single charged residue is 25–30 kcal/mol greater than that predicted with our model. For charged membrane proteins and membrane-associated proteins, Generalized Born models that treat the membrane dielectric as a uniform slab could give rise to incorrect results. However, for proteins that are not highly charged, such methods may be sufficient.
Some voltage sensor segments are stable in the membrane
Voltage-gated potassium (Papazian et al., 1987), sodium (Noda et al., 1984), and proton channels (Ramsey et al., 2006; Sasaki et al., 2006), as well as voltage-gated phosphatases (Murata et al., 2005), all contain 4–7 charged residues in their fourth TM segment that are critical for their ability to sense changes in membrane potential. How such highly charged segments stably incorporate into the membrane is an outstanding question in membrane biophysics, and many researchers believe that other TM segments are required for incorporation (Catterall, 1986b; Guy and Seetharamulu, 1986). In contrast, both experiment (Hessa et al., 2005b) and simulations (Freites et al., 2005; Wee et al., 2011) suggest that some S4 segments favorably adopt a transmembrane configuration.
To further explore these controversial results, we performed calculations on idealized helices with sequences corresponding to the S4 segments from the Kv1.2 Shaker-like potassium channel from rat and the KvAP archaebacterial channel. Note that the sequence of Kv1.2’s S4 segment is identical to Shaker. In the absence of membrane bending, both helical segments are highly unstable in the membrane, with transfer free energies of +72 kcal/mol and +99 kcal/mol for Kv1.2 and KvAP, respectively. Remarkably, when we allow the membrane to bend, our model predicts that both segments are quite stable in the membrane with the S4 from Kv1.2 stabilized by −31 kcal/mol and the segment from KvAP stabilized by −33 kcal/mol (Fig. 5). This drastic reduction in energy is brought about by relatively modest distortions in the membrane as can be seen from the minimum energy configurations, also pictured in the figure. It is not surprising that the S4 segment from KvAP is 2 kcal/mol more stable than the Kv1.2 segment because it has only 5 positive charges, whereas Kv1.2 has 6 charges.
With any new approach, it is useful to have a benchmark with which to validate the model. Fully atomistic MD simulations are the gold standard in this case, but detailed free energy calculations, even on a single pass TM, in the presence of a membrane are extremely demanding, and there is a large potential for sampling error. Fortunately, the Sansom group has performed free energy calculations on both of these S4 segments using a more tractable, coarse-grained model of the system (Bond et al., 2008). They report insertion energy values of −36 kcal/mol and −38 kcal/mol for the Kv1.2 and KvAP S4 segments, respectively (Wee et al., 2011), which is in strikingly good agreement with our absolute values and our predicted relative stability between both voltage-sensor segments. The Sansom group used their coarse-grained model results as a starting point to carry out fully atomistic free energy calculations on the Kv1.2 S4 segment using the GROMOS96 force field (Wee et al., 2011). This resulted in a minimum energy configuration that was −45 kcal/mol more stable than the segment in water, which is 9 kcal/mol lower than the coarse-grained model results and 14 kcal/mol more stable than our results (Wee et al., 2011). Thus, although the membrane deformations and energies in our model are similar to the deformations from coarse-grained and fully atomistic simulations, our predicted energy values are slightly higher. Although there are differences in the three energy potentials that may account for these discrepancies, we believe that the most obvious deficiencies in our model are the lack of protein tilt and side chain reorientation. Adding these two additional degrees of freedom will reduce the minimum energies that our model will produce, and hopefully bring our values into closer absolute agreement with other calculations. Additionally, as explored below, the material properties of the bilayer can significantly impact the insertion energy of the helix, and it is not clear what the bilayer parameters in our model should be to most closely approximate the properties of the coarse-grained membrane.
It is thought that the RXXR spacing of residues in many S4 segments is important to their stability. To examine this, we created 18 mutated Kv1.2 sequences preserving the total charge but disrupting the spacing, and we found every mutation led to a higher insertion energy value. Interestingly, the S4 segment from the hyperpolarization-activated potassium channel KAT1 has an unconventional charge spacing with two adjacent arginines, and there is experimental evidence that this segment will not insert into the membrane when isolated from the rest of the channel (Zhang et al., 2007). We applied our model to the KAT1 S4 segment and observed a stabilizing transfer free energy of only −24 kcal/mol, which is 7 kcal/mol less stable than either the Kv1.2 or KvAP S4 segments. Therefore, we believe that this lower insertion energy may be related to the charge spacing on KAT1; however, our results still predict that the S4 from KAT1 should be stable in the membrane, which is at odds with the finding that the S3 segment is also required for membrane stabilization (Zhang et al., 2007). To specifically explore the importance of charge spacing, we systematically varied this spacing and measured the stability of single pass TM segments. Our results suggest that charged residues spaced two apart, for example RXXR, are 2–4 kcal/mol more stable than those spaced by 0 or 1 uncharged residues (Table 1). Visualization of the minimum energy configuration for each case shows that the membrane need only bend on one side of the helix to expose charged residues separated by two intervening nonpolar residues, whereas the membrane undergoes much more extensive distortions to expose charged residues with different spacings (unpublished data).
Recent x-ray structures of voltage-gated ion channels suggest that the S4 helix exists, at least partially, in a 310 helix rather than an ideal α-helix (Long et al., 2007; Clayton et al., 2008; Vieira-Pires and Morais-Cabral, 2010). The 310 configuration places all the charges with an RXXR spacing on one face, which localizes the membrane bending to one side of the helix and may reduce the bending energy. We explored this possibility by creating 310 helices of the two sequences in Table 1 that have RXXR motifs, one harbors 2 arginines (indicated by * in the table) and the second harbors 3 arginines (indicated by ** in the table). The insertion energy values are higher when these sequences adopt a 310 conformation versus an α-helical conformation by +5.3 kcal/mol (sequence with 2 arginine) and +5.1 kcal/mol (sequence with 3 arginine). Our calculations indicate that the α-helix configuration is slightly more stable because of an increased nonpolar stabilization; the α-helix is more compact and buries more surface area in the membrane. Thus, we believe that the propensity for portions of the S4 helix to adopt a 310 configuration may be determined by local interactions with the rest of the channel rather than energetic interactions with the lipid membrane.
In our previous manuscript, we showed that the membrane dipole only moderately influenced membrane protein stability because the majority of the amino acids between the upper and lower leaflets were hydrophobic and therefore neutral (Choe et al., 2008). However, this was the case for helices harboring a single charged residue, and it is possible that with many basic residues more charge becomes buried between the leaflets and experiences the significant positive electrostatic potential created by the lipid headgroups. To explore this scenario, we added a thin layer of positive and negative charge to the upper and lower leaflets, with the negative layer closer to the membrane–water interface, as described in our previous work (Choe et al., 2008). The net charge sums to zero, and the separation length and charge values were chosen to create an interior membrane potential of +300 mV, which is near the peak dipole potential value of +260 mV measured for phosphatidylcholine headgroups with ester linkages to the tail (Wang et al., 2006). For the configurations pictured in Fig. 5, the dipole potential destabilizes the TM helix by 2.6, 2.3, and 1.6 kcal/mol for the KvAP, Kv1.2, and KAT1 S4 segments, respectively. Although larger than previously observed for helices with a single arginine, the destabilization is not very large, because the membrane bends to keep most charged groups out of the core.
Membrane protein stability depends on bilayer stiffness
Studies on outer membrane proteins have shown that the elastic properties of the membrane influence protein stability in the membrane (Kleinschmidt and Tamm, 2002; Hong and Tamm, 2004). We therefore investigated the effect of increasing bilayer stiffness on the insertion energetics by varying the bilayer compression modulus (Ka). In all calculations, we used the Ka value of 142.5 pN/nm measured experimentally by White (1978) and employed by Nielsen et al. (1998) in their mattress models. However, this value is at the low end of the physiological range, and in the presence of cholesterol Ka can be as high as 1,200 pN/nm (Tierney et al., 2005). We used this later value as an upper value, and we calculated the insertion energies for each of the voltage sensor S4 helices over the entire range, as shown in Fig. 6. In each case, we performed the full search procedure to identify the optimal shape that minimizes insertion. As Ka increases, the total insertion energy values increase significantly due to the increase in the membrane deformation energy. Importantly, even moderate increases in the compression modulus destabilize the KAT1 S4 segment, whereas segments from KvAP and Kv1.2 remain stable. This observation is in very good agreement with the experimental observation that the S4 from KAT1 is not stable in the membrane alone (Zhang et al., 2007), but that the isolated voltage sensor helix from KvAP readily inserts into the membrane (Hessa et al., 2005b).
We used our fast continuum method for determining the insertion energy of membrane proteins to explore several outstanding questions in membrane protein biophysics. Linking the three numeric solvers of our method and adding the search algorithm permits us to determine arbitrary distortions in the membrane, which is essential for understanding the true energetics of embedded proteins. Although some implicit membrane models account for membrane flexibility (Tang et al., 2006; Zhou et al., 2010), many treat it as a rigid slab (Im et al., 2003a; Lomize et al., 2006; Ulmschneider et al., 2007). As shown in Fig. 3, our algorithm readily identifies the putative membrane-spanning region of membrane proteins and moves to bury the hydrophobic residues. As the protein is translated away from the center of the bilayer, the membrane undergoes a large, low energy conformational change to minimize the energy of the system. We believe that such distortions are important for understanding the shapes and energies of membrane proteins attached to cytoskeletal and extracellular elements such as integrins and cadherins. Similarly, the interaction between the membrane and the protein has been shown to regulate the function of some proteins, such as the stretch-activated channel MscL. By systematically varying the equilibrium length of the membrane, we predict an optimal equilibrium value that minimizes the total system energy, which includes not only strain in the membrane but also electrostatics and nonpolar effects caused by hydrophobic mismatch (Fig. 4). Based on these calculations, we predict membrane thinning at the edge of the channel for a range of equilibrium thicknesses, and the extent of this compression is in good agreement with experiments (Powl et al., 2005).
Our model suggests that amino acid insertion energies are nonadditive. This is most important when considering the placement of multiple charged residues into the TM domain as we show that placing a second arginine into a TM already containing one may cost as little as 1 kcal/mol. In part because of this effect, we show that S4 voltage sensor segments can be stable in the membrane, as already suggested by experiment (Hessa et al., 2005b), qualitative MD simulations (Freites et al., 2005), and quantitative free energy calculations (Wee et al., 2011). Our method provides a simple mechanical explanation for this nonadditivity—once the membrane bends to accommodate one charged residue it no longer needs to bend for the next one. This feature is an integral component of our method, but it is missing from other implicit membrane models (Im et al., 2003a; Lomize et al., 2006; Ulmschneider et al., 2007), as we explicitly demonstrate in Table 2. However, our method still predicts insertion energies for single charged amino acids that are 4–8 kcal/mol larger than those predicted by the translocon scale (Hessa et al., 2005a), and 6 kcal/mol larger for arginine compared with the porin folding scale (Moon and Fleming, 2011), but a nearly identical value for lysine compared with the porin folding scale (Moon and Fleming, 2011; Fig. S2). In general, our larger values may result from limitations of our system discussed below; however, there are open questions concerning the interpretation of the translocon studies, including whether the H-segment is actually centered in the membrane (Dorairaj and Allen, 2007) and the role of the two additional TM segments that may alter the stability of the central residues (Shental-Bechor et al., 2006).
Previously, we determined that inserting the KvAP S4-S3 helix–turn–helix motif into a flat membrane is energetically unfavorable (Grabe et al., 2004), but here we show that the S4 helix is stable if the membrane is deformable. We incorporated membrane bending into our solvation model by using classical elastostatics to describe the equilibrium shape of the membrane initially proposed by Helfrich (1973) and expanded to include mean bending, bilayer compression, and surface tension by Nielsen et al. (1998). Recently, a multiscale modeling approach was developed that used a similar continuum model of the membrane to quantify the energetics of membrane deformations observed in fully atomistic membrane protein simulations (Mondal et al., 2011). This model used a simple finite difference method on a square grid to solve for the membrane shape, and with a relatively fine spatial grid, they were able to compute shapes and energies for noncylindrical G protein–coupled receptors. As discussed below, we intend to develop a more general membrane model that can handle arbitrarily shaped membrane proteins. Even when proteins are hydrophobically matched to the width of the membrane, it has been shown that noncylindrically shaped proteins can induce strain in the membrane (Kim et al., 1998), and this strain can influence protein function when bilayer leaflets contain lipids favoring spontaneous curvature (Dan and Safran, 1998). These considerations can be incorporated into our model through modifications to the membrane energy density in Eq. 2, as detailed by Dan and Safran (1998). At the same time, we would like to incorporate more atomistic detail of the lipid headgroups into our membrane model. For instance, charged lipids may unevenly accumulate in certain regions near the embedded protein, and the Weinstein laboratory has developed a mean field theory for dealing with this phenomena (Khelashvili et al., 2009), which is well suited for our continuum method. Additionally, we could assign spatially dependent material properties that could attempt to account for different lipid densities or types in the upper versus the lower leaflet, which has been shown to affect opening of the MscL channel using simulation (Jeon and Voth, 2008).
In summary, we have a fast, predictive method for determining the stability of proteins in the membrane that does not rely on scales or the assumption of additivity. Our results agree well with coarse-grained models and fully atomistic simulations; however, we estimate our method to be roughly 600 times faster than comparable coarse-grained MD calculations (Wee et al., 2011) and 40,000 times faster than fully atomistic calculations (Dorairaj and Allen, 2007). In the near future, we intend to include three additional features into our model. First, we intend to adopt a 3-dimensional model of the membrane that explicitly represents the strain field between the upper and lower leaflets, and we will employ edge detection to accurately represent the protein–membrane boundary. This will allow us to incorporate protein tilt into our model, which has been shown to be important for single TM WALP peptides (van der Wel et al., 2002; Strandberg et al., 2004; Ozdirekcan et al., 2005), and it will also allow us to investigate proteins of arbitrary shape. Second, MD simulations have shown that charged residue side chains “snorkel” to interact with the polar headgroups and solvent (Dorairaj and Allen, 2007). To account for these changes, which can impact insertion energies, we will incorporate a rotamer library search on the charged residues to minimize membrane distortions and electrostatic penalties. Third, we use a simple SASA model for the nonpolar energy. We will explore different continuum models for this energy that are more specific for water–membrane transfer free energies, as well as those that account for dispersive solute–solvent interactions (Wagoner and Baker, 2006). Finally, we intend to integrate the membrane deformation algorithm presented in this manuscript into our software APBSmem (Callenberg et al., 2010) to help orient proteins in the membrane for use in interpreting experiments and carrying out MD simulations in a similar manner to the Orientations of Proteins in Membranes database (Lomize et al., 2006).
We would like to thank Charles Wolgemeth (University of Connecticut), Yann Chemla (University of Illinois, Urbana Champaign), and Sukrit Suksombat (University of Illinois, Urbana Champaign) for fruitful discussions, and Patrick van der Wel (University of Pittsburgh) for helpful comments on the manuscript.
This work was supported by NSF CAREER award MCB0845286 to M. Grabe and a Howard Hughes Medical Institute Undergraduate Research Fellowship to N.R. Latorraca.
Kenton J. Swartz served as editor.