Membrane environment modulates the pKa values of transmembrane helices.

In this work, we apply the recently developed constant pH molecular dynamics technique to study protonation equilibria of titratable side chains in the context of simple transmembrane (TM) helices and explore the effect of pH on their configurations in membrane bilayers. We observe that, despite a significant shift toward neutral states, considerable population of different side chains stay in the charged state that give rise to pKa values around 9.6 for Asp and Glu and 4.5 to 6 for His and Lys side chains, respectively. These charged states are highly stabilized by favorable interactions between head groups, water molecules, and the charged side chains that are facilitated by substantial changes in the configuration of the peptides. The pH dependent configurations and the measured pKa values are in good agreement with relatively recent solid state NMR measurements. Our results presented here demonstrate that all-atom constant pH molecular dynamics can be applied to membrane proteins and peptides to obtain reliable pKa values and pH dependent behavior for these systems.

Introduction

Membrane proteins play a significant role in myriad biological processes and account for ∼30% of all proteins in the cell.[1−3] The transmembrane (TM) segments of these proteins are usually composed of hydrophobic helices whose sequence determines their orientation and position of in membranes.[4−6] Despite the hydrophobic nature of the TM segments, occurrences of charged residues are believed to be important for the function of membrane proteins. In integrins, for instance, basic amino acid side chains that are buried in the membrane regulate transmembrane signaling,[7] and the ion selectivity of nicotinic-type receptors is determined by the presence or absence of a pore facing carboxylate ring in the membrane bilayer.[8,9] Recent studies suggest that the microenvironment of these buried residues, for example, lipids and phosphate head groups, has a profound effect on their charged states, which would influence the orientation of the titratable side chains and their helical configurations in the membrane.[6,8,10−12] Therefore, detailed information about the pKa values of buried residues in hydrophobic environments can shed light on the extent of microenvironment effects on the population of the charged and neutral states and, thus, pH-modulated biological function of membrane proteins.

The role of buried ionizable residues in determining helix orientation and the effect of the membrane bilayers in modulating their pKa values have been investigated by solid state NMR (SSNMR) in simple model peptides such as the WALP[13] and GWALP[14] series using quadrupolar splitting of labeled Ala residues. These TM peptides have provided the opportunity to capture peptide–membrane interactions in great detail and predict possible scenarios for protein–membrane interactions in more complicated systems.[10,11,15−18] Vostrikov and Sansom et al. have shown that buried Arg residues at different positions in the GWALP23 peptides gives rise to different peptide behavior in membrane bilayers. The GWALP23 peptide with an Arg in position 14 has an average tilt angle of 16.2° which is almost 10° greater than the GWALP23. In addition, the GWALP23 peptides with one Arg in position 12 tend to exit the lipid bilayer and do not form stable TM orientations in the membrane.[11] Using the same peptide models, Koeppe et al. revealed that Arg maintains its positive charge in the membrane bilayer while Lys containing peptides titrate and change their tilt angles. For instance, the pKa of a Lys residue in the 14th position of Lys-containing GWALP23 peptides is shifted down by ∼4 pH units (6.2 at 323 K and an estimated value of 6.8 at 298 K) compared to the standard pKa of Lys in aqueous solution (10.4). These peptides also adopt a tilt angle of 15°, which is ∼10° larger than the tilt angle of the GWALP23 peptide.[10] Similar to Arg-containing peptides,[11] the GWALP23 with Lys in the 12th position shows a smaller tendency to remain in a single membrane bound configuration, which makes it difficult to determine the pKa of this residue.[10]

These examples and several more[8,19−22] are strong evidence of the effect of the microenvironment on the protonation equilibria of protein side chains. Thus, along with experimental techniques such as NMR, several computational methods have been developed to calculate the pKa values of titratable residues in different environments. Among these techniques, constant pH molecular dynamics (CPHMD) approaches are particularly devised to study pH dependent behavior of proteins for which little information about the charged states of the key residues are available. Unlike Poisson–Boltzmann based approaches,[23] in the CPHMD methods, protonation states of a side chain are coupled with the protein dynamics and are allowed to propagate with time following degrees of freedom in the system. There are two distinct CPHMD approaches that are in use. In discrete CPHMD methods, first developed by van Gunsteren et al.,[24−26] the charge states of titratable residues are determined by performing Monte Carlo (MC) steps at specific intervals of regular molecular dynamics simulations. In this approach, the abrupt change in the protonation states of the residues (and overall charge of the system) can lead to discontinuity of the potential energy and possible side effects after MC moves along the titration coordinate. This issue is alleviated in the continuous CPHMD methods pioneered by Brooks et al. In this category of CPHMD techniques, the titration coordinates propagate in time using a λ dynamics approach,[27−29] which provides a continuous switch between different protonation states and prevents sudden jumps in the potential energy and its derivatives. We note that λ dynamics is related in its underlying formulation to techniques devised by Pettitt and co-workers for the addition and deletion of particles within a grand canonical ensemble to achieve constant chemical potential.[30,31]

