Structure and Dynamics of Phospholamban in Solution and in

Jan 20, 2005 - ACS eBooks; C&EN Global Enterprise .... Differences between the solution and micelle structures of PLB are mainly in the orientation of...
0 downloads 0 Views 468KB Size
1780

Biochemistry 2005, 44, 1780-1792

Structure and Dynamics of Phospholamban in Solution and in Membrane Bilayer: Computer Simulations† Yao Houndonougbo,‡ Krzysztof Kuczera,*,# and Gouri S. Jas*,‡ Departments of Chemistry and Molecular Biosciences, UniVersity of Kansas, 1251 Wescoe Hall DriVe, Lawrence, Kansas 66045, and Higuchi Biosciences Center, UniVersity of Kansas, 2095 Constant AVenue, Lawrence, Kansas 66047 ReceiVed June 5, 2004; ReVised Manuscript ReceiVed NoVember 5, 2004

ABSTRACT: We have performed molecular dynamics simulations of phospholamban (PLB), a 52-residue integral membrane protein that inhibits calcium ATPase in the cardiac sarcoplasmic reticulum. We present a microscopic description of the structure and dynamics of PLB in solution and membrane environments, based on 10 ns molecular dynamics simulations of PLB in lipid bilayer and 5 ns simulations in methanol and water, and a water-soluble model of PLB in water. Throughout the simulations, PLB retains its “L” shape, with two well-defined helical domains at the N- and C-termini. In the simulations of PLB in methanol and water, the helices were almost perpendicular, with average interhelix angles of 54 ( 13° and 63 ( 15°, respectively. In the lipid bilayer trajectory, both the interhelix angle and its fluctuations were larger, with an average of 130 ( 19° and with the transmembrane C-terminal approximately perpendicular to the bilayer plane. The internal dynamics of phospholamban is characterized by large amplitude collective motions of the two helical domains: hinge bending, twisting of both N- and C-terminal helices, and flexing of the C-terminal helix. The central linker of PLB is highly flexible, due mostly to elastic deformations of this region. The simulation results are in good agreement with NMR data on PLB secondary structure and helix orientations in solution, micelles, and lipid bilayers, as well as fluorescence measurements of interdomain distances. Our most interesting findings involve the details of the PLB dynamics, which are difficult to obtain by experimental approaches. Two kinds of motions of the helical domains found in the simulations can clearly have functional roles. The population of conformations with relatively open interdomain angles, as well as large fluctuations of this coordinate in the bilayer, allows the N-terminal helix to come into contact with the PLB binding site on the calcium ATPase, while the presence of twisting motions around its axis enables the helix to orient the correct face to the binding site.

Phospholamban (PLB) is a 52-residue integral membrane protein that inhibits calcium ATPase in the cardiac sarcoplasmic reticulum (SR) (1, 2). Inhibition of the calcium ATPase is removed upon PLB phosphorylation at serine 16 and/or threonine 17. Recent NMR solution studies in a mixture of chloroform/methanol of the C41F mutant of pig PLB showed that the protein takes on an “L-shape”, with mainly R-helical structure in the cytoplasmic domain (residues 1-16) and transmembrane domain (residues 21-52), and a type III β-turn conformation in the cytoplasmic connecting region (residues 17-20) (3). These results confirm previous NMR studies on the overall shape and secondary structure of PLB in mixed solvents (4, 5). An analogous structure, with an overall “L-shape” and predominantly R-helical in the two domains was also found in NMR experiments on monomeric phospholamban in dodecylphos† This work was supported by American Heart Association Grant 0225485Z providing a postdoctoral fellowship for Y.H. and Scientist Development Grant 0360061Z for G.S.J. * To whom correspondence should be addressed. Phone: (785) 8645060. Fax: (785) 864-5396. E-mail: [email protected] or gourijas@ ku.edu. ‡ Higuchi Biosciences Center. # Department of Chemistry and Department of Molecular Biosciences.

phocholine (DPC) micelles (6). Differences between the solution and micelle structures of PLB are mainly in the orientation of the cytoplasmic helix. The two phosphorylation sites (Ser16 and Thr17) and the two arginine residues (Arg9 and Arg13) of the cytoplasmic domain face the transmembrane domain in the structure in chloroform/methanol solution (3), whereas these four residues are oriented toward the bulk solvent in the DPC micelle structure (6). Experimental studies also provide information about the orientations of the PLB helical domains. The measured interhelical angle is 68 ( 23° in organic solvent (3), 80 ( 22° in DPC micelles (6), and in the 80-120° range in oriented lipid bilayers (7, 8). Solid-state NMR (7) and infrared (8) measurements indicate that the transmembrane domain is perpendicular while the cytosolic domain is approximatively parallel to the membrane plane. It has been suggested that, to interact with calcium ATPase, PLB must adopt an extended structure (3, 9). In some studies, PLB has been represented as a continuous R-helix (interhelix angle 180° in our convention) (4, 10). PLB can be found in different oligomeric states. Sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDSPAGE) experiments suggested that it has a strong tendency for forming stable pentamers (11-13). The specific residues

10.1021/bi0488404 CCC: $30.25 © 2005 American Chemical Society Published on Web 01/20/2005

Computer Simulations of Phospholamban that are important in maintaining the pentameric structure are confined to the hydrophobic transmembrane domain (1416). The pentamer is noncovalent and has been modeled as a left-handed symmetrical coiled-coil, with two different relative orientations of the five helices (10, 12, 17). The model by Simmerman et al. (18), in which the pentamer is stabilized by a leucine zipper-like arrangement, is consistent with the labeling experiments in SDS (19) and recent studies in a lipid bilayer (17). The residues of the membrane-spanning helix are hydrophobic, so PLB is difficult to solubilize. To make PLB more amenable to experimental study, a water-soluble model of phospholamban (WSPLB) has recently been produced by mutating several hydrophobic residues (20). In that study, it was shown that WSPLB has the same properties as regular PLB, that is, helical structure, oligomer formation, and oligomer stabilization upon phosphorylation. Other watersoluble PLB models have been produced previously (21, 22). Thus far, phospholamban has mostly been the object of experimental studies, as the lack of atomic-resolution structural information limited the application of molecular modeling approaches (3, 6). These methods, and especially molecular dynamics simulations, have become a powerful tool in the study of biological molecules, providing detailed, realistic information about their structure, motion, and interactions in different environments (23-26). The power of computer simulations lies in their atomic detail, high temporal resolution, and the complementarity of results to experimental data. Thus, the optimal situation is when both experimental and simulation data are available, enabling verification of calculations by comparison with observations and insight into the microscopic mechanisms of observed effects through models. In this work, we present results of molecular dynamics simulations that model the behavior of PLB in methanol, aqueous solution, and a membrane bilayer, as well as the water-soluble version of the protein (WSPLB) in aqueous solution, using explicit all-atom representations of both solute and solvent. We analyze the hydrogen bonding patterns, secondary structure, domain motions, helix orientations, solvent accessible surface areas, and translational and rotational diffusion coefficients. We find that PLB tends to retain the mostly R-helical structure of its transmembrane and cytosolic domains in all studied environments. Also, the protein tends to undergo motions involving significant reorientations of the two helical domains. Interestingly, we find larger values of both the average and fluctuations of the interhelical angle of PLB in the membrane environment, compared to solution. Thus, the membrane environment tends to induce a PLB structure that is more elongated and more flexible, which is compatible with the requirements for PLB interactions with the calcium ATPase. METHODS Simulations of PLB in Methanol and Water. The initial coordinates of monomeric phospholamban, in which Cys (41) is replaced by Phe (C41F mutant), were taken from the Protein Data Bank file 1FJK, generated from a highresolution NMR study in a mixture of chloroform and methanol (3). The CHARMM program was used to add hydrogen atoms needed in the all-atom model, and the

Biochemistry, Vol. 44, No. 6, 2005 1781 positions of the constructed atoms were subjected to a brief energy minimization in a vacuum, with the rest of the protein kept fixed (27). PLB has a net charge of +3e in normal charge states of ionizable groups corresponding to pH 7. To neutralize the system at biologically realistic conditions, six sodium and nine chlorine counterions were added using the program Solvate 1.0 by Helmut Grubmueller (28). The protein and counterions were placed in truncated orthorhombic cell with dimensions 64.0 × 64.0 × 84.0 Å. For PLB in methanol, the cell was overlaid with 2239 methanol molecules, generating a simulation system of 14 346 atoms. For PLB in water the cell was overlaid with 5443 TIP3P water molecules (29), with the final simulation system of 17 241 atoms. For each cell, a brief energy minimization followed by a 50 ps MD equilibration were applied with the protein fixed, followed by 50 ps MD with the whole protein-solvent system being allowed to move. Finally, 5 ns MD trajectories of the unconstrained systems were generated. In both equilibration and trajectories generation, the leapfrog algorithm was used for integration of equations of motion, the Hoover method was employed to maintain constant temperature at 300 K (30), and the Langevin piston method was used for maintaining constant pressure at 1 atm (31). Periodic boundary conditions were applied. SHAKE constraints were used to fix all the bonds involving hydrogen atoms (32) to integrate the equations of motion with a 2 fs time step (33). For van der Waals terms, standard nonbond cutoffs were employed, with an atombased 12.0 Å nonbonded cutoff distance and a switching function between 10.0 and 12.0 Å. To avoid truncation of the long-range electrostatic interactions, the particle-mesh Ewald (PME) method was used, with a charge mesh of 64 × 64 × 90 grid points, grid spacing of approximately 1 Å, the Ewald parameter κ set to 0.34 Å-1, and direct interactions truncated at 12.0 Å. Simulation of WSPLB in Water. The WSPLB coordinates were built from the file 1FJK (3). The coordinates for the side chains of the mutated residues (20) were constructed in extended conformations using the program CHARMM (27). WSPLB was placed in a 63 × 63 × 82 Å truncated orthorhombic cell, and overlaid with 5636 TIP3P waters, creating a system of 17 805 atoms. The trajectory setup and generation procedure for the 5 ns trajectory of WSPLB in water were the same as described above for PLB. Simulation of PLB in Membrane Bilayer. For the simulation of PLB in the membrane bilayer, the initial coordinates of monomeric phospholamban were taken from the first of the 20 structures in Protein Data Bank file 1N7L (6). To construct the initial configuration of the protein-membrane complex, we used the CHARMM membrane tools based on the method developed by Woolf and Roux (34, 35). The cross-section area of one DPPC phospholipid is 62.9 Å2 (36) and approximatively 180.6 Å2 for PLB. The bilayer was constructed with 56 DPPC molecules for the top layer and 59 DPPC molecules for the bottom layer, producing a system with dimensions 61.0 × 61.0 Å in the XY (membrane) plane. To create the protein-lipid system starting configuration, we randomly selected preequilibrated and prehydrated DPPC molecules from a library of 2000 structures sampled from a mean-field simulation (37). The center of the membrane was placed at Z ) 0, and the transmembrane domain of the protein was oriented along the Z axis perpendicularly to the

1782 Biochemistry, Vol. 44, No. 6, 2005 plane of the bilayer, with the cytoplasmic helix laid on the top layer surface. Systematic rigid body rotations and translations of the lipids with their primary TIP3P water molecules were used to remove unfavorable core-core overlaps present in the random initial configuration. The system was then energy-minimized with a stepwise increase of the lipid van der Waals radii to their full values, while the protein was kept fixed. The DPPC/PLB complex has a net charge of +4e in normal charge states of ionizable groups corresponding to pH 7. Thus, we added two sodium and six chlorine ions to neutralize the system at close to physiological ionic strength of 0.15 M, using a modified version of the program SOLVATE 1.0 (28). The membrane-protein complex was overlaid with a preequilabrated TIP3P water box of size 61.0 × 61.0 × 85.0 Å, and those water molecules that overlapped the protein, membrane, previously placed primary waters, or counterions were removed. The resulting system comprises of 31 058 atoms, and includes 5198 water molecules, 115 DPPC molecules, two Na+, and six Cl- ions. To relax the system, a brief energy minimization was followed by a 130 ps equilibration with the temperature kept at 330 K using Langevin dynamics (31). During this stage, the protein backbone, the centers of mass of the phospholipid headgroups and waters were restrained by harmonic forces which were gradually released over the first 100 ps. Further equilibration involved 25 ps of Langevin dynamics at 330 K and 25 ps of velocity scaling to adjust the temperature to 300 K. In this phase, the drift of the bilayer and PLB was prevented by constraining the center of mass of the protein to remain close to the Z axis and the center of mass of the lipids to remain close to the XY plane. Finally, 10 ns MD trajectories of the system were generated with the same cylindrical and planar constraints described above. The Hoover method was employed to maintain a constant temperature at 300 K (30) and the Langevin piston method (31) was used to maintain constant pressure at 1 atm. Leapfrog integration was used (38), with periodic boundary conditions. All bonds involving hydrogen atoms were fixed using SHAKE constraints, enabling a 2 fs time step (33). An atom-based 12.0 Å nonbonded cutoff distance with a switching function between 10.0 and 12.0 Å was used for van der Waals terms. The PME method was used to avoid truncation of the electrostatic interactions. The direct interactions were truncated at 12.0 Å. A charge mesh of 72 × 72 × 96 grid points was used, giving a grid spacing of approximately 1 Å, and the Ewald parameter κ was set to 0.34 Å-1. An instantaneous configuration of the membrane system is shown in Figure 1. Programs and Force Fields. The program CHARMM version 30 was used to perform all the simulations. The CHARMM version 22 all-atom topology was employed (27, 39). For PLB, the CMAP correction to the peptide backbone force field was used, designed to accurately reproduce the quantum mechanical φ/ψ map of dipeptides (40). Generation of a 5 ns trajectory took about 1000 CPU-hours for the solution and 1300 CPU-hours for the DPPC bilayer systems, running in parallel on 2.4 GHz dual-processor Pentium IV workstations. Trajectory Analysis. In each of the four simulations, frames were saved every 0.5 ps. To follow internal motions, oriented trajectories were generated by overlaying the backbone atoms

Houndonougbo et al.

FIGURE 1: Molecular graphics view of phospholamban (red) in DPPC bilayer membrane in aqueous solution of 0.15 M [Na+,Cl-]. The membrane bilayer is represented by carbon chains (gray) and phosphate atom (green). The chlorine ions are in orange, the sodium ions are in yellow, and the water molecules are in light blue.

of the trajectory structures on starting coordinates to remove the effects of the overall translation and rotation. These oriented trajectories were used to calculate root-mean-square deviations (RMSDs) from the starting structures and rootmean-square fluctuations (RMSFs) within the trajectories and to visualize the dynamics. We analyzed the helicities and hydrogen bonding patterns in the simulations. The helical hydrogen bonding information was extracted from the trajectories using the following criteria. A hydrogen bond (A‚H-D) was defined by an A‚‚‚H distance less than 3.5 Å and an A-H-D angle greater than 120°. The helicities of the two domains were calculated as the ratio of the number of R-helical residues in a domain to the number of R-helical residues in the starting PDB file for that domain. We follow the definition of Kabsch and Sander (41) and define residues i + 1, i + 2, i + 3 as helical if the hydrogen bond (i, i + 4) exists. To examine the interactions of PLB with the lipid membrane, we measured the number of contacts between the protein and the lipid headgroups as a function of time. A contact occurs when a lipid-protein interatomic distance is less than 3.5 Å. We have performed ED analysis of the trajectories to analyze the large amplitude collective motions of PLB (42, 43). The mass-weighted covariance matrix Cij of position fluctuations of the backbone CR, C, and N atoms, was calculated according to

Cij ) xmimj 〈(xi - 〈xi〉)(xj - 〈xj〉)〉

(1)

where xi and mi are individual atomic Cartesian coordinates and masses, and 〈‚‚‚〉 is an average over the trajectory in which the effects of the overall translations and rotations of the protein have been removed. The eigenvectors with the largest eigenvalues give models of collective coordinates that best describe the calculated structural fluctuations (42, 43),

Computer Simulations of Phospholamban

Biochemistry, Vol. 44, No. 6, 2005 1783

while the eigenvalues may be used to calculate the quasiharmonic conformational entropy of the system (44). The translational diffusion coefficient Dt of the protein in solution was calculated from the mean square displacement 〈∆r2(t)〉 of the center of mass of the protein as a function of time using the expression (45):

Dt )

〈∆r2(t)〉 6t

(2)

where ∆r(t) is the displacement of the center of mass in time t and 〈‚‚‚〉 denotes an average over starting structures. To calculate the rotational diffusion coefficients, we assumed that the principal axes of the diffusion tensor and the moment of inertia tensor coincide. For a molecule undergoing rotational diffusion, the autocorrelation function for axis i is (46):

bi(0)‚n bi(t)〉 ) e- (3Dr-Di)t i ) x, y, z C1, i(t) ) 〈n

