Effect of Stacking Interactions on the Thermodynamics and Kinetics of

Article pubs.acs.org/JPCA

Effect of Stacking Interactions on the Thermodynamics and Kinetics of Lumiflavin: A Study with Improved Density Functionals and Density Functional Tight-Binding Protocol Caitlin G. Bresnahan,† Clorice R. Reinhardt,† Thomas G. Bartholow, John P. Rumpel, Michael North, and Sudeep Bhattacharyya* Department of Chemistry, University of WisconsinEau Claire, Eau Claire, Wisconsin 54702, United States S Supporting Information *

ABSTRACT: The π−π stacking interaction between lumiflavin and a number of π-electron-rich molecules has been studied by density functional theory using several new-generation density functionals. Six known lumiflavin-aromatic adducts were used and the models were evaluated by comparing the geometry and energetics with experimental results. The study found that dispersion-corrected and hybrid functionals with larger (>50%) Hartree−Fock exchanges produced superior results in modeling thermodynamic characteristics of these complexes. The functional producing the best energetics for these model systems was used to study the stacking interactions of lumiflavin with biologically relevant aromatic groups. Additionally, the reduction of flavinin the presence of both a hydride donor and a nondonor π-electronic system was also studied. Weak interactions were observed in the stacked lumiflavin complexes of benzene, phenol, and indole, mimicking phenyl alanine, tryptophan, and tyrosine side chains, respectively, of an enzyme. The stacked complex of naphthalene and flavin showed little change in flavin’s redox potential indicating insignificant effect on the thermodynamics of the hydride transfer reaction. In contrast, the hydride transfer reaction with the hydride donor N-methyl nicotinamide tells a different story, as the transition state was found to be strongly impacted by the stacking interactions. A comparison of performance between the density functional theory (DFT) and the computationally less expensive dispersion-corrected self-consistent density functional tight-binding (SCC-DFTB-D) theory revealed that the latter produces consistent energetics for this hydride transfer reaction and additional DFT-computed perturbative corrections could significantly improve these results. the flavin ring causes significant electronic change producing butterfly bending.20,21 These π-stacking interactions are also important for catalysis in enzymes such as quinone oxidoreductases (QRs).3,5,6,10,22 For example, the enzyme QR2 (being studied in our laboratory23,24) catalyzes hydride transfer reactions of quinones amidst the formation of stacked donor−acceptor complexes.6,8,23−25 The active site of QR uses a two-substrate-based ping-pong kinetics that involves substrates shuttling and catalytic hydride transfer steps. These processes are mediated by certain alterations in weak forces (van der Waals interactions) between the flavin ring and the substrate (or the product).4,5 These weak forces fall in the category of medium-long-range (3−5 Å and >5 Å) interactions and are difficult to assess26−28 because of theoretical limitations to treat electron correlations. Theoretical studies in recent years have helped to understand redox21,23,24,29,30 and photochemical31−35 reactions involving flavin. Recent efforts in our lab provided insights into the effects

1. INTRODUCTION Flavin is the primary electron−proton mediatory group in many enzyme families.1−5 These tricyclic heteronuclear isoalloxazine ring systems occur primarily as two cofactors: flavin mono nucleotide and flavin adenine dinucleotide. Occurring predominantly as noncovalently attached to the enzyme matrix,3,6−9 they catalyze a wide variety of chemical reactions.1,10,11 The catalytic diversity of flavin arises in part due to its ability to undergo coupled electron and proton transfers,12 which can occur either by two sequential 1e−/ 1H+ additions or a single-step 2e−/2H+ transfer. Additionally, the versatility is also shaped by flavin’s redox potential, which is modulated through its noncovalent interactions with active site amino acid residues. These interactions include hydrogen bonding with flavin ring atoms,13 electrostatic effects14 with polar side chains, and π−π stacking interactions4,5,15 with hydrophobic aromatic substrates and residues. The flavin ring possesses a large delocalized π-electron density,4 which plays a key role in the binding of aromatic substrates.5,16,17 Several recent studies suggest that van der Waals interactions play crucial role in modulating the optical/ electronic properties of flavin.18,19 It is known that reduction of © XXXX American Chemical Society

Received: October 3, 2014 Revised: December 2, 2014

A

dx.doi.org/10.1021/jp510020v | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

of ring substitution,21 hydrogen bonding, and active site electrostatics23,24 on flavin’s geometry and energetics. Nevertheless, the effects of aromatic interactions on the catalytic redox changes of flavin was only partially revealed in these studies. In the present study, we have investigated the flavin chemistry with density functional theory (DFT) keeping our focus on the stacking interactions. This is the first report of a detailed DFT-based study aiming to unravel the impact of π−π stacking on the thermodynamics and kinetics of flavin’s redox reactions. Although, the popular Becke’s 3-parameter Lee−Yang− Parr36 (B3LYP) functional has been used to study the oxidation−reduction processes of flavin,13,37 it is unsuitable to study redox events dominated by dispersion/van der Waals interactions.38 Therefore, in the present study, we have used improved density functionals to probe the geometry and thermodynamics of lumiflavin (hereafter will be abbreviated as LF) adducts of six planar polycyclic aromatic molecules (Figure 1). In particular, we examined the performances of several

continuum solvation using the Polarizable Continuum Model (PCM)57 within Gaussian 09.50 In all these calculations, the harmonic oscillator-rigid rotor approximation was used to obtain frequencies and moments of inertia. Hybrid quantum/classical simulations were carried out using CHARMM program suite51 with the CHARMM all-atom force field.58 The three-point-charge TIP3P model59 was used to treat water molecules. Nonbonded interactions were truncated using a switching function between 11 and 12 Å. Bond lengths and bond angles of water molecules were constrained by the SHAKE algorithm.60 In all molecular dynamics simulations, a time step of 1 fs was used in the leapfrog Verlet algorithm for integration.61,62 All necessary conditions for simulations including the partitioning scheme were maintained as described in the previous studies.23,24 Briefly, the stacked LF-nicotinamide complex was treated with SCC-DFTB-D,52,53 while embedded in a 30 Å thick layer consisting of molecular mechanically treated solvent atoms. Stochastic boundary conditions63 were employed throughout to simulate solvent effects. 2.1. Energetics of the Binding of Flavin-Aromatic Stacked Complexes. The standard-state free energy for the association of LF-aromatic complexes was determined in Scheme 1

Figure 1. Six polycyclic aromatic compounds (A−F) used in the binding study with flavins are coded as (A) 1,2-benzanthrancene, (B) 1,2-benzofluorene, (C) triphenylene, (D) fluoranthene, (E) 2,3benzofluorene, and (F) phenanthrene.

ethanol following Scheme 1. Representing the aromatic molecule (listed in Figure 1) as X

functionals built with the second derivative of the electron density (i.e., kinetic energy density) in addition to the local density and its first derivative (gradient), under the metageneral gradient approximation (meta-GGA). Four meta-GGA density functionals of varied formulations have been used: pure meta-GGA (TPSSTPSS-D339), meta-GGA hybrids with a constant amount of Hartree−Fock (HF) exchange (mPW1PW91-D3,40 M06, M06-2X, and M06-HF41,42), rangeseparated hybrid functionals without short-range HF exchange (LC-wPBE43), and range-separated hybrid functionals with short-range HF exchange (wB97X-D44). In addition, selfconsistent charge-density functional tight binding with dispersion corrections (SCC-DFTB-D),45,46 a less computerintensive approximate DFT formalism, was also used as a comparison. This has previously shown great promise especially in its capabilities in modeling proton transfer in enzymes47−49 as well as redox changes of flavoenzymes23,24,37,47 including QR223,24 through the use of a hybrid quantum mechanical/ molecular mechanical setup. Therefore, the study was extended to probe the reduction and hydride transfer reactions of flavin.

