This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Huber, T. Articles by Brown, M. F. Search for Related Content PubMed PubMed Citation Articles by Huber, T. Articles by Brown, M. F.

Membrane Model for the G-Protein-Coupled Receptor Rhodopsin: Hydrophobic Interface and Dynamical Structure

Thomas Huber ^* , Ana V. Botelho , Klaus Beyer  and Michael F. Brown  

^* Laboratory of Molecular Biology and Biochemistry, Howard Hughes Medical Institute and Rockefeller University, New York, New York 10021 USA; Department of Chemistry and Department of Biochemistry & Molecular Biophysics, University of Arizona, Tucson, Arizona 85721 USA; and Lehrstuhl für Stoffwechselbiochemie der Universität München, D-80336 Munich, Germany

Correspondence: Address reprint requests to Thomas Huber, Tel.: 212-327-8676; Fax: 212-327-7904; E-mail: hubert{at}mail.rockefeller.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION COMPUTATIONAL AND EXPERIMENTAL... RESULTS DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Rhodopsin is the only member of the pharmacologically important superfamily of G-protein-coupled receptors with a known structure at atomic resolution. A molecular dynamics model of rhodopsin in a POPC phospholipid bilayer was simulated for 15 ns, revealing a conformation significantly different from the recent crystal structures. The structure of the bilayer compared with a protein-free POPC control indicated hydrophobic matching with the nonpolar interface of the receptor, in agreement with deuterium NMR experiments. A new generalized molecular surface method, based on a three-dimensional Voronoi cell construction for atoms with different radii, was developed to quantify cross-sectional area profiles for the protein, lipid acyl chains and headgroups, and water. Thus, it was possible to investigate the bilayer deformation due to curvature of the individual lipid monolayers. Moreover, the generalized molecular surface derived hydrophobic interface allowed benchmarking of the hydropathy sequence analysis, an important structural genomics tool. Five water molecules diffused into internal hydration sites during the simulation, yielding a total of 12 internal waters. The cytoplasmic loops and the C-terminal tail, containing the G-protein recognition and protein sorting sequences, exhibited a high mobility, in marked contrast to the extracellular and transmembrane domains. The proposed functional coupling of the highly conserved ERY motif to the lipid-water interface via the cytoplasmic loops provides insight into lipid effects on G-protein-coupled receptor activation in terms of a flexible surface model, involving the spontaneous monolayer curvature.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL AND EXPERIMENTAL... RESULTS DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Heterotrimeric guanine nucleotide-binding protein (G-protein)-coupled receptors (GPCRs) constitute a large family of sensory proteins, which respond to a variety of different external stimuli such as hormones, neurotransmitters, neuromodulators, odors, and light (Strader et al., 1995; Ballesteros et al., 2001). They act as the interface between a signal and the cellular response mediated by molecular switches, the family of regulatory G-proteins. Cell membrane receptors, many of them GPCRs, are pharmacologically very important, as they comprise 45% of all molecular drug targets of current therapies (Drews, 2000). Rhodopsin is the only GPCR whose structure is known at atomic resolution (Palczewski et al., 2000); it is a member of the largest receptor Family A, with 569 possible candidates identified in the human genome (Lander et al., 2001). The opsins and the protein, peptide, and cationic amine receptor subfamilies, constitute 90% of all known GPCRs (Palczewski et al., 2000). They share a conserved structure of seven transmembrane (TM) helices, and the sequences of many other GPCRs are used to build homology models based on rhodopsin (Baldwin et al., 1997; Lomize et al., 1999; Stenkamp et al., 2002). Light-induced photoisomerization of the covalently bound ligand of rhodopsin, 11-cis-retinal, to all-trans-retinal is followed by a series of thermal relaxations known in kinetic detail (Lewis and Kliger, 2000), leading to a tertiary structural rearrangement of the protein (Hubbell et al., 2003). With respect to conformational transitions there is a high similarity among the different GPCRs (Ghanouni et al., 2001). The subsequent steps of the rhodopsin-dependent signal transduction cascade (Ebrey and Koutalos, 2001) involve binding of the heterotrimeric G-protein transducin, resulting in GDP–GTP exchange, and dissociation of the subunits. Activation of a phosphodiesterase by the G·GTP subunit is followed by subsequent enzymatic processes, which regulate the cytosolic levels of second messengers. These lead ultimately to closing of cation channels and hyperpolarization of the cellular plasma membrane, giving the visual signal. The structural transitions of the primary photoproduct of rhodopsin culminate in the signaling state, involving the metastable intermediates metarhodopsin I (meta I) and metarhodopsin II (meta II), where the latter is the activated receptor conformation (Hofmann et al., 1995).

The thermodynamic equilibrium between the meta I and meta II conformations of photoactivated rhodopsin is strongly dependent on the lipid composition of the membrane (Beach et al., 1984), and a native-like lipid mixture results in optimal function (Wiedmann et al., 1988). Evidence has increasingly accumulated that the functions of other membrane proteins are similarly modulated by the lipid composition (Jensen and Schutzbach, 1984; Navarro et al., 1984; Epand, 1990; Lee, 1998; Perozo et al., 2002). One possibility entails specific protein and lipid interactions that are related to the domain organization of biomembranes (rafts) (Simons and Ikonen, 2000). Yet in case of rhodopsin the membrane lipids are entirely in the fluid, liquid-crystalline (L) state at physiological temperature (Deese et al., 1981; Miljanich et al., 1985; Ellena et al., 1986). An alternative is that nonspecific lipid/protein interactions may influence function due to biophysical properties of the membrane lipids, as described by a flexible surface model (FSM) (Gibson and Brown, 1993; Botelho et al., 2002). It has been argued that the microscopic organization of rod outer segment membranes leads to an active control or regulation of rhodopsin function by the local lipid environment (Beach et al., 1984; Baldwin and Hubbell, 1985; Wiedmann et al., 1988; Mitchell et al., 1992; Brown, 1997). Therefore, the standard fluid-mosaic model for cellular membranes, in which fluid-like phospholipid bilayers constitute a matrix for integral and peripheral membrane proteins, has been extended to account for modulation of membrane protein function due to the composition of the lipid matrix (Botelho et al., 2002). Such lipid regulation of protein function has clear implications for human diseases, as revealed by studies of the influences of dietary intake of essential and polyunsaturated fatty acids (Neuringer, 2000), and the well-known effects of cholesterol in cardiovascular dysfunction and various neurological disorders.