(3)

where b ni is the unit axis vector, Dr ) (Dr,x + Dr,y + Dr,z)/3 is the average rotational diffusion coefficient, and Dr,x, Dr,y, and Dr,z are the rotational diffusion coefficients around the individual axes. The individual diffusion coefficients Dr,x, Dr,z, and Dr,z were calculated by first fitting to single exponential decays the autocorrelation functions for the three individual axes x, y, z and then by using eq 3. The five rotational correlation times τi, i ) 1, ..., 5 of the secondorder correlation function of C2(t) ) 〈3(n bi(0)‚n bi(t))2 - 1〉/2 of an arbitrary molecule-fixed vector b n were calculated from the values of the three diffusion coefficients found in the trajectories (46, 47). For PLB domains, we have calculated the diffusion coefficients D⊥ for reorientation of the helical axis (tumbling) and D| for reorientations along the helix axis (spinning). To calculate D⊥, we evaluate the autocorrelation function C1,⊥(t) ) e-2D⊥t of the helix axis vector. Next, we compute D| from the correlation functions of the two small axes of the domain moment of inertia and eq 3. The helix axes were determined from the CR coordinates using the method of Aqvist (48). The statistical errors of the translational and rotational diffusion coefficients were estimated from the difference of the calculated values obtained by splitting the trajectory into two equal parts and by taking into account the t-coefficients for the 0.95 confidence level (49). RESULTS OVerall Structure and Dynamics. The final and starting structures of all simulations are compared in Figure 2. For the simulations of PLB in water and methanol, the time evolution of the RMSDs from the NMR starting structure are shown in Figures 3-5. In methanol and water, the average overall RMSD is 3.6 Å. The average RMSDs within the N-terminal (cytoplasmic domain, residues 1-16) and the C-terminal (transmembrane domain, residues 2152) are in the 1.0-2.0 Å range. The RMSDs of individual residues (Figure 6) mostly follow the expected pattern of lower deviations within helices and higher deviations in the connector and termini. The average atomic RMSFs of PLB N- and C-terminal domains (Figure 7) were 0.4 and 0.9 Å in methanol, and 0.4 and 1.1 Å in water, respectively. Thus,

