Molecular Dynamics Simulation Study on Stabilities and Reactivities of

Apr 16, 2008 - Toshio Asada,*, Shigeru Nagase,, Kichisuke Nishimoto, andShiro Koseki ... Shigayuki Yagi , Hiroyuki Nakazumi , and Takeshi Matsushita...
1 downloads 0 Views 517KB Size
Download PDF
5718

J. Phys. Chem. B 2008, 112, 5718-5727

Molecular Dynamics Simulation Study on Stabilities and Reactivities of NADH Cytochrome B5 Reductase Toshio Asada,*,† Shigeru Nagase,‡ Kichisuke Nishimoto,§ and Shiro Koseki† Department of Chemistry, Faculty of Science, Osaka Prefecture UniVersity, 1-1 Gakuen-cho, Naka-ku, Sakai, Osaka 599-8531, Japan, Department of Theoretical and Computational Molecular Science, Institute for Molecular Science, Myodaiji, Okazaki 444-8585, Japan, and Department of Chemistry, Faculty of Science, Osaka City UniVersity, Sumiyoshi-Ku, Osaka 558-8585, Japan ReceiVed: October 19, 2007; In Final Form: February 16, 2008

Binding free energies between coenzyme (FAD and NADH) and the apoenzyme of NADH-cytochrome b5 reductase (b5R) were estimated by applying the continuum Poisson-Boltzmann (PB) model to structures sampled from molecular dynamics simulations in explicit water molecules. Important residues for the enzymatic catalysis were clarified using a computational alanine scanning method. The binding free energies calculated by applying an alanine scanning method can successfully reproduce the trends of the measured steady-state enzymatic activities kcatNADH/KmNADH. Significant decreases in the binding free energy are expected when one of the four residues Arg91, Lys110, Ser127, and Thr181 is mutated into Ala. According to the results of the molecular dynamics simulation, Thr181 is considered to be one of the key residues that helps NADH to approach the isoalloxazine in FAD. Finally, we have constructed very simplified model systems and carried out density functional theory calculations using B3LYP/LANL2DZ//ROHF(or RHF)/LANL2DZ level of theory in order to elucidate a realistic and feasible mechanism of the hydride-ion transfer from NADH to FAD affected by HEME(Fe3+) as an electron acceptor. Our calculated results suggest that the electron and/or hydrideion transfer reaction from NADH to FAD can be accelerated in the presence of HEME(Fe3+).

1. Introduction Flavoproteins (flavoenzymes) catalyze various important biochemical reactions, such as oxidation (oxidase, oxygenase, dehydrogenase),1,2 reduction (reductase),3,4 and electron transfer (both one-electron and two-electron transfers). Figure 1 depicts the structure of the active site in NADH-cytochrome b5 reductase (b5R), which is a member of the flavoenzyme family5,6 of dehydrogenases-electron transferases that mediates the transfer of electron from NADH to cytochrome b5 (b5). Crystallographic analysis has revealed that b5R consists of two domains that are separated by a large interdomain cleft.7,8 The N-terminal and C-terminal domains comprise the FAD and NADH binding sites, respectively. It is suggested that the large interdomain cleft permits the nicotinamide moiety of NADH to access the re side of the isoalloxazine ring of the flavin without inducing a conformational change in the C-terminal region.7 An outline of the catalytic cycle of the solubilized catalytic domain of pig b5R (pb5R) has already been proposed experimentally by Iyanagi et al.4 In their mechanism, two electrons are transferred from NADH to FAD by hydride-ion (H-) transfer at the initial stage of the catalytic cycle. The two-electron reduced enzyme-NAD+ complex (E-FADH--NAD+) then transfers two electrons to two one-electron acceptors one by one, after which the reduced enzyme returns to the oxidized state. Furthermore, Strittmatter has suggested that the reduction of FAD by NADH is the rate-limiting step in electron transfer * Corresponding author. E-mail: [email protected]. † Osaka Prefecture University. ‡ Institute for Molecular Science. § Osaka City University.

Figure 1. Structure of the active site in b5R. Nicotine amide moiety in NADH and isoalloxazine ring in FAD are represented by ball and stick model. The red circle indicates the hydrogen atom that moves from NADH to FAD in the catalytic cycle. The gray circle indicates the C-terminal of b5R.

catalyzed by b5R.9-11 However, the experimentally determined crystal structure of pb5R was inconsistent with its biochemistry; for example, Lys83 and Cys245 were not located at the active site. Later, Bewley et al. determined the structure of the soluble domain of rat b5R at a resolution of 2.3 Å.8 In the rat model, Lys110 (Lys83 in the pig model) lies in front of the isoalloxazine

10.1021/jp7101513 CCC: $40.75 © 2008 American Chemical Society Published on Web 04/16/2008