This study centers on a biomembrane model of a protein in a native-like phospholipid bilayer environment under physiological conditions, comprising dark-state rhodopsin and POPC at full hydration. The molecular dynamics (MD) model was based on the crystallographic membrane protein structure (Palczewski et al., 2000; Teller et al., 2001; Okada et al., 2002), supplemented with experimental solid-state ²H NMR data for the phospholipid component. The microscopic structure of the rhodopsin-containing membrane model is founded on a 15-ns molecular dynamics computer simulation in an all-atom representation (Fig. 1), and is interpreted in terms of the hydrophobic protein/lipid interface. The MD simulations show that once rhodopsin is liberated from crystal-packing constraints and embedded in a native-like phospholipid bilayer environment, it relaxes to an ensemble of conformations different from the crystal structure. Information is obtained for additional internal water molecules and for parts of the protein that are crystallographically unresolved (Palczewski et al., 2000; Teller et al., 2001; Okada et al., 2002). The definition of a new generalized molecular surface (GMS) approach has allowed us to quantify the contact interfaces between the molecular components of a biomembrane simulation in terms of hydropathy analysis. The average lipid monolayer curvature is inferred to be negative, which is relevant to the native membrane environment of rhodopsin, with a high degree of nonlamellar, reverse hexagonal (H_II) phase-promoting phospholipids (Deese et al., 1981). We also suggest that association of the membrane-water interface with the cytoplasmic loop C2 and the adjacent ERY sequence motif of rhodopsin can impose stress/strain on the protein, thereby affecting its function. A remarkable stability of the transmembrane (TM) helix bundle in the bilayer is found in comparison to the extramembranous loops, which pertains to the physiological role of the receptor as a single quantum detector (Baylor, 1996). Finally, the current molecular dynamics model serves as a platform for future structure-function studies of rhodopsin and other homologous GPCRs (Huber et al., 2003). Preliminary results of this study have been communicated (Huber et al., 2002a).

View larger version (124K): [in this window] [in a new window]   FIGURE 1  Simulated system demonstrating the location of rhodopsin in the bilayer. A side view of rhodopsin in a POPC bilayer membrane is shown (protein/lipid ratio 1:103); two adjacent unit cells indicate the protein with the extracellular (intradiskal) side up and the cytosolic side down. Transmembrane helices are depicted as rods, and the extramembranous loops as strings. White spheres indicate the lipid acyl carbonyl carbons representing the hydrophobic boundary, and the phosphates representing the headgroups are in orange. Red and blue spheres illustrate the ionic atmosphere of the aqueous region, corresponding to Na+ and Cl-, respectively. The palmitoylation of Cys322 and Cys323 is shown by the white space-filling chains penetrating deep into the lower monolayer. The retinal chromophore is depicted in red. A total of 45,000 atoms was used to describe the system, which was simulated for 15 ns using a 1-fs time step and full electrostatic calculations (cf. text). The system is shown at an early state of the simulations, after 3 ns of equilibration time, in which the cytoplasmic monolayer is not yet fully adapted to the protein.

	COMPUTATIONAL AND EXPERIMENTAL METHODS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL AND EXPERIMENTAL... RESULTS DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

The MD simulation involved an all-atom model of the phospholipid/rhodopsin system (including the covalently bound retinal, the N-terminal acetylation and glycosylations, the C-terminal palmitoylations, and the disulfide bridge) with excess water and physiological salt concentration. To account for the experimentally determined protonation state of specific amino acids, changes of their ionization were performed during refinement of the protein models, as summarized in Table 1. The models were built on the basis of the recent crystal structures of rhodopsin (Palczewski et al., 2000; Teller et al., 2001), in which unresolved parts were added by distance geometry modeling methods (Blaney et al., 1995; Huber et al., 2002b). Briefly, missing residues of the crystal structure (cf. Table 1) were built in an extended conformation. The atoms represented in the crystal structure were treated in a rigid body approximation, fixing all internal degrees of freedom. The preprocessing algorithm of DGEOM95 (Blaney et al., 1995) generated a distance matrix of upper and lower distance limits from the covalent geometry and chirality of the fragments and their atomic van der Waals radii. Additional constraints were provided, describing the covalent bonds linking the fragments and the crystallographic-defined rhodopsin polypeptide chain. The distance geometry calculations generated a set of 20 protein models scored by pseudoenergy values derived from the distance constraints. The lowest energy conformations were free of constraint violations and were inspected for proper geometry and packing. Subsequently, the models were energy minimized with the CHARMM potential energy function, and inserted in an equilibrated bilayer model of POPC. The initial model of the POPC bilayer for the rhodopsin-containing simulation was generated by duplicating a previous, 72-lipid-containing POPC bilayer model (Huber et al., 2002b). Some of the water molecules and 20 lipids in each leaflet of the bilayer, whose excluded volumes overlapped with the rhodopsin molecule, were removed. Subsequently another lipid was taken out of the extracellular side, giving a total of 103 lipids per rhodopsin. The resulting defects in the bilayer structure were annealed with all the rhodopsin atoms constrained to fluctuate about the initial structure. Multiple time-step integration of the equations of motion allowed for a fully flexible covalent bond treatment with a 1.0-fs time step for the fastest degrees of freedom. The thermodynamic ensemble was a constant particle number, normal pressure, surface tension, and temperature (NPT) ensemble applying Nosé-Hoover dynamics (Feller and Pastor, 1996). The target pressure was 1.0 bar, and target surface tensions for the simulations are given in the Appendix. As a consequence of the periodic boundary conditions, the systems were oriented multibilayer stacks of the pure lipid or the rhodopsin-containing bilayer at full hydration in the liquid-crystalline state. The temperature for all simulated systems was 27°C (300 K), with one exception noted in Table 2. To account for the long-range nature of the Coulombic interactions, a full electrostatics calculation with particle-mesh Ewald lattice summation (Darden et al., 1993) was applied. The three-dimensional (3-D) fast Fourier transform grid size of 80 x 80 x 96 for the rhodopsin-containing systems replaced the 64 x 64 x 64 grid for the reference POPC bilayer system. All other parameters for the MD integration were identical to previous membrane simulations (Huber et al., 2002b).

