Resolving the hydride transfer pathway in oxidative conversion of proline to pyrrole
Thiotemplated pyrrole is a prevailing intermediate in the synthesis of numerous natural products where the pyrrole is tethered to a carrier protein (CP). Biosynthesis of the pyrrole requires oxidation of an L-proline side chain. Herein, we investigated the biocatalytic mechanism of proline-to-pyrrole synthesis using the recently reported (Thapa et al., Biochemistry, 2019, 58, 918) structure of a Type II non-ribosomal protein synthetase (NRPS) Bmp3-Bmp1 (Oxidase-CP) complex. The substrate (L-proline side chain) is attached to the Bmp1(CP) and the catalytic site is located inside the flavin-dependent oxidase (Bmp3). Interestingly, the FAD molecule is free (unbound) within the Bmp3 catalytic site. We show that the FAD isoalloxazine ring is stabilized in the catalytic site of Bmp3 by strong hydrogen bonding with Asn123, Ile125, Ser126, and Thr158. The stability of the dimeric Bmp3 system including one FAD molecule in our simulation suggests that the tetrameric Bmp3 assembly, as found in the crystal structure, is not functionally important. After the initial deprotonation followed by an enamine-imine tautomerization, oxidation of either the C2-C3 bond or the C2-N1 bond, through a hydride transfer (either from C3 or N1), is required for the pyrrole synthesis. Using molecular dynamics simulations, quantum mechanics/molecular mechanics simulations, and electronic structure calculations, we conclude that the hydride transfer is more likely to occur from C3 than N1. Additionally, we demonstrated the elasticity in the oxidase active site through enzymatic synthesis of proline derivatives.
Introduction<!>Model description and the substrate bound state<!>Similarity to acyl-CoA dehydrogenases and prior mutational analysis suggests that the<!>Hydride transfer in the CS-II intermediate<!>Experimentally querying active site flexibility<!>Transition state and the free-energy barrier of hydride transfer<!>Summary and Conclusions<!>System preparation and MD simulations<!>QM/MM simulations<!>QM calculations in a cluster model<!>Conflicts of interest
<p>The pyrrole heterocycle is widespread in biological chemistry. In primary metabolism, pyrroles serve to bind metal ions in porphyrins. The metal binding capacity of pyrroles is exploited by bacterial secondary metabolites as well. 1 Pyrroles intercalate DNA and serve as a core synthon around which many pharmaceutical molecules are constructed. 2 Pyrroles in several secondary metabolite biosynthetic schemes are constructed via a type II non-ribosomal peptide synthetase (NRPS) adenylation-oxidation cascade which involves the 4e − oxidation of the pyrrolidine side chain of L-proline (Fig. 1a). 3 Prior to oxidation, the L-proline carboxylate is activated by adenylation, catalyzed by an adenyltransferase, which then transfers L-proline to a carrier protein (CP) phosphopantetheinyl thiol. [4][5][6][7] The phosphopantetheinyl thiol is derived from coenzyme A (CoA) and is appended post translationally to a serine side chain hydroxyl of the apo-CP to create the holo-CP. 8 An oxidase then catalyzes the conversion of the pyrrolidine to pyrrole in a sequence of two 2e − oxidations with electrons funneled to molecular oxygen via a flavin cofactor. 9,10 The thiotemplated pyrrole thus developed can be delivered to type I polyketide synthases for chain elongation, 11 to type I NRPSs for ligation to other amino acids, 12 and to other modification enzymes such as halogenases 13 or methyltransferases 14 for further modification.</p><p>The biosynthesis of thiotemplated pyrroles using the abovementioned type II NRPS scheme relies on intramolecular protein-protein interactions, of the CP with the adenyltransferase and of the CP with the oxidase, as well as catalysis in the respective adenyltransferase and oxidase active sites. Biochemically validated atomistic detail for both sets of protein-protein interactions are now available. 5,10 Catalysis in the adenyltransferase active site follows the well accepted two-step process involving the ATP-mediated activation of the amino acid carboxylate followed by thioesterification of proline to the CP phosphopantetheinyl thiol. Conversion of this thiotemplated proline-derived pyrrolidine to pyrrole via a pyrroline intermediate 9 in the oxidase active site remains less well resolved. The structure and overall mechanism of pyrroli(di)ne oxidases resembles the mitochondrial fatty acid acyl-CoA dehydrogenases (Fig. 1b). Catalysis in fatty acid dehydrogenases proceeds by proton abstraction from the fatty acyl C α by a catalytic glutamate residue followed by hydride transfer from the C β to the flavin cofactor leading to the C α -C β oxidation. 15,16 Pyrrolidine and pyrroline thioesters are unstable which has prevented a structural description of the substrate-or the intermediate-bound states of the pyrroli(di)ne oxidases. 10,17,18 Furthermore, a lack of insight into substrate and intermediate binding in the oxidase active site has precluded deciphering the site of hydride transfer from the pyrroli(di)ne to the flavin cofactor. Two mechanistic routes are conceivable, which converge to identical outcomes with interspersed imine-enamine tautomerizations (Fig. 1c). 19 First is that the enzyme acts as an amine oxidase which would involve hydride transfer from prolyl N1 to flavin N5, as has been proposed previously. 17,20 The second possibility is that the pyrroli(di)ne oxidase acts as a C α -C β desaturase with the hydride being transferred from the prolyl C3 to flavin N5. At present, these two routes are indistinguishable.</p><p>Herein, we investigate the first two steps of pyrrole biosynthesis using molecular dynamics (MD), quantum mechanics/molecular mechanics (QM/MM) simulations, and electronic structure calculations. In particular, we assess the feasibility of the Glu242 residue to initiate the biosynthesis. More importantly, we provide mechanistic insights into the hydride transfer from the substrate to the FAD bound to the Bmp3 enzyme. We investigate the energy barriers for the unique hydride transfer process from the C3 atom of the substrate instead of from the N1 atom. Furthermore, by enzymatic synthesis of proline derivatives, we demonstrate elasticity in the oxidase active site, which supports the computational description of substrate binding.</p><!><p>Prior genetic and enzymological description of a thiotemplated pyrrole biosynthetic pathway from marine bacteria had identified the L-proline-derived pyrrolidine/pyrrolyl carrier protein (CP) to be localized to the N-terminus of the Bmp1 polypeptide (henceforth referred to as Bmp1(CP)) and the pyrrolidine oxidase as the Bmp3 polypeptide. 7 Subsequent crystallographic description of the Bmp1(CP)/Bmp3 complex identified a tetrameric Bmp3 core similar to that of acyl-CoA dehydrogenases. 10 Each Bmp3 monomer binds one Bmp1(CP), thus organizing a 4:4 Bmp1(CP)/Bmp3 stoichiometric complex. Each Bmp1(CP) contacts two Bmp3 monomers in this tetrameric complex. One of these Bmp3 monomers is the one to whose active site the phosphopantetheinylated molecular cargo is delivered to by Bmp1(CP).</p><p>The second Bmp3 monomer only makes contact with the Bmp1(CP) without participating in catalysis. For computational analysis of catalysis occurring in the Bmp3 active site, our model system comprises a Bmp1(CP) peptide in complex with two Bmp3 peptides with the flavin cofactor preserved in the catalytically relevant Bmp3 monomer only. The model henceforth used for computation is illustrated in Figure 2a.</p><p>Structure of the pyrrolidine oxidase Bmp3 in complex with Bmp1(CP) was described in the product bound state. 10 The substrate bound state, referred to as catalytic state-I (CS-I, Figure 1) was inaccessible. Other crystal structures of pyrrolidine oxidases were similarly not available with bound thiotemplated substrates. 18,21 Hence, we first sought to recover the physiological substrate-bound CS-I state for Bmp3. To describe the CS-I state, we equilibrated the system where the thioesterified pyrrolidine heterocycle was modelled in the Bmp3 active site adjacent to the cofactor FAD isoalloxazine, followed by extended MD simulations. The initial position of the substrate in the CS-I state was set as found in the structure of the pyrrole in the product-bound state. We analyzed two replicas of 500-ns MD simulations of the CS-I model. For the occupancy analysis, we pooled data from both replicas.</p><p>In this CS-I model, we observed that the FAD isoalloxazine ring is stabilized by hydrogen bonding with Asn123, Ile125, Ser126, and Thr158 (Figures S2a-e). The percentage of all hydrogen bonding interactions between FAD, Bmp3, and water is provided in Table S1 and binding interactions of the FAD are illustrated in Figure S1. The pyrrolidine ring forms a hydrogen bond with the Asn123 side chain (Figure S2a), while the carbonyl group forms hydrogen bonds with the Ser363 backbone (Figure S2b,c). The substrate forms hydrogen bonds with Bmp3 predominantly through the Ppant moiety (Figure S2c-i), which attaches the pyrrolidine ring with the Bmp1(CP). In summary, the FAD isoalloxazine ring and the substrate form several hydrogen bonds with the Bmp3 protein to stabilize the catalytic site. Glu242 side chain carboxylate acts as the catalytic base for proton abstraction from the pyrrolidine side chain to initiate catalysis (Figure 1). We, 10 and others 18 had postulated that Glu242 abstracts the pyrrolidine C2-proton ( pyrrolidine C2-H). Therefore, the distance between the Glu242 side chain carboxylate and pyrrolidine C2 and the pyrrolidine C2-H atoms is critical for catalysis. We calculated these distances along the MD trajectories (Figure S3).</p><!><p>We observed that the pyrrolidine C2 and pyrrolidine C2-H atoms stay close to the Glu242 side chain in MD simulations (Figure S3b). Since rotation is allowed around the pyrrolidine C2-carbonyl carbon bond, the pyrrolidine ring can adopt two dominant conformations. To quantify the population distribution, we analyzed the differences between pyrrolidine C2−Glu242 and pyrrolidine C2-H−Glu242 distances as ∆d = d1−d2 (Figure 2b). One conformation brings the pyrrolidine C2-H proximal to the Glu242 carboxylate while the other conformation places the pyrrolidine C2-H distal to the Glu242 (Figure 2c). Since the two conformations can interconvert by rotation around a single bond, the rotational energy barrier should be small. This rotation does not change the pyrrolidine C2−Glu242 distance. If the C2-H atom is pointed towards Glu242, then ∆d>0. In contrast, ∆d<0 when the C2-H atom is pointed away from the Glu242 sidechain.</p><p>The combined data from two independent simulations show that the relative abundance of two rotamers are comparable (Figure 2d). The catalytically proficient conformation is represented by ∆d>0, since the deprotonated carboxylate of Glu242 side chain is geometrically close to the pyrrolidine C2-H atom. We observed that the catalytically proficient conformation was favored in the second simulation, while the catalytically deficient conformation was favored in the first simulation (Figure 2c,d).</p><p>We performed additional QM/MM simulations on the CS-I state to obtain more accurate dynamics of the catalytic site. Details of the QM/MM simulations are described in the Methods. From the QM/MM simulation, we observed that the most probable pyrrolidine C2−Glu242 distance is larger than the pyrrolidine C2-H−Glu242 distance, while both substrate atoms remain close to the Glu242 side chain (also see Figure S5a,b). Note that we do not observe the catalytically deficient ∆d<0 state from QM/MM simulation of the CS-I state, possibly because of the shorter time-scale (15 ps) of the QM/MM simulation. Within that time scale, QM/MM simulation only sampled conformations conducive to the proton transfer from the substrate to the Glu242 side chain (Figure S5a,b). While the substrate can bind in the active site in either conformation, catalysis can only proceed in one conformation and not in the other. Based on these results, we posit that the ∆d>0 state illustrated in Figure 2c is the catalytically proficient substrate bound state for pyrrolidine oxidase Bmp3. Mutations of Glu242 in Bmp3 and the corresponding residue Glu244 in pyrrolidine oxidase AnaB render the enzymes inactive. 10,17 Taken together, these results indicate the feasibility of the Glu242 protonation by proton transfer from the C2 atom of the substrate. More importantly, the active site is open enough for free rotation around the C-C2 bond allowing the pyrrolidine ring to fully rotate.</p><!><p>After thiotemplated substrate binding and abstraction of the pyrrolidine C2-H by the Glu242 side chain carboxylate, the catalytic state-II (CS-II, Figure 1) is defined by the protonated Glu242 side chain (Figure 1c). Pyrrolidine oxidation can proceed by tautomerization followed by hydride transfer from the pyrrolidine C3 to the flavin N5, or via hydride transfer from pyrrolidine N1 to flavin N5 (Figure 1c). Both of these routes have been proposed in literature, but evidence supporting either of the two routes has been lacking as the enolate intermediate was not accessible experimentally. To address this knowledge gap, we performed six independent 20-25 ps QM/MM simulations for the CS-II intermediate state using the same QM/MM scheme adopted for the CS-I simulations. In the CS-I state, all carbon and nitrogen atoms of the pyrrolidine ring are sp 3 hybridized. In contrast, the pyrrolidine C2 and pyrrolidine C3 atoms (for C2-C3 oxidation) of the substrate are in sp 2 hybridization states (Figure 1). In the case of plausible C2-N1 oxidation, C2 and N1 will become sp 2 hybridized. Both mechanisms have been proposed before in the flavin-catalyzed biosynthesis processes. 9,17,20 We investigated which mechanism is more likely to occur using a total of 145 ps of QM/MM simulations along with transition state calculations.</p><p>We analyzed the pyrrolidine C3− flavin N5 and pyrrolidine N1− flavin N5 distances to identify the closer hydride donor between pyrrolidine N1 and flavin C3 (Figure 3a). We observed that the most probable pyrrolidine C3− flavin N5 distance was shorter (∼ 3.0 Å) than the pyrrolidine N1− flavin N5 (∼ 4.2 Å) (Figure 3a). Additionally, we calculated the distance of the flavin N5 atom from the pyrrolidine C3-H and pyrrolidine N1-H as illustrated in Figure 3b. The most probable pyrrolidine C3-H− flavin N5 distances were shorter than the pyrrolidine N1-H− flavin N5 (1.0−2.0 Å and 4.0−4.5 Å, respectively) as shown in Figure 3b. Therefore, pyrrolidine C3 and pyrrolidine C3-H atoms are closer to the FAD cofactor as compared to pyrrolidine N1 and pyrrolidine N1-H, respectively. These results allow us to posit that the enolate intermediate tautomerizes such that the hydride transfer to the flavin cofactor occurs from pyrrolidine C3 and not from pyrrolidine N1. This then sets up the CS-III state in which the pyrrolidine C2−C3 bond should be unsaturated (Figure 1).</p><p>In support of the aforementioned assertion, in one of the QM/MM simulations, we captured the hydride transfer event from the CS-II pyrrolidine C3 to the flavin N5 (Figures 3c,d).</p><p>Therefore, the pyrrolidine C3-H− flavin N5 distance distribution has two peaks for one of the simulations as illustrated in Figure 3b. We further queried the structural perturbation accompanying the hydride transfer in this simulation. We observed that the pyrrolidine C2−C3 bond length reduced from 1.50 Å to 1.37 Å, while the carbonyl C− pyrrolidine C2 bond length increased from 1.42 Å to 1.50 Å with the hydride transfer from the C3 atom (Figure 3e).</p><p>We optimized a model substrate in the CS-II and the CS-III states in a low-dielectric solvent ( = 6.2528) using the conductor-like polarizable continuum model (CPCM) implicit solvent model to obtain ideal bond lengths for pyrrolidine C2−C3 and pyrrolidine C−C2. The optimization was performed at the CAM-B3LYP/6-31+G(d) level of theory. The optimized pyrrolidine C2−C3 and carbonyl C− pyrrolidine C2 bond lengths are 1.51 Å and 1.35 Å in the CS-II state, respectively, and 1.34 Å and 1.35 Å in the CS-III state, respectively (Figure S7). The changes in these bond lengths match the probable bond lengths obtained from the QM/MM simulation of the states before and after the hydride transfer event. In contrast, pyrrolidine C2−N1 and carbonyl C−O distances remain similar before and after the hydride transfer event. Therefore, we observed only one peak in the most probable pyrrolidine C2−N1 distance (Figure S6c), indicating that pyrrolidine N1 did not participate in this hydride transfer event. Taken together, these results indicate that pyrrolidine C2−C3 oxidation during the hydride transfer is more likely to occur compared to pyrrolidine C2−N1 oxidation.</p><!><p>Data presented above points towards the Bmp3 active site accommodating a large degree of flexibility in substrate binding with catalytically proficient and catalytically non-proficient modes of substrate binding being energetically equivalent. To further explore substrate tolerance of the Bmp3 active site, we enzymatically generated the pyrrolidine C4-Me diastereomers and explored if they could be processed in the Bmp3 active site. Using previously described enzymatic synthesis schemes, [22][23][24] L-leucine oxidative cyclization by the LdoA/NosE/NosF enzymes furnished the (2S,4S)4-methylproline (1, Figure 4) and the GriE/GriF/ProC enzymes furnished the (2S,4S)4-methylproline (2) diastereomer. The 4-methylproline isomers were purified from preparative scale enzyme reactions and confirmed by liquid chromatography/mass spectrometry after derivatization of the carboxylate (Figures S12a,b). Molecules 1 and 2 were used as substrates in a coupled adenyltransfer-pyrrolidine oxidation reaction in which the adenyltransferase Bmp4 transferred the proline derivatives to the phosphopantetheinyl thiol of holo-Bmp1(CP) in vitro followed by oxidation by wild type Bmp3.  As monitored by the formation of the characteristic cyclo-acylpantetheine fragmentation ejection ion, 25 full conversion of holo-Bmp1(CP) to 4-methylpyrrole-Bmp1(CP) was achieved in this reaction. While the activity of the adenyltransferase Bmp4 would indeed be compromised as compared to activation of the native substrate L-proline, this result demonstrates elasticity in the Bmp3 active site, such that it can successfully incorporate diverse substrates.</p><!><p>Although QM/MM simulations indicate the feasibility of the C2-C3 oxidation, they do not provide any information about the energy required for the hydride transfer event. To obtain additional insight into the energy required for hydride transfer, we performed transition state calculations using a cluster model. We discovered that the CS-II and CS-III intermediates are separated by 16.6 kcal/mol, with CS-III being lower in energy. Additionally, we observed that while the FAD isoalloxazine ring is planar in the CS-II intermediate, upon hydride transfer, the isoalloxazine ring bends out-of-plane around the short axis of the ring, leading to a reduction in conjugation in the CS-III intermediate (Figure 5). The bending of the isoalloxazine ring has been observed in prior studies 26 when FAD is reduced to FADH − . Furthermore, we observed that the HOMO and LUMO in the CS-II state are located on the FAD and the substrate (Figure 5), with a slightly higher contribution from the substrate.</p><p>The overlap of the FAD and the substrate orbitals originates from the π-π * overlap between them. However, the loss of planarity of the isoalloxazine ring in the CS-III intermediate leads to the abrogation of the π-π * overlap. Hence, the HOMO and LUMO of the CS-III state are located only on FAD and the substrate, respectively. The localization of the HOMO on FAD in the CS-III intermediate confirms that the negative charge of the hydride is located on the FAD molecule in the CS-III intermediate.</p><p>The transition state (TS) is located 6.7 kcal/mol above the CS-II intermediate state.</p><p>The TS is associated with a single imaginary frequency centered on the C3-H atoms, which captures the hydride transfer from C3 to N5 as indicated by the red arrow on the TS (Figure 5). Therefore, we conclude that the free energy barrier for the hydride transfer is 6.7 kcal/mol followed by structural rearrangements, which include slight bending of the isoalloxazine ring.</p><p>In line with abovementioned findings, here, we additionally detect that the direct hydride transfer from N1 is geometrically constrained, and a few structural rearrangements must occur before the transfer. For example, the substrate must rotate around the C−C2 bond (see Figures 1 and S7 for definitions) to come closer to the flavin N5. A similar rotation around the C−C2 bond was observed in the CS-I state (Figures 2c,d). In contrast to the CS-I state, C−C2 is a double bond in the CS-II state. Therefore, the rotational energy barrier is expected to be much higher in this case. We computed a potential energy surface (PES) the O-C-C2-N1 dihedral angle (Figure S8) in a CPCM solvent model using a low dielectric-constant ( = 6.2528) to mimic the protein environment. The geometry of the molecule was optimized along the dihedral angle at 10 • intervals starting from the optimized CS-II state geometry of the substrate. Consistent with other calculations using the cluster model, the relaxed PES scan was performed at the CAM-B3LYP/6-31+G(d) level of theory.</p><p>We found that the aforementioned rotation has a 33.6 kcal/mol energy barrier, which is much higher than the free-energy barrier for FAD reduction by hydride transfer from the C3 atom. Although a complete rotation may not be required for C2-N1 oxidation, even a ∼70 • rotation is ∼18 kcal/mol higher in energy than the initial state (Figure S8). To estimate the free energy of activation for this rotation, we also performed a TS calculation at the same level of theory. The TS is confirmed by a single imaginary frequency, which represents the rotation around the dihedral angle. We observed that the free energy barrier for this rotation is 33.5 kcal/mol while the two end-states differ by only 0.5 kcal/mol. The free-energy barrier of the rotation is ∼5-fold higher than the hydride transfer from C3, which does not require additional structural rearrangements. Therefore, we conclude that the C2-C3 oxidation of the substrate is more likely to occur compared to C2-N1 oxidation during hydride transfer to the FAD.</p><!><p>We have simulated the initial state in the pyrrole biosynthesis in which the substrate protein is bound to an active and an inactive Bmp3-enzyme. We demonstrated that the tetrameric Bmp1-Bmp3 complex is not functionally important for the enzyme catalysis. One Bmp1-Bmp3 can perform the chemical reaction in the presence of a second inactive Bmp3. Furthermore, we demonstrated that the active site of Bmp3 can stabilize FAD and proline substrate. We showed that the Glu242 can initiate the biosynthesis process by capturing a proton from the substrate using MD and QM/MM simulations. Additionally, our results show a significant population of the substrate conformations that favors proton uptake by Glu242. Our simulations show that the C3 atom remains closer to the FAD hydride acceptor site than the N1 atom and one QM/MM simulation captured the hydride transfer event from C3 to FAD within 12 ps without any applied bias. We also determined that hydride transfer to FAD is coupled with C2-C3 oxidation of the substrate and must overcome a very small (6.7 kcal/mol) energy barrier. In contrast, the alternative mechanism, via C2-N1 oxidation, requires rotation around the C=C bond, which has a 33.5 kcal/mol activation free energy barrier. Our findings provide mechanistic insights for a reaction that is well represented in biological chemistry, but, has evaded experimental characterization due to instability of the substrate and intermediate species.</p><!><p>We used the tetrameric complex structure (PDB ID: 6CY8) of the holo Bmp3-Bmp1 system at 2.7-Å resolution. 10 To test the stability and the functional relevance of the Bmp3-Bmp1 complex, we prepared a system containing the substrate, a single bound Bmp1 (chain A), and a dimer of Bmp3 (chains B and C), where only chain B contains the catalytic FAD molecule (Figure 1A). Therefore, we call chains B and C active and inactive Bmp3, respectively.</p><p>Bmp3 structures have eight missing residues (residues 24 to 31), which were modeled using SWISS model. 27 The Bmp1 structure is missing five N-terminal residues and 1 C-terminal residue, which were added using the VMD PSFgen plugin. 28 We solvated this complex in an explicit TIP3P water box 29 with at least 12 Å of solvent padding around the protein. The system was neutralized by adding 28 Na + ions. We used the CHARMM36m force field 30 for protein residues. Initial force-field parameters for the FAD molecule and the substrate were taken from CGenFF, 31 which assigns them by analogy to existing CHARMM parameters.</p><p>Charges on the isoalloxazine ring of FAD and the pyrrolidine ring of the substrate (in the CS-I state) were further improved using the Force Field Toolkit (ffTK), as implemented in VMD. 32,33 Parameters for the dihedral angles, where the pyrrolidine rings connect with the phosphopantetheine appendage, were also improved using ffTK (see Supporting Information for details). MD simulations were performed in the NPT ensemble using periodic boundary conditions at a constant temperature of 310 K and a constant pressure of 1 atm. Long-range electrostatics were computed using the Particle-mesh Ewald (PME) method. 34 We used a 12-Å cutoff with a switching function between 10-12 Å. All MD simulations were performed using a 2-fs time step.</p><p>We first minimized the modeled missing pieces followed by minimization of the entire system. We equilibrated the solvation box for 0.5 ns with the proteins restrained. After that, we equilibrated the system for another 0.5 ns with a 2 kcal/mol/ Å2 restraint on protein backbone atoms, followed by a 50-ns equilibration without any restraints. The final structure from this equilibration was used to initiate all QM/MM simulations. Finally, we ran two independent production simulations for 500 ns, from which we collected trajectory data every 20 ps for subsequent analysis.</p><!><p>The protein system was extracted from the 50-ns equilibration in the CS-I state. Since the active-site charge distribution changes significantly during the reaction, each system was resolvated and neutralized before QM/MM simulations. The QM region included all residues directly involved in the reaction. The water molecule next to FAD is also included in the QM region to investigate its role in each step of the reaction. The resultant QM region included: FAD, substrate, Glu242, and one water molecule (Figure 6a). The rest of the system was treated as point charges. The link-atom scheme was used to satisfy the valency of the terminal QM atom, where a hydrogen atom was added as the link atom. The redistributed charge and dipole (RCD) method was used to modify charges and dipoles in the QM-MM boundary region to avoid over-polarization. Additionally, we used the switch function to modify point charges and enforced net-zero electrostatic interactions between charges at the cutoff distance. The net-charge of the system was also set to zero. The QM part of the QM/MM calculation was performed at BP86/def2-SVP level of theory using ORCA. 35 The same level of theory was successfully used to model enzymatic hydride transfer from flavin mononucleotide (FMN) using QM/MM methodology previously. 36 Becke-Johnson damping (D3BJ) was used for dispersion corrections to the DFT. 37 QM/MM simulations were performed using the NAMD/ORCA interface 38 and a 0.5-fs time step. After minimization, we raised the temperature of the system to 310 K followed by equilibration for 100 fs with backbone restraints.</p><p>Finally, we ran one production simulation for the CS-II state using the NPT ensemble at 310 K and 1 atm. The same procedure was followed for QM/MM simulations of the CS-I and CS-II states. The difference between the two states was included in the QM part of the QM/MM simulation.</p><!><p>The cluster model was constructed using the QM/MM equilibrated structure of the CS-II system. In addition to residues already included in the QM region, we added Ser363, Asn123, Thr158, and Phe85 residues, which are located close to the active site (Figure 6b). Phe85 and Thr158 can play important roles since the hydrogen bond of FAD with Thr158 and hydrophobicity near the active site are required to tune the catalysis. 39 In total, the cluster model included five residues, substrate, FAD, and one water molecule. All C α -atom positions were frozen during subsequent QM optimizations to preserve the overall structure of the residues. Optimization of CS-II and CS-III states and the TS were performed at the CAM-B3LYP/6-31+G(d) with D3BJ dispersion corrections 37 using Gaussian16. 40 While the B3LYP functional has been used before to describe electronic states of FAD, 26 usage of additional long-range corrections is usually recommended for molecules with fused rings. 41 Therefore, we used the CAM-B3LYP functional, along with an implicit dispersion correction, to correctly describe non-bonded interactions within the cluster model. The average dielectric constant value inside a protein is usually in the 6-7 range. 42 While some studies have estimated a slightly lower dielectric constant (3-4) inside a protein, 43,44 it is unlikely that implicit solvents with a dielectric constant in 3-7 range will show discernible differences in molecular properties. Thus, we mimic the protein environment by using the CPCM implicit solvent model 45 with a low dielectric constant ( = 6.2528).</p><!><p>There are no conflicts to declare.</p>