INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Bacteriorhodopsin (bR) is the light-driven proton pump protein from the purple membrane of the bacterium Halobacterium salinarium (Oesterhelt and Stoeckenius, 1971). BR contains a retinal chromophore (see Fig. 1) covalently linked to Lys²¹⁶ via a protonated Schiff base. After absorption by bR of a 568-nm photon, the retinal undergoes a conformational change that leads to the transfer of a proton from the intracellular to the extracellular side of the membrane.

View larger version (5K): [in this window] [in a new window]   FIGURE 1   Retinal bonded to Lys216. Arrows indicate the isomerized dihedral angles in the dark-adapted state.

The light-adapted form of bR contains ~100% all-trans retinal (in which all of the double bonds of the polyene are in the trans conformation). However, after several minutes in the dark, the protein reaches a dark-adapted state. This state has an absorption maximum at 558 nm and contains a mixture of two isomers of retinal (Harbison et al., 1984): all-trans and C13==C14 cis, C15==N16 syn, abbreviated here to (13,15)cis, (in which the dihedral angles C13==C14 and C15==N16 are both cis). The two corresponding forms of the protein, denoted bR₅₆₈ (all-trans retinal) and bR₅₄₈ ((13,15)cis retinal), have absorption maxima at 568 nm and 548 nm, respectively. Experiments on retinal extraction followed by high-pressure liquid chromatography have led to the suggestion that the bR₅₄₈ form populates about two-thirds of the total in the dark-adapted state (Scherrer et al., 1989; Song et al., 1995). The population ratio is modified by changes in temperature, pH, amino acid sequence of bR (Song et al., 1995), and pressure (Schulte and Bradley, 1995; Schulte et al., 1995). The retinal isomer ratio has also been determined on membranes treated with Triton X-100, a detergent that produces monomers of bR (Massote and Aghion, 1991), giving a population of 71% of bR₅₄₈ at 277 K (Scherrer et al., 1989). The relative populations of the two forms in the dark-adapted bR suggest that the difference in their free energies is k_BT (where k_B is the Boltzmann constant and T is the temperature).

Important questions remain concerning dark-adapted bR, including the isomerization pathway from (all-trans) to (13,15)cis retinal, and the factors influencing the conformational equilibrium. Calculations based on atomic models can be used to address these questions. Molecular dynamics simulations have been performed to examine other aspects of bR, including steps along the photocycle (Humphrey et al., 1994, 1998; Xu et al., 1995, 1996; Edholm et al., 1995; Ben-Nun et al., 1998; Hermone and Kuczera, 1998). A simulation model of bR₅₆₈ has been presented (Ferrand et al., 1993b) and the thermodynamic stability of internal water molecules examined (Nina et al., 1993, 1995; Roux et al., 1996). Recently a structure was proposed for (13,15)cis bR, and its photocycle was simulated (Logunov et al., 1995). The results of this study provided insight into why the photocycle of (13,15)cis does not pump protons; the Schiff base nitrogen points to the intracellular side after photoisomerization of retinal. The isomerization pathway and the dark-adaptation kinetics have also been examined (Logunov and Schulten, 1996).

In preliminary work on the dark-adapted conformational equilibrium, we used quantum chemical calculations and free energy simulations to examine the relative energies of the all-trans and (13,15)cis conformers of isolated retinal molecules, subject to harmonic restraints designed such that the conformational flexibility of the chromophore in the protein was approximately reproduced (Baudry et al., 1997). The all-trans form was found to be ~2.1 kcal/mol lower in free energy than the (13,15)cis form. This is in contrast to the higher proportion of the (13,15)cis form observed in dark-adapted bR, suggesting that the equilibrium in bR is significantly affected by interactions between the protein and the retinal and/or the effect of the retinal on the protein-protein interactions. Here the dark-adaptation pathway for the conformational change from all-trans to (13,15)cis retinal is examined. The (13,15)cis-all-trans free energy difference is calculated using umbrella sampling with the Weighted Histogram Analysis Method (Kumar et al., 1992). Agreement with the experiment is found to within ~k_BT. Calculations on modified bR allow factors strongly influencing the conformational equilibrium to be identified.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Molecular mechanics force field