First implemented in an implicit solvent framework[32] and later improved to include proton tautomerism, continuous CPHMD[33] has been shown to be successful in calculating the pKa values of titratable residues for a diverse set of proteins.[33−38] In spite of numerous successes, this method suffers from the well-known disadvantages of the generalized Born implicit solvent models. For instance, the overestimation of the desolvation energy for the buried charges in the protein core,[39,40] or the compactness of the protein structures,[38,41] can induce systematic errors in pKa calculations. More significant errors can occur for membrane proteins with buried charge residues as a result of assumptions used in the implicit membrane models, such as the rigid nature of the membrane bilayer,[42−44] the insertion of a charged residue in the hydrophobic core of the membrane requires complete loss of its solvation shell. However, as evidenced in detailed explicit membrane/solvent simulations, charged residues in the membrane are surrounded by a few water molecules and possibly phosphate head groups.[45,46] Partial solvation of the charged residues helps to compensate for the low dielectric environment of the membrane and significantly reduces the insertion energy, while total loss of the solvation shell in implicit membrane models leads to systematic overestimation of side chain pKa’s.[47,48] Recently, Brooks and co-workers developed a novel explicit solvent CPHMD approach based on a new multisite λ dynamics framework (CPHMDMSλD)[49,50] and applied it to study the pH-dependent dynamics of several RNA[51,52] and protein[53] systems. In addition, to address the convergence problems of earlier CPHMD variants,[51,54] pH-based replica exchange may be employed to enhance the sampling of the spatial and titration coordinates.[52]

Inspired by the recent pKa measurements of the buried Lys residues in 1,2-di(9Z-octadecenoyl)-sn-glycero-3-phosphocholine (DOPC) bilayers,[10] we were motivated to utilize the CPHMDMSλD method on several TM-GWALP23 peptides containing titratable residues at position 12 or 14 to examine the influence of the membrane on the structure and pKa of the ionizable residues in these peptides. Our results indicate that the CPHMDMSλD method reproduces the pKa’s of Lys residues in the Lys-containing GWALP23 peptides inserted in the membrane bilayer, in good accord with the pKa measurements from quadrupolar splitting (QS) of deuterated Ala residues.[10] In addition, we observe that the peptide tilt angle is affected by the protonation state of the titratable residues: for example, the peptides are significantly more tilted in charged states, and this general trend of the tilt angle variation is in agreement with the SSNMR data. In addition, we directly predict the pKa variations of His in an ether-linked DOPC bilayer and Glu and Asp in an ether-linked 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC) in the context of TM-GWALP23 peptides. Overall, our results indicate that both protonation states and tilt angles in the peptides containing titratable residues are considerably influenced by their surrounding environment.

Methods

The list of the peptides simulated in this paper along with their sequence and the membrane they are inserted in are provided in Table 1. Insertion of the peptides along with the equilibration steps are described in more detail in the Supporting Information.

Table 1

Sequences, Lipid Bilayer Composition, Number of Replicas, pH range, and Calculated and Experimental Values for Each Simulated TM Peptide

name	sequence	membrane	no. of replicas (pH range)	pKa
K12	acetyl-GGALY(LA)3K(AL)3WLAGA-amide	DOPC	8(2–9)	5.7 ± 0.2(<7)b
K14	acetyl-GGALY(LA)4K(AL)2WLAGA-amide	DOPC	8(2–9)	5.6 ± 0.2(6.8)b
H12	acetyl-GGALY(LA)3H(AL)3WLAGA-amide	DOPCa	9(0–8)	4.0 ± 0.1
H14	acetyl-GGALY(LA)4H(AL)2WLAGA-amide	DOPCa	9(0–8)	4.5 ± 0.3
E12	acetyl-GGALY(LA)3E(AL)3WLAGA-amide	DLPCa	9(5–13)	9.8 ± 0.1
E14	acetyl-GGALY(LA)4E(AL)2WLAGA-amide	DLPCa	9(5–13)	9.7 ± 0.1
D12	acetyl-GGALY(LA)3D(AL)3WLAGA-amide	DLPCa	9(5–13)	9.4 ± 0.2
D14	acetyl-GGALY(LA)4D(AL)2WLAGA-amide	DLPCa	9(5–13)	9.4 ± 0.2

