A Modified Arrhenius Approach to Thermodynamically Study Regioselectivity in Cytochrome P450‐Catalyzed Substrate Conversion

Abstract The regio‐ (and stereo‐)selectivity and specific activity of cytochrome P450s are determined by the accessibility of potential sites of metabolism (SOMs) of the bound substrate relative to the heme, and the activation barrier of the regioselective oxidation reaction(s). The accessibility of potential SOMs depends on the relative binding free energy (ΔΔG bind) of the catalytically active substrate‐binding poses, and the probability of the substrate to adopt a transition‐state geometry. An established experimental method to measure activation energies of enzymatic reactions is the analysis of reaction rate constants at different temperatures and the construction of Arrhenius plots. This is a challenge for multistep P450‐catalyzed processes that involve redox partners. We introduce a modified Arrhenius approach to overcome the limitations in studying P450 selectivity, which can be applied in multiproduct enzyme catalysis. Our approach gives combined information on relative activation energies, ΔΔG bind values, and collision entropies, yielding direct insight into the basis of selectivity in substrate conversion.


Introduction
CytochromeP 450 enzymes (P450s) form as uperfamily of heme-containing proteins that playa ni mportant role in the oxidative metabolism of many lipophilic xenobiotics,a sw ella s in the biosynthesis and catabolism of endogenous compounds. [1] P450sa re consideredt ob et he catalytically most diversee nzymes in nature and because of their versatility they have many potential biotechnological applications. [2][3][4][5][6][7][8][9][10][11][12][13][14][15] So far, sequences of over 300 000 isoformsh ave been determined in all domains of life. [16] In humans,P 450sa re responsible for approximately 75 %o fp hase Im etabolism of currently marketed drugs and are involved in the activation of several prodrugs and toxicants. [17] Therefore, there is great interest in predictive tools to determine the metabolic properties of P450s. So far more than thirty different types of reactions are described for P450s. [18À20] The predominantly occurring P450 reactions include C-hydroxylation, heteroatom dealkylation,e poxidation and heteroatom oxidation. Figure 1s hows the catalytic cycle for P450-mediated hydroxylation reactions consisting of 1) substrate binding,2 )one-electron reduction of the ferric iron, 3) binding of molecular oxygen to the ferrous iron,4 )a second one-electronr eduction,5 )protonation of the Fe 2 + OO À , 6) heterolytic cleavage of the hydroperoxyl bond to yield FeO 3 + ,7 )hydrogen abstraction of CÀH-bond, 8) rebound of hydroxy group, and 9) releaseo ft he product. Which step is rate limiting appears to depend on the specific combination of P450 isoform and substrate involved in the hydroxylation reaction. [21][22][23][24][25][26][27][28][29][30] For severals ubstrates the rate-limiting nature of hydrogen abstraction has been demonstrated by the kinetic isotope effect (KIE) observed after deuterium substitution. [23,27,30] The regio-(and stereo-)selectivity and specific activity of cyto-chromeP 450s are determined by the accessibility of potential sites of metabolism (SOMs) of the bound substrater elativet o the heme, and the activation barrier of the regioselective oxidation reaction(s). The accessibility of potential SOMs depends on the relative bindingf ree energy (DDG bind )ofthe catalytically active substrate-binding poses, and the probability of the substrate to adopt at ransition-state geometry.A ne stablished experimental method to measure activation energies of enzymat-ic reactions is the analysis of reactionr ate constants at different temperatures and the constructiono fA rrhenius plots. This is ac hallenge for multistep P450-catalyzed processes that involve redox partners. We introduce am odified Arrhenius approach to overcome the limitations in studying P450 selectivity,w hich can be applied in multiproduct enzyme catalysis. Our approachg ives combinedi nformation on relative activation energies, DDG bind values, and collision entropies, yielding direct insight into the basis of selectivity in substrate conversion.  It is generally accepted that the regio-and stereoselectivity of P450s are governed by 1) the preference and probability of the substrate to bind in ar eactive orientation relative to the activatedo xygen species, and 2) the activation energies of the specific oxidation reactionsa tt he exposed sites of metabolism (SOMs). [31] To predict and/or rationalize regioselectivem etabolite formation,s everal in silicoa pproaches have been used to study orientations in substrate binding and/ort he activation energy (E a )o ft he oxidation reactions involved. Computational approaches to study the orientation(s) and dynamics of substrate binding include dockingm ethods, molecular dynamics (MD) simulations, and/or binding free energy computation. [12,13,[32][33][34] Furthermore,q uantum mechanical methods such as density functional theory (DFT), [35][36][37][38][39][40][41][42] albeit in combination with molecular mechanics (MM) techniques, have been used to calculate values for E a .T here is still limited experimental data of preferred binding modes anda ctivation energies to validate the predictivity of these computational approaches. Possible binding orientations of substrates in the active site of P450s have been studied experimentallyu sing co-crystallography [43][44][45] and spin-relaxation studies. [46][47][48] However,i nc o-crystallography studies the bound substrate sometimes appears too distant from the active center to be catalytically accessible. [44,49] In such cases, rearrangement of the substrate in the active site, which can be triggered by heme-iron reduction, may well be required to adopt ap roductive complex. [50,51] In addition, resolving the electron density can be difficult when a substrate is able to bind in multiple orientations, which is not unusualfor P450s.
Ac ommone xperimental methodt od etermine activation energieso fc hemical and enzymatic reactions is the quantification of the reaction rate constant for substrate-to-productc onversion (k cat )a td ifferent temperatures, and by subsequently constructing logarithmic plots of (ln k cat )v ersus the reciprocal absolutet emperature (1/T)a ccording to the Arrhenius equation, [52,53] According to Equation (1) the slope of the linear plot obtained equals the negative value of the activation energy E a divided by the gas constant R. The pre-exponential factor A comprises the frequency or collisione fficiencya tw hich the activated enzyme·substrate complex is formed. As such, it can be considered an entropic measure for the probability to form the transition state out of the enzyme·substrate complex.
Until now,s everale xamples of Arrhenius plots of P450 catalyzed reactions have been reported. [18,[54][55][56][57][58][59] These studies show that direct application of this or related approaches to P450s has several limitations. First, the slopes of the Arrhenius plots do not necessarily represent the activation energy of the oxidation reaction because steps prior to the oxidation reaction, such as reduction by the NADPH cytochrome P450 reductase (P450 reductase) and/or cytochrome b 5 reductase (steps 2a nd 4, Figure 1), may also be rate limiting for the overall reaction. [21] Second, nonlinear Arrhenius plots were found with reactions catalyzed by microsomalP 450s with ad iscontinuitya ta pproxi-mately 20 8C, which wasa ttributed to at ransition of membrane fluiditya ffecting the interaction between P450 and P450 reductase. [57,58] Furthermore, usually only as mall range of temperatures is used because above the optimal temperature, the reactionr ates decrease againd ue to enzymed enaturation. Finally,P 450 catalysis often leads to differentp roducts at different ratios. Experimentally, k cat values are usually determined by dividing values for the maximal velocityi ns ubstrate conversion as obtained from enzyme kinetic studies (V max )b yt he total enzyme concentration [E] total (k cat = V max /[E] total ), assuming that at maximal enzyme activity all enzymes are occupied by the substrate in ar eactive binding pose. However, in the case of parallel reactions the individual V max values cannot be divided by [E] total but should be divided by the concentration of the enzyme-substrate complexeso ft he corresponding reactive binding poses, designated [ES 1 ]a nd [ES 2 ]i nS cheme 1f or possible formation of two products (P 1 and P 2 ). No direct experimental methods are available to accurately determine ratios of different bound conformations of ag iven enzyme·substrate complex.
In the present study we proposea nd evaluate am odified Arrhenius approach to overcome several of the above-mentioned intrinsic limitations of thermodynamic studies on P450catalyzed reactions. In this modified approach the temperature dependence of the ratio of V max values of parallel reactions is analyzed rather than studying the temperature dependence of kinetic parameters for individual pathways. Dividing the Arrhenius-equations of the competing reactions in Scheme 1( i.e.,o f reactions 1a nd 2w ith maximal velocities V max,1 and V max,2 and rate constants k cat,1 and k cat,2 ,r espectively)a nd using V max = k cat [ES] max gives [Eq. (2)]: The subscript max for concentrationso ft he enzyme-substrate complexes indicatesE S, ES 1 or ES 2 concentrations at maximal enzymea ctivity.U nder as teady-state approximation the ratio between [ES 1 ]a nd [ES 2 ]c an be related to with DDG bind = DG bind,1 ÀDG bind,2 ,w hich is the (possible)d ifference in binding free energies DG bind between substrate binding poses associated with formation of products 1a nd 2, respectively.T he steady-state approximation of Equation (3) is valid when substrate concentration [S] @ [E] total ,a nd when the constantso fs ubstrate binding and unbinding( k b and k u in Scheme1)a re substantially higher than those for formation of the enzyme-product (k p,1 and k p,2 in Scheme 1) and/or alternatively,w hen rapid interconversion between ES 1 and ES 2 is possible (with k i and k -i in Scheme 1b eing highert han the k p values). Hence, when correlating the natural logarithm of the ratio V max,1 /V max,2 with the inverse absolutet emperature (in a modified Arrhenius plot), as traight line is expected of which the slope of the plot will represent the sum (D)o ft he differences in activation energy (E a )a nd DG bind of the parallel reactions (divided by the gas constant R), Equations (4) and (5) ln where D = DDG bind + DE a ,a nd DE a = E a,1 ÀE a,2. Because we assumet hat the rates and temperature dependence of steps 2-6 of the catalytic cycle in Figure 1w ill be similar for the parallelp athways, these factors will cancelo ut when evaluating the ratios of product formation and therefore will not contributet od ifferences in D in Equation ( 5). In addition, changes in membrane fluidity( in case of microsomal P450s), suboptimal interaction between P450 and P450 reductases, and protein denaturation at increased temperatures are expected to affect V max values of both pathways to as imilar extent,w hich will allow studies over al arger temperature range.
In previous studies, ratios of product formation obtained at one incubation temperature have been used to estimate the overall difference in activation free energies( DDG overall )f or competing reactions by P450s, Figure 2. [25,60] Also in this case, possible ratel imiting factorsp rior to the oxidation reaction (steps 2-6, Figure 1) are expected to cancelo ut. According to the Curtin-Hammett principle, [61] the product ratio of two competing reactions is governed by the difference DDG overall between the free energies of the correspondingt ransition states Figure 2) when the barrier to interconversion between reactive binding poses ES 1 and ES 2 (either direct or via enzyme/substrate unbinding) is much smaller than the barriert op roduct formation. The differences in the free energy of formingt hese transition states dependo nt he one hand on differences in free energies of binding of the pro-ductive binding poses ES 1 and ES 2 (DDG bind )a nd secondly by the differencesb etween free energies of the binding poses and the transitions tates involved in product formation,w ith DDG°= DG 1°-DG 2°, Figure 2A [Eq (6)]. ln ½product 1 Previously,H iggins et al. assumed that when different humanP 450s show both similark inetic isotope effects and product ratios, the product ratios observed are determined by the differences in activation free energyo ft he transition states DDG°[Eq (7)]: [60] ln ½product 1 This may wellb ev alid for small substrates that can rapidly adopt multiple binding poses (high k i and k -i in Scheme1)a nd that have only as mall differencei nf ree energy of binding Figure 2. Gibbs free energy profile for acatalyzed reaction with two possible products A) with ab inding free energy difference and B) without abinding free energy difference between the binding poses leadingtot he different products P 1 and P 2 ,s tarting from enzyme and substrate (E + S). DG°, x comprises the sum of an activation energy( E a,x )a nd ac ollision entropy (TDS°, x ) term.
(DDG bind ), cf. Figure 2B.I ns uch cases, DDG°can be directly estimatedf rom ratios in product formation using Equation (7). In the case of substrates with high molecular weight and/or P450s with restrictive active sites, next to differencesi nf ree energieso factivation, steric factorsmay also play an important role in the regioselectivity of P450 reactions, and DDG bind in Equation (6) cannot be neglected apriori. This is exemplified for example by recent bindingf ree energy calculations [32] for one of the pairs of product formation considered in the current work, and Equation (6) is in such cases to be used instead of Equation (7). The differencei nu sing our modified Arrhenius approachc ompared to direct use of the Curtin-Hammett formalism is that the entropic contribution to differencesi nt he barrierf or forming the transition state from the enzyme-substrate complex (i.e., the ratio A1/A2 in Equation (5)) can be separated from other contributions (DDG bind and DE a ). Therefore, our approachcan be of direct helpinvalidating (thecombined use of) free energy,q uantum chemical, and MD studies on preferred modes of substrate binding, activation energies, and/or probabilities to adopt catalytically activeb inding orientations, respectively.
In this study,t he human P450 isoform1 A2 (CYP1A2) and a drug metabolizing mutant of bacterialP 450 BM3 (CYP102A1), that is, BM3 M11, are used to evaluate the applicabilitya nd to illustrate the value of our approach to analyze thermodynamic determinants of selectivity in P450-catalyzed product formation. For that purpose, we determinedt he temperature dependence of product ratios for pairs of different substrate conversions as catalyzed by the same isoform. Mefenamic acid (MF) and testosterone (TE) were selected as substrates. MF is oxidizedb yC YP1A2 and P450 BM3 M11t ot wo or three metabolites, respectively,w hile TE conversion catalyzed by BM3 M11l eads to three different products as well, Figure 3. [13,62] Recombinant CYP1A2 was selected as model for am embranebound P450, which dependso nc o-expressed NADPH cyto-chromeP 450 oxidoreductase as redox partner.P 450 BM3 M11 was used as model for as oluble P450. Wildtype P450 BM3 is a natural fusion protein betweenaP 450 domain and P450 oxidoreductase domain and is often used for mechanistic studies of P450. [4] Because it has the highest turnover recorded for any P450, it also has promising biotechnological perspective for biosynthesis of fine chemicals. [11] Mutant P450 BM3 M11w as developed by combination of site-directed and random mutagenesisa nd catalyzes oxidation reactions of aw idev ariety of pharmaceuticals and other chemicals. [6] To determine D and relative collisione fficiencies ln (A 1 /A 2 )i nE quation (5) for the multiple substrate conversions catalyzed by CYP1A2 or P450 BM3 M11, enzyme kinetic parameters were determined for each reactiona td ifferent incubation temperatures. In support of our steady-statea pproximation in Equation (3), we also measured kinetic isotope effects for the pair of product formation (i.e.,c onversion of MF to either 3'-OH-MF or 4'-OH-MF by BM3 M11, Figure 3) for which the corresponding ES binding poses were previously reported to be similar, [32] and hence may well rapidly interconvert. In addition, molecular dynamics (MD) computer simulations of selected isoform-substrate combinations were carried out to quantify the probability of the substrates to adopt different catalytically active binding poses, using geometric criteria for transition state formationb ased on combined QM/MM studies by Mulholland and co-workers. [63] The resultsw ere compared with the relative collision efficiencies as determined from the intercepts of our Arrheniusp lots, and with differences between our estimated D values and corresponding Curtin-Hammett estimates for relative activation barriers, as measures for possible differencesi nt he entropyo f transition state formation.T of urtheri nterpret and cross-validate our modifiedA rrhenius and in silico analyses we also computed differences in activation energies E a (using the SMARTCyp web server) [64] and/oro btainedt hem from literature, and where possible we combined these estimates with DDG bind values reported in literature for ad irect comparison with our values for D.