FIGURE 2: Superposition of the CR atoms of the starting structure (black) and the final structure (gray) the simulations.

FIGURE 3: Instantaneous backbone RMSDs (for residues 4-49) from the starting structure 1FJK for PLB in methanol and water and WSPLB in water.

in methanol and water, the two domains of PLB tend to remain stable, close to their experimental structures. The relatively large RMSDs from the starting NMR structure are due to rigid-body motions of the helices. The average RMSDs from starting structure and RMSFs around the trajectory average tend to be lower in the cytoplasmic domain than in the transmembrane domain. The simulation results for WSPLB in aqueous solution are similar to our calculations for the regular PLB in methanol and water. Although the overall WSPLB RMSD from the 1FJK starting structure reaches up to 8 Å (Figure 3), the RMSDs within the two domains are smaller, in the 1-4 Å range. The average atomic RMSFs in WSPLB were 0.5 Å in the N-terminal domain and 1.3 Å in the C-terminal domain. The RMSDs from starting structure and trajectory

1784 Biochemistry, Vol. 44, No. 6, 2005

Houndonougbo et al.

FIGURE 4: Instantaneous backbone RMSDs from the starting structure 1FJK for PLB in methanol and water and WSPLB in water. Left: N-terminal domain, residues 4-16. Right: C-terminal domain, residues 21-49.

FIGURE 5: Instantaneous backbone RMSDs from the starting structure 1N7L for PLB in DPPC bilayer. Top: whole protein, residues 2-50. Middle: N-terminal domain, residues 2-16. Bottom: C-terminal domain, residues 22-50.

average RMSFs from individual residues (Figures 6 and 7) are similar as for regular PLB. The simulation of PLB in the DPPC membrane bilayer also exhibits features that are similar to the results found for PLB in solution. While the overall RMSD from the 1N7L structure reaches up to 9 Å (Figure 5), the RMSDs within the two domains remain small, 0.5-2.6 Å. The average atomic RMSFs in DPPC were 0.7 Å in the N-terminal

domain and 0.8 Å in the C-terminal domain. Except for the region centered at residue 30, the RMSDs from the starting structure and trajectory average RMSFs mostly follow the expected pattern, with the largest values at the helix ends in the central linker region (Figures 6 and 7). The overall structure and dynamics seen in the simulations thus agree qualitatively with experimental data. The simulations predict stable structures of the N- and C-terminal domains for the regular PLB in water, methanol, and DPPC bilayer, and for WSPLB in water. The fact that deviations from starting structures are much larger for the whole protein than for its domains indicates the existence of large amplitude domain motions for PLB in all studied environments. Interestingly, the structural fluctuations appear to be systematically higher in the C-terminal domain compared to the N-terminal in all our PLB simulations, suggesting that the former domain is more conformationally labile. Hydrogen Bonds. To analyze the secondary structure of PLB in the different environments, we examined the helicity of the two domains using the hydrogen bonding patterns found in our trajectories (see Methods). In the NMR starting structures, the two helices involve residues 4-16 and 2149 in chloroform/methanol solution (3) and residues 2-16 and 22-50 in DPC micelles (6). The experimentally observed helicities of the N-terminal domain are thus 81 and 93%, in solution and micelles, respectively. For the Cterminal domain, the observed helicity is 90% in both environments. The domain helicities found in the simulations are displayed in Figure 8 and Table 1. Two trends may be seen in the simulation results. The first is that, while the average trajectory helicities tended to be somewhat lower than corresponding experimental data, the phospholamban N- and C-terminal domains tend to retain their R-helical structure throughout the simulations. The

Computer Simulations of Phospholamban

Biochemistry, Vol. 44, No. 6, 2005 1785

FIGURE 6: Distribution along the polypeptide chain of backbone atoms RMSDs from the simulation starting structures of PLB. Left: deviations of PLB simulations in water and methanol, and WSPLB in water from structure 1FJK. Right: deviations of PLB simulations in DPPC bilayer from structure 1N7L.

FIGURE 7: Distribution along the polypeptide chain of backbone atom RMSFs around trajectory average structures. Left: simulations of PLB in water and methanol, and WSPLB in water. Right: simulations of PLB in DPPC bilayer.