Carboxyl group is substituted with ether.

Experimental values in parentheses taken from ref (10).

The pH-replica exchange simulations (pH-REX) on the membrane-inserted peptides were performed in the following manner. The structures obtained from the equilibration runs of protonated and unprotonated side chains were used as the initial structures of the pH-REX calculations with each equilibrated structure used as alternative replicas. The number of replicas and the pH range for each system is listed in Table 1. Replica exchange in the pH variable has previously been shown to improve the convergence of pKa calculations.[55] All simulations were performed using CPHMDMSλD with the BLOCK facility and the λ functional form with a coefficient of 5.5.[53] Expressing λ in this particular form in terms of virtual particles that propagate through the course of the simulation (θα) has been shown to increase energy stability. More detailed discussion on the effect of functional form of λ versus θα is provided in ref.[49] All nonbonded energy terms were linearly scaled by λ while all bond, angle and dihedral terms were unscaled. λ values were saved every 10 steps. A nonbonded cutoff of 12 Å was used with an electrostatic force-shifting function and a van der Waals switching function. A fictitious mass of 12 amu·Å2 was assigned to all θα. In the pH-REX simulations, exchange attempts were made every 1 ps. Each pH-REX simulation was repeated 3 times for 5 ns per replica, from which the last 3 ns was used for analysis. The average pKa’s are reported in Table 1. All parameters for the pH-Rex simulations were kept as default values as listed in ref (53).

The fraction of unprotonated state at each pH (Sunprot) was determined from the total population of the unprotonated (Nunprot) and protonated states (Nprot) states.These populations are defined as the total number of times in the trajectory that the condition λ > 0.8 for eitnher the unprotonated or protonated state is satisfied.

To obtain the pKa of the system of interest, Sunprot is fitted versus pH using the Henderson–Hasselbalch equation.

The configuration of the peptides at different pH values was also examined. The tilt angle of the peptides was measured as the angle between the membrane normal (z-axis) and the peptide axis, which was defined as the vector connecting the Cα atoms. The rotation angle is defined as the angle between the vector that connects the helical axis to the Cα atom of residue Leu 14 and is perpendicular to the helix axis, and the projection vector of the z-axis onto the plane made by the second and third principal axes. For a more detailed definition of the tilt and rotation angles, please see ref (56).

Results and Discussion

We have simulated eight TM peptides, each having one titratable residue (i.e., K, H, E, or D) at position 12 or 14 in the peptide. The purpose of this study is to investigate the effect of the membrane bilayers on the pKa values of these residues, as well as to examine the pH-dependent variation of the peptide configurations in the membranes. In addition, we compare the pKa values obtained from the CPHMDMSλD simulations with those of Lys containing residues from the SSNMR measurements[10] and make direct predictions for the other titratable side chains.

The pKa values for different side chains, averaged over three runs, are listed in Table 1. The corresponding titration curves are displayed in Figure 1. Our results show that the Lys residue in the K14 peptide has a pKa of 5.6 ± 0.2, in good agreement with the estimated experimental value of 6.8 at 300 K[10] and within the error bars observed for insertion free energies of different side chains in membrane bilayers.[57] The difference we see here may be due to the steep gradient of the membrane dielectric across its normal. Small changes in the position of the side chain across the membrane normal is accompanied by large variation in the apparent dielectric constant of the membrane, which can cause significant changes in the side chains pKa.[42,58] The QS measurements of Ala residues for the K12 peptide only determined an upper bound of 7 for the pKa value of the Lys side chain.[10] Our calculated pKa of 5.7 ± 0.2 is consistent with this bound. Thus, we find that, consistent with the SSNMR experiment, the pKa of Lys residues at the 12th or 14th position of GWALP peptides are shifted by ∼4.4 pH units. Interestingly, similar pKa shifts were previously observed for Lys residues in engineered mutants of staphylococcal nuclease.[21] Among 25 single mutations, 4 mutants with Lys residues located in the inner core of the protein showed depressed pKa values between 5.3 and 6.5,[21] which were accompanied by partial or global unfolding of the protein at low pH values.