Solid-state deuterium NMR spectroscopy
Deuterium NMR spectroscopy of the lipid dispersions at 11.7 T employed a Bruker AMX-500 spectrometer using the solid-state quadrupolar echo method with composite pulses, as described (Huber et al., 2002b). Briefly, the ²H NMR spectra of the randomly oriented bilayer membrane dispersions of POPC-d₃₁ in the presence or absence of rhodopsin were numerically deconvoluted (de-Paked). The C–²H segmental order parameters were calculated from the observed quadrupolar splittings, using a static quadrupolar coupling constant of 170 kHz for the aliphatic C–²H bonds. The assignments were made in reversed order from the methyl terminus toward the plateau region. The effect of rhodopsin on the lipid order parameter profiles was investigated by ²H NMR at pH 7 over the temperature range from -43 to 37°C for each sample, and the data determined at 27°C were compared with the MD simulations at the same temperature.

	RESULTS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL AND EXPERIMENTAL... RESULTS DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
All-atom molecular dynamics model of membrane-bound rhodopsin
In what follows, we present a refined model for dark-adapted rhodopsin in the native-like environment of a bilayer membrane comprising POPC in the liquid-crystalline (L) state at full hydration. The model is based on a 15-ns molecular dynamics computer simulation, utilizing an all-atom representation of the system. During the first 5 ns of the MD simulations, several modifications were applied to the initial model to obtain an accurate description of the real system. In particular, care had to be been taken to choose the correct protonation state of titratable amino acids, to solvate internal cavities with water, to adjust the number of lipids in each of the monolayers, and to find an area per lipid that reproduces the experimental observables. The model was equilibrated for a further 5 ns, and the final 5 ns of the 15-ns trajectory, designated segment P6 (cf. Appendix) were used for the analysis of the protein structure, the lipid order parameters, and the generalized molecular surface of the various interfaces of the biomembrane system. Likewise, the protein-free POPC bilayer was simulated for 39 ns, and a5-ns segment, designated L4, with a similar area per lipid was used for comparison of the order parameters. The details of the protein models are described in Table 1, and a complete description of the bilayer systems and equilibration can be found in the Appendix.

Comparison of calculated and experimental order parameter profiles for rhodopsin/POPC membranes
Fig. 2 a shows the results of theoretical MD simulations for the sn-1 acyl chain order parameters of POPC in the rhodopsin-containing membrane model (segment P6), and the protein-free lipid control (segment L4) (cf. Appendix). The average area per lipid for the rhodopsin simulation was 62.6 Ų, with some uncertainty introduced by the vase-shape protein inclusion (vide infra). The lipid control simulation segment was performed at a constant area per lipid of 62.4 Ų. As can be seen in Fig. 2 a, the orientational order of the lipid acyl chain segments C₈–C₁₅ is reduced due to rhodopsin incorporation in the MD model system. How do the MD simulations correspond to the experimental effects of rhodopsin on POPC acyl chain order parameters? In this regard, solid-state deuterium nuclear magnetic resonance (²H NMR) allows investigation of the bilayer structure (Petrache et al., 2001) under physiological conditions (hydration, ionic strength, and temperature). The acyl chain order parameters obtained by ²H NMR spectroscopy of deuterated POPC-d₃₁ bilayers in the presence of rhodopsin can be compared with the order parameters calculated from the MD simulations (Huber et al., 2002b). The order parameter profiles are sensitive to the average bilayer structure (Petrache et al., 2001), and can help to guide the choice of simulation parameters, such as the applied surface tension, as well as judge the quality of the final MD model.

View larger version (21K): [in this window] [in a new window]   FIGURE 2  Comparison of POPC acyl chain order parameters from MD simulations and solid-state deuterium (2H) NMR experiments. (a) Simulated order parameter profile for sn-1 chain of POPC at 27°C: (red --) presence of rhodopsin at zero surface tension (simulation segment P6; cf. Table 2) with an average cross-sectional area per lipid in the acyl chain region of 62.6 Å2; (black --) order parameter profile in absence of rhodopsin at a constant area per lipid of 62.4 Å2 (simulation segment L4). (b) Experimental 2H NMR sn-1 chain order parameter profile of POPC-d31 in the presence (red --) and absence (black --) of rhodopsin at 27°C and pH 7, with a lipid-to-protein ratio of 100:1. Insert: numerically deconvoluted solid-state 2H NMR spectra of POPC-d31 in the absence (black line) and presence of rhodopsin (red line). The order parameters nearest the aqueous interface show little effect of rhodopsin, whereas the carbon segments C8–C15 are more disordered.

Formulation of a generalized molecular surface based on Voronoi polyhedra
Currently there is no generally accepted standard method to describe the molecular packing of the protein, lipids, and water in biomembrane simulations. Hence, we developed a novel volumetric method based on a generalization of Voronoi polyhedra to differently sized atoms (Richards, 1985) to compute the molecular cross-sectional area and the solvent contact interface, the generalized molecular surface (cf. Appendix). The GMS quantifies the interface of rhodopsin with lipid and water, and allows for further analysis in terms of a hydropathy profile. One can then consider how the shape of rhodopsin in the membrane corresponds to the partial cross-sectional areas of the lipids (cf. Fig. 3), that is to say, how the transmembrane protein is embedded in the lipid bilayer. A related aspect is how rhodopsin is related to the lipid composition of its native membranes. For instance, in the case of a rodlike or cylinder-shaped protein inclusion in a bilayer membrane, the cross-sectional area profile as a function of depth (z) within the membrane would be uniform over the length of the protein, and zero elsewhere. A second case would be a globular or ellipsoidal protein having a parabolic cross-sectional area distribution over the length of the protein, with a maximum at the center of the protein. The third case is that of a vase-shaped protein inclusion, with larger cross-sectional area at the height of the lipid headgroups than in the center of the bilayer. Below we consider how application of the GMS approach can be used to quantify the various interfaces and describe the bilayer deformation due to rhodopsin.