ΔGo(ethanol, LF ··· X) = ΔGo(g , LF ··· X) + ΔGS(LF ··· X) − ΔGS(LF) − ΔGS(X)

(1)

where ΔG (g, LF···X) is the gas-phase binding free energy of the LF−aromatic complex and ΔGS(LF···X), ΔGS(LF), and ΔGS(X), are the solvation free energies of the optimized complex, free LF, and the free aromatic molecule, respectively, in ethanol. For the study of lumiflavin complexes of three biologically ring systems (benzene, indole, and phenol), an implicit solvation model with a dielectric constant of 4.2 was used64 to approximate the hydrophobic environment of a protein’s active site. 2.2. Reduction of Flavin−Aromatic Complexes. The aqueous-phase free energy of the hydride addition reaction to the LF stacked with an aromatic nonhydride donor molecule was computed using the thermodynamic diagram of Scheme 2 o

ΔGo(aq) = ΔGo(g ) + ΔGS(LFH‐ ··· X) − ΔGS(LF ··· X)

2. METHODOLOGY Electronic structure calculations were carried out with the Gaussian 0950 and CHARMM51 computer programs. The structure and energy of the flavin and flavin bound aromatic systems were calculated using two different levels of theory: SCC-DFTB52,53 and DFT.54,55 The 6-31+G(d,p) basis set56 was used in all DFT calculations. All calculations were carried out at 298 K. Solvation effects were modeled implicitly with

− ΔGS(H+)

(2)

where ΔG (g) is the gas-phase free energy to transfer a hydride onto the N5 atom of the LF-aromatic complex (Scheme 3a) and ΔGS(LFH−···X) and ΔGS(LF···X) are the solvation free energies of the optimized hydride transferred product and the LF-aromatic complex, respectively, in water. In these condensed phase calculations, the electron has been conceived o

B



Article

C---H and N---H, for breaking and formation of bonds, respectively (Figure 2). On the basis of the optimized structure

Scheme 2

Scheme 3

Figure 2. Representation of the reaction coordinate defined in the free energy calculation for the reaction of lumiflavin and N-methyl nicotinamide using SCC-DFTB-D/MM.

of the reactant and product systems, the range of the reaction coordinate (zi) was chosen to be from −1.5 to +1.5 Å. The complete range of the reaction coordinate was divided into 11 segments of 0.4 Å called “windows”. In each of these windows, independent simulations were run by applying a harmonic biasing potential. This potential represents a restraining force centered at the midpoint (zoi ) of that particular window. The biasing potential was chosen based on a trial-and-error method so that the populations generated in a specific window produce a Gaussian distribution around zoi . In each window, the system was equilibrated initially for 50 ps, and then configurations along z were sampled for 500 ps to yield the PMF. Quantum Mechanical Corrections to the Classical Free Energy Surface. Two corrections were applied to the SCCDFTB-D-computed energetics: zero-point corrections (ΔZPE) and a high-level correction (ΔHLC) to account for the difference of thermodynamics of association between SCCDFTB-D and DFT.

of as an isolated particle in the gas-phase. As free energy is a state function, the present treatment is adequate and helps to avoid calculations of multiple reaction steps involving solvated electrons, which are highly reactive and can form radical-anion complex intermediates.65 ΔGoS(H+) is the free energy of solvation of an isolated proton and is equal to −264.0 kcal/ mol.66 The reaction in Scheme 3a can also be described as a halfreaction, where a hydride (2e−/1H+) is added onto LF, in the presence of the stacking interactions due to a nonhydride donor π-electronic system. An example is given in Scheme 3a, where the nonhydride donor aromatic molecule is naphthalene. The change in the Gibbs’ free energy, ΔGo(aq), can be used to compute the standard reduction potentials:

o o ΔGcorr = ΔGuncorr + ΔZPE + ΔHLC

ΔGo(aq) − EHo (3) nF where n is the number of electrons, F is the Faraday constant (23.06 kcal/mol), and EoH is the normal hydrogen electrode,67 whose reported value is 4.28 V.66 2.3. Hydride Transfer between Flavin and an Aromatic Donor. Using implicit solvation and DFT, the electron transfer reaction within the stacked complex of Nmethyl nicotinamide (NMH) and LF was studied (Scheme 3b). In parallel, a comparative study of the same reaction was carried out using hybrid SCC-DFTB-D/MM potentials with explicit solvation treatment. The activation barrier for the aqueous-state hydride transfer reaction was calculated independently using both implicit and explicit solvation techniques. Intrinsic Reaction Coordinate (IRC) Calculation. In the implicit solvation treatment, DFT-based electronic structure calculations were carried out for reactant and product states. The activation barrier was computed by following the IRC, which is defined as the minimum energy reaction pathway in mass-weighted Cartesian coordinates between the reactant and the product. Potential of Mean Force (PMF) Calculation. The activation energies for the hydride transfer reaction in the aqueous-state was calculated using explicit solvation treatment. MD simulations were performed using hybrid SCC-DFTB-D/MM potentials and the free energy change was obtained by plotting the potentials of mean force68 (PMF) as a function of the reaction coordinate. The reaction coordinate is defined as z = rbroken − rformed, where r represents the internuclear distances, Eo = −

(4)

where ΔZPE represents the difference of M06-2X-computed zero-point corrections between the transition state and the reactant for computing ΔG#(aq) or the zero-point corrections between the product and the reactant for computing ΔGo(aq). The perturbative correction, HLC for a particular state refers to the gas-phase Born−Oppenheimer energy differences between SCC-DFTB-D and DFT (M06-2X). In order to obtain the ΔHLC, the difference of the two perturbative corrections were calculated: one, between the transition state and the reactant for correcting ΔG#(aq) and the other, between the product and the reactant for correcting the ΔGo(aq). ΔHLC was computed using the average of 10 conformations for each of the reactant, the transition state, and the product. These conformations were chosen from the product/reaction coordinate corresponding to these states: within 0.01 Å of zrc for the reactant state, 0.005 Å of zts for the transition state, and 0.01 Å of zpr for the product state.

3. RESULTS AND DISCUSSION The presentation of results is done in three parts illustrating the three phases of this study. In the first phase, we studied the stacking interactions involving flavin and several known aromatic polycyclic molecules using various functionals. Next, we explored the energetics of reduction (hydride addition) of flavin stacked with a nonhydride donor aromatic molecules. In the third and final phase, the intracomplex hydride transfer in a model stacked complex (LF and N-methyl nicotinamide) was examined. C



Article

Table 1. Inter-Planar Distances (Å) Observed in the Optimized Adduct Complexes of Lumiflavin and the Polycyclic Aromatic Compounds shown in Figure 1 Using Various Density Functionals stacked complexes

TPSSTPSS-D3

mPW1PW91-D3

LC-wPBE

M06

M06-2X

M06-HF

wB97X-D

SCC-DFTB-D

LF--A LF--B LF--C LF--D LF--E LF--F average MUEa

3.42 3.46 3.45 3.39 3.39 3.38 3.42 0.06

3.30 3.29 3.26 3.24 3.25 3.30 3.27 0.14

3.82 3.81 3.61 3.60 3.80 3.53 3.70 0.30

3.58 3.41 3.49 3.41 3.41 3.53 3.47 0.08

3.26 3.26 3.29 3.25 3.28 3.35 3.28 0.13