Results and Discussion
Temperature-dependent mefenamic acid hydroxylation catalyzedbyB M3 M11 As described previously, [33] mefenamic acid was metabolized by P450 BM3 M11t ot he three regioisomeric hydroxy metabolites shown in Figure 3. At all incubation temperatures 4'-hydroxymefenamic acid (4'-OH-MF) was the major product, followed by 3'-hydroxymethylmefenamic acid (3'-OH-MF), and 5hydroxymefenamic acid( 5-OH-MF) as relatively minor product, Figure 4A.A ss ummarized in Ta ble 1, the catalytic efficiency (V max /K M )f or all three pathways increased from 4t o2 5 8C. At higher temperatures the catalytic efficiency decreased again. Also for the V max values the lowest values were obtained at the lowest and highest incubation temperatures.A saresult, the plots of ln V max versus 1000/T were strongly nonlinear when analyzing the kinetics of each metabolite individually (data not shown). Thisw as expected based on previousA rrhenius studies on P450 catalyzed reactions. [56][57][58][59] Before applyingo ur modifiedA rrhenius analysis to the (three)p airs of BM3 M11c atalyzed product formations (Table 2), we measured relative kinetici sotope effects in an attemptt oe xplicitly verify the steady-statea pproximation taken in Equation (3) for the ratio of MF conversion to 3'-OH-MF and 4'-OH-MF.R epresentative extracted ion chromatograms of the mixed deuterated and non-deuterated hydroxy metabolites in Figure4Bs howt hat the relative peak area of the 3'-OH-MF  Table 3s ummarizes that deuteration of mefenamic acid resultsi na na pproximately fourfold decrease in the rate of 3'-methyl hydroxylation but to a5 0% increase of the 4'-hydroxylation pathway,a nd no change in 5-OH mefenamic acid formation. These resultsi ndicate that the hydrogen abstractiono ft he 3'-methyl group is rate limiting to amore significant extent than aromatic hydroxylation reactions. Thisi si nl ine with previous studies in which substantially larger KIEs were observed for aliphatic than for aromatic hydroxylation by P450s. [30] The increasein4'-OH-MF formationa fter deuteration may be well explained by metabolic switching resulting from the strongly decreased 3'-hydroxylation. Thus, our kinetic isotope effect measurements for these reactionss uggestt hat rapid interconversion between catalytically active poses for 3'-OH-MF and 4'-OH-MF formation is possible and accordingly,that k i and k Ài values for binding-poseinterchange and/or k b and k u values are probably highert han k p for the corresponding product formations. Therefore, even in  Table 1. Te mperature dependence of enzymek inetic parametersf or the regioselective hydroxylation of mefenamica cid by P450 BM3 M11a nd recombinant human CYP1A2 (i.e.,f or the formation of the products 3'-hydroxymethylmefenamic acid, 4'-hydroxymefenamica cid and 5-hydroxymefenamic acid).