View larger version (23K): [in this window] [in a new window]   FIGURE 3  Lateral deformation of the bilayer membrane due to the protein inclusion. (a) The GMS method was applied to determine the partial contributions of the different components of the membrane system to the cross-sectional area as a function of position along the bilayer normal. The center of mass of rhodopsin defines the origin of the coordinate system; positive z values are on the extracellular side and negative z values are on the cytoplasmic side. The grayscale of the transverse planes indicates the different components (from light to dark: anions, cations, water, glycerophosphocholine, acyl ester, acyl chain hydrocarbon, and protein atoms) of the Voronoi regions. Rhodopsin is depicted in terms of its backbone secondary structure elements, and the internal water molecules are shown as van der Waals sphere (CPK) models. The central slice shows only three components; viz. atoms from the protein, acyl chain hydrocarbon, and a single internal water. The slices through the headgroup regions illustrate the heterogenity in this region with all components present. (b) Quantitative results from the GMS analysis. The cross-sectional area as a function of slice position along the z axis in the unit cell is plotted for the four principal components: aqueous solvent (blue; AW(z), all water and salt atoms), hydrophobic part of the phospholipids (pale yellow; AH(z), all aliphatic acyl chain atoms), polar part of the phospholipids (cyan; AP(z), all lipid atoms except the aliphatic acyl chain atoms), and protein (black; AR(z), including posttranslational modifications). Note that the bilayer is shifted toward the extracellular domain of rhodopsin (cyan and yellow curve), where the protein has a vase-like shape, with larger cross-sectional area values close to the lipid/water interfaces than in the bilayer center.

An important aspect is that the area distributions in Fig. 3 allow one to calculate the position of the bilayer relative to the rhodopsin molecule. In Fig. 3 b, integration of the sum of the area distributions of the lipid components, A_P(z) and A_H(z), gives the volume of the bilayer, which is centered atz = +1.7 Å relative to the center of mass of rhodopsin, taking the different numbers of lipids for each monolayer into account. The following calculations are based on integrations for z values <+1.7 Å for the cytosolic monolayer, and for the extracellular monolayer the integration is for largerz values. The cross-sectional area distributions of the acyl chains A_H(z) (pale yellow), Fig. 3 b, are centered at +8.9 Å and -5.7 Å for the extracellular and cytosolic monolayers, respectively. The center of the distribution for the lipid polar headgroups A_P(z) (cyan), Fig. 3 b, in the extracellular monolayer is at a z position of +20.4 Å, and the corresponding value for the cytosolic leaflet is at -16.9 Å. The bilayer thickness from the headgroup distributions is 37.3 Å.

What is the average cross-sectional area of rhodopsin for the different regions of the bilayer, viz. corresponding to the polar headgroups and the acyl chains? The weighted average of the cross-sectional area of rhodopsin A_R(z), using the partial areas of the lipid polar headgroups A_P(z) as weights, gives the cross-sectional protein area in the vicinity of membrane surface. Its value for the extracellular monolayer is 870 Ų, and it is 960 Ų for the cytoplasmic monolayer. The corresponding area for the acyl chain region of the extracellular monolayer is 905 Ų, whereas it is 858 Ų for the cytoplasmic monolayer. These weighted average cross-sectional protein area values are subtracted from the unit cell area to obtain the monolayer area values for the different regions of the bilayer. In the headgroup region, the calculations yield a value of 63.3 Ų for the area per lipid in the extracellular leaflet, and 60.3 Ų in the cytoplasmic leaflet, indicating that the cytoplasmic monolayer is more condensed. This is perhaps not surprising, because it contains one additional lipid compared to the extracellular monolayer. By contrast, the calculation of the average area per lipid in the hydrophobic lipid acyl chain regions for each monolayer gives almost identical values of 62.3 and 62.6 Ų for the cytosolic and extracellular leaflets, respectively.

The cross-sectional area profiles for the different components, viz. protein, lipid acyl chains and headgroups, and aqueous solvent, allow one to express the bilayer deformation in terms of the monolayer curvature. We can then gain insight into the relation of the rhodopsin structure to the composition of the native lipid environment (Brown, 1997; Botelho et al., 2002). Briefly, the curvature is related to the different cross-sectional areas in the headgroup and acyl chain regions of a monolayer (Thurmond et al., 1993; Rand and Parsegian, 1997). A concave monolayer, when viewed on the headgroups, corresponds to a negative curvature, and a flat monolayer to zero curvature, i.e., an infinite radius of curvature. How can one define the monolayer curvature for a planar bilayer membrane? This apparent paradox is resolved when we consider the protein interface with water to be in effect positively curved, as in case of a vase-shaped bilayer inclusion, and the lipid component to be negatively curved. Hence a mean curvature of zero is obtained for both components together. The monolayer curvature can be estimated from the ratio of the average area per lipid in the headgroup region and in the acyl chain region, together with the average positions of the polar and hydrophobic parts of the monolayers. The results are expressed as the radius of curvature of the so-called pivotal surface, which in our calculations divides the positions of the centers of the headgroup and acyl chain distributions in a ratio of 1/3 to 2/3, that is, the pivotal surface is closer to the headgroups than the acyl chains. In this way a radius of curvature of -692 Å is obtained for the cytoplasmic monolayer, and +2144 Å for the extracellular monolayer. This is in agreement with the lipid composition of the native rod disk membranes, in which a preferential distribution of nonlamellar phase-promoting lipids is found within the cytoplasmic monolayer, with a tendency to form negatively curved, reverse hexagonal (H_II) phases (Deese et al., 1981; Miljanich et al., 1981).