average helicities were above 80% for both domains in most of the simulations. The exceptions were the PLB and WSPLB N-terminal helices in water, which had average helicities of 75 and 77%, respectively. The second trend was for the occurrence of rare structural fluctuations involving instantaneous helix breaking, seen in all the trajectories. In the simulation in membrane bilayer, several hydrogen bonds involving PLB side chains were found. These hydrogen bonds involved residues Gln-5‚‚‚Arg-9, Ser-10‚‚‚Arg14, Val-4‚‚‚Thr-8, Arg-13‚‚‚Thr-17, Ala-25‚‚‚Gln-22, Gln26‚‚‚Asp-30, Gln-29‚‚‚Asn-34, Leu-37‚‚‚Cys-41, and Arg14‚‚‚Gln-22. Domain Motions. We used three coordinates to measure relative motions of the domains, the interhelix angle, and distances between CR atoms of residues Tyr6 and Ala24 and between CR’s of residues Met1 and Leu52. Trajectory averages of these coordinates are presented in Table 2. Time evolutions of the interhelix angles and distances are shown in Figures 9 and 10. The average values of the interhelical angles vary between 54 and 130° in the four simulations.

The corresponding average interdomain distances were within a 18-26 Å range for Tyr 6‚‚‚Ala24 and 42-68 Å range for Met 1‚‚‚Leu52. For the three simulations of regular PLB, the average interhelix angles are 54° in methanol, 63° in water, and 130° in the membrane bilayer. Thus, the domains tend to remain close to perpendicular to each other. Additionally, there is a marked difference of domain orientations between solution and membrane environments. In the methanol and water solution simulations, the helices tend to be in more “closed” conformations, with the average interhelical angle below 90°, while in the membrane the helices assume a more “open” orientation, with the average interhelix angle of 130°. Additionally, the helix angle fluctuations are larger in the membrane environment as in the methanol and water solutions. We note that in solution and membrane the interhelix angle appears to fluctuate around a well-defined average value. The average interhelix distances in the PLB trajectories track the trend of the helix angles. The average Tyr 6‚‚‚Ala24 distances are 18-23 Å in the solution simula-

1786 Biochemistry, Vol. 44, No. 6, 2005

Houndonougbo et al.

FIGURE 8: Evolution of helicity of PLB, C- and N-terminal helix domains, during molecular dynamics simulations. All helix percentages were calculated with respect to the number of R-helix residues as defined in the PDB file of each starting structure. Table 1: Average Helicities of Phospholamban Domainsa

Table 2: Average Domain Angles and Distancesa

simulations

N-terminal helix (%)

C-terminal helix (%)

PLB in methanol PLB in water WSPLB in water PLB in DPPC

