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

Article abstract of DOI:10.1002/cbic.201900751

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.

Full text of DOI:10.1002/cbic.201900751

Accepted Article
Title: A modified Arrhenius approach to thermodynamically study
regioselectivity in Cytochrome P450 catalyzed substrate
conversion
Authors: Rosa A. Luirink, Marlies C.A. Verkade-Vreeker, Jan N.M.
Commandeur, and Daan Geerke
This manuscript has been accepted after peer review and appears as an
Accepted Article online prior to editing, proofing, and formal publication
of the final Version of Record (VoR). This work is currently citable by
using the Digital Object Identifier (DOI) given below. The VoR will be
published online in Early View as soon as possible and may be different
to this Accepted Article as a result of editing. Readers should obtain
the VoR from the journal website shown below when it is published
to ensure accuracy of information. The authors are responsible for the
content of this Accepted Article.
To be cited as: ChemBioChem 10.1002/cbic.201900751
Link to VoR: http://dx.doi.org/10.1002/cbic.201900751
A Journal of
www.chembiochem.org

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
A modified Arrhenius approach to thermodynamically study
regioselectivity in Cytochrome P450 catalyzed substrate
conversion
Rosa A. Luirink[a] [b], Marlies C.A. Verkade-Vreeker[a] [b], Jan N.M Commandeur[a] and Daan P.
Geerke*[a]
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,
epoxidation and heteroatom oxidation. Figure 1 shows 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-electron reduction, (5) protonation of the
Fe2+OO-, (6) heterolytic cleavage of the hydroperoxyl bond to
yield FeO3+, (7) hydrogen abstraction of C-H-bond, (8) rebound
of hydroxy-group, and (9) release of the product. Which step is
rate limiting appears to depend on the specific combination of
P450 isoform and substrate involved in the hydroxylation
reaction.[21–30] For several substrates the rate-limiting nature of
hydrogen abstraction has been demonstrated by the kinetic
isotope effect (KIE) observed after deuterium-substitution.[23,27,30]
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. The accessibility of potential SOMs depends on the relative
binding free energy ∆∆Gbind 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 multi-step
P450-catalyzed processes which involve redox partners. We
introduce a modified Arrhenius approach to overcome the limitations
in studying P450 selectivity, which can be applied in multi-product
enzyme catalysis. Our approach gives combined information on
relative activation energies, ∆∆Gbind values and collision entropies,
yielding direct insights into the basis of selectivity in substrate
conversion.
Keywords: Cytochrome P450 catalysis
•
selectivity
•
thermodynamics • validation • Molecular Dynamics simulation
Introduction
Cytochrome P450 enzymes (P450s) form
a
superfamily of heme-containing proteins that play an important
role in the oxidative metabolism of many lipophilic xenobiotics,
as well as in the biosynthesis and catabolism of endogenous
compounds.[1] P450s are considered to be the catalytically most
diverse enzymes in nature and because of their versatility they
have many potential biotechnological applications.[2-15] So far,
sequences of over 300,000 isoforms have been determined in
all domains of life.[16] In humans, P450s are responsible for
approximately 75% of phase I metabolism 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.
Figure 1. Catalytic cycle for Cytochrome P450 catalyzed hydroxylation of
substrate RH.
[a]
[b]
R.A. Luirink, M.C.A Verkade-Vreeker, dr. J.N.M Commandeur, dr.
D.P. Geerke
AIMMS Division of Molecular Toxicology
Vrije Universiteit
De Boelelaan 1108, 1081 HZ Amsterdam, the Netherlands
E-mail: d.p.geerke@vu.nl
Contributed equally
It is generally accepted that the regio- and stereo-
selectivity of P450s are governed by (a) the preference and
probability of the substrate to bind in a reactive orientation
relative to the activated oxygen-species, and (b) the activation
energies of the specific oxidation reactions at the exposed sites-
Supporting information for this article is given via a link at the end of
the document.
1
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
of-metabolism (SOMs).[31] To predict and/or rationalize
regioselective metabolite formation, several in silico approaches
have been used to study orientations in substrate binding and/or
the activation energy Ea of the oxidation reactions involved.
Computational approaches to study the orientation(s) and
dynamics of substrate binding include docking methods,
molecular dynamics (MD) simulations, and/or binding free
energy computation.[13,32–34] Furthermore, quantum mechanical
methods such as density functional theory (DFT),[35–42] albeit in
combination with molecular mechanics (MM) techniques, have
been used to calculate values for Ea. There is still limited
experimental data of preferred binding modes and activation
energies to validate the predictivity of these computational
approaches. Possible binding orientations of substrates in the
active site of P450s have been studied experimentally using co-
crystallography[43–45] and spin-relaxation studies.[46–48] However,
in co-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 a productive 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
unusual for P450s.
Until now, several examples of Arrhenius plots of P450
catalyzed reactions have been reported.[18,54–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 b5 reductase (steps 2 and
4, Figure 1), may also be rate limiting for the overall reaction.[21]
Second, non-linear Arrhenius plots were found with reactions
catalyzed by microsomal P450s with
a discontinuity at
approximately 20 oC, which was attributed to a transition of
membrane fluidity affecting the interaction between P450 and
P450 reductase.[57,58] Furthermore, usually only a small range of
temperatures is used since above the optimal temperature, the
reaction rates decrease again due to enzyme denaturation.
Finally, P450 catalysis often leads to different products at
different ratios. Experimentally, kcat values are usually
determined by dividing values for the maximal velocity in
substrate conversion as obtained from enzyme kinetic studies
(Vmax) by the total enzyme concentration [E]total (kcat = Vmax
/
[E]total), assuming that at maximal enzyme activity all enzymes
are occupied by the substrate in a reactive binding pose.
However, in the case of parallel reactions the individual Vmax
values cannot be divided by [E]total but should be divided by the
concentration of the enzyme-substrate complexes of the
corresponding reactive binding poses, designated [ES1] and
[ES2] in Scheme 1 for possible formation of two products (P1 and
P2). No direct experimental methods are available to accurately
determine ratios of different bound conformations of a given
enzyme-substrate complex.
A
common experimental method to determine
activation energies of chemical and enzymatic reactions is the
quantification of the reaction rate constant for substrate-to-
product conversion (kcat) at different temperatures, and by
subsequently constructing logarithmic plots of (ln kcat) versus the
reciprocal absolute temperature (1/T) according to the Arrhenius
equation,[52,53]
þ
ċč
Ĕ
͟
Ɣ ̻͙
.
(1)
ꢀꢁ/
According to Equation (1) the slope of the linear plot obtained
equals the negative value of the activation energy Ea divided by
the gas constant R. The pre-exponential factor A comprises the
frequency or collision efficiency at which 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.
Scheme 1. Kinetic scheme for the catalytic conversion of substrate S to two
possible products, P1 and P2. The associated rate contants kcat,1 and kcat,2
depend on rate constants for the individual steps of the catalytic conversion,
which comprise binding to (kb) and unbinding of (ku) the enzyme-substrate
complex ES, interconversion between ES1 and ES2 (ki, k-i), formation of the
enzyme-product complex EP out of ES (kp), and unbinding of the EP complex
(kr).
2
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
In the present study we propose and evaluate a
modified Arrhenius approach to overcome several of the above-
mentioned intrinsic limitations of thermodynamic studies on
P450 catalyzed reactions. In this modified approach the
temperature dependence of the ratio of Vmax 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., of reactions 1 and 2 with maximal velocities
Vmax,1 and Vmax,2 and rate constants kcat,1 and kcat,2, respectively)
and using Vmax= kcat • [ES]max gives:
or
ln ƴꢂꢂĠĔī,uƸ Ɣ ln ʠꢅ ʡ Ǝ ꢇꢈ
(5)
ꢅ
∆
u
ĠĔī,v
v
where ∆ = ∆∆Gbind + ∆Ea , and ∆Ea = Ea,1 – Ea,2
.
Because we assume that the rates and temperature
dependence of steps 2-6 of the catalytic cycle in Figure 1 will be
similar for the parallel pathways, these factors will cancel out
when evaluating the ratios of product formation and therefore
will not contribute to differences 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 Vmax values of both pathways to a similar extent, which will
allow studies over a larger temperature range.
þ
þ
ċč
Ĕ,u Ĕ,v
ꢂ
&
&
∙ʞꢃꢄ ʟ
u ĠĔī
ʞꢃꢄ ʟ
u ĠĔī
∙ꢅ
∙ꢅ
ĠĔī,u
ĖĔħ,u
u
ͯ
Ɣ
Ɣ ʞꢃꢄ ʟ
∙ ͙
(2)
ꢂ
∙ʞꢃꢄ ʟ
v ĠĔī
ĠĔī,v
ĖĔħ,v
v ĠĔī
v
The subscript max for concentrations of the enzyme-substrate
complexes indicates ES, ES1 or ES2 concentrations at maximal
enzyme activity. Under a steady-state approximation the ratio
between [ES1] and [ES2] can be related to
A
∆∆Ā
ċč
ĕĜġė
ʞꢃꢄ ʟ
u ĠĔī
ͯ
Ɣ ͙
(3)
ʞꢃꢄ ʟ
v ĠĔī
with ∆∆Gbind
=
∆Gbind,1
–
∆Gbind,2, which is the (possible)
difference in binding free energies ∆Gbind between substrate
binding poses associated with formation of products 1 and 2,
respectively. The steady-state approximation of Equation (3) is
valid when substrate concentration [S] >> [E]total, and when the
constants of substrate binding and unbinding (kb and ku in
Scheme 1) are substantially higher than those for formation of
the enzyme-product (kp,1 and kp,2 in Scheme 1) and/or
alternatively, when rapid interconversion between ES1 and ES2
is possible (with ki and k-i in Scheme 1 being higher than the kp
values).
B
Hence, when correlating the natural logarithm of the ratio
Vmax,1/Vmax,2 with the inverse absolute temperature (in a modified
Arrhenius plot), a straight line is expected of which the slope of
the plot will represent the sum (∆) of the differences in activation
energy (Ea) and ∆Gbind of the parallel reactions (divided by the
gas constant R), Equations (4) and (5).
Figure 2. Gibbs free energy profile for a catalyzed reaction with two possible
products with A) a binding free energy difference and B) without a binding free
energy difference between the binding poses leading to the different products
ln ƴꢂꢂĠĔī,uƸ Ɣ ln ʠꢅ ʡ Ǝ
(4)
ꢅ
∆∆ꢆ
ͮ ∆ꢃ
Ĕ
u
ĕĜġė
ꢇꢈ
P1 and P2, starting from enzyme and substrate (E + S). ∆G‡ comprises the
,x
ĠĔī,v
v
sum of an activation energy (Ea,x) and a collision entropy (T∆S‡,x) term.
3
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
In previous studies, ratios of product formation
obtained at one incubation temperature have been used to
estimate the overall difference in activation free energies
(∆∆Goverall) for competing reactions by P450s, Figure 2.[25,60] Also
in this case, possible rate limiting factors prior to the oxidation
reaction (steps 2-6, Figure 1) are expected to cancel out.
According to the Curtin-Hammett principle,[61] the product ratio of
two competing reactions is governed by the difference ∆∆Goverall
between the free energies of the corresponding transition states
such cases to be used instead of Equation (7). The difference in
employing our modified Arrhenius approach compared to direct
use of the Curtin-Hammett formalism is that the entropic
contribution to differences in the barrier for forming the transition
state from the enzyme-substrate complex (i.e., the ratio A1/A2 in
Equation (5)) can be separated from other contributions (∆∆Gbind
and ∆Ea). Therefore, our approach can be of direct help in
validating (the combined use) of free energy, quantum chemical,
and MD studies on preferred modes of substrate binding,
activation energies, and/or probabilities to adopt catalytically-
active binding orientations, respectively.
‡
‡
([ES-O]1 and [ES-O]2 in Figure 2) when the barrier to
interconversion between reactive binding poses ES1 and ES2
(either direct or via enzyme/substrate unbinding) is much smaller
than the barrier to product formation. The free energies of
forming these transition states depend on the one hand on
differences in free energies of binding of the productive binding
poses ES1 and ES2 (∆∆Gbind) and secondly by the differences
between free energies of the binding poses and the transition
In the present study the human P450 isoform 1A2
(CYP1A2) and a drug metabolizing mutant of bacterial P450
BM3 (CYP102A1), i.e., BM3 M11, are used to evaluate the
applicability and illustrate the value of our approach to analyze
thermodynamic determinants of selectivity in P450-catalyzed
product formation. For that purpose, we determined the
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 oxidized by CYP1A2 and P450 BM3 M11 to
two or three metabolites, respectively, while TE conversion
catalyzed by BM3 M11 leads to three different products as well,
Figure 3.[62][13] Recombinant CYP1A2 was selected as model for
a membrane-bound P450, which depends on co-expressed
NADPH cytochrome P450 oxidoreductase as redox partner.
P450 BM3 M11 was used as model for a soluble P450. Wildtype
P450 BM3 is a natural fusion protein between a P450-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 M11 was developed by
combination of site-directed and random mutagenesis and
‡
‡
states involved in product formation, with ∆∆G‡ = ∆G1 -∆G2 ,
Figure 2A.
‡
ʞΞΠΝΒΣΑ΢ ʟ
ͮ ∆∆ͻ
u ʡ Ɣ − ∆∆ͻ-´£°Ÿªª Ɣ − ∆∆ͻ
(6)
§¬¢
ln ʠ
ʞΞΠΝΒΣΑ΢ ʟ
ΆΈ
ΆΈ
v
Previously, Higgins et al. assumed that when different human
P450s show both similar kinetic isotope effects and product
ratios, the product ratios observed are determined by the
differences in activation free energy of the transition states
∆∆G‡:[60]
‡
∆∆ꢆ
ꢇꢈ
ln ʠʞʞΞΞΠΠΝΝΒΒΣΣΑΑ΢΢ ʟʟʡ Ɩ Ǝ
(7)
u
v
This may well be valid for small substrates that can rapidly adopt
multiple binding poses (high ki and k-i in Scheme 1) and that
have only a small difference in free energy of binding (∆∆Gbind),
cf. Figure 2B. In such cases, ∆∆G‡ can be directly estimated
from ratios in product formation using Equation (7). In case of
substrates with high molecular weight and/or P450s with
restrictive active sites, next to differences in free energies of
activation, also steric factors may play an important role in the
regioselectivity of P450 reactions, and ∆∆Gbind in Equation (6)
cannot be neglected a priori. This is exemplified e.g. by recent
binding free energy calculations[32] for one of the pairs of product
formation considered in the current work, and Equation (6) is in
catalyzes oxidation reactions of
a
wide variety of
pharmaceuticals and other chemicals.[6] To determine ∆ and
relative collision efficiencies ln(A1/A2) in Equation (5) for the
multiple substrate conversions catalyzed by CYP1A2 or P450
BM3 M11, enzyme kinetic parameters were determined for each
reaction at different incubation temperatures. In support of our
steady-state approximation in Equation (3), we also measured
4
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
Figure 3. Structures of the compounds used in this study. Arrows indicate the sites of metabolism as experimentally observed for P450 BM3 mutant M11.
Abbreviations refer to products formed by the regiospecific oxidation reactions. Note that mefenamic acid hydroxylation as catalyzed by P450 isoform 1A2 leads
to formation of 4’-OH-MF and 5-OH-MF only.
kinetic isotope effects for the pair of product formation (i.e.
conversion 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 formation based on combined QM/MM
studies by Mulholland and co-workers.[63] The results were
compared with the relative collision efficiencies as determined
from the intercepts of our Arrhenius plots, and with differences
between our estimated ∆ values and corresponding Curtin-
Hammett estimates for relative activation barriers, as measures
for possible differences in the entropy of transition state
formation. To further interpret and cross-validate our modified
Arrhenius and in silico analyses we also computed differences in
activation energies Ea (using the SMARTCyp web server)[64]
and/or obtained them from literature, and where possi-ble we
combined these estimates with ∆∆Gbind values reported in
literature for a direct comparison with our values for ∆.
minor product, Figure 4A. As summarized in Table 1, the
catalytic efficiency (Vmax/Km) for all three pathways increased
from 4 oC to 25 oC. At higher temperatures the catalytic
efficiency decreased again. Also for the Vmax-values the lowest
values were obtained at the lowest and highest incubation
temperatures. As a result, the plots of ln(Vmax) versus 1000/T
were strongly nonlinear when analyzing the kinetics of each
metabolite individually (data not shown). This was expected
based on previous Arrhenius studies on P450 catalyzed
reactions.[56–59]
Figure 4. Extracted ion chromatograms of A) non-deuterated (d0) and B)
deuterated (d12
) hydroxy metabolites of mefenamic acid, formed after
Results and Discussion
incubations of P450 BM3 M11 with a 52/48 mixture of deuterated (d13) and
non-deuterated mefenamic acid (total concentration 750 µM), respectively.
Temperature dependent mefenamic acid hydroxylation
catalyzed by BM3 M11
Before applying our modified Arrhenius analysis to the
(three) pairs of BM3 M11 catalyzed product formations (Table 2),
we measured relative kinetic isotope effects in an attempt to
explicitly verify the steady-state approximation taken in Equation
(3) for the ratio of MF conversion to 3’-OH-MF and 4’-OH-MF.
Representative extracted ion chromatograms of the mixed
deuterated and non-deuterated hydroxy metabolites in Figure 4B
As described previously,[33] mefenamic acid was
metabolized by P450 BM3 M11 to the 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 5-hydroxymefenamic acid (5-OH-MF) as relatively
5
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
Table 1. Temperature dependence of enzyme kinetic parameters for the regioselective hydroxylation of mefenamic acid by P450 BM3 M11 and recombinant
human CYP1A2 (i.e., for the formation of the products 3’-hydroxymethylmefenamic acid, 4'-hydroxymefenamic acid and 5-hydroxymefenamic acid).
3'-hydroxymethylmefenamic acid
4'-hydroxymefenamic acid
5-hydroxymefenamic acid
a
b
a
b
a
b
T (K)
Km (µM)
Vmax
Vmax /Km
Km (µM)
Vmax
Vmax /Km
Km (µM)
Vmax
Vmax /Km
P450 BM3 M11
277.0
283.4
290.2
292.7
298.6
304.4
313.3
317.6
72 16
135 10
1875
1514
3484
3434
5447
2281
644
76 17
230 19
343 12
859 82
736 52
608 50
787 41
431 37
208 14
3026
2271
5053
4628
6909
2223
714
77 22
24
34
85
72
64
81
42
21
3
2
8
4
7
7
5
2
312
274
574
383
744
198
109
27.6
140 18
161 14
145 45
85 20
266 56
556 200
388 184
212
9
151 15
170 17
159 52
88 23
354 34
603 150
538 34
124 19
148 15
188 37
86 25
409 165
384 202
762 118
561 44
498 33
463 33
607 53
358 28
194 19
500
387
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
CYP1A2
279.4
287.7
299.4
309.8
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
N.A.
204 118
210 92
144 68
182 33
0.56 0.13 2.75
0.99 0.17 4.71
0.91 0.18 6.88
6
2
0.50 0.02 83.3
0.78 0.09 16.6
0.67 0.03 16.0
0.68 0.02 14.5
47 21
42
47
7
4
1.0 0.09
5.49
a.
b.
Unit: nmol product/min/nmol enzyme;
Unit: µL/min/nmol enzyme.
N.A. Not applicable.
Table 2. Comparison of thermodynamic properties for the oxidation reactions of mefenamic acid and testosterone catalyzed by cytochrome P450 BM3 M11 and
recombinant human CYP1A2, as determined using our modified Arrhenius approach (Equation (5)), the Curtin-Hammett principle (Equation (6)), or SMARTCyp
(to calculate differences in activiation energies, ∆Ea).
Comparison
for the
formation of:
ln(Vmax,1/Vmax,2) vs. 1000/T
Modified
Arrhenius
∆ (kJ mol-1)
Curtin-Hammett SMARTCyp
Enzyme
∆∆Goverall
∆Ea
(kJ mol-1) a
Intercept
Slope (K) R2
(kJ mol-1)
Mefenamic Acid
3’-OH-MF vs. 4’-OH-MF
3’-OH-MF vs. 5-OH-MF
4’-OH-MF vs. 5-OH-MF
BM3 M11
BM3 M11
BM3 M11
3.0 0.2
5.3 0.2
2.3 0.2
-0.97 0.05 0.97
-0.97 0.07 0.98
0.00 0.07 0.14
8.1 0.5
8.1 0.6
0.0 0.6
0.7
-4.9
-5.6
- 1.8
- 1.8
0
4’-OH-MF vs. 5-OH-MF
CYP1A2
2.8 0.3
-0.74 0.10 0.97
6.1 0.8
-0.8
0
Testosterone
15ß-OH-T vs. 16ß-OH-T BM3 M11
2ß-OH-T vs. 16ß-OH-T BM3 M11
2ß-OH-T vs. 15ß-OH-T BM3 M11
7.6 0.5
9.1 0.5
1.4 1.0
-1.95 0.14 0.99
-2.82 0.16 0.99
-0.87 0.29 0.82
16.7 1.1
23.5 1.3
7.2 2.5
-2.8
0.8
3.6
0
-9.5
-9.5
a.
SMARTCyp activation energies: 3’-OH-MF, 66.4 kJ mol-1; 4’-OH-MF, 68.2 kJ mol-1; 5-OH-MF, 68.2 kJ mol-1; 2β -OH-T, 66.4 kJ mol-1; 15β -OH-T, 75.9
kJ mol-1; 16β -OH-T, 75.9 kJ mol-1.
show that the relative peak area of the 3’-OH-MF product (as
compared to the peak areas for 4’-OH-MF or 5-OH-MF)
decreases significantly when going from non-deuterated to
deuteration may be well explained by metabolic switching
resulting from the strongly decreased 3’-hydroxylation. Thus, our
kinetic isotope effect measurements for these reactions suggest
that rapid interconversion between catalytically active poses for
3'-OH-MF and 4'-OH-MF formation is possible and accordingly,
that ki and k-i values for binding-pose interchange and/or kb and
ku values are probably higher than kp for the corresponding
product formations. Therefore, even in the unexpected case that
kbind << kcat, the steady-state approximation used in Equations
(3)-(5) can still be assumed to be valid.
deuterated mefenamic acid. Table
3
summarizes that
deuteration of mefenamic acid results in an approximately 4-fold
reduction in the rate of 3’-methyl hydroxylation but to a 50%
increase of the 4’-hydroxylation pathway, and no change in 5-OH
mefenamic acid formation. These results indicate that the
hydrogen-abstraction of the 3’-methyl-group is rate-limiting to a
more significant extent than aromatic hydroxylation reactions.
This is in line with previous studies in which substantially larger
KIEs were observed for aliphatic than for aromatic hydroxylation
by P450s.[30] The increase in 4’-OH-MF formation after
6
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
Table 3. Kinetic isotope effects for hydroxylation of mefenamic acid by P450
BM3 M11: ratios of nondeuterated (d0) vs deuterated (d12) hydroxy metabolites
formed in incubations of P450 BM3 M11 with a 52/48% mixture of deuterated
and nondeuterated mefenamic acid (total concentration 750 µM).
Kinetic isotope effect (ratio d0 vs d12 OH-MF)
Temperature (K)
3’-OH-MF
4’-OH-MF
5-OH-MF
278
288
298
308
3.84 0.20
4.35 0.20
3.78 0.48
3.58 0.55
0.65 0.002
0.66 0.001
0.65 0.002
0.79 0.006
1.08 0.02
1.02 0.03
1.06 0.01
1.00 0.01
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 Vmax
values of mefenamic acid metabolites against 1000/T, linear
curves were obtained, Figure 5A. This observed correlation
supports our assumptions that the temperature dependence of
steps 2-6 of the catalytic cycle in Figure 1 will be similar for the
different hydroxylation paths and cancel when relating Vmax
ratios to the inverse temperature. The slopes (∆ = ∆∆Gbind + ∆Ea)
and intercepts of the thus obtained modified Arrhenius plots are
summarized in Table
2
and discussed below in our
thermodynamic analysis of the observed regio-specificity in MF
conversion. In addition, ∆∆Goverall values as derived (using the
Curtin-Hammett principle and Equation (6)) from pairs of Vmax
values at 300K are also reported in Table 2, as well as
differences in activation energies ∆Ea calculated using
SMARTCyp. Note that we used two versions of SMARTCyp
(versions 2 and 3), which gave identical results for Ea. The small
differences in activation barrier for all three hydroxylation
reactions are in line with the similar values for Ea calculated at
the B3LYP level of Density Functional Theory (DFT) by Leth et
al.[65] Leth modeled Compound I as a porphyrin moiety without
Figure 5. Modified Arrhenius plots of product ratios of mefenamic acid formed
by A) P450 BM3 M11 and B) recombinant human CYP1A2, for the ratio
ln(Vmax,4’OH-MF/Vmax,5-OH-MF) (triangles), ln(Vmax,3’OH-MF/Vmax,5-OH-MF) (squares), and
ln(Vmax,3’OH-MF/Vmax,4’-OH-MF) (circles).
−
side chains and with axial coordinating O2 and CH3S− ligands,
and showed e.g. a difference in the range of −6 to 3 kJ mol-1 in
activation barrier when comparing 3’-methyl-OH-MF and 4’-OH-
MF formation. This is to be compared with the corresponding
SMARTCyp value of −2 kJ mol-1, Table 2.
Despite the lower binding affinity of MF to BM3 M11 in
its pose that is catalytically active for 4’-OH-MF formation, the
difference with 3’-methyl hydroxylation in the overall activation
free energy ∆∆Goverall is close to zero from our Curtin-Hammett
analysis (0.7 kJ mol-1, Table 2). Thus, the difference in ∆Gbind is
for a large part counterbalanced by a lower entropy penalty to
form the transition state for 3'-methyl-OH-MF formation out of
the corresponding ES complex, as reflected by the higher
collision efficiency for 3’-methylhydroxylation (with a value of 3.0
+/- 0.2 for the intercept of the modified Arrhenius plot, Table 2).
We could cross-validate these findings with MD simulations of
BM3 M11 in complex with MF, in which we compared the
frequencies of occurrence of MF binding poses that can
potentially adopt transition state geometries for either 3’-Me or 4’
hydroxylation. Indeed, our simulations showed higher
frequencies for substrate orientations that are in line with the
For the ratio between 3'-methyl-OH-MF and 4'-OH-MF
formation the modified Arrhenius plot (Figure 5A) shows a slope
of -0.97 ± 0.05 K which corresponds to a value for ∆ (in Equation
(5)) of 8.1 ± 0.5 kJ/mol, with a lower sum of ∆Gbind and Ea for
hydroxylation at the 4’ aromatic SOM, Table 2. As stated
previously, ∆Ea is close to 0 (and even slightly negative),
therefore the positive ∆ value should be interpreted to result
from a more negative (favorable) binding free energy for the
catalytic binding orientation for 4’ hydroxylation of mefenamic
acid as compared to 3’-methyl hydroxylation (∆∆Gbind = ∆ − ∆Ea =
8.1 − − 1.8 = 9.9 kJ mol-1), in agreement with and confirming the
corresponding ∆∆Gbind values previously computed by us that
range between 9.3 and 11.6 kJ mol-1.[32]
7
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
simulations with mefenamic acid in the 3’methyl hydroxylation
pose (Figure S1), which was not observed in the 5-hydroxylation
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 M11 is overall favorable over 3’methyl hydroxylation,
despite the higher entropic cost for transition state formation and
the slightly higher activation energy, should thus be understood
in terms of the lower binding affinity for substrate binding in a
pose that leads to 4’-OH-MF formation. Our previous detailed
free energy perturbation computation study[32] on selectivity in
MF hydroxylation came to this conclusion as well and could thus
be verified with our modified Arrhenius approach.
pose (Figure S2). Furthermore, we indeed observe
a
substantially higher frequency in MD simulations of substrate
orientations corresponding to transition state formation for the 3’-
methyl hydroxylation compared to 5 hydroxylation (Table 4).
For the ratio between 4’-OH-MF and 5-OH-MF
formation by BM3 M11, the modified Arrhenius plot (Figure 5A)
shows a slope of 0.00 ± 0.07 K which corresponds to a value for
∆ (in Equation (5)) of 0.0 ± 0.6 kJ mol-1, Table 2. The predicted
identical Ea values for both pathways (Table 2) thus suggest a
similar binding free energy for the corresponding catalytically
active poses (∆∆Gbind = ∆−∆Εa) whereas the lower B3LYP-D3
value of Leth et al. (by 11 kJ mol-1) [65] hints at preferred binding
in the pose enabling 4’ hydroxylation. The preference of 4’ over 5
hydroxylation can again be understood in terms of the higher
probability (lower entropic cost) of transition state formation for 4’
hydroxylation compared to 5-OH product formation (Table 2),
probably also due to hydrogen bonding with Serine 72 in the
catalytic-active binding pose in the former case (Figure S3).
For the ratio between 3'-methyl-OH-MF and 5-OH-MF
formation the modified Arrhenius plot (Figure 5A) shows a slope
of -0.97 ± 0.07 K which corresponds to a value for ∆ (in Equation
(5)) of 8.1 ± 0.6 kJ/mol, with a lower sum of ∆Gbind and Ea for
hydroxylation at the 5 aromatic SOM, Table 2. When considering
the similar activation energy for 5-OH-MF as predicted by
SMARTCyp (with ∆Ea = −1.8 kJ mol-1, Table 2), the binding free
energy should be higher for the 3’methyl binding pose,
suggesting that binding in an orientation that can lead to 5-OH-
MF formation is more favorable. As an alternative explanation,
the difference in ∆ may be due to a lower Ea value for 5’-
hydroxylation, as indicated 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 Ea value for 5-OH-MF formation by 14 kJ mol-1.[65] In any
case, from Table 2 the preference of 3’-methyl over 5-OH
mefenamic acid formation by P450 BM3 M11 can be understood
in terms of the higher probability (lower entropic cost) of
transition state formation for 3’-methyl hydroxylation, as
indicated by the higher intercept of our modified Arrhenius plots
and the observed difference between ∆ and the Curtin–Hammett
estimate for ∆∆Goverall. This observed difference is equal to
T∆∆S‡ under the assumption that trends in enthalpy and energy
are equal, Cf Equation (8) and Figure 2.
Table 4.
Percentages of substrate binding orientations during two
independent MD simulations of mefenamic acid bound to BM3 M11 (A and B)
that are suitable for transition state formation for hydroxylation of mefenamic
acid at its 3’methyl, 4’ or 5 position.
3'methyl
79.2
4'
66.1
46.7
5
29.0
16.0
A
B
76.6
Table 5.
Percentages of substrate binding orientations during two
independent MD simulations of mefenamic acid bound to CYP1A2 (A and B)
that are suitable for transition state formation for hydroxylation of mefenamic
acid at its 4’- or 5 position.
4'
5
A
B
33.1
65.1
1.9
7.2
ʚ
ʛ
ʚ
‡ʛ
∆ Ǝ ∆∆GΝΤΓΠΏΚΚ Ɣ ∆∆GΐΗΜΒ ƍ ∆EΏ Ǝ ∆∆GΐΗΜΒ ƍ ∆∆G
Ɣ ∆EΏ Ǝ ʚ∆∆H‡ Ǝ ͎∆∆S‡ʛ ≅ ͎∆∆S‡
(8)
Table 6.
Percentages of substrate binding orientations during three
This is in accord with our in silico data from references [32] and
[33] showing strong hydrogen bonding interaction between
mefenamic acid’s carboxylate group and BM3 M11’s Ser72
residue, which directs substrate-binding orientations for 3’-
methyl hydroxylation to adopt a catalytically active pose. Such
an anchoring hydrogen bond is not present when bound in a
pose enabling hydroxylation at the 5 position. Also in our MD
simulations this particular hydrogen bond was observed in the
independent MD simulations of testosterone bound to BM3 M11 (A, B and C)
that are suitable for transition state formation for hydroxylation of testosterone
at its 15β or 16β position.
15 β
16 β
A
B
C
49.7
4.5
49.7
40.3
16.2
53.1
8
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
Temperature dependent mefenamic acid hydroxylation
catalyzed by CYP1A2
Oxidation of mefenamic acid by recombinant CYP1A2
resulted in formation of 4’-OH-MF and 5-OH-MF. At lower
substrate concentrations 5-OH-MF was the major metabolite, as
indicated by the higher Vmax/Km-values, Table 1. Vmax-values of
the 4’-hydroxylation pathway were slightly higher than for 5-
hydroxylation. Non-linear Arrhenius plots are obtained when
plotting ln(Vmax) versus 1000/T for the metabolites individually
(data not shown). As for the pairs of BM3 M11 mediated product
formations, the modified Arrhenius plot of the ratio of Vmax values
for the CYP1A2 catalyzed pathways as plotted against 1000/T
showed linear behavior, Figure 5B.
Figure 6.
Modified Arrhenius plots of product ratios for testosterone
conversion by BM3 M11, for the ratios ln(Vmax,15β-OH-T/Vmax, 16β-OH-T) (triangles),
ln(Vmax,2β-OH-T/Vmax, 16β-OH-T) (squares), and ln(Vmax, 2β-OH-T/Vmax, 15β-OH-T) (circles).
A higher collision frequency for formation of the
transition state for 4’ hydroxylation (compared to 5-hydroxylation,
cf. the positive intercept and positive difference between ∆ and
∆∆Goverall in Table 2) suggests an entropically more favorable
transition state formation for 4’-OH-MF formation. This is line
with the higher frequency observed in MD of substrate
orientations corresponding with transition state formation for 4’
hydroxylation (Table 5). During the simulations, we observe a
hydrogen bond between mefenamic acid and Thr469 when
binding in the 5-hydroxylation position (Figure S4), whereas we
do not observe any stabilizing or positioning hydrogen bonds in
the 4’-hydroxylation binding pose (Figure S5). The slight
preference of 5-OH mefenamic acid formation by CYP1A2 (as
reflected by the negative ∆∆Goverall estimate in Table 2) can in
this case be associated to a lower binding free energy for and
preferred binding in the corresponding catalytically-active pose.
Furthermore, the Ea value for 5’-hydroxylation may be lower (see
above).
Using the slopes of these curves, ∆ values for the
three different pairs of product formation were obtained, Table 2.
The minor pathway leading to 2β-OH-T was found to have the
highest sum of ∆Gbind and Ea, with ∆ = 23.5 kJmol-1 for the ratio
with 16β-OH-T, and 7.2 kJ/mol for the ratio with 15β-OH-T. Thus,
∆ is significantly higher for 2β-OH-T when compared to the other
products, whereas the activation energy for formation of the
corresponding transition state is lowest. This is apparent from
the value for Ea, which was predicted by SMARTCyp to be 9.5
kJ mol-1 lower than for 15β and 16β hydroxylation (Table 2),
probably due to the adjacent C=O moiety beside of the C2
carbon. From previous comparative DFT calculations on the
energy of C-H bond breaking we even found differences of more
than 20 kJ mol-1 in favor of C-H bond activation at position 2
compared to positions 15 and 16.[66] The higher ∆ and lower Ea
values for 2β-OH-T imply a significantly higher binding free
energy (i.e., lower binding affinity) for TE binding in a pose
compatible with 2β 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 M11 with TE bound
in this binding pose.
Temperature
dependent
testosterone
hydroxylation
catalyzed by P450 BM3 M11
As illustrated in Figure 3, testosterone can be hydroxylated by
P450 BM3 M11 at three positions, leading in order of their
relative amounts to 15β-OH-T, 16β-OH-T or 2β-OH-T formation,
respectively.[46,62] The enzyme kinetic parameters of the
reactions performed at temperatures ranging from 6 to 36 oC are
shown in supplemental Table S1. As was observed for
mefenamic acid, no linear Arrhenius plots are obtained when
plotting ln(Vmax) versus 1000/T for the metabolites individually
(data not shown). Encouragingly, the modified Arrhenius plots of
the ratios of Vmax values versus 1000/T showed again linear
behavior for all three combinations of pathways, Figure 6.
The ∆ value for the ratio of formation of the major
metabolite 15β-OH-T and the less abundant 16β-OH-T is 16.7
kJ mol-1, Table 2. Ea is probably higher for 15β hydroxylation
than for 16β hydroxylation due to the substrate’s hydroxyl group
at C17, which was also indicated by our previous finding that the
energy cost of breaking the C15-H and C16-H aliphatic bonds is 5
kJ mol-1 higher for the former.[66] In this particular case, the
difference in Ea may be the only contribution to ∆ as our MD
simulations suggest that interconversion between binding
orientations suited for 15β- and 16β-hydroxylation can occur on
9
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
the ns time scale. This is illustrated in Table 6 which shows that
the geometric criteria for adopting the transition states for the
15β and 16β hydroxylation pathways can be both fulfilled during
a single simulation. Therefore, information on the frequency of
MD configurations consistent with transition state formation was
for both positions taken from the same set of (three) simulations.
The fact that 15β-OH-T formation is prevalent over 16β
hydroxylation can be explained in terms of a higher collision
frequency and entropy of transition state formation out of the
enzyme-substrate complex, as can be seen from the relatively
large intercept in the modified Arrhenius plot for the ratio of
15β/16β hydroxylation and from the according difference
between ∆ and ∆∆Goverall, Table 2. These results are in line with
the higher frequency observed in MD of substrate orientations
corresponding with transition state formation for 15β
hydroxylation (Table 6).
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 M11
was constructed as described earlier.[6] All other chemicals were
of analytical grade and obtained from standard suppliers.
Methods
Expression of Cytochrome P450 BM3 M11 and Cytochrome
P450 1A2
His-tagged Cytochrome P450 BM3 M11 was
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]
Conclusions
We presented a method that makes it possible to obtain
experimental estimates for thermodynamic determinants of
regio- (and/or stereo-)selectivity in P450 catalyzed substrate
conversion, by means of studying temperature dependent ratios
of pairs of metabolite formation as catalyzed by a single P450
isoform. We illustrated the use of this modified Arrhenius
approach by studying the determinants of the regioselectivity in
mefenamic acid and testosterone hydroxylation by P450 BM3
M11 and CYP1A2. For the selected P450-substrate
combinations, our approach gave insight into the basis of
selectivity by giving combined information on relative activation
energies, ∆∆Gbind values, and collision entropy differences. We
cross validated the observed collision entropy differences with
molecular dynamics simulations, and we were able to verify
previous computational free energy calculations. The obtained
agreement suggests that the presented method can also be
applied to other combinations of P450 isoforms and substrates
that involve formation of two or more products. The methods and
experiments described here are useful tools for future research
on regio- and stereo-selectivity of P450 catalysis to ultimately
improve biocatalysts, and the data that can be obtained with our
method allow to validate results from computational models to
understand and predict selectivity.
A bicistronic plasmid containing the cDNA of human
CYP1A2 cDNA and human NADPH Cytochrome P450
reductase was transformed into Escherichia coli strain DH5α. A
300 mL terrific broth (TB) supplemented with
1 mM δ-
aminolevulinic acid, 0.5 mM thiamine, 400 µL/L trace elements,
100 µg/mL ampicillin, 1 mM isopropyl-β-d-thiogalactopyranoside
(IPTG), 0.5 mM FeCl3 was inoculated with a 7,5 mL pre-culture
grown from a single colony. The cells were allowed to grow for
40 h at 28 °C and 125 rpm. Escherichia coli-cells were collected
by centrifugation (4000 × g, 4 °C, 15 min) and resuspended in
20 mL 0.1 M potassium phosphate buffer at pH 7.4 containing
20% glycerol, v/v, 0,25 mM ethylenediaminetetraacetic acid
(EDTA) and 0,1 mM dithiothreitol (DTT)). The cells were treated
with 0.5 mg/mL lysozyme for one hour at 4°C and 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 and 4 °C.
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°C until use.
Assessment of enzyme kinetic parameters at different
temperatures
All incubations were performed at enzyme
concentration and incubation times for which the product
10
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
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 5 mM
MgCl2 and 2 mM EDTA. Seven substrate concentrations ranging
from 10 to 750 µM were used in a final incubation volume of 100
µL. Reaction mixtures were preincubated in a shaking water
bath set at the incubation temperature for 10 minutes. Reactions
were initiated by addition of 10% (v/v) of a prewarmed solution
containing NADPH regenerating system; final concentrations
were 0.5 mM NADPH, 10 mM glucose 6-phosphate, and 0.4
U/mL glucose-6-phosphate dehydrogenase. The reaction was
allowed to proceed for 4 minutes at different temperatures and
then stopped by the addition of 100 µL ice-cold methanol. The
denatured enzyme fractions were precipitated by centrifugation
for 20 minutes at 14,000 rpm. The supernatants were isolated
and analyzed by HPLC or LC-MS as described below.
localized at
a
position not corresponding to the SOMs.
Therefore, the metabolic incubations were performed using the
mixture of d13-mefenamic acid and d12-mefenamic acid.
Because the deuterated compound contained low levels of
unidentified side-products, the kinetic isotope effects were
determined in a competitive intermolecular experiment with
equimolar mixtures of mefenamic acid and d12,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 5 mM MgCl2 and 2 mM EDTA.
Total substrate concentrations of the mixture of labelled and
unlabeled mefenamic acid were 75 and 750 µM. Reaction
mixtures (total volume 100 µL) were pre-warmed at the
incubation temperature for 10 minutes, before initiating the
reaction by addition of a NADPH regenerating system (final
concentrations of 0.5 mM NADPH, 10 mM glucose 6-phosphate,
and 0.4 unit/mL glucose-6-phosphate dehydrogenase). The
reaction was allowed to proceed for 4 minutes at different
temperatures and then stopped by the addition of 100 µL ice-
cold methanol. The protein was removed by centrifugation for 20
minutes at 14,000 rpm. The supernatants were analyzed on an
Agilent 1200 series rapid resolution LC equipped with a TOF
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 nondeuterated 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 a single concentration, no Arrhenius-plots were
constructed since activities did not represent Vmax-values.
For mefenamic acid and M11 enzyme kinetic
parameters for the formation of the three metabolites were
determined at eight temperatures ranging from 4 °C to 45 °C, to
investigate the linearity in more detail. The other substrates were
incubated at four temperatures ranging from 4 to 35 °C. The
enzyme kinetic parameters Vmax and Km were determined by
nonlinear regression according to the Michaelis-Menten
equation using Graphpad Prism 7.0 software (GraphPad, San
Diego USA).
∆ and ln(A1/A2) in Equation (5) for competing enzyme
reactions were determined by plotting the logarithms of the ratio
of the Vmax values measured at different temperatures against
the inverse of the absolute temperature T. According to Equation
(5), the slope of this curve, when linear, corresponds to the sum
∆ of ∆∆Gbind and ∆Ea divided by the negative gas constant R
(8.3145 J mol-1 K-1). The slopes and intercepts of the modified
Arrhenius plots were determined by linear regression using the
Graphpad Prism 7.0 software (GraphPad, San Diego USA).
Determination of competitive intermolecular kinetic isotope
effects of hydroxylation of mefenamic acid
A method 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). Repeated cycles of
microwave-assisted H/D exchange in presence of platinum and
palladium catalysts resulted in a 52%/48% mixture of d13-
mefenamic acid and d12-mefenamic acid. Because all three
metabolites showed a 52%/48% ratios of d12- and d11-labelling
degree, it was concluded that the remaining hydrogen-atom was
Analytical methods
The analyses of metabolites were performed by
reversed-phase liquid chromatography using a Shimadzu HPLC
equipped with two LC- 20AD pumps, a SIL20AC autosampler
and a SPD20A UV detector. Lab 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
11
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
for the possible detrimental effect of hydrogen interposition on C-
H activation by the ferryl oxygen.[32] Thus, a given enzyme-
substrate conformation was identified as a catalytically active
pose and suitable for transition state formation when fulfilling the
criteria summarized in Table 7.
1200 series rapid resolution LC was used which was connected
to an Agilent 6230 time-of-flight mass spectrometer equipped
with an electrospray ionization (ESI) source operating in positive
ion mode. A capillary voltage of 3500 V was used, and nitrogen
was used both as drying gas (10 L/min) and nebulizing gas
(pressure 50 psig) at a constant gas temperature of 350 °C.
1000 MS spectra/s were acquired and analysis was performed
using Agilent MassHunter Qualitative analysis software (version
2.0).
Table 7. Geometric criteria to identify conformations from molecular dynamics
simulations as catalytically-active binding poses for hydroxylation of aromatic
or aliphatic C-H sites of metabolism (SOMs). Unless noted otherwise, these
criteria were derived from combined QM/MM studies of Mulholland and co-
workers.[63]
Type of
hydroxylation
Aromatic
Distance criteria
Angle criteria
For all compounds, a Luna 5µm C18 column (4.6x150
mm) was used as stationary phase and gradients were
constructed by using two mobile phases: eluent A (0.8 %
acetonitrile, 99% water, 0.2% formic acid) and eluent B (0.8%
water, 99% acetonitrile, 0.2% formic acid).
Distance carbon of site
of metabolism to ferryl
oxygen (CSOM-OFe) <
0.35 nm
Distance of hydrogen
of site of metabolism to
CSOM-HSOM-OFe angle
should not be
between 140 and 220
degrees[a]
Aliphatic
HSOM-OFe-Fe angle
should be between
110 and 130 degrees
ferryl oxygen (HSOM
-
OFe) < 0.35 nm
For the analysis of metabolites of mefenamic acid, the
first 5 min was isocratic at 40% eluent B. From 5 until 30 min,
the concentration of eluent B was increased linearly to 100%,
followed by linear decrease back to 40% between 30- and 30.5-
min. Isocratic re-equilibration at 40% eluent B was maintained
until 45 min. The flow rate was 0.5 mL/min. UV-VIS detection
was performed at 254 nm.
[a]Adapted from reference [32]
.
For BM3 M11, chain B of the crystal structure of the heme
domain of mutant BM3 M11 (PDB ID 5E9Z)[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 2HI4 was
used.[70] To obtain protein-binding 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. with a ferryl-oxygen dummy atom). Prior to docking,
initial (steepest-descent) energy minimization of MF and TE with
the MMFF94 force field was performed using MOE.[73] After
docking and MD with the heme modeled in the resting state,
enzyme-substrate conformations consistent with 3’methyl
hydroxylation (mefenamic acid in BM3 M11), 4’ hydroxylation
and 5 hydroxylation (mefenamic acid in BM3 M11 and CYP1A2),
and 16β hydroxylation (testosterone in BM3 M11) were selected
to start MD simulations from, in which the heme group was
modelled in its compound I state. For this purpose, two starting
poses for MD were selected per combination of mutant and
product formation using the geometric criteria in Table 7.[32] For
MD simulations of testosterone in its catalytically-active pose for
15β/16β-OH-testosterone formation, a single starting pose was
selected from which three independent MD simulations were
started. All MD simulations (including thermal equilibration and
100 ns production simulations) were performed using identical
For the analysis of metabolites of testosterone, the first
1 min was isocratic at 50% eluent B. From 1 to 20 min the
percentage of eluent B was increased linearly to 99%; and from
20 to 20.5 min linearly decreased to 50% B and maintained at
50% for re-equilibration until 30 min. The flow rate was 0.5
mL/min. UV-VIS detection was performed at 254 nm.
Molecular dynamics (MD) simulations of substrate binding
to P450 BM3 M11 and CYP1A2
MD simulations 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 M11 or CYP1A2. To define occurrence of binding poses
that are suitable for hydroxylation of either an aliphatic or
aromatic C-H moiety, we used 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 a rule to exclude assignment of conformations to
be catalytically active for aromatic hydroxylation, in case the
angle between the C-H site-of-metabolism bond and the vector
connecting the corresponding hydrogen with the ferryl oxygen
was between 140 and 220 degrees (Table 7), in order to account
12
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
V. B. Urlacher, This J. is 2011, 47, 944–946.
simulation settings and force-field parameters for the protein and
heme group (either in the resting or the compound I state) as
described in reference [74]. Atomic coordinates were written out
to disk every 100 ps. Partial atomic charges of MF and TE for
use in MD were obtained with GAMESS (Version 1 May 2012)[75]
at the Hartree-Fock level using the 6-31G* basis set. Other
interaction parameters for MF and TE were used from the
General Amber Force Field for organic molecules version 1.7.[76]
After MD, protein-ligand interaction profiles during simulation
[10]
S. T. Jung, R. Lauchli, F. H. Arnold, Curr. Opin. Biotechnol.
2011, 22, 809–817.
[11]
[12]
L.-L. Wong, ChemBioChem 2011, 12, 2537–2539.
J. Kirchmair, M. J. Williamson, J. D. Tyzack, L. Tan, P. J.
Bond, A. Bender, R. C. Glen, 2012, DOI
10.1021/ci200542m.
[13]
[14]
H. Venkataraman, S. B. A. de Beer, L. A. H. van Bergen, N.
van Essen, D. P. Geerke, N. P. E. Vermeulen, J. N. M.
Commandeur, ChemBioChem 2012, 13, 520–523.
H. Venkataraman, M. W. den Braver, N. P. E. Vermeulen, J.
N. M. Commandeur, Chem. Res. Toxicol. 2014, 27, 2071–
2081.
were
analyzed
in
terms
of
protein
residue-ligand
interaction frequencies using the dedicated Python-based
biomolecular analysis library MDInteract, which is freely
available at https://github.com/MD-Studio/MDInteract.[77]
[15]
[16]
M. Tavanti, J. L. Porter, S. Sabatini, N. J. Turner, S. L.
Flitsch, ChemCatChem 2018, 10, 1042–1051.
D. R. Nelson, Biochim. Biophys. Acta - Proteins Proteomics
2018, 1866, 141–154.
Acknowledgements
[17]
[18]
[19]
F. P. Guengerich, Chem. Res. Toxicol. 2008, 21, 70–83.
F. P. Guengerich, Chem. Res. Toxicol. 2001, 14, 611–650.
E. M. Isin, F. P. Guengerich, Biochim. Biophys. Acta - Gen.
Subj. 2007, 1770, 314–329.
This work was supported by the Netherlands Organization for
Scientific Research (NWO, VIDI grant 723.012.105). We thank
Marieke Heijink for her contribution.
[20]
[21]
[22]
X. Zhang, S. Li, Nat. Prod. Rep. 2017, 34, 1061–1089.
F. P. Guengerich, Biol. Chem. 2002, 383, 1553–1564.
J. A. Krauser, F. P. Guengerich, J. Biol. Chem. 2005, 280,
19496–506.
References
[1]
[2]
P. R. Ortiz de Montellano, Cytochrome P450: Structure,
Mechanism, and Biochemistry, Springer Science &
Business Media, 2005.
[23]
[24]
F. P. Guengerich, J. A. Krauser, W. W. Johnson, Am.
Chem. Soc. 2004, 43, 10775–10788.
M. A. Noble, C. S. Miles, S. K. Chapman, D. A. Lysek, A. C.
Mackay, G. A. Reid, R. P. Hanzlik, A. W. Munro, Roles of
Key Active-Site Residues in Flavocytochrome P450 BM3,
1999.
V. L. Arcus, E. J. Prentice, J. K. Hobbs, A. J. Mulholland, M.
W. Van der Kamp, C. R. Pudney, E. J. Parker, L. A.
Schipper, Biochemistry 2016, 55, 1681–1688.
D. N. Kim, K.-H. Cho, W. S. Oh, C. J. Lee, S. K. Lee, J.
Jung, K. T. No, J. Chem. Inf. Model. 2009, 49, 1643–1654.
J. I. Seeman, Chem. Rev. 1983, 83, 83–134.
K.-H. Kim, E. M. Isin, C.-H. Yun, D.-H. Kim, F. P.
Guengerich, FEBS J. 2006, 273, 2223–2231.
L. C. Bell, F. P. Guengerich, J. Biol. Chem. 1997, 272,
29643–29651.
[25]
[3]
[4]
T. W. . Ost, C. S. Miles, J. Murdoch, Y.-F. Cheung, G. A.
Reid, S. K. Chapman, A. W. Munro, FEBS Lett. 2000, 486,
173–177.
[26]
[27]
A. W. Munro, D. G. Leys, K. J. McLean, K. R. Marshall, T.
W. . Ost, S. Daff, C. S. Miles, S. K. Chapman, D. A. Lysek,
C. C. Moser, C. C. Page, P. L. Dutton, Trends Biochem. Sci.
2002, 27, 250–257.
[28]
[29]
K.-H. J. Ling, R. P. Hanzlik, Biochem. Biophys. Res.
Commun. 1989, 160, 844–849.
[5]
A. J. Warman, O. Roitel, R. Neeli, H. M. Girvan, H. E.
Seward, S. A. Murray, K. J. McLean, M. G. Joyce, H.
Toogood, R. A. Holt, D. Leys, N. S. Scrutton, A. W. Munro,
Biochem. Soc. Trans. 2005, 33, 747–753.
[30]
[31]
F. P. Guengerich, Methods Enzymol. 2017, 596, 217–238.
C. F. Butler, C. Peet, A. E. Mason, M. W. Voice, D. Leys, A.
W. Munro, J. Biol. Chem. 2013, 288, 25387–25399.
L. Capoferri, R. Leth, E. ter Haar, A. K. Mohanty, P. D. J.
Grootenhuis, E. Vottero, J. N. M. Commandeur, N. P. E.
Vermeulen, F. S. Jørgensen, L. Olsen, D. P. Geerke,
Proteins Struct. Funct. Bioinforma. 2016, 84, 383–396.
H. Venkataraman, M. C. A. Verkade-Vreeker, L. Capoferri,
D. P. Geerke, N. P. E. Vermeulen, J. N. M. Commandeur,
Bioorg. Med. Chem. 2014, 22, 5613–5620.
[6]
[7]
B. M. A. van Vugt-Lussenburg, E. Stjernschantz, J.
Lastdrager, C. Oostenbrink, N. P. E. Vermeulen, J. N. M.
Commandeur, J. Med. Chem. 2007, 50, 455–461.
C. J. C. Whitehouse, S. G. Bell, H. G. Tufton, R. J. P.
Kenny, L. C. I. Ogilvie, L.-L. Wong, Chem. Commun. 2008,
8, 966–968.
[32]
[33]
[34]
[8]
[9]
H. Arnold, J. C. Lewis, F. H. Arnold, Chimia (Aarau). 2009,
63, 309–312.
P. Rydberg, P. Vasanthanathan, C. Oostenbrink, L. Olsen,
E. Weber, A. Seifert, M. Antonovici, C. Geinitz, J. J. Pleiss,
13
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
L. Higgins, K. R. Korzekwa, S. Rao, M. Shou, J. P. Jones,
Arch. Biochem. Biophys. 2001, 385, 220–230.
D. Y. Curtin, Rec. Chem. Prog 1954, 15, 110–128.
E. Vottero, V. Rea, J. Lastdrager, M. Honing, N. P. E.
Vermeulen, J. N. M. Commandeur, JBIC J. Biol. Inorg.
Chem. 2011, 16, 899–912.
ChemMedChem 2009, 4, 2070–2079.
[60]
[35]
[36]
M. Kazemi, J. Åqvist, Nat. Commun. 2015, 6, 7293.
P. Rydberg, R. Lonsdale, J. N. Harvey, A. J. Mulholland, L.
Olsen, J. Mol. Graph. Model. 2014, 52, 30–35.
P. Rydberg, U. Ryde, L. Olsen, J. Chem. Theory Comput.
2008, 4, 1369–1377.
[61]
[62]
[37]
[38]
[39]
[40]
[41]
[42]
[43]
[44]
L. Olsen, P. Rydberg, T. H. Rod, U. Ryde, J. Med. Chem.
2006, 49, 6489–6499.
[63]
[64]
[65]
[66]
J. Oláh, A. J. Mulholland, J. N. Harvey, Proc. Natl. Acad.
Sci. 2011, 108, 6050–6055.
P. Rydberg, U. Ryde, L. Olsen, J. Phys. Chem. A 2008, 112,
13058–13065.
P. Rydberg, D. E. Gloriam, L. Olsen, Bioinformatics 2010,
26, 2988–2989.
T. S. Dowers, D. A. Rock, D. A. Rock, B. N. S. Perkins, J. P.
Jones, Drug Metab. Dispos. 2004, 32, 328–332.
B. Meunier, S. P. De Visser, S. Shaik, Chem. Rev. 2004,
104, 3947–3980.
R. Leth, B. Ercig, L. Olsen, F. S. Jørgensen, J. Chem. Inf.
Model. 2019, 59, 743–753.
S. B. A. de Beer, L. A. H. van Bergen, K. Keijzer, V. Rea, H.
Venkataraman, C. Fonseca Guerra, F. M. Bickelhaupt, N. P.
E. Vermeulen, J. N. M. Commandeur, D. P. Geerke, Curr.
Drug Metab. 2012, 13, 155–166.
S. Shaik, D. Kumar, S. P. De Visser, A. Altun, W. Thiel,
Chem. Rev. 2005, 105, 2279–2328.
A. Wang, U. Savas, M.-H. Hsu, C. D. Stout, E. F. Johnson,
J. Biol. Chem. 2012, 287, 10834–10843.
[67]
[68]
R. Omura, T., Sato, J. Biol. Chem. 1964, 239, 2379–85.
V. Derdau, J. Atzrodt, J. Zimmermann, C. Kroll, F. Brückner,
Chem. Eur. J. 2009, 15, 10397–10404.
P. A. Williams, J. Cosme, D. M. Vinkovic, A. Ward, H. C.
Angove, P. J. Day, C. Vonrhein, I. J. Tickle, H. Jhoti,
Science 2004, 305, 683–686.
[69]
[70]
A. Sali, Mol. Med. Today 1995, 1, 270–277.
S. Sansen, J. K. Yano, R. L. Reynald, G. A. Schoch, K. J.
Griffin, C. D. Stout, E. F. Johnson, J. Biol. Chem. 2007, 282,
14348–14355.
[45]
[46]
P. A. Williams, J. Cosme, A. Ward, H. C. Angove, D. Matak
Vinković, H. Jhoti, Nature 2003, 424, 464–468.
V. Rea, A. J. Kolkman, E. Vottero, E. J. Stronks, K. A. M.
Ampt, M. Honing, N. P. E. Vermeulen, S. S. Wijmenga, J. N.
M. Commandeur, Biochemistry 2012, 51, 750–760.
S. Poli-Scaife, R. Attias, P. M. Dansette, D. Mansuy,
Biochemistry 1997, 36, 12672–12682.
[71]
[72]
[73]
O. Korb, T. Stützle, T. E. Exner, Swarm Intell. 2007, 1,
115–134.
O. Korb, T. Stützle, T. E. Exner, J. Chem. Inf. Model. 2009,
49, 84–96.
[47]
[48]
“Molecular Operating Environment (MOE) software;
Chemical Computing Group Inc.
K. A. Regal, S. D. Nelson, Arch. Biochem. Biophys. 2000,
384, 47–58.
http://www.chemcomp.com,” 2012.
[49]
[50]
H. Li, T. L. Poulos, Nat. Struct. Biol. 1997, 4, 140–6.
S. Modi, M. J. Sutcliffe, W. U. Primrose, L.-Y. Lian, G. C. K.
Roberts, Nat. Struct. Biol. 1996, 3, 414–417.
I. Schlichting, J. Berendzen, K. Chu, A. M. Stock, S. A.
Maves, D. E. Benson, R. M. Sweet, D. Ringe, G. A. Petsko,
S. G. Sligar, Science 2000, 287, 1615–22.
[74]
[75]
R. A. Luirink, S. J. Dekker, L. Capoferri, L. F. H. Janssen, C.
L. Kuiper, M. E. Ari, N. P. E. Vermeulen, J. C. Vos, J. N. M.
Commandeur, D. P. Geerke, J. Inorg. Biochem. 2018, 184,
115–122.
[51]
M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M.
S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A.
Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery,
J. Comput. Chem. 1993, 14, 1347–1363.
[52]
[53]
[54]
S. Arrhenius, Zeitschrift für Phys. Chemie 1889, 4, 96–116.
S. Arrhenius, Zeitschrift für Phys. Chemie 1889, 4, 226–248.
X. Yuan, Q. Wang, J. H. Horner, X. Sheng, M. Newcomb,
Biochemistry 2009, 48, 9140–9146.
[76]
[77]
J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A.
Case, J. Comput. Chem. 2004, 25, 1157–1174.
E. A. Rifai, M. Van Dijk, N. P. E. Vermeulen, A. Yanuar, D.
P. Geerke, J. Chem. Inf. Model. 2019, 59, 4018–4033.
[55]
[56]
[57]
[58]
[59]
S. Puntarulo, A. I. Cederbaum, Arch. Biochem. Biophys.
1989, 269, 569–575.
J. A. Peterson, R. E. Ebel, D. H. O’keeffe, T. Matsubara, R.
W. Estabrook, J. Biol. Chem. 1976, 251, 4010–4016.
W. Duppel, V. Ullrich, Biochim. Biophys. Acta - Biomembr.
1976, 426, 399–407.
C. S. Yang, F. S. Strickhart, L. P. Kicha, Biochim. Biophys.
Acta - Biomembr. 1977, 465, 362–370.
D. M. Enoch, Zeitschrift für Naturforsch. C 1994, 49, 763–
771.
14
This article is protected by copyright. All rights reserved.

10.1002/cbic.201900751
ChemBioChem
FULL PAPER
Entry for the Table of Contents
FULL PAPER
The current work presents
modified Arrhenius approach to
enable studying the thermodynamic
a
R.A. Luirink, M.C.A Verkade-Vreeker,
dr. J.N.M Commandeur, dr. D.P.
Geerke
determinants
Cytochrome
of
P450
selectivity
in
catalyzed
Page No. – Page No.
substrate conversion.
A modified Arrhenius approach to
thermodynamically study regio-
selectivity in Cytochrome P450
catalyzed substrate conversion
15
This article is protected by copyright. All rights reserved.