Membrane protein generalized molecular surface compared with hydropathy analysis
Let us next consider the relative exposure of the membrane-embedded rhodopsin molecule to the hydrophobic and polar regions of its environment. The total contact or exposed areas with water or lipid were calculated from the GMS for each residue as described in the Appendix. The GMS area values of the ith residue of rhodopsin at the lipid/protein hydrophobic interface a_H(i), the lipid/protein polar interface a_P(i), and the protein/water interface a_W(i) are shown as a function of sequence position in Fig. 4, a–c, respectively. These interface area functions represent the accessibility of the explicit environment molecules to rhodopsin as obtained from the all-atom MD simulations. Note especially the periodicity of the lipid/protein hydrophobic interface area function a_H(i) in the transmembrane helical regions of rhodopsin, Fig. 4 a, which is characteristic for transmembrane helix bundles. The helical domains of the rhodopsin MD structure are indicated by a gray background. The periodicity is also found to extend to the headgroup region of the lipid/protein polar interface a_P(i), Fig. 4 b. However it is not as pronounced in the aqueous exposure profile of the protein/water interface a_W(i), Fig. 4 c, which is mainly outside of the helical secondary structure in the loop and turn regions.

View larger version (61K): [in this window] [in a new window]   FIGURE 4  Hydropathy analysis compared with GMS analysis of the MD model of rhodopsin in the POPC bilayer. (a) Contribution aH(i) to area of lipid/protein hydrophobic interface versus residue position (i). The exposure to the hydrophobic bilayer interior is calculated as the GMS of rhodopsin with the hydrocarbon acyl chain atoms of POPC for each amino acid residue. Posttranslational modifications are not included, for example, the palmitoylations of Cys322 and Cys323 do not appear as peaks. (b) Residue contributions aP(i) to area of lipid/protein polar interface, computed from the GMS of rhodopsin with all lipid atoms excluding the hydrocarbon part of the acyl chains. (c) Contribution aW(i) of individual amino acid residues to area of protein/water interface, calculated as the aqueous solvent-exposed GMS of rhodopsin with internal and external water molecules, as well as with the salt ions that give a very low contribution. (d) Gibbs free energy of partitioning 19-residue-long sequence windows into the membrane, obtained from the ratio of the lipid-exposed GMS and the water-exposed GMS of rhodopsin (solid line). The hydropathy index from the GMS was calculated from the protein/water interface, and the sum of the lipid/protein hydrophobic and polar interfaces (cf. text). The results are compared with a hydropathy analysis (White and Wimley, 1999) of the protein amino acid sequence (dotted line). Note the striking similarities of the two plots, in which the null crossover points of both functions are very close to one another, and the locations of the local minima and maxima are similar. (e) The hydrophobic and polar interfaces mapped onto the residue sequence, where the sequence topology scheme of rhodopsin illustrates the protein residues associated with the different GMS interfaces (cf. text). The extracellular e-side and cytoplasmic c-side of the bilayer are depicted, where the light and dark gray regions denote the headgroup and acyl chain regions, respectively. Helices are indicated as cylinders, where the first and last residues are labeled by the residue index. The lipid interface type is denoted by the filling color of the residue symbols, whereas the symbol shape marks the water interfaces. Residues exposed to the bilayer hydrocarbon interior (the lipid/protein hydrophobic interface) are shown in yellow; exposure to the lipid headgroups (the lipid/protein polar interface) is indicated by blue; and any contact with the lipid/protein hydrophobic and polar interfaces, is depicted by green. Protein residues not exposed to water are shown as circles, whereas residues exposed to internal or external water are displayed as squares. The remaining symbols indicate posttranslational modifications, and the charge states of the residues in the MD model are designated by red or blue symbol frames. An important detail is the exposure of the C2 loop residues to both the lipid/protein polar and hydrophobic interfaces, which is related to the model of lipid effects on rhodopsin activation.

(1)

The calculations were performed for the same ionization of the residues in rhodopsin as in the MD model Rh_C (Table 1). Note that due to the smoothing, the functional form does not provide a value for the hydropathy index function for the first and last nine residues of a protein. Peaks in the hydropathy plot of greater than a predefined cutoff level, zero for these whole-residue hydropathy analysis methods, are indicative of potential transmembrane helices (White and Wimley, 1999).

As distinct from the sequence-based hydropathy analysis, the GMS can be used to define a hydropathy index corresponding to the actual interfacial areas of the rhodopsin-containing MD model. The free energy of transfer from membrane to water of a 19-residue sequence window, i.e., the hydropathy index, is related to the probability to find these residues exposed to a membrane environment as compared to the probability for water exposure. Assuming the area of the molecular surface is independent of the type of environment, viz. membrane or water, then the probabilities for membrane and water exposure are proportional to the exposed molecular surface areas. The GMS-derived area exposure functions for the lipid/protein interface, a_H(i) and a_P(i), and for the protein/water interface, a_W(i), are summed independently over the window from residue to residue , as in the sequence-based analysis. The ratio of these 19-residue area exposures is interpreted as a mol fraction partition coefficient between the phospholipid and water compartments, which is then related to the standard free energy of transfer from the phospholipid bilayer to water:

(2)

It should be noted that the identification of the GMS interface with partition coefficients representing free energies of transfer is conceptually different from the relation of molecular surface exposure to the hydrophobic effect. In the latter case, the expulsion of nonpolar compounds from water is considered in terms of the water-to-hydrocarbon transfer free energy, which is related to the interfacial surface tension and the accessible surface area (Dill, 1997; White and Wimley, 1999). On the other hand the GMS-based hydropathy analysis describes partitioning of the molecular components resulting from the MD simulations of rhodopsin. When one compares the partitioning free-energy profile computed in this way (solid line) with the sequence-based hydropathy analysis (dotted line), Fig. 4 d, one can see clear similarities between both. The results clearly reflect the transmembrane helical and loop topology of rhodopsin. A positive value of the hydrophobic moment function corresponds to the transmembrane region, whereas a negative value relates to the aqueous region. Note that the zero crossover points of the hydrophobic moment functions, i.e., the MD-calculated hydropathy profile and the sequence-based hydropathy profile, are very similar, confirming that the actual position of rhodopsin in the biomembrane model is defined by the amino acid sequence. This conclusion is significant to homology-based modeling of GPCRs and bioinformatics in general. The lipid/protein and protein/water interfacial contacts as a function of the rhodopsin sequence are schematically summarized in the form of a two-dimensional topology model in Fig. 4 e. It can be seen that many residues of the transmembrane domain of helices H3 and H7 (indicated in white) are not exposed to the lipid bilayer, resulting in small hydrophobic moments relative to the other transmembrane helices. Moreover, helices H4 and H6 have the shortest transmembrane-spanning sequences corresponding to local thickness variations of the bilayer, a possible origin of rhodopsin dimerization (Fotiadis et al., 2003). Note the association of the cytoplasmic loops C2 and C3 with the lipid/water interface, which is significant for conformational changes due to the membrane environment (vide infra). The organization of the transmembrane helix bundle, which is conserved among the members of the GPCR superfamily (Baldwin et al., 1997), gives rise to a complicated pattern of protein/water and lipid/protein interfaces as described herein.

