Mechanistic Studies of a Flavin Monooxygenase: Sulfur Oxidation of

Research Article Cite This: ACS Catal. 2018, 8, 9298−9311

pubs.acs.org/acscatalysis

Mechanistic Studies of a Flavin Monooxygenase: Sulfur Oxidation of Dibenzothiophenes by DszC Ana C. C. Barbosa,‡ Rui P. P. Neves,‡ Seŕ gio F. Sousa,§ Maria J. Ramos, and Pedro A. Fernandes* UCIBIO,REQUIMTE, Departamento de Química e Bioquímica, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, s/n, 4169-007 Porto, Portugal

ACS Catal. Downloaded from pubs.acs.org by UNIV OF SUNDERLAND on 09/07/18. For personal use only.

S Supporting Information *

ABSTRACT: Flavin monoxygenases (FMOs) are enzymes of increasing biotechnological (e.g., crude oil biodesulfurization) and pharmacological (e.g., drug metabolism) interest that perform the oxidation of soft nucleophiles and play key roles in the excretion of xenobiotics or in sulfur amino acid metabolism. DszC is a key FMO involved in sulfur oxidation of dibenzothiophenes (DBTs) through the 4S metabolic pathway of some bacteria. This pathway can be a cheaper and greener alternative for sulfur removal, as DBTs are the major source of crude oil sulfur. Here, we investigate the reaction mechanism of DszC with quantum mechanics/molecular mechanics methods (SCS-MP2/def2-TZVPP:ff10// B3LYP/6-31G(d):ff10). We observe that the reaction mechanism of DBT oxidation occurs in three stages: (1) spin-forbidden formation of a C4a-hydroperoxyflavin (C4aOOH) intermediate; (2) oxidation of DBT to DBTO, upon nucleophilic attack of the DBT-sulfur on the distal oxygen of C4aOOH; and (3) proton transfer from the N5H of the flavin group to the His92-imidazole via Ser163-hydroxyl, releasing a water molecule and oxidized flavin mononucleotide. The overall reaction is computed to be exergonic (−38.7 kcal·mol−1), and the rate-limiting step is the oxidation of DBT to DBTO (ΔG‡ = 19.7 kcal·mol−1, consistent with the experimental turnover of 1.6 min−1). We observe that oxygen activation is a nearly spontaneous process that occurs through a proton-coupled electron transfer to produce a hydroperoxyl radical, followed by a triplet-singlet spin-forbidden inversion to form the C4aOOH intermediate. In agreement with other studies, His391 is a key acid to activate O2 and form the covalent bond. Further clarifying previous mutagenesis results, we also propose that His92 and Tyr96 are key residues for the mechanism: His92 acts as acid to deprotonate N5H in flavin via Ser163; and Tyr96 enhances the oxidation of DBT-sulfur by anchoring the proximal oxygen of C4aOOH, and acts as acid to form the water byproduct and regenerate the flavin cofactor. These are important results to clarify the chemistry of flavin monoxygenases and to open doors for the rational design of DszC mutants with improved catalytic activity. KEYWORDS: biodesulfurization, oxygen activation, sulfur oxidation, flavin monoxygenase, quantum mechanics/molecular mechanics

■

INTRODUCTION Desulfurization of Organic Sulfur Compounds by the Oil Industry. Organic sulfur compounds, such as mercaptans, thiols, sulfides, disulfides, and thiophenes, are found in crude oil, in amounts that can range from 0.03 to 7.89% (w/w).1,2 The combustion of these compounds leads to the emission of sulfur oxides (SOX), with known environmental and health problems.3,4 These sulfur-containing compounds have also been related to several complications in petroleum refining, such as catalyst poisoning or corrosion in refinery equipment.2,5 Currently, one of the main challenges of this industry is to achieve increasingly lower sulfur content in transportation fuels, as imposed by regulatory agencies.5−8 Hydrodesulfurization has been the more commonly employed chemical approach; however, it has major drawbacks, such as high costs, massive CO2 emissions, reduction of the fuel’s energetic value, and the difficulty in removing sulfur from dibenzothio© XXXX American Chemical Society

phene (DBT) and its alkylated derivatives (DBTCs), unless very high temperature/pressure conditions are used.1,9 Unfortunately, DBTCs are the main organic sulfur compounds found in crude oil, accounting for up to 70% of these1,10 As a result, the development of alternative or complementary desulfurization techniques that can overcome the inefficiency of hydrodesulfurization and produce cleaner fossil fuels at lower cost and with less CO2 emissions is a research area of great interest.5 Biodesulfurization: The 4S Pathway and the Role of DszC. In the last decades, biodesulfurization has drawn wide attention1,11,12 because it is a much greener technique for sulfur removal in fossil fuels. In particular, the metabolism of Received: May 15, 2018 Revised: August 9, 2018

9298

DOI: 10.1021/acscatal.8b01877 ACS Catal. 2018, 8, 9298−9311

Research Article

ACS Catalysis

Figure 1. (A) DszC tetramer (PDB ID: 3X0Y), highlighting the two homodimers (chains AF in the front, chains GH at the back); the FMN cofactor is shown in van der Waals representation. (B) The catalytic site of one monomer of DszC (chain A) in ball-and-stick representation, with the most relevant distances identified for the reaction mechanism of DszC. (C) The inner and outer chambers of one monomer of DszC (chain A); the binding pose of DBT (as proposed by Zhang et al.) is represented, with DBT in white ball-and-stick representation.

Scheme 1. Proposed General Catalytic Cycle of DszC of Zhang and Co-workers, Based on the Catalytic Mechanism Proposed for the Similar Oxygenase Component (C2) of p-Hydroxyphenylacetate Hydroxylase (HPAH):22,31 (A) Mechanism of Oxygen Activation through a Proton-Coupled Electron Transfer Mechanism Followed by TS Spin Transition; (B) DBT Oxygenation by Distal Oxygen (Od) of C4a-Hydroperoxyflavin Intermediate to DBT-Sulfur; (C) N5-Deprotonation and Water Elimination To Regenerate FMN

efficiency, 1.3 μM−1 min−1 (kcat = 1.6 ± 0.3 min−1, KM = 1.4 ± 0.3 μM) and 1.1 μM−1 min−1 (kcat = 1.7 ± 0.2 min−1, KM = 1.3 ± 0.3 μM) .12 DszC has a key role in the 4S pathway, because it has been referred to as one of its kinetic bottlenecks.12 It is the first enzyme of the pathway, and it catalyzes the successive oxidation of DBT to DBT sulfoxide (DBTO) and after to DBT sulfone (DBTO2). This reaction requires FMNH2 as a cofactor (provided by DszD) and molecular oxygen. X-ray crystallographic studies on its apo and holo forms have provided essential structural and mechanistic insights. DszC consists of a dimer of homodimers, with two similar catalytic sites per homodimer, separated by 50 Å (Figure 1A). The catalysis takes place close to the interface of the homodimer and involves the isoalloxazine ring of FMNH2, the DBT substrate, and residues from one monomer (the conserved His92, Tyr96, Asn129, Ser163 His388, and His391) Figure 1B. The phosphate moiety of FMNH2 is stabilized by a network of hydrogen bonds with the second monomer in the homodimer (Arg282, Ala369, and Arg370).21−23 No X-ray structure of the Rhodococcus erythropolis strain IGTS8 DszC complexed with DBT was obtained up to now.21−24

organic sulfur compounds by the bacterium Rhodococcus erythropolis strain IGTS8 has been regarded with great interest. This strain has the ability to utilize DBTCs as a source of sulfur, without degrading the DBTCs carbon skeleton, thus preserving the energetic value of the fuel.13−15 This highly specific metabolic pathway is often referred to as the 4S pathway and involves four enzymes: two flavin mononucleotide (FMN)H2-dependent monooxygenases, DszA and DszC, the desulfinase, DszB, and the NADH-dependent FMN oxidoreductase, DszD. In the 4S pathway, DBT is converted into 2-hydroxybiphenyl (HBP) and sulfite (SO32−). This pathway occurs under mild conditions, it yields no undesirable end products, and requires low capital and operating costs; hence, it has been suggested as a green alternative for the improvement of the quality of fossil fuels.14−20 However, the rate of biodesulfurization is slow compared to hydrodesulfurization, thus limiting its implementation in largescale processes.8 Its rate has to be increased at least 500 fold to be of use in oil refining plants.18 An increase of the enzymatic activity of the enzymes in the 4S pathway is a crucial step in this regard. In particular, DszC and DszB exhibit a low catalytic 9299


Research Article

ACS Catalysis

active site? (5) what is the rate-limiting step? (6) what roles do the conserved His92, Tyr96, Asn129, and His391 play in the reaction? There is a lot to wonder about what concerns the complex mechanism of oxidation of DBT by DszC. In an effort to bring some light to these issues, we have performed hybrid quantum mechanics/molecular mechanics (QM/MM), full QM calculations, and classical molecular dynamics simulations. The results have substantially expanded the knowledge on the reaction mechanism of DszC and the mechanism of oxygen activation by flavin-dependent monoxygenases. This knowledge will be of future use in the development of efficient biocatalysts for desulfurization, either by mutagenesis or through de novo synthesis.

DszC likely follows an ordered sequential mechanism, comprising three major steps:22 (1) the activation of molecular oxygen by FMNH2 forming a C4a-(hydro)peroxyflavin intermediate; (2) the oxidation of DBT to DBTO; and (3) the protonation of the C4a-hydroxyflavin (C4aOH) to yield oxidized FMN (FMNox) and a water molecule (Scheme 1).16,25 The FMNox cofactor will then be reduced to FMNH2 by DszD, restoring DszC for a new catalytic cycle. Mechanism of Oxygen Activation in DszC. The mechanism of oxygen activation in flavin-dependent enzymes is present in a wide variety of biological reactions, and its chemistry is a subject of ongoing research.26−28 The original reaction pathway was proposed by Bruice et al.29 and Massey et al.,30 and it comprised the formation of a C4a-peroxyflavin (C4aOO−) intermediate, which would be protonated to produce C4a-hydroperoxyflavin (C4aOOH). Recently, Wongnate et al. reported a new viewpoint on this activation mechanism in pyranose 2-oxidase (P2O).27 In this study, a positively charged His, located near the C4a-carbon of the isoalloxazine ring, is involved in a proton-coupled electron transfer (PCET) mechanism that reduces the molecular oxygen (O2) to a •OOH radical. This is followed by the formation of the C4aOOH intermediate. Later calculations from Visitsatthawong et al. on the oxygenase component (C2) of p-hydroxyphenylacetate hydroxylase (HPAH) enforced that His396 was similarly placed to interact with the molecular oxygen and also proposed that the formation of the C4aOOH would take place via the PCET mechanism.28 Similarly to the C2 of HPAH, the DszC from Rhodococcus erythropolis (ReDszC) has a conserved His391 near the C4a-carbon of the isoalloxazine ring (the Nε of His391 is 4.23 Å away from the C4a-flavin, Figure 1B), suggesting that DszC may follow a similar PCET mechanism for the activation of molecular oxygen (formation of C4aOOH intermediate, steps 1−3 in Scheme 1).21 However, it is not clear if the necessary spin inversion leading to the formation of the C 4a OOH intermediate occurs vertically or through the cross-section of the singlet and triplet surfaces describing the formation of the C4a−OpOH bond.27,28 Oxidation of DBT by DszC. The C2 of HPAH shows an active site similar to that of DszC (Figure S2 in SI) and has been thoroughly used as a reference to predict the reactivity of DszC.21,22 In analogy to the proposed mechanism for the C2 of HPAH,28,31,32 the oxidation of DBT to DBTO in DszC should occur through nucleophilic attack of the DBT sulfur to the Od of the C4aOOH, as shown in step 4 of Scheme 1. Similarly to Ser171 in C2 of HPAH, the side chain of Ser163 makes a hydrogen-bond with the N5-nitrogen of the isoalloxazine ring (2.82 Å, Figure 1B), suggesting that FMN can also be restored from the deprotonation of N5-nitrogen by the leaving C4ahydroxide, via the Ser163,27,28 leading to the elimination of a water molecule and formation of oxidized FMN, which can be recycled by DszD to undergo a new cycle (formation of oxidized FMN in Scheme 1). Nevertheless, there is a large lack of chemical insight (structural, chemical and thermodynamic) on how all these transformations take place, and a set of unanswered questions, such as (1) what is the mechanism of molecular oxygen activation in DszC? (2) during molecular oxygen activation, does the triplet-to-singlet spin inversion occurs vertically or is there a minimum energy crossing point (MECP) in which it occurs? (3) how is DBT placed in the active site? (4) does the mechanism of oxygen activation by DszC require DBT at the

■

METHODOLOGY Molecular Models for the DszC Complex. We modeled the system starting from an X-ray structure of FMN-bound DszC from Rhodococcus erythropolis (PDB ID: 3X0Y, 2.3 Å resolution).23 This crystal structure contained two tetramers in the asymmetric unit, with four equivalent catalytic sites. The association of the tetramers did not seem to be important for the catalysis.23 Each tetramer was composed by two catalytically competent homodimers with 417 amino acid residues per monomer. Hence, to run molecular mechanics calculations, our initial model of DszC contained all the residues of one catalytically competent homodimer (chain A and F), and all crystallographic water molecules within a range of 6 Å from the catalytic site. We modeled the three states of FMN (FMNH2, FMNH−, and the C4a-hydroperoxyflavin intermediate, which is C4aOOH) bound to the catalytic site using the oxidized FMN molecule present in the X-ray structure of DszC. The fact that we did the initial modeling using the oxidized form of FMN should not constitute a problem, as the monomer that we picked to model the active site (chain A) displayed a welldefined electronic density for FMN and a closed-loop that was described to tightly lock FMN in the active site of DszC. As such, and since the isoalloxazine ring of FMN was tightly bound through hydrogen bonds, it can be safely assumed that the binding mode of reduced FMN and C4aOOH (aside from the pyrimidin-like ring of isoalloxazine) are very close to of the crystallographic pose of oxidized FMN. As the reaction of DszC with DBT was likely to require modification of the FMN cofactor before DBT binding, we modeled the C4aOOH intermediate with and without DBT in the catalytic site to study the mechanism of oxygen activation in both cases. There was no X-ray structure of the Rhodococcus erythropolis strain IGTS8 DszC complexed with FMN and DBT. Experiments suggested two likely binding poses for DBT,21−23 even though the one proposed by Zhang et al. was clearly the more adequate to accommodate DBT (Figure 1C).23 In Zhang’s proposal, the binding pocket was composed of two chambers: an inner chamber composed of hydrophobic residues and an outer chamber composed mostly by hydrophilic residues (Figure 1C).22 It was a large pocket where DBT binds to the DBT-sulfur (SDBT) close to the C4a-carbon. Given the similarity between the active site and the fold of DszC and the C2 of HPAH, we modeled DBT in the pose of the hydroxyphenyl acetate (HPA) substrate (Figure S2 in SI) of the C2 of HPAH from Acinetobacter baumanni (PDB ID: 2JBT, 2.8 Å resolution).31 This enzyme showed a similar binding pocket to that proposed by Zhang et al., and it exhibited a 9300


Research Article

ACS Catalysis

Figure 2. Poses for DBT modeled from the C2 of HPAH (on the left) and from the ortholog TdsC (on the right). Distances between the nucleophilic sulfur of DBT and the distal oxygen of C4aOOH are indicated in Å.

AMBER force field. Coordinates for FMNH2, FMNH−, and C4aOOH were determined from a geometry optimization in vacuum at the HF/6-31G(d,p) level, where only the isolloxazine region was freely optimized. Coordinates for DBT were also obtained from an optimization in vacuum at the HF/6-31G(d,p) level. The atomic point charges for FMNH2, FMNH−, C4aOOH, and DBT were calculated in vacuum at the HF/6-31G(d,p) level with the restrained electrostatic potential (RESP) fitting,37 using the Antechamber tool. All parameters are included in SI. The Xleap module of AMBER12 was used to protonate the complex, assuming physiological protonation states for all residues (as predicted with Propka3.1).38 The His92, His388, and His391 were protonated at the Nδ (refer to section B in SI). Thirty-three Na+ ions were added to neutralize the system (two N-terminal residues Ala-Asp were present only in chain A). The complex was surrounded by a rectangular box of TIP3P water molecules39 within a minimum distance of 12 Å from the protein. A minimization and an equilibration protocol was applied to prepare the system for MD simulations. All MD productions were performed in the NPT ensemble (310 K and 1 bar). Details follow in section A of SI. The QM calculations were performed with Gaussian09.40 The AMBER 12 software41 was used to perform the classical mechanics calculations and run the conventional molecular dynamics (MD) simulations. QM/MM Model. Several computational methods have been successfully employed in the modeling of enzyme-catalyzed reactions.42−47 We have been able to clarify several enzyme reaction mechanisms using a combination of quantum mechanics/molecular mechanics (QM/MM) modeling and the N-layered integrated molecular orbital and molecular mechanics (ONIOM) methodology.48−50 Three aspects were carefully taken into account during the construction of the QM/MM model: (i) the initial protonation state of bound FMN; (ii) the position of O2 in relation to the C4a-carbon of FMN; (iii) the possibility of oxygen activation to occur prior or after DBT binding to the active site. We modeled two complexes to evaluate the most probable coupled protonation state of FMN and of the nearby His391, which is unclear from the X-ray structure and will have an impact in the mechanism: DszC:FMNH2 with His391protonated at the Nδ (δ-His391) and DszC:FMNH− with His391 positive (p-His391). Our calculations have shown that pHis391 was the most likely (refer to section B in SI). This

similar folding to DszC. The root-mean-square deviation (rmsd) for the backbone of the homodimeric complex was 3.1 Å (over 3617 pairs of heavy atoms) and 0.71 Å regarding residues in the active site (considering 38 pairs of heavy atoms−backbone and Cβ of His92(Trp112 in C2), Ser163(Ser171 in C2) and His391(His396 in C2), and heavy atoms of FMN). On a side note, a structure of an enzyme similar to DszC (TdsC from Paenibacillus sp., with 64% identity and identical active site residues) with DBT bound was very recently published (PDB ID: 5XDE).33 The pose of DBT in this structure showed some differences in relation to the X-ray structure of the parent enzyme HPAH from Acinetobacter baumanni used for the modeling here (see Figure 2). As there are no structures of any flavin-dependent monoxygenase complexed with both C4aOOH and substrate, choosing the correct pose of DBT in any DszC full complex just by direct comparison with the incomplete X-ray structures becomes doubtful. Furthermore, an inspection of the electron density map obtained by Hino and co-workers, with a resolution of 1.6 Å, showed an electron density for DBT less defined than that of other residues in the active site of TdsC, suggesting that the binding pose of DBT should not be very rigid and well-defined. To be on the safe side, we repeated the DBT oxidation reaction step in DszC, with DBT bound in a position modeled from the one found in TdsC (following the same modeling procedure as described). The plane of the DBT pose from TdsC (Figure 2, on the right) was approximately perpendicular to the one modeled from HPA and resulted in a configuration that seemed more favorable to carry the SN2 nucleophilic attack by the DBT-sulfur −3.09 Å from the distal oxygen of C4aOOH (Od) and an angle of 164.1° for the nucleophilic attack. Aside from that, the active site configuration was similar for both DBT poses (refer to section C in SI for details). Our energies calculated at the RI-SCSMP2/def2-TZVPP:ff10 level showed that this pose was in fact more favorable for the oxidation of DBT than the one modeled from HPA. However, we observed very similar geometry changes during the reaction, regardless of the starting pose of DBT. Classical MD Simulations of the DBT Complex. The ff10 force field34 was employed to parametrize the enzyme, while the DBT, FMNH2, FMNH− and C4aOOH parameters were taken from the General AMBER Force Field (GAFF).35 Parameters for oxidized FMN (FMNox) were taken from the literature.36 We followed the parametrization procedure of the 9301


Research Article

ACS Catalysis

Figure 3. Representation of the QM/MM model built to study the mechanism of DBT oxidation by DszC (on the left), where red spheres represent the water molecules within 6 Å of the active catalytic site; and the DFT layer employed in the calculations (on the right), where the thinwhite lines represent the bonds that were capped with hydrogen link atoms.

active site (total of 112 atoms in DFT layer, with zero overall charge, in the singlet spin state). The MM layer comprised all the remaining atoms of the model. The system was further divided into a frozen region and a free region, keeping in the free region all residues within a radius of 15 Å from the catalytic site. No water molecules were in close contact with the catalytic site. Additionally, to evaluate the effect of the frozen region in the rate-limiting step of the catalytic mechanism of DszC, we built two more models for the DszC:C4aOOH:DBT complex: (1) a model with all residues kept free, except those in the second monomer beyond a radius of 15 Å from the catalytic site; and (2) a model in which all residues in the second monomer beyond a radius of 15 Å from the catalytic site, and a 3 Å cap of water molecules surrounding the protein were kept frozen, while the remaining system was kept free. The input coordinate files for each model are included in SI. QM/MM Calculations. All QM/MM calculations were carried out at the B3LYP/6-31G(d):ff10 level of theory,51−59 with the Coulomb QM/MM interactions calculated with the electrostatic embedding scheme, as implemented in Gaussian 09.40 The B3LYP/6-31G(d) level has been successfully employed in other QM/MM studies of the catalytic mechanism of enzymes.49,60−62 In particular, to study the triplet-singlet spin inversion taking place upon formation of the C4aOOH, we allowed breakage of the symmetry of the singlet configuration by accepting the inclusion of imaginary coefficients and/or mixed-spin orbitals in the wave function.63−65 We used hydrogens as link atoms to complete the valence of the bonds in the boundary of the two layers.66−69 All stationary points of the reaction mechanism were retrieved from full geometry optimization. Vibrational analysis confirmed the nature of each stationary point, with a single imaginary frequency for transition states and all real vibrational frequencies for minima. The harmonic zero-point energy (ZPE), the entropy and thermal energy at 298.15 K were used to estimate the Gibbs energy within the single conformation approach. The translational, rotational and vibrational contributions for the Gibbs energy were determined with the particle-in-a-box, rigid rotor and the combined free-rotor/ harmonic oscillator approximations. We only considered vibrational temperatures larger than 100 K for the calculation of the vibrational entropy, including 14 750 vibrations for minima (14 545 in the absence of DBT), and 14 749 for TSs.

observation is also supported by previous studies on C2:HPAH, which indicated that upon mutation of this residue, there was a lower concentration of C4aOOH and lower hydroxylation efficiency. Hence, we assumed that FMNH2 is deprotonated by His391, after it enters the active site. On the second aspect of the modeling, as O2 is a small molecule, it was difficult to predict its precise position within the active site. Hence, we built the QM/MM models starting from the C4a-hydroperoxyflavin-intermediate (C4aOOH) intermediate, which had the O2 molecule already bound to FMN and therefore with a position that was obvious and clearly assignable, and we studied the mechanism of activation of O2 in the reverse direction. Finally, to study the mechanism of oxygen activation, we built two QM/MM models for the reaction of formation of C4aOOH: a first with DszC complexed with C4aOOH intermediate only (DszC:C4aOOH) and a second with DszC complexed with the C4aOOH and DBT(DszC:C4aOOH:DBT). Our QM/MM model was cut from the previously modeled homodimer, after energy minimization at the MM level, as it has been shown that the use of the X-ray conformation is usually the best option for single-conformation high-level calculations.50 We kept in the model one monomer binding FMN and DBT (except when testing the mechanism of oxygen activation in the absence of DBT), all residues from the remaining monomers within a radius of 18 Å from the catalytic site (constituted by His92, Tyr96, Asn129, Ser163, His388 and His391, FMN, and DBT), and the water molecules at least 6 Å away from the catalytic site (Figure 3). We have divided the system in two layers: the QM layer, treated at the DFT level, and the MM layer, treated at the classical level with the ff10 force field.34 The QM layer includes all atoms from the DBT substrate (when present), the isoalloxazine ring (except for the two methyl groups in the phenyl moiety), the N10-bound hydroxymethyl of the ribitol tail of the C4a-hydroperoxyflavin intermediate (C4aOOH), and the side chains of the residues involved directly or indirectly in the reaction (His92, Tyr96, Asn129, Ser163, His388, and His391). Two slightly different QM/MM models were used: to study the formation of the C4aOOH intermediate, we did not include DBT (total of 91 atoms in the DFT layer, with zero overall charge, in the singlet spin state); to study the oxidation of DBT, we included both C4aOOH and DBT in the 9302


Research Article

ACS Catalysis

Figure 4. (A) DFT region of the metastable [REACTO2]* with FMNH• and the superoxide anion (O2•−) upon deprotonation of His391 by the superoxide anion (the point used corresponds to a N−H bond of His391 around 1.0 Å); (B) DFT region for the REACTOOH state and (C) INTC4aOOH state. The distances are indicated in Å. Blue-to-red shades in atoms represent the most relevant changes in atomic charge relative to the [REACTO2]* (from 0.07 to 0.60 au). Blue stands for a more positive atomic charge, and red stands for a more negative atomic charge.

singlet and triplet configurations and the singlet−triplet/ triplet-singlet (S-T/T-S) vertical transitions along these PESs. Then, we retrieved the energy crossing points (ECP) in these vertical transitions to choose starting geometries for our DFT cluster model. These ECPs correspond to points where the geometry and energy of the QM/MM model are the same for the singlet and triplet configurations. That said, they may not correspond to the minimum energy crossing point (MECP) that characterizes triplet-singlet spin inversion, but they are points in the seam line formed by the intersection of the singlet and triplet PESs. All valences for the capped atoms of the DFT layer of the QM/MM model were filled with hydrogen atoms. The DFT cluster comprised 106 atoms and zero overall charge. For each of the ECPs determined, we used the code developed by Jeremy Harvey89,90 to search for the MECP for the spin inversion from triplet to singlet configuration, in the formation of the C4aOOH. The MECP was characterized at the B3LYP/ 6-31G(d) level.

The free-rotor and harmonic oscillator approximations were combined through a frequency-dependent Head-Gordon damping function,70 for a reference vibrational temperature of 120 K (as previously reported by Grimme).71 We carried out single-point energy calculations with the larger triple-ζ basis set 6-311+G(2d,2p)72−74 and the B3LYP density functional. We have also tested other density functionals widely used for precise activation and reaction energy calculations (e.g., mPWB1K,75 M06,76 M06-2X,76 or BMK77); however, we observed significant energy differences for the barrier of DBT oxidation as calculated by these density functionals (between 3 to 13 kcal·mol−1 comparatively to B3LYP). Consequently, we performed single-point energy calculations78,79 on DFT clusters taken from the DFT layer of the QM/MM model (all link atoms were replaced with −CH3) with the resolution of the identity spin component scaled MP280,81 (RI-SCS-MP2) and B3LYP methods, and the def2TZVPP basis set82 combined with the Coulomb-fitting auxiliary basis set def2/J83 to employ the resolution of the identity approximation (RIJCOSX84 for RI-SCS-MP2 and RIJONX85 for B3LYP), with the ORCA 4.0.1 software.78,79 The final Gibbs energy profile was obtained at the SCS-MP2/ def2-TZVPP:ff10//B3LYP/6-31G(d):ff10 level of theory. Atomic charges were calculated at the B3LYP/ 6-311+G(2d,2p) level, in all the stationary points with the HirshfeldEE extension of the Hirshfeld scheme,86−88 as implemented in Gaussian 09. In particular, singly occupied α/β orbitals were calculated at the B3LYP/6-31G(d) level of theory for the stationary points of the formation of the C4ahydroperoxyflavin, upon biorthogonalization of the calculated molecular orbitals. Natural bond orbital analysis (NBO) was also conducted for the DFT layer of all the stationary points of the mechanism of reaction of DBT oxidation by DszC at the B3LYP/6-31G+(d,p) level of theory. DFT Calculations on Cluster Models. To address the study of the mechanism of oxygen activation, we built a cluster model from the DFT layer of the QM/MM model without DBT. Our goal was to address the possible spin inversion mechanism for the formation of the C4a-hydroperoxyflavin intermediate in DszC. We studied the QM/MM potential energy profiles (along the putative reaction coordinate) for the

■

RESULTS AND DISCUSSION First Stage: Formation of the C4a-Hydroperoxyflavin Intermediate. Steps 1, 2, and 3 in Scheme 1. A triplet spin configuration was the most favorable for the system after O2 was modeled inside the active site ([REACTO2]* in Figure 4A). However, this state was not a stationary point and, as soon as O2 entered the active site, the deprotonation of His391 by O2 occurred with no energy barrier. This resulted in a neutral OOH species (REACTOOH in Figure 4B). An analysis of the total Hirshfeld charge in FMNH, His391 and O2 (Table 1) before (preREACT) and after O2 entered Table 1. Total Hirshfeld Charge (a.u.) for His391, FMNH, and O2 for the Formation of the C4a-Hydroperoxyflavina preREACT [REACTO2]* REACTOOH INTC4aOOH a

9303

His391

FMNH

O2

0.75 0.60 0.11 0.14

-0.98 -0.33 −0.26 −0.38

− -0.39 -0.03 −

The most relevant changes in charge are highlighted in bold. DOI: 10.1021/acscatal.8b01877 ACS Catal. 2018, 8, 9298−9311

Research Article

ACS Catalysis

Figure 5. (A) Potential energy surface (PES) for the C4a···Op distance for the triplet state (Tαα) and the singlet state (S) configurations, in colorcontoured lines. The PES for the lowest energy path for the formation of the C4a···Op bond is represented in bold gray lines. The geometries optimized in the Tαα configuration are represented in blue, and the geometries optimized in the S configuration are represented in orange. The PES for the vertical spin transition between S and Tαα in the S geometry (STαα) and TααS and between Tαα and S in the Tαα geometry (TααS) are represented in dotted lines; the corresponding energy crossing points (ECPS and ECPTαα) are highlighted with blue and orange circles. (B) Optimized minimum energy crossing point (MECP) obtained from geometry optimization of the DFT cluster from ECPS (orange) and from ECPTαα (cyan blue), calculated at the B3LYP/6-31G(d) level. Average distances are indicated in Å.

the active site ([REACTO2]*) showed that the negative charge in FMNH was partially transferred to O2, with O2 exhibiting a partial superoxide character (−0.39 au), and FMNH showing a small partial negative charge (−0.33 au). The transfer of electron density to O2 was particularly evident when comparing the negative charge of FMNH in the preREACT and REACTOOH complexes (−0.98 and −0.33 au respectively). The Od−Op bond length (1.33 Å) was also consistent with the single bond character of the superoxide anion (1.35 Å), with Od and Op being the distal and the proximal oxygen to the C4a of FMNH (refer to Figure 4C). This proton-coupled electron transfer mechanism (PCET) has been often observed in the literature for similar flavindependent (mono)oxygenases;27,28 however, these earlier studies indicate that the electron transfer occurs after the protonation of O2. Alternatively, we propose that the partial electron transfer from FMNH to O2 preceded the deprotonation of His391 by its distal oxygen (Od). An analysis of the highest occupied molecular orbitals at [REACTO2]* provided additional evidence that the triplet state is a result of two doublet radicals in two different molecules: there were two semioccupied orbitals with highest occupancy for Op (0.50) and Od (0.47) in O2, and for N5 (0.13) and C4a (0.12) in FMNH. The same analysis at the REACTOOH state (after the formation of OOH) also revealed two semioccupied orbitals: one encompassed mostly OOH (occupancy of 0.67 for Op and 0.38 for Od), and the other presented similar contributions from the N5 and C4a of the isoalloxazine moiety in FMNH (occupancy of 0.23 and 0.28). Moreover, the overall charge of OOH and FMNH were also closer to neutral. Hence, OOH and FMNH are best described as two doublet radicals. In addition, the hydrogen bond between the imidazole of His391 and OOH (1.76 Å) anchored the OOH in a position close to the isoalloxazine of FMNH (Figure 4B). We calculated 1D potential energy surfaces (PES) for the subsequent formation of the C4a−Op covalent bond in the singlet (S) and triplet (Tαα) configurations (Figure 5A). The analysis of the spin state of the expected reactants and products

showed that the REACTOOH state was more stable in the triplet configuration (ΔESTαα = −13.2 kcal·mol−1), while the INTC4aOOH state was much more stable in the singlet state (ΔESTαα = +59.2 kcal·mol−1). Additionally, in the singlet state, REACTOOH was converted to INTC4aOOH in a barrierless process (see orange-contoured line in Figure 5A), but in the triplet state the REACTOOH state was a stationary point. Therefore, a spin-forbidden reaction, with an inversion from the triplet to the singlet state, was due to occur. To understand the process, we calculated the vertical transition from S to Tαα (STαα) and from Tαα to S (TααS), as well as the PES in which contribution of the unpaired antiparallel spin state (Tαβ) was accounted for. The Tαβ state corresponded to the existence of two doublet spin states, one in each of the two reacting molecules. A vertical TααS spin inversion at REACTOOH or a vertical STαα spin inversion at INTC4aOOH would be unlikely (ΔE > 20 kcal·mol−1). Instead, the lowest-energy path for the formation of the C4a···Op bond showed that a potential energy barrierless spin inversion could occur between the unpaired antiparallel spin state (Tαβ) and the closed-shell singlet state (S) through an intersystem crossing point (ICP), at a C4a···Op distance close to 2.16 Å (see Tαβ and S PESs in Figure 5A). The results are in line with those proposed in recent literature for C2:HPAH, which have calculated an ICP from the triplet to the closed-singlet state occurring at ca. 2.2 Å with an activation enthalpy of 1.4 −1.28 Furthermore, both the TααS and STαα energy profiles suggested that a TααS spin inversion could occur in between two energy crossing points (ECPS at 2.36 Å and ECPTαα at 2.76 Å, Figure 5A). We have optimized these ECPs, performing DFT calculations on a cluster model comprising the DFT layer of our QM/MM model, at the B3LYP/ 6-31G(d) level, using a code developed by Jeremy Harvey.89,90 A minimum energy crossing point (MECP)90,91 was located at a C4a···Op distance of 2.43 Å (Figure 5B), very close to the intersection of the S and Tαα PESs along the putative reaction coordinate. This MECP should correspond to an energy lower than 6 kcal·mol−1 (energy difference from the ECP with lowest 9304


Research Article

ACS Catalysis

Figure 6. DFT layer of the REACT, TS1, and INT1 states of the oxidation of DBT by C4aOOH, for the two poses of DBT modeled from HPA of C2:HPAH (top) and DBT of TdsC (bottom). The distances are indicated in Å. Blue-to-red shades in atoms represent the most relevant changes in atomic charge relative to the REACT state (from 0.07 to 0.30 au). Blue stands for increase in atomic charge, and red stands for decrease in atomic charge.

while the bottom TS1 showed to be slightly earlier and more asymmetric (Od was positioned at 1.91 Å from Op and 2.08 Å from the SDBT). Both TS1 exhibited an almost linear SDBT··· Od···Op angle (172.7° and 173.4°, respectively), typical for an SN2 reaction. The most favorable Gibbs activation energy was 19.7 kcal·mol−1, corresponding to the bottom DBT pose (DBTTdsC), in good agreement with the experimentally derived one (a turnover of 1.6 min−1, corresponding to an overall ratelimiting step of 20.4 kcal·mol−1).12 An analysis of the vibrational frequencies in TS1 provided a sole imaginary frequency (313.7−343.9i cm−1) that was dominated by the antisymmetric stretching of the Od−Op and Od−SDBT bonds. A NBO analysis suggested that the interaction between Op···SDBT was dominated by natural atomic orbitals (NAOs) with p character. All these observations again pointed out that the DBT oxidation to DBTO follows a SN2 mechanism. Upon comparison of the DBTHPAH and DBTTdsC poses, we also observed that the kinetics of DBT oxidation depended more on the initial orientation of SDBT relatively to Od (ΔG‡ of 29.3 −1 for an SDBT···Od···Op angle of 146.5°, and 19.7 −1 for an SDBT···Od···Op angle of 164.1°) than on the SDBT−Od distance (refer to section C in SI), which is also in line with previous results on SN2 substitutions.48 From REACT to TS1, we noted that the hydrogen bond between Od and the Nε of His391 seemed stronger for the DBTTdsC pose (0.06 Å increase) than for DBTHPAH (0.32 Å increase), again suggesting that the DBTTdsC pose was better placed for the nucleophilic attack of SDBT to Od; while the hydrogen bond between Op and the Tyr96-hydroxyl became shorter in about 0.58 and 0.31 Å, in agreement with the increase in electron density in Op (−0.08 to −0.24 au). An analysis of the atomic charges at the REACT and TS1 also showed that the transfer in electron density occurred mainly from the SDBT (δq ≈ + 0.20 au) to Op (δq ≈ − 0.15 au). Previous mutagenesis studies on DszC did not led to clear conclusions about the role of Tyr96: Zhang et al. indicated that mutation of this residue to Ala decreased DBT consumption

energy, ECPTαα), and it implied that the mechanism of spin inversion through this MECP was also more favorable than the vertical spin inversion at the reactant state (∼21 kcal·mol−1). Upon characterization of the triplet REACTOOH state and the singlet INTC4aOOH, we calculated a Gibbs reaction energy of 68.8 kcal·mol−1 for the mechanism of oxygen activation. We then repeated the same procedure having the DBT in the active site. In the presence of DBT, the energy gap between the singlet and triplet states was negligible, even though the relative stabilization of INT C4aOOH comparatively to REACTOOH was alike in both cases (ΔE ∼ 30 kcal·mol−1). As a result, the TS vertical transition mechanism seemed to be preferred over the MECP mechanism, but could still be competitive with the ICP mechanism (Figure S9 in SI). Because either of the discussed possibilities was equally fast, our results emphasized that it is more likely that DBT should enter the active site once C4aOOH is formed. Second Stage: Oxidation of the Dibenzothiophene. Step 4 in Scheme 1. An X-ray structure of a DszC ortholog (TdsC) was very recently obtained with FMN and DBT present in the active site. Therefore,33 we decided to characterize two reactant states, with DBT modeled from HPA in C2:HPAH structure (DBTHPAH), and from the DBT in the TdsC structure (DBTTdsC) (both shown in Figure 2). While DBTHPAH was bent 69.1° relatively to the SDBTOd axis, with the nucleophilic DBT-sulfur (SDBT) 3.76 Å away from the C4aOOH−hydroperoxyl Od atom (Figure 6, top), DBTTdsC was bent 105.7° along the same axis, with SDBT 3.09 Å away from the Od atom (Figure 6, bottom). In both structures, C4aOOH established a short hydrogen bond between Od and the Nε of His391 (1.89−1.92 Å), as well as with the Asn129carbamide and Ser163-hydroxyl, through the isoalloxazine moiety (see REACT in Figure 6). A transition state (TS1, Figure 6) for the nucleophilic attack of the DBT-sulfur (SDBT) to Od was optimized for both REACT states. The top TS1 exhibited a symmetric nature (Od was positioned at 2.00 Å from Op and 1.97 Å from the SDBT), 9305


Research Article

ACS Catalysis

Figure 7. DFT layer of the INT1, TS2, INT2, TS3, PROD1, and PROD2 stationary points for the regeneration of the FMN cofactor. The distances are indicated in Å. Blue-to-red shades in atoms represent the most relevant changes in atomic charge relative to the INT1 state (from 0.07 to 0.30 au). Blue stands for increase in atomic charge, and red stands for decrease in atomic charge.

to the Oγ of Ser163 (1.85 Å apart), concerted with the deprotonation of the Ser163-hydroxyl by the Nε of His92 (1.71 Å apart). There were no major geometric changes in the active site during this proton transfer. The transition state (TS2) was characterized by a frequency of 1021.2i cm−1, with the N5hydrogen at 1.32 Å from N5 and 1.20 Å from Oγ. The deprotonation of Ser163 by His92 was slightly asynchronous with the late proton transfer from N5. At TS2, His92 already received the proton from Ser163 (TS2, Figure 7). The Gibbs activation energy of the process was 12.6 kcal·mol−1, well within the range of barriers calculated from similar reactions in flavin-dependent oxidases.27,93−95 We could only determine a stationary state at the end of this step at 0 K (INT2, Figure 7); the Gibbs energy profile showed that this INT2 should be a transient state instead. After TS2, the C4a−Op bond increased (from 1.47 to1.97 Å), and the isoalloxazine group of the C4aOH intermediate adopted a more planar pose, characteristic of the oxidized form of FMN (TS3 in Figure 7): the N5−C4a single bond shortened (1.42 Å to 1.29 Å), indicating a pronounced double bond character, and the N5−C4a-C10a and C4a-C10a-N1 planes changed from 41° to 22°. An NBO analysis at TS3 confirmed the sp2 character of the natural atomic orbitals characterizing the N5−C4a bond, even though it was not a double bond. Upon TS3 characterization, we also observed that the hydrogen bond formed between the Ser163-hydroxyl and the O4 of the isoalloxazine ring (1.89 Å) promoted the exit of the OpH leaving group (reorientation of Ser163-hydroxyl from N5 to O4 decreased the energy of TS3 in ca. 5 kcal·mol−1).The analysis of the geometry changes, along with the atomic charges and overall charge of C4aOH (Table S4 in SI), suggested that the elongation of the C4a−Op bond occurred with concomitant sp2-rehybridization of N5 and C4a, which stabilized the increase in electron density at N5 resulting from the proton transfer to Ser163 (INT2 and PROD1, Figure 7, and Table S4 in SI). The sp2 character of the N5−C4a bond was also validated by the NBO analysis, which provided two NBOs for the bond: one with high sp2 character, and the other involving mostly p orbitals from N5 and C4a.

by 25−50%,22 while Liu et al. indicated that the same mutation abolished enzyme activity.21 Our results indicated that the hydrogen bond between Tyr96 and Op anchored the latter Op during the nucleophilic attack of SDBT, and prevented steric hindering by Tyr96 on DBT, which could compromise the DszC activity (refer to section B in SI). From TS1 to INT1, the Od−SDBT bond fully formed and Od underwent a late deprotonation by Op, forming DBTO and C4a-hydroxyflavin (C4aOH). The SDBT−Od bond length was 1.53 Å, and it was neither a characteristic single bond (1.66 Å) nor a double bond (1.47 Å). NBO analysis showed that this bond was best represented by a single bond with largest contributions from p NAOs, as opposed to the predicted double bond.12 His391 lost the hydrogen bond to the (former) OOH group of C4aOOH, while DBTO became hydrogen bonded to the OpH. There were no significant geometry changes from TS1 to INT1, other than the position of Od and Hd. The changes in atomic charges from REACT to INT1 were consistent with the oxidation of the DBT-sulfur (δq ≈ + 0.30 au) and the reduction of the Op (δq ≈ − 0.15 au) and Od (δq ≈ − 0.20 au) that took place during the reaction. The process was very exoergonic, with a Gibbs reaction energy of −32.6 kcal·mol−1. Third Stage: Formation of Oxidized Flavin. Steps 5 and 6 in Scheme 1. In the last step of the catalytic mechanism of DszC, the FMN cofactor was regenerated from the C4aOH intermediate. For this to occur, the C4a−OpH bond had to be cleaved. The results indicated that this cleavage is preceded by the deprotonation of N5 of C4aOH (Figure 7), because a preliminary PES for the direct cleavage of the C4a− OpH led OpH to bond to C10a, with a barrier above 25 kcal· mol−1. In agreement with previous mutagenesis studies,21,92 His92 was critical for this step. The hydrogen bond between Ser163 and His92 was critical for the deprotonation of the N5 in C4aOH, because it was only in the presence of this interaction that we observed a transition state for the reaction; in its absence, Ser163 was hydrogen bonded to the O4 of C4aOH and the reaction did not occur (refer to section B in SI). The most favorable pathway comprised a proton transfer from N5 9306


Research Article

ACS Catalysis TS3 exhibited one imaginary frequency (147.2i cm−1) corresponding to an antisymmetric stretching of the C4a−Op bond and the Op−H bond from the Tyr96-hydroxyl, a pronounced out-of-plane displacement of the N5−C4a bond in C4aOH, and a rotation of the Ser163-hydroxyl toward C4carbonyl of C4aOH. The Gibbs activation energy was 16.3 kcal·mol−1 relatively to the INT1 state (also in agreement with experimental studies performed in flavin-dependent oxidases that suggest barriers in between 13.5−16.2 kcal·mol−1 for this step);27,93−95 it was very close to the Gibbs energy of TS2 and INT2, indicating that the releasing of the OpH-hydroxyl was coupled with the proton transfer from N5 to the Ser163-hydroxyl. The atomic charges throughout INT1 to TS3 showed that there was an increase in electron density in FMN and a decrease in the His92, via the Ser163 (overall charge was similar in INT1, INT2 and TS3, Table S4 in SI). Since His92 was close to Tyr96, the positive charge built in His92 (qHis92 ≈ + 0.43) could support the asynchronous deprotonation of Tyr96-hydroxyl by OpH, after TS3, and stabilize the negatively charged Tyr96 (qTyr96 ≈ −0.65). The product of the reaction (PROD1) was oxidized FMN (FMNox) and a water molecule, which was hydrogen bonded to the Tyr96-phenoxide (1.63 Å) and to the DBTO-sulfoxide (2.11 Å). The Gibbs energy of reaction was +8.1 kcal·mol−1, meaning that the recycling of FMN would be an endoergonic process (PROD1 in Figure 8). The release of the water

This emphasized that restoring the uncharged state of Tyr96 should be key to prepare the active site of DszC to unbind the DBTO substrate. Gibbs Energy Profile for the Reaction Mechanism of DszC. Overall, the determined Gibbs energy profile is in very good agreement with the rate constant experimentally determined for DszC, which is ca. 1.6 min−1 (corresponding to an observed barrier of 20.4 kcal.mol−1).12 Figure 8 depicts the Gibbs energy profile at 298 K for the chemical reaction catalyzed by DszC. An analysis of the Gibbs energy profile for the reaction clearly defined two critical steps in this mechanism, the oxidation of DBT (ΔG‡ of 19.7 kcal· mol−1) and the dehydration/regeneration of FMN (ΔG‡ of 16.3 kcal·mol−1). The overall chemical step of DszC was very exergonic (−38.7 kcal·mol−1). In light of these results, we propose the oxidation of DBT to DBTO as the rate-limiting step of the chemical cycle of DszC. Zero-point energy (ZPE) and entropy corrections made a small but non-negligible contribution to the ΔG profile of the chemical reaction by DszC. The reaction is mostly enthalpydriven, with entropy contributions below 2 kcal·mol−1 for most steps of the reaction. In transition states involving proton transfer reactions (TS2 and TS4), ZPE corrections lowered the energy of the enzyme complex in ca. 3 kcal·mol−1 relatively to the closest minimum (INT1 and PROD1); while in TS3 and PROD1, the ZPE corrections from the release of OpH (TS3 and PROD1) lowered the complex’s energy in about 2 kcal· mol−1 relatively to the REACT state (Table S6 in SI). Regarding entropy contributions, in TS2, vibrational entropy contributions led to an increase in the Gibbs activation energy of 1.7 kcal·mol−1. This is most likely related to the stiffness of the active site, since from INT1 to INT2, there are mainly hydrogen bond rearrangements, and the more negatively charged C4aOH is likely to favor binding of the cofactor at the active site. On the other hand, upon protonation of Tyr96 by the positively charged His92, the vibrational entropy of the product state (PROD2) lowered the Gibbs energy in 1.3 kcal· mol−1 (relative to PROD1). This observation emphasizes what we described to be a more flexible active site in the absence of salt interactions between His92 and Tyr96. Our theoretical proposal of the complete reaction mechanism of DBT oxidation by DszC is depicted in Scheme 2. The main mechanistic aspects proposed by this study are highlighted in blue. The radial chart in light gray details the Gibbs activation energy for each step of the reaction mechanism relatively to the proposed rate-limiting step (19.7 kcal·mol−1, corresponding to the oxidation of DBT, refer to Figure 8).

Figure 8. Gibbs energy profile at the SCS-MP2/def2-TZVPP:ff10// B3LYP/6-31G(d):ff10 level of theory for the reaction mechanism catalyzed by DszC, discriminating the ΔG, ΔH, and −TΔS contributions at 298 K. ΔG and ΔH are shown on the left y-axis, and −TΔS is shown on the right y-axis. The most favorable Gibbs energy barrier for the rate-limiting step (TS1) from the DBTTdsC pose is emphasized over the less favorable Gibbs barrier.

■

molecule from the active site of DszC should be entropically favored, thus completing the first catalytic cycle of DszC (oxidation of DBT to DBTO). We also studied the protonation of Tyr96 by His92. In the reaction, the Tyr96phenoxide deprotonated the His92-imidazolium, overcoming a Gibbs activation energy of +13.0 kcal·mol−1. The characterization of the transition state resulted in a single imaginary frequency (1026.1i cm−1) dominated by the antisymmetric stretching of the O−H and Nδ−H bonds of Tyr96 and His92, with the Nδ-hydrogen 1.20 Å far from the Nδ-nitrogen of His92 and at 1.36 Å from the Tyr96-phenoxide. Upon the protonation of Tyr96 by His92, the water molecular byproduct ended up less tightly bonded to the active site, because it established a shorter hydrogen bond with the DBTO-sulfoxide. The product state (PROD2) was very exoergonic, with a Gibbs reaction energy of −16.6 kcal·mol−1 relatively to INT1.

CONCLUSIONS We studied the reaction mechanism of DBT oxidation by the flavin monoxygenase (FMO) DszC, which has been regarded as an important target for biodesulfurization of crude oil through the 4S pathway of Rhodococcus erythropolis. We provide detailed atomistic and thermodynamic results regarding both the mechanism of molecular oxygen activation common to many FMOs, and the chemistry of DBT oxidation, which is a major sulfur carrier in crude oil. The reaction mechanism by DszC proceeds in three stages: (1) molecular oxygen activation through proton-coupled electron transfer from FMNH2 to dioxygen, via His391, followed by a triplet-to-singlet spin inversion during the formation of the hydroperoxyflavin-intermediate (C4aOOH); 9307


Research Article

ACS Catalysis Scheme 2. Summary of Our Proposal for the Reaction Mechanism of DBT Oxidation by DszC

bond with the Ser163 that allows deprotonation of N5 in C4aOH and further promotes dehydration of C4aOH by Tyr96. Our insight on the molecular and mechanistic aspects of the catalysis by DszC and other FMOs is a significant step in the pursuit of the detailed molecular picture for this class of enzymes and provides a full assessment of their potential in rational drug design and bioengineering campaigns. In particular, the clear identification of key residues for the oxidation of DBT derivates further paves the way for the needed rational enzyme engineering of DszC, so that its optimization may become of technological use in biodesulfurization of crude oil.

(2) oxidation of DBT to DBTO through nucleophilic attack of the DBT-sulfur to the distal oxygen (Od) in C4aOOH; and (3) dehydration of FMN cofactor through a proton transfer from the N5 to the His92-imidazole via the Ser163-hydroxyl, and cleavage of the C4a−Op bond by Tyr96. The overall reaction is computed to be exoergonic, and entropy contributions do not play a major role in the chemical reaction by DszC. An in-depth study of the triplet-to-singlet inversion mechanism during oxygen activation by FMOs enforced that the proton-coupled electron transfer mechanism to form the hydroperoxyl radical (OOH·) involves the positively charged His391 (proton donor) and FMNH (electron donor). The resulting FMNH• and OOH• radicals covalently bond through a barrierless antiferromagnetic-coupled doublet state, or through a minimum energy crossing point mechanism (with a barrier lower than 6 kcal·mol−1) in which Op is 2.43 Å far from the C4a-carbon. DBT oxidation is proposed to be the rate-limiting step of the reaction and exhibits a Gibbs activation energy of 19.7 kcal· mol−1, which is in line with current available experimental data that estimate a barrier of ca. 20 kcal·mol−1. During DBT oxidation, a hydrogen bond between Tyr96 and the proximal oxygen (Op) of C4aOOH anchors it for the nucleophilic attack by DBT-sulfur and prevents steric impediment on DBT from the bulky side chain of Tyr96. Furthermore, we proposed that the regeneration of FMN occurs in a single step, in which His92 and Tyr96 have important direct roles. His92 forms a mandatory hydrogen

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acscatal.8b01877. Modeling of the DszC:FMN:DBT complex and details on MD simulations (section A); structural results regarding the protonation of FMN and His391 from MD simulations (section B); QM/MM results of the DBT oxidation for the different protonations of the His92-His388-His391 triad (section B); QM/MM results of the DBT oxidation for different models of the DszC:FMN:DBT complex (section C); PES for the mechanism of oxygen activation in the presence and absence of DBT (section D); relevant atomic charge variations for the reaction mechanism by DszC (section 9308


Research Article

ACS Catalysis

■

(2) Kilbane, J. J.; Le Borgne, S. Petroleum biorefining: the selective removal of sulfur, nitrogen, and metals. In Studies in Surface Science and Catalysis; Vazquez-Duhalt, R., Quintero-Ramirez, R., Eds.: Elsevier: Amsterdam, 2004; Vol. 151, pp 29−65. (3) Harrison, R. M.; Yin, J. X. Particulate matter in the atmosphere: which particle properties are important for its effects on health? Sci. Total Environ. 2000, 249, 85−101. (4) Pope, C. A., III; Burnett, R. T.; Thun, M. J.; Calle, E. E.; Krewski, D.; Ito, K.; Thurston, G. D. Lung cancer, cardiopulmonary mortality, and long-term exposure to fine particulate air pollution. JAMA, J. Am. Med. Assoc. 2002, 287, 1132−1141. (5) Boniek, D.; Figueiredo, D.; dos Santos, A. F. B.; de Resende Stoianoff, M. A. Biodesulfurization: a mini review about the immediate search for the future technology. Clean Technol. Environ. Policy 2015, 17, 29−37. (6) Bachmann, R. T.; Johnson, A. C.; Edyvean, R. G. J. Biotechnology in the petroleum industry: An overview. Int. Biodeterior. Biodegrad. 2014, 86, 225−237. (7) Kulkarni, P. S.; Afonso, C. A. M. Deep desulfurization of diesel fuel using ionic liquids: current status and future challenges. Green Chem. 2010, 12, 1139−1149. (8) Srivastava, V. C. An evaluation of desulfurization technologies for sulfur removal from liquid fuels. RSC Adv. 2012, 2, 759−783. (9) Shafi, R.; Hutchings, G. J. Hydrodesulfurization of hindered dibenzothiophenes: an overview. Catal. Today 2000, 59, 423−442. (10) Vazquez-Duhalt, R.; Torres, E.; Valderrama, B.; Le Borgne, S. Will biochemical catalysis impact the petroleum refining industry? Energy Fuels 2002, 16, 1239−1250. (11) Izumi, Y.; Ohshiro, T. Purification and characterization of enzymes involved in desulfurization of dibezothiophene in fossil fuels. J. Mol. Catal. B: Enzym. 2001, 11, 1061−1064. (12) Abin-Fuentes, A.; Mohamed, M. E.; Wang, D. I. C.; Prather, K. L. J. Exploring the Mechanism of Biocatalyst Inhibition in Microbial Desulfurization. Appl. Environ. Microbiol. 2013, 79, 7807−7817. (13) Gallagher, J. R.; Olson, E. S.; Stanley, D. C. Microbial desulfurization of dibenzothiophene: a sulfur-specific pathway. FEMS Microbiol. Lett. 1993, 107, 31−35. (14) Gray, K. A.; Mrachko, G. T.; Squires, C. H. Biodesulfurization of fossil fuels. Curr. Opin. Microbiol. 2003, 6, 229−235. (15) Gupta, N.; Roychoudhury, P. K.; Deb, J. K. Biotechnology of desulfurization of diesel: prospects and challenges. Appl. Microbiol. Biotechnol. 2005, 66, 356−366. (16) Denome, S. A.; Oldfield, C.; Nash, L. J.; Young, K. D. Characterization of the Desulfurization Genes from Rhodococcus Sp Strain Igts8. J. Bacteriol. 1994, 176, 6707−6716. (17) Gray, K. A.; Pogrebinsky, O. S.; Mrachko, G. T.; Xi, L.; Monticello, D. J.; Squires, C. H. Molecular mechanisms of biocatalytic desulfurization of fossil fuels. Nat. Biotechnol. 1996, 14, 1705−1709. (18) Kilbane, J. J. Microbial biocatalyst developments to upgrade fossil fuels. Curr. Opin. Biotechnol. 2006, 17, 305−314. (19) Le Borgne, S.; Quintero, R. Biotechnological processes for the refining of petroleum. Fuel Process. Technol. 2003, 81, 155−169. (20) Monticello, D. J. Biodesulfurization and the upgrading of petroleum distillates. Curr. Opin. Biotechnol. 2000, 11, 540−546. (21) Liu, S. H.; Zhang, C. G.; Su, T. T.; Wei, T. D.; Zhu, D. Y.; Wang, K.; Huang, Y.; Dong, Y. H.; Yin, K.; Xu, S. J.; Xu, P.; Gu, L. C. Crystal structure of DszC from Rhodococcus sp. XP at 1.79 angstrom. Proteins: Struct., Funct., Genet. 2014, 82, 1708−1720. (22) Zhang, L.; Duan, X. L.; Zhou, D. M.; Dong, Z.; Ji, K. H.; Meng, W. Y.; Li, G. Q.; Li, X.; Yang, H. T.; Ma, T.; Rao, Z. H. Structural insights into the stabilization of active, tetrameric DszC by its Cterminus. Proteins: Struct., Funct., Genet. 2014, 82, 2733−2743. (23) Guan, L. J.; Lee, W. C.; Wang, S.; Ohshiro, T.; Izumi, Y.; Ohtsuka, J.; Tanokura, M. Crystal structures of apo-DszC and FMNbound DszC from Rhodococcus erythropolis D-1. FEBS J. 2015, 282, 3126−3135. (24) Duan, X.; Zhang, L.; Zhou, D.; Ji, K.; Ma, T.; Shui, W.; Li, G.; Li, X. Crystallization and preliminary structural analysis of dibenzothiophene monooxygenase (DszC) from Rhodococcus

D); and thermochemistry results for the reaction mechanism by DszC (section D) (PDF) Gaussian input files for all stationary states of the reaction mechanism of DszC (ZIP) PDB coordinate files for all stationary states of the reaction mechanism of DszC (ZIP) Orca input files for the DFT and SCS-MP2 calculations for all stationary states of the reaction mechanism of DszC (ZIP) Animation 1 for the reaction mechanism of oxygen activation by DszC (MPG) Animation 2 for entire reaction mechanism of DBT oxidation by DszC (MPG) Molecular dynamics parameters for the DszC:FMN and DszC:FMN:DBT complexes (ZIP) Gaussian input files for stationary states in section B; PDB coordinate files for stationary states in section B; Gaussian input files for stationary states in section C; PDB coordinate files for stationary states in section C (ZIP, ZIP)

AUTHOR INFORMATION

Corresponding Author

*E-mail for P.A.F.: [email protected]. ORCID

Rui P. P. Neves: 0000-0003-2032-9308 Maria J. Ramos: 0000-0002-7554-8324 Present Address

§ UCIBIO@REQUIMTE, BioSIM, Departamento de Biomedicina, Faculdade de Medicina da Universidade do Porto, Alameda Professor Hernâni Monteiro, 4200-319, Porto, Portugal

Author Contributions ‡

(A.C.C.B. and R.P.P.N.) These authors contributed equally.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge financial support from European Union (FEDER funds POCI/01/0145/FEDER/007728) and National Funds (Fundaçaõ para a Ciência e Tecnologia and Ministério da Educaçaõ e Ciência − FCT/MEC) under the Partnership Agreement PT2020, UID/MULTI/04378/2013; NORTE-01-0145-FEDER-000024, supported by Norte Portugal Regional Operational Programme (NORTE 2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF); Fundaçaõ para a Ciência e a Tecnologia (FCT) through projects EXCL/ QEQ-COM/0394/2012 and PTDC/QUI-QFI/28714/2017, IF/00052/2014.

■

ABBREVIATIONS FMN, flavin mononucleotide; DBT, dibenzotiophene; ONIOM, our own N-layered integrated molecular orbital and molecular mechanics;; QM/MM, quantum mechanics/ molecular mechanics

■

REFERENCES

(1) Monticello, D. J.; Finnerty, W. R. Microbial Desulfurization of Fossil-Fuels. Annu. Rev. Microbiol. 1985, 39, 371−389. 9309


Research Article

ACS Catalysis erythropolis. Acta Crystallogr., Sect. F: Struct. Biol. Cryst. Commun. 2013, 69, 597−601. (25) Lei, B. F.; Tu, S. C. Gene overexpression, purification, and identification of a desulfurization enzyme from Rhodococcus sp strain IGTS8 as a sulfide/sulfoxide monooxygenase. J. Bacteriol. 1996, 178, 5699−5705. (26) Huijbers, M. M.; Montersino, S.; Westphal, A. H.; Tischler, D.; van Berkel, W. J. Flavin dependent monooxygenases. Arch. Biochem. Biophys. 2014, 544, 2−17. (27) Wongnate, T.; Surawatanawong, P.; Visitsatthawong, S.; Sucharitakul, J.; Scrutton, N. S.; Chaiyen, P. Proton-coupled electron transfer and adduct configuration are important for C4a-hydroperoxyflavin formation and stabilization in a flavoenzyme. J. Am. Chem. Soc. 2014, 136, 241−253. (28) Visitsatthawong, S.; Chenprakhon, P.; Chaiyen, P.; Surawatanawong, P. Mechanism of Oxygen Activation in a FlavinDependent Monooxygenase: A Nearly Barrier less Formation of C4aHydroperoxyflavin via Proton-Coupled Electron Transfer. J. Am. Chem. Soc. 2015, 137, 9363−9374. (29) Bruice, T. C. Oxygen-Flavin Chemistry. Isr. J. Chem. 1984, 24, 54−61. (30) Massey, V. Activation of Molecular-Oxygen by Flavins and Flavoproteins. J. Biol. Chem. 1994, 269, 22459−22462. (31) Alfieri, A.; Fersini, F.; Ruangchan, N.; Prongjit, M.; Chaiyen, P.; Mattevi, A. Structure of the monooxygenase component of a twocomponent flavoprotein monooxygenase. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 1177−1182. (32) Thotsaporn, K.; Chenprakhon, P.; Sucharitakul, J.; Mattevi, A.; Chaiyen, P. Stabilization of C4a-Hydroperoxyflavin in a Twocomponent Flavin-dependent Monooxygenase Is Achieved through Interactions at Flavin N5 and C4a Atoms. J. Biol. Chem. 2011, 286, 28170−28180. (33) Hino, T.; Hamamoto, H.; Suzuki, H.; Yagi, H.; Ohshiro, T.; Nagano, S. Crystal structures of TdsC, a dibenzothiophene monooxygenase from the thermophile Paenibacillus sp A11−2, reveal potential for expanding its substrate selectivity. J. Biol. Chem. 2017, 292, 15804−15813. (34) Zgarbova, M.; Otyepka, M.; Sponer, J.; Mladek, A.; Banas, P.; Cheatham, T. E.; Jurecka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7, 2886−2902. (35) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157−1174. (36) Schneider, C.; Suhnel, J. A molecular dynamics simulation of the flavin mononucleotide-RNA aptamer complex. Biopolymers 1999, 50, 287−302. (37) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269−10280. (38) Sondergaard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H. Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values. J. Chem. Theory Comput. 2011, 7, 2284−2295. (39) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (40) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene,

M.; Knox, J. E.; 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.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (41) Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Goetz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M. J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. Amber 12; University of California: San Francisco, 2012. (42) Swiderek, K.; Marti, S.; Tunon, I.; Moliner, V.; Bertran, J. Peptide Bond Formation Mechanism Catalyzed by Ribosome. J. Am. Chem. Soc. 2015, 137, 12024−12034. (43) Himo, F. Recent Trends in Quantum Chemical Modeling of Enzymatic Reactions. J. Am. Chem. Soc. 2017, 139, 6780−6786. (44) van Keulen, S. C.; Rothlisberger, U. Exploring the inhibition mechanism of adenylyl cyclase type 5 by n-terminal myristoylated G alpha(i1). PLoS Comput. Biol. 2017, 13, e1005673. (45) Dubey, K. D.; Wang, B. J.; Vajpai, M.; Shaik, S. MD simulations and QM/MM calculations show that single-site mutations of cytochrome P450(BM3) alter the active site’s complexity and the chemoselectivity of oxidation without changing the active species. Chem. Sci. 2017, 8, 5335−5344. (46) Prejanò, M.; Marino, T.; Rizzuto, C.; Madrid Madrid, J. C.; Russo, N.; Toscano, M. Reaction Mechanism of Low-Spin Iron(III)and Cobalt(III)-Containing Nitrile Hydratases: A Quantum Mechanics Investigation. Inorg. Chem. 2017, 56, 13390−13400. (47) Jindal, G.; Ramachandran, B.; Bora, R. P.; Warshel, A. Exploring the Development of Ground-State Destabilization and Transition-State Stabilization in Two Directed Evolution Paths of Kemp Eliminases. ACS Catal. 2017, 7, 3301−3305. (48) Neves, R. P. P.; Fernandes, P. A.; Ramos, M. J. Mechanistic insights on the reduction of glutathione disulfide by protein disulfide isomerase. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, E4724−E4733. (49) Neves, R. P. P.; Fernandes, P. A.; Ramos, M. J. Unveiling the Catalytic Mechanism of NADP+-Dependent Isocitrate Dehydrogenase with QM/MM Calculations. ACS Catal. 2016, 6, 357−368. (50) Sousa, S. F.; Ribeiro, A. J. M.; Neves, R. P. P.; Bras, N. F.; Cerqueira, N. M. F. S. A.; Fernandes, P. A.; Ramos, M. J. Application of quantum mechanics/molecular mechanics methods in the study of enzymatic reaction mechanisms. WIREs Comput. Mol. Sci. 2017, 7, e1281. (51) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic-Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (52) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (53) Becke, A. D. Density-functional thermochemistry. IV. A new dynamical correlation functional and implications for exact-exchange mixing. J. Chem. Phys. 1996, 104, 1040−1046. (54) Hariharan, Pc; Pople, J. A. Influence of Polarization Functions on Molecular-Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213−222. (55) Hariharan, P. C.; Pople, J. A. Accuracy of Ah Equilibrium Geometries by Single Determinant Molecular-Orbital Theory. Mol. Phys. 1974, 27, 209−214. (56) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular-Orbital Methods 0.12. Further Extensions of GaussianType Basis Sets for Use in Molecular-Orbital Studies of OrganicMolecules. J. Chem. Phys. 1972, 56, 2257−2261. (57) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. 9310


Research Article

ACS Catalysis (58) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Results Obtained with the Correlation-Energy Density Functionals of Becke and Lee, Yang and Parr. Chem. Phys. Lett. 1989, 157, 200−206. (59) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin-Density Calculations - a Critical Analysis. Can. J. Phys. 1980, 58, 1200−1211. (60) Calixto, A. R.; Brás, N. F.; Fernandes, P. A.; Ramos, M. J. Reaction Mechanism of Human Renin Studied by Quantum Mechanics/Molecular Mechanics (QM/MM) Calculations. ACS Catal. 2014, 4, 3869−3876. (61) Ribeiro, A. J. M.; Yang, L.; Ramos, M. J.; Fernandes, P. A.; Liang, Z.-X.; Hirao, H. Insight into Enzymatic Nitrile Reduction: QM/MM Study of the Catalytic Mechanism of QueF Nitrile Reductase. ACS Catal. 2015, 5, 3740−3751. (62) Gesto, D. S.; Cerqueira, N. M. F. S. A.; Fernandes, P. A.; Ramos, M. J. Unraveling the Enigmatic Mechanism of L-Asparaginase II with QM/QM Calculations. J. Am. Chem. Soc. 2013, 135, 7146− 7158. (63) Seeger, R.; Pople, J. A. Self-Consistent Molecular-Orbital Methods 0.18. Constraints and Stability in Hartree-Fock Theory. J. Chem. Phys. 1977, 66, 3045−3050. (64) Bauernschmitt, R.; Ahlrichs, R. Stability analysis for solutions of the closed shell Kohn-Sham equation. J. Chem. Phys. 1996, 104, 9047−9052. (65) Schlegel, H. B.; McDouall, J. J. Do You Have SCF Stability and Convergence Problems? In Computational Advances in Organic Chemistry; Ö gretir, C., Csizmadia, I. G., Eds.; Kluwer Academic: The Netherlands, 1991; Vol. 330, pp 167−185. (66) Maseras, F.; Morokuma, K. Imomm - a New Integrated AbInitio Plus Molecular Mechanics Geometry Optimization Scheme of Equilibrium Structures and Transition-States. J. Comput. Chem. 1995, 16, 1170−1179. (67) Singh, U. C.; Kollman, P. A. A Combined Ab initio QuantumMechanical and Molecular Mechanical Method for Carrying out Simulations on Complex Molecular-Systems - Applications to the CH3Cl + Cl− Exchange-Reaction and Gas-Phase Protonation of Polyethers. J. Comput. Chem. 1986, 7, 718−730. (68) Field, M. J.; Bash, P. A.; Karplus, M. A Combined QuantumMechanical and Molecular Mechanical Potential for MolecularDynamics Simulations. J. Comput. Chem. 1990, 11, 700−733. (69) Das, D.; Eurenius, K. P.; Billings, E. M.; Sherwood, P.; Chatfield, D. C.; Hodoscek, M.; Brooks, B. R. Optimization of quantum mechanical molecular mechanical partitioning schemes: Gaussian delocalization of molecular mechanical charges and the double link atom method. J. Chem. Phys. 2002, 117, 10534−10547. (70) Chai, J.-D.; Head-Gordon, M. Systematic optimization of longrange corrected hybrid density functionals. J. Chem. Phys. 2008, 128, 084106. (71) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem. - Eur. J. 2012, 18, 9955−9964. (72) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. Efficient Diffuse Function-Augmented Basis Sets for Anion Calculations. Iii. The 3-21+G Basis Set for First-Row Elements, Li-F. J. Comput. Chem. 1983, 4, 294−301. (73) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-Consistent Molecular-Orbital Methods 0.25. Supplementary Functions for Gaussian-Basis Sets. J. Chem. Phys. 1984, 80, 3265−3269. (74) Mclean, A. D.; Chandler, G. S. Contracted Gaussian-Basis Sets for Molecular Calculations 0.1. 2nd Row Atoms, Z = 11−18. J. Chem. Phys. 1980, 72, 5639−5648. (75) Zhao, Y.; Truhlar, D. G. Hybrid Meta Density Functional Theory Methods for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions: The MPW1B95 and MPWB1K Models and Comparative Assessments for Hydrogen Bonding and van der Waals Interactions. J. Phys. Chem. A 2004, 108, 6908−6918. (76) Zhao, Y.; Truhlar, D. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new

functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (77) Boese, A. D.; Martin, J. M. L. Development of density functionals for thermochemical kinetics. J. Chem. Phys. 2004, 121, 3405−3416. (78) Neese, F. The ORCA program system. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (79) Neese, F. Software update: the ORCA program system, version 4.0. WIREs Comput. Mol. Sci. 2018, 8, e1327. (80) Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (81) Grimme, S. Improved second-order M?ller-Plesset perturbation theory by separate scaling of parallel- and antiparallel-spin pair correlation energies. J. Chem. Phys. 2003, 118, 9095−9102. (82) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (83) Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (84) Kossmann, S.; Neese, F. Efficient Structure Optimization with Second-Order Many-Body Perturbation Theory: The RIJCOSXMP2Method. J. Chem. Theory Comput. 2010, 6, 2325−2338. (85) Neese, F.; Olbrich, G. Efficient use of the resolution of the identity approximation in time-dependent density functional calculations with hybrid density functionals. Chem. Phys. Lett. 2002, 362, 170−178. (86) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge-Densities. Theor. Chim. Acta 1977, 44, 129−138. (87) Ritchie, J. P. Electron-Density Distribution Analysis for Nitromethane, Nitromethide, and Nitramide. J. Am. Chem. Soc. 1985, 107, 1829−1837. (88) Ritchie, J. P.; Bachrach, S. M. Some Methods and Applications of Electron-Density Distribution Analysis. J. Comput. Chem. 1987, 8, 499−509. (89) Harvey, J. N. Understanding the kinetics of spin-forbidden chemical reactions. Phys. Chem. Chem. Phys. 2007, 9, 331−343. (90) Harvey, J. N. Spin-forbidden reactions: computational insight into mechanisms and kinetics. WIREs Comput. Mol. Sci. 2014, 4, 1− 14. (91) Lykhin, A. O.; Kaliakin, D. S.; Depolo, G. E.; Kuzubov, A. A.; Varganov, S. A. Nonadiabatic transition state theory: Application to intersystem crossings in the active sites of metal-sulfur proteins. Int. J. Quantum Chem. 2016, 116, 750−761. (92) Zhang, X.; Tang, Y. Y.; Qu, S. Q.; Da, J. W.; Hao, Z. P. H2SSelective Catalytic Oxidation: Catalysts and Processes. ACS Catal. 2015, 5, 1053−1067. (93) Witt, S.; Wohlfahrt, G.; Schomburg, D.; Hecht, H. J.; Kalisz, H. M. Conserved arginine-516 of Penicillium amagasakiense glucose oxidase is essential for the efficient binding of beta-D-glucose. Biochem. J. 2000, 347, 553−559. (94) Roth, J. P.; Klinman, J. P. Catalysis of electron transfer during activation of O-2 by the flavoprotein glucose oxidase. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 62−67. (95) Sucharitakul, J.; Wongnate, T.; Chaiyen, P. Hydrogen Peroxide Elimination from C4a-hydroperoxyflavin in a Flavoprotein Oxidase Occurs through a Single Proton Transfer from Flavin N5 to a Peroxide Leaving Group. J. Biol. Chem. 2011, 286, 16900−16909.

9311


Mechanistic Studies of a Flavin Monooxygenase: Sulfur Oxidation of

Recommend Documents