German Cavelier

PROTEINS: Structure, Function, and Genetics 43:420 – 432 (2001) Mechanism of NAD(P)H:Quinone Reductase: Ab Initio Studies of Reduced Flavin German Cavelier1,2 and L. Mario Amzel1* 1 Department of Biophysics and Biophysical Chemistry, Johns Hopkins University School of Medicine, Baltimore, Maryland 2 Program in Molecular Biophysics, Johns Hopkins University, Baltimore, Maryland ABSTRACT NAD(P)H:quinone oxidoreductase type 1 (QR1, NQO1, formerly DT-diaphorase; EC 1.6.99.2) is an FAD-containing enzyme that catalyzes the nicotinamide nucleotide-dependent reduction of quinones, quinoneimines, azo dyes, and nitro groups. Animal cells are protected by QR1 from the toxic and neoplastic effects of quinones and other electrophiles. Alternatively, in tumor cells QR can activate a number of cancer chemotherapeutic agents such as mitomycins and aziridylbenzoquinones. Thus, the same enzyme that protects the organism from the deleterious effects of quinones can activate cytotoxic chemotherapeutic prodrugs and cause cancer cell death. The catalytic mechanism of QR includes an important initial step in which FAD is reduced by NAD(P)H. The unfavorable charge separation that results must be stabilized by the protein. The details of this charge stabilization step are inaccessible to easy experimental verification but can be studied by quantum chemistry methods. Here we report ab initio quantum mechanical calculations in and around the active site of the enzyme that provide information about the fine details of the contribution of the protein to the stabilization of the reduced flavin. The results show that (1) protein interactions provide approximately 2 kcal/ mol to stabilize the planar conformation of the reduced flavin isoalloxazine ring observed in the X-ray structure; (2) the charge separation present in the reduced planar form of the flavin is stabilized by interactions with groups of the protein; (3) even after stabilization, the reduction potential of the cofactor remains more negative than that of the free flavin, making it a better reductant for a larger variety of quinones; and (4) the more negative reduction potential may also result in faster kinetics for the quinone reduction step. Proteins 2001;43:420 – 432. effects of quinones, whose redox cycling is avoided by the obligatory two-electron reduction carried out by QR1. Surprisingly, the same reaction that protects against the toxic effects of quinones can be lethal for tumor cells by activating cytotoxic quinone anticancer drugs such as mitomycins and aziridylbenzoquinones. Because levels of QR1 are elevated in some tumors, the enzyme provides an opportunity to develop improved chemotherapeutic agents.2–7 QR1 is a homodimer (273 residues per monomer) with two identical catalytic sites at two equivalent locations on the dimer interface (Fig. 1). Each binding site, formed by residues of both monomers, contains an FAD prosthetic group that remains noncovalently bound during catalysis. NAD(P)H and NAD(P)⫹ cycle in and out of the catalytic site of the enzyme, providing reducing equivalents for the reaction. The kinetics of the reaction is ping-pong: the oxidized cofactor NAD(P)⫹ must be released before the quinone substrate can bind. The structure of rat QR1 (rQR1), as determined by X-ray diffraction methods, shows that each monomer is formed by two domains: an N-terminal catalytic domain and a C-terminal domain involved in dimerization.5 The major portion of each monomer, the catalytic domain, has an overall fold characteristic of other flavoproteins: a twisted central parallel ␤-sheet surrounded on both sides by connecting helices. The structures of rQR1 determined in complexes with NAD(P)⫹ and with a substrate and an inhibitor have provided a detailed description of the catalytic site and important insights into the mechanism of the enzyme.5 Recently, the structures of apo human and mouse QR1 have been reported8 as well as the structure of another cytosolic quinone reductase, QR2.9 These structures show that the isoalloxazine moiety of FAD interacts with resi- © 2001 Wiley-Liss, Inc. Key words: density functional; enzymes; flavoproteins; DT-diaphorase; isoalloxazine; quantum chemistry INTRODUCTION NAD(P)H:quinone oxidoreductase type 1 (QR1, EC 1.6.99.2; for a review see Ref. 1) is an FAD-containing protein that catalyzes the nicotinamide nucleotide-dependent reduction of quinones, quinoneimines, azo dyes, and nitro groups.1–5 It protects against the toxic and neoplastic © 2001 WILEY-LISS, INC. Abbreviations: NAD(P)H, nicotinamide adenine dinucleotide (phosphate), reduced; NAD(P)⫹, nicotinamide adenine dinucleotide (phosphate), oxidized; FAD, flavin adenine dinucleotide, also oxidized form; FADH2, flavin adenine dinucleotide, reduced; SCF, self-consistent field method. Grant sponsor: National Institute of General Medical Sciences; Grant numbers: GM45540, 5T32GM-. *Correspondence to: L. Mario Amzel, Department of Biophysics and Biophysical Chemistry, Johns Hopkins University School of Medicine, 725 North Wolfe Street, Baltimore, MD-. E-mail:-Received 22 May 2000; Accepted 1 February 2001 Published online 00 Month 2001 AB INITIO STUDY OF NAD(P)H:QUINONE REDUCTASE Fig. 1. Ribbon diagram of the X-ray structure of rQR1. The active sites between the two domains of the homodimer show the noncovalently bound flavin cofactor. dues from loops of both monomers. Aromatic residues Tyr-104, Trp-105, and Phe-106 and the main chain of Leu-103 interact directly with the rings of the isoalloxazine, anchoring it. The two oxygen atoms of the flavin ring (O2F1 and O4F) form hydrogen bonds with mainchain NH groups of the protein: O4 with Phe-106 and O2 with Gly-150. Nitrogens in the rings also form hydrogen bonds with NH groups of the protein: N1F with Gly-149 and N5F with Trp-105. The two methyl groups of the isoalloxazine ring are seated in a pocket formed by two residues of one monomer (Tyr-104 and Trp-105) and four residues of the other monomer (Ile-50, Tyr-67, Pro-68, and the main chain of Glu-117). The ribitol, phosphates, and adenine ring of FAD interact noncovalently with several loops and helices, anchoring the FAD cofactor to the enzyme.5 The mechanism proposed for QR1 has two distinct steps: (1) hydride transfer from NAD(P)H to the enzyme-bound flavin followed by the release of NAD(P)⫹ and (2) hydride transfer from the flavin to the quinone followed by the release of the hydroquinone. In the first half of the reaction cycle, the nicotinamide of NAD(P)H is in an ideal position for a direct hydride transfer to FAD: C4 of the nicotinamide is 4.0 Å from N5 of flavin.5 As mentioned previously, the N1F is hydrogen-bonded to the NH of Gly-149. If a negative charge develops on N1F upon FAD reduction, there are no groups that can donate a proton to compensate the charge, so the most likely tautomer of the FAD is the enolic (FADH2) form, with part of the negative charge at O2F. This oxygen is already an acceptor of a hydrogen 1 The usual atom nomenclature for each cofactor is followed with a final capital letter indicating the cofactor: F, flavin; A, adenosine of FAD; and N, nicotinamide-ribose. 421 bond from Tyr-155, which can in turn be stabilized by the positive charge of (or the transfer of a proton from) His-161. This process is reversed5 when the hydride is transferred, following the opposite direction, to the quinone that replaced the leaving NAD(P)⫹. The structures provide a simple rationale for the observed ping-pong kinetics of the reaction; the binding of the substrate cannot occur until NAD(P)⫹ is released because the quinone and nicotinamide share the same site. In this scheme, after NAD(P)⫹ leaves, the quinone arrives and binds in an ideal orientation to receive the hydride coming back from FADH2. The quinone is reduced by the hydride to the singly ionized hydroquinone (hydroquinolate), and the isoalloxazine is oxidized to the quinonoid form. The proton on O2F is transferred back to the OH of Tyr-155. The imidazole ring of His-161 becomes fully protonated again and can either transfer a proton to the hydroquinolate or simply stabilize its negative charge. Therefore, in addition to hydride transfer, the net effect of the second half of the reaction is to transfer a proton from O2F (accepted in the first half of the reaction) to the hydroquinolate.5 This mechanism explains the obligatory two-electron reduction promoted by QR1: both halves of the reaction involve a hydride transfer (not treated in this article), first from NAD(P)H to FAD and then from FADH2 to the quinone. The charge relay formed by Tyr-155 and His-161 allows the reaction to proceed without unfavorable charge separations.5 The catalytic mechanism of FAD reduction and subsequent oxidation in its active pocket involve a fine-tuned succession of steps. Charge separation and subsequent charge stabilization by the enzyme are at the core of this mechanism. For a mechanism of this kind, theoretical quantum chemical studies of the electronic configuration of the free FAD and the FAD in the binding environment of the catalytic site can provide direct information about the distribution of the charges during the reaction. In this article, we report results of quantum mechanical ab initio calculations aimed at probing the QR1 catalytic mechanism, particularly the electrostatic stabilization of the FAD after reduction, and the energetic contributions of the cofactor and the protein to the completed process. Hydride transfers are not treated in this article; only the states before and after the reduction of the FAD are treated. Reaction mechanisms in enzymes have generally been studied and rationalized on the basis of structural and experimental data.4,10 –12 Significant advances have been made in developing methods that can give a quantum mechanical description of the electronic structure of entire macromolecules.13–16 An important goal that has yet to be accomplished is the detailed description of enzymatic reaction mechanisms from first principles. Enzymatic mechanisms involving electron transfer, hydride transfer, and proton tunneling might be explained in greater detail on the basis of precise electronic structures. For example, recent experimental findings17–19 point to the possible involvement of hydrogen tunneling in enzyme catalysis, a process that would be better understood with the knowl- 422 G. CAVELIER AND L.M. AMZEL edge of the electronic structure of reactants, products, and intermediates calculated quantum chemically in the vicinity of the enzyme active site. Quantum chemical calculations on coenzymes and other enzyme ligands, in isolation or interacting with a limited number of atoms, have provided additional insight into their function and properties during enzyme catalysis.20 –24 Recently, quantum chemistry has been used to elucidate the precise events occurring at the active site of ribonuclease A.25 Calculations based on the extended Hückel method have been reported for the enzyme DNA photolyase.26 Combined quantum mechanics/molecular mechanics is a well-established method in which the quantum mechanics is applied only to a restricted part of the molecule and the rest studied with molecular mechanics.27–29 Other approaches involve model systems to study with ab initio methods the transition-state structure.30 Other recent methods include the use of density functional theory (DFT) with continuum dielectric theory31 and the application of DFT to optimize geometries of the transition state in a model system.32 Flavin-dependent reductases and oxidases comprise an important family of enzymes.10,33 Several mechanisms have been proposed and studied with different experimental and theoretical methods in reactions involving this cofactor.4,34,35 In addition, quantum mechanical calculations have been used to explore the electronic structure of the isolated isoalloxazine ring in the reduced and oxidized states.36 –39 However, no calculations have been reported for flavins in a protein environment. QR1 is an excellent system for this type of calculation. First, many of the mechanistic questions involve understanding the distribution of charges and their stabilization, exactly the kind of questions that quantum mechanical calculations can answer. Second, the reaction takes place in a cavity within the protein that is not directly connected to a bulk solvent. Thus, the calculations do not require the inclusion of complex solvent models. In general, it can be said that some redox enzymes, those having properties similar to QR1, are amenable to studies such as the one presented here. MATERIALS AND METHODS Programs All quantum mechanical calculations were carried out with the program GAUSSIAN 9840 (Gaussian Inc.) run on J90 and SV1 CRAY computers at the Advanced Biomedical Computing Center of the National Cancer Institute (Frederick, MD). The three-dimensional visualizations were performed with the graphic capabilities of the software UNICHEM (Oxford Molecular, Inc.) implemented on SGI workstations. The initial structure of the cofactor was built on an SGI workstation with the program QUANTA (Molecular Simulations, Inc.). The method of convergence used in all SCF calculations was DIIS (Pulay’s Direct Inversion in the Iterative Subspace extrapolation method) with a tight (10⫺8) convergence limit. Additionally, all geometry optimizations used a tight (10⫺8) convergence limit. Isolated Isoalloxazine Ring The isoalloxazine ring was used with a methyl termination at N10F instead of the ribitol-adenine portion of riboflavin. The initial structure constructed with QUANTA was optimized with GAUSSIAN 98. Initial optimization with the MNDO semiempirical method was followed by a full refinement optimization at the RHF/6-31G level.40 The optimized structures were used for all further calculations. Initial exploratory calculations were performed with the RHF/STO-3G model chemistry. To test the importance of polarization effects in the calculations, we repeated several RHF calculations with the 6-31G(d) basis set. Calculations with this basis set were also performed with the hybrid DFT method B3LYP (Becke’s41 style three-parameter DFT with LYP gradient-corrected correlation functional of Lee et al.42) to include some of the effects of electronic correlation. Based on the proposed mechanism, calculations were carried out for the following configurations of the isoalloxazine ring (Fig. 2): (A) oxidized form, (B) reduced form I (anion, charge on N1F, O2F in keto form), (C) reduced form II (anion, charge on O2F enolate), (D) resonance form of the configurations in A and B, (E) reduced form III (neutral, O2F as OH), and (F) reduced form IV (neutral, O2F in keto form, NH at N1F). The RHF/STO-3G calculations used 309 –318 basis functions. The central processing unit (CPU) time employed for each calculation was typically 15–22 h. The B3LYP/6-31G(d) calculations used around 310 basis functions and CPU times of 6 –7 h. The numerical values computed comprised the energy of each configuration, the Mulliken charge of each atom, and the three-dimensional representations of electrostatic potential, the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). Isoalloxazine Ring in the Protein Environment We performed equivalent calculations with the isoalloxazine ring anchored in the environment of the catalytic site of QR1 by including a number of selected residues around the ring with the coordinates of the QR1 X-ray structure. The electrostatic potential and charges around the flavin molecule, as well as the hydrogen bonds the flavin makes with active site residues, are thought to have a predominant role in the catalytic mechanism.5,20,23,24,34,43,44 On the basis of the analysis of the reaction mechanism previously summarized,5 several amino acids in the environment of the flavin cofactor were selected for the calculations because (1) they form the hydrophobic binding site, (2) they make hydrogen bonds with the isoalloxazine ring, or, more importantly, (3) they are part of the proposed charge-relay system. On the basis of these criteria, the following amino acids were selected: Trp-105, Phe-106, Gly-149, Gly-150, Tyr-155, and His-161. All individual amino acids and amino acid pairs (Gly-149 & Gly-150, Trp-105 & Phe-106) were terminated at their carboxy terminus with a methyl group. In the calculations, the X-ray coordinates of the flavin were replaced in turn by the coordinates of each of the optimized isoalloxazine rings described previously, main- AB INITIO STUDY OF NAD(P)H:QUINONE REDUCTASE 423 Fig. 2. Forms of the isoalloxazine ring used in this study: (A) oxidized form before reduction, which represents the state of the FAD moiety at the beginning of the catalytic cycle immediately after the release of the reduced quinone; (B) reduced form I, with a negative charge on N1F and with O2F in the keto form, which depicts the standard way of representing the reduced FAD; (C) reduced form II, with a negative charge initially located at the O2F atom, which corresponds to the ionized enolate form of the isoalloxazine ring; (D) resonance form of the configurations in parts B and C (only this configuration exists for this combination of charges and atoms); (E) reduced form III, in which the charge of the reduced FAD was neutralized by the addition of an H⫹ to O2F, which is, therefore, protonated in the enol form; and (F) reduced form IV, with O2F in the keto form and the proton on N1F. taining the orientation and position of the ring in the experimental X-ray rQR1 structure. The three structures shown in Figure 3(A–C) for the isolated rings were used in this portion of the work, with proper allowance for the variations induced by the protein. Because of the size of the analyzed structures, exploratory calculations were carried out with a minimal basis set with restricted Hartree–Fock theory (RHF/STO-3G; ⬃150 atoms, ⬃800 basis functions). Final calculations for selected structures were carried out with the 6-31G(d) basis set with DFT [B3LYP/6-31G(d)]. These calculations used 1380 basis functions and took 40 h to 14 days of CPU time. Effects of Geometry and Dielectric To analyze the effect of the protein on the geometry of the reduced isoalloxazine, we performed single-point energy calculations for bent configurations of the cofactor, isolated and anchored in the catalytic site of QR at several dihedral angles. The coordinates used in these calculations were generated first with the bent optimized geometry found for the reduced form of the free isoalloxazine ring. The different conformations were generated with the program QUANTA by changes in the dihedral angle in 5° intervals around the virtual bond between N5F and N10F. The values for these intervals were begun at 145° (a more bent configuration than the equilibrium configuration of the isolated ring, 151.3°) and increased past the planar configuration (180°) to 210°. To release the strain in the isoalloxazine ring created by the changes, we optimized the coordinates of the substituents on N10F and N5F at each new dihedral angle. These restrained optimizations of the rings, isolated and in the protein environment, used the semiempirical MNDO 94 program with the AM1 Hamiltonian as implemented in UNICHEM. Energies for 424 G. CAVELIER AND L.M. AMZEL Fig. 3. Initial exploratory electronic structure of the RHF/STO-3G results. The simplified, compressed representations use the following color scheme: electrostatic potentials are represented in red (negative) and blue (positive), LUMOs are represented by solid green surfaces, and HOMOs are represented by solid gold. Only one phase of the molecular orbitals is shown for clarity; the lobes with the phase of the opposite sign are located symmetrically around the corresponding atoms. The Mulliken charges of selected atoms are indicated as fractions of the electron charge. The free isoalloxazine ring is shown here in three of the configurations of Figure 2: (A) oxidized form of FAD before reduction, (B) initial effect of the reduction by the hydride (note the replacement of LUMOs by HOMOs and the appearance of an unfavorable charge separation, very negative around N1F and O2F), and (C) enolic form of the isoalloxazine ring upon protonation of O2F (note the neutralization of the previous charge separation and the symmetrical redistribution of HOMOs and LUMOs). these optimized structures were calculated with the RHF/ STO-3G model chemistry. This level of theory was considered adequate for these calculations because the resulting Fig. 4. DFT B3LYP/6-31G(d) electronic structure results. The three free isoalloxazine rings shown correspond to the configurations of Figure 3(A–C). energies were used only to estimate differences with respect to the energy at the minimum. To partially account for the effects of the solvent, we repeated calculations for the isolated FAD with the self-consistent-reaction-field (SCRF) method45 (as implemented in GAUSSIAN46), with a relative dielectric constant of 80. Another way of partially accounting for the effects of the solvent is to satisfy the hydrogen-bonding potential of N1F, O2F, N3F, O4F, and N5F with individual water molecules. To this effect, water molecules were placed in hydrogen-bonding positions around the isoalloxazine ring of the structure with a proton in N1F [Fig. 2(F)]. The geometry of this hydrogenbonded isoalloxazine ring was then optimized at the AB INITIO STUDY OF NAD(P)H:QUINONE REDUCTASE RHF/STO-3G level of theory, and electronic structure calculations were carried out with the 6-31G(d) basis set as described and shown later in Figure 7(A–B). RESULTS AND DISCUSSION Isolated Isoalloxazine Ring Quantum mechanical calculations for the chosen configurations of the isoalloxazine ring [Fig. 2(A–F)] were carried out as described previously. The oxidized form [Fig. 2(A)] represents the state of the FAD moiety at the beginning of the catalytic cycle, immediately after the release of the reduced quinone. Structural optimization shows that in this form the isoalloxazine ring adopts a planar conformation. A comparison of the calculated optimized structure with the experimental structure from the Cambridge Crystallographic Database shows an excellent agreement (root-mean-square deviation ⫽ 0.08 Å). Reduced form I [Fig. 2(B)], with a negative charge on N1F and with O2F in the keto form, depicts the standard way of representing the reduced FAD. Reduced form II [Fig. 2(C)] corresponds to the form with the O2F in the enolate form (with a negative charge). Upon optimization, however, the two configurations give almost exactly the same coordinates (root-mean-square deviation ⫽ 0.0009 Å) and atomic Mulliken charges, as expected for true resonance forms. Therefore, only one configuration exists for such a combination of charges and atoms [Fig. 2(D]. The optimized structure shows that the isoalloxazine adopts a bent, butterfly-like conformation with a dihedral angle of 151.3° between the two halves of the molecule. Two forms were calculated in which the charge of the reduced FAD was neutralized by the addition of an H⫹: reduced form III [Fig. 2(E)] with O2F protonated in the enol form and reduced form IV [Fig. 2(F)] with O2F in the keto form and the proton on N1F. The dihedral angles of these conformations are 155.3 and 148.4° respectively. This conformation is in excellent agreement with the experimentally determined structure of lumiflavin47 (dihedral angle ⫽ 144.5°). A previous quantum chemical study39 also found a bent conformation, but the bend angle (164°) was significantly farther from the experimental value. Initial exploratory RHF/STO-3G results of the quantum mechanical calculations performed after geometry optimization are shown in Figure 3(A–C). The simplified, compressed representation is described in the figure. The Mulliken charges of selected atoms are indicated as fractions of the electron charge. Only three of the configurations in Figure 2(A–F) are shown [Fig. 3(A–C)]: (1) the optimized oxidized FAD, corresponding to the form in Figure 2(A); (2) the optimized reduced FADH⫺, corresponding to the resonance form in Figure 2(D) that arises from the optimizations of the forms in Figure 2(B,C); and (3) the optimized reduced FADH2, corresponding to the form in Figure 2(E). The optimized reduced form IV [Fig. 2(F)] is not shown, as it was not used in the comparisons done in the protein environment because it had a higher energy and a less favorable charge separation than the one in Figure 3(C) (discussed later). 425 The initial exploratory RHF/STO-3G calculations for the oxidized isoalloxazine ring [Fig. 3(A)] show a positive electrostatic potential almost everywhere in the molecule, except for slightly electronegative regions around the O2F and O4F oxygens and to a lesser extent around the N1F and N5F nitrogens, all of which have partial negative charges. The most negatively charged atoms are N1F (Mulliken charge ⫽ ⫺0.333) and N3F (charge ⫽ ⫺0.388; it is compensated by a covalently bound hydrogen with a charge of ⫹0.220). The oxygens are charged more or less equally, OO2F with ⫺0.263 and O4F with ⫺0.243. The electrostatic potential is nonetheless more negative around the oxygens than around the nitrogens because the oxygens, in addition to being more electronegative, are bound to only one atom. By contrast, the nitrogens are more closely surrounded by carbon and hydrogen atoms. The dimethylbenzene ring has an almost uniform charge distribution when hydrogens are taken into account. The unoccupied electron-acceptor orbital (LUMO) in the hydride acceptor N5F has high overlap with the LUMO in C4F, C4aF, and C10aF, which in turn overlaps the LUMO in N1F, C2F, and N3F [Fig. 3(A)], providing a molecular orbital path for the delocalization of the electrons after reduction. Therefore, after reduction of the isolated isoalloxazine ring to the reduced anionic form [Fig. 3(B)], the most important changes in atomic Mulliken charges occur in N1F (from ⫺0.333 to ⫺0.423), O2F (from ⫺0.263 to ⫺0.383), and O4F (from ⫺0.243 to ⫺0.392). The rest of the atoms of the isoalloxazine ring experience smaller changes in charge. Accordingly, the electrostatic potential becomes much more electronegative in the vicinity of N1F, O2F, and O4F [Fig. 3(B)], particularly in the region around O2F. In the absence of a neutralizing proton, the negatively charged N1F shares its negative electrostatic potential with O2F [Fig. 3(B)]. N5F is also more electronegative than before reduction, but the hydride hydrogen neutralizes its negative potential. The HOMOs in the reduced isoalloxazine ring [Fig. 3(B)] are in the regions where the LUMOs were in the oxidized ring [Fig. 3(A)]: the most significant HOMOs occur on N5F, C4F, C4aF, and C10aF; on N1F, C2F and N3F; and on O2F and O4F. The most important difference, however, is that after reduction the planar conformation of the oxidized isoalloxazine ring [Fig. 3(A)] changes to a bent conformation (dihedral angle ⫽ 151.3°) when it becomes reduced [because of the orientation of the drawing, the bent structure is not very apparent in Fig. 3(B)]. In the mechanism proposed for QR1, the reduced isoalloxazine ring is in the enolic form protonated at O2F. For this protonated form, the isolated isoalloxazine ring has the charges redistributed in a symmetrical and nonlocalized manner [Fig. 3(C)]. The negative electrostatic potential produced by charge separation at N1F and O2F in the unprotonated form is almost absent in the uncharged form. Although the two electrons of the reduction are still there, their distribution is less polarized. The overall effect of the protonation of O2F is, therefore, a smoothing of the charge separation and a decrease in the charge asymmetry caused by the reduction. Protonation also has a major 426 G. CAVELIER AND L.M. AMZEL effect on the molecular orbitals: the HOMOs are redistributed toward the central ring in a symmetrical fashion, whereas the LUMOs are redistributed more evenly along the whole isoalloxazine ring. Finally, another possible structure is one protonated at N1F [Fig. 2(F)]. The effect of this protonation (not shown in Figs. 3– 4) is similar to that of protonation at O2F, but the potential at O2F and in O4F is more negative, and the charges remain somewhat separated and more polarized. That is, the charge imbalance is not as uniformly distributed as for the protonation of O2F. There is also some strain in the molecule because the methyl in the N10F bends toward the back of the isoalloxazine ring. For these reasons and because N1F cannot be protonated when bound to rQR1, this form was not used in the comparisons with the ring in the protein environment. Calculations with DFT methods and a polarized basis set that includes d functions for the heavy atoms [B3LYP/ 6-31G(d); Fig. 4(A–C), which parallels Fig. 3(A–C)] give highly similar results. The main difference of including d functions is evidenced as expected in the Mulliken charges that are all larger in the calculation with B3LYP/6-31G(d). Although this is a more energetically accurate calculation, the Mulliken charges calculated this way are probably unrealistic because the higher delocalization afforded by the use of a polarized basis set makes it less appropriate to use the resulting electron density to estimate point charges. Isoalloxazine Ring in the Protein Environment The structure of the isoalloxazine ring refined as part of the crystal structure of QR1 has slight differences in bond lengths and bond angles compared with the structure of the optimized free flavins calculated previously [Fig. 2(A– F)]. Small errors in the experimental coordinates of the QR1 structure, determined at a 2.1-Å resolution, can account for these differences. Therefore, for these calculations hybrid coordinate sets were used: the coordinates of the protein atoms were directly those of the X-ray structure, and the coordinates of the flavin ring were replaced in turn by those of each of the optimized isoalloxazine rings described previously. The orientation and position of the ring were adjusted to be those in the X-ray structure. From energetic analysis (shown later in Fig. 8) and restrained optimizations of the reduced ring in the protein environment, it is clear that the reduced ring inside the protein is more planar than the free reduced ring, adopting a dihedral of around 171–175° in the protein environment instead of the dihedral of 151.3° found for the free ring. The planar form (dihedral angle ⫽ 171–175°) was then used for the final calculations for the reduced isoalloxazine ring. Calculations carried out with the optimized structures of the isoalloxazine ring anchored in the protein environment gave the results shown in Figure 5(A–C). They correspond to the isolated ring structures in Figure 3(A– C). Exploratory RHF/STO-3G calculations show that charges in the isoalloxazine ring have tendencies similar to those found in the isolated ring, but the influence of the protein environment on the values of the charges and the electrostatic potential is significant. Placing the oxidized isoalloxazine ring in the protein environment results in small changes in the Mulliken charges [Fig. 5(A)]. The general trend is toward a higher polarization of the charges [compared with the isolated ring in Fig. 3(A)]: the negative charges of the oxygens and nitrogens increase, and the positive charges of the carbon atoms adjacent to them also increase. This effect helps to anchor the isoalloxazine ring in its place by intensifying the strength of the corresponding hydrogen bonds. The electrostatic potential at and around the atoms of the oxidized isoalloxazine ring is only slightly negative with respect to the surrounding protein. The LUMOs are very similar to those observed in the free isoalloxazine. There is an unoccupied overlapping path from the LUMO on N5F to the LUMO on C4F, C4aF, and C10aF. This second LUMO in turn overlaps the LUMO on N1F, which approaches that of O2F. Therefore, the same unoccupied path is present for the arriving reducing electrons as for the isolated oxidized isoalloxazine ring, from N5F through C4aF and C10aF to N1F and from there to O2F. The influence of the surrounding amino acids on the flavin leaves the oxidized isoalloxazine ring in the protein environment mainly with acceptor (LUMO) orbitals, an optimal situation for receiving the reducing electrons through hydride transfer. The configuration of the anionic form of the reduced isoalloxazine ring in the protein environment is shown in Figure 5(B). The geometry of the free, isolated reduced isoalloxazine ring is the bent conformation described previously (dihedral angle ⫽ 151.3°). However, all X-ray structures of rQR1, as well as most flavoproteins, show a rather planar isoalloxazine ring in the reduced form.48 This conformation must be a consequence of the protein environment, which imposes a planar structure on the reduced isoalloxazine ring. Quantum mechanical calculations confirm this explanation. A restrained optimization (with MNDO 94, AM1 Hamiltonian; see the Materials and Methods section) of the dihedral of this form of the reduced isoalloxazine ring in the protein environment produced a configuration with a dihedral of 171.8°. In addition, an analysis of the bend angle based on single-point energy calculations (shown later in Fig. 8) also shows that a conformation with a dihedral of 171.1° is the most energetically favorable for this form of the reduced flavin in the protein environment, partly because it leads to more favorable hydrogen bonds. Probably, the bent reduced conformation (dihedral angle ⫽ 151.3°) of the free isoalloxazine ring never occurs in the protein environment. The protein environment seems to be ideally suited to adjust the charges, whereas the geometry and the electrostatic potential seem to stabilize a more planar conformation of the isoalloxazine ring.48 Figure 5(B) shows that upon reduction the negative charges of the N1F, O2F, and O4F atoms is increased by the charge of the two reducing electrons migrating through the unoccupied path present in the oxidized FAD [LUMOs in Fig. 5(A)]. The occupied molecular orbitals in the reduced form [HOMOs in Fig. 5(B)] however, have a similar configuration as in the isolated reduced isoallox- AB INITIO STUDY OF NAD(P)H:QUINONE REDUCTASE 427 Fig. 5. Initial exploratory RHF/STO-3G results for the electronic configurations of the isoalloxazine ring in the protein environment (cf. Fig. 3): (A) oxidized isoalloxazine ring in the enzyme catalytic site; (B) reduced isoalloxazine ring in the enzyme catalytic site, the initial step of the charge relay [cf. Fig. 3(B)]; (Bⴕ) reduced isoalloxazine ring in the enzyme catalytic site, the intermediate step in the charge-relay mechanism [cf. Fig. 3(B,C)], with a proton transferred from Tyr-155 to O2F; and (C) reduced isoalloxazine ring in the enzyme catalytic site, the completion of the charge relay [cf. parts B and B⬘ and Fig. 3(C)]. azine ring, with occupancy around the whole molecule. However, the negative charges and potential show a tendency to increase with respect to the oxidized form, mainly in N1F, O2F, and O4F [Fig. 5(B)]. The large charge separation that would exist if this configuration remained as shown (i.e., no charge relay) can be seen in a comparison of the charges and the negative electrostatic potential of this form with those of the reduced ring in Figure 5(C). At this stage, the result of the reduction is the migration of the charge of the reducing electrons from N5F to N1F, 428 G. CAVELIER AND L.M. AMZEL Fig. 6. DFT B3LYP/6-31G(d) electronic structure results for the isoalloxazine ring in the protein environment. The three configurations shown correspond to those of Figure 5(A–C). O2F, and O4F. Despite the ring planarity imposed by the protein, there is still a charge separation, but it is only temporary, and further compensating mechanisms are put in place by the protein environment. Note in Figure 5(B) that a LUMO is present at residue His-161, a residue postulated to cooperate with Tyr-155 in the establishment of a charge-relay system that helps to neutralize the initial charge separation.5 The charge stabilization is brought about by the migration of the proton of the OH of Tyr-155 to O2F; stabilization of the negative charge left in Tyr-155 is accomplished by either proton migration from or direct charge stabilization by His-161. AB INITIO STUDY OF NAD(P)H:QUINONE REDUCTASE Fig. 7. Effects of geometry and dielectric: (A) hydrogen-bonded waters partially account for the effect of the solvent around the isoalloxazine ring and (B) the dielectric effect is taken into account with the SCRF method. Compare the effect of the waters or the SCRF with the effect of the protein in Figure 5(C). An intermediate step in the charge-relay mechanism can be observed in Figure 5(B⬘). This result correspond to the case in which the proton from Tyr-155 has been transferred to the O2F in the ring, but the proton of His-161 has not been transferred to the Tyr-155 phenolic oxygen, which remains negatively charged. Note in Figure 5(B⬘) that there is an additional charge in O2F, N3F, O4F, and related carbons and hydrogens, with respect to the full charge relay shown in Figure 5(C). However, this intermediate step shows that the sole transfer of a proton from Tyr-155 to O2F is enough to cause a neutralization of the initial charge separation in the ring [see the corresponding charges in Fig. 5(B)]. Moreover, occupied HOMOs migrate from the isoalloxazine ring to Tyr-155, with a corresponding negative potential around the side chain of this residue: its oxygen has a charge of ⫺0.531. As a result of the charge relay [Fig. 5(C)], a proton is transferred from Tyr-155 to O2F and stabilizes the enol form. In turn, His-161 stabilizes the charge in Tyr-155. Partial charges of the isoalloxazine atoms [Fig. 5(C)] are changed with respect to those of the oxidized, planar form [Fig. 5(A)] and much more changed with respect to those of the configuration obtained before the charge relay [Fig. 5(B)]. The difference in the calculated energies between the structure before [Fig. 5(B)] and after [Fig. 5(C)] the charge relay is ⫺9.9 kcal/mol, indicating that the charge 429 relay has a large stabilizing effect on the reduced state. A comparison of Figure 5(C) with the case of the free reduced isoalloxazine ring [Fig. 3(C)] shows greater negative charges in N1F, O2F, N3F, and O4F that may result in strong hydrogen bonds. Similar charges to the ones found in the isolated ring of Figure 3(C) are found when the electronic structure is calculated for the unstable, planar configuration of Figure 5(C) taken outside of the protein (not shown). Therefore, the increased charges of the bound cofactor probably reflect the planar configuration of the isoalloxazine ring stabilized by the binding energy provided by the protein [Fig. 5(C), dihedral angle ⫽ 175°]. This planar configuration may in turn contribute to the enhanced reducing power of the FADH2 in the protein and the rate of the second step of the QR1 reaction (FADH2 reducing the quinones; see the Stabilization Energy section later). Comparing Figure 5(A,B) with Figure 5(C) suggests that the energetically favorable interactions result from the electrostatic potential around atoms N5F, N1F, O2F, and O4F. As a result of the charge-relay mechanism, a corresponding equalization (spreading) of the HOMO in the isoalloxazine ring is observed. The negative charges of N1F, O2F, and O4F are lower in magnitude than before the changes associated with the charge relay [Fig. 5(B)], and the adjoining carbon atoms are more positive, such that the cofactor shows mainly a weakly positive electrostatic potential relative to the rest of the environment [Fig. 5(C)]. The HOMOs are distributed evenly throughout the isoalloxazine; in particular, there is a strong HOMO at the N5F atom overlapping with the other HOMOs of the ring. These HOMOs hold the electronic charge for the quinone reduction that can occur after NAD(P)⫹ leaves and the quinone arrives in its place. The LUMOs are no longer at His-161; they are now in the isoalloxazine ring. As in the free isoalloxazine ring, calculations at the B3LYP/6-31G(d) level of theory give very similar results [Fig. 6(A–C)] to those in Figure 5(A–C). Like the other cases discussed previously (Figs. 3 and 4), the main difference is evidenced in the Mulliken charges. Effects of Geometry and Dielectric To study the energetic influence of the protein environment on the conformation of the isoalloxazine ring, singlepoint energy calculations were made for different bent configurations of the ring at several dihedral angles. They were performed at the RHF/STO-3G level for the isolated ring and for the ring anchored in the catalytic site of the enzyme. Although this is a minimal basis set and thus the computed energies may have a bias, it has been found that energy differences from RHF/STO-3G calculations may have the required accuracy.46 Therefore, in this analysis only differences in energy with respect to a minimum were used. To represent more closely the aqueous environment for the free isoalloxazine ring (not in the protein environment), we performed calculations with an approximate method that employs a bulk dielectric constant to represent the solvent: the SCRF method.46 A relative dielectric constant of 80 was used in calculations with the SCRF 430 G. CAVELIER AND L.M. AMZEL Fig. 8. Stabilization energy. The bent configuration of the free isoalloxazine ring (circles) is destabilized in the protein environment, and the planar configuration of the reduced isoalloxazine ring (black squares) is stabilized by the protein environment. Finally, the more planar configuration of the protonated enol form in the catalytic site (black triangles) is stabilized by the protein. method. The relaxed optimized configurations of the rings found before were the basis for these single-point energy calculations at the RHF/STO-3G level with the SCRF method. These results are shown in Figure 7(A–B) and are further discussed later. Calculations with hydrogen-bonded waters, as described in the Materials and Methods section, were also performed as an alternative method to partially account for the effect of the solvent around the free isoalloxazine ring outside the protein. The results with hydrogen-bonded waters [corresponding to the reduced structure of Figure 2(F) with a proton in N1F, the most likely protonated solution form] are shown in Figure 7(A) and can be compared to the visualizations obtained with SCRF electronic structure calculations on the isolated reduced isoalloxazine ring [Fig. 7(B)]. The hydrogen-bonded waters provide an environment similar to the hydrogen bonds of the protein environment: the charges, molecular orbitals, and electrostatic potential are similar to those of Figure 5(C), the result of the charge relay inside of the protein, where the protonation is only possible at O2F, not at N1F. Some charges and negative potentials are much higher than inside the protein because of the negative charge of the oxygens in the waters [Fig. 7(A)]. In general, the trend from in vacuo (not shown) to the SCRF calculation [Fig. 7(B)] to the water hydrogen-bonded structure [Fig. 7(A)] is toward an increase in the charges of negative atoms such as N1F, N5F, O2F, and O4F, together with an increase in the positive charge of positive atoms surrounding them, such as the hydrogen in N1F, the carbon in C2F, and the hydrogen at N5F. The molecular orbitals in the three cases are very similar to the one in Figure 5(C) (protein environment, protonated at O2F). The calculations performed with surrounding waters and those performed with SCRF show similar trends and magnitudes of the changes in the charges. Stabilization Energy The energies associated with a number of conformations of the reduced isoalloxazine ring were calculated for the free ring and the ring in the protein environment before [Figs. 3(B)–5(B)] and after the charge relay [Figs. 3(C)– 5(C)]. The total energy was calculated for different degrees of bending (the angle between the planes of the two outer isoalloxazine rings) to estimate the energetic contribution that the protein makes to catalysis by imposing a planar form to the reduced isoalloxazine (Fig. 8). The most stable conformations of the bound reduced cofactor have the following bending angles: reduced form before the charge relay, 171°, and after the charge relay, 175° (Fig. 8). In the free isoalloxazine, these conformations have energies of 1.9 and 2.5 kcal/mol, respectively, above the energy of the most stable conformation (151°). Thus, interactions with the protein provide close to 2 kcal/mol for the stabilization of the planar reduced form. Equally important, the stable conformation of the free isoalloxazine (bent angle ⫽ 151°) is destabilized in the protein environment by 14.8 kcal/mol (Fig. 8). The interaction of the isoalloxazine ring with rQR1 also has an effect on the energy of reduction. Thus, the analysis can be complemented by the estimation of the difference in energy between the reduction of the isolated/solvated and 431 AB INITIO STUDY OF NAD(P)H:QUINONE REDUCTASE TABLE II. Energies in the Thermodynamic Cycle: Stabilization Energy Reduced FAD configuration FADH2 FADH2 (SCRF) Fig. 9. Stabilization energy and reduction energy of the flavin: the thermodynamic cycle considered for FAD reduction isolated and in the QR environment. This thermodynamic cycle permits us to calculate an experimentally inaccessible quantity, ⌬E2 ⫺ ⌬E4, by the results of the quantum calculations, because ⌬E2 ⫺ ⌬E4 is equal to ⌬E1 ⫺ ⌬E3, and the latter can be found from the calculations. Species QR is the enzyme environment alone, with the oxidized flavin isoalloxazine ring deleted. Species QR⬘ is the same enzyme environment but with the reduced flavin isoalloxazine ring deleted. It is different from QR because Tyr-155 in the protein environment has transferred a proton to the reduced flavin as part of the charge-transfer relay. TABLE I. Species Considered in the Thermodynamic Cycle Reduced (R) Oxidized (OX) Isolated FAD FAD in QR environment QR environment alone FAD (R) FAD (OX) FAD (R) 䡠 QR⬘ FAD (OX) 䡠 QR QR⬘ QR protein-bound isoalloxazine rings (⌬Ebound ⫺ ⌬Efree). This energy difference can be estimated by the study of the thermodynamic cycle shown in Figure 9. Calculations based on this cycle with the species shown in Table I overcome the problem of calculating energy differences between structures that have a different number of atoms. The structures correspond to the cases shown in Figures 3(A,C) and 5(A,C). The energies of the thermodynamic cycle were calculated with the various model chemistries explained in the Materials and Methods section. It is clear from Figure 9 that ⌬E1 ⫺ ⌬E3 ⫽ ⌬E2 ⫺ ⌬E4 where ⌬E1 ⫽ EFAD(R)䡠QR⬘ ⫺ (EQR⬘ ⫹ EFAD(R)) ⌬E3 ⫽ EFAD(OX)䡠QR ⫺ (EQR ⫹ EFAD(OX)) QR⬘ refers to the protein environment QR without the proton, which has been transferred in the first step of the charge relay to O2F. The results of the calculations in vacuo and with an SCRF are shown in Table II. The contribution of the enzyme to the FAD two-electron obligatory reduction with respect to the reduction of the isolated FAD (i.e., ⌬E2 ⫺ ⌬E4, which is equal to ⌬E1 ⫺ ⌬E3) is to make the reduction less favorable by 2.39 – 4.74 kcal/mol (Table II). That means that the reduction potential of the rQR1 flavin is 50 –100 mV more negative than that of the ⌬E1 (kcal/mol) ⌬E3 (kcal/mol) ⌬E1 ⫺ ⌬E3 ⫽ ⌬E2 ⫺ ⌬E4 (kcal/mol) - ⫺2.48 8.18 4.74 2.39 free flavin. This change may be very important for the function of the enzyme. Even with this change of potential, o⬘ the bound flavin can be reduced by NAD(P)H (⌬Efree 䡠 flavin o⬘ ⫽ ⫺180 mV; ⌬Efree 䡠 NAD(P) H ⫽ ⫺320 mV;). However, having a more negative reduction potential, the rQR1bound flavin has a greater driving force for the reduction of quinones, allowing the enzyme to reduce a large variety of quinones, even those with a large negative reduction potential. This conclusion agrees with the observed broad range of quinones that can be reduced by QR1. In addition, a greater driving force may result in a more efficient hydride transfer or higher catalytic rates for that step. This study, in addition to providing valuable insights into the mechanism of QR, presents an excellent example of the contribution that quantum mechanical calculations can make to understanding enzymatic mechanisms. ACKNOWLEDGEMENTS The authors thank Mario Bianchet for his always useful and generous collaboration and B. Brooks and his group for helpful discussions and suggestions. Thanks are also due to the Advanced Biomedical Computing Center of the National Cancer Institute (Frederick, MD) for the Cray computer time and some of the software used and a very efficient and friendly system administration. G. Cavelier’s work was supported in part by an NIGMS/NIH-National Research Service Award. REFERENCES 1. Ross D, Siegel D, Beall H, Prakash AS, Mulcahy RT, Gibson NW. DT-diaphorase in activation and detoxification of quinones. Bioreductive activation of mitomycin C. Cancer Metastasis Rev 1993;12: 83–101. 2. Knox RJ, Boland MP, Friedlos F, Coles B, Southan C, Roberts JJ. The nitroreductase enzyme in Walker cells that activates 5-(aziridin-1-yl)-2,4-dinitrobenzamide (CB 1954) to 5-(aziridin-1-yl)-4hydroxylamino-2-nitrobenzamide is a form of NAD(P)H dehydrogenase (quinone) (EC 1.6.99.2). Biochem Pharmacol 1988;37:-. Knox RJ, Friedlos F, Jarman M, Roberts JJ. A new cytotoxic, DNA interstrand crosslinking agent, 5-(aziridin-1-yl)-4-hydroxylamino2-nitrobenzamide, is formed from 5-(aziridin-1-yl)-2,4-dinitrobenzamide (CB 1954) by a nitroreductase enzyme in Walker carcinoma cells. Biochem Pharmacol 1988;37:-. Tedeschi G, Chen S, Massey V. DT-diaphorase. Redox potential, steady-state, and rapid reaction studies. J Biol Chem 1995;270: 1198 –1204. 5. Li R, Bianchet MA, Talalay P, Amzel LM. The three-dimensional structure of NAD(P)H:quinone reductase, a flavoprotein involved in cancer chemoprotection and chemotherapy: mechanism of the two-electron reduction. Proc Natl Acad Sci U S A 1995;92:8846 – 8850. 6. Talalay P, De Long MJ, Prochaska HJ. Identification of a common chemical signal regulating the induction of enzymes that protect against chemical carcinogenesis. Proc Natl Acad Sci U S A 1988;85:-. 432 G. CAVELIER AND L.M. AMZEL 7. Skelly JV, Sanderson MR, Suter DA, Baumann U, Read M, Gregory DS, Bennett M, Hobbs SM, Neidle S. Crystal structure of human DT-diaphorase: a model for interaction with the cytotoxic prodrug 5-(aziridin-1-yl)-2,4-dinitrobenzamide (CB1954). J Med Chem 1999;42:-. Faig M, Bianchet MA, Talalay P, Chen S, Winski S, Ross D, Amzel LM. Structures of recombinant human and mouse NAD(P)H: quinone oxidoreductases: species comparison and structural changes with substrate binding and release. Proc Natl Acad Sci U S A 2000;97:-. Foster CE, Bianchet MA, Talalay P, Zhao Q, Amzel LM. Crystal structure of human quinone reductase type 2, a metalloflavoprotein. Biochemistry 1999;38:-. Walsh CT. Enzymatic reaction mechanisms. San Francisco: Freeman;- p. 11. Bugg T. An introduction to enzyme and coenzyme chemistry. Oxford: Blackwell Science;- p. 12. Silverman RB. The organic chemistry of enzyme-catalyzed reactions. San Diego: Academic;- p. 13. Lee T-S, York DM, Yang W. Linear-scaling semiempirical quantum calculations for macromolecules. J Chem Phys 1996;105:2744 – 2750. 14. Sánchez-Portal D, Ordejón P, Artacho E, Soler JM. Densityfunctional method for very large systems with LCAO basis sets. Int J Quantum Chem 1997;65:-. Ordejón P, Artacho E, Soler JM. Self-consistent order-N densityfunctional calculations for very large systems. Phys Rev B 1996;53: R10441–R-. Ordejón P, Drabold DA, Martin RM, Grumbach MP. Linear system-size scaling methods for electronic-structure calculations. Phys Rev B 1995;51:1456 –1476. 17. Kohen A, Klinman JP. Enzyme catalysis: beyond classical paradigms. Acc Chem Res 1998;31:-. Kohen A, Cannio R, Bartolucci S, Klinman JP. Enzyme dynamics and hydrogen tunneling in a thermophilic alcohol dehydrogenase. Nature 1999;399:496 – 499. 19. Ringe D, Petsko GA. Tunnel vision. Nature 1999;399:-. Yagi K. Hydrogen bonding in flavoproteins. Chem Scr A 1987;27: 27–30. 21. Mestres J, Duran M, Bertran J. Characterization of the transition state for the hydride transfer in a model of the flavoprotein reductase class of enzymes. Bioorg Chem 1996;24:69 – 80. 22. Zheng Y-J, Ornstein RL. A theoretical study of the structures of flavin in different oxidation and protonation states. J Am Chem Soc 1996;118:-. Niemz A, Imbriglio J, Rotello VM. Model systems for flavoenzyme activity: one- and two-electron reduction of flavins in aprotic hydrophobic environments. J Am Chem Soc 1997;119:-. Breinlinger EC, Keenan CJ, Rotello VM. Modulation of flavin recognition and redox properties through donor atom interactions. J Am Chem Soc 1998;120:8606 – 8609. 25. Vishveshwara S, Jacob R, Nadig G, Maizel JV Jr. The role of lysine-41 in RNase A catalysis. A quantum chemical study on the active site–ligand complex. J Mol Struct 1998;471:1–11. 26. Antony J, Medvedev DM, Stuchebrukhov A. Theoretical study of electron transfer between the photolyase catalytic cofactor FADH⫺ and DNA thymine dimer. J Am Chem Soc 2000;122:-. Várnai P, Richards WG, Lyne PD. Modelling the catalytic reaction in human aldose reductase. Proteins 1999;37:218 –227. 28. Chatfield DC, Eurenius KP, Brooks BR. HIV-1 protease cleavage 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. mechanism: a theoretical investigation based on classical MD simulation and reaction path calculations using a hybrid QM/MM potential. J Mol Struct 1998;423:79 –92. Liu H, Müller-Plathe F, Gunsteren WFv. A combined quantum/ classical molecular dynamics study of the catalytic mechanism of HIV protease. J Mol Biol 1996;261:454 – 469. Tapia O, Andrés J, Safont VS. Transition structures in vacuo and the theory of enzyme catalysis. Rubisco’s catalytic mechanism: a paradigmatic case? J Mol Struct 1995;342:131–140. Li J, Nelson MR, Peng CY, Bashford D, Noodleman L. Incorporating protein environments in density functional theory: a selfconsistent reaction field calculation of redox potentials of [2Fe2S] clusters in ferredoxin and phthalate dioxygenase reductase. J Phys Chem A 1998;102:-. Rick SW, Abashkin YG, Hilderbrandt RL, Burt SK. Computational studies of the domain movement and the catalytic mechanism of thymidine phosphorylase. Proteins 1999;37:242–252. Walsh CT. Flavin coenzymes: at the crossroads of biological redox chemistry. Acc Chem Res 1980;13:148 –155. Ghisla S, Massey V. Mechanisms of flavoprotein-catalyzed reactions. Eur J Biochem 1989;181:1–17. Murthy YV, Massey V. Synthesis and properties of 8-CN-flavin nucleotide analogs and studies with flavoproteins. J Biol Chem 1998;273:-. Meyer M. Ab initio study of flavonoids. Int J Quantum Chem 2000;76:724 –732. Meyer M. Density functional study of isoalloxazine and C4ahydroperoxidihydroisoalloxazine. J Mol Struct 1997;417:163–168. Wouters J, Durant F, Champagne B, André J-M. Electronic properties of flavins: implications on the reactivity and absorption properties of flavoproteins. Int J Quantum Chem 1997;64:721– 733. Nakai S, Yoneda F, Yamabe T. Theoretical study on the lowestfrequency mode of the flavin ring. Theor Chem Acc 1999;103:109 – 116. Frisch MJ, Trucks GW, Schlegel HB, et al. Gaussian 98. Pittsburgh, PA: Gaussian;- p. Becke AD. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys Rev A 1988;38:3098 –3100. Lee C, Yang W, Parr RG. Development of the Colle–Salvetti correlation-energy formula into a functional of the electron density. Phys Rev B 1988;37:785–789. Bruice TC. Mechanisms of flavin catalysis. Acc Chem Res 1980;13: 256 –262. Nishimoto K. Ab-initio MO study of the redox function of flavin. In: Nagata C, editor. Biomolecules. Tokyo: Japan Scientific Societies Press; 1985. p 9 –19. Wong MW, Frisch MJ, Wiberg KB. Solvent effects 1. The mediation of electrostatic effects by solvents. J Am Chem Soc 1991;113: 4776 – 4782. Foresman JB, Frisch A. Exploring chemistry with electronic structure methods. Pittsburgh, PA: Gaussian;- p. Werner P-E, Ronnquist O. The crystal structure of 5-acetyl-9bromo-1,3,7,8,10-pentamethyl-1,5-dihydroisoalloxazine. Acta Chem Scand 1970;24:-. Lennon BW, Williams CH, Ludwig ML. Crystal structure of reduced thioredoxin reductase from Escherichia coli: structural flexibility in the isoalloxazine ring of the flavin adenine dinucleotide cofactor. Protein Sci 1999;8:2366 –2379.