Identification of internal water molecules
An important aspect of the refined rhodopsin structure is the identification of four internal water molecules together with further internal cavities, which are not accessible to the protein surface (Teller et al., 2001). The presence of voids and internal cavities may be significant for distinguishing regions of potential mobility implicated in the functional mechanism of rhodopsin activation. How can one identify such internal water molecules? The detection of internal water molecules constitutes an important application of our GMS method. In this context, a clear definition of the term internal is necessary: an internal volume has to be surrounded by only the protein. By this definition any lipid molecule of the surrounding membrane is considered as external. Moreover, any nonprotein molecule in contact with any previously assigned external molecule can likewise be assigned as external. Here a molecule was designated in contact with another molecule if a nonzero contact area of the GMS was formed by both molecules. The consequent application of this iterative definition for external molecules results in a negative selection of the protein and any internal molecules, where the internal and external space is given by the connectivity of the network. Initially, eight molecules buried in the rhodopsin structure, four of them corresponding to the locations proposed from analysis of the crystal structure (Teller et al., 2001), and four additional water molecules occupying empty cavities in the protein, were placed in the refined crystal structure (Teller et al., 2001) utilizing a combination of the programs SOLVATE (Grubmüller, 1996) and DOWSER (Zhang and Hermans, 1996). These internal water molecules were introduced into the simulation at segment P4 together with model Rh_B, cf. Appendix. The hydrophilicity of the internal cavities of rhodopsin was then further explored using the GMS approach. One of the initially placed water molecules dissociated from its cavity, and five additional water molecules entered the protein from the bulk solvent. Thus, in the equilibrated rhodopsin MD model, a total of 12 internal water molecules could be identified by this protocol. The distribution of these internal waters within the rhodopsin molecule is visualized by their GMS with the protein in Fig. 5 a. Together with the water-exposed residues found throughout the transmembrane region (Fig. 4), a more or less continuous network of residues in contact with water is evident, even though the water molecules are not arranged in an uninterrupted chain.

View larger version (71K): [in this window] [in a new window]   FIGURE 5  The principal interfaces of a membrane protein as a result of the generalized molecular surface analysis. Each of the four views is rotated about the vertical axis by 90°. (a) The GMS interface of rhodopsin with the internal water molecules identified by the simulations, where the protein is indicated by a simplified secondary structure depiction (cf. text). (b) The GMS interface of rhodopsin with the acyl chain hydrocarbon layer (red), the acyl chain ester group atoms (gray), the glycerophosphocholine headgroup atoms (orange), and the external, or bulk, water molecules (yellow). The protein/water interface (yellow) of the extracellular and cytoplasmic domains of rhodopsin comprises 62% of the total area, which is twice as large as the lipid/protein hydrophobic interface (red) with 31% of the area. The remaining 7% corresponds to the lipid/protein polar interface (gray and orange areas). (c) The GMS area values for the lipid/protein hydrophobic interface mapped in red/white shading onto a Connolly-type molecular surface with basis color blue. For the rest of the MD simulation unit cell the generalized Voronoi regions are shown in cross sections. The Voronoi regions are gray-coded for the different components of the environment of rhodopsin (from dark to light: hydrocarbon, ester, and glycerophosphocholine headgroup atoms of the phospholipids, water, anions, and cations). (d) Simplified representation of the lipid/water interfaces calculated as described in the text (gray- and cyan-colored surfaces, the acyl chain carbonyl atoms and the headgroup phosphorus and nitrogen atoms, respectively). Rhodopsin is depicted by a simplified secondary structure representation, and its internal water molecules are shown in a CPK-type van der Waals sphere representation. The GMS method allows one to quantify the protein/water and lipid/protein interfaces at virtually any level of detail.

How is the structure of the surrounding bilayer affected by the above rhodopsin interfaces? In Fig. 5 c planar slices through the system are shown, where the gray maps depict the Voronoi regions of the components utilized for the GMS computation. It can clearly be seen that the hydrocarbon part of the lipid bilayer forms a sharp interface with the aqueous solvent and the polar lipid headgroups; it is essentially free of polar components such as water or lipid headgroups. In what follows, we refer to this as the hydrophobic lipid/water interface, which also marks the relatively sharp transmembranous lipid/protein hydrophobic interface. The latter is illustrated by the GMS mapped on the Connolly-type molecular surface (vide supra). We note that the lipid/water interfaces, as identified by the GMS method, can be compared with other definitions of the bilayer membrane geometry. For instance, the Luzzati thickness of a bilayer is the lipid head-head separation derived from electron density profiles, and the hydrophobic thickness is relevant for hydrophobic matching of proteins. Implicit in these definitions are two different lipid/water interfaces, which are approximated by the average positions of the headgroup atoms or the acyl chain carbonyl atoms, respectively (Huber et al., 2002b). We refer to these interfaces as the PN-type and COO-type lipid/water interfaces. They are calculated for each monolayer from the positions of the phosphorus and nitrogen atoms of the headgroups for the PN-type, and the carbonyl carbons of the acyl chain ester groups for the COO-type.

With the above definition, the MD model can be analyzed in terms of the PN-type and COO-type lipid/water interfaces, using the distance-weighted averages of the corresponding atom positions in each monolayer. Briefly, a surface z(x, y) as a function of xy position approximating N individual atoms, having coordinates (x_i, y_i, z_i) for the ith atom, is given by:

(3)