3'-Hydroxymethylmefenamic acid
4'-Hydroxymefenamic acid 5-Hydroxymefenamic acid In the next step we used Equation (5) to study the temperature dependence of the ratios in MF conversion by BM3 M11. When plotting natural logarithms of the ratios of V max values of mefenamic acid metabolites against 1000/T,l inear curves were obtained, Figure 5A.T his observed correlation supports our assumption that the temperature dependence of steps 2-6 of the catalytic cycle in Figure 1w ill be similarf or the different hydroxylation paths andc ancelw hen relating V max ratios to the inverse temperature. The slopes (D = DDG bind + DE a )a nd intercepts of the thus obtainedm odified Arrhenius plots are summarizedi nT able 2a nd discussed below in our thermodynamic analysiso ft he observed regio-specificity in MF conversion.I n addition, DDG overall values as derived (using the Curtin-Hammett principle and Equation (6)) from pairs of V max values at 300 Ka re also reported in Table 2, as wella sd ifferences in activation energies DE a calculated by using SMARTCyp.N ote that we used two versions of SMARTCyp (versions2and 3), which gave identicalr esultsf or E a .T he small differences in activation barrierf or all three hydroxylationr eactions are in line with the similar values for E a calculated at the B3LYP level of Density Functional Theory (DFT) by Leth et al. [65] Leth modeled compound Ia sap orphyrin moietyw ithout side chains and with axial coordinating O 2 À and CH 3 S À ligands, and showed for example ad ifferencei nt he range of À6t o3kJ mol À1 in activation barrier when comparing 3'-methyl-OH-MF and 4'-OH-MF formation.T his is to be compared with the corresponding SMARTCyp value of À2kJmol À1 (Table 2).
For the ratio between 3'-methyl-OH-MF and 4'-OH-MF formationt he modified Arrhenius plot ( Figure 5A)s hows as lope of À0.97 AE 0.05 Kw hich corresponds to av alue for D [in Eq. (5)] of 8.1 AE 0.5 kJ mol À1 ,w ith al ower sum of DG bind and E a for hydroxylation at the 4' aromatic SOM, Ta ble 2. As stated    previously, DE a is close to 0( and even slightly negative); therefore the positive D value should be interpreted to result from amore negative (favorable) binding free energy for the catalytic binding orientation for 4' hydroxylation of mefenamic acid as compared to 3'-methyl hydroxylation (DDG bind = DÀDE a = 8.1À1.8 = 9.9 kJ mol À1 ), in agreement with and confirming the corresponding DDG bind values previously computed by us that range between 9.3 and 11.6 kJ mol À1 . [32] Despite the lower binding affinity of MF to BM3 M11i ni ts pose that is catalytically active for 4'-OH-MF formation,t he differencew ith 3'-methyl hydroxylation in the overall activation free energy DDG overall is close to zero from our Curtin-Hammett analysis( 0.7 kJ mol À1 ,T able 2). Thus, the difference in DG bind is for al arge part counterbalanced by al ower entropyp enalty to form the transition state for 3'-methyl-OH-MF formationo ut of the corresponding ES complex, as reflected by the higherc ollision efficiency for 3'-methylhydroxylation (with av alue of 3.0 AE 0.2 for the intercepto ft he modified Arrhenius plot, Table 2). We could cross-validate these findings with MD simulations of BM3 M11i nc omplex with MF,i nw hich we compared the frequencies of occurrenceo fM Fb inding poses that can potentially adopt transition state geometriesf or either 3'-Me or 4' hydroxylation. Indeed, our simulations showedh igher frequencies for substrate orientations that are in line with the transition state geometries for 3'-Me-OH-MF than for 4'-OH-MF product formation, Table 4. The fact that 4' hydroxylation of MF by BM3 M11i soverall favorable over 3'-methyl hydroxylation, despite the highere ntropic cost for transition state formation and the slightly higher activation energy,s hould thus be understoodi nt erms of the lower binding affinity for substrate binding in ap ose that leads to 4'-OH-MF formation. Our previous detailedf ree energy perturbation computations tudy [32] on selectivity in MF hydroxylation came to this conclusion as well and could thus be verifiedw ith our modified Arrhenius approach.
For the ratio between 3'-methyl-OH-MF and 5-OH-MF formation the modified Arrheniusp lot (Figure5A) shows as lope of À0.97 AE 0.07 Kw hich corresponds to av alue for D (in Equation (5)) of 8.1 AE 0.6 kJ mol À1 ,w ith al ower sum of DG bind and E a for hydroxylationa tt he 5a romatic SOM, Ta ble 2. When considering the similara ctivation energy for 5-OH-MF as predicted by SMARTCyp (with DE a = À1.8 kJ mol À1 ,T able 2), the bindingf ree energy should be higher for the 3'-methyl bindingp ose, suggesting that binding in an orientation that can lead to 5-OH-MF formation is more favorable. As an alternative explanation, the difference in D might be due to al ower E a value for 5'-hy-droxylation,a si ndicated by additional DFT calculations of Leth et al. in which dispersion corrections were explicitly included (i.e. using the B3LYP-D3 level of theory), resulting in lower E a value for 5-OH-MF formation by 14 kJ mol À1 . [65] In any case, from Table 2t he preference of 3'-methyl over 5-OH mefenamic acid formationb yP 450 BM3 M11c an be understood in terms of the higherp robability (lowere ntropic cost) of transition state formation for 3'-methyl hydroxylation,a si ndicated by the higher intercepto fo ur modified Arrhenius plots and the observed difference between D and the Curtin-Hammett estimate for DDG overall .T his observed difference is equal to TDDS°u nder the assumptiont hat trends in enthalpy and energy are equal,c f. Equation (8) and Figure2.
This is in accord with our in silico data from refs. [32] and [33] showing strong hydrogen bondingi nteraction between mefenamic acid's carboxylate group and BM3 M11'sS er72 residue, which directss ubstrate-bindingo rientations for 3'-methyl hydroxylation to adopt ac atalytically active pose. Such an anchoringh ydrogen bond is not present when bound in ap ose enablingh ydroxylation at the 5p osition. Also in our MD simulationst his particular hydrogen bond waso bserved in the simulations with mefenamic acid in the 3'-methyl hydroxylation pose ( Figure S1), which was not observed in the 5-hydroxylation pose (FigureS2). Furthermore, we indeed observe as ubstantially higher frequencyi nM Ds imulations of substrate orientations corresponding to transition statef ormation for the 3'-methyl hydroxylationc ompared to 5h ydroxylation (Table 4).
For the ratio between 4'-OH-MFa nd 5-OH-MF formation by BM3 M11, the modified Arrhenius plot (Figure5A) shows a slope of 0.00 AE 0.07 Kw hich corresponds to av alue for D [in Eq. (5)] of 0.0 AE 0.6 kJ mol À1 ,T able 2. The predicted identical E a valuesf or both pathways ( Table 2) thus suggestasimilar binding free energy for the corresponding catalytically active poses (DDG bind = DÀDE a )w hereas the lower B3LYP-D3 value of Leth et al. (by 11 kJ mol À1 ) [65] hints at preferred binding in the pose enabling4 ' hydroxylation.T he preference of 4' over 5h ydroxylation can again be understoodi nt erms of the higherp robability (lower entropic cost) of transition state formation for 4' hydroxylation compared to 5-OH product formation ( Table 2), probablya lso due to hydrogen bonding with Serine 72 in the catalytic-active binding pose in the formercase ( Figure S3).