Figure 1

Titration curves for three trials of pH-REX for different peptides listed in Table 1. The unprotonated fractions obtained from pH-REX are shown with circles, and the lines represent the Henderson–Hasselbach fit. In some panels, due to overlap of the curves and data points, only two out of three colors are visible.

Since no measurements yet exist for substitutions of His, Asp, or Glu at these same positions in the GWALP peptides, we decided to make direct predictions of the pKa values and the corresponding peptide configurations. For all substituted residues, significant pKa shifts with respect to the standard reference pKa values are observed. In the case of the His-containing peptides, H12 and H14 in Table 1, pKa shifts of ∼2 and ∼1.5 units toward the neutral form are observed, respectively. Glu residues at these positions (E12 and E14) yield calculated pKa values of 9.8 ± 0.1 and 9.7 ± 0.1. For the Asp residue at position 12 (D12) and 14 (D14), pKa values of 9.4 ± 0.2 are obtained. We note that in a different study using fluorescence emission of a Trp residue, a pKa of 8.7 was reported for an Asp residue in the transmembrane helix formed by the peptide K2GL7DLWL9K2A (termed pL(D11)) indicating a ∼4.7 pH unit shift in favor of the neutral state in a DOPC bilayer.[12] While, due to different peptide sequences and membrane compositions, direct comparison between this value and the pKa observed in our simulations is difficult, the difference can be attributed to the interfacial orientation of the peptide that is suggested to be adopted in the latter study.[12] The significant shift of the pKa values of E and D at positions 12 and 14 suggest that the ratio of protonated versus unprotonated side chains in membranes is severely shifted toward the latter and, perhaps, in the absence of any other stabilizing interactions (e.g., salt bridges with neighboring residues), these side chains tend to stay protonated in the membrane within the biologically relevant pH ranges.

In order to investigate the effect of pH on the charge state and eventually the configurations of these peptides in membrane bilayers, we examined the tilt and rotation angles of the peptides at different pH values. The potential of mean forces (PMF) for basic (K and H) and acidic (E and D) side chains along their tilt and rotation angles were calculated and are shown in Figures 2 and 3, and the dominant structures extracted from the replicas with specified pH values are depicted in Figures 4 and 5, respectively. For positively charged side chains, depicted in the left column of Figure 2, the peptides adopt tilt angles between 30° and 50°, and rotation angles between 20° and 100°. The neutral peptides, which are populated at high pH values, on the other hand adopt smaller tilt angles roughly form 5° to 25° with rotation angles that are not significantly affected (Figure 2 right column). The experimental values for the K14 peptide at high and low pH values and for K12 at high pH obtained from the QS of deuterated Ala residues[10] are marked in the figure with red circles. Due to complexity of the QS data for K12 at low pH, no tilt and rotation angles were reported. It can be seen that the tilt and rotation angles obtained experimentally for K12 and K14 at high pH values are part of diverse ensembles that are sampled by these peptides in our simulations. For positively charged K14, on the other hand, the tilt and rotation angles that are adopted by the peptide in our simulations are significantly larger.

Figure 2

2D PMFs of rotation angle vs tilt angle for TM helices containing K and H side chains for low and high pH values where the titratable side chains are predominantly protonated (positively charged) and unprotonated (neutral), respectively. The experimental values are marked with red circles.[10] The color bars are in kcal/mol.

Figure 3

2D PMFs of rotation angle vs tilt angle for TM helices containing E and D side chains for low and high pH values where the titratable side chains are predominantly protonated (neutral) and unprotonated (negatively charged) respectively. The color bars are in kcal/mol.

Figure 4

Representative structures of TM peptides containing one basic side chain obtained from the replica with specified pH value. The structures were rendered using VMD.[63] The peptide backbone is depicted with green cartoon representation, while the titratable side chains are shown with blue van der Waals representations. The bilayer phosphate head groups are shown with orange spheres. The lipid tails are not shown to avoid complexity. Water molecules whose center of masses are within 5 Å of the side chains are illustrated with van der Waals spheres and red color.

Figure 5

Representative structures of peptides obtained from the replica with the specified pH. The color coding and representation is the same as described in Figure 4. The structures were rendered using VMD.[63]