The CHARMM (Brooks et al., 1983) potential energy function was employed in all of the molecular mechanics and dynamics calculations and has the following form:

(1)

The potential energy function includes bonded interactions, comprising bond stretches, bond angle bends and dihedral angle contributions, and nonbonded interactions between pairs (i, j) of atoms. In Eq. 1, b, u, , and  are the bond lengths, Urey-Bradley 1:3 distances, bond angles, and improper dihedral angles in any given configuration, and b₀, u₀, ₀, and ₀ are the reference values for these properties. The associated force constants are k, k_u, k, and k. The improper dihedral contributions are used to represent out-of-plane deformations of sp² groups. For the intrinsic dihedral angles , k is the force constant, n is the symmetry of the rotor (e.g., 3 for a methyl group), and  is the phase angle.

The nonbonded interactions are included between pairs i,j of atoms on different molecules and on the same molecules separated by three or more bonds. They consist of a Lennard-Jones term, with parameters _i,j and _i,j, and a Coulombic electrostatic term between partial charges q_i, q_j. The dielectric constant  = ₀. _r was set to  = ₀, i.e., _r = 1. Hydrogen bonds are described by the nonbonded terms in the energy function. In all of the calculations the long-range electrostatic terms were smoothly brought to zero at a cutoff of 12 Å by multiplication by a cubic switching function between 10 Å and 12 Å. Pairs of atoms on the same molecule separated by only two bonds may interact via a Urey-Bradley term harmonic in the distance between atoms i and j.

Derivation of an adiabatic potential energy map

The Automatic Map Refinement Procedure (AMRP) was used to derive adiabatic potential energy maps for rotation about the C13==C14 and C15==N16 dihedral angles of retinal. This method, described in detail by Baudry et al. (1997), uses iterative energy minimizations in all of the wells in the potential energy surface examined. Smooth interpolation between the points obtained was performed using Akima's quintic polynomials implemented in the IDL package.

Free energy calculations on retinal in bR

The goal of the free energy calculations was to compute the free energy difference between the two conformers in question. In theory this is independent of the pathway between them. However, the pathway is important in practice, as it influences the statistical accuracy of the calculations. In the present work we performed umbrella sampling along the "bicycle-pedal" pathway (Warshel, 1976), defined as ₁ = ₂, where ₁ is the C12-C13==C14-C15 dihedral and ₂ is that of C14-C15==N16-C. This pathway was chosen because it preserves the direction of the retinal chain and is accompanied by only a small change in its shape, thus minimizing the environmental perturbation at each step. It has been suggested to be the pathway taken by retinal in bR during dark adaptation (Orlandi and Schulten, 1979; Tavan et al., 1985; Seltzer, 1987; Logunov and Schulten, 1996).

The starting structure for the MD free energy calculations on these molecules was obtained from Roux et al. (1996), based on the crystallographic structure of Grigorieff et al. (1996). Water molecules buried in the protein were placed according to thermodynamic criteria (Nina et al., 1995; Roux et al. 1996). These essentially reproduce the hydration patterns suggested by the newer x-ray crystallographic analyses of Pebay-Peyroula et al. (1997) and Luecke et al. (1998). Asp⁸⁵ is involved in hydrogen bonds with two water molecules that are in very good agreement with those of the model of Luecke et al. (1998). Moreover, Pebay-Peyroula et al. suggest that four water molecules may lie between Asp⁹⁶ and the retinal, as in the model of Roux et al. (1996). However, we have repeated the basic conformational equilibrium calculation discussed in the present paper on the wild-type bR with the newer structure of Pebay-Peyroula et al. (1997). The result obtained was very similar to that obtained with the older, electron crystallographic structure.

To investigate the protein-retinal or water-retinal interactions, five systems were studied, four of which include some modification of the wild-type:

Wild-type bR, including five water molecules placed by Roux et al. (1996), i.e., four water molecules in the channel between H15 and Asp⁹⁶, and a water near the Schiff base forming a hydrogen bond between the Schiff base, Asp⁸⁵ and Asp²¹². We call this "hydrated" bR.

A bR in which these five water molecules close to the Schiff base are removed (called here the "dehydrated retinal binding site").

A bR in which Asp⁸⁵ is mutated to Ala (called "D85A").

A bR in which Asp²¹² is mutated to Ala (called "D212A").

A bR in which both Asp⁸⁵ and Asp²¹² are mutated to Ala (called here the "double mutant").

The retinal force field of Nina et al. (1995), modified as in Baudry et al. (1997), was used in the MD calculations to evaluate the ((13,15)cis-all-trans) potential energy and free energy differences. For the protein part, version 22 of the CHARMM force field was used (MacKerell et al., 1998). The umbrella sampling method of Valleau and Torrie (1977) was used to calculate the two-dimensional potential of mean force (PMF) along ₁ and ₂. In this method, N_W biased distributions are generated using harmonic functions of the form

(2)

centered on successive values of ₁ⁱ and ₂ⁱ. The umbrella sampled distributions were analyzed using the Weighted Histogram Analysis Method (WHAM) (Kumar et al., 1992), which has been extended to calculations along more than one degree of freedom (Roux, 1995). The WHAM equations express the optimal estimate for the unbiased distribution function as a weighted sum over the N_W individual biased distribution functions (Kumar et al., 1992; Roux, 1995) In the case of two variables, ₁ and ₂, if n_i is the number of independent data points used to construct each biased distribution function, we have (Roux, 1995)

(3)

In Eq. 3, (₁, ₂) is the distribution function of the dihedral angles ₁, ₂, and w_i and w_j are the harmonic restraint potentials used in umbrella sampling, where i and j denote the two constrained degrees of freedom. The undetermined constants F_j are defined from

(4)

Because the distribution function itself depends on the set of constants F_j, the WHAM equations 3 and 4 must be solved self-consistently. In practice this is achieved through an iterative procedure. Self-consistency was typically reached after 2500 iterations, corresponding to a maximum change in F_j of 2 × 10⁴ kcal/mol between successive iterations. The PMF, W() along a coordinate  = (₁, ₂), is defined from the average distribution function ():

(5)

where * and W(*) are arbitrary constants. The version of WHAM used takes advantage of the periodicity of the reaction coordinate to remove hysteresis.

The potential of mean force was calculated for ₁ (= ₂) varying from +180° to 180° in steps of 30°, using the following protocol:

Step i. ₁ and ₂ were restrained to the desired values (i.e., ₁ = 180°, ₂ = 180°; ₁ = 150°, ₂ = 150°; ... ; ₁ = 180°, ₂ = 180°) with a force constant of 15 kcal/mol/degree², giving a total of 13 windows of simulation.

Step ii. A 1-ps equilibration run was performed.

Step iii. Langevin production dynamics runs were performed, using the Langevin parameters listed in Roux et al. (1996).

For the simulations on mutated bR (i.e., mutations D85A, D212A, double mutant) and dehydrated binding-site bR, the time lengths were 20 ps per production window, giving a total production time of 260 ps for each modified protein. For the simulations on wild-type, hydrated bR, two simulations were performed:

A simulation with ₁ varying from +180° to 180°, for which the production runs were performed for 18 ps. An additional 18-ps production run was performed in a window centered on ₁ = 160°, ₂ = 160° to increase the sampling on this part of the diagonal, over which the free energy was found to have a relatively steep gradient. The total production time was 252 ps. This simulation is called the "forward" simulation.

A simulation with ₁ varying from 180° to +180°, for which the production runs were 20 ps long. The total production time was thus 260 ps. This simulation is called the "backward" simulation.

For all of the simulations the potentials of mean force were calculated with 81 bins of width 5° for values of ₁ and ₂, ranging from +202.5° to 202.5°. For the simulations on wild-type, hydrated bR, the potential of mean force was calculated from the forward and backward simulations independently, and from the combined umbrella-sampled distributions from both the forward and backward simulations, which is equivalent to a production time of 512 ps.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Wild-type bacteriorhodopsin