The nonlocal functional form, a sum of Gaussian-type radial weighting functions with width r, enables one to effectively extrapolate the lipid/water interfaces, as approximated by these surfaces, into the protein interior. Note that the surfaces are calculated under periodic boundary conditions. In Fig. 5 d, the approximation of the lipid/water PN-type interface is shown as the cyan-colored surfaces, and the lipid/water COO-type interface is depicted by the gray-shaded surfaces. These interfaces can be compared with the boundaries of the Voronoi regions from the GMS calculations, as illustrated by the slices through the bilayer in Fig. 5 c. Consequently, the PN-type and COO-type interfaces approximate the all-atom GMS in terms of a simple continuum approximation, and can be utilized for visualization of local membrane deformations.

Structural fluctuations of rhodopsin
In the previous sections, our focus has been on the influences of rhodopsin on the phospholipid bilayer. How is the protein structure affected due to the bilayer environment? Here one should remember that the x-ray structure of rhodopsin corresponds to a nonnative medium, for example its photoactivation leads to disruption of the protein crystals (Okada and Palczewski, 2001). The mobility of the equilibrated structure of rhodopsin in the membrane was studied by computing the rms fluctuations of the C backbone structure versus the average equilibrated structure for the last 5 ns of the simulation (segment P6; cf. Appendix). In Fig. 6 a, the rms fluctuation values from the MD model (shown in red) are juxtaposed with those calculated from the crystallographic Debye-Waller B-factors (black), which are related to the variances of the atom positions. The B-factors were averaged for all six available crystallographic rhodopsin structures (Palczewski et al., 2000; Teller et al., 2001; Okada et al., 2002). The rms fluctuations are measures of mobility that can be compared for the MD-generated structure and for the crystallographic structure. Note that the rms fluctuations calculated from the crystallographic B-factor for a particular atom are only valid in the absence of static disorder, so that a constant offset is zsometimes subtracted (Bernèche and Roux, 2000). However, subtraction of a constant would imply a nonuniform static disorder over the whole molecule, so that such a correction does not seem theoretically well justified. Consequently the data shown in Fig. 6 a are not offset corrected. The general trends of the rms fluctuations from the crystallographic B-factors and the MD-generated ensemble are comparable, as shown in Fig. 6 a. The rms plot from the B-factors has some gaps, as some of the rhodopsin atoms in the crystal structures are either missing, or no B-factor was assigned (Palczewski et al., 2000; Teller et al., 2001; Okada et al., 2002). It is interesting that the helical domains of rhodopsin exhibit relatively small fluctuations, whereas the extramembranous domains, especially the cytoplasmic side of rhodopsin, show higher mobility. Next we shall analyze the rms deviations of these structures relative to the crystal structure.

View larger version (50K): [in this window] [in a new window]   FIGURE 6  Quantitative analysis of the rhodopsin MD structure as compared with the crystal structure, indicating backbone fluctuations and deviations, and secondary structure stability. (a) Polypeptide backbone mobility shown by root mean square (rms) fluctuation of C carbons from the equilibrated MD ensemble for the last 5-ns segment (P6) of 15-ns simulation (red --), and calculated from the crystallographic B-factors (black line; cf. text). For clarity, the helical domain structure in parts a and b is indicated by a gray background in the plots. (b) Backbone C atom rms deviations for simulation segment P6 versus the crystal structure (blue --) in comparison with rms fluctuations from MD (red line) and B-factors (gray line). (c) Secondary structure analysis calculated with STRIDE (Frishman and Argos, 1995) for the rhodopsin polypeptide chain versus simulation time. The color code indicates the secondary structures: purple, -helical (62%); pink, 310-helical (0.11%); yellow, ß-sheet (3.8%); white, coil (16%); green, turn (18%); brown, isolated bridge (0.004%); and red,-helix (0%). Various segments of the trajectory are indicated (cf. Appendix). Note the difference of the secondary structure for model RhA in the region of cytoplasmic loop C3 compared to the later models RhB and RhC. In the second half of the C3 loop, the MD model shows fluctuations between coil and turn structures, and the rms fluctuations have the largest values after the C-terminal tail, which also oscillates between coil and turn conformations. Other regions with large rms deviations from the crystal structure appear relatively stable in terms of their secondary structure. The helical domains remain for the length of the simulations, except for ends and breaks in helices H5 and H7, which persist over the time course of the simulations.

It is noteworthy that the MD simulation of rhodopsin in a native-like membrane environment yields a new equilibrium conformation versus the crystal structure. This MD model can be considered as a refined structure for rhodopsin in a native-like bilayer environment. The packing in the crystal unit cell results in a large contact area among neighboring proteins. Moreover a number of detergent molecules (ß-nonylglucoside), are present, together with the crystallization adjuvant heptane-1,2,3-triol, and several heavy atoms such as zinc and mercury (Palczewski et al., 2000; Teller et al., 2001; Okada et al., 2002). The large rms deviations involving the loop regions might reflect the removal of these packing constraints and the heavy atoms found in the crystal together with newly formed interactions with the phospholipid bilayer. One might argue that the observed conformational changes versus the crystal structure are an artifact due to model building. Especially the conformation of the missing residues 236–240 in loop C2 and 330–334 in the C-terminal tail had to be surmised by a distance geometry approach. A particular problem was that the available crystal structures (Palczewski et al., 2000; Teller et al., 2001) came to conflicting conclusions about the direction of the turn from loop C3 into helix H6 around residue 241, hence the resulting models of residues 236–240 were on opposite sides of helix H6. It was possible to close these gaps in the polypeptide chain without sterical overlap of the modeled structure for all models. However, the MD simulations of the protein models in the bilayer environment resulted in a significant structural rearrangement of the modeled fragments and the adjacent residues on a timescale of several nanoseconds. Because the computational cost of the simulations prohibited exhaustive testing of the equilibration properties of alternative starting configurations, we could not investigate the effect of different starting structures. The MD simulations starting with model Rh_A based on the first structure (Palczewski et al., 2000) were not long enough to obtain an equilibrium conformation for this region. The subsequent models Rh_B and Rh_C were based on the second structure (Teller et al., 2001), and the MD simulations yielded an equilibrium position of this region approximately half-way between both crystal structures. The fluctuations of this turn are quite large, which might be relevant for the rotation of helix H6 upon receptor activation (Hubbell et al., 2003). On the other hand, a detailed inspection of the resulting model gave a plausible explanation for the structural changes: several residues increased dramatically their interaction with the phospholipid/water interface. For example, the initial portion of the C-terminal tail bends toward the bilayer, leading to a deeper insertion in the bilayer of the palmitoyl chain on residue Cys³²³ and to a salt-bridge formation of residue Lys³²⁵ with a lipid phosphate group. Several residues away, the side chain of Lys³³⁹ also forms another salt bridge to an adjacent phospholipid, but the residue remains close to the position in the crystal structure. The peptide loop between these two anchor points exhibits a new conformational equilibrium with small rms fluctuations but large rms deviations relative to the crystal structure (cf. Fig. 6 b). As the cysteine palmitoylations can be viewed as an anchor for the amphipathic helix H8, the presence of the bilayer also leads to appreciable deviations of this helix. The C-terminal tail beyond residue Lys³³⁹ is more or less unstructured and very flexible, in agreement with its function as a protein trafficking signal (Sung et al., 1994). The strong interaction of Lys³³⁹ with a phospholipid is remarkable given that recent solution-state ¹⁵N NMR spectroscopy data for rhodopsin show high mobility of this residue in dodecylmaltoside micelles, which lack phosphate headgroups (Klein-Seetharaman et al., 2002).