We also examined the distribution of tilt and rotation angles for the acidic side chains at both high and low pH (see Figure 3). At low pH, where the side chains are protonated and neutral (Figure 3 left column), the peptides tend to have smaller tilt and rotation angles compared to high pH (Figure 3 right column), where the side chains are negatively charged.

While the pH dependent variations of the tilt and rotation angles of K12 and K14 agree well with the experimental values, the values we observe are larger than those reported from the experimental measurements. Similar apparent disagreements between tilt and rotation angles obtained from molecular dynamics simulations and the values calculated from SSNMR observables have been noticed before for several TM peptides including the WALP series.[17,18,43,59] These discrepancies are believed to stem from orientational fluctuations of the peptides in membrane bilayers. Recently, extensive attempts have been made to capture the complex motions of TM helices in different models in order to extract more representative tilt and rotation angles from SSNMR QS measurements. However, SSNMR observables are time and ensemble averaged data, which may not fully represent a single helix orientation. Therefore, more complicated methods are needed to obtain the tilt and rotation angles from these measurements.[60]

Finally, the effect of pH on the orientation of the peptides and the structure of the membrane aqueous interface can be highlighted with further details by looking at representative structures from our simulations. The dominant structures of each peptide at high and low pH values were selected as noted in the computational methodology and the results for positive and negative side chains are depicted in Figures 4 and 5, respectively. For the basic side chains, as illustrated in Figure 4, at low pH (left column), the peptides tend to have highly tilted structures with more interfacial exposure of the side chain. The tilting of the peptides accompanied by the deformation of the bilayer, which allows solvation of the charged state with water molecules and head-groups, leads to considerable stabilization of the charged states.[45,46] However, increasing the pH (Figure 4 right column) causes the solvation shell to be lost and the tilt angle to decrease. The same effect can be observed in Figure 5 for acidic side chains. At low pH values, the peptides tend to have smaller tilt angles while at high pH the tilt angles are increased and the solvation shell around the side chains is formed (compare left and right columns of Figure 5). In addition to the improved charge solvation provided by helix tilting, the enhanced tilting of the peptide in the presence of the charged side chain can also be explained by the effect of partial solvation on the membrane thickness. The water molecules and head groups in the solvation shell that surround the charged side chain and penetrate into the hydrophobic core of the membrane, cause the local membrane thickness around the side chain to decrease. The thinning of the membrane causes the helix to tilt to maximize the interactions between membrane core and the hydrophobic peptide. In addition, tilting of the helix when the side chain is charged provides more exposure of the side chain to the membrane exterior where it incurs a lesser desolvation penalty. The effect of peptide sequence and membrane thickness fluctuation on the peptide tilt angle have been discussed in more detail in the following references: (45, 47, 48, and 61) and the references therein.

Our simulations suggest that the observed pKa values are tightly coupled to the conformation of the peptide and membrane/water interface. Detailed explicit solvent simulations of insertion of proteins and small molecules into lipid bilayer[57,62] have revealed that these fluctuations are typically kinetically hindered and long time scale simulations (possibly ms) may be required to ensure convergence of the thermodynamics properties. This may be reflected in the pKa values reported here. In principle, CPHMDMSλD simulations can be combined with umbrella sampling along rotation and tilt angles of the peptide or advanced enhanced sampling techniques such as the ensemble dynamics approach[56,60] to strengthen the sampling along these two reaction coordinates. However, our generally good agreement with experiment in the case of the lysine based peptides and the similar charge-state dependence of the peptide orientation and membrane interface structure for both basic and acidic peptides suggest that our findings are at least semiquantitatively reasonable.

Conclusions

In this study, we have applied for the first time constant pH molecular dynamics techniques to calculate the pKa values of titratable side chains in the context of TM helices and investigated the effect of pH on the configuration of the studied peptides. Our calculated pKa values for the Lys-containing peptides and the values previously obtained using SSNMR QS of deutrated Ala residues[10] are in good agreement. We also observed that the configuration of the various peptides and the arrangement of the lipid molecules in the vicinity of the side chains are coupled to the external pH. At pH values where the population of the charged states dominates, the peptides tend to have more tilted structures with notable membrane deformation. On the other hand, when the pH is such that the population of the neutral side chains is pronounced, the tilt angles are decreased and the membrane deformation is less significant.

The overall good agreement between the calculated pKa values and the experimental values especially for the case of Lys containing peptides[10] and the variation of the peptide conformation with pH suggest that the CPHMDMSλD captures the essence of pH dependency of peptide membrane interactions.