84 ( 5 81 ( 7 77 ( 3 86 ( 3

88 ( 6 90 ( 5 87 ( 2 85 ( 8

a For simulations, trajectory averages and standard deviations of data are shown. For details see Methods.

Met1‚‚‚

tions and 26 Å in the bilayer. The corresponding Leu52 average distances are 42-62 and 68 Å. Thus, larger distances roughly correlate with larger angles, consistent with the picture of the helical domains executing mainly rigidbody motions. WSPLB exhibits an average interhelix angle of 100°, and average distances of 23 and 63 Å for Tyr6‚‚‚Ala24 and Met1‚‚‚Leu52, respectively. Thus, this system tends to occupy an “open” conformation, intermediate between those seen for regular PLB in solution and membrane bilayer. The domain orientations and distances exhibit the largest fluctuations for this system (Table 2), consistent with the large atomic fluctuations (see OVerall Structure and Dynamics). The time evolution of the interhelix angles contains a large contribution from fluctuations on a time scale of about 1 ns,

simulations

interhelical angle (°)

Tyr-6‚‚‚Ala-24 distance (Å)

Res-1‚‚‚Res-52 distance (Å)

PLB in methanol PLB in water WSPLB in water PLB in DPPC NMRb

54 ( 13 63 ( 15 100 ( 27 130 ( 19 80 ( 22

18 ( 1 19 ( 1 23 ( 3 25 ( 2 26 ( 2

42 ( 5 47 ( 7 63 ( 8 68 ( 4 51 ( 9

a Angles are calculated using the helical limits defined in PDB Files 1FJK for solution and 1N7L for DPPC simulations. The distances are between residue CR atoms. The errors are the standard deviations of the data. b Calculated from the 20 structures (1N7L) determined in DPC micelle.6

as may be seen from the plots angle time series (Figure 9), as well as analysis of their power spectra. To study the motions of domains relative to the membrane plane, we calculated the angles between the PLB helices and the bilayer plane normal (Z axis) in the DPPC trajectory. As shown in Figure 11, the angle between the bilayer normal and the N-terminal helix varies between 75 and 177° with an average value of 132°. This corresponds to tilts relative to the membrane plane ranging from -15 to 87°, with an average of 42°. The angle between the C-terminal helix and

Computer Simulations of Phospholamban

FIGURE 9: Instantaneous interhelical angles of PLB and WSPLB. Angles are calculated using the helical limits defined in files 1FJK.pdb and 1N7L.pdb. The helix axis is the normalized vector defining the principal axis obtained from diagonalizing the moment of inertia tensor for the selected atoms.

the membrane plane normal varies in the 0-17° range with an average value of 8°. The domain orientations and distances found in our trajectories are in good agreement with experimental data. The interhelix angle calculated from our trajectory of PLB in methanol, 54 ( 14°, explores a similar range as the 68 ( 23° seen in NMR studies in organic solvents (3, 5). The interhelical angle of 130 ( 19° found in our DPPC bilayer simulation tends to be higher than experimentally observed values in lipid environments, 80 ( 22° in DPC micelles (6) and 80-120° range in lipid bilayers (7). The observed and calculated angle distributions exhibit partial overlap, especially with the bilayer data. The experimental trend of more open conformations in lipids relative to solution is reproduced by the simulations, although the calculated effect is larger. The mainly perpendicular orientation of the C-terminal helix relative to the membrane plane found in the DPPC simulation agrees with solid-state NMR results (7, 9). Fluorescence resonance energy transfer (FRET) experiments provide a distance of 21 ( 1 Å between fluorescent probes attached to residues 6 and 24 for PLB in lipid vesicles (50). This distance is in good agreement with our calculated Tyr6‚‚‚Ala24 distance of 26 ( 3 Å. Qualitatively, the large domain motions found in the simulations agree with the conformation heterogeneity of PLB observed in the FRET study (50). We note here that we follow the convention that helices aligned along one line, such as a continuous, extended R-helical structure of PLB, correspond to an interhelix angle of 180°. This is consistent with the angles reported for the solution and micelle NMR structures (3, 6). To simplify comparisons, we have converted the solid-state NMR results (7) to conform with this convention. Essential Dynamics. We have performed Essential Dynamics (ED) analysis of our trajectories to further describe

Biochemistry, Vol. 44, No. 6, 2005 1787 the concerted motions of PLB. As in previous studies (42, 51, 52), analysis of our MD trajectories has shown that most of the trajectory fluctuations can be described by a few largescale collective motions that correspond to the eigenvectors of the covariance matrix with largest eigenvalues, eq 1. The first three modes represent between 70 and 90% of the fluctuations in the different trajectories. To represent 95% of the fluctuations, it takes 7 modes of PLB in methanol, 6 modes of PLB in water, 5 modes of PLB in membrane, and 12 modes of WSPLB in water. Molecular graphics analysis of the ED modes with the largest eigenvalues shows the presence of three main kinds of collective motions of PLB. The first is a hinge bending motion in which flexing of the linker region changes the interhelix angle, with the two helical domains reorienting as rigid bodies. The second type involves twisting of the C-terminal and N-terminal domains around their helical axes. The third kind of motion is a flexing of the longer C-terminal helix, in which this domain fluctuates between a straight and curved conformation while retaining its helical structure. In the hinge bending and helix twisting motions, the two helical domains behave essentially as rigid bodies. In contrast, the C-terminal helix flexing mode represents an elastic deformation of the internal structure of that domain. This kind of deformation is present both in the solution and lipid bilayer simulations, although it has a lower amplitude in the latter environment. The helix flexing motion appears to be present only in the C-terminal helix, the longer of the two helical domains of PLB. This might be the reason the C-terminal domain tends to exhibit systematically larger structural fluctuations than the N-terminal domain, in all the modeled environments. Animations of selected ED modes may be found at our website (53). Contacts and Interactions with Lipids. To assess the interactions between PLB and the surrounding phospholipid molecules, we have estimated the number of interatomic contacts between the PLB N-terminal domain (residues 1-21) and the lipid headgroups. Most of these contacts involved residues Arg-14, Met-20, and Pro-21. These contacts existed during 9, 49, and 4% of the 10 ns trajectory for Arg-14, Met-20, and Pro-21, respectively. Thus, most of the interactions of the PLB N-terminal domain with the lipid headgroups involved Met-20, located in the interdomain linker. The contacts involved mainly the methionine side chain. Interdomain Linker. To analyze the dynamics of the interdomain linker (residues 17-20), we have calculated the time evolution of the backbone dihedrals of the residues involved, as well as the linker end-to-end distance (CR17...CR20 distance) and linker pseudodihedral (Cβ17-CR17-CR20-Cβ20 dihedral). The last coordinates are presented in Figure 12. The PLB interdomain linkers show significant flexibilities in the four simulations. The linker end-to-end distance varies in the range 5-6 Å in the solution simulations and 8-10 Å in DPPC. Thus, the more open interhelix angle in DPPC corresponds also to elongated structures of the linker. The linker pseudodihedrals fluctuate over a wide range of values, spanning almost the full possible orientation range. Interestingly, there are few structural transitions in the linker region. Thus, the linker appears to undergo fluctuations around a well-defined average structure, undergoing primarily elastic deformations, without internal structural changes.

1788 Biochemistry, Vol. 44, No. 6, 2005

Houndonougbo et al.

FIGURE 10: Instantaneous end-to-end distance (left) and Tyr-6‚‚‚Ala-24 distance (right) during the simulations of PLB and WSPLB.

FIGURE 11: Time evolution of the N- and C-terminal helix angle relative normal to the membrane plane.

SolVent Accessible Surface Areas. Using a probe radius of 1.6 Å, we computed the solvent accessible surface areas (SASA) of PLB in our simulations. The results for the whole protein, domains, and phosphorylation sites are compared with corresponding values from the starting NMR structures in Table 3. The average SASAs of whole PLB from trajectories in methanol, water, and DPPC were 5645, 5558, and 5641 Å2, respectively. Thus, the simulated structures in solution tend to have systematically larger SASAs than the NMR starting

coordinates, which exhibited total SASA of 5406 Å2 in organic solvent and 5066 Å2 in micelles. Similar effects have been seen previously in our simulations of calmodulin (24, 52, 54). As the domain internal structures sampled in the trajectories tend to be similar to the starting NMR coordinates in terms of RMSDs and secondary structure, this effect is probably due to subtle rearrangements of the molecular surface and may be traced to the use of vacuum force field calculations in NMR structure refinement (55). The whole protein tends to have a lower SASA in the lipid environment compared to organic solvent, due to the assumption of a more compact structure by the hydrophobic transmembrane domain. A similar trend is found in the experimental data, Table 3. In our simulations, the SASA of the N-terminal domain tends to be larger in the lipid environment than in solution. This is consistent with the more open conformation of the two helices in micelles and bilayers (see Table 2). For the engineered WSPLB in aqueous solution, the average total SASA was larger than for the regular PLB in the same environment. This is due to a larger SASA of the C-terminal helix of WSPLB compared to PLB. Such an effect is consistent with the fact that WSPLB was created by introducing several hydrophobic to hydrophilic mutations in the C-terminal helix of PLB (20). For the two phosphorylation sites, Ser16 and Thr17, the general trend is for the residues to remain solvent accessible, with larger SASA in the lipid environment relative to solution. This trend is in qualitative agreement with the experimental structures. The calculated SASAs of Ser and Thr in tripeptide models, Ala-Ser-Ala and Ala-Thr-Ala built in extended conformation with the program CHARMM, are 128 and 155 Å2, respectively. Thus, the trajectory average

Computer Simulations of Phospholamban

Biochemistry, Vol. 44, No. 6, 2005 1789

FIGURE 12: Evolution of CR17‚‚‚CR20 distance (left) and pseudodihedral Cβ17-CR17-CR20-Cβ20 (right) of the hinge region during the MD simulations. Table 3: Trajectory Average SASAs of Phospholamban Compared with Experimental Dataa simulations

whole protein

N-terminal domain (res 1-16)

C-terminal domain (res 21-52)

loop (res 17-20)

Ser16

Thr17

PLB in methanol PLB in water WSPLB in water PLB in DPPC NMRb NMRc

5645 ( 25 5558 ( 35 5760 ( 21 5641 ( 17 5406 5066

1725 ( 9 1711 ( 24 1708 ( 13 1881 ( 9 1652 1623

3445 ( 18 3392 ( 14 3641 ( 10 3319 ( 13 3351 3051

475 ( 4 455 ( 7 411 ( 6 441 ( 6 403 392

43 ( 3 35 ( 3 52 ( 4 94 ( 1 37 79

63 ( 3 67 ( 2 67 ( 6 102 ( 4 83 86

a

Units: Å2. b In organic solvent, structure 1FJK.3

c

In micelles, structure 1N7L.6

SASA in the DPPC bilayer correspond to 73% of maximum possible exposure of Ser16, and 65% for Thr17. Additionally, large SASA fluctuations are found for these sites over the course of the trajectories. In the water simulation, there were on average five water molecules within 4.5 Å of each of the side chain hydroxyl groups of Ser16 and Thr17. In the DPPC bilayer, these numbers were seven for Ser16 and four for Thr17. Rotational Diffusion. To analyze reorientational motions of PLB, we calculated the reorientation rates around the individual molecular axes and their average for the whole protein, as well as the rotational diffusion coefficients of the helical axis, D⊥, and around the helix axis, D|, for the two helical domains. The results are presented in Table 4. The average rotational diffusion coefficients of PLB were Dr ) 0.04 ns-1 in methanol and Dr ) 0.08 ns-1 in water, corresponding to the average (second-order) rotational relaxation time τr ) 1/6Dr of 4 and 2 ns, respectively. The diffusion coefficient values exhibit large statistical errors, because the simulation time is relatively short, comparable to τr. Thus, it is not productive to analyze the overall rotation results in too much detail. However, the calculated reorientation rates of the N-terminal and C-terminal domains of PLB in solution exhibit some interesting features that rise above the noise level.

Table 4: Reorientation Motions of PLB and WSPLB in Solutiona PLB in methanol

PLB in water

WSPLB in water

Dr, 1/ns τr, ns Dx, 1/ns Dy, 1/ns Dz, 1/ns τi, ns

Whole Protein 0.04 ( 0.01 0.08 ( 0.02 4(1 2(1 0.05 ( 0.04 0.18 ( 0.03 0.04 ( 0.04 0.05 ( 0.03 0.03 ( 0.04 0.01 ( 0.03 4-6 1-6

0.04 ( 0.03 4(2 0.11 ( 0.04 0.01 ( 0.04 0.01 ( 0.04 2-14

D⊥, 1/ns D||, 1/ns

C-Terminal Helix 0.04 ( 0.02 0.04 ( 0.01 1.9 ( 0.1 1.9 ( 0.3

0.08 ( 0.01 0.8 ( 0.5

D⊥, 1/ns D||, 1/ns

N-Terminal Helix 0.11 ( 0.03 0.27 ( 0.03 1.2 ( 0.4 0.5 ( 0.1

0.16 ( 0.03 0.3 ( 0.1

a The principal axis (x) is roughly along the C-terminal helix (see starting structure in Figure 2), the axis (y) is perpendicular to (x) in the figure plane, and the z axis is normal to the figure plane. The helix boundaries are taken as defined in the PDB files.

The C-terminal domain tumbling proceeds at a rate similar to overall protein reorientation, with diffusion coefficients D⊥ of 0.04-0.08 ns-1. The tumbling motion is somewhat faster in the N-terminal domain, with D⊥ of 0.11-0.27 ns-1. We note that the rate is faster in the more opened conforma-

1790 Biochemistry, Vol. 44, No. 6, 2005 tion. The time scale of these motions, τ⊥ ) 1/2D⊥, is 6-12 ns for the C-terminal and 4-6 ns for the N-terminal. The most interesting finding is that the rates for domain reorientations around the helix axes (spinning) D| were 0.5-1.9 ns-1, an order of magnitude faster than the tumbling. The time scale of the spinning motions, τ| ) 1/D| , is in the range of 0.5-2 ns. Faster rotation around the long axis of a freely moving prolate object is in accord with hydrodynamic theory (56). For WSPLB in aqueous solution, the rotational results are very similar to those for PLB. Overall tumbling of WSPLB proceeded at an average rate of Dr ) 0.04 ns-1, corresponding to a correlation time τr of 4 ns. The nanosecond spinning motions around the domain helical axes were also found in the WSPLB trajectory. Translational Diffusion. The translational diffusion coefficients Dt were 0.5 × 10-6 cm2/s for PLB in methanol, 1.2 × 10-6 cm2/s for PLB in water, and 3.6 × 10-6 cm2/s for WSPLB in water. These results are in the range expected for small proteins.56 For a rigid body moving in a continuous solvent, we would expect diffusion to be about 60% faster in methanol than in water due to the difference in viscosity (0.544 cP for methanol vs 0.890 cP for water at 300 K (57)). The opposite is found in the calculated values, with the diffusion coefficient of PLB in methanol being about 50% lower than in water. A possible reason for this effect is the exaggerated self-diffusion of the TIP3 water model used in our simulations (29). The simulation values of Dt correlate with the average interhelix angle (Table 2), exhibiting a trend for faster translational diffusion of systems with more open angles, i.e., more elongated structures. This trend is in qualitative agreement with hydrodynamic theory (56). DISCUSSION Our simulation results are consistent with a wide range of experimental data. The overall structure of PLB, in the form of an “L”-shaped arrangement of two R-helical domains, is conserved throughout the simulations, in agreement with the NMR starting structures (3, 6). The largescale domain motions, reflected in the variation of the interhelical angle, are also found in both the simulations and NMR data (3, 6). The distribution of interhelix angles for PLB in our methanol simulation, 54 ( 13°, exhibits a large overlap with the 68 ( 23° found from NMR studies in organic solvents (3). In the DPPC bilayer simulation, the interhelical angle distribution of 130 ( 19° was centered above the 80 ( 22° found in DPC micelles (6) and 80120° range found in lipid bilayers (7). This last discrepancy may be explained by the fact that the solid-state NMR signal was modeled assuming the helices of PLB are rigid cylinders (7). The simulations show significant flexing of the Cterminal domain, which should lead to a wider distribution of interhelical angles consistent with the NMR data on individual N-H bond orientations (7). The observed and calculated angle distributions exhibit partial overlap, especially in the case of the bilayer data, and the experimental trend of more open conformations in lipids relative to organic solvents is qualitatively reproduced by the modeling results. The mainly perpendicular orientation of the C-terminal helix relative to the membrane plane found in the DPPC simulation agrees with solid-state NMR results (7, 9).

Houndonougbo et al. Another measure of domain motions is the Tyr 6‚‚‚Ala24 distance of 26 ( 2 Å in the DPPC bilayer trajectory. This is in good agreement with the FRET determination of 21 ( 1 Å (50), given that the experimental data describe the distance of attached probes, while the simulations refer to the CR distance. Qualitatively, the large domain motions found in the simulations agree with the conformation heterogeneity of PLB observed in the FRET study (50). Another aspect on which the simulations and experimental data agree is the adoption of a more compact structure of the hydrophobic C-terminal domain in the lipid environment relative to organic solvent (3, 6). Although PLB is known to be insoluble in water and thus no experimental data are available in this case, we have simulated the dynamics of PLB in aqueous solution to provide a reference point for comparison. On the basis of the 5 ns trajectory, it appears that PLB retains its helical domains and “L” shape in water, with interhelix angles intermediate between those found in methanol and DPPC bilayer. Our results thus indicate that the insolubility of PLB in water might not be due to unfolding, but rather to the hydrophobic nature and aggregation of the transmembrane domain. Such aggregation was not possible in our simulations that involved a single PLB molecule. The simulation results for WSPLB tended to be generally similar to those found for the regular PLB system in solution and membrane bilayer. The WSPLB interhelix angle of 100 ( 27° falls in the range between those explored by PLB in water and membrane trajectories. Thus, the simulations predict that the structure and dynamics of WSPLB should be similar to PLB. This is in accord with experimental data, which show that WSPLB exhibits similar functional properties to the original protein (20). Interestingly, the helix angle fluctuations for WSPLB are the largest of all our four trajectories. Thus, the water-soluble model might be more labile than the regular PLB. Our modeling techniques involved molecular dynamics with the CHARMM force field, including the CMAP correction term for peptide backbone, with standard nonbond cutoffs for van der Waals forces and Ewald summations of electrostatic interactions. Thanks to the good overall agreement with experimental data, we believe that our simulations provide a reliable description of the behavior of PLB in different environments. Our most interesting results are the identification of several kinds of collective motions of PLB that occur on the nanosecond time scale. The major modes of collective, large-amplitude motion involve hinge bending, twisting of both helical domains, and flexing of the Cterminal helix. The first two kinds of motions of the helical domains found in the simulations can clearly have functional roles. The presence of the hinge bending motion leads to large fluctuations in the interhelix angle, and thus the population of elongated conformations which allow the N-terminal helix to come into contact with the PLB binding site on the calcium ATPase (58). The twisting motions around the helical axis enable the N-terminal helix to orient the correct face to its binding site on the ATPase. The flexing and twisting of the C-terminal domain, which is found in both solution and bilayer environments, may be related to its ability to aggregate into stable pentamers, which have been modeled as coiled-coils (18). Interestingly, the central linker region of PLB tends to retain its turn structure in the

Computer Simulations of Phospholamban trajectories. Thus, the large-scale helical domain motions mostly involve elastic deformations of this region rather than structural transitions. Another interesting effect are the significant trajectory average values and large fluctuations of solvent accessibility of the phosphorylation sites Ser16 and Thr17. Thus, these sites are essentially always accessible to solvent, and presumably also to phosphorylation. Overall, we were able to successfully model PLB in organic solvent and bilayer environments. This opens the way to simulating other processes involving PLB, such as phosphorylation, aggregation, and interaction with the calcium ATPase. CONCLUSIONS We present a microscopic description of the structure and dynamics of PLB in solution and membrane environments based on molecular dynamics simulations. The simulations, with 5 ns length in solution and 10 ns length in the membrane, show the presence of two main features of PLB structure and dynamics: the presence of two well-defined helical domains at the N- and C-termini, and large-amplitude rigid-body motions of these domains. In all our simulations, the N-terminal and C-terminal domains of PLB remained close to their experimental starting coordinates, with a high content of R-helical structure. In PLB trajectories in methanol and water, the two helices were almost perpendicular, with average interhelix angles of 54 ( 13° and 63 ( 15°, respectively. In the lipid bilayer trajectory, both the interhelix angle and its fluctuations were larger, with an average of 130 ( 19°. As the transmembrane, C-terminal, helix remained approximately perpendicular to the bilayer plane, the N-terminal helix tended to be in an open, or extended conformation, with an average tilt of 132° with respect to the bilayer normal. The large amplitude collective motions of the helical domains fell into three types: hinge bending, twisting of both N- and C-terminal helices, and flexing of the C-terminal helix. The central linker of PLB exhibited high flexibility, due mostly to elastic deformations and not to any significant structural transitions. Many features observed in the simulations are in agreement with experimental data. These include the mostly helical secondary structure of the protein, the range explored by the interhelix angle in methanol, as well as the interhelix distance and C-terminal helix orientation in the DPPC bilayer. The calculated interhelix angles in the DPPC bilayer tend to be larger than observed values in lipid environments. However, the simulation and experimental data exhibit significant overlap, and the observed effect of opening up of the PLB interhelix angle in the lipid environment relative to solution is qualitatively reproduced in the simulations. Thanks to the good overall agreement with a wide range of experimental data, we believe that our simulations provide a reliable description of the behavior of PLB in solution and lipid bilayer environments. This allows us to draw some conclusions about the microscopic details of PLB structure and dynamics that may be difficult to obtain by experimental approaches alone, and to relate these details to PLB biological function. Two kinds of motions of the helical domains found in the simulations can clearly have functional roles. The population of conformations with relatively open interdomain

Biochemistry, Vol. 44, No. 6, 2005 1791 angles, as well as large fluctuations of this coordinate in DPPC bilayers, allow the N-terminal helix to come into contact with the PLB binding site on the calcium ATPase (58). Additionally, the presence of twisting motions around the helical axis enables the helix to orient the correct face to the binding site. Generally, the presence of the flexible central linker is crucial for the correct positioning and orientation of the N-terminal helix for its regulatory interactions with the calcium ATPase. Another interesting effect is that the phosphorylation sites Ser16 and Thr17 are essentially always accessible to solvent, and presumably also to phosphorylation. Overall, we were able to successfully model PLB in organic solvent and bilayer environments. This opens the way to modeling other processes involving PLB, such as phosphorylation, aggregation, and interaction with the calcium ATPase. We are currently pursuing molecular dynamics simulations of pentameric and phosphorylated PLB to understand how aggregation and phosphorylation modify the properties of this protein to modulate its interactions with the calcium ATPase. REFERENCES 1. Simmerman, H. K. B., Jones, L. R., and Amgen, T. 0. (1998) Phospholamban: protein structure, mechanism of action, and role in cardiac function, Physiol. ReV. 78, 921. 2. Brittsan, A. G, and Kranias, E. G. (2000) Phospholamban and cardiac contractile function, J. Mol. Cell. Cardiol. 32, 2131. 3. Lamberth, S., Schmid, H., Muenchbach, M., Vorherr, T., Krebs, J., Crafoli, E., and Griesinger, C. (2000) NMR solution structure of phospholamban, HelV. Chim. Acta 83, 2141. 4. Maslennikov, V., Sobol, A. G., Anagli, J., James, P., Vorherr, T., Arseniev, A. S., and Carafoli, E. (1995) The secondary structure of phospholamban: a two-dimensional NMR study, Biochem. Biophys. Res. Commun. 217, 1200. 5. Pollesello, P., Annila, A., and Ovaska, M. (1999) Structure of the 1-36 amino-terminal fragment of human phospholamban by nuclear magnetic resonance and modeling of the phospholamban pentamer, Biophys. J. 76, 1784. 6. Zamoon, J., Mascioni, A., Thomas, D. D., and Veglia, G. (2003) NMR solution structure and topological orientation of monomeric phospholamban in dodecylphosphocholine micelles, Biophys. J. 85, 2589-2598. 7. Mascioni, A., Karim, C., Zamoon, J., Thomas, D. D., and Veglia, G. (2002) Solid-state NMR and rigid body molecular dynamics to determine domain orientations of monomeric phospholamban, J. Am. Chem. Soc. 124, 9392. 8. Tatulian, S. A., Jones, L. R., Reddy, L. G., Stokes, D. L., and Tamm, L. K. (1995) Secondary structure and orientation of phospholamban reconstituted in supported bilayers from polarized attenuated total reflection FTIR spectroscopy, Biochemistry 34, 4448. 9. Toyoshima, C., Asahi, M., Sugita, Y., Khanna, R., Tsuda, T., and MacLennan, D. H. (2003) Modeling of the inhibitory interaction of phospholamban with the Ca2+ atpase, Proc. Natl. Acad. Sci. U.S.A. 100, 467. 10. . Arkin, T., Rothman, M., Ludlam, C. F., Aimoto, S., Engelman, D. M., Rothschild, K. J., and Smith, S. O. (1995) Structural organization of the pentameric transmembrane alpha-helices of phospholamban, a cardiac ion channel, J. Mol. Biol. 248, 824834. 11. Jones, L. R. (1985) Sarcolemmal enzymes mediating 3-adrenergic effects on the heart, Curr. Top. Membr. Transp. 25, 11-41. 12. Wegener, A. D., and Jones, L. R. (1984) Phosphorylation-induced mobility shift in phospholamban in sodium dodecyl sulfatepolyacrylamide gels. Evidence for a protein structure consisting of multiple identical phosphorylatable subunits, J. Biol. Chem. 259, 1834-1841. 13. Fuji, J., Adoma, M., Tada, M., Tada, H., and Sakiyama, F. (1986) Characterization of structure unit of phospholamban by amino acid sequencing and electrophoretic analysis, Biochem. Biophys. Res. Commun. 138, 1044-1050.

1792 Biochemistry, Vol. 44, No. 6, 2005 14. Simmerman, H. K. B., Collins, J. H., Theibert, J. L., Wegener, A. D., and Jones, L. R. (1986) Sequence analysis of phospholamban. Identification of phosphorylation sites and two major structural domains, J. Biol. Chem. 261, 13333-13341. 15. I.. Arkin, T., Adams, P. D., MacKenzie, K. R., Lemmon, M. A., Brunger, A. T., and Engelman D M. (1994) Structural organization of the pentameric transmembrane alpha-helices of phospholamban, a cardiac ion channel, EMBO J. 13, 4757-4764. 16. Wegener, A. D., Simmerman, H. K., Liepnieks, J., and Jones, L. R. (1986) Proteolytic cleavage of phospholamban purified from canine cardiac sarcoplasmic reticulum vesicles. Generation of a low resolution model of phospholamban structure, J. Biol. Chem. 261, 5154-5159. 17. Torres, J., Adams, P. D., and Arkin, I. T. (2000) Use of a new label, 13Cd18O, in the determination of a structural model of phospholamban in a lipid bilayer. Spatial restraints resolve the ambiguity arising from interpretations of mutagenesis data, J. Mol. Biol. 300, 677-685. 18. H. K. B., Heather, K. B., Kobayashi, Y. M., Autry, J. M., and Jones, L. R. (1996) A leucine zipper stabilizes the pentameric membrane domain of phospholamban and forms a coiled-coil pore structure, J. Biol. Chem. 271, 5941-5946. 19. Karim, C., Stamm, J., Karim, J., Jones, L., and Thomas, D. (1998) Cysteine reactivity and oligomeric structures of phospholamban and its mutants, Biochemistry 37, 12074-12081. 20. Slovic, A. M., Summa, C. M., Lear, J. D., and DeGrado, W. F. (2003) Computational design of a water-soluble analog of phospholamban, Protein Sci. 12, 337-348. 21. Frank, S., Kammerer, R. A., Hellstern, S., Pegoraro, S., Stetefeld, J., Lustig, A., Moroder, L., and Engel, J. (2000) Toward a highresolution structure of phospholamban: Design of soluble transmembrane domain mutants, Biochemistry 39, 6825-6831. 22. Li, H., Cocco, M. J., Steitz, T. A., and Engelman, D. M. (2001) Conversion of phospholamban into a soluble pentameric helical bundle, Biochemistry 40, 6636-6645. 23. Mahadevan, J., Xu, C., and Kuczera, K. (2002) Molecular dynamics simulations of conformational behavior of linear RGD peptidomimetics and cyclic pro- drugs in aqueous and octane solutions, J. Biomol. Struct. Dyn. 19, 775. 24. Yang, C., Jas, G. S., and Kuczera, K. (2004) Structure, dynamics and interaction with kinase taxgets: computer simulations of calmodulin, Biochim. Biophys. Acta 1697, 289-300. 25. Wong, T. C. (2003) Membrane structure of the human immunodeficiency virus gp4l fusion peptide by molecular dynamics simulation. ii. The glycine mutants,. Biochim. Biophys. Acta 1609, 45-54. 26. Faraldo-Gomez, J. D., Smith, G. R., and Sansom, M. S. (2003) Molecular dynamics simulations of the bacterial outer membrane protein fhua: a comparative study of the ferrichrome-free and bound states, Biophys J. 85, 1406-1420. 27. Brooks, B. R., Bruccoleri, R., Olafson, B., States, D., Swaminathan, S., and Karplus, M. (1983) CHARMM: A program for macromolecular energy, minimization and dynamics calculations, J. Comput. Chem. 4, 187-217. 28. Grubmueller, H., Ed., SOLVATE 1.0 program, Munich, LudwigMaximilians University, 1996. 29. Jorgensen, W. L., Chandrasekhar, J., Impey, R. W., Madura, J. D., and Klein, W. L. (1983) Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79, 926. 30. Hoover, W. G. (1985) Phys. ReV. A 31, 1695. 31. Feller, S. E., Zhang, Y., Pastor, R. W., and Brooks, B. R. (1995) J. Chem. Phys. 103, 4613-4621. 32. Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. C. (1977) Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes, J. Comput. Phys. 23, 327. 33. Barth, E., Kuczera, K., Leimkuhler, B., and Skeel, R. D. (1995) Algorithms for constrained molecular dynamics, J. Comput. Chem. 16, 1192-1209. 34. Woolf, T. B., and Roux, B. (1994) Molecular dynamics simulation of the gramicidin channel in a phospholipid bilayer, Proc. Natl. Acad. Sci., U.S.A. 91, 11631-11635. 35. Woolf, T. B., and Roux, B. Molecular dynamics of PfI coat protein in a phospholipid bilayer in Biological Membranes: A Molecular PerpectiVe from Computation and Experiment (Merz, K. M., and Roux, B., Eds.) Berkhauser, Boston, 1996. 36. Feller, S. E., Venable, R. M., and Pastor, R. W. (1997) Computer simulation of a dppc phospholipid bilayer: Structural changes as a function of molecular surface area, Langmuir 13, 6555.

Houndonougbo et al. 37. Venable, R. M., Zhang, Y., Hardy, B. J., and Pastor, R. W. (1993) Molecular dynamics simulations of a lipid bilayer and of hexadecane: an investigation of membrane fluidity, Science 262, 223226. 38. Frenkel, D., and Smit, B. Understanding Molecular Simulation, Academic Press, New York, 1996. 39. MacKerell, A. D., Jr., Bashford, D., Bellott, R. L., Dunbrack, R. L., Jr., Evanseck, J. D., Field, M. J., Fischer, S., Gao, J., Guo, H., Ha, S., Joseph-McCaxthy, D., Kuchnir, L., Kuczera, K., Lau, F. T. K., Mattos, C., Michnick, S., Ngo, T., Nguyen, D. T., Prodhom, B., Reiher, III, W. E., Roux, B., Schlenkrich, M., Smith, J. C., Stote, R., Straub, J., Watanabe, M., Wiorkiewicz-Kuczera, J., Yin, D., and Karplus, M. (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins, J. Phys. Chem. B. 102, 586-616. 40. MacKerell, A. D., Feig, M., and Brooks, C. L. (2004) Improved treatment of the protein backbone in empirical force fields, J. Am. Chem. Soc. 126, 698-699. 41. Kabsch, W., and Sander, C. (1983) Dictionary of protein secondary structure: Pattern recognition of hydrogen-bonded and geometrical features, Biopolymers 22, 2577. 42. Garcia, A. E. (1992) Nonlinear dynamics of proteins, Phys. ReV. Lett. 68, 2696-2699. 43. Amadei, A., Linssen, A. B. M., Van Aalten, D. M. F., and Berendsen, H. J. C. (1996) An efficient method for sampling the essential subspace of proteins, J. Biomol. Struct. Dyn. 13, 615626. 44. Andricioaei, I., and Karplus, A. (2001) J. Chem. Phys. 115, 62896292. 45. McQuarrie, D. A. Statistical Mechanics, Harper and Row, New York, 1976. 46. Ehrenberg, M., and Rigler, R. (1972) Polarized fluorescence and rotational brownian motion, Chem. Phys. Lett. 14, 539-544. 47. Chuang, T. J., and Eisenthal, K. B. (1972) Theory of fluorescence depolarization anisotropic rotational diffusion, J. Chem. Phys. 57, 5094-5097. 48. Aqvist, J. (1986) Comput. Chem. 10, 97-99. 49. Shoemaker, D. P., Garland, C. W., and Nibler, J. W. Experiments in Physical Chemistry, 6th ed., The McGraw-Hill Companies, Inc., New York, 1996. 50. Li, J., Bigelow, D. J., and Squier, T. C. (2003) Phosphorylation by cAMP-dependent protein kinase modulates the structural coupling between the transmembrane and cytosolic domains of phospholamban, Biochemistry 42, 10674-10682. 51. Berendsen, J. C. (1999) Molecular dynamics simulations: The limits and beyond in Computational Molecular Dynamics: Challenges, Methods, Ideas (Deuflhard, P., Hermans, J., Leimkuhler, B., Mark, A. E., and Reich, S., and Skeel, R. D. Eds.) pp 3-26, Springer, Berlin. 52. Yang, C., Jas, G. S., and Kuczera, K. (2001) Structure and dynamics of calcium-activated calmodulin in solution, J. Biomol. Struct. Dyn. 19, 247-271. 53. Houndonougbo, Y. A., and Kuczera, K. (2004) Images and animations of phospholamban dynamics in solution an membrane bilayer. http://129.237.102.17/ kuczera/plb.html. 54. Yang, C., and Kuczera, K. (2002) Molecular dynamics simulations of a calmodulin-peptide complex in solution, J. Biomol. Struct. Dyn. 20, 179-197. 55. Xia, B., Tsui, V., Case, D. A., Dyson, H. J., and Wright, P. E. (2002) Comparison of protein solution structures refined by molecular dynamics simulations in a vacuum, with a generalized Born model and with explicit water, J. Biomol. NMR 22, 317331. 56. Cantor, C. R., and Schimmel, P. R. Biophysical Chemistry, Vol. H. Techniques for the Study of Biological Structure and Function, W. H. Freeman, New York, 1980. 57. Riddick, J. A., Bunger, W. B., and Sakand, T. K. Organic SolVents: Physical Properties and Methods of Purification, Wiley, New York, 1986. 58. Toyofuku, T., Kurzydlowski, K., Tada, M., and MacLennan, D. H. (1994) Amino acids Glu2 to Ile18 in the cytoplasmic domain of phospholamban are essential for functional association with the Ca(2+)-atpase of sarcoplasmic reticulum, J. Biol. Chem. 269, 3088-3094. BI0488404