3.22 3.22 3.20 3.30 3.34 3.27 3.26 0.15

3.41 3.39 3.35 3.47 3.35 3.45 3.40 0.05

3.21 3.19 3.15 3.19 3.30 3.24 3.21 0.20

a

The lowest row indicates the mean unsigned error (MUE) over all six complexes calculated from the reference LF--naphthalene complex, where the interplanar distance is 3.41 Å69

Table 2. Comparison of the Computed Gibbs Binding Free Energies (ΔG°comp(s)) in kcal/mol for the Lumiflavin−Aromatic adducts (LFX, where X = A − F in Figure 1) Using Various Density Functionals as Well as SCC-DFTB-D and Their Experimentally Determined Analogues19 for the Riboflavin Adducts, ΔG°expt(s) stacked complexes

TPSSTPSS-D3

mPW1PW91-D3

M06

M06-2X

M06-HF

wB97XD

SCC-DFTB-D

exptl19

LF--A LF--B LF--C LF--D LF--E LF--F

−2.4 −2.3 −2.8 −1.6 −1.9 −0.4

−11.0 −10.4 −9.6 −9.4 −9.5 −8.2

3.1 4.3 2.8 2.5 1.9 3.7

−2.1 −1.5 −1.3 −0.6 −0.5 −0.1

−1.0 −1.4 −2.0 −1.7 −2.3 −1.5

−5.2 −3.1 −3.3 −3.1 −3.7 −2.2

−2.1 −1.7 −2.0 −1.8 −1.7 −1.7

−2.6 −2.4 −2.3 −2.1 −2.2 −2.0

3.1. Flavin−Aromatic Stacked Complexes. LF forms aromatic stacked complexes indicating the presence of attractive van der Waals interactions between aromatic molecules and the flavin ring. In the present study, three classes of functionals have been used to probe these noncovalent interactions. Three dispersion-corrected functionals, namely, TPSSTPSS-D3,39 mPW1PW91-D3,40 and wB97X-D44 and one long-range corrected functional LCwPBE43 were used. Separately, Minnesota functionals M06, M06-2X, and M06-HF41,42 and as a comparison, the computationally less-expensive dispersion-corrected SCC-DFTB-D45 were also used. 3.1.1. Geometry. The interplanar separations of the gasphase optimized LF adduct complexes of the six polycyclic molecules (Figure 1) are listed in Table 1. The structure of the LF−naphthalene complex, whose X-ray structure has been previously reported,69 was used as a reference for comparison. The interplanar separation between the aromatic rings was overestimated in LC-wPBE (by 0.3 Å) when compared to the experimentally observed value of 3.41 Å69 (Table 1). For the remaining 6 functionals, namely, TPSSTPSS-D3, mPW1PW91D3, M06, M06-2X, M06-HF, and wB97X-D, the computed interplanar distances are within 0.1 Å of the experimental value indicating that they were able to correctly reproduce the geometry. The semiempirical SCC-DFTB-D calculations produced geometric results similar to M06 group of functionals (Table 1). 3.1.2. Energetics. As illustrated in Scheme 1, the binding free energies of six LF−aromatic adducts in ethanol were computed using SCC-DFTB-D and six functionals, namely, TPSSTPSSD3, mPW1PW91-D3, M06, M06-2X, M06-HF, and wB97X-D (Table 2). In the case of SCC-DFTB-D, the Gibbs’ free energies of binding were obtained from the SCC-DFTB-D gasphase potential energy by including the M06-2X-computed thermal, zero-point corrections, and solvation energies. All computed Gibbs’ energies, ΔG°comp(ethanol) and the exper-