Temperature-dependent mefenamic acid hydroxylation catalyzed by CYP1A2
Oxidation of mefenamic acidb yr ecombinant CYP1A2 resulted in formation of 4'-OH-MF and 5-OH-MF.Atlower substrate con-centrations5-OH-MF wasthe major metabolite, as indicated by the higher V max /K M values, Table 1. V max values of the 4'-hydroxylation pathway were slightly highert han for 5-hydroxylation. Nonlinear Arrhenius plotsa re obtained when plotting ln V max versus1 000/T for the metabolites individually (data not shown). As for the pairs of BM3 M11m ediated product forma-  Figure 5B. Ah igherc ollisionf requency for formationo ft he transition state for 4' hydroxylation (compared to 5-hydroxylation,c f. the positive intercepta nd positive difference between D and DDG overall in Ta ble 2) suggestsa ne ntropically more favorable transition state formationf or 4'-OH-MF formation. This is in line with the higherfrequency observed in MD of substrate orientations corresponding to transition state formation for 4' hydroxylation (Table 5). Duringt he simulations, we observe ah ydrogen bond between mefenamic acid and Thr469w hen binding in the 5-hydroxylation position( Figure S4), whereas we do not observe any stabilizing or positioning hydrogen bonds in the 4'-hydroxylation bindingp ose ( Figure S5). The slightpreference of 5-OH mefenamic acid formation by CYP1A2 (asr eflected by the slightly negative DDG overall estimate in Ta ble 2) can in this case be associated to al ower binding free energy for and preferred binding in the corresponding catalytically active pose. Furthermore, the E a value for 5'-hydroxylation may be lower (see above).