Isomerization pathway

View larger version (80K): [in this window] [in a new window]   FIGURE 2   Adiabatic potential energy map for rotation around C13==C14 (1) and C15==N16 (2). (a) In vacuo retinal (Baudry et al., 1997), calculated with the parameters of Table 1, set B. The dashed line indicates the bicycle-pedal diagonal. (b) Retinal in bR, calculated with the parameters of Table 1, set B.

1

2

1

2

1

2

1

2

View larger version (77K): [in this window] [in a new window]   FIGURE 3   Adiabatic potential energy map for rotation around C13==C14 (1) and C15==N16 (2) retinal in bR. (a) Calculated with the parameters of Table 1, set A. The arrow indicates the lowest energy pathway for trans to cis isomerization of the two double bonds. (b) Calculated with the intrinsic dihedral force constants of 1.88 kcal/mol for these dihedrals. Arrows indicate the lowest energy pathways for trans to cis isomerization of the two double bonds.

Decomposition of the (13,15)-all-trans potential energy difference

Potential of mean force along the bicycle-pedal pathway

View larger version (28K): [in this window] [in a new window]   FIGURE 4   (a) Free energy surface along the bicycle-pedal diagonal, calculated with the parameters of Table 1, set B. The 22 kcal/mol energy contour is represented by a dashed line. (b) Free energy profile along the bicycle-pedal diagonal for retinal in vacuo constrained to have the same (all-trans) average structure as retinal in bR (see Baudry et al., 1997). Only the portion of the two-dimensional free energy map corresponding to the strict diagonal-pedal (i.e., 1 = 2) is shown here, to facilitate the comparison with a. (c) Free energy surface along the bicycle-pedal diagonal, calculated with the parameters of Table 1, set A.

View this table: [in this window] [in a new window]   TABLE 2   Free energy (A) and potential energy (E) differences

Modified bacteriorhodopsin

Calculations were performed of the (13,15)cis-all-trans free energy difference for the modified bR structures described in Methods. The results are summarized in Table 2.

The interaction energy calculations discussed above (under Wild-Type Bacteriorhodopsin) suggest that internal water molecules might stabilize the all-trans conformer in the protein. To determine if this result is also true in free energy terms, the (13,15)cis-all-trans free energy difference for dehydrated binding-site bR, with the internal water molecules removed from retinal binding pocket, was calculated. This led to the (13,15)cis being the form of lowest free energy, as in the hydrated wild type. However, the free energy of the the all-trans, dehydrated retinal binding pocket species was found to be 4.2 kcal/mol higher than the (13,15)cis species (1.1 kcal/mol in the case of the hydrated species). Therefore these calculations confirm the stabilization of the all-trans species by the internal buried water molecules in bR. Experimental dehydration has been shown to destabilize the all-trans chromphore (Korenstein and Hess, 1977). However, the water molecules removed in the present calculations are the five water molecules placed by Roux et al. (1996) between residues Asp⁹⁶ and Asp⁸⁵. This dehydration probably does not correspond to that experimentally performed (Zaccaï, 1987; Ferrand et al., 1993a; Lehnert et al., 1998), in which the water molecules removed by dehydration are mainly localized around the lipid headgroups. The dehydration simulated here corresponds to a selective dehydration of the core of bR, in the retinal binding site.

The interaction energy calculations also suggested that the effect of Asp²¹² might be to stabilize the (13,15)cis isomer, the opposite of the effect of the water molecules. The calculation of A for the D212A mutant led to the all-trans species being 0.8 kcal/mol lower in free energy than the (13,15)cis species. This result again concurs with the interaction energy calculations, in that Asp²¹² significantly stabilizes the (13,15)cis isomer.

The interaction energies showed no significant influence of the Asp⁸⁵ residue on the cis-trans interaction energy difference. This was also found for the calculations of A for the D85A mutant, which gave the result that the (13,15)cis species is 1.6 kcal/mol lower in free energy than the all-trans species. This A is not significantly different from the wild type, again in accord with the interaction energy calculations. There is experimental evidence that the protonation of Asp⁸⁵ affects the equilibrium, although not very strongly (Fischer and Oesterhelt, 1979; Turner et al., 1993; Balashov et al., 1996). However, the accuracy of the present free energy calculations is certainly not better than k_BT, and consequently only large population changes are amenable to investigation by simulation techniques.

Finally, to investigate the effect of complete removal of the negative charge around the Schiff base, free energy calculations were performed on the double mutant D85A and D212A. The (13,15)cis-all-trans free energy difference obtained was 1.9 kcal/mol, the all-trans species being the more stable, as in the case of the single mutation D212A. This result suggests that mutation of these aspartic acids to alanine, which renders impossible the formation of hydrogen bonds between retinal and these side chains, leads to a stabilization of the all-trans species.

Table 2 also gives the (13,15)cis-all-trans adiabatic potential energy differences (E) for the above molecular systems. The values of E and A are generally quite close, with a maximum difference of ~1 kcal/mol. This suggests that the calculated free energy differences are essentially enthalpic.

Structures of bR in the wild-type, hydrated protein

View larger version (28K): [in this window] [in a new window]   FIGURE 5   (a) Average structure of the retinal region in the all-trans form. Labels A-E for the water molecules are taken from Roux et al. (1996). (b) Average structure of the retinal region in the (13-15)cis form. This figure was made with the program VMD (Humphrey et al., 1996).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The present work involved investigating factors influencing the (13,15)cis-all-trans conformational transition and thermodynamic equilibrium in dark-adapted bacteriorhodopsin. Errors in the calculated free energies are unlikely to be smaller than ~k_BT (i.e., 0.6 kcal/mol at 300 K), even for long simulations on systems of a limited number of atoms (Baudry et al., 1997). Comparison of the wild-type hydrated bR free energy difference calculated from the forward, backward, and forward + backward simulations shows results that are in agreement with each other to within ~k_BT. Moreover, our best estimate of the wild-type A of 1.1 kcal/mol is also within ~k_BT of experiment. Consequently, we consider the present results to be satisfactory, given the inherent limitations of the methods.

The calculations identify two groups of atoms that strongly influence the conformational equilibrium. The double mutation D85A/D212A modifies the free energy differ-ence by ~3 kcal/mol, in favor of the all-trans species. In contrast, removal of the internal water molecules stabilizes the (13,15)cis species by ~3 kcal/mol. Thus the Asp residues located near the Schiff base and the internal water molecules have effects on the stability of the two species that are opposite but of similar magnitude. The present calculations suggest that the opposing actions of these two groups of atoms approximately cancel in wild-type bR, leading to a similar free energy of the two species to within ~k_BT. The effects on the free energy difference of the protein modifications examined here are ~5k_BT and are likely to be above the error level. The results suggest that if it were possible experimentally to remove the internal water molecules, or to perform the double Asp  Ala mutation, the cis/trans equilibrium should be measurably affected. To our knowledge, these experiments have not been yet been performed.

The bR model used in this paper is a monomeric model of the protein, in which the average effect of the environment is simulated approximately by the use of weak harmonic restraints on the position of -carbons located at the surface of bR (Ferrand et al., 1993b; Roux et al., 1996). Repetition of the present calculations using a model including explicitly trimeric bR in a lipid bilayer would be of interest, although computationally expensive. However, it has been found that the 1:2 population ratio between the two conformers of retinal in dark-adapted bR remains on dissociation into monomeric bR (Scherrer et al., 1989; Song et al., 1995).

The pathway for conformational change from (13,15)cis to all-trans in dark-adapted bR is currently unknown. Our calculations separate the factors influencing the pathway into two: the intrinsic torsional term and the rest (the "environment"). Whereas the intrinsic terms favor sequential rotation, the environment favors a bicycle-pedal mechanism. We find that the lowest-energy pathway found depends critically on the balance between these effects. As quantum-chemical and molecular-mechanical calculations become more precise, it should be possible in the near future to use theory to help decide which of the two pathways is favored in bR.