imental binding free energies (ΔG°expt(ethanol) with riboflavin, last column) are listed in Table 2. The positive ΔG°comp(ethanol) values computed using the M06 functional indicates repulsion between the aromatic ring systems. In contrast, the ΔG°comp(ethanol) values were underestimated in mPW1PW91-D3 calculations by a significant amount (6−8 kcal/mol). As shown in Table 2, the Gibbs’ binding free energies computed by SCC-DFTB-D, TPSSTPSS-D3, M06-2X, M06-HF, and wB97X-D produced the best results comparable to the experiment. In order to compare the trend of the computed results, the computed binding free energies (ΔGocomp(ethanol,LF···X)) were plotted against known experimental values (ΔGoexpt(ethanol,LF···X, Figure 3). The SCCDFTB-D and M06-HF-computed energetics shows no apparent correlation with the experimental results. However, the computed energetics displayed a somewhat consistent trend for mPW1PW91-D3, TPSSTPSS-D3, M06-2X, and wB97X-D (Figure 3). The linearity of the plotted data suggests that these functionals were able to distinguish the subtle changes in the stacking interactions as the linear trend shows a R2 value ranging between 0.8 and 0.9. 3.1.3. Functional Dependencies. The π−π stacking interactions involve dispersion forces between the delocalized π-clouds over the two ring systems: the flavin ring and the aromatic molecule (Figure 1).70 Dispersion forces arise due to instantaneous multipole-induced multipole interactions necessitating a correct way to compute the motion of electrons, which are correlated. In normal Hartree−Fock treatment, the electron’s motion is derived using the average field created by other electrons and hence completely misses the local correlation effects. With Kohn−Sham DFT, the electron correlation is approximated by exchange-correlation functionals that are primarily based on local density approximations (LDA) under the uniform electron gas model. These local density functionals fail to appreciate medium-to-long-range forces (>3 Å) and are quite unreliable. D



Article

o Figure 3. A plot of Gibbs’ free energy of binding, ΔG comp (ethanol,LF···X), for the stacked lumiflavin complexes computed by M06-2X, wB97X-D, TPSSTPSS-D3, and mPW1PW91-D3 vs the comparable values (ΔGocomp (ethanol,LF···X)) obtained in experiment as given in Table 2. The x-coordinate for each data corresponds to a particular lumiflavin-aromatic complex as indicated at the top.

Figure 4. Potential energy surface of the complex LF--A in the direction perpendicular to the flavin ring as computed by various functionals.

(Figure 5) reveal low correlation coefficients for the variation of energies computed by TPSSTPSS-D3 and mPW1PW91-D3. In

In order to correctly model the medium-to-long-range (3−8 Å) behavior, several approaches have been taken in the past. First, the use of the generalized gradient approximation (GGA) by adding the gradient (first derivative) of exchange-correlation density and later the meta-GGA functionals, by adding the Laplacian (second derivative of the density) as corrections. Second, instead of pure local density a hybrid density is computed utilizing a portion of the nonlocal Hartree−Fock (HF) exchange admixed to a semi local density functional exchange. Third, in some functionals, the long-range behavior of the functional is attenuated by a range-separator parameter, γ, which defines the linear range beyond which the density functional approximation will be replaced with a 100% HF exchange. Fourth, following Grimme’s hybrid approach, classically computed dispersion energies are introduced into the optimization after the wave function of the molecule has been optimized and then use the aforementioned optimized energy for reoptimizing the molecular geometry. As observed in subsections 3.1.1 and 3.1.2, the two functionals that produced the most reliable geometry and energetics (Tables 1 and 2) are M06-2X and wB97X-D. The M06-2X functional possesses a constant amount (54%) of HF exchange. In contrast, wB97X-D has Becke’s 97 density functional with 22% short change HF exchange, a dispersion correction, but a 100% long-range HF exchange effected by a lower range-separator parameter, γ = 0.2 bohr−1. 3.1.4. The Potential Energy Well. In order to study how various functionals treat the stacking interactions, we computed potential energy surface (Figure 4) for a representative complex (with 1,2 benzantracene, LF---A, Figure 1). A large depth of 30 kcal/mol was obtained for mPW1PW91-D3, while M06 calculations produced a shallow well of only 16 kcal/mol. The depth of the well computed by TPSSTPSS-D3, M06-2X, wB97X-D, and SCC-DFTB-D are similar differing by only ∼2 kcal/mol (Figure 4). The decaying tail of the potential energies (in the direction perpendicular to the ring) were plotted against 1/r6, where r is the interplanar separation; the type of functional dependency is characteristic of dispersion forces. The fitted lines in the plot

Figure 5. Decaying tail of the potential energies for the LF---A plotted against 1/r6 (r = interplanar distance, in the range of 4.4−7.9 Å) using various theoretical protocol: (a) SCC-DFTB-D, (b) TPSSTPS-D3, (c) M06-2X, (d) wb97X-D, (e) mPW1PW91-D3, and (f) M06.

contrast, SCC-DFTB-D- and wB97X-D-computed potential energies displayed moderate correlation with 1/r6, whereas excellent correlations were obtained for M06 group of functionals (Figure 5). As a control study, the potential energy surface of a stacked benzene dimer was also computed using M06-2X (Figure S1 of Supporting Information), which is consistent to the findings using a higher level of theory.71 Thus, the calculations unequivocally demonstrate that M06-2X group provides a better estimate of the medium-range interactions in terms of the depth and form of the potential energy well. 3.1.5. Comparison of M06-2X and SCC-DFTB-D. As evident in previous studies, the use of SCC-DFTB-D in modeling the redox chemistry of flavin and flavoenzymes provided reasonably accurate energetics of flavin reduction.24,37,47 We evaluated its performance in predicting energetics of flavin-aromatic association using dispersion-corrections.45 The gas-phase calculation shows that geometry obtained for the optimized E



Article

ability to discern the small variations in the stacking interactions between LF and different aromatic molecules (Figure 3). Although SCC-DFTB-D produced reasonably accurate absolute binding free energy (column 1, Table 3), it was less impressive in reproducing the subtle difference in binding energies of the various LF−aromatic adducts (Table 2). 3.2. Interactions with Biologically Relevant Aromatic Rings. In enzyme active sites, the flavin ring interacts with hydrophobic side chains such as phenyl alanine, tryptophan, and tyrosine. The stacking interactions of the three biologically relevant aromatic ring systems, namely, benzene (for phenyl alanine), indole (for tryptophan), and phenol (for tyrosine) were investigated using the M06-2X functional. The Gibbs’ free energies for the stacked complexes (LF--benzene, LF--phenol, and LF--indole) in gas-phase as well as in protein interior region were computed (Table 4). 3.2.1. Characterization of the Binding Complexes. The calculated gas-phase Gibbs’ binding free energies (first row, Table 4) varied between +1.9 to −2.5 kcal/mol indicating that the binding of these biological substituents are weak. The LF-indole complex exhibited the strongest stacking interactions, while the LF--benzene had the lowest affinity with a positive (+1.9 kcal/mol) value for the free energy of binding. The average intermolecular distance of 3.25 Å was observed in LF-indole, while hovering above the pyrazine and uracil rings of the flavin (Figure 5). Comparatively, the benzene and phenol rested above the xylene and pyrazine with intermolecular distances of 3.40 and 3.35 Å, respectively (Figure 5). Using implicit solvation model, the binding affinities for these complexes were also determined in an environment that mimics the protein interior region. A dielectric constant of 4.2 was chosen, for continuum electrostatics calculation, which has been found to account for the electronic polarization and minor backbone fluctuations in protein active site.64 The decomposition of the free energy changes indicates that upon complexation there was a small reduction of the solvation free energies, presumably because of reduction of solvent accessible surface area (Table 4). This is similar to what was observed for larger polycyclic aromatics, indicating that solvation favors dissociation of these adducts. The observed weak interactions are consistent with the versatile chemistry exhibited by flavoenzymes. The reactivity of these enzymes stems out of flavin-containing cofactors that catalyze electron transfer reactions with diverse substrate molecules. The electron transfer process requires a substrate to come within the van der Waals contact of the plane of the flavin ring. If common aromatic groups of hydrophobic side chains in a protein had a stronger affinity for the flavin ring, its binding would be competitive with that of the natural substrate

complexes are comparable to the DFT results (Table 1). In addition, the SCC-DFTB-D-computed binding free energies are within 2 kcal/mol compared to that computed by the M062X functional in DFT (Table 3). These results show that Table 3. Decomposition of the Binding Free Energies (ΔG°(g)) of the Lumiflavin−Aromatic Adducts Computed by SCC-DFTB-Da stacked complexes LF--A LF--B LF--C LF--D LF--E LF--F

ΔG°(g) −6.2 −7.4 −6.6 −5.3 −9.2 −5.0

(−5.7) (−6.3) (−5.6) (−3.8) (−7.5) (−3.9)

ΔGS(LF···X) − ΔGS(LF) − ΔGSX) 3.6 4.8 4.3 3.2 7.0 3.0

a

As a comparison, the M06-2X-computed energies are given in parentheses. The thermal and zero-point energy corrections in the gasphase were computed using M06-2X density functional. The solvation energy differences between the product and the reactant states were computed using PCM and the M06-2X. 6-31+G(d,p) basis set were used throughout the density functional calculations. All energies are given in kcal/mol.

although SCC-DFTB-D can be used to model the geometry of these association complexes, a small perturbative correction with M06-2X-computed energies would increase the accuracy of the assessed energetics. This approach could be of profound significance in predicting the energetics of flavoenzymes’ catalysis. In the above calculations, the contribution of solvation in the energetics was computed by M06-2X (Table 3). These data indicate that solvation does not favor the complexation process. As observed, for each of the stacked complexes, the solvation energy of adducts (column 3, Table 3) are less than the sum of the individual solvation free energies of LF and the aromatic molecule. Thus, it is apparent that the complexation involves two opposite effects: the stacking interactions and solvation. The first decreases the energy and favors complexation, while the solvation effect tends to resist it. Therefore, the study shows that the attractive stacking interaction is the predominant contributory factor for these stacking complexes. Taken together, it appears that M06-2X produced the most accurate energetics and exhibited impressive ability to discriminate various polycyclic aromatic molecules in the condensed-phase (Table 2, Figure 3). In contrast, mPW1PW91-D3, TPSSTPSS-D3, and wB97X-D did underestimate the binding free energy by various amounts (Table 2). Nevertheless, all three functionals exhibited moderate-to-strong

Table 4. Decomposition of the Gibbs’ Free Energies of Binding (in kcal/mol) for the Stacked Lumiflavin Complexes of Benzene, Indole, and Phenol (Represented as “S” in Column 1) Mimicking the Protein Residues Phenyl Alanine, Tryptophan, and Tyrosinea free energy changes association reaction LF + S → LF--S

a

ΔGo(g) ΔGS(LF) ΔGS(A) ΔGS(LF--A) ΔGo(protein)

LF--benzene

LFindole

LF--phenol

1.9 −11.0 −1.3 −10.3 3.9

−2.6 −11.0 −3.4 −9.6 2.2

−1.5 −11.0 −3.2 −10.3 2.3

The protein interior was approximated by using implicit solvation with a dielectric constant of 4.2. F



Article

Table 5. Comparison of the Born−Oppenheimer Potential Energies of Association for Biologically Relevant Lumiflavin− Aromatic Complexes Using Various Basis Setsa Born−Oppenheimer potential energy of association

a