Temperature-dependent testosterone hydroxylation catalyzed by P450 BM3 M11
As illustrated in Figure3,t estosteronec an be hydroxylated by P450 BM3 M11a tt hree positions, leading in order of their relative amounts to 15b-OH-T,1 6 b-OH-T or 2b-OH-T formation,r espectively. [46,62] The enzyme kinetic parameters of the reactions performed at temperatures rangingf rom 6t o3 68Ca re shown in Supporting Information Ta ble S1.A sw as observed for mefenamic acid, no linear Arrheniusp lots are obtained when plotting ln V max versus 1000/T for the metabolites individually (data not shown). Encouragingly,t he modified Arrhenius plots of the ratios of V max values versus 1000/T showeda gain linear behavior for all three combinationsofp athways, Figure 6.
Using the slopes of these curves, D values for the three different pairs of product formation were obtained, Table 2. The minor pathway leading to 2b-OH-T was found to have the highest sum of DG bind and E a, with D = 23.5 kJ mol À1 for the ratio with 16b-OH-T,a nd 7.2 kJ mol À1 for the ratio with 15b-OH-T.T hus, D is significantly higher for2 b-OH-T when compared with the other products,w hereas the activation energy for formation of the corresponding transition state is lowest. This is apparent from the value for E a ,w hich was predicted by SMARTCyp to be 9.5 kJ mol À1 lower than for 15b and 16b hydroxylation (Table 2), probablyd ue to the adjacent C=O moiety beside of the C 2 carbon. From previousc omparative DFT calculations on the energy of CÀHb ond breaking we even found differences of more than 20 kJ mol À1 in favor of CÀH bond activation at position 2c ompared with positions 15 and 16. [66] The higher D and lower E a values for 2b-OH-T imply a significantly higher bindingf ree energy (i.e.,l ower binding affinity) for TE binding in ap ose compatible with 2b hydroxylation than for the other pathways, which may explain why it was difficult to find suitable starting poses from docking to start MD from of BM3 M11w ith TE bound in this binding pose.
The D value for the ratio of formation of the major metabolite 15b-OH-T and the less abundant 16b-OH-T is 16.7 kJ mol À1 , Ta ble 2. E a is probably higherf or 15b hydroxylation than for 16b hydroxylation due to the substrate's hydroxy group at C17, which was also indicated by our previousf inding thatt he energy cost of breaking the C15-H and C16-H aliphatic bonds is 5kJmol À1 higher for the former. [66] In this particular case, the difference in E a might be the only contribution to D as our MD simulations suggest that interconversion between binding orientations suited for 15b-a nd 16b-hydroxylation can occur on the ns time scale. This is illustrated in Ta ble 6w hich shows that the geometric criteria for adopting the transition states fort he 15b and 16b hydroxylation pathways can be both fulfilled during as ingle simulation. Therefore, information on the frequencyo fM Dc onfigurations consistent with transition state formation wasf or both positions taken from the same set of (three)s imulations. The fact that 15b-OH-T formation is prevalent over 16b hydroxylation can be explained in terms of a higherc ollisionf requency and entropy of transition state formationo ut of the enzyme-substrate complex, as can be ob- Table 5. Percentages of substrate binding orientations during two independent MD simulations of mefenamica cid bound to CYP1A2 (A and B) that are suitable for transition state formation for hydroxylation of mefenamic acid at its 4'-o r5 -position.   served from the relatively largei ntercept in the modified Arrheniusp lot for the ratio of 15b/16b hydroxylationa nd from the according differenceb etween D and DDG overall ,T able 2. These resultsa re in line with the higher frequency observed in MD of substrate orientationsc orresponding to transition state formation for 15b hydroxylation (Table 6).