Similar rearrangements in the cytoplasmic loops C2 and C3 and in the adjacent helices result in optimized hydrophobic and electrostatic interactions with the adjacent phospholipids. In particular, Lys¹⁴¹, Phe¹⁴⁶, His¹⁵², Tyr²²³, Phe²²⁸, Lys²³¹, and Arg²⁵² exhibit movement toward the lipid/water interface, whereas Glu¹⁵⁰ and Glu²³² shift away from the bilayer toward the aqueous space. This demonstrates that the conformations of the cytoplasmic loops C2 and C3 in the crystal structure are not optimized for the presence of a bilayer. It is particularly interesting that Arg²⁵² exchanges its counterion Glu²³² in the crystal structure in favor of the negative partial charge of the phosphate group in the zwitterionic lipid POPC. By contrast, the anionic amino acyl residues exhibit no pronounced interaction with the corresponding positive partial charge of the rather hydrophobic choline group. In summary, we suggest that the observed structural changes in the functionally important cytoplasmic domain of rhodopsin are related to the transfer from the crystal packing into the phospholipid bilayer environment, and that the residue hydrophobicity (White and Wimley, 1999) provides the driving free energy.

Stability of secondary structure elements
As described above, the crystal structure and the MD-equilibrated model exhibit different conformations as a function of the sequence of rhodopsin. In addition, the MD simulations allow one to characterize the trajectories of the secondary structure fluctuations. The dependence of the secondary structure for each residue as a function of time is shown in Fig. 6 c. The relation of the simulation time to the different rhodopsin models (Rh_A, Rh_B, and Rh_C; cf. Table 1) and the different segments of the simulations (P1–P6; cf. Appendix) are indicated in the figure. In Fig. 6 c, the stability or time invariance of the -helical, the ß-strand or extended, and several of the turn and coil structures can be inferred from the noninterrupted band patterns. However, there are several regions with fluctuations between two or more secondary structures. For example, the ends of some -helices, especially H3, H4, and H7, are dynamic and depart partially from a helical backbone conformation. Furthermore, helices H5 and H7 exhibit a break around His²¹¹ and Lys²⁹⁶, respectively, which is possibly related to the nearby proline residues Pro²¹⁵ and Pro³⁰³. This can be seen from the interrupted -helical patterns at sequence positions 211 and 296. The two antiparallel ß-strands in the N-terminal tail are slightly more mobile than the second pair of ß-strands in the E2 loop. The structure of the loop C3 and the C-terminal tail is very important, and it is interesting to follow the equilibration of the models in the regions of the crystallographically unresolved residues. The initial models were generated by fitting the missing peptide fragments into the open gaps in the rhodopsin crystal structure by distance geometry, and were very similar to results calculated with the homology modeling program MODELER (Baker and Sali, 2001) (not shown). It is clear that for these regions (residues 236–239 and 328–333) a new conformational equilibrium emerged after several nanoseconds of simulation. Such large structural rearrangements do not typically result from the standard homology modeling approach for membrane proteins (Stenkamp et al., 2002). Rather, only MD simulations allow for such structural relaxations within a meaningful explicit representation of the molecular membrane environment. This observation is important with regard to GPCR structures based on homology models derived from the rhodopsin structure (Huber et al., 2003).

Equilibrium structure of rhodopsin
The rather substantial differences in the mobility of the structural domains of rhodopsin led to a further analysis of its structure. The ensemble of rhodopsin conformations from the last 5 ns of the MD model (segment P6) was aligned together with the different crystal structures that are currently available to identify the structural changes. A representative set of six conformations from the MD model and the six protein molecules from the three crystal structures, having two monomers per unit cell, are shown in Fig. 7. The spectral color of the MD-generated structures indicates the position in the sequence, starting from blue for the N-terminus to red for the C-terminus, with the crystal structures shown in gray for comparison. The 11-cis-retinal ligand can be identified in the side view, Fig. 7 a, and extracellular view, Fig. 7 b. Visual inspection of the rhodopsin structural fluctuations compared to the crystal structures reveals a high stability of the transmembrane helical fold. By contrast, significant mobility and deviations from the crystal structures are evident in the extramembranous loop portions, with larger deviations and fluctuations in the cytoplasmic regions, Fig. 7 c, versus the extracellular regions, Fig. 7 b. Appreciable differences are found in cytoplasmic loops C2 and C3, in the C-terminal tail, in the cytoplasmic ends of the transmembrane helices H4, H5, H6, and H7, and in the cytoplasmic helix H8. The light arrows indicate the cytoplasmic loops C2 and C3, as well as the C-terminal tail, which exhibit the largest deviations relative to the crystal structures.

View larger version (66K): [in this window] [in a new window]   FIGURE 7  Conformational snapshots from the MD model compared with the crystal structures. A set of five representative conformations from the MD simul