stacked complexes

6-31+G(d,p)

6-311+G(d,p)

6-311+G(2d,p)

6-311+G(2d,2p)

aug-ccpVTZ

uncertainty

LF--benzene LF--indole LF--phenol

−9.1 −13.7 −12.8

−9.8 −14.6 −13.9

−8.5 −12.9 −12.0

−8.3 −12.6 −11.8

−7.8 −12.0 −11.0

0.8 1.0 0.9

The last column shows the standard deviation for a specific adduct across all basis sets.

Table 6. Changes in the Computed Geometric and Energetic Properties for the Hydride-Transfer Reaction of Flavin in the Free State, in Complex with Non-Hydride Donor π-Electronic Systems, and with a Hydride Donor π-Electronic Moleculea interplanar distanceb (Å) reactions −

−

LF+ H → LFH LF-- NAPTc,d + H− → LFH−--NAPT LF--NMH → LFH−--NM+

bend angles (deg)

R

TS

P

R

TS

P

ΔGo(aq) (Eo)

− 3.31 3.26

− − 2.83

− 3.33 3.25

180 180 180

− − 174

152 155 144

−190.0 (−160 mV) −190.2 (−156 mV) −5.6

a Naphthalene, and N-methyl nicotinamide are abbreviated as NAPT, and NM respectively. For a particular reaction, R, TS, and P refer to the reactant, the transition state, and the product of that reaction. The mid-point potentials, Eo (in parentheses) are calculated using eq 3. The optimized geometry and zero-point corrected Gibbs’ free energies were computed using M06-2X functional and 6-31+G(d,p) basis set. All energies are in kcal/ mol. bThe interplanar distance for the lumiflavin-benzoic acid complex is 3.36 Å72 cThe experimental value for the interplanar distance in LF--NAPT complex is 3.41 Å69 dThe difference of the midpoint reduction potential of the two electron reduction of the LF--NAPT complex and a lone flavin ring is ±20 mV69

molecules. This would be catastrophic as the catalytic efficiency of the flavoenzyme would decrease significantly. 3.2.2. Variation of Basis Set Sizes. The gas-phase Born− Oppenheimer potential energies of association were computed using M06-2X (Table 5) for the three stacked complexes using basis sets of varying sizes. The larger basis sets used are 6311+G(d,p), 6-311+G(2d,p), 6-311+G(2d,2p), and augccpVTZ. For a specific stacked complex, the potential energy of association fluctuated by only 1 kcal/mol (Table 5) indicating that the calculation of energetics throughout this study using 6-31+G(d,p) basis set is quite reliable. 3.3. Hydride Transfer Reactions with LF. One of the important aspects of flavin chemistry is its ability to mediate redox reactions in the form of a hydride transfer. The reduction energetics of flavin was computed following the usual thermodynamic diagram Scheme 2. As noted by Walsh and Miller, during the reductive half cycle, the hydride is added to the N5 position of the isoalloxazine ring (shown in Scheme 3a) resulting in a change of electron distribution on the flavin ring, which in turn causes butterfly bending.20 Our previous work with the local functional M06-L was also able to reproduce this butterfly bending effect due to flavin’s reduction.21 As M06-2X produced better results with association energetics, we studied the electron-transfer induced effects on the geometry and energetics of the π−π stacked complexes of flavin in the presence of a nonhydride donor, as well as hydride donor πelectron rich substrate. 3.3.1. Electron Transfer-Induced Effects on Flavin Complexed with a Nonhydride Donor π-Electron System. We studied the effect of redox-mediated changes on the geometric and energetic features of the flavin-aromatic stacked complexes. In particular, we focused on a known naphthalene (NAPT) adduct of flavin because of the availability of experimental data. In each case, the reduction reaction is defined by a hydride addition to the flavin N5 atom as shown in Scheme 3a. The distances between the rings and the butterfly bend angles of flavin are presented in Table 6. The interplanar distance was calculated by taking the average of the three closest separations

between the two optimized ring systems. The computed interplanar distances for LF NAPT complexes of oxidized flavin was ∼3.3 Å and are quite consistent with the X-ray crystallographic results (Table 6).72 As evident from the interplanar distances, the separation between two planes does not change upon reduction. However, the butterfly angle values in Table 4 indicate that the reduction causes ∼30° bending of the flavin ring. The impact of stacking interactions on the energetics of these two π−π complexes were analyzed by computing the Gibbs’ free energies of the 2e− /1H+ transfer reaction. Table 6 shows that the aqueous free energies of the reduction (ΔGo(aq)) of the two stacked complexes exhibit no perceptible difference from that of the free LF. This is consistent with the cyclic voltammetric results as the difference between free flavin and the derivatized flavin-naphtalene complex69 was only 20 mV, which is less than 1 kcal/mol. 3.3.2. Electron Transfer-Induced Effects of Flavin Complexed with a Hydride Donor π-Electronic System. N-Methyl nicotinamide (NMH) is a known reductant of LF (Scheme 3b). The reactant, the transition state, and the product of the hydride transfer reaction between LF and NMH in aqueous state were optimized using M06-2X. The geometric analysis demonstrates a major change in the adduct occurring in the transition state. The interplanar distances between the flavin and nicotinamide ring systems are ∼3.3 Å for both the reactant and the product states. In contrast, the interplanar separation is only 2.8 Å for the transition state (Table 6). The 0.5 Å reduction is quite significant and speaks about the importance of the stacking interactions to shape the activation barrier. The reduced flavin is significantly bent with a butterfly bend angle of 144°. The computed barrier was found to be 17.3 kcal/mol (Figure 6a, Table 7), which is consistent to the experimental73 kcat value of 17 s−1 (an estimated value of 16.7 kcal/mol using a transmission factor of 1). The obtained barrier is consistent to the hydride transfer barrier obtained with acyl-CoA dehydrogenase using AM174 (with specific reaction parameters), where G



Article

Figure 6. Stacked lumiflavin complexes of (a) benzene, (b) indole, and (c) phenol. The left panel has the side view of the complexes, while the right panel shows a view from the top. The average interplanar distance observed for these complexes are 3.35 Å for LF--benzene, 3.25 Å for LF--indole, and 3.40 Å for LF--phenol.

Table 7. Energetics of the Hydride Transfer Reaction between Lumiflavin and N-Methyl Nicotinamide Using Two Different Theoretical Formalisms: DFT (M06-2X) with Implicit Solvation and SCC-DFTB-D/MM with Explicit Solvationa calculation method DFT(M062X)

SCC-DFTBD/MM

solvation treatment

components of free energy

ΔGo(aq) (kcal/mol)

ΔG#(aq) (kcal/mol)

implicit

ΔGuncorr.

−6.3

20.5

explicit

ΔZPE ΔGcorr. ΔGuncorr.

0.7 −5.6 −1.0

−3.2 17.3 16.0

ΔZPE ΔHLC ΔGcorr.

0.7 −4.0 −4.3

−3.2 5.1 17.9

Figure 7. Free energy profile for the hydride transfer reactions between lumiflavin and N-methyl nicotinamide modeled using (a) DFT with M06-2X functional and implicit solvation and (b) SCCDFTB-D with explicit molecular mechanically treated solvent molecules. The SCC-DFTB-D energies are not corrected for zeropoint energy and high-level-correction as given in Table 7.

interactions play a key role in shaping the transition state of flavin’s hydride transfer reaction.

a