Stabilities and Reactivities of NADH Cytochrome B5 Reductase ring, where it is oriented toward the re face of the cofactor. They have mentioned that this side chain is free to interact with an incoming NADH molecule as suggested by its biochemistry. Despite the presence of many experimental reports regarding b5R, theoretical investigations of this enzymatic reaction have yet to be carried out. To this end, we have calculated the binding free energies between coenzyme (cofactor FAD and NADH) and the apoenzyme of b5R by applying the continuum PoissonBoltzmann (PB) model12 to structures sampled from molecular dynamics (MD) simulations in explicit solvents. First, we will discuss the binding free energies in the wild type of rat b5R structure. Some important residues for the binding free energies will then be clarified using the computational alanine scanning method proposed by Massova et al.13 Although their original method can usually reproduce observed experimental trends for the binding free energies, it sometimes fails when applied to the mutation of buried bulky residues that the solvent is not accessible. We used a modified method, which improves upon their method. Significant decreases in the binding free energy are expected when one of the four residues Arg91, Lys110, Ser127, and Thr181 is mutated into Ala. Finally, the possible mechanism of H- transfer will be discussed by the results of the density functional theory (DFT) calculations using two simplified model systems: the first is from nicotine amide (NAH), which is the model of NADH, to 10-methyl isoalloxazine (10-MIA), which is the model of FAD; the second is from NAH to 10-MIA with (HEME(Fe3+) + imidazole (IMI)), which is an electron acceptor. In the present paper, we obtained important results for the activation energy of the H- transfer reaction using model systems. That is, when we calculate NAH+10-MIA system, the activation energy is 20.9 kcal/mol. But when we calculate NAH+10-MIA+(HEME+IMI), the activation energy drastically decreases to be 11.3 kcal/mol. 2. Method 2.1. Model Structure Setting. The X-ray crystal structure of rat b5R8 reported in the Protein Data Bank (PDB) was used as the starting structure for the molecular dynamics (MD) simulations. The PDB ID code is 1ib0, with a resolution of 2.3 Å. All crystal water molecules have been included in our calculations. Hydrogen atoms were added to the heavy atoms of the crystal structure using the Leap module of the Amber 7.0 package.14 The numbering of the atoms and the atom types of FAD and NADH is shown in Figure 2 and in Figure 1S in the Supporting Information, respectively. Partial atomic charges of FAD and NADH were calculated using RESP software in the Amber package, in which atomic charges are fitted to reproduce the electrostatic potential surrounding the molecule. It is identical to the approach used to derive molecular mechanical electrostatic charges in AMBER99 force field.15,16 The total atomic charges of the substrate were limited to -2.0e for both FAD and NADH using Lagrangian constraints. We cannot determine whether these phosphates are completely deprotonated or not in the protein environment from the present available X-ray results. There are some cases for protonation of diphosphate moiety. At the present, we assumed that both FAD and NADH contain diphosphates with a minus two charge based on the experimentally proposed model by Kimura et al.17 All RESP charges and bond angle parameters for FAD and NADH not included in the AMBER99 force field parameters are listed in Tables 1S and 2S in the Supporting Information. 2.2. Molecular Dynamics Simulations. Molecular dynamics simulations were performed using the SANDER program in the Amber 7.0 package. The complex b5R was surrounded by

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5719

Figure 2. Numbering of atoms: (a) FAD, (b) NADH.

TIP3P18 cap water molecules within a radius of 35.0 Å from the geometrical center of the isoalloxazine ring of FAD. The geometry of each system was initially optimized with the steepest descent followed by the application of the conjugate gradient energy minimization algorithm. The long-range nonbonded interactions were truncated at 20 Å. In our simulations, only residues within 18 Å from N5 atom of FAD, which is called the belly region, were allowed to move, and the coordinates of the rest of the residues were frozen. The SHAKE algorithm was used to fix all bond lengths at equilibrium distances and to remove bond stretching, which is the fastest of all atomic displacements. The system was gradually heated to 300 K and equilibrated for an additional 100 ps with the a 1.0 fs time step, followed by a 2 ns period of dynamics simulation and data collection. The temperature control was carried out using the Berendsen coupling scheme, which uses separate scaling factors for atoms of the solute and the solvent to overcome the “cold solute/hot solvent” problem.19 2.3. Binding Free Energy Calculation. The binding free energies were evaluated according to the strategy of MM-PB/ SA,13,20,21 where MM and SA stand for molecular mechanics and surface area approximation, respectively. The outlines are summarized below. To use MM-PB/SA for the binding free energy analysis, 80 snapshots were taken at 25 ps intervals from MD trajectories of the b5R complex. To evaluate free energies, the conformational averaging of these snapshots is used. In this method, all solvent and crystal water molecules are not included explicitly. The binding free energy ∆Gbind is written as

∆Gbind ) ∆EMM + ∆Gsolv - T∆S

(1)

5720 J. Phys. Chem. B, Vol. 112, No. 18, 2008

Asada et al.

Figure 3. Temporal development of the root-mean-square displacement, RMSD; RMSD for protein backbone atoms of b5R in the belly region is shown in green color. Thr181 in red. R helix from Gly182 to Leu194 in black. β sheet ranging from Val167 to Gly180 in blue. The interatomic distance between N5 atom in FAD and H14 atom in NADH, R(N5-H14), is also shown for comparison in purple color.

where EMM, Gsolv, T, and S are the molecular mechanical energy, the solvation free energy, the absolute temperature, and the entropy of the solute, respectively. As Massova has pointed out,20 the relative contribution of the change in conformational entropy to ∆Gbind is negligible for the mutations. Thus, we will neglect the contribution of the last term in eq 1 in the following discussions since we are interested in the relative values of ∆Gbind for the wild type and some mutants of b5R. The molecular mechanical energy EMM can be approximated as the summation of several physical terms

EMM ) Ebind + Eangle + Edihed + Eelectrostatic + EvdW (2) where Ebind, Eangle, and Edihed are bond, angle, and dihedral angle energies, and Eelectrostatic and EvdW are electrostatic and van der Waals (VDW) interaction energies, respectively. Nonbonded cutoff was not used in the energy component analysis. Further, the solvation free energy Gsolv consists of the electrostatic contribution GPB and the nonpolar contribution Gnp. The electrostatic contribution GPB was calculated by solving the Poisson-Boltzmann equation using the finite-difference method proposed by Nicholls et al.12 We have used the PARSE parameter set22 and Cornell et al. charges15,16 with a dielectric constant of 80 for water22 and 1 for the solute.23 This value of the dielectric constant for the solute is consistent with our use of a nonpolarizable molecular mechanics force field; in this model, the solute response to charge fluctuations is estimated through explicit averaging of the conformations, rather than by use of an interior dielectric constant greater than unity.24 The solvent dielectric constant of 80 is widely used for roomtemperature water. The contribution of the nonpolar term Gnp was calculated from