Conclusions
We have presented am ethodt hat makes it possible to obtain experimentale stimates for thermodynamic determinantso f regio-(and/or stereo-)selectivity in P450 catalyzeds ubstrate conversion,b ym eans of studying temperature dependent ratios of pairs of metabolite formation as catalyzed by as ingle P450 isoform.W ei llustrated the use of this modified Arrhenius approachb ys tudying the determinants of the regioselectivity in mefenamic acid and testosteroneh ydroxylationb yP 450 BM3 M11a nd CYP1A2.F or the selected P450-substrate combinations, our approach gave insight into the basis of selectivity by giving combined information on relative activation energies, DDG bind values, andc ollisione ntropy differences. We cross validated the observed collision entropyd ifferences with molecular dynamics simulations, and we were able to verify previousc omputational free energy calculations.T he obtained agreement suggestst hat the presented method can also be appliedt oo ther combinationso fP 450 isoformsa nd substrates that involve formation of two or more products.T he methods and experiments described here are useful tools for future research on regio-and stereoselectivity of P450 catalysis to ultimately improve biocatalysts, and the data that can be obtained with our method allow to validate resultsf rom computationalmodelstou nderstand and predict selectivity.

Experimental Section
Materials: Mefenamic acid, testosterone, NADPH, glucose 6-phosphate and glucose-6-phosphate dehydrogenase were purchased from Sigma-Aldrich (Steinheim, Germany). Supersomes containing recombinant human P450s were obtained from BD Biosciences (Breda, Netherlands). The plasmid containing P450 BM3 M11w as constructed as described earlier. [6] All other chemicals were of analytical grade and obtained from standard suppliers.
Expression of cytochromes P450 BM3 M11a nd P450 1A2: Histagged cytochrome P450 BM3 M11w as expressed by transforming competent Escherichia coli BL21 cells with the corresponding pET28a + vector and purified using nickel affinity chromatography as described previously. [62] The method of Omura and Sato was used to determine the cytochrome P450 concentration. [67] Ab icistronic plasmid containing the cDNA of human CYP1A2 cDNA and human NADPH cytochrome P450 reductase was transformed into E. coli strain DH5a.A300 mL terrific broth (TB) supplemented with 1mm d-aminolevulinic acid, 0.5 mm thiamine, 400 mLL À1 trace elements, 100 mgmL À1 ampicillin, 1mm isopropylb-d-thiogalactopyranoside (IPTG), 0.5 mm FeCl 3 was inoculated with a7 .5 mL pre-culture grown from as ingle colony.T he cells were allowed to grow for 40 ha t2 8 8Ca nd 125 rpm. E. coli cells were collected by centrifugation (4000 g,4 8C, 15 min) and resuspended in 20 mL 0.1 m potassium phosphate buffer at pH 7.4 con-taining 20 %g lycerol, v/v,0 .25 mm ethylenediaminetetraacetic acid (EDTA) and 0.1 mm dithiothreitol (DTT). The cells were treated with 0.5 mg mL À1 lysozyme for one hour at 4 8Ca nd subsequently disrupted by three cycles of disruption by Emulsiflex C3 emulsifier. The membranes containing the CYP1A2 were isolated by ultracentrifugation for 75 min at 40 000 rpm (169 936 g)a nd 4 8C. The pellet was resuspended in the potassium phosphate-glycerol buffer and subsequently homogenized by Potter-Elvehjem. The concentration of CYP1A2 was determined using the method of Omura and Sato [67] and the enzyme was stored at À80 8Cu ntil use.
Assessment of enzyme kinetic parameters at different temperatures: All incubations were performed at enzyme concentration and incubation times for which the product formation was linear in time and proportioned to enzyme concentration (data not shown). The incubation mixtures contained 50 nm P450 BM3 M11 or CYP1A2 in 100 mm potassium phosphate buffer (pH 7.4) supplemented with 5mm MgCl 2 and 2mm EDTA. Seven substrate concentrations ranging from 10 to 750 mm were used in af inal incubation volume of 100 mL. Reaction mixtures were preincubated in a shaking water bath set at the incubation temperature for 10 min. Reactions were initiated by addition of 10 %( v/v)o faprewarmed solution containing NADPH regenerating system;f inal concentrations were 0.5 mm NADPH, 10 mm glucose 6-phosphate, and 0.4 units mL À1 glucose-6-phosphate dehydrogenase. The reaction was allowed to proceed for 4min at different temperatures and then stopped by the addition of 100 mLi ce-cold methanol. The denatured enzyme fractions were precipitated by centrifugation for 20 min at 14 000 rpm (20 817 g). The supernatants were isolated and analyzed by HPLC or LC-MS as described below.
For mefenamic acid and M11e nzyme kinetic parameters for the formation of the three metabolites were determined at eight temperatures ranging from 4t o4 5 8C, to investigate the linearity in more detail. The other substrates were incubated at four temperatures ranging from 4t o3 58C. The enzyme kinetic parameters V max and K M were determined by nonlinear regression according to the Michaelis-Menten equation using GraphPad Prism 7.0 software (GraphPad, San Diego, CA, USA).
D and ln (A 1 /A 2 )i nE quation (5) for competing enzyme reactions were determined by plotting the logarithms of the ratio of the V max values measured at different temperatures against the inverse of the absolute temperature T. According to Equation (5), the slope of this curve, when linear,c orresponds to the sum D of DDG bind and DE a divided by the negative gas constant R (8.3145 Jmol À1 K À1 ). The slopes and intercepts of the modified Arrhenius plots were determined by linear regression using the GraphPad Prism 7.0 software.
Determination of competitive intermolecular kinetic isotope effects of hydroxylation of mefenamic acid: Amethod for full deuteration of mefenamic acid was described previously [68] and used to synthesize deuterated mefenamic acid and to study kinetic isotope effects for its conversion by BM3 M11 ( Figure 3).  13 ]mefenamic acid, so that eventual effects of these side products on the reactions were applicable to both labelled and unlabeled mefenamic acid to the same extent. [30] All incubations were performed with 100 nm P450 in 100 mm potassium phosphate buffer (pH 7.4) supplemented with 5mm MgCl 2 and 2mm EDTA. To tal substrate concentrations of the mixture of labelled and unlabeled mefenamic acid were 75 and 750 mm.R eaction mixtures (total volume 100 mL) were pre-warmed at the incubation temperature for 10 min, before initiating the reaction by addition of aN ADPH regenerating system (final concentrations of 0.5 mm NADPH, 10 mm glucose 6-phosphate, and 0.4 unit mL À1 glucose-6-phosphate dehydrogenase). The reaction was allowed to proceed for 4min at different temperatures and then stopped by the addition of 100 mLi ce-cold methanol. The protein was removed by centrifugation for 20 min at 14 000 rpm (20 817 g). The supernatants were analyzed on an Agilent 1200 series rapid resolution LC equipped with aT OF Agilent 6230 mass spectrometer (Agilent technologies, Waldbronn, Germany). Data processing was performed with the Mass Hunter Qualitative Analysis software package (version B.06.00). Note that the extracted ion chromatograms of the deuterated hydroxy metabolites and substrate showed slightly shorter retention times than their non-deuterated counterparts, indicating that the lipophilicity was slightly reduced upon deuteration. Assuming that the deuteration of the aromatic rings and methylene group does not affect the ionization efficiency of the metabolites, the kinetic isotope effect of full deuteration were for all three metabolites directly calculated from the peak areas in the ion chromatograms. Because kinetic isotope effects were studied at as ingle concentration, no Arrhenius plots were constructed because activities did not represent V max values.
Analytical methods: The analyses of metabolites were performed by reversed-phase liquid chromatography using aS himadzu HPLC equipped with two LC-20AD pumps, aS IL20AC autosampler and a SPD20A UV detector.L ab Solution software of Shimadzu was used to control the HPLC-system, data acquisition and data analysis. For metabolite identification and quantification of isotope ratios in the competitive isotope experiment, an Agilent 1200 series rapid resolution LC was used which was connected to an Agilent 6230 timeof-flight mass spectrometer equipped with an electrospray ionization (ESI) source operating in positive ion mode. Acapillary voltage of 3500 Vw as used, and nitrogen was used both as drying gas (10 Lmin À1 )a nd nebulizing gas (pressure 50 psig) at ac onstant gas temperature of 350 8C; 1000 MS spectra per second were acquired and analysis was performed using Agilent MassHunter Qualitative analysis software (version 2.0).
For the analysis of metabolites of mefenamic acid, the first 5min was isocratic at 40 %e luent B. From 5u ntil 30 min, the concentration of eluent Bw as increased linearly to 100 %, followed by linear decrease back to 40 %b etween 30 and 30.5 min. Isocratic re-equilibration at 40 %e luent Bw as maintained until 45 min. The flow rate was 0.5 mL min À1 .U V/Vis detection was performed at 254 nm.
For the analysis of metabolites of testosterone, the first 1min was isocratic at 50 %e luent B. From 1t o2 0min the percentage of eluent Bw as increased linearly to 99 %; and from 20 to 20.5 min linearly decreased to 50 %Band maintained at 50 %f or re-equilibration until 30 min. The flow rate was 0.5 mL min À1 .U V/Vis detection was performed at 254 nm.
Molecular dynamics (MD) simulations of substrate binding to P450 BM3 M11a nd CYP1A2:M Ds imulations were carried out to quantify the occurrences of different catalytically active binding poses over time, as an (entropic) measure of the frequency or efficiency of collisions allowing for transition state formation. Simulations were carried out both for mefenamic acid and testosterone, either bound to BM3 M11orCYP1A2. To define occurrence of binding poses that are suitable for hydroxylation of either an aliphatic or aromatic CÀHm oiety,w eu sed geometric criteria for transition state formation as reported by Mulholland and co-workers. [63] These criteria are similar as we used before and were as before extended with ar ule to exclude assignment of conformations to be catalytically active for aromatic hydroxylation, in case the angle between the CÀHs ite-of-metabolism bond and the vector connecting the corresponding hydrogen with the ferryl oxygen was between 1408 and 2208 (Table 7), in order to account for the possible detrimental effect of hydrogen interposition on CÀHa ctivation by the ferryl oxygen. [32] Thus, ag iven enzyme-substrate conformation was identified as ac atalytically active pose and suitable for transition state formation when fulfilling the criteria summarized in Ta ble 7.
For BM3 M11, chain Bo ft he crystal structure of the heme domain of mutant BM3 M11( PDB ID:5 E9Z) [32] was used as template for docking and subsequent MD simulations. Missing residue Q73 and missing atoms of residues K31, Q73, K94, K97, Q109, Q110, D136, K187, K218, Q229, T245, R255, Q288, K306, K449 were added with Modeller 9.3. [69] For docking and MD simulations with CYP1A2, the crystal structure from PDB ID:2 HI4 was used. [70] To obtain proteinbinding poses for MF and TE to start MD simulations from, they were docked into the protein templates (using the PLANTS docking software, version 1.2 [71] and the ChemPLP scoring function [72] ) and equilibrated in MD simulations in which the heme group was described in its resting state (i.e.,w ith af erryl-oxygen dummy atom). Prior to docking, initial (steepest-descent) energy minimiza- Table 7. Geometric criteria to identify conformationsf rom moleculard ynamics simulations as catalytically active binding poses for hydroxylation of aromatic or aliphatic CÀHs iteso fm etabolism (SOMs). Unlessn oted otherwise, these criteria were derived from combined QM/MMs tudies of Mulholland and co-workers. [63] Type of hydroxylation