4. CONCLUSION Density functionals are known for their poor ability to model van der Waals interactions.38 In recent years, a number of density functionals have been developed that can efficiently deal with noncovalent interactions beyond short-range (3−5 Å and >5 Å) distances. In this study we have studied the performances of a number of new-generation density functionals as well as SCC-DFTB-D to describe the stacking interactions between flavin and π-electron rich molecules. This is the first theoretical report delineating the critical role the π−π stacking has on the thermodynamics and kinetics of flavin’s redox chemistry. Our results show that dispersion correction provides a major boost in predicting adduct formation. The association of six aromatic polycyclic molecules (Figure 1) with LF was studied. The geometric and energetic results (Tables 1 and 2) obtained using the semiempirical formalism SCC-DFTB-D as well as density functionals TPSSTPSS-D3 and WB97X-D were quite consistent to the experimental results (Table 2). The M06-2Xcomputed Gibbs’ binding free energies were within 0−2 kcal/ mol and display a strong correlation to experimental results (R2 = 0.9) (Figure 3).

The corrective terms used in QM/MM calculations are described in subsection 2.3.

a barrier of 18 kcal/mol was computed for flavin’s hydride transfer reaction.75,76 The activation energy barrier for the reaction was also computed using explicitly treated water molecules and SCC-DFTB-D (Table 7). Potentials of mean force calculation was carried out for the hydride transfer reaction, whose reaction coordinate is shown in Figure 6b. The observed activation barrier without zero-point correction was 16.0 kcal/mol. The corrected activation barrier using eq 4 was found to be 17.9 (Table 7) kcal/mol, which is consistent with the DFT-computed value73 of 17.3 (Table 7) kcal/mol. Analysis of molecular orbitals indicates that the HOMO is localized on the nicotinamide ring in the reactant state but on the flavin ring in the product state consistent with the hydride transfer (Figure 7). The HOMO is shared between the two planes in the transition state with significant overlap of the molecular orbitals (Figure 7). This observation is also consistent to the significant reduction of interplanar distance (∼0.5 Å) and provides a solid evidence that the aromatic H



Article

The variation of the potential energy (Figure 4) with the interplanar separation was also studied for one of the representative complexes (LF---A). It appears that the depth of the potential energy well are similar for SCC-DFTB-D, TPSSTPSS-D3, WB97X-D, and M06-2X. The decaying tail of the potential energies, shown in Figure 4, were found to be varying as a function of 1/r6 (r = interplanar separation) indicating predominantly dispersive nature of the interaction (Figure 5). The correlations were found to be rather low for TPSSTPSS-D3 and mPW1PW91-D3. The computed energetics by SCC-DFTB-D and wB97X-D exhibited medium correlations. In contrast, M06-2X (and M06) displayed an excellent correlation (R2 = 0.99), indicating the most reliable treatment of the dispersive interactions. These results clearly shows the superiority of M06-2X over other functionals in assessing the weak interacting forces between flavin and aromatic molecules. The decomposition of the interaction energies (Table 3 and 4) suggests that the complexation involves two opposite thermodynamic effects: the stacking interactions stabilizes the adduct formation while the solvation disfavors it. The attractive stacking interaction outweighs the solvation effect and therefore is the predominant contributory factor for these stacked complexes of flavin. The stacking interactions of the three biologically relevant aromatic ring systems, namely, benzene (for phenyl alanine), indole (for tryptophan), and phenol (for tyrosine) were also investigated. The study shows positive Gibbs’ free energies of binding (Table 4) demonstrating weak stacking interactions between flavin ring and these aromatic groups, emanating from the hydrophobic side chains of proteins. This observation is supportive of the relative ease of the substrate’s accessibility of the flavin ring during electron transfer reaction. The accuracy of the energetics for these complexes was evaluated using larger basis sets. The calculated Born−Oppenheimer potential energies of association fluctuated within 1 kcal/mol across various basis sets (Table 5). This indicates that the chosen basis set 6-31+G(d,p) is suitable for studying these association complexes. The effect of stacking interactions on the electron transfer of the LF ring was studied using M06-2X. The reduction potential of the naphthalene-bound LF calculated using M06-2X shows little change (16 mV) from that of the free LF (Table 6), which is quite consistent with the 20 mV difference observed by cyclic voltammetric experiment.69 This also indicates an insignificant effect of the weak stacking interactions on flavin reduction. Similar to the free LF state, a significant butterfly bending was observed for the reduced flavins (Table 6). The computed activation barrier for the hydride transfer reaction between LF and a known hydride donor substrate, NMH was found to be within 0.5 kcal/mol of the experimental result (Table 7).73 The resultant negative Gibbs’ free energy of the reaction computed in the study demonstrated a favorable hydride transfer reaction (Table 7). The hydride transfer reaction presents an interesting geometric change during the transition state. Analysis of the structures suggests that the two parallel rings moved closer by ∼0.5 Å during the transition state with an interplanar separation of 2.8 Å. Analysis of HOMO in the transition state further revealed extensive electron delocalization indicating greater role of the π-stacking interactions during the hydride transfer reaction (Figure 8). As evident in previous studies, the use of SCC-DFTB-D in modeling the redox chemistry of flavin and flavoenzymes

Figure 8. Relative position of lumiflavin and N-methyl nicotinamide in the reactant, transition state, and product in the hydride transfer reaction. The top panel shows the HOMOs obtained in the DFT/ M06-2X calculation using implicit solvation, while the bottom panel exhibits the conformations obtained in SCC-DFTB-D/MM simulations using a stochastic boundary with an explicitly treated solvent sphere of 30 Å.

provided reasonably accurate energetics of flavin reduction.24,37,47 We evaluated its performance in predicting energetics of flavin-aromatic association using dispersioncorrections.45 The gas-phase calculation shows that geometry obtained for the optimized complexes are comparable to the DFT results (Table 1). The SCC-DFTB-D-computed binding free energies are also quite consistent (within 1 kcal/mol) to the experimental results (Table 2). Closer scrutiny of the energy components revealed that the gas-phase binding free energies in SCC-DFTB-D are within 2 kcal/mol compared to the M06-2X-computed values (Table 3). These results show that although SCC-DFTB-D can be used to model the geometry of these association complexes, a small perturbative correction would be required for assessing energetics. This approach could be of profound significance in predicting the energetics of flavoenzymes’ catalysis. Computation of the hydride transfer reaction using SCCDFTB-D/MM with explicit water molecules shows that the activation barrier is underestimated by ∼2 kcal/mol compared to the value obtained using DFT (M06-2X) with implicit solvation (Table 7). To account for this difference, perturbative corrections using M06-2X were applied (Table 7). SCC-DFTB has been previously used for flavoenzymes to model redox changes of flavin. The present study shows that the combined SCC-DFTB-D/MM simulation followed by M06-2X-based perturbative corrections would be a promising tool to model flavin’s hydride transfer reactions, which is dominated by effects due to delocalized π-electron density.

■

ASSOCIATED CONTENT

S Supporting Information *

Additional figure and discussion for the potential energy variations for the dibenzene complex. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*(S.B.) E-mail: [email protected]. Author Contributions †

Contributed equally to the study

Notes

The authors declare no competing financial interest. I



■

Article