Gnp ) γSA + b

(3)

where γ ) 0.00542 kcal/[mol Å2], b ) 0.92 kcal/mol using the surface area estimation of the program MSMS.25 The probe radius used for calculatating the solvent accessible surface area was set to the conventional 1.4 Å.

2.4. Alanine Scanning Method. The alanine scanning method is a convenient method to evaluate the component of stabilization effects for the binding free energies between surrounded residues and the active site. We have applied two models to analyze the effects of point mutation on the binding free energy between the apoenzyme and both FAD and NADH. One is the conventional alanine scanning method,13,24 which is referred to as model A in the following discussion, in which the Cγ atom of a given residue is simply replaced by a hydrogen atom that is located at a distance of 0.94 Å from the Cβ atom in the direction from the Cβ atom to the Cγ atom. Accordingly, the topological parameters have been changed to Ala. The present study uses a modified model, which is referred to as model B, in which the relative dielectric constant of the cavity produced by the mutation to Ala is considered to be equivalent to the internal dielectric constant (1 in our calculation) for the PB calculation. In other words, when an original residue is mutated to Ala, the space occupied by a truncated side-chain is considered to be filled by the continuum medium with a dielectric constant equal to the internal dielectric constant in model B, regardless of the fact that the dielectric constant of the solvent is used (80 for water) in model A. In both models, the single trajectory obtained from the original MD simulation was used for residue mutations, after which MM-PB/SA calculations were carried out for these mutated trajectories to calculate the binding free energies for each mutant. If we carry out different MD runs for each mutant, the large computational time is required and component analysis for the binding energy become more complicated because of respective thermal fluctuations. 2.5. Density Functional Theory (DFT) Calculations. To elucidate a realistic and feasible reaction mechanism of Htransfer from NADH to FAD using DFT calculations, we used 10-methyl isoalloxazine (10-MIA), nicotinamide (NAH), HEME(Fe2+/Fe3+) and imidazole (IMI) instead of FAD, NADH, protoheme IX and His, respectively. All DFT calculations were carried out by using the Gaussian 03 program.26 Owing to the limitations of our computational resources, the discussion for the present model of reactive systems is based on results obtained by using B3LYP/LANL2DZ//ROHF(or RHF)/ LANL2DZ level of theory. 3. Results and Discussion 3.1. Structure and Molecular Dynamics Simulation. MD simulations of the solvated b5R were performed with a spherical cap of explicit water centered at isoalloxazine of FAD. Figure 3 shows the temporal development of the rmsd for all protein backbone atoms in b5R, Thr181, the β sheet ranging from Val167 to Gly180, and the R helix ranging from Gly182 to Leu194, as well as the interatomic distance between N5 atom in FAD and H14 atom in NADH, R(N5-H14). Stable trajectories of the complexes could be obtained as demonstrated by the steady rmsd for protein backbone atoms in b5R in the belly region. In the initial structure at t ) 0.0 ps, the direct H- transfer from NADH to FAD is apparently prevented by the long R(N5H14) distance of 10.4 Å. In the 2 ns trajectory, R(N5-H14) significantly fluctuates over the range from 6 to 12 Å, whereby the shortest R(N5-H14) distance is fluctuating around 7.08.0 Å. It is interesting that the R(N5-H14) fluctuation correlates with the displacement of the residue Thr181. Thr181 is located at the hinge region between the R helix and the β sheet, and it has a hydrogen bond with the phosphate moiety in NADH (See Figure 4b). The rmsd of the R helix is almost 0.5 Å, which is

Stabilities and Reactivities of NADH Cytochrome B5 Reductase

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5721

Figure 4. Schematic representation of the binding site in (a) FAD and (b) NADH, as given by the snap shot after 2 ns MD simulation. Rectangles and ellipsoids represent residues involved in hydrogen bonding (broken) and vdW interactions (wide broken), respectively.

on the same order of magnitude as the one measured for the protein backbone. Although the β sheet shows that the rmsd is large in comparison with the R helix, its fluctuation is more steady than that of the R helix. Apparently, the rmsd of Thr181 is significantly larger than that of the other elements. It can be

assumed that the small fluctuations, such as rotation, of the R helix connected with Thr181 causes destabilization of the hydrogen bond (with strong electrostatic interaction) between Thr181 and the negatively charged phosphate moiety in NADH. Thus, the displacement of NADH correlates with that of Thr181.

5722 J. Phys. Chem. B, Vol. 112, No. 18, 2008

Asada et al.

TABLE 1: Energetic Analysis of the Binding Free Energy for the Wild Type Complex: Mean Values and Standard Errors ∆Eelectrostatic ∆EvdW ∆EMM ∆Gnp ∆GPB ∆Gsolv ∆Gelec,tot ∆Gtotal

mean (kcal/mol)

se

-921.0 -147.7 -1068.7 -13.0 1023.9 1010.9 102.9 -57.8

2.7 0.6 2.6 0.0 2.5 2.5 1.0 0.8

TABLE 2: Computational Alanine-Scanning Mutagenesis Results of b5R (∆∆ ) ∆mutant - ∆wildtype) for Important Eight Residues: Mean Values and Standard Errorsa contribution

Arg91Ala

Tyr93Ala

Thr94Ala

Lys110Ala

mean

mean

mean

mean

std

se

se

se

∆∆Eelectrostatic 211.4 0.7 1.8 0.2 1.5 0.1 223.5 1.0 ∆∆EvdW 3.4 0.2 6.4 0.1 0.3 0.0 0.1 0.2 ∆∆EMM 214.8 0.7 8.2 0.2 1.8 0.1 223.6 0.9 ∆∆Gnp 0.2 0.0 -0.3 0.2 0.0 0.1 -0.1 0.0 ∆∆GAPBb -206.6 0.6 -11.7 0.3 -2.0 0.1 -207.3 0.7 ∆∆GBPBc -190.8 0.5 -5.1 0.1 -2.0 0.1 -192.1 0.6 ∆∆GAtotalb 8.4 0.2 -3.8 0.3 -0.2 0.1 16.2 0.3 ∆∆GBtotalc 24.2 0.3 2.8 0.2 -0.2 0.1 31.4 0.4 Tyr112Ala

Lys125Ala

Ser127Ala

Thr181Ala

mean

mean

se

mean

std

mean

∆∆Eelectrostatic 20.6 0.3 138.5 1.1 ∆∆EvdW 8.7 0.2 0.7 0.0 ∆∆EMM 29.3 0.2 139.2 1.2 ∆∆Gnp 0.2 0.0 0.0 0.0 ∆∆GAPBb -26.7 0.2 -137.1 1.2 ∆∆GBPBc -15.5 0.1 -136.3 1.1 ∆∆GAtotalb 2.8 0.3 2.1 0.1 ∆∆GBtotalc 14.0 0.2 2.9 0.1

17.0 -1.5 15.5 0.0 -8.7 -9.6 6.8 5.9

0.2 16.6 0.3 0.2 2.9 0.2 0.2 19.5 0.3 0.0 -0.1 0.0 0.1 -11.2 0.3 0.1 -11.3 0.2 0.1 8.2 0.4 0.1 8.1 0.3

contribution

se

se

a Results for less important residues are summarized in the Supporting Information, Table 3S. Values are in kcal/mol. b Model A:  ) 80. c Model B:  ) 1.

We could not find statistically meaningful correlations between the distance R(N5-H14) and structural fluctuations for other residues. Figure 4 illustrates the schematic representation of the active sites corresponding to the snapshot following the 2 ns MD simulation. The negatively charged diphosphate moiety of FAD forms the bulk of the favorable electrostatic interactions between FAD and the surrounding residues. Indeed, this diphosphate group interacts with two positively charged residues, Arg91 and Lys125, with an additional hydrogen bond with Ser127. On the other hand, the diphosphate moiety of NADH is bounded with only one positively charged residue of Lys110 through electrostatic interaction, with additional hydrogen bonds with Tyr112, Thr181, and Gln210. Therefore, the rmsd of NADH becomes larger than that of FAD in our MD simulations (data not shown). As Kimura et al. have pointed out,17 Arg91 and Lys125 neutralize the negative charge on the phosphate groups of FAD, and assist in the binding of NADH by reducing the electrostatic repulsion between the negative charges of the phosphates of NAD. Thus, the negatively charged phosphate of the bound NADH is positioned near the phosphate of FAD. It should be noted that there is also the strong electrostatic interaction between Lys110 and the diphosphate group of NADH. There is also a hydrogen bond between the backbone NH moiety of Lys110 and the O12 atom in FAD, as seen in Figure 4a. Thus, Lys110 is the important bridging residue between FAD and

NADH, which is located in the FAD domain in b5R. In fact, in the crystal structure of the rat b5R, Lys110 lies in front of the isoalloxazine ring, where it is oriented toward the re face of the cofactor. Therefore, this side chain is free to interact with the NADH molecule, as suggested by its biochemistry. 3.2. Energy Analysis of b5R using Alanine Scanning Method. To investigate the binding free energies ∆Gbind between the apoenzyme and the coenzyme (cofactor FAD and NADH), we collected the Cartesian coordinates of all molecules over 2 ns periods at a temperature of T ) 300 K, after which MM-PB/SA and alanine scanning analyses were performed. The configuration used for the MM-PB/SA calculations includes 80 equally spaced snapshots taken over the span of one MD simulation. All energy analyses were performed for a single MD trajectory on the empty enzyme and on the enzyme with both FAD and NADH bound. The calculated results for each component of ∆Gbind for the wild type are listed in Table 1. The favorable formation of the complex, (∆Gtotal ) -57.8kcal/ mol), is determined by the contribution of the VDW forces (∆EvdW ) -147.7kcal/mol) and the nonpolar contribution to solvation (∆Gnp ) -13.0kcal/mol). Both terms are nonelectrostatic components. On the other hand, the electrostatic contribution (∆Gelec,tot ) 102.9 kcal/mol) acts as the unfavorable factor for the binding free energy as demonstrated by numerous studies,24,27-31 because of the fact that this contribution comes from both MM electrostatic energy (∆Eelectrostatic ) -921.0 kcal/ mol) and the electrostatic contribution to solvation (∆GPB ) 1023.9 kcal/mol). The unfavorable change in the electrostatics of the solvation is mostly, but not fully, compensated by the favorable electrostatics within the resulting complex. Although binding free energies have not yet been experimentally measured, the values of the catalytic rate constant kcat and the Michaelis constant Km for NADH have been reported by Kimura et al.17,32 and by Bewley et al.8 using the site-directed mutagenesis approach for some residues surrounded FAD. These experimental data become a good reference for understanding the role each residue plays in the formation of the catalytic properties in the tertiary structure of the solubilized catalytic domain in b5R. Thus, we have carried out the computational alanine scanning analysis for 35 surrounded residues within 8 Å from NADH and FAD in order to clarify the importance of the residues for the binding free energy in the complex. Table 2 lists the calculated components of the relative binding free energies based on the wild-type enzyme for eight residues. The less important results of the rest are summarized in the Supporting Information, Table 3S. The negative numbers corresponding to ∆∆GA/B total indicate the favorable substitution, and the positive numbers indicate the unfavorable substitution. These results show significant decreases in the binding free energy are expected when one of the four residues Arg91, Lys110, Ser127, and Thr181 is mutated to Ala. Tyr112 is characterized as an important residue in a model B since ∆∆ GBtotal is large in contrast with a model A. The positively charged residue Arg91 determines the favorable electrostatic interaction with FAD and NADH, which corresponds to remarkable decreases in the binding free energy for the Arg91Ala mutant. As seen in Figure 4, the negative charge of the diphosphate moiety of FAD interacts with the two positive residues of Arg91 and Lys125; however, Arg91 is considered to be a more important residue than Lys125 in our calculations. In the crystal structure, since Lys125 is more exposed to the solvent region than Arg91, screening effects of electrostatic interactions become larger for Lys125 than for Arg91. Indeed, Kimura et al. have suggested, on the basis of experiment results,

Stabilities and Reactivities of NADH Cytochrome B5 Reductase

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5723

TABLE 3: Comparison between Calculated Relative Binding Free Energies (∆∆ ) ∆mutant - ∆wildtype) and Experimental Steady-State Enzymatic Activities for Six Mutants enzyme WT Thr94Ala Lys125Ala Tyr93Ala Ser127Ala Arg91Ala Lys110Ala a

KmNADH (µM)

KcatNADH (s-1)

KcatNADH/KmNADH (µM-1 s-1)

∆∆GAtotal (kcal/mol)

∆∆GBtotal (kcal/mol)

3.1 ( 0.3a 6 ( 1b 2.9 ( 0.2c 4.4 ( 0.3a 4.2 ( 0.4a 26 ( 2a 104 ( 18a 380 ( 25b

1100 ( 84a 800 ( 8b 984 ( 30c 1060 ( 19a 180 ( 3a 435 ( 8a 403 ( 27a 270 ( 25b

354 ( 61a 133 ( 10b 339c 240 ( 21a 42.9 ( 4.8a 16.7 ( 1.6a 3.9 ( 0.9a 0.7 ( 0.5b

-0.2 2.1 -3.8 6.8 8.4 16.2

-0.2 2.9 2.8 5.9 24.2 31.4

See ref 17. b See ref 8. c See ref 32.

Figure 5. Numbering of atoms: (a) 10-MIA, (b) NAH.

that Arg91 is critical for the affinity of Pb5R to NADH and NAD+ in the electrostatic interactions. Because the positively charged Lys110 is the important residue as shown above between FAD and the two negatively charged phosphate moieties in NADH, the contribution to the binding free energy ∆∆GA/B total of Lys110 is largest in the 35 mutants we have investigated. Lys110 is considered to be one of the most important residues for the stabilization of this enzyme, which is co-incident with the experimental site-directed mutagenesis.8 In addition, the two residues Ser127 and Thr181 have a hydroxyl group that makes favorable hydrogen bonds with the phosphate moieties in FAD and NADH, thus the changes in the binding free energies ∆∆GA/B total in their mutants have also large positive values. Tyr112 has the stable VDW interaction with adenine (RA3 in Figure 4a) in FAD and also has a hydrogen bond with the diphosphate moiety in NADH. Thus, this residue is considered to be an important one. In this study, we have examined two extreme cases for alanine mutation by using models A and B. In the protein, the value of  in the new cavity might be very small, which resembles to results of a model B. The steady-state enzymatic activities kcatNADH/KmNADH are compared to our alanine scanning results in Table 3, where the smaller kcatNADH/KmNADH value for the mutant than for the wildtype is considered to be an important residue for the enzymatic catalysis. The experimental data suggest that the Lys110Ala mutant entails a significant decrease in the enzymatic activities, followed by the Arg91Ala, Ser127Ala, Tyr93Ala, Lys125Ala,

and Thr94Ala mutants, whose trends were successfully reproduced within an error for ∆∆GBtotal by using the computational alanine scanning analyses. On the other hand, using the conventional model A with the alanine scanning method, Tyr93Ala mutant becomes the favorable mutation with ∆ ∆GAtotal ) - 3.8 kcal/mol, in which the negative change in the solvation free energy ∆∆GAPB overcomes the unfavorable electrostatic and VDW interactions, ∆∆Eelectrostatic and ∆∆EvdW, which contradicts with the experimental results. Kimura et al.17 have pointed out that Tyr93 makes an aromatic contact with the si-face of the isoalloxazine ring through a direct hydrogen bond to the 4′-hydroxyl group of the ribityl moiety of the FAD cofactor (see Figure 4a). In the case of the Tyr93Ala mutant, the weakening of the VDW interaction, ∆∆EvdW ) 6.4 kcal/ mol, indicates that the VDW contact between the phenyl ring of Tyr93 and the isoalloxazine ring of FAD in the wild-type enzyme plays an important role in the stability of the complex as suggested by Kimura et al.17 Usually, we can easily calculate the binding free energies for mutants by means of the alanine scanning approach, which uses a single trajectory obtained for a wild-type MD simulation to generate the structures of monomers and dimers. However, it should be noted that the alanine scanning technique used here cannot explain the conformational changes caused by the mutation, whereas the experimental site-directed mutagenesis leads to both the change in the inter-residue interactions and the conformational changes. Because Tyr93 has the buried bulky side chain, it is a reasonable assumption that the solvent water molecules are not accessible to the vacant cavity produced by Ala mutation. In addition, ∆∆GAPB as resulting from model A is too far down on the negative scale to be accepted for comparison with the loss of the electrostatic interaction ∆∆Eeletrostatic in the Tyr93Ala mutant, because the PB term is an electrostatic contribution to the solvation free energy and the absolute value of ∆∆GPB is usually comparable in magnitude to ∆∆Eelectrostatic. Thus we can conclude a model B is more appropriate model than a model A for Try93Ala mutant. In this analysis, the new important residue Tyr112 and Thr181 are found for the binding free energy using a model B. Considering the favorable binding free energy and the correlation between rmsd of Thr181 and the temporal development of R(N5-H14) distance, we can conclude that Thr181 is one of the important key residues which helps NADH to approach the isoalloxazine in FAD. In contrast, the correlation between rmsd of Tyr112 and R(N5-H14) are not observed in this study. 3.3. Ab Initio Calculation of the Possible H- Transfer and Electron-Transfer Processes using a Very Simplified Model System. In a previous paper,2 we proposed, on the basis of DFT calculations, a mechanism for the flavin-catalyzed dehydrogenation of amino acids via the protonation at O12 in FAD through a hydrogen bonding between the apoenzyme and FAD (we refer

5724 J. Phys. Chem. B, Vol. 112, No. 18, 2008

Asada et al.

Figure 6. Two perpendicular views of optimized structures with HEME(Fe)+IMI for (1a and 1b) reactant, (2a and 2b) model transition state, and (3a and 3b) product using B3LYP/LanL2DZ. In addition, two perpendicular views of optimized structures without HEME(Fe)+IMI for (4a and 4b) reactant, (5a and 5b) model transition state, and (6a and 6b) product using B3LYP/LanL2DZ are given.

Stabilities and Reactivities of NADH Cytochrome B5 Reductase

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5725

TABLE 4: Calculated Energies for All Model Complexes Related with Our Model Reaction Mechanism molecule reactant

1, 2

transition state

1, 2

product

1, 2

reactant

0, 1

transition state

0, 1

product

0, 1

10-MIA

0, 1

10-MIAH

0, 2

NA+

1, 1

NAH

0, 1

HEME(Fe2+

)+IMI

a

energy (hartree)

∆E (kcal mol-1)

NAH+10-MIA+HEM+IMI ROHF/LANL2DZ B3LYP/LANL2DZ ROHF/LANL2DZ B3LYP/LANL2DZ ROHF/LANL2DZ B3LYP/LANL2DZ

-2650.09768 -2667.05983 -2650.03125 -2667.04188 -2650.13998 -2667.09186

0.0 0.0 41.7 11.3 -26.5 -20.1

NAH+10-MIA RHF/LANL2DZ B3LYP/LANL2DZ RHF/LANL2DZ B3LYP/LANL2DZ RHF/LANL2DZ B3LYP/LANL2DZ

-1321.08209 -1329.37115 -1321.00228 -1329.33786 -1321.08512 -1329.37768

0.0 0.0 50.1 20.9 -1.9 -4.1

Miscellaneous RHF/LANL2DZ B3LYP/LANL2DZ ROHF/LANL2DZ B3LYP/LANL2DZ RHF/LANL2DZ B3LYP/LANL2DZ RHF/LANL2DZ B3LYP/LANL2DZ RHF/LANL2DZ B3LYP/LANL2DZ

-866.55652 -871.94905 -867.14414 -872.56110 -453.74989 -456.59869 -454.51623 -457.41054 -1329.22467 -1337.90958

computational levela

charge, spin multiplicity

0, 1

The procedures used for the geometry optimizations were ROHF/LanL2DZ for all molecules.

SCHEME 1: Calculated Mulliken and Natural (parenthesis) Charges on Each Molecule along the Assumed Reaction Path

to this as the proton transport channel for driving H- transfer from the substrate to FAD). When O12 in FAD is protonated, it becomes a good H- acceptor. This protonation also determines the decreasing of the activation energy of the reaction. In the case of a b5R-b5 complex, it is reasonable to assume that the positively charged HEME(Fe3+) approaches the isoalloxazine ring in FAD and becomes a trigger of H- and/or the electron-transfer reactions. As no experimental structural data of b5R-b5 complexes have been reported so far, we decided to construct some very simplified model systems with 10-MIA,

NAH, HEME(Fe3+) and IMI in order to elucidate the possible H- transfer reaction mechanism from NAH to 10-MIA. In the present study, we have restricted some structural parameters of the relative conformations between NAH and 10-MIA during optimizations, because both of them are rather strongly bounded with the apoenzyme. The numbering of the atoms for 10-MIA and NAH are shown in Figure 5. The bond distances R[C17(NAH)-N10(10-MIA)] are set to 5.6 Å, and the bond angle A[C17(NAH)-N10(10-MIA)-N5(10-MIA)] is set to 80.0 degrees. In addition, in order to construct the transition state model structure, the distances R[N5(10-MIA)-H14(NAH)] and R[C4(NAH)-H14(NAH)] are restricted to 1.24 and 1.38 Å, respectively, whereby those parameters are extracted from the transition state of the NAH/10-MIA complex by B3LYP/6-31G(d,p) calculation, which has one imaginary normal-mode frequency of 2390i cm-1. The optimized structures and the calculated electronic energies are summarized in Figure 6 and Table 4 at the B3LYP/ LANL2DZ//ROHF(or RHF)/LANL2DZ level of theory. The optimized distances between N5 in 10-MIA and Fe in HEME, R(Fe-N5) were calculated to be 4.75, 4.04, and 3.73 Å for complexes 1, 2, and 3, respectively. The relative energy of the product was found to be -20.1 kcal/mol for the complex of NAH‚‚‚10-MIA‚‚‚HEME+IMI (see Table 4), which is lower than that of the reactant, in contrast with the complex of NAH‚ ‚‚10-MIA, for which the relative energy of the product is -4.1 kcal/mol. Table 4 shows that the presence of HEME apparently reduces the activation energy of H- transfer reaction to +11.3 kcal/mol and stabilizes the product. Thus, we propose a favorable reaction mechanism illustrated in Scheme 1 together with the calculated total Mulliken and natural populations on each molecule along the proposed reaction path. As soon as the reactant is prepared, the electron in NAH is immediately transferred to HEME(Fe3+) to produce the reduced form HEME(Fe2+). Three kinds of electronic structures of the model systems are possible. The iron in the HEME could be Fe3+(d5) (S ) 1/2, 3/2, and 5/2) in the initial state, which becomes Fe2+(d6) (S )

5726 J. Phys. Chem. B, Vol. 112, No. 18, 2008

Asada et al. Quantum mechanics/molecular mechanics (QM/MM) calculations of the b5R+b5 docking complex are doing on a parallel computer system, on the basis of the results of the present simplified model system. The results will be reported in our forthcoming paper. 4. Conclusion

Figure 7. Calculated potential energy curves based on the energy of the reactant at R(Fe-N5) ) 5.0 Å. The circle, square, and diamond mean reactant, transition state, and product, respectively.

0, 1, and 2) after the one-electron transfer reaction is completed. On the other hand, the NAH‚‚‚10-MIA could have S ) 0 and S ) 1/2 for the initial and the final state, respectively. According to the spin combinations, three spin states, S ) 1/2 (doublet), 3/2 (quartet), and 5/2 (sextet), are possible. However, Swart et al.33 showed that the lowest energy spin state for Fe(III) with a 6-coordinated system is doublet, which is co-incident with the experimental data. Thus, we have used doublet states for NAH10-MIA complexes in this study. Figure 7 shows calculated potential energy curves along the distance between Fe atom in HEME and N5 atom in 10-MIA, R(Fe-N5). Only R(Fe-N5) has been changed, whereas the relative conformations between NAH/10-MIA and HEME+IMI are restricted to be same structural parameters of each optimized structure. The energy difference between the reactant and the transition state becomes 16.5 kcal/mol at R(Fe-N5) ) 5.0 Å. If the conformational relaxation is taken into account, the activation energies of the model complexes may take an even lower value of 11.3 kcal/ mol as listed in Table 4. Our theoretical calculations indicate that H- transfer from NAH to 10-MIA is an exergonic reaction. In the present calculation, H- transfer is assumed to be induced by the formation of a b5R-b5 complex, in which HEME is in direct contact with 10-MIA. Thus, there is room for further discussion about the possibility of reactions taking place when the HEME is located at a long distance from FAD. The important conclusion from our model calculations is that the electron and/or H- transfers from NADH to FAD are obviously accelerated in b5R+HEME (in b5) complexes. The reaction mechanism proposed by Iyanagi et al. implies that two electrons are transferred from NADH to FAD by H- transfer at the initial stage of the catalytic cycle. Then the two-electron reduced enzyme-NAD+ complex (E-FADH--NAD+) transfers two electrons to two one-electron acceptors one by one, after which the reduced enzyme returns to the oxidized state. In our calculations, their reaction mechanism is rather difficult because the barrier is too high. However, upon docking of the b5, the H- transfer takes place easily (see Table 4). Our reaction mechanism proposes that b5 approaches to b5R at the initial stage of the catalytic cycle through favorable electrostatic interactions. Concerted with the H- transfer is a reduction of the HEME moiety of b5, decreasing the oxidation state from Fe(III) to Fe(II). Thus the first electron has been transferred to the b5. After the first electron transfer is realized, the proton may be released as mentioned by Iyanagi et al.4 Then, the other b5 approaches to b5R, the second electron transfer will proceed and reduced FAD returns to the oxidized state. However, we have not carried out calculations after the release steps of the first incoming b5 yet.

We have carried out MD simulations of rat b5R in explicit water molecules at room temperature to clarify the role of each surrounded residue for the enzymatic catalysis. The alanine scanning method proposed by Massova et al.13 was modified to reproduce the results of the experimental steady-state enzymatic activities kcatNADH/KmNADH. Although the conventional method predicts the favorable mutation for Tyr93Ala mutant of b5R in contrast with the unfavorable experimental measurement, the modified model can reproduce the correct trends of experimental enzymatic activities within an error. Both approaches were applied for 35 surrounding residues within 8 Å from NADH and FAD, and significant decreases in binding free energy are expected when one of the four residues Arg91, Lys110, Ser127, and Thr181 is mutated to Ala. There is a favorable electrostatic interaction between the positively charged residue Arg91 and the dinegative diphosphate moiety in FAD. Lys110 is considered to be one of the most important residues for the stabilization of the complex. Lys110 is the important bridging residue between FAD and NADH, which is located in the FAD domain in b5R. Ser127 and Thr181 have a hydroxyl group that makes favorable hydrogen bonds with the negatively charged diphosphate moiety in FAD and NADH. Tyr112 is estimated to be important using a model B. As there are no experimental structural data on the b5R-b5 complex to date, we have constructed very simplified model systems in order to elucidate the possible H- transfer mechanism from NADH to FAD affected by the addition of HEME(Fe3+). All DFT calculations have been carried out using B3LYP/ LANL2DZ//ROHF(or RHF)/LANL2DZ level of theory. Our calculated results suggest that the electron and/or H- transfer reaction from NADH to FAD can be noticeably accelerated in the presence of HEME(Fe3+). We have proposed that b5 approaches to b5R at the initial stage of the catalytic cycle through favorable electrostatic interactions. Concerted with the H- transfer is a reduction of the HEME moiety of b5, decreasing the oxidation state from Fe(III) to Fe(II). Thus the first electron has been transferred to the b5. Acknowledgment. This work was supported by the Next Generation Super Computer Project, Nanoscience Program, MEXT, Japan. The first author acknowledges with sincere appreciation the financial support provided by the Core Research for Evolutionary Science and Technology (CREST) “High Performance Computing for Multi-scale and Multi-physics Phenomena” from the Japan Science and Technology Agency. A part of this work was supported by the “Nanotechnology Support Project” of the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. Research described in this article was also supported in part by Grants-in-Aid for scientific research for Priority Area “Molecular Theory for Real Systems” from the Japanese Ministry of Education, Culture, Sports, Science and Technology. Supporting Information Available: Figure 1S shows the atomic type of Amber99 force-field parameter sets used in our molecular dynamics simulations for FAD and NADH. All RESP charges and bond angle parameters for FAD and NADH not

Stabilities and Reactivities of NADH Cytochrome B5 Reductase included in the AMBER99 force field parameters are listed in Tables 1S and 2S. The computational alanine-scanning mutagenesis results of b5R for less important residues are given in Table 3S (PDF). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Hildebrandt, A.; Estabrook, R. W. Arch. Biochem. Biophys. 1971, 143, 66. (2) Nishimoto, K.; Higashimura, K.; Asada, T. Theor. Chem. Acc. 1999, 102, 355. (3) Spatz, L.; Strittmatter, P. J. Biol. Chem. 1973, 248, 793. (4) Iyanagi, T.; Watanabe, S.; Anan, K. F. Biochemistry 1984, 23, 1418. (5) Karplus, P. A.; Daniels, M. J.; Herriott, J. R. Science 1991, 251, 60. (6) Correll, C. C.; Batie, C. J.; Ballow, D. P.; Ludwig, M. L. Science 1992, 258, 1604. (7) Nishida, H.; Inaka, K.; Yamanaka, M.; Kaida, S.; Kobayashi, K.; Miki, K. Biochemistry 1995, 34, 2763. (8) Bewley, M. C.; Marohnic, C. C.; Barber, M. J. Biochemistry 2001, 40, 13574. (9) Strittmatter, P. J. Biol. Chem. 1962, 237, 3250. (10) Strittmatter, P. J. Biol. Chem. 1963, 238, 2213. (11) Strittmatter, P. J. Biol. Chem. 1965, 240, 4481. (12) Nicholls, A.; Honig, B. J. Comput. Chem. 1991, 12, 435. (13) Massova, I.; Kollman, P. A. J. Am. Chem. Soc. 1999, 121, 8133. (14) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, W. J., III; Thomas, E.; Ross, W. S.; Simmerling, C.; Darden, T.; Merz, K. M.; Stanton, R. V.; Cheng, A.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. A. Amber-7; University of California at San Francisco: San Francisco, 2002. (15) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (16) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049. (17) Kimura, S.; Nishida, H.; Iyanagi, T. J. Biochem. 2001, 130, 481. (18) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926.

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5727 (19) Guenot, J.; Kollman, P. A. J. Comput. Chem. 1993, 14, 295. (20) Massova, I.; Kollman, P. A. Perpect. Drug DiscoVery Des. 2000, 18, 113. (21) Srinivasan, J.; Miller, J.; Kollman, P. A.; Case, D. A. J. Biomol. Struct. Dyn. 1998, 16, 671. (22) Sitkoff, D.; Sharpe, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978. (23) Kuhn, B.; Kollman, P. A. J. Am. Chem. Soc. 2000, 122, 3909. (24) Asada, T.; Gouda, H.; Kollman, P. A. J. Am. Chem. Soc. 2002, 124, 12535. (25) Sanner, M. F.; Olson, A. J.; Spehner, J. C. Biopolymers 1996, 38, 305-320. (26) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J. Gaussian 03; Gaussian, Inc.: Wallingford, CT, 2004. (27) Novotny, J.; Bruccoleri, R. E.; Davis, M.; Sharp, K. A. J. Mol. Biol. 1997, 268, 401. (28) Misra, V. K.; Sharp, K. A.; Friedman, R. A.; Honig, B. J. Mol. Biol. 1994, 238, 245. (29) Misra, V. K.; Hecht, J. L.; Sharp, K. A.; Friedman, R. A.; Honig, B. J. Mol. Biol. 1994, 238, 264. (30) Sharp, K. A. Biophys. Chem. 1996, 61, 37. (31) Bruccoleri, R. E.; Novotny, J.; Davis, M. E. J. Comput. Chem. 1997, 18, 268. (32) Kimura, S.; Kawamura, M.; Iyanagi, T. J. Biol. Chem. 2003, 278, 3580. (33) Swart, M.; Groenhof, A. R.; Ehlers, A. W.; Lammertsma, K. J. Phys. Chem. A 2004, 108, 5479.