Computational Analysis of the Proton Translocation from Asp96 to

Oct 24, 2006 - Yoshiharu Sato,† Masayuki Hata, Saburo Neya, and Tyuji Hoshino*,‡. Graduate School of Pharmaceutical Sciences, Chiba UniVersity, Ch...
0 downloads 0 Views 672KB Size
Download PDF
22804

J. Phys. Chem. B 2006, 110, 22804-22812

Computational Analysis of the Proton Translocation from Asp96 to Schiff Base in Bacteriorhodopsin Yoshiharu Sato,† Masayuki Hata, Saburo Neya, and Tyuji Hoshino*,‡ Graduate School of Pharmaceutical Sciences, Chiba UniVersity, Chiba 263-8522, Japan ReceiVed: May 25, 2006; In Final Form: August 16, 2006

The potential energy change during the M f N process in bacteriorhodopsin has been evaluated by ab initio quantum chemical and advanced quantum chemical calculations following molecular dynamics (MD) simulations. Many previous experimental studies have suggested that the proton transfer from Asp96 to the Schiff base occurs under the following two conditions: (1) the hydrogen bond between Thr46 and Asp96 breaks and Thr46 is detached from Asp96 and (2) a stable chain of four water molecules spans an area from Asp96 f Schiff base. In this work, we successfully reproduced the proton-transfer process occurring under these two conditions by molecular dynamics and quantum chemical calculations. The quantum chemical computation revealed that the proton transfer from Asp96 to Shiff base occurs in two-step reactions via an intermediate in which an H3O+ appears around Ala215. The activation energy for the proton transfer in the first reaction was calculated to be 9.7 kcal/mol, which enables fast and efficient proton pump action. Further QM/MM (quantum mechanical/molecular mechanical) and FMO (fragment molecular orbital) calculations revealed that the potential energy change during the proton transfer is tightly regulated by the composition and the geometry of the surrounding amino acid residues of bacteriorhodopsin. Here, we report in detail the Asp96 f Schiff base proton translocation mechanism of bacteriorhodopsin. Additionally, we discuss the effectiveness of combining quantum chemical calculations with truncated cluster models followed by advanced quantum chemical calculations applied to a whole protein to elucidate its reaction mechanism.

1. Introduction Bacteriorhodopsin, a 27 kDa membrane protein, is a lightdriven proton pump found in a membrane of the extreme halophile Halobacterium salinarum.1 Bacteriorhodopsin consists of 248 amino residues and has very stable2 seven transmembrane helices to make a pocket for retinal chromophore. The retinal is covalently bound to Lys216 through a protonated Shiff base. Light-induced isomerization of the chromophore from the alltrans to the 13-cis form triggers a cyclic reaction that comprises a series of intermediates named I, J, K, L, M, N, and O. These intermediates can be identified by the difference in infrared (IR) spectra, particularly labeled by BR568,I480,J625,K610,L550,M412,N550,O640.3-5 The proton pump of bacteriorhodopsin is achieved by the following five reactions: (1) proton transfer from the Schiff base to Asp85, (2) proton release from Glu204 or Glu194 to the extracellular side during the L f M conversion, (3) proton transfer from Asp96 to the Schiff base during the M f N conversion, (4) proton uptake from the cytoplasmic side during the N f O conversion, and (5) proton translocation from Asp85 to Glu204 or Glu194 in the O f BR conversion.6 Water molecules are assumed to play an essential role in the unidirectional proton transfer and the switching mechanism in bacteriorhodopsin. The role of water molecules in the extracellular region at the early reaction stages has been studied in detail by several experimental6,7 and theoretical approaches.8 * Corresponding author. Phone: +81-43-290-2926. Fax: +81-43-2902925. E-mail: [email protected]. † Present address: Japan Biological Information Research Center (JBIRC) AIST Waterfront Bio-IT Research Bldg. 7F 2-42 Aomi, Koto-ku Tokyo 135-0064. ‡ PREST, Japan Science and Technology Agency, 4-1-8 Honcho Kawaguchi, Saitama, Japan.

However, the presence and the role of water molecules in the cytoplasmic half side are still unclear. It has been commonly accepted that water molecules are required for proton transfer from Asp96 to the Schiff base;9,10 however, there remain two possibilities for the mechanism and pathway of the proton transfer. One is that the reprotonation of the Schiff base from Asp96 is due to the Brownian motions of water molecules and amino acid side chains. The other is that a single linkage water chain between Asp96 and the Schiff base is formed, leading to the proton translocation. Our previous studies using ab initio quantum chemical calculation suggested that the proton transport in the L-M-N process was achieved by the presence of a water chain,11-13 although such a single-linkage water chain had not been found by any crystallographic or NMR structure analysis. Recently, the validity of this proposal has been supported by the N′ state X-ray crystallographic structure14 and by the results of a hybrid QM/MD simulation study.15 In the present study, we focused on the unidirectional proton translocation mechanism in the M f N conversion, and we revealed the precise mechanism of proton transfer in the cytoplasmic half side using several computational methods. First, we performed MD simulation in the M state of bacteriorhodopsin in order to confirm the stability of the single water chain between Asp96 and the Schiff base. MD simulation has been widely used for elucidating the biophysical phenomena such as signal transduction16 and molecular recognition.17 The dynamics and function of water molecules in a protein have been successfully elucidated by MD simulation,18,19 while the dynamics and geometries of water molecules are usually difficult to reveal by experimental methods. This is a fine advantage of MD simulation. However, MD simulation does not allow covalent

10.1021/jp0632081 CCC: $33.50 © 2006 American Chemical Society Published on Web 10/24/2006

Proton Translocation in Bacteriorhodopsin

J. Phys. Chem. B, Vol. 110, No. 45, 2006 22805

Figure 1. MD simulation model for bacteriorhodopsin. (a) Crystal structure of bacteriorhodopsin. Bacteriorhodopsin is solvated spherically in (b) and is solvated with water molecules in a rectangular box in (c). The lipid bilayer-water-protein model system is shown in (d), and a snapshot of the lipid bilayer and protein viewed from the cytoplasmic side is shown in (e).

bonds to break or form and therefore cannot be applied to the analysis of a chemical reaction process. For observation of proton transfer, we performed ab initio quantum chemical calculation using a truncated model system constructed from the structure obtained by MD simulation, and we elucidated the pathway of proton transfer and its potential energy change. The results obtained from quantum chemical calculation using the truncated model system contain some degree of inaccuracy in geometry and energy. Hence, we performed additional energy calculation with IMOMM20 and FMO methods21,22 in order to refine the estimation of energy and geometry of the proposed intermediates and to assess the results of the truncated quantum chemical calculation. 2. Materials and Methods 2.1. Molecular Dynamics Simulations with the ProteinWater Solvent System. For the late-M state MD simulation, the initial structure was constructed from an X-ray crystal structure (PDB code 1CWQ23) (Figure 1a). We performed two kinds of MD simulations for the M state bacteriorhodopsin: a no-cutoff (NCO) method and a particle mesh Ewald (PME) method.24 These two methods differ in the treatment of longrange electrostatic and van der Waals interactions. The former treats the long-range interactions without approximation, but it is only applicable to isolated systems. Hence, the NCO model was spherically solvated by about 10 000 water molecules (Figure 1b). On the other hand, the PME method, which has been used in most of the recent MD studies, treats long-range interaction by introducing a cutoff radius to calculate the real space term and the reciprocal contribution separately. The PME model was solvated with a rectangular box of 50.7 Å × 64.7 Å × 81.8 Å in size (Figure 1c). The different treatment of surrounding electrostatic effects may cause some difference in the estimation of the dynamics and the fluctuations of a protein

and water molecules. By executing two different models’ simulations, we closely examined the dynamics of the hydrogen bond network (HBN) and the water chain. We also examined the validity of the existence and the stability of the water chain by comparing the results. 2.2. Molecular Dynamics Simulation with the ProteinLipid Bilayer-Water Solvent System. The model of MD simulation in the lipid bilayer condition for the M state bacteriorhodopsin was constructed as shown in Figure 1, parts d and e. In our previous study, we successfully performed MD simulations of sensory rhodopsin in the lipid bilayer condition.16 We used the same procedure for the construction of the model system. For the component of lipid bilayer, we selected phosphatidylglycerophosphate monomethyl ester (PGP-Me), that is a main component of purple membrane (PM).25 Bacteriorhodopsin was surrounded by 56 lipid molecules, and 10 442 water molecules were generated at the upper and lower sides of the bacteriorhodopsin-lipid bilayer complex. 2.3. Computational Details of the Molecular Dynamics Simulation. MD simulations and analyses were carried out using the AMBER7 program package26 together with the Amber ff99 force field27 except for retinal. Atomic charges and stable conformation of the retinal Schiff base were deduced from ab initio quantum chemical calculations, in which density functional theory (DFT)28,29 was applied to a whole retinal Schiff base structure using GAUSSIAN98.30 The hybrid B3LYP method and 6-31G** basis set were used for the DFT calculations, and charges were computed by the RESP method31 based on the DFT results. The force field was assigned for the retinal by using the antechamber module in the AMBER7 program package, except for the force parameter for C15-N torsion. The torsional potential around this C15-N bond was set to 30.0.32 An integration time step of the simulation was set to 1 fs. In the PME calculation, a cutoff distance of 10 Å was set

22806 J. Phys. Chem. B, Vol. 110, No. 45, 2006 for nonbonded interactions and a dielectric constant of  ) 1.0 was employed. A periodic boundary condition was applied, and the pressure was kept constant for the whole system. In the NCO calculation, the harmonic boundary potential was applied to maintain the water sphere and keep the volume of the system constant. Hence, the simulation was carried out under the NVT ensemble condition. The force constant of the boundary potential was 1.5 kcal/mol/Å. Coulomb interaction and van der Waals term were explicitly calculated without any approximation by using MDGRAPE-2.33,34 The protein-lipid bilayer-water solvent system calculation was done by the periodic boundary condition. A cutoff distance of 12 Å was used, and a harmonic boundary constraint was applied to protein residues for keeping their appropriate positions during the first 100 ps of MD simulation to prevent the structural collision due to the manual docking of the protein and lipid bilayer and to flex the interactions between them. Other conditions were same as those used in the PME calculation with the water solvent system. 2.4. Ab Initio Quantum Chemical Calculation with the Truncated Model System. Theoretical approaches using quantum chemical calculation provide us a variety of information such as structures of intermediates, potential energy surfaces, and a distribution of electron density. For example, we can reveal the activation energy and the direction of the reaction by estimating the potential energy surface. It is especially hard to experimentally reveal the transition state structures or the intermediates in a rapid reaction. Quantum chemical calculations can predict these structures and clarify the reaction mechanism. This is a fine advantage of this theoretical approach. The protontransfer mechanisms of bacteriorhodopsin have been investigated in detail using this approach. 2.4.1. Construction of the Model Reaction System. The cluster model for the quantum chemical calculation was constructed on the basis of the structure obtained by the NCO model MD simulation. Since it was revealed by our MD simulation that Asp96, retinal, Lys216, and the main chain atoms of Thr46 and Ala215 are essentially involved in the formation of a water chain, these residues and five water molecules were extracted from the final structure of the 1 ns NCO MD simulation. To construct the computational model, the carbon and nitrogen atoms located at the truncating boundary were replaced with H atoms. The total number of atoms was 77. The positions of the carbon and hydrogen atoms terminating covalent bonds at the boundary were fixed during the geometry optimization based on the assumption that these atoms do not move flexibly during the proton-transfer reaction. 2.4.2. Calculation Method. Computations for geometry optimization were performed by the restricted Hartree-Fock (HF) method with the 3-21G** basis set. The 6-31G(d,p) basis set was used for the refinement of the energy and the geometry estimated on each stationary point (late-M, MH, N). To reproduce the environment inside the protein by taking the influence of surroundings into account, the self-consistent reaction field calculation was carried out by the Onsager method.35 The radius for the cavity was set at 6.27 Å, and the dielectric constant  was set at 20.0. The computational program used was Gaussian98.30 2.5. Advanced Quantum Chemical Calculation. We performed additional QM/MM calculations by the IMOMM method for the purpose of reevaluation of the surrounding effects and the potential energy change of the reaction. The protein structure was obtained by superimposing the respective optimized structure that was obtained by the above quantum chemical calculation onto the final structure of NCO MD simulation that

Sato et al. had been used as the template for the truncated model of the quantum chemical calculation. The active-site atoms included in the above truncated quantum chemical calculations were also treated quantum mechanically with the GAMESS36 code, while the remaining part of the protein was treated molecular mechanically using the TINKER37 code with the AMBER ff99 force field. The 3-21G** basis set was used for energy and geometry estimations of the QM region. Fragment molecular orbital (FMO) calculations were performed using the Abinit-MP package.21,22 Single-point energy calculations in the HF/STO-3G level were performed for the late-M and N state structures obtained by the IMOMM calculation. The unit of fragment was set to be one residue, except for retinal and Lys216. We treated the complex of Lys216 with retinal as one fragment. Figures of the structure in this paper were drawn by molecular visualizing software RasMol,38 VMD,39 and MacMolPlt.40 3. Results 3.1. Molecular Dynamics Simulation of the Late-M State Bacteriorhodopsin. Molecular dynamics simulation for the M state bacteriorhodopsin was performed in order to confirm whether the water chain from Asp96 to Schiff base is stably present in the late-M state. The presence of this water chain in the N′ state (the N state can be divided into early and late stages,41 and the late stage of the N state is called the N′ state) was experimentally detected by a crystallographic study.14 If we assume that the proton transfer from Asp96 to Schiff base occurs directly, a single water chain must exist in the hydrophobic channel. However, such a water chain has not been detected in the M state. We performed MD simulation on the M state by using different two methods. The M state crystallographic structure determined by Sass et al.23 is used as an initial structure of bacteriorhodopsin in both models. Two water molecules were added to the calculation model at the inside of the cytoplasmic hydrophobic channel between Asp96 and Schiff base so that four water molecules formed a single linkage of water molecules like a chain. In both NCO and PME MD calculations, the water chain from Asp96 to Schiff base is formed with four water molecules at the end of the 1 ns MD simulation (Figure 2). However, the dynamics of water molecules and the change in the hydrogen bond network (HBN) pattern during the simulation are different in the two methods. It is important for a direct proton translocation that a single linkage of water molecules is formed in such a manner that Asp96 is a proton donor and the Schiff base is an acceptor. In the NCO calculation, the hydroxyl hydrogen atom of Asp96 has a direct hydrogen bond with Thr46 (Figure 2a) and cannot work as a donor for the Schiff base. Accordingly, a rearrangement of HBN is required to initiate the proton translocation. To investigate the energy change with the rearrangement of HBN and the reaction process of proton transfer, we used the snapshot structure of the NCO simulation at 1 ns as a template of the truncated model for the following ab initio quantum chemical calculation. 3.2. Ab Initio Quantum Chemical Calculation with the Truncated Model. 3.2.1. Potential Energy Change with Recombination of HBN. Geometry optimization was performed by using the cluster model that was constructed from the snapshot structure at the last point of the 1 ns NCO simulation. The optimized structure was obtained as shown in Figure 3 and is denoted by “S1”. We prepared four additional truncated model structures (S2-S5) in which HBN patterns were different and changed in sequence from S1 to S5. Geometry optimization was

Proton Translocation in Bacteriorhodopsin

J. Phys. Chem. B, Vol. 110, No. 45, 2006 22807

Figure 2. Snapshot structures at 0 ps (left column) and average structures for the last 3 ps in the 1 ns MD simulation. Dotted lines indicate the hydrogen bond interactions. (a) NCO MD simulation. Four water molecules form a water chain connecting Asp96 and the Schiff base. The pattern of the hydrogen bond network (HBN) at the end of the 1 ns simulation is different from the snapshot structure at 0 ps. (b) PME MD simulation. Four water molecules form a water chain, and the alignment of water molecules in the snapshot structure at 0 ps is still kept in the average structure for the last 3 ps.

performed for each structure, and each potential energy was estimated. The transition state structures and their energies were also calculated by the QSTN method implemented in the Gaussian98 package. Asp96 works as a donor, and the Shiff base works as its acceptor in structures S3, S4, and S5. The patterns of hydrogen bonds from the hydroxyl side chain of Thr46 are different in these three structures. From the viewpoint of potential energy, S5 is the most favorable among the five different structures in HBN patterns, and the S5 structure indeed appears in the average structure during the last 3 ps of the PME calculation (Figure 2b). The activation energy level of each rearrangement reaction is considerably low, indicating that these rearrangements of the HBN easily occur in physiological conditions. From these results, we regarded the S5 structure as being ready for the proton translocation from Asp96 to the Schiff base, and we employed it as the starting structure in the calculation of potential energy change with proton translocation.

3.2.2. Proton Translocation from Asp96 to Schiff Base. The starting structure obtained by the previous calculation is shown in Figure 4 and is denoted by “a”. This represents the theoretically determined structure for the late-M intermediate. We obtained the potential energy profiles along the proton translocation by the quantum chemical calculation. By closely monitoring the potential energy curve with proton transfer from Asp96 to the Schiff base, we confirmed the existence of an intermediate in which an H3O+ molecule appears around Ala215 (Figure 4b); that is, the proton translocation from Asp96 to the Schiff base occurs in two-step reactions. We call this intermediate the MH intermediate here. The carboxyl oxygen atom of Ala215 stabilizes the H3O+ molecule, and Ala215 also has an essential role in stabilization of the water chain from Asp96 to the Schiff base as observed in our MD simulation on the late-M state. Comparison of the potential energies of the late-M and MH intermediates indicated that the MH intermediate was more stable than the late-M intermediate by 6.7 kcal/mol. Computation

22808 J. Phys. Chem. B, Vol. 110, No. 45, 2006

Sato et al.

Figure 3. Structural and potential energy changes following the recombination of the HBN. Structures S1-S5 are different in the pattern of the hydrogen bonds. The potential energy of the respective stable structure and the transition state is presented, relative to structure S1. Hydrogen bonds that are different from the previous structure are shown in red.

Figure 4. Structural changes following the proton transfer from Asp96 to the Schiff base. Numerals are interatomic distances in units of angstrom. Hydrogen bonds that are different in the late-M state and the N state are emphasized with thick dotted lines.

was performed to search for a saddle point on the potential energy surface. As a result, the optimized structure (see the Supporting Information for the optimized structure), which is the transition state of the late-M f MH conversion reaction (TS1), was obtained. A vibrational analysis at the transition state indicated the appearance of one vibration mode that has an imaginary frequency. The steepest decent reaction paths were calculated for both the forward and reverse directions following the vibration mode corresponding to the imaginary frequency. This routine eventually provided the lowest potential energy curve along the intrinsic reaction coordinate (IRC) that connected the reactant and the product (Figure 5). It is suggested

Figure 5. Energy diagram for the late-M f N proton transport reaction in bacteriorhodopsin calculated by the truncated model quantum chemical calculation; (a-c) correspond to the structures shown in Figure 4. The late-M f MH conversion appears to be the rate-limiting step of the reaction. The abscissa is the intrinsic reaction coordinate (IRC) distance measured from (a).

from Figure 5 that the activation energy for the late-M f MH conversion is 9.7 kcal/mol. The structure for the N intermediate was determined through the geometry optimization, shifting the location of protons from the MH intermediate. The optimized N state structure is shown in Figure 4c. An energy comparison indicates that the N intermediate is more stable than the MH intermediate by 14.8 kcal/mol (Figure 5). The structure for the transition state was also obtained by the geometry optimization (see the Supporting

Proton Translocation in Bacteriorhodopsin

J. Phys. Chem. B, Vol. 110, No. 45, 2006 22809

TABLE 1: Comparison of the Potential Energy Change among Different Dielectric Constants and Different Modelsa method

late-M state

MH state

N state

QM(HF3-21G**,  ) 20.0) QM(HF6-31G(d,p),  ) 1.0) QM(HF6-31G(d,p),  ) 4.0) QM(HF6-31G(d,p),  ) 20.0) IMOMM (HF3-21G**,  ) 1.0) IMOMM (HF3-21G**,  ) 4.0) IMOMM (HF3-21G**,  ) 20.0)

0 0 0 0 0 0 0

-6.7 0b 19.1 8.3 -4.4 -11.9 -20.9

-21.5 43.2 2.5 -21.5 23.2 -17.9 -37.0

a The energies are presented relative to the late-M state in units of kcal/mol. b Atom geometry was converged to that of the late-M intermediate when optimization was performed starting from this structure.

Information for the optimized structure). A vibrational analysis at the transition state showed a vibration mode corresponding to an imaginary frequency. Following the forward and reverse directions of the vibration mode, the lowest energy curve connecting the MH state and the N state is shown in Figure 5. This curve indicates that the activation energy for MH f N conversion is 3.4 kcal/mol. IRC computations have successfully converged to the late-M, MH, and N intermediates, because the last points of IRC computations coincide with the structures separately obtained by geometry optimization. Hence, it is confirmed that proton transfer from Asp96 to Schiff base occurs in two-step reactions via the intermediate that generates an H3O+ molecule around the Ala215 backbone carboxyl oxygen. The hydrogen bond network changes slightly with the proton transfer. In the late-M state, the hydrogen bonds near Asp96 are short and stronger than those in the vicinity of the Schiff base. On the other hand, the hydrogen bonds near the Schiff base become shorter in the N state as shown in Figure 4 by thick dot lines. 3.3. MD Simulation with the Lipid Bilayer-ProteinWater System and Investigation of Surrounding Effects. We performed additional 1 ns MD simulations with the lipid bilayer-protein-water system (Figure 1, parts d and e) in order to investigate the effect of the lipid bilayer on stabilization of the water chain. A water chain from Asp96 to the Schiff base observed in NCO calculation is stably maintained during the 1 ns MD simulation (see the Supporting Information). This result supports the presence of a water chain from Asp96 to the Schiff base in the M state structure, which has been deduced from an N′ state X-ray crystallographic experiment. Furthermore, the consistency of the results of the two MD simulations with the lipid bilayer and the water solvent systems indicates that the lipid bilayer has little influence on stabilization of the water chain, although it has been reported that some compounds of the lipid bilayer have important roles in the photocycle of bacteriorhodopsin.42 We performed more accurate quantum chemical calculations with the truncated model system constructed from the average structure during the final 3 ps of the lipid-bilayer MD simulation in order to investigate the effect of surrounding amino acid residues. The potential energies of optimized structures for three kinds of intermediates were calculated with three different dielectric constants,  ) 1.0, 4.0, and 20.0 (Table 1). The 6-31G(d,p) basis set was used for the geometry optimization and energy estimation. In the case of  ) 1.0, the atom geometry starting from the MH structure converges to the late-M structure in optimization, and the potential energy level of the N state structure was much higher than that of the late-M intermediate. In the case of  ) 4.0, the potential energy of the N state structure is almost the same as that of the late-M state structure, and the potential energy of the MH state is high. We also obtained the transition structures of both the late-M f MH and

the MH f N conversions with  ) 4.0 and the 6-31(d,p) basis set by executing the optimized calculations in a similar way as described in the previous section. Activation energies are estimated to be 20.8 kcal/mol for the late-M f MH conversion and 4.6 kcal/mol for the MH f N conversion reaction. In the case of  ) 20.0, the change in energy is similar to that obtained by 3-21G** quantum chemical calculations. The potential energy level of the N state is lower than that of the late-M state. These results indicate that the surrounding electrostatic effect tightly regulates the transition reaction of the late-M-MH-N conversion and that the stability of the N state structure relative to that of the late-M state increases in proportion to the surrounding electrostatic intensity. This tendency can be easily understood if attention is paid to the electrostatic change in the proton acceptor and donor during proton transfer. Both Asp96 and the Schiff base are neutral in the late-M state, but they are charged in the N state. That is, the N state becomes advantageous in a high dielectric condition. Hence, we performed advanced quantum chemical calculations to assess the results of the quantum chemical calculations with the truncated model and to reveal which amino acids have important roles in the control of proton transfer. 3.4. Reevaluation with Advanced Quantum Chemical Calculations. A combination of classical molecular dynamics simulation and truncated model quantum chemical calculation has been widely used to elucidate a biological reaction mechanism. However, the error caused by truncation has been problematic, and several methods have recently been developed to overcome this problem. One of the promising methods for exact treatment of surroundings is the hybrid quantum mechanics/molecular mechanics (QM/MM) method. In this study, we used the integrated molecular orbital molecular mechanics (IMOMM) method implemented in the GAMESS package in order to investigate the surrounding effects and reevaluate the potential energy of each intermediate. First, we executed IMOMM calculation by using the same boundary as that used for the truncated model QM calculation. The rmsd values from the structures obtained by the truncated QM calculation in the late-M, MH, and N states are 0.48, 0.50, and 0.47 Å, respectively, indicating that there is no discrepancy between the truncated QM and the IMOMM calculations and that no QM atom overlaps at the position where side chains of the omitted amino acid residues exist. Hence, the truncated model structure used in this study is appropriate for investigating the potential energy surface of the proton transfer. Furthermore, the potential energies of three intermediates estimated by the IMOMM method with several different dielectric constants () are listed in Table 1. In the IMOMM calculation using 1.0 for , the potential energy of the N state is less stable than that of the late-M state. That is, the potential energy estimation is incompatible with the truncated QM calculation. In the IMOMM calculation with  ) 4.0 or 20.0, the order of the energetic stability among the three states is consistent with the results obtained by the truncated QM calculation. There are two possible reasons for the incompatibility in the case of  ) 1.0. One is that the electrostatic interaction is underestimated in the IMOMM calculation. Because the electrostatic interaction between the QM region and the MM region is calculated molecular mechanically, the charges on atoms in the MM region may not be sufficiently reflected in the QM wave function. One solution for this problem is to include all significant charged residues in the QM region. The other possible reason is that the dielectric constant of the inside of this protein is not so low or some factors

22810 J. Phys. Chem. B, Vol. 110, No. 45, 2006

Sato et al. (Figure 6c), the interaction being prominent with the acidic amino residues in the extracellular side (Glu194,204,Asp212). It is clear from the FMO results that Glu194, Glu204, and Asp212 interact with the active site electrically and that these residues affect the potential energy change along the proton transfer. 4. Discussion

Figure 6. Changes in the interresidue interaction energy due to the proton transfer from Asp96 to the Schiff base. (a) Interaction of Thr46 with other residues in bacteriorhodopsin. (b) Interaction of Asp96 with other residues. (c) Interaction of retinal + Lys216 with other residues. The horizontal axis represents the residue number, and vertical axis represents the difference in interaction energy of the N state and the late-M state.

other than the protein make the dielectric constant high, e.g., high ionic concentration in the channel. To clarify the important elements not taken into consideration in the IMOMM model, we executed single-point FMO calculation on the late-M and the N state structures. We used the atom geometry obtained by the IMOMM method as an input for the FMO calculation. The potential energy of the N state structure was found to be more stable than the late-M state by 47.8 kcal/ mol. Furthermore, the interresidue interaction energies were compared between two intermediates. Figure 6a shows the deviation of the interaction energy of Thr46 with each of the other residues in bacteriorhodopsin between the late-M and the N intermediates. The deviation was estimated by the equation ∆G ) G(N state) - G(late-M state); that is, a negative value means that the interaction with Thr46 becomes larger in the N state than in the late-M state. Figure 6, parts b and c, shows the interaction energies with Asp96 and retinal-Lys216, respectively. It was confirmed from this analysis that the interaction of Thr46 with Asp96 significantly stabilizes the N state compared to the late-M state. The interaction of Asp96 with other residues is obviously changed due to the proton transfer, but the prominent energy gain from Ala98 is compensated by the loss from other residues, and totally the sum of interaction energy is almost constant for Asp96. The interaction of retinal-Lys216 with other residues is also drastically changed due to the proton transfer

4.1. Proton Pathway in the Cytoplasmic Half Side. The water chain in the M state proposed here has not been detected in previous experimental studies.23,43,44 However, the locations of crystal water molecules near the Schiff base in those studies were different; that is, a crystal water molecule detected in some crystallographic studies was not observed in another study. In a recent X-ray crystallographic experiment on the N′ state,14 the presence of two water molecules in hydrophobic regions near the Shiff base was detected, and the existence of at least one water molecule directly making a hydrogen bond with the Schiff base was also detected in an IR spectroscopic study.45 In our theoretical approach, two water molecules were added to the region near the Schiff base to make the water chain. This assumption is plausible according to the detection of water molecules in the above experimental studies. Furthermore, it was revealed from our calculations that the proton translocation from Asp96 to Schiff base occurs via an intermediate in which an H3O+ molecule appears around Ala215. This intermediate was not detected in a recent QM/MD simulation study,15 in which Ala215 was treated as an MM region. Hence, it is important to treat Ala215 and Thr46 in the same level of calculation because both of them directly interact with the water chain in the cytoplasmic region. Our MD simulation further suggests the importance of Ala215 for the formation of the water chain, in which the carboxyl group of the main chain of Ala215 firmly supports the second water molecule from the Schiff base. This water molecule interacted with carboxyl oxygen atom at the distance of 2.7 Å on average (between oxygen atoms) in the MD simulation of the lipid bilayer system and corresponds to the WAT506 interacting with the carboxyl oxygen atom of Ala215 at the distance of 3.1 Å in the N′ state crystal structure.14 The water molecules forming the single water chain have not been observed in the crystallographic experiments on the ground state structures; thus, those water molecules are expected to enter from the cytoplasmic solvent before the late-M state. Furthermore, those water molecules stably formed the single water chain after the reprotonation of the Schiff base as seen in the optimized structure by the quantum chemical calculation (Figure 4c). Considering the consistency of our simulated N state structure with the N′ state crystal structure, the single water chain is assumed to be maintained after the reprotonation of Asp96 as well. Activation energies for the proton translocation from Asp96 to the Schiff base were calculated to be 9.7 for the first reaction step and 3.5 for the second reaction step (Figure 5), suggesting that the second reaction step occurs quickly, making it difficult for the H3O+ molecule to be detected in experiments. The existence of this intermediate will decrease the activation energy for the proton translocation and enhance the proton transfer from Asp96 to the Schiff base. To the best of our knowledge, this study is the first report to succeed in computational reproduction of the water chain inside a bacteriorhodopsin by MD simulation as shown in Figure 2. In both of NCO and PME MD simulations, the water chain was stably maintained in the cytoplasmic region. The water chain was also stably maintained in the following MD simulation

Proton Translocation in Bacteriorhodopsin with lipid-bilayer system. The existence and the stability of the water chain were reproduced by three different models’ MD simulations. It is assumed that the dynamics of the water chain is not strongly affected by the surrounding electrostatic effects, which is attributable to the considerably stable seven helices of bacteriorhodopsin. On the other hand, the patterns of the HBN in the two MD simulations were different, indicating that the HBN is affected by the surrounding electrostatic effects and easily rearranged. The low activation energy for the HBN rearrangement was confirmed by the ab initio quantum chemical calculations. As an initial three-dimensional structure of this study, we selected the structure presented by Sass et al.23 Several crystal structures have already been obtained by different research groups through X-ray crystallographic studies on the M-state of bacteriorhodopsin.23,43,44 Because those experimental crystal structures have a certain degree of difference among them and this causes controversy in the understanding of the proton transport mechanism, it will be important to perform MD simulation using each of the crystal data as an initial structure for the calculation and to compare the stabilities of the water chain inside the bacteriorhodopsin. Therefore, an investigation of several different crystal structures by precise MD simulations will be an important subject in our future study. 4.2. Effects of Surroundings. The treatment of surrounding effects has been controversial in the field of computational chemistry. The accuracy of results obtained by simulations depends on the model structure constructed and on the surrounding effects taken into account in computation. In this study, we elucidated the mechanism of proton translocation from Asp96 to the Shiff base in bacteriorhodopsin by the following multiple-step computations. (1) First, MD simulation was executed to observe the dynamics of water molecules in the cytoplasmic side and the pattern of interaction of water molecules with amino acid residues in the channel, which provided geometric data for the best choice of the truncated model employed in the following quantum chemical calculation. (2) A model structure was constructed on the basis of the atom geometry obtained by the first step, and then the potential energy changes with the HBN rearrangement and the following proton transfer were estimated by quantum chemical calculation. The dielectric constant  was set to 20.0 based on the results obtained by our previous ab initio quantum chemical calculation studies.11-13 (3) For reevaluation of the surrounding effects that had been treated with a continuum model in the truncated quantum chemical calculation, we executed additional advanced quantum chemical calculations, IMOMM and FMO calculations, which treat the effects explicitly. The effect of surroundings was included by the self-consistent reaction field method in the truncated model quantum chemical calculation. Computational results of the truncated model quantum chemical calculation, IMOMM calculation, and FMO calculation were qualitatively consistent. The potential energy decreased with the translocation of proton from Asp96 to the Schiff base. In computation by the FMO method, the difference between potential energy in the late-M state and that in the N state was calculated to be 47.8 kcal/mol. This value is too large and seems to be overestimated, judging from the experimental measurement that the excessive enthalpy at the K state is 12 kcal/mol compared to the ground state.46 This result is also inconsistent with previous kinetics experiments in which the biphasic decay of the M state was explained well with the assumption of steady equilibration of M and N followed by unidirectional decay of N to the initial state.47,48 Judging from the compatibility between the experimental observation and the potential energy estimation

J. Phys. Chem. B, Vol. 110, No. 45, 2006 22811 obtained by the truncated model quantum chemical calculation with  ) 4.0 and 6-31G(d,p) basis set (Table 1), the actual dielectric constant inside the protein is around  ) 4.0 and the energy difference of the intermediates is overestimated by FMO calculation due to the low level of the basis set and to the singlepoint calculation without geometry optimization. However, the following two observations are common between quantum chemical calculations with  ) 4.0 and 20.0: (1) proton transfer from Asp96 to the Schiff base occurs in a two-step reaction, and there exists an MH state that generates an H3O+ molecule around the Ala215 backbone atom and (2) the rate-limiting step is the late-M f MH conversion, which causes detachment of the hydrogen atom from Asp96 hydroxyl oxygen, and the activation energy for the MH f N conversion is considerably small. Despite the energetic overestimation, the FMO method and the truncated model quantum chemical calculations provided the important findings about the interactions among amino acid residues in bacteriorhodopsin. First, it is revealed by the FMO calculation that Thr46 has a hydrogen bond with Asp96 and stabilizes the N state. Because Asp96 is ionic in the N state, the hydrogen bond becomes stronger in the N state than that in the late-M state. This ionic effect, however, is expected to be canceled by the following reprotonation of Asp96 in the N′ state. The experimental findings that the decay of M state is more rapid in the T46V mutant than in the wild type47 is compatible with the temporary binding effect due to the hydrogen bond between Thr46 and Asp96, and the other factors such as the change of hydrophobic or steric interactions may have stronger effects on the stabilization of the N and the N′ states. Second, it is also revealed by the FMO calculation that three acidic residues, Glu194, Glu204, and Asp212 contribute to the stabilization of the N state. It is shown that the proton release to the extracellular region occurs before the Schiff base reprotonation in the E204Q mutant49 or under the condition at low pH < 6.50 Thus, the ionic state of Glu194 or Glu204 is considered not to be an essential factor for the reprotonation of the Schiff base, and the contribution of these residues to the stabilization of the N state will be minor. The proton-transfer reaction becomes 10-fold slower in the E204 mutant than in the wild type and the rate-limiting step also changes due to the mutation.49 These indicate that the ionic state of the acidic residue affect the energy profile of the enzymatic reaction by bacteriorhodopsin during the M f N f O conversion. Furthermore, our FMO results indicate that these three residues induce the unidirectional proton transfer from the cytoplasmic to the extracellular side cooperatively. The interactions among amino acid residues are considered to have an important effect to the energy dynamics of the reaction. In the present study, the interaction energy was estimated based on the results of FMO calculation. We mainly examined the contribution of amino acid residues in terms of electrostatic interaction. The other interactions such as hydrophobic interactions and steric interactions also have important effects to the protein function. Thus, it is very interesting to estimate these interaction energies separately and clearly reveal the effect of each component. The highly accurate and costdemanding methods such as QM/MM and FMO will be more feasible in the future for the investigation of reaction mechanisms and to refine the energy estimation. In addition to this, a combinatorial approach using classical MD simulation, quantum chemical calculation, and advanced methods will also be a promising strategy because of its efficiency within limited computational resources and more importantly from the point

22812 J. Phys. Chem. B, Vol. 110, No. 45, 2006 of view that the reliability of computational analysis increases by integrating a wide range of information obtained from different methods as described here. 5. Conclusions Several experimental studies have suggested that proton transfer from Asp96 to the Schiff base occurs as follows: (1) a hydrogen bond between Asp96 and Thr46 breaks and (2) a proton is transferred via the water chain of four water molecules from Asp96 to the Schiff base.44,51-53 We have successfully reproduced these observations by the MD simulations and the following truncated model QM simulations. We further propose the presence of an intermediate having a stable H3O+ ion near the carboxyl oxygen of Ala215 in the proton-transfer process. The rate-limiting step is the breaking of the OH bond of Asp96, which corresponds to the late-M - MH conversion. Examinations of the effects of surrounding amino acid residues by using advanced quantum chemical calculations have further revealed that the energy change with the proton transfer in bacteriorhodopsin is strongly affected by the electrostatic environment created by charged amino acid residues, which are mainly the extracellular titrating amino acid residues, Glu194, Glu204, and Asp212. Acknowledgment. The authors thank the staff of the Research Center for Computational Science, Okazaki, and Institute of Media and Information Technology, Chiba University. This work was partly supported by a Grant-in-Aid for Center Of Excellence (COE) research from the Ministry of Education, Science, Sport, and Culture, Japan. This work was also supported by the Japan Science and Technology Agency. Supporting Information Available: Figure of the optimized structures of two transition states (Figure S1) and the results of molecular dynamics simulation of bacteriorhodopsin with the lipid bilayer-protein-water system (Figure S2). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Oesterhelt, D.; Stoeckenius, W. Proc. Natl. Acad. Sci. U.S.A. 1973, 70, 2853. (2) Heyes, C. D.; El-Sayed, M. A. J. Phys. Chem. B 2003, 107 (44), 12045. (3) Mathies, R. A.; Lin, S. W.; Ames, J. B.; Pollard, W. T. Annu. ReV. Biophys. Biophys. Chem. 1991, 20, 491. (4) Lanyi, J. K. Nature (London) 1995, 375, 461. (5) Lanyi, J. K. J. Struct. Biol. 1998, 124, 164. (6) Kandori, H. Biochim. Biophys. Acta 2004, 1658, 72. (7) Shibata, M.; Tanimoto, T.; Kandori, H. J. Am. Chem. Soc. 2003, 125, 13312. (8) Hayashi, S.; Ohmine, I. J. Phys. Chem. B 2000, 104 (45), 10678. (9) Cao, Y.; Varo, G.; Chang, M.; Ni, B.; Needleman, R.; Lanyi, J. K. Biochemistry 1991, 30, 10972. (10) Dencher, N. A.; Sass, H. J.; Buldt, G. Biochim. Biophys. Acta 2000, 1460, 192. (11) Murata, K.; Fujii, Y.; Enomoto, N.; Hata, M.; Hoshino, T.; Tsuda, M. Biophys. J. 2000, 79, 982. (12) Murata, K.; Hoshino, T.; Sato, Y.; Hata, M.; Tsuda, M. Chems Bio Inf. J. 2002, 2, 97. (13) Murata, K.; Hoshino, T.; Sato, Y.; Hata, M.; Tsuda, M. J. Mol. Struct.: THEOCHEM 2003, 664/665, 125. (14) Schobert, B.; Brown, L. S.; Lanyi, J. K. J. Mol. Biol. 2003, 330, 553. (15) Lee, Y. S.; Krauss, M. J. Am. Chem. Soc. 2004, 126, 2225. (16) Sato, Y.; Hata, M.; Neya, S.; Hoshino, T. Protein Sci. 2005, 14, 183. (17) Ode, H.; Ota, M.; Neya, S.; Hata, M.; Sugiura, W.; Hoshino, T. J. Phys. Chem. B 2005, 109, 565. (18) Fujii, Y.; Okimoto, N.; Hata, M.; Narumi, T.; Yasuoka, K.; Susukita, R.; Suenaga, A.; Futatsugi, N.; Koishi, T.; Furusawa, H.; Kawai, A.; Ebisuzaki, T.; Neya, S.; Hoshino, T. J. Phys. Chem. B 2003, 107, 10274.

Sato et al. (19) Yamanaka, K.; Okimoto, N.; Neya, S.; Hata, M.; Hoshino, T. J. Mol. Struct.: THEOCHEM 2006, 758, 97. (20) Shoemaker, J. R.; Burggraf, L. W.; Gordon, M. S. J. Phys. Chem. A 1999, 103, 3245. (21) Kitaura, K.; Sawai, T.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 312, 319. (22) Inadomi, Y.; Nakano, T.; Kitaura, K.; Nagashima, U. Chem. Phys. Lett. 2002, 364, 139. (23) Sass, H. J.; Buldt, G.; Gessenich, R.; Hehn, D.; Neff, D.; Schlesinger, R.; Berendzen, J.; Ormos, P. Nature (London) 2000, 406, 649. (24) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (25) Kates, M.; Moldoveanu, N.; Stewart, L. C. Biochim. Biophys. Acta 1993, 1169, 46. (26) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Wang, J.; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER7; University of California: San Francisco, CA, 2002. (27) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049. (28) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (29) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B: Condens. Matter 1988, 37, 785. (30) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A. J.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; AlLaham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98; Gaussian, Inc.: Pittsburgh, PA, 1998. (31) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1357. (32) Tajkhorshid, E.; Baudry, J.; Schulten, K.; Suhai, S. Biophys. J. 2000, 78, 683. (33) Narumi, T.; Susukita, R.; Ebisuzaki, T.; McNiven, G.; Elmegreen, B. Mol. Simul. 1991, 21, 401. (34) Narumi, T.; Susukita, R.; Koishi, T.; Yasuoka, K.; Furusawa, H.; Kawai, A.; Ebisuzaki, T. Proceeding of the 2000 ACM/IEEE Conference on Supercomputing, Dallas, Texas, Nov 4-10, 2000, p 54. (35) Onsager, L. J. Am. Chem. Soc. 1936, 58, 1486. (36) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (37) Ponder, J. W.; Richards, F. M. J. Comput. Chem. 1987, 8, 1016. (38) Sayle, R. A.; Milner-White, E. J. Trends Biochem. Sci. 1995, 20, 374. (39) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33. (40) Bode, B. M.; Gordon, M. S. J. Mol. Graphics Modell. 1998, 16, 133. (41) Dioumaev, A. K.; Brown, L. S.; Needleman, R.; Lanyi, J. K. Biochemistry 2001, 40, 11308. (42) Joshi, M. K.; Dracheva, S.; Mukhopadhyay, A. K.; Bose, S.; Hendler, R. W. Biochemistry 1998, 37, 14463. (43) Takeda, K.; Matsui, Y.; Kamiya, N.; Adachi, S.; Okumura, H.; Kouyama, T. J. Mol. Biol. 2004, 341, 1023. (44) Lanyi, J.; Schobert, B. J. Mol. Biol. 2002, 321, 727. (45) Fischer, W. B.; Sonar, S.; Marti, T.; Khorana, H. G.; Rothschild, K. J. Biochemistry 1994, 33, 12757. (46) Birge, R. R.; Cooper, T. M.; Lawrence, A. F.; Masthay, M. B.; Zhang, C. F.; Zidovetzki, R. J. Am. Chem. Soc. 1991, 113 (11), 4327. (47) Brown, L. S.; Zimanyi, L.; Needleman, R.; Ottolenghi, M.; Lanyi, J. K. Biochemistry 1993, 32, 7679. (48) Fukuda, K.; Kouyama, T. Biochemistry 1992, 31, 11740. (49) Misra, S.; Govindjee, R.; Ebrey, T. G.; Chen, N.; Ma, J. X.; Crouch, R. K. Biochemistry 1997, 36, 4875. (50) Cao, Y.; Brown, L. S.; Needleman, R.; Lanyi, J. K. Biochemistry 1993, 32, 10239. (51) Luecke, H.; Schobert, B.; Cartailler, J. P.; Richter, H. T.; Rosengarth, A.; Needleman, R.; Lanyi, J. K. J. Mol. Biol. 2000, 300, 1237. (52) Lanyi, J. K.; Schobert, B. J. Mol. Biol. 2003, 328, 439. (53) Lanyi, J. K. Biochim. Biophys. Acta 2005.