(19) Ju, S. Y.; Papadimitrakopoulos, F. Synthesis and redox behavior of flavin mononucleotide-functionalized single-walled carbon nanotubes. J. Am. Chem. Soc. 2008, 130, 655. (20) Walsh, J. D.; Miller, A.-F. Flavin reduction potential tuning by substitution and bending. J. Mol. Struct. (Theochem) 2003, 623, 185. (21) North, M. A.; Bhattacharyya, S.; Truhlar, D. G. Improved Density Functional Description of the Electrochemistry and StructureProperty Descriptors of Substituted Flavins. J. Phys. Chem. B 2010, 114, 14907. (22) Sollner, S.; Deller, S.; Macheroux, P.; Palfey, B. A. Mechanism of flavin reduction and oxidation in the redox-sensing quinone reductase Lot6p from Saccharomyces cerevisiae. Biochemistry 2009, 48, 8636. (23) Mueller, R. M.; North, M. A.; Yang, C.; Hati, S.; Bhattacharyya, S. Interplay of Flavin’s Redox States and Protein Dynamics: An Insight from QM/MM Simulations of Dihydronicotinamide Riboside Quinone Oxidoreductase 2. J. Phys. Chem. B 2011, 115, 3632. (24) Rauschnot, J. C. J.; Yang, C.; Yang, V.; Bhattacharyya, S. Theoretical determination of the redox potentials of NRH:quinone oxidoreductase 2 using quantum mechanical/molecular mechanical simulations. J. Phys. Chem. B 2009, 113, 8149. (25) McCorrmick, D. B.; Li, H.-C.; MacKenzie, R. E. Spectral evidence for the interaction oi riboflavin with aromatic hydrocarbons. Spectrochim. Acta 1967, 23A, 2353. (26) Johnson, E. R.; Keinan, S.; Mori-SaIǹ chez, P.; ContrerasGarciIa,̀ J.; Cohen, A. J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498. (27) Antony, J.; Grimme, S. Is Spin-Component Scaled SecondOrder M?ller-Plesset Perturbation Theory an Appropriate Method for the Study of Noncovalent Interactions in Molecules? J. Phys.Chem. A 2007, 111, 4862. (28) Sinnokrot, M. O.; Sherrill, C. D. High-Accuracy Quantum Mechanical Studies of pi-pi Interactions in Benzene Dimers. J. Phys. Chem. A 2006, 110, 10656. (29) Orville, A. M.; Lountos, G. T.; Finnegan, S.; Gadda, G.; Prabhakar, R. Crystallographic, spectroscopic, and computational analysis of a flavin C4a-oxygen adduct in choline oxidase. Biochemistry 2009, 48, 720. (30) Wang, X. L.; Quan, J. M. Intermediate-assisted multifunctional catalysis in the conversion of flavin to 5,6-dimethylbenzimidazole by BluB: a density functional theory study. J. Am. Chem. Soc. 2011, 133, 4079. (31) Salzmann, S.; Martinez-Junza, V.; Zorn, B.; Braslavsky, S. E.; Mansurova, M.; Marian, C. M.; Gartner, W. Photophysical properties of structurally and electronically modified flavin derivatives determined by spectroscopy and theoretical calculations. J. Phys. Chem. A 2009, 113, 9365. (32) Khrenova, M. G.; Domratcheva, T.; Schlichting, I.; Grigorenko, B. L.; Nemukhin, A. V. Computational characterization of reaction intermediates in the photocycle of the sensory domain of the AppA blue light photoreceptor. Photochem. Photobiol. 2010, 87, 564. (33) Ai, Y. J.; Tian, G.; Liao, R. Z.; Zhang, Q.; Fang, W. H.; Luo, Y. Intrinsic property of flavin mononucleotide controls its optical spectra in three redox states. Chemphyschem. 2011, 12, 2899. (34) Rieff, B.; Mathias, G.; Bauer, S.; Tavan, P. Density functional theory combined with molecular mechanics: the infrared spectra of flavin in solution. Photochem. Photobiol. 2011, 87, 511. (35) Rieff, B.; Bauer, S.; Mathias, G.; Tavan, P. DFT/MM description of flavin IR spectra in BLUF domains. J. Phys. Chem. B 2011, 115, 11239. (36) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (37) Bhattacharyya, S.; Stankovich, M. T.; Truhlar, D. G.; Gao, J. Combined quantum mechanical and molecular mechanical simulations of one- and two-electron reduction potentials of flavin cofactor in water, medium-chain acyl-CoA dehydrogenase, and cholesterol oxidase. J. Phys. Chem. A 2007, 111, 5729. (38) Zhao, Y.; Truhlar, D. G. Density functionals with broad applicability in chemistry. Acc. Chem. Res. 2008, 41, 157.

ACKNOWLEDGMENTS This work was supported by XSEDE Grant (CHE-110018), NSF Grant (CHE-1229354), and the Office of Research and Sponsored Programs, University of WisconsinEau Claire. We gratefully acknowledge the computational support from the XSEDE-supported Trestles cluster of the University of California at San Diego and the in-house Blugold Supercomputing Cluster, University of WisconsinEau Claire. We thank Learning and Technology Services of the University of WisconsinEau Claire for their assistance in the computing facilities.

■

REFERENCES

(1) Joosten, V.; van Berkel, W. J. Flavoenzymes. Curr. Opin. Chem. Biol. 2007, 11, 195. (2) Chaiyen, P. Flavoenzymes catalyzing oxidative aromatic ringcleavage reactions. Arch. Biochem. Biophys. 2009. (3) Macheroux, P.; Kappes, B.; Ealick, S. E. Flavogenomics–a genomic and structural view of flavin-dependent proteins. FEBS J. 2011, 278, 2625. (4) Nandwana, V.; Samuel, I.; Cooke, G.; Rotello, V. M. Aromatic stacking interactions in flavin model systems. Acc. Chem. Res. 2013, 46, 1000. (5) Leung, K. K.; Shilton, B. H. Chloroquine binding reveals flavin redox switch function of quinone reductase 2. J. Biol. Chem. 2013, 288, 11242. (6) Deller, S.; Macheroux, P.; Sollner, S. Flavin-dependent quinone reductases. Cell. Mol. Life Sci. 2008, 65, 141. (7) Venkataram, U. V.; Bruice, T. C. On the Mechanism of FlavinCatalyzed Dehydrogenation a,b to an Acyl Function. The Mechanism of 1,5-Dihydroflavin Reduction of Maleimides. J. Am. Chem. Soc. 1984, 1984, 5703. (8) Wu, K.; Knox, R.; Sun, X. Z.; Joseph, P.; Jaiswal, A. K.; Zhang, D.; Deng, P. S.; Chen, S. Catalytic properties of NAD(P)H:quinone oxidoreductase-2 (NQO2), a dihydronicotinamide riboside dependent oxidoreductase. Arch. Biochem. Biophys. 1997, 347, 221. (9) Legrand, Y. M.; Gray, M.; Cooke, G.; Rotello, V. M. Model systems for flavoenzyme activity: relationships between cofactor structure, binding and redox properties. J. Am. Chem. Soc. 2003, 125, 15789. (10) Mansoorabadi, S. O.; Thibodeaux, C. J.; Liu, H. W. The diverse roles of flavin coenzymes–nature’s most versatile thespians. J. Org. Chem. 2007, 72, 6329. (11) Miura, R. Versatility and specificity in flavoenzymes: control mechanisms of flavin reactivity. Chem. Rec. 2001, 1, 183. (12) Cukier, R. I. Proton-Coupled Electron Transfer Reactions: Evaluation of Rate Constants. J. Phys. Chem. 1996, 100, 15428. (13) Cuello, A. O.; McIntosh, C. M.; Rotello, V. M. Model Systems ̂ Hydrogen Bonding in for Flavoenzyme Activity. The Role of N(3)â’H Flavin Redox Processes. J. Am. Chem. Soc. 2000, 122, 3517. (14) Goodman, A. J.; Breinlinger, E. C.; McIntosh, C. M.; Grimaldi, L. N.; Rotello, V. M. Model systems for flavoenzyme activity. Control of flavin recognition via specific electrostatic interactions. Org. Lett. 2001, 3, 1531. (15) Sharifi, R.; Samaraweera, M.; Gascon, J. A.; Papadimitrakopoulos, F. Thermodynamics of the quasi-epitaxial flavin assembly around various-chirality carbon nanotubes. J. Am. Chem. Soc. 2014, 136, 7452. (16) Breinlinger, E. C.; Rotello, V. M. Model Systems for Flavoenzyme Activity. Modulation of Flavin Redox Potentials through ̈ l€-Stacking Interactions. J. Am. Chem. Soc. 1997, 119, 1165. (17) Breinlinger, E. C.; Keenan, C. J.; Rotello, V. M. Modulation of Flavin Recognition and Redox Properties through Donor Atom-π Interactions. J. Am. Chem. Soc. 1998, 120, 8606. (18) Ju, S. Y.; Kopcha, W. P.; Papadimitrakopoulos, F. Brightly fluorescent single-walled carbon nanotubes via an oxygen-excluding surfactant organization. Science 2009, 323, 1319. J



Article

(39) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: nonempirical meta-generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 146401. (40) Perdew, J. P.; Burke, K.; Wang, Y. Generalized gradient approximation for the exchange-correlation hole of a many-electron system. Phys. Rev. B 1998, 54, 16533. (41) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theor. Comp. 2006, 2, 364. (42) 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. (43) Vydrov, O. A.; Scuseria, G. E. Assessment of a long-range corrected hybrid functional. J. Chem. Phys. 2006, 125, 234109. (44) Chai, J. D.; Head-Gordon, M. Systematic optimization of longrange corrected hybrid density functionals. J. Chem. Phys. 2008, 128, 084106. (45) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen bonding and stacking interactions of nucleic acid base pairs: a density-functional-theory based treatment. J. Chem. Phys. 2001, 114, 5149. (46) Elstner, M.; Porezag, D.; Jungnickel, G.; Elstner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge densityfunctional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260. (47) Li, G.; Zhang, X.; Cui, Q. Free energy perturbation calculations with combined QM/MM potentials − complications, simplifications and applications to redox potential calculations. J. Phys.Chem. B 2003, 107, 8643. (48) Konig, P. H.; Hoffmann, M.; Frauenheim, T.; Cui, Q. A critical evaluation of different QM/MM frontier treatments with SCC-DFTB as the QM method. J. Phys. Chem. B 2005, 109, 9082. (49) Bartholow, T. G.; Sanford, B. L.; Cao, B.; Schmit, H. L.; Johnson, J. M.; Meitzner, J.; Bhattacharyya, S.; Musier-Forsyth, K.; Hati, S. Strictly conserved lysine of prolyl-tRNA Synthetase editing domain facilitates binding and positioning of misacylated tRNA(Pro.). Biochemistry 2014, 53, 1059. (50) 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.; et al. Gaussian 09, Revision A.1; Gaussian, Inc.: Wallingford CT, 2009. (51) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S. J. CHARMm: A program for macromolecular energy, minimization and dynamics calculations. J. Comput. Chem. 1983, 4, 187. (52) Elstner, M.; Cui, Q.; Munih, P.; Kaxiras, E.; Frauenheim, T.; Karplus, M. Modeling zinc in biomolecules with the self consistent charge-density functional tight binding (SCC-DFTB) method: applications to structural and energetic analysis. J. Comput. Chem. 2003, 24, 565. (53) Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M. A QM/MM Implementation of the Self-Consistent Charge Density Functional Tight Binding (SCC-DFTB) Method. J. Phys. Chem. B 2001, 105, 569. (54) Hohenberg, P.; Walter, K. Inhomogeneous Electron Gas. Phys. Rev. B 1964, 136, 864. (55) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. A 1965, 140, 1133. (56) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213. (57) Cossi, M.; Barone, V. Solvent effect on vertical electronic transitions by the polarizable continuum model. J. Chem. Phys. 2000, 112, 2427.

(58) MacKerell, A. D. J.; Bashford, D.; Bellott, M.; Dunbrack, R. L. J.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Gou, J.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586. (59) 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. (60) Ryckaert, J. P.; Ciotti, G.; Berensden, H. J. C. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327. (61) Hockney, R. W. The potential calculation and some applications. Methods Comput. Phys. 1970, 9, 136. (62) Verlet, L. Computer ″″experiments″ on classic fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev. 1967, 159, 98. (63) Brooks, C. L. I.; Karplus, M. Deformable stochastic boundaries in molecular dynamics. J. Chem. Phys. 1983, 79, 6312. (64) Li, L.; Li, C.; Zhang, Z.; Alexov, E. On the Dielectric “Constant” of Proteins: Smooth Dielectric Function for Macromolecular Modeling and Its Implementation in DelPhi. J. Chem. Theory Comput. 2013, 9, 2126. (65) Abel, B.; Buck, U.; Sobolewski, A. L.; Domcke, W. On the nature and signatures of the solvated electron in water. Phys. Chem. Chem. Phys. 2012, 14, 22. (66) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous solvation free energies of ions and ion-water clusters based on an accurate value for the absolute aqueous solvation free energy of the proton. J. Phys. Chem. B 2006, 110, 16066. (67) Lewis, A.; Bumpus, J. A.; Truhlar, D. G.; Cramer, C. J. Molecular Modeling of Environmentally Important Processes: Reduction Potentials. J. Chem. Educ. 2004, 81, 596. (68) Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300. (69) Accorsi, G.; F, B.; Farrán, A.; Herranz, F.; Claramunt, R. M.; Marcaccio, M.; Valentic, G.; Paoluccic, F.; Pinilla, E.; Torres, M. R. Intramolecular interactions and photoinduced electron transfer in isoalloxazine-naphthalene bichromophores. J. Photochem. Photobiol. A: Chem. 2009, 203, 166. (70) Grimme, S. Density functional theory with London dispersion corrections. WIRES Comp. Mol. Sci. 2011, 1, 211. (71) Sinnokrot, M. O.; Sherrill, C. D. Substituent effects in pi-pi interactions: sandwich and T-shaped configurations. J. Am. Chem. Soc. 2004, 126, 7690. (72) Shieh, H.-S.; Ghisla, S.; Hanson, L. K.; Ludwig, M. L.; Nordman, C. E. Molecular complex of lumiflavin and 2-aminobenzoic acid: crystal structure, crystal spectra, and solution properties. Biochemistry 1981, 20, 4766. (73) Blankenhorn, G. Intermolecular Complexes between N- Methyl1,4-dihydronicotinamide and Flavines. The Influence of Steric and Electronic Factors on Complex Formation and the Rate of FlavineDependent Dihydronicotinamide Dehydrogenation. Biochemistry 1975, 14, 3172. (74) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J. Am. Chem. Soc. 1985, 107, 3902. (75) Bhattacharyya, S.; Ma, S.; Stankovich, M. T.; Truhlar, D. G.; Gao, J. Potential of mean force calculation for the proton and hydride transfer reactions catalyzed by medium-chain acyl-CoA dehydrogenase: effect of mutations on enzyme catalysis. Biochemistry 2005, 44, 16549. (76) Poulsen, T. D.; Garcia-Viloca, M.; Gao, J.; Truhlar, D. G. Free energy surface, reaction paths, and kinetic isotope effecr of short-chain acyl-coa dehydrogenase. J. Phys. Chem. B 2003, 107, 9567.

K


Effect of Stacking Interactions on the Thermodynamics and Kinetics of

Recommend Documents