Article Information

• PDF (530 kb)
• Full Text with Large Figures
• Cited by in Scopus (13)

PubMed

• Articles by Lintao Bu
• Articles by Charles L. Brooks

Related Articles

• Structure of the 1–36 Amino-Terminal Fragment of Human Phosp...

  Structure of the 1–36 Amino-Terminal Fragment of Human Phospholamban by Nuclear Magnetic Resonance and Modeling of the Phospholamban Pentamer Biophysical Journal, Volume 76, Issue 4, 1 April 1999, Pages 1784-1795 Piero Pollesello, Arto Annila and Martti Ovaska Abstract The structure of a 36-amino-acid-long amino-terminal fragment of phospholamban (phospholamban[1–36]) in aqueous solution containing 30% trifluoroethanol was determined by nuclear magnetic resonance. The peptide, which comprises the cytoplasmic domain and six residues of the transmembrane domain of phospholamban, assumes a conformation characterized by two -helices connected by a turn. The residues of the turn are Ile18, Glu19, Met20, and Pro21, which are adjacent to the two phosphorylation sites Ser16 and Thr17. The proline is in a conformation. The helix comprising amino acids 22–36 is well determined (the root mean square deviation for the backbone atoms, calculated for a family of 18 nuclear magnetic resonance structures is 0.57Å). Recently, two molecular models of the transmembrane domain of phospholamban were proposed in which a symmetric homopentamer is composed of a left-handed coiled coil of -helices. The two models differ by the relative orientation of the helices. The model proposed by Simmerman et al. (H.K. Simmerman, Y.M. Kobayashi, J.M. Autry, and L.R. Jones, 1996, J. Biol. Chem. 271:5941–5946), in which the coiled coil is stabilized by a leucine-isoleucine zipper, is similar to the transmembrane pentamer structure of the cartilage oligomeric membrane protein determined recently by x-ray (V. Malashkevich, R. Kammerer, V. Efimov, T. Schulthess, and J. Engel, 1996, Science 274:761–765). In the model proposed by Adams et al. (P.D. Adams, I.T. Arkin, D.M. Engelman, and A.T. Brunger, 1995, Nature Struct. Biol. 2:154–162), the helices in the coiled coil have a different relative orientation, i.e., are rotated clockwise by ∼50°. It was possible to overlap and connect the structure of phospholamban[1–36] derived in the present study to the two transmembrane pentamer models proposed. In this way two models of the whole phospholamban in its pentameric form were generated. When our structure was connected to the leucine-isoleucine zipper model, the inner side of the cytoplasmic domain of the pentamer (where the helices face one another) was lined by polar residues (Gln23, Gln26, and Asn30), whereas the five Arg25 side chains were on the outer side. On the contrary, when our structure was connected to the other transmembrane model, in the inner side of the cytoplasmic domain of the pentamer, the five Arg25 residues formed a highly charged cluster. Abstract | Full Text | PDF (609 kb)

• A Fluorescence Energy Transfer Method for Analyzing Protein ...

  A Fluorescence Energy Transfer Method for Analyzing Protein Oligomeric Structure: Application to Phospholamban Biophysical Journal, Volume 76, Issue 5, 1 May 1999, Pages 2587-2599 Ming Li, Laxma G. Reddy, Roberta Bennett, Norberto D. Silva, Larry R. Jones and David D. Thomas Abstract We have developed a method using fluorescence energy transfer (FET) to analyze protein oligomeric structure. Two populations of a protein are labeled with fluorescent donor and acceptor, respectively, then mixed at a defined donor/acceptor ratio. A theoretical simulation, assuming random mixing and association among protein subunits in a ring-shaped homo-oligomer, was used to determine the dependence of FET on the number of subunits, the distance between labeled sites on different subunits, and the fraction of subunits remaining monomeric. By measuring FET as a function of the donor/acceptor ratio, the above parameters of the oligomeric structure can be resolved over a substantial range of their values. We used this approach to investigate the oligomeric structure of phospholamban (PLB), a 52-amino acid protein in cardiac sarcoplasmic reticulum (SR). Phosphorylation of PLB regulates the SR Ca-ATPase. Because PLB exists primarily as a homopentamer on sodium dodecyl sulfate polyacrylamide gel electrophoresis, it has been proposed that the pentameric structure of PLB is important for its regulatory function. However, this hypothesis must be tested by determining directly the oligomeric structure of PLB in the lipid membrane. To accomplish this goal, PLB was labeled at Lys-3 in the cytoplasmic domain, with two different amine-reactive donor/acceptor pairs, which gave very similar FET results. In detergent solutions, FET was not observed unless the sample was first boiled to facilitate subunit mixing. In lipid bilayers, FET was observed at 25°C without boiling, indicating a dynamic equilibrium among PLB subunits in the membrane. Analysis of the FET data indicated that the dye-labeled PLB is predominantly in oligomers having at least 8 subunits, that 7–23% of the PLB subunits are monomeric, and that the distance between dyes on adjacent PLB subunits is about 10Å. A point mutation of PLB (L37A) that runs as monomer on SDS-PAGE showed no energy transfer, confirming its monomeric state in the membrane. We conclude that FET is a powerful approach for analyzing the oligomeric structure of PLB, and this method is applicable to other oligomeric proteins. Abstract | Full Text | PDF (230 kb)

• Interactions between Ca-ATPase and the Pentameric Form of Ph...

  Interactions between Ca-ATPase and the Pentameric Form of Phospholamban in Two-Dimensional Co-Crystals Biophysical Journal, Volume 90, Issue 11, 1 June 2006, Pages 4213-4223 David L. Stokes, Andrew J. Pomfret, William J. Rice, John Paul Glaves and Howard S. Young Abstract Phospholamban (PLB) physically interacts with Ca-ATPase and regulates contractility of the heart. We have studied this interaction using electron microscopy of large two-dimensional co-crystals of Ca-ATPase and the I40A mutant of PLB. Crystallization conditions were derived from those previously used for thin, helical crystals, but the addition of a 10-fold higher concentration of magnesium had a dramatic effect on the crystal morphology and packing. Two types of crystals were observed, and were characterized both by standard crystallographic methods and by electron tomography. The two crystal types had the same underlying lattice, which comprised antiparallel dimer ribbons of Ca-ATPase molecules previously seen in thin, helical crystals, but packed into a novel lattice with p222 symmetry. One crystal type was single-layered, whereas the other was a flattened tube and therefore double-layered. Additional features were observed between the dimer ribbons, which were substantially farther apart than in previous helical crystals. We attributed these additional densities to PLB, and built a three-dimensional model to show potential interactions with Ca-ATPase. These densities are most consistent with the pentameric form of PLB, despite the use of the presumed monomeric I40A mutant. Furthermore, our results indicate that this pentameric form of PLB is capable of a direct interaction with Ca-ATPase. Abstract | Full Text | PDF (875 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 3, 854-863, 1 February 2007

doi:10.1529/biophysj.106.095216

Biophysical Theory and Modeling

Membrane Assembly of Simple Helix Homo-Oligomers Studied via Molecular Dynamics Simulations

Lintao Bu, Wonpil Im and Charles L. Brooks^,  

Department of Molecular Biology (TPC6) and Center for Theoretical Biological Physics, The Scripps Research Institute, La Jolla, California

Address reprint requests to C. L. Brooks, Tel.: 858-784-8035.

Integral membrane proteins account for 30% of all proteins in the cell and play key roles in communication between the cell and its environment 1. Biological activity is clearly linked to protein folding, with misfolding leading to malfunction and disease for both membrane and non-membrane-associated proteins. However, in contrast to the wealth of available information regarding the folding of water-soluble proteins, relatively little is known about how membrane proteins fold to their native states. One idea, the two-stage model of integral membrane protein folding proposed by Popot and Engelman more than a decade ago 2,3, suggests a mechanism for helix-bundle membrane protein folding: the insertion of helix into the membrane (stage 1) and the assembly of the inserted helices in the membrane (stage 2). In earlier studies of insertion and folding, we explored aspects of stage 1 and observed a rather general mechanism governing these processes 4. Exploring the mechanism of assembly of membrane proteins is fundamentally related to our understanding of the biological functions of these proteins. Although the preponderance of transmembrane helical structure makes the prediction of membrane protein structures in one sense simpler than that of water-soluble proteins, the prediction of helix assembly in these systems remains an outstanding problem in computational biology because it requires a detailed structural and thermodynamic understanding of protein-protein and protein-lipid interactions 5.

The experimental determination of three-dimensional structures of membrane proteins is extremely difficult. Among the ∼30,000 protein structures found in the Protein Data Bank (PDB) 6, only 0.2% are of membrane proteins. Considering their biological importance and significant presence in genomes, a challenge to theory and computational biology is to assist experiment in understanding the structure and function of membrane proteins.

Several other methods have been used to explore the interfacial structures of transmembrane helices based on molecular dynamics (MD) simulation or energy minimization methods 7,8,9,10,11, using additional experimental information to identify the near-native structures. Engelman and co-workers developed a computational search algorithm to explore the interfacial structures of transmembrane helix homo-oligomers. They found that the van der Waals interactions alone provide sufficiently stabilizing forces to determine the specific helix association in phospholamban 7, glycophorin A 8, and synaptobrevin 9. Kukol et al. performed an exhaustive molecular dynamics global search protocol to obtain a structure of the M2 protein from the influenza A virus using the orientational data derived from site-directly infrared dichroism spectra as an unbiased refinement energy term 10. Torres et al. explored the interfacial structures of glycophorin A, the M2 protein, and phospholamban using global searching molecular dynamics simulations and helix tile as restraints 11. Ponder and co-workers performed an ab initio prediction of the glycophorin A structure using a novel potential smoothing and search algorithm 12. Helms and co-workers developed a novel scoring function for modeling structures of oligomers of transmembrane helices assuming that van der Waals interaction dominates in the packing of transmembrane helices 13. Kokubo and Okamoto used a replica-exchange Monte Carlo simulation method to study the structures of transmembrane helices of bacteriorhodopsin 14 and glycophorin A 15,16. Recently, Bowie and co-workers proposed a simple Monte Carlo method to study the association of helices using only sequence and native oligomerization state information 17,18. These approaches usually ignored the heterogeneous membrane/solvent environment and incorporated specific information from about the systems of interest from experiment, and consequently may not generalize to other cases. In this study we demonstrate that with only sequence and oligomerization state information we are able to assemble conformational ensembles that are in excellent agreement with experiment for three transmembrane assemblies, suggesting that the combination of a more physical model for the aqueous/membrane interface and enhanced sampling methods provide a more broadly applicable approach to predicting and modeling transmembrane assemblies.

An explicit membrane/solvent model provides the most detailed information to molecular modeling and represents the most accurate model 19,20,21,22. However, due to the increase in computing resources needed as the system size increases, significant efforts have been directed to the development of implicit membrane models. In general, continuum electrostatics can be used to define the electrostatic potential and the electrostatic solvation energy of a solute with arbitrary shape by solving the Poisson-Boltzmann equation using finite-difference methods 23,24,25. Unfortunately, the cost of solving the Poisson-Boltzmann equations has limited its application in molecular dynamics simulations 20,26. Alternatively, implicit membrane models based on generalized Born (GB) theories and dielectric screening functions have been used quite successfully to estimate the electrostatic solvation energy. Spassov et al. first extended the GB method to include an implicit membrane. They proposed an empirical approach to model the membrane within the context of a pairwise additive GB model 27. Lazaridis used an effective energy function approach to model protein solvation 28. Im et al. proposed an improved GB method based on a smooth dielectric boundary to study the structure, stability, and interactions of membrane proteins 29,30. For more information, see recent reviews by Brooks and co-workers 31,32. More recently, Feig and co-workers devised an implicit membrane model based on GB theories developed in the Brooks group 33.

Im and Brooks studied the interfacial folding and membrane insertion of designed peptides 4, using their implicit membrane GB model 29,30 and replica-exchange (REX) 34,35 molecular dynamics (MD). Their results demonstrated the mechanism of stage 1 of the two-stage model, and the success of using an implicit membrane model combined with advanced sampling methods to simulate biological membranes. In this article, we focus on the second stage, the assembly of transmembrane helices. Starting from an idealized helix, we sampled the conformational space of various oligomerization states using the implicit membrane GB model of Im et al. 29 REX/MD simulations, and the imposition of rotational symmetry to define the extent of oligomerization. We applied our method to predict the transmembrane structures of three peptides: glycophorin A (GpA), the M2 proton channel (M2-TMP), and phospholamban (PLB), which are experimentally known to form dimeric, tetrameric, and pentameric structures, respectively. We first explored the structures of each peptide in the native oligomeric state. We compared the predicted structures of GpA dimer, the M2-TMP tetramer, and PLB pentamer with experimental structures to examine our prediction with the native oligomerization state information as the structural constraint. Furthermore, we compared the potential energy between different oligomerization states for each peptide to address the challenging question of whether one can predict the native state energetically. In other words, whether one can predict the structures of helix homo-oligomers in membrane without using any experimental information.

We began our calculations with idealized α-helices, i.e., ϕ=−65° and ψ=−40°, for each peptide using the sequences given in Table 1. Each structure was oriented along the membrane normal in the membrane, and then rotated by 22.5° around the Z axis to produce 16 replicas. Each replica was then translated by 20Å from the symmetry axis in the X,Y plane, and these were taken as the initial structures of the monomers in our REX/MD simulation. We imposed m-fold rotational symmetry using the IMAGE facility in CHARMM 36 to provide putative oligomers of order m. For glycophorin A, which is a dimer in the native state, we imposed twofold, threefold, fourfold, and fivefold symmetries on the single peptide to simulate the structure of a dimeric, trimeric, tetrameric, and pentameric assembly, respectively. The M2 protein forms a tetramer and phospholamban forms a pentamer in their native states. We imposed twofold, threefold, fourfold, fivefold, and sixfold symmetries to simulate the structure of a dimeric, trimeric, tetrameric, pentameric, and hexameric peptide oligomer, for these two proteins. As a reference, a simulation of each peptide itself was also carried out.

Our studies were performed using the GBSW (a Generalized Born model with a simple SWitching function) module 29,30 in the CHARMM program 36. All MD simulations used a time-step of 2 fs and no cutoff for the nonbonded energy evaluation. The all-hydrogen parameter set PARAM22 37 of the CHARMM force field was used. The physical parameters representing the membrane in our GB model are 0.03kcal/(mol×Ų) for the surface tension coefficient (representing the nonpolar solvation energy), 25Å for the thickness of the membrane hydrophobic core, and 1Å for a membrane smoothing length over which the hydrophobic region is gradually changed to the solvent region. The planar membrane is perpendicular to the Z axis and centered at Z=0.

The MMTSB Tool Set 38 was used to control the REX simulations. We used 16 replicas that were distributed over an exponentially spaced temperature range from 300K to 600K. Langevin dynamics with a friction coefficient of 5.0ps⁻¹ for heavy atoms was used. A cylindrical harmonic restraint with a 25Å radius and a force constant of 1.0kcal/(mol×Ų) was applied to prevent the peptides from drifting radially away from each other, i.e., away from the symmetry axis. (Note that this is much larger than the radius of any of the assemblies we studied.) The REX/MD simulations were carried out for 10ns for each oligomerization state of each peptide. Every 1ps, a replica exchange was attempted and the coordinates were saved for further analyses. The pairwise exchange ratio was ∼40% for each run.

Using the CLUSTER facility in MMTSB Tool Set 38, we clustered the sampled structures in the native state of each peptide. Due to the size limitation of the ensemble of structures used in the CLUSTER facility, we collected every other structure during the last 7ns of the REX/MD simulation providing 3500 structures to be used in clustering stage. We chose the structure located at the center of the largest cluster as the predicted structure.

Glycophorin A (GpA) is one of the most abundant proteins located on the surface of red blood cells; however, despite its ubiquitous presence, its function remains unknown. It is also one of the most well-studied model systems in the field of helix-helix interactions in membranes 39,40. The NMR structure of glycophorin A in micelles (PDB:1AFO) 41 shows that it forms a right-handed helical dimer with the packing motif LIXXGVXXGV. The two Gly residues form a flat surface to facilitate tight packing of the backbone atoms 42. The two Val residues play the key role in the van der Waals interaction between the transmembrane helices.

The 3500 sampled structures formed two clusters with group size of 1865 (53%) and 1635 (47%) structures, respectively. The representative structure from the largest cluster has a C_α root mean-square derivation (RMSD) value of 2.2Å relative to the experimental structure, whereas the representative structure from the other cluster has a C_α RMSD value of 7.6Å.

Figure 1AB, shows the interhelical crossing angle of the simulated dimeric structures and C_α RMSD of the dimeric structures relative to the native structure (PDB:1AFO) as a function of time. For clarity, only five trajectories out of 16 during the last 7ns are shown. We can see a few transitions between the two configurations (left-handed dimer and right-handed dimer) occurring at high temperatures, indicating the sampling efficiency of REX/MD simulation. The RMSD is well-correlated with the interhelical crossing angle. Figure 1CD, shows the distribution of crossing angles and C_α RMSD of the structures sampled at the lowest temperature (300K) during the MD simulation. Based on the distribution of crossing angles, the helices could be clustered into two distinct families of conformations: a right-handed dimer (crossing angle at −50°), and a left-handed dimer (crossing angle at 50°). The right-handed dimer has a most probable RMSD value of 2.2Å, whereas the most probable RMSD value of the left-handed dimer is 7.8Å. The solid line in Figure 1CD, shows the integrated population. While we see both conformations occurring with some probability, the native right-handed dimer occupies >60% of the total conformations sampled. These results are relatively consistent with the clustering results using the CLUSTER facility.

Display large version of this figure

Figure 1

(A) Interhelical crossing angle of GpA dimeric structures during REX/MD simulation as a function of time. The crossing angle is measured by the angle between two principal axes defined by the backbone heavy atoms of each monomer. The negative sign means that it forms a right-handed dimer and the positive sign stands for a left-handed dimer. For clarity, only five trajectories are shown. (B) RMSD of the C_α atoms of GpA sampled structures relative to the native structure (1AFO) as a function of time. (C) The distribution of crossing angle of GpA dimer. The solid curve represents the integrated population. (D) The distribution of RMSD of C_α atoms of sampled structures relative to the native structure.

Figure 2A shows the structure of a representative dimer model of glycophorin A derived from our simulation. As shown in Figure 2BC, the comparison of contact map for the C_α-C_α distances between our model and the NMR structure reveals that the interfacial residues of our model, including Leu⁷⁵, Ile⁷⁶, Gly⁷⁹, Val⁸⁰, Gly⁸³, and Val⁸⁴, are identical with those of the solution NMR structure.

Display large version of this figure

Figure 2

(A) The representative dimer model of GpA derived from our simulation. The stick representations show the interfacial residues. Contact map of the C_α atoms in GpA NMR structure model (B) and our predicted model (C). Color is coded by the distance between two C_α atoms.

The M2 protein from Influenza A contains 97 residues and is a proton selective ion channel that forms a left-handed tetrameric helical domain 43. The structure of a 25-residue (from Ser²² to Leu⁴⁶) peptide, which is also called M2-TMP, was recently determined in a DMPC bilayer using rotational echo double-resonance solid-state NMR 44.

The 3500 sampled structures from the REX/MD simulation formed five clusters with group size of 1373 (39%), 919 (26%), 849 (24%), 226 (6%), and 133 (4%) structures. The representative structure from the largest cluster has a C_α RMSD value of 2.7Å relative to the experimental structure (PDB:1NYJ). The representative structures from other clusters have a C_α RMSD value of 8.5, 5.1, 3.9, and 4.1Å relative to the experimental structure.

Figure 3AB, show the interhelical crossing angle of the sampled tetrameric structures and C_α RMSD of the sampled tetrameric structures relative to the native structure as a function of simulation time. Again, the RMSD is well-correlated with the interhelical crossing angle. The distribution of RMSD in the sampled structures at the lowest temperature in Figure 3D, showing the existence of five clusters, is consistent with the clustering results using the CLUSTER facility. Based on the distribution of crossing angles in Figure 3C, the helices could be clustered into three families of conformations: two right-handed tetramers (crossing angle at −25° and −5°), and a left-handed tetramer (crossing angle at 35°). The left-handed tetrameric state, which is also the native state, has the largest population of 50%. The population of the two right-handed tetramers is ∼30% and 20%, respectively. The left-handed tetramer has a most probable RMSD value of 2.6Å at the lowest temperature, whereas the most probable RMSD value of the right-handed tetramer is 8.9Å.

Display large version of this figure

Figure 3

(A) Interhelical crossing angle of M2-TMP tetramer structures derived from simulation as a function of time. For clarity, only five trajectories are shown. (B) RMSD of the C_α atoms of M2-TMP sampled structures relative to the native structure (1NYJ) as a function of time. (C) The distribution of crossing angle of M2-TMP tetrameric structures. (D) The distribution of RMSD of C_α atoms of sampled structures relative to the native structure.

Figure 4A shows the structure of a representative tetrameric model of M2-TMP derived from our simulation. Figure 4BC, shows that the interfacial residues, including Val²⁷, Ser³¹, Gly³⁴, His³⁷, Leu³⁸, and Trp⁴¹, are identical with those of the experimental structure derived from solid-state NMR.

Display large version of this figure

Figure 4

(A) The representative tetramer model of M2-TMP derived from our simulation. The stick representations show the interfacial residues. Contact map of the C_α-C_α in the M2-TMP NMR structure model (B) and our predicted model (C).

Located in the membrane of the cardiac sarcoplasmic reticulum, phospholamban (PLB) is involved in regulation of a Ca²⁺ pump 45. We compare our predicted model with a model structure (PDB:1PLN) 46, which was created by the direct structure modeling of mutagenesis data. Structures have been determined for the helical monomer in the solid state using rotational echo double-resonance 47, and more recently for the pentameric structure in micelles via solution NMR methods 48. The model structure differs very little from the more recent experimental structure.

The 3500 sampled structures formed four clusters with group size of 1653 (47%), 1233 (35%), 498 (14%), and 116 (4%) structures. The representative structure from the largest cluster has a C_α RMSD value of 0.62Å relative to the model structure (PDB:1PLN). Our predicted model shows similar agreement with the pentameric NMR structure (PDB:1ZLL)48 which, for the TM region, is in the range of 0.71–0.94Å compared with the 20 NMR structures and has an average C_α RMSD of 0.84Å. The representative structures from other clusters have a C_α RMSD value of 3.4, 1.6, and 4.6Å, respectively, relative to the model structure (PDB:1PLN).

Figure 5A shows the distribution of crossing angle of phospholamban pentameric structures at 300K, which suggests that the helices only form a left-handed pentamer (crossing angle at 19°). Figure 5B shows the distribution of C_α RMSD of phospholamban sampled structures at 300K relative to PDB:1PLN during the MD simulation. The distribution of RMSD, which is characterized by the existence of four clusters with the most probable RMSD values of 0.62, 1.7, 3.3, and 4.6Å, is consistent with the clustering results using the CLUSTER facility.

Display large version of this figure

Figure 5

(A) The distribution of crossing angle of the sampled PLB pentameric structures. (B) The distribution of RMSD of C_α atoms of the sampled structures relative to the model structure in model structure (PDB:1PLN).

Figure 6A shows the structure of the representative pentameric model of PLB derived from our simulation. Figure 6BC, illustrates that the positions of the interfacial residues, including Leu³⁷, Ile⁴⁰, Leu⁴⁴, Ile⁴⁷, and Leu⁵¹, are identical to PDB:1PLN.

Display large version of this figure

Figure 6

(A) The representative pentamer model of PLB derived from our simulation. The stick representations show the interfacial residues. Contact map of the C_α-C_α in PLB model derived from experimental data (B) and our predicted model (C).

In this section, we compare the potential energy between the different oligomerization states of each peptide. The question we would like to address is whether we can predict the native oligomerization state of each peptide energetically. We assume that the native state not only has the lowest free energy, but also the lowest potential energy.

Figure 7A shows the potential energy profile of each peptide for various oligomerization states at the lowest temperature (300K), averaged over the last 7ns of the REX/MD simulations, relative to the corresponding monomeric states. The potential energy of each oligomerization state of each peptide converged after the initial 2ns of simulation (data not shown). In the case of GpA, the dimeric state does not have the lowest potential energy. Table 2 shows the decomposition of the potential energy. The differences in potential energy between monomeric/dimeric state and dimeric/trimeric state are −21.1kcal/mol and −15.3kcal/mol, respectively. The differences are dominated by the differences in van der Waals interaction between monomeric/dimeric state and dimeric/trimeric state, which are −21.8kcal/mol and −13.3kcal/mol, respectively. Clearly, van der Waals interactions between the interfacial residues play the key role in the packing of helices in our model.

Display large version of this figure

Figure 7

Potential energy (A) and free energy (B) profiles of GpA, M2-TMP, and PLB oligomers at the lowest temperature (300K) during the REX/MD simulations, relative to the corresponding monomers.

In the case of M2-TMP, as shown in Figure 7A, the tetrameric state, which is the native state, does not have the lowest potential energy. Table 3 shows the decomposition of the energy. The potential energy differences between monomeric/dimeric state, dimeric/trimeric state, and tetrameric/pentameric state are all dominated by the differences in interhelical van der Waals interactions. Interestingly, the difference of interhelical van der Waals interaction in two adjacent oligomerization states becomes smaller as the oligomerization number increases. This trend is also seen in glycophorin A and phospholamban. The potential energy difference between trimeric state and tetrameric state is relatively small, compared to the other two adjacent states. The differences in interhelical van der Waals interaction between trimeric/tetrameric state are canceled by other unfavorable energy differences, such as internal energy, which is 418.6±0.5kcal/mol for monomer, dimer, and trimer, whereas it is 421.7±0.2kcal/mol for tetramer, pentamer, and hexamer. The van der Waals interaction between the helices in the pentameric state and hexameric state is identical, whereas the electrostatic interaction dominates the potential energy difference between pentameric state and hexameric state. This is perhaps due to the electrostatic interaction between the polar residues located in the interhelical interface, mainly Ser³¹, His³⁷, and Trp⁴¹. As demonstrated in the solid-state NMR structure, the distance between His³⁷ and Trp⁴¹ is <3.9Å, which suggests that the interaction between His³⁷ and Trp⁴¹ plays the key role in sterically closing the channel 44. We observed that the close packing of pentamer and hexamer does not have the correct handedness. As shown in Fig. 8, the sampled structures at pentameric and hexameric states are mostly right-handed, whereas the native structure should be left-handed.

Display large version of this figure

Figure 8

The distribution of crossing angle of M2-TMP sampled structures at tetrameric, pentameric, and hexameric states.

In the case of PLB, as shown in Figure 7A, the pentameric state, which is the native state, has the lowest potential energy. As shown in Table 4, the difference between interhelical van der Waals interactions again dominates the potential energy difference between adjacent oligomerization states. Since the interfacial residues in PLB are all hydrophobic residues, the electrostatic interaction between the helices does not contribute significantly to the stabilization of the oligomers.

We have investigated the membrane assembly of GpA, M2-TMP, and PLB peptide, using REX/MD and an implicit membrane GB model. Our approach is quite successful in predicting the structures of homo-oligomers, using only the native oligomerization state as a structural constraint. It is noteworthy that this property can often be gleaned from measurements utilizing analytical ultracentrifugation, equilibrium dialysis, and other biochemical approaches without the necessity of atomic level structural information 41,49,50,51. For our predicted models, we find the RMSD value of C_α atoms relative to the corresponding experimental and model structures are 2.2Å (GpA), 2.7Å (M2-TMP), and 0.62Å (PLB), respectively. Also of interest is the observation that a distribution of conformations appear to be present in each case. Whether this is a true reflection of some level of conformational heterogeneity or a limitation of our model remains to be investigated.

Using only the peptide sequence we do not always predict the native oligomerization state as predominant based on energetic criteria. We successfully predicted the native oligomerization state for PLB, but not for GpA and M2-TMP. One explanation for this may be that we did not consider the entropy loss during helix association. Shown in Table 5 is an estimation of translational, rotational, and conformational entropy. The translational entropy and rotational entropy were calculated based on principal RMS fluctuations of the center of mass or Euler angles 52. The translational entropy can be expressed as

(1)

(2)

More accurate evaluation of entropy loss may be needed to improve the first-principles calculation of folding and oligomerization equilibria for these peptides. However, we also find that the interhelical van der Waals interaction dominates in the packing of the GpA, M2-TMP, and PLB peptides. Thus, we anticipate that the interhelical van der Waals interaction is overestimated in our GBSW model, since we did not include peptide-lipid dispersion interactions that should compete with the peptide-peptide interactions. Currently, we are extending our implicit membrane model to include interactions between protein and lipid. The optimization of parameters in our GBSW model is ongoing.