Mechanism of Cytochrome P450 17A1-Catalyzed Hydroxylase and

Apr 7, 2017 - Department of Drug Design and Pharmacology, Faculty of Health and Medical Sciences, University of Copenhagen, Universitetsparken 2, ...
0 downloads 0 Views 1MB Size
Subscriber access provided by HACETTEPE UNIVERSITESI KUTUPHANESI

Article

Mechanism of Cytochrome P450 17A1 Catalyzed Hydroxylase and Lyase Reactions Silvia Bonomo, Flemming Steen Jørgensen, and Lars Olsen J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.6b00759 • Publication Date (Web): 07 Apr 2017 Downloaded from http://pubs.acs.org on April 9, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Mechanism of Cytochrome P450 17A1-Catalyzed Hydroxylase and Lyase Reactions Silvia Bonomo†, Flemming Steen Jørgensen†, Lars Olsen†* †

Department of Drug Design and Pharmacology, Faculty of Health and Medical Sciences,

University of Copenhagen, Universitetsparken 2, DK-2100, Copenhagen Ø, Denmark.

Keywords: CYP17A1, Hydroxylase, Lyase, Compound I, Peroxoanion, DFT, Docking, Molecular Dynamics

ACS Paragon Plus Environment

1

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 36

ABSTRACT: Cytochrome P450 17A1 (CYP17A1) catalyzes C17-hydroxylation of pregnenolone and progesterone and the subsequent C17-C20 bond cleavage, lyase reaction, to form androgen precursors. Compound I (Cpd I) and peroxoanion (POA) are the heme-reactive species underlying the two reactions. We have characterized the reaction path for both the hydroxylase and lyase reactions using density functional theory (DFT) calculations and the enzyme-substrate interactions by molecular dynamics (MD) simulations. Activation barriers for positions subject to hydroxylase reaction have values close to each other and span from 54 to 60 kJ·mol-1 with a small preference for 17α-hydroxylation in agreement with experimental observations. For the lyase reaction, two different types of mechanisms, concerted and stepwise, were identified with identical activation energies (87 kJ·mol-1). Embedding the DFT optimized transition states (TS) for the two reactions into the active site of CYP17A1 show that the TS for the C17hydroxylation needs to be distorted by 13 kJ·mol-1, whereas the TS for the 17,20-lyase reaction easily can be accommodated in the protein. Finally, differences in the hydrogen bond pattern of the substrates were detected both in the CYP17A1-Cpd I and CYP17A1-POA complexes, with the former found to be pivotal for the hydroxylation site than the latter, suggesting a possible explanation for the slower conversion rate of CYP17A1 for 17α-hydroxyprogesterone over 17α-hydroxypregnenolone. The results support the concept that selectivity of the steroidogenic CYPs is ruled by direct interactions to the enzyme in contrast to the selectivity of drug metabolizing CYPs, where reactivity of the substrates dominate.

ACS Paragon Plus Environment

2

Page 3 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

INTRODUCTION Among the enzymes involved in the steroidogenic processes, cytochrome P450 17A1 (CYP17A1) produces precursors of corticosteroids and sexual hormones by converting pregnenolone and progesterone into their 17α-hydroxyl derivatives and, subsequently, into their 17-oxo products dehydroepiandrostenedione (DHEA) and androstenedione (cf. Scheme 1).1 Conversions take place in the same binding pocket and start with the hydroxylase step, where the radical cation compound I (Cpd I) initiates the C-H homolytic scission on the steroidal scaffold, followed by the C-OH bond formation to release the hydroxylated

products.2-5

While

pregnenolone

is

selectively

converted

to

17α-hydroxypregnenolone, progesterone is hydroxylated in position 16α, 17α, and 21, although the selectivity can be regained by the A105L mutation present in the baboon enzyme, which exclusively yields 17α-hydroxyprogesterone.6 A double binding mode for progesterone in human CYP17A1 has been proposed by combining structural and functional data to explain its lack of selectivity.7 However, the orientation of progesterone, together with the other endogenous substrates, has been elucidated via X-ray crystallography only in the A105L mutant of CYP17A1.7 Deacetylation on position 17 and oxidation of the hydroxyl group constitutes the lyase step necessary for the production of androgens and estrogens (cf. Scheme 1). Thus, CYP17A1 is a validated target in hormone-dependent cancer therapy, especially against prostate cancer.8 Notably, inhibition of the lyase reaction is preferred over the hydroxylase reaction to avoid side effects from overproduction of mineralocorticoids.9 However, no structural differences have yet been identified to account for the hydroxylation versus lyase reactions, preventing structure-based approaches to design inhibitors selective for the lyase step. Elucidation of the reaction mechanism, including characterization of the transition

ACS Paragon Plus Environment

3

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 36

state (TS) for each reaction, is expected to provide information important for understanding and controlling selectivity. The cytochrome reactive species responsible for the lyase reaction has been under debate for the past two decades.6, 10 Experimental data from different research groups suggest that this reaction could proceed via either a radical mechanism catalyzed by Cpd I11 or a nucleophilic attack on C20 of the steroidal scaffold promoted by a peroxoanion (POA) .1214

However, the steps underlying the formation of 17-oxo products from peroxohemiketal

intermediates remain unknown. Furthermore, the 17,20-lyase reaction is faster for 17α-hydroxypregnenolone than for 17α-hydroxyprogesterone,15 suggesting different interactions of the substrates with enzyme pocket and POA.6, 16 Oxidative C-C bond cleavage is a reaction carried out primarily by CYP11A1, CYP19A1 and CYP51 in addition to CYP17A1.1 However, the lyase reaction being the object of our study is more similar to the conversion of nabumetone into its active metabolite 6-methoxy-2-naphtyl acetic acid performed by CYP1A2.17 For both enzymes, the catalytic mechanism is initiated by a Cpd I-mediated hydrogen abstraction from an alkyl carbon (cf. Scheme 2), followed by a C-C bond cleavage of the hydroxylated nabumetone to an aldehyde, which finally is oxidized to its correspondent acid.17 These last two steps are combined in the CYP17A1 lyase reaction. Recently, cryoradiolysis and Raman spectroscopy techniques have been used to characterize two structures involved in the earlier stages of lyase catalysis.14 A directional hydrogen bond between the 17α-hydroxyl group of the substrate and the proximal oxygen (Op) of POA, namely the atom bound to the ferric ion, has been found in both structures. Notably, one of the two complexes is poised for nucleophilic attack on C20 by the POA terminal oxygen (Ot), whereas the second structure is the so-called peroxohemiketal intermediate which decays into products.14 Furthermore, Sen and Hackett studied the

ACS Paragon Plus Environment

4

Page 5 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

deformylation reaction of Mycobacterium tuberculosis CYP51 with computational chemistry methods and identified a stepwise mechanism with the C–C bond cleavage from the deprotonated peroxohemiketal species being the energetically most favored.18 In this paper, we have investigated how CYP17A1 catalyzes the hydroxylase and lyase reactions. Density functional theory (DFT) calculations, docking studies, and molecular dynamics (MD) simulations were applied to characterize the most important species in the hydroxylase and lyase catalytic mechanisms. DFT calculations were used to determine the transition state structures and energies for each possible position of the steroidal scaffold undergoing catalysis by CYP17A1. Docking and MD simulations linked the obtained geometries with the actual conformations in the protein pocket and provided insights into substrate selectivity. The unveiled features are important for understanding differences between the steroidogenic and drug metabolizing CYPs as well as for the development of inhibitors with an increased selectivity for the lyase reaction over the hydroxylase reaction.19

COMPUTATIONAL METHODS DFT Calculations. Activation barriers were calculated at the DFT/B3LYP20-22 and B3LYP-D323 level of theory with a protocol similar to the one adopted by Leth et al.24. Cpd I was modeled as a Fe+5 ion embedded in a porphyrin ring without side chains and with SCH3- and O2- as axial ligands.4, 24 The iron ion of POA was formally a Fe+2 and the formal charge of the system was set to -2 according to the data published by Rydberg et al.25 A simplified substrate formed only by rings C and D was used in the calculations mimicking the full steroidal scaffold. Geometry optimizations, vibrational analyses and solvent calculations were performed in the gas phase using the 6-31G(d) basis set for all atoms26 except iron, for which the double-ζ basis set of Schäfer et al. enhanced with a p

ACS Paragon Plus Environment

5

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 36

function27 was applied. The continuum conductor-like screening model (COSMO)28 with an effective dielectric constant of 4 was applied to calculate the effect of the solvent. Final energies for the hydroxylase mechanism (cf. Table S1 and S2, Supporting Information) were calculated at the B3LYP(-D3)/6-311++G(2d,2p) level – or bs2 - with the double-ζ basis set of Schäfer et al.27 enhanced with s, p, d, and f functions on the iron29 and corrected with the 6-31G(d) level-calculated zero-point vibrational energies (ZPE). The frequency calculations also verified that the structures represented true minima or transition states. Reported final energies for the lyase mechanism were calculated at the B3LYP/def2-TZVPP basis,30 because inconsistency between the wave functions from 631G(d) and bs2 basis set was observed and also def2-TZVPP describes the electronic state of the system more accurately.30 Energies were rectified with ZPE and solvation effects corrected for the charge calculated both at the 6-31G(d) level of theory. Detailed calculations are listed in Table S3 of the Supporting Information. Spin and charge distributions were obtained by Natural Bond Order analysis from the bs2 or def2-TZVPP basis set calculations.31 All the calculations were performed in the doublet and quartet spin states using the software package TURBOMOLE.32 Docking Studies. The crystal structure of CYP17A1 in complex with galeterone (PDB ID 3SWZ) was used as starting point for docking and MD experiments. Firstly, water molecules were deleted from the original PDB structures and then Protein Preparation Wizard protocol implemented in Maestro33 was used to add hydrogen atoms and zero-order bonds and to cap the N- and C-terminal of the protein. The hydrogen-bonding network of the protein was optimized at pH 7 and a restrained minimization allowing heavy atoms to move max. 0.3 Å was carried out using the OPLS2005 force field. Protonation states and tautomers of the ligands were generated at pH 7 ± 2 using Epik.

34

Docking was carried out with heme, Cpd I and POA in the binding pocket and it was

ACS Paragon Plus Environment

6

Page 7 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

performed with Gold Suite version 5.2.235 with the heme-tailored ChemScore scoring function.36 Cpd I and POA were modeled according to the data published by Shahrokh et al.37 A 15 Å box was set around the center of mass of the co-crystallized ligand and 50 poses for each compound were produced, with the exception of 17α-hydroxypregnenolone for which 100 poses were generated. The protein was kept rigid during each docking. A heavy atom RMSD value of 0.46 Å between the docked and the co-crystallized galeterone in the heme-containing CYP17A1 validated our protocol. Molecular Dynamics Simulations. Every representative binding pose obtained from docking was used as starting conformation for MDs. Simulations were performed and analyzed using GROMACS 5.0.2,38 although the antechamber and t-Leap modules of AMBER12 were used to assign the Generalized Amber Force Field (GAFF) parameters to the ligands and the ff99SB force field parameters to the proteins and to neutralize each complex.39, 40 Cpd I and POA were covalently bound to a cysteinate residue via their iron ion. Force field parameters and charges for the two reactive species were retrieved from a recent

publication.37

ACPYPE

converted

topologies

and

coordinates

into

a

GROMACS-readable format.41 Each complex was then solvated with a TIP3P truncated octahedron water box setting a distance between the protein and the box wall to 7 Å and used as starting geometry for the molecular dynamics (MD) simulation under periodic boundary conditions. One sodium ion was added to each of the POA-based systems to neutralize the system, whereas the Cpd I-containing systems were already neutral. The complexes underwent 10000 steps energy minimization of water and counter ions, whereas protein, ligand and compound I or peroxoanion were restrained. Three NVT simulations of 160 ps each at 100K, 200K, and 300 K, respectively, with the ligand positionally restrained with harmonic potentials of force constants of 10000, 5000, and 50 kJ nm-2 mol-1, respectively, allowed the systems to gradually reach 300 K. A 1 ns unrestrained NPT

ACS Paragon Plus Environment

7

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 36

simulation at 300 K and 1.01325 bar closed the equilibration process. Production runs between 30 and 50 ns of extension at 300 K and constant volume followed. The final 20 ns of each trajectory, corresponding to 4000 frames, were analyzed in terms of CYP17A1 reactive species-substrate distances, angles and hydrogen bonds formed. Cutoff values were obtained by averaging the DFT-based TS geometries for the hydroxylase reaction adding a tolerance limit of +2 Å in case of the distance (cf. Figure 1). For the lyase step, the cutoffs were based on the highest values retrieved from the DFT-based TS geometries (cf. Figure 2). All bonds were constrained using the LINCS algorithm42, the time-step for the simulations was set to 2 fs and temperature and pressure regulations were achieved by the Berendsen coupling algorithm.43 A similar protocol was successfully applied by Capoferri et al. in a recent publication,44 albeit we decided to randomly assign molecular velocities five times for each starting complex to increase the efficiency of the conformational sampling.

RESULTS DFT Calculations. Density functional theory (DFT) has previously been employed to investigate the reaction mechanism of several enzymes, including cytochrome P450s.5, 45 In our study, substrates differ only for the substitutions on the A-B ring system, which do not affect the intrinsic reactivity of positions that are relevant for the hydroxylase and lyase reactions. Thus, DFT calculations were performed on steroidal scaffolds containing only ring C and D, as described in the Experimental section. Hydroxylase Reaction. Experimental data for CYP17A1-mediated hydroxylation indicate a preference for the formation of products in position 16α, 17α and 21 (cf. Scheme 1).46, 47 In addition, it has been proposed that position 16β is accessible for

ACS Paragon Plus Environment

8

Page 9 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

H-abstraction while the subsequent OH rebound forms the 16α-OH substrate. 6 The intrinsic reactivity of each of these positions was therefore investigated to compare calculated activation barriers with experimental data. Aliphatic hydroxylation has been extensively studied with DFT calculations.2,

3, 5

The

Cpd I-mediated homolytic C–H bond cleavage has been identified as the rate-determining step of the reaction followed by a rebound step where the C–OH bond is formed.2 In agreement with these data, the hydrogen abstraction barriers calculated for the quartet spin state range from 54 kJ·mol-1 for position 17α, to 65 kJ·mol-1 for position 16β (cf. Table 1). The energies for the doublet state range from 56 kJ·mol-1 to 65 kJ·mol-1 for the position 17α and 21, respectively (cf. Figure 1 and Table S1, Supporting Information). Thus, the activation energies at the different positions and different spin states do not differ significantly, with the lowest activation energy either in the doublet or quartet spin state for the four different positions being 57.3 ± 3.0 kJ·mol-1. The reaction complexes are 5-9 kJ·mol-1 more stable than the separated species indicating a small intermolecular stabilization in the reaction complexes. The B3LYP functional does not consider dispersion effects and this has been related to a decrease in accuracy when predicting barrier heights for hydroxylation reactions.48,

49

Thus, the B3LYP-D3 functional has been used to determine dispersion-corrected activation barriers for the aliphatic hydroxylations of our system. B3LYP-D3 calculated activation energies relative to the reactant complex (45-73 kJ·mol-1, cf. Table 1) are comparable to the activation energies calculated without dispersion (63-70 kJ·mol-1 relative to the reactant complex in the quartet state), although the lower barriers for 16β and 17α hydroxylations relative to 16α and 21 hydroxylations suggest that dispersion is more pronounced in the TS for 16β and 17α hydroxylation. Inspection of the TS geometries confirms this. For example, in the 17α−TS structure, the ligand is almost parallel to the

ACS Paragon Plus Environment

9

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 36

porphyrin ring, maximizing the contribution of non-covalent interactions between the two species,23 which probably causes a decrease in activation energy. Conversely, the ligand-porphyrin angle in the 16α-TS complex is almost orthogonal to the porphyrin ring and, accordingly, no dispersion effect is present, leading to smaller differences in activation energies. Thus, the energies of some of the B3LYP-D3 TS complexes may be artificially too low. The B3LYP-D3 calculated activation energies relative to the separated species are typically 50 kJ·mol-1 lower, but the ordering for the four sites is similar. This suggests that the effect of forming the reaction complex from the separated species is comparable for the four different sites. The energy profiles determined with both B3LYP and B3LYP-D3 functionals show a preference for position 17α as the most favored site for H-abstraction, but the differences between the sites indicate that, from an energetic point of view, hydroxylation is possible in all four positions of the steroidal scaffold. Lyase Reaction. One of the proposed mechanisms for the lyase reaction involves the Cpd I-mediated H-abstraction from the newly added hydroxyl group on position 17α.6,

11

However, scanning the H transfer from the 17-OH group to the O of Cpd I resulted in a linear energy increase, and thus, no TS could be identified (data not shown). Therefore, the steps underlying the POA-based lyase mechanism were investigated with the 17α-hydroxylated derivative of the simplified steroidal scaffold used for the hydroxylase reaction. Calculations were performed in the doublet spin state of POA, since this was found to be energetically more favorable than the quartet spin state, as also reported previously.18, 25 In order to identify the transition state for the lyase mechanism, the Op–Ot, Op–H17, Ot-C20 and C17–C20 distances were explored as potential reaction coordinates. Scanning

ACS Paragon Plus Environment

10

Page 11 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

the distance between POA–Ot and substrate–C20 influences the rest of the complex to a higher extent than changing other bond lengths, highlighting its role in the product formation. The initial structures used to investigate the feasibility of a stepwise mechanism were generated with Ot–C20 = 1.70 Å and C17–C20 = 1.63 Å for complex 3 shown in Figure 2 and 2.10 Å and 1.45 Å, respectively, for complex 4, in accordance with previous findings published for the CYP51-catalyzed deformylation.18 Scanning the Ot-C20 and C17-C20 distances as reaction coordinates led to a flat potential surface for the nucleophilic attack on C20. Thus, no TS structure was identified from 1 to the deprotonated peroxohemiketal complex 3, which subsequently undergoes homolytic C17– C20 bond cleavage to yield 4 (Figure 2 and Table S4) as described for CYP51.18 The TS 4 rearranges into the intermediate structure 5, which is a radical on C17. Finally, 5 decays into the products 6 by a barrierless mechanism involving the charge-driven transfer of H17 from the substrate to POA-Op (Figure 2). We also identified the TS for a concerted pathway (cf. 2, Figure 2) with the same energy as the TS for the stepwise pathway (cf. 4, Figure 2). The two TS complexes differ mainly for the Op–H17 and C17–C20 bond lengths denoted d1 and d3, respectively (cf. Figure 2). The TS 2 is characterized by shorter d1 and d3 distances than TS 4, which allows the simultaneous cleavage of these bonds to release the product 6. Natural population analysis showed a higher spin density on the Fe ion in 2 compared to 4. The spin distribution suggests a non-radical C17–C20 and Op–Ot bond scission from the concerted TS 2 with the simultaneous formation of Op–H17 and C17=O17 bonds as shown in Figure S2, Supporting Information. Based on the solvation-corrected potential energy profiles for the lyase reaction it seems that both the concerted and stepwise mechanisms are equally accessible, although the solvation-free energy values (cf. Table S3, Supporting Information) indicate the stepwise

ACS Paragon Plus Environment

11

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 36

mechanism as the preferred one for the lyase reaction. Moreover, the reaction is driven by a radical C17–C20 bond cleavage for the stepwise mechanism, whereas this is not observed in the concerted mechanism. Docking and Molecular Dynamics. Substrate orientations and specificity have been studied with X-ray crystallography in the past years using the A105L mutant of the human CYP17A1.7 Here, the steric hindrance of Leu105 is reducing the volume of the binding pocket, thus stabilizing ligands, especially progesterone, in only one possible conformation which favor the hydroxylation in position 17.7 To avoid such bias, we have investigated substrate orientations in wild type CYP17A1.50 Binding of Hydroxylase Substrates. Pregnenolone and progesterone are the main substrates for the Cpd I-mediated hydroxylase reaction (cf. Scheme 1).6 Figure 3 shows the results from docking pregnenolone and progesterone into Cpd I-containing CYP17A1. Two orientations of pregnenolone and progesterone were identified: one with the substituent in position 3 of the steroidal scaffold hydrogen bonded to Arg239, herein called the A orientation, and one where this substituent was closer to Asn202, the B orientation. After 10 ns MD simulation, these two binding modes of pregnenolone converged into a B-like orientation, in contrast to progesterone, which maintained two distinct binding modes with conformation A moving closer to Asn202, which also rotated its terminal amide group to better interact with the 3-oxo group of the substrate (cf. Figure 3) as argued by Petrunak et al.7 Therefore, extended conformational sampling for at least 30 ns MD simulation (cf. Experimental section) was undertaken for the three distinct binding modes (pregnenolone and progesterone A and B) to assess interactions with CYP17A1 and differences in accessibility of positions subject to hydroxylase activity. In Figure 4 we have summarized Cpd I-substrate distances and angles from the last 20 ns of each MD simulation, while average values are reported in Table S5 of the Supporting

ACS Paragon Plus Environment

12

Page 13 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Information. H17α is clearly the most accessible position for both pregnenolone and progesterone during the H-abstraction step of the hydroxylase reaction. Notably, the I helix in CYP17A1 forces the Fe=O---H17α angle to increase by about 30° compared to the DFT determined TS structure. The penalty associated with this reorientation amounts to 13 kJ·mol-1 (cf. Figure S3, Supporting Information) estimated by distorting the Fe=O---H17α angle of the DFT-obtained TS structure. Both the H16α and H21 positions are accessible, although to a lesser extent than H17α with the H16α–Cpd I distance and the Fe=O---H21 angle, respectively, being the limiting factors. In the MD, the H16β position is not accessible, primarily due to the H16β–Cpd I distance never getting close enough to the oxoiron group. This is in contrast to the observation by Yoshimoto and Auchus that Habstraction of H16β can occur. 6 Analysis of substrate-protein hydrogen bond interactions over the MD simulations is reported in Figure S4 of Supporting Information. Both pregnenolone and progesterone make water-mediated hydrogen bonds to the protein with their O3 and O20, while only pregnenolone interacts directly with Asn202 in the protein in almost all the simulations analyzed. During the simulations we observed that the amide group of Asn202 rotated to optimize the interaction with progesterone. However, only two of the simulations starting with progesterone in the A conformations have about 50% of the frames hydrogen bonded to either Arg239 or Asn202 (cf. system 1 and 4, Figure S4, Supporting Information). Trajectories from the simulations starting from conformation B of progesterone have less than 10% of frames with progesterone directly interacting with CYP17A1 residues. Binding of Lyase Substrates. Compounds produced in the hydroxylase reaction undergo further modifications in the lyase reaction, which have POA as reactive species (cf. Scheme 1). 17α−Hydroxypregnenolone is converted faster than 17α−hydroxyprogesterone

ACS Paragon Plus Environment

13

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 36

in human CYP17A1 which has been linked to the different hydrogen bonding in the binding pocket.16 Our docking results show only one binding mode for each substrate with 17α−hydroxypregnenolone hydrogen-bonded to Asn202 and with the steroid part forming the same angle with POA as described in literature (cf. Figure S5, Supporting Information).7 Figure 5 shows the accessibility of the hydroxyl group in position 17α and C20, representing the preferred sites for attack in the lyase mechanism, while average geometric values are collected in Table S6, Supporting Information. For both lyase substrates, four of the five simulations showed that the hydroxyl group on position 17α adopts a Fe-Op---H distance, i.e. d1 of Figure 2, within the cutoff value in all the frames analyzed. Overall, the three-dimensional conformation of CYP17A1 is maintained in presence of both lyase substrates. Asn202 rotates its amide side chain in order to make a hydrogen bond between its amino group and 17α-hydroxyprogesterone, as also noticed for the Cpd I-based simulations in presence of progesterone. 17α-Hydroxypregnenolone makes watermediated hydrogen bonds in all the five MD simulations, whereas only half as many hydrogen bonds are present in the MD simulations of 17α-hydroxyprogesterone. A distinct difference

was

observed

for

the

hydrophilic

substrate-protein

interactions:

17α-hydroxypregnenolone interacts with Asn202, in agreement with the experimental data,16 whereas 17α-hydroxyprogesterone establishes hydrophilic contacts in only one MD simulation and for less than 25% of the time analyzed (cf. Figure S6, Supporting Information).

DISCUSSION

ACS Paragon Plus Environment

14

Page 15 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

CYP17A1 catalyzes two subsequent reactions on its two substrates, namely progesterone and pregnenolone (cf. Scheme 1). In the first reaction, a hydroxyl group is added to the substrates primarily in position 17α, but traces of hydroxylation in position 16α and 21 have been detected for progesterone. Subsequently, the 17α−hydroxylated substrates are converted into the 17-oxo products in the lyase reaction, which is mandatory for the production of steroidal hormones (cf. Scheme 1). Furthermore, CYP17A1 shows substrate specificity especially in the lyase reaction, with 17α−hydroxypregnenolone being converted faster than 17α−hydroxyprogesterone by the human enzyme.6, 12 We initially investigated both mechanisms of reactions with DFT calculations to map the reaction path, to determine activation energies and to identify stable TS geometries. Subsequently, docking and MD simulations were used to evaluate, if the DFT obtained geometries could be accommodated in the protein and to unveil important substrate-protein interactions. B3LYP-barrier heights for the H-abstraction are similar for positions 16α, 16β, 17, and 21 (cf. Figure 1), suggesting that they can all be subject to hydroxylation. Moreover, the activation energies are comparable with previously determined activation energies for aliphatic oxidations.4 The non-covalent interaction term parameterized in B3LYP-D3 lowered the energies of the TS complexes with the substrate and the porphyrin moiety adopting a parallel arrangement, whereas energies for TS with more orthogonal orientations were not affected to the same extent (cf. Figure S1, Supporting Information). This preference for structures maximizing attractive van der Waals interactions is well known for small model systems,48 but may be considerably reduced in an enzyme cavity where additional and competing attractive van der Waals interactions are present. Pregnenolone and progesterone initially had the same two binding modes when docked into CYP17A1, with the substituent in position 3 interacting with Arg239 or Asn202

ACS Paragon Plus Environment

15

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 36

(reported in Figure 3 as orientation A and B, respectively). After a short MD simulation, binding mode A of pregnenolone converged into orientation B indicating that this molecule prefers only one binding mode, orientation B, in CYP17A1. Conversely, both binding modes of progesterone were maintained, as also mentioned previously.7 Position 21 of the substrate was the closest to Cpd I in the A orientation. Binding mode B, on the other hand, had all hydroxylase positions close to Cpd I (cf. Figure 3). To investigate which of the reactions are possible in the protein, we examined whether the TS geometry from the DFT calculations could be mimicked in the binding pocket. Thus, both distance and angle criteria should be satisfied. H17α is the most accessible position for Cpd I among the systems analyzed (cf. Figure 4 and Table S4, Supporting Information). However, the presence of the I helix in CYP17A1 forces pregnenolone and progesterone to bend their Fe=O---H17α angle to values close to 160° (cf. Figure S3, Supporting Information). The energy penalty for distorting the DFT-optimized TS geometry to the protein-induced conformation was only 13 kJ·mol-1, and the increased activation barrier of position H17α is still within the acceptable range of values for hydroxylations.4 Direct interactions between substrate and amino acids in the binding cleft may partially compensate for this minor increase in energy. However, only MD simulations with pregnenolone and progesterone in the A conformation were characterized by direct hydrogen bonds with Asn202 or Arg239 (cf. Figure S4), whereas the B orientation of progesterone exclusively made water-mediated interactions with the protein. This may explain the differences in product distribution between pregnenolone and progesterone during the hydroxylase reaction. The MD simulations allowed us to determine the accessibility of the different sites for hydroxylation (cf. Figure 4). For both pregnenolone and progesterone (A as well as B) the H17α position was clearly the most accessible. Position H16β was not accessible,

ACS Paragon Plus Environment

16

Page 17 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

primarily because the distance from H16β to the heme Fe atom was too large. Although experimental data on 16-deuteroprogesterone have shown that H16β also undergoes H-abstraction, the final product detected was 16α-hydroxyprogesterone, thereby confirming the hydroxylation on the β-face of position 16 is unlikely to occur.6,

47

The

remaining two sites, H16α and H21, were considerably less accessible compared to H17α. The actual accessibilities in percent depend on the cut-off values, but the relative ordering for the four sites remains constant as H17α > Η21 ~ H16α > H16β. Thus, a more quantitative interpretation of the propensities is probably not justified. The above ordering is also in agreement with the experimental data showing that 17α represents the preferred site for hydroxylation for both pregnenolone and progesterone. Whereas there is a clear preference for 17α hydroxylation in the simulation starting from the A conformation of progesterone, it seems that both H16α and H21 hydroxylation is more likely to occur when considering the simulation starting from the B conformations (cf. Figure 4). Thus, it seems that the direct hydrogen bond interactions between CYP17A1 and pregnenolone drive the hydroxylation specifically and efficiently towards position H17α, whereas the water-mediated interactions and the double binding mode characteristic for progesterone favors more promiscuous reactions and, thus, different hydroxylated products are formed. The 17α-hydroxylated compounds produced in the initial hydroxylase reaction are further modified in the lyase reaction to release androgen and estrogen precursors (cf. Scheme 1). With this reaction, the steroid is deacetylated on position 17 while its hydroxyl group is oxidized to an oxo moiety using POA as the catalytic species. A Cpd I-mediated mechanism similar to the H-abstraction has recently been proposed.6, 11 However, when we investigated this possibility by DFT calculations, no TS geometries could be identified.

ACS Paragon Plus Environment

17

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 36

Thus, the Cpd I-mediated mechanism of action were not further pursued, which might agree with the possible involvement of POA in the reaction. 12, 13 Based on experimental data, it has been proposed that the lyase reaction starts with a POA-Ot-mediated nucleophilic attack on the steroidal C20 to form a peroxohemiketal intermediate, which then decays to products.14 In this step, a directional hydrogen bond between POA-Op and the hydroxyl group of the steroidal scaffold is required in order to get the optimal POA-substrate orientation needed for the nucleophilic attack to occur (cf. complex 1, Figure 2).14 Moving further along the reaction pathway, a flat potential energy surface was found for the formation of the deprotonated peroxohemiketal complex 3 (cf. Figure 2) as also reported for the CYP51-mediated deformylation mechanism.18 The peroxohemiketal 3 decayed into products by a stepwise mechanism characterized by the homolytic C20-C17 bond cleavage of the TS structure 4. Our DFT calculations identified also a concerted pathway with a TS (cf. 2, Figure 2) with the same energy as the TS (4) for the stepwise pathway Substrate-POA orientations, namely angle ρ and σ reported in Figure 2, did not significantly differ between 2 and 4, whereas the distance between POA-Op and the hydroxyl group in position 17 of the steroid was considerably longer for 4 than for 2. This can be the reason for the slightly higher energy of 2 when solvation effects are not considered in the calculations (cf. Table S3, Supporting Information). Recently, a peroxohemiketal structure has been described experimentally.14 This complex has a Fe-Op-Ot-C20 bond similar to our DFT-identified structure 3 (cf. Figure 2). However, the experimentally determined peroxohemiketal is a neutral species, whereas our structure is negatively charged (cf. Table S4, Supporting Information). Compared to what was observed for the hydroxylase reaction, DFT-obtained TS geometries for the lyase reaction were more accessible in the CYP17A1 active site, suggesting a minor distortion from the protein (cf. Figure 5). 17α-Hydroxypregnenolone

ACS Paragon Plus Environment

18

Page 19 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

forms stable hydrogen bonds with both the protein, especially Asn202, and water molecules present in the binding pocket (cf. Figure S6, Supporting Information). For the 17α-hydroxyprogesterone-CYP17A1 MD simulations, water-mediated contacts were reduced by 50% compared to the corresponding 17α-hydroxypregnenolone-CYP17A1 simulations and direct protein-substrate hydrogen bonds were basically absent (cf. Figure S6, Supporting Information). The reduced amount of hydrophilic interactions could indicate that the affinity for the substrate is lower, and thus, 17α-hydroxyprogesterone is less efficiently converted than 17α-hydroxypregnenolone. This interpretation is supported by

the

different

conversion

rates

for

17α-hydroxypregnenolone

and

17α-

hydroxyprogesterone determined experimentally.15 Whereas the steroidogenic CYPs are selective enzymes catalyzing specific reactions on structurally similar substrates, the drug metabolizing CYPs are promiscuous. All the isoforms of the drug metabolizing CYPs convert structurally diverse compounds, especially the CYP3A4 isoform. For these CYPs, the selectivity of reaction can be explained to a large extent by considering the DFT calculated energy profiles of the various substrates.5,

51, 52

Also, the formation of different metabolites from the same

molecule, a characteristic feature of drug metabolizing CYPs, suggests that the binding mode plays a less pronounced role in the catalytic process. Conversely, to catalyze site-specific reactions, steroidogenic CYPs have evolved to be highly complementary and provide strong binding with their substrates.53 This work has shown that both energetic and geometric features need to be considered with the latter being most important in order to explain the enzymatic reactions.

CONCLUSIONS

ACS Paragon Plus Environment

19

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 36

CYP17A1 hydroxylase and lyase mechanisms of reactions have been elucidated by DFT calculations and MD simulations. The heme cofactor was replaced by Cpd I and POA, representing the hydroxylase and lyase reactive species, respectively, in the docking studies and the MD simulations to better mimic the enzymatic processes. For the hydroxylase reaction DFT calculated activation barriers were similar for the four positions studied and, accordingly, all sites were possible hydroxylation sites from an energetical point of view. The MD simulations revealed that for both pregnenolone and progesterone H17α is clearly the most accessible position, H16α and H21 display intermediate accessibility, and H16β is practically inaccessible due the distance to the oxo group being slightly longer. Furthermore, the MD simulations showed that the I helix in CYP17A1 would distort the DFT-based optimal TS geometry for the hydroxylation on position H17α, whereas H16α and H21 did not suffer such a strain. The presence of two binding modes of progesterone in Cpd I-containing CYP17A1 has also been confirmed. Substrate-protein hydrogen bond interactions were crucial for the choice of the hydroxylation site, especially in the case of progesterone. This molecule seems to be hydroxylated in position 17α when directly interacting with CYP17A1, whereas hydroxylation may take place on position H16α and H21 when this interaction is absent. Thus, geometric features like accessibility and hydrogen-bonding seem to be more important for explaining the selectivity than energy considerations. The lyase reaction can take place via both a stepwise and a concerted mechanism, since our DFT-calculations show that they are equally possible from an activation energy point of view and that both geometries easily can be accommodated in the protein. Hence, a minor stabilizing effect by CYP17A1 was observed for the lyase reaction compared to the hydroxylase reaction. Furthermore, the lower number of hydrogen bond interactions

ACS Paragon Plus Environment

20

Page 21 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

between 17α-hydroxyprogesterone and CYP17A1 compared to the one measured for 17α-hydroxypregnenolone can play a role for the different substrate affinities reported in literature. It may be a general trend that the specificity of drug metabolizing CYPs primarily is determined by energetic requirements whereas the specificity of the steroidogenic CYPs are determined by geometrical features and that further work on both CYP types may provide more evidence for this trend.

ACS Paragon Plus Environment

21

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 36

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]. Phone (+45) 35336305. Author contributions All authors have given approval to the final version of the manuscript. Notes The authors declare no competing financial interests. ACKNOWLEDGMENT The authors are thankful to Ulf Ryde from the Lund University, Lund, Sweden, for the helpful discussion.

ASSOCIATED CONTENT Supporting Information. Detailed DFT calculated and dispersion-corrected reaction path for hydroxylation and lyase, binding modes of lyase substrates in CYP17A1, analysis of molecular dynamics trajectories, and spin and charge distributions for the lyase reaction are available free of charge (file type PDF). This information is available free of charge via the Internet at http://pubs.acs.org

ABBREVIATIONS CYP17A1, cytochrome P450 isoform 17A1; DFT, density functional theory; TS, transition state; Cpd I, compound I; POA, peroxoanion; MD, molecular dynamics; H-bond, hydrogen bond; RMSD, root mean square deviation.

ACS Paragon Plus Environment

22

Page 23 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

TABLE OF CONTENTS Scheme 1. Overview of the processes carried out by CYP17A1. For the catalytic species, Cpd I and POA, the porphyrin ring is simplified by two thick lines. Scheme 2. Comparison between proposed reaction mechanisms for CYP1A2 (A) and CYP17A1 (B). Figure 1. Summary of the DFT-calculated transition states geometries for the H-abstraction from position H16α, H16β, H17α and H21. Figure 2. Proposed concerted mechanism (solid arrows) and stepwise mechanism (dashed arrows) for the lyase reaction. Figure 3. Conformational sampling of pregnenolone (top) and progesterone (bottom) in CYP17A1 from docking and after 10 ns of MD. Figure 4. Accessibilities from the MD simulations for the positions subject of hydroxylase. Figure 5. Accessibilities from the MD simulations for positions undergoing lyase reaction. Table 1. B3LYP energies relative to the sum of Cpd I and the substrate for the hydroxylase reaction with Cpd I in the quartet spin state.

ACS Paragon Plus Environment

23

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 36

Scheme 1. Overview of the processes carried out by CYP17A1. For the catalytic species, Cpd I and POA, the porphyrin ring is simplified by two thick lines. Adapted from Yoshimoto and Auchus.6

ACS Paragon Plus Environment

24

Page 25 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Scheme 2. Comparison between proposed reaction mechanisms for CYP1A2 (A) and CYP17A1 (B).a

a

Oxygen atoms transferred from Cpd I and POA to the substrates are marked with a star or reported in bold, respectively.

ACS Paragon Plus Environment

25

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 36

Figure 1. Summary of the DFT-calculated transition states geometries for the H-abstraction

from

position

H16α,

H16β,

H17α

and

H21.

ZPE-corrected

B3LYP/6-311++G(2d,2p) energies are compared to the sum of Cpd I, in the quartet spin state, and the substrate energies.

ACS Paragon Plus Environment

26

Page 27 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 2. Proposed concerted mechanism (solid arrows) and stepwise mechanism (dashed arrows) for the lyase reaction. Energies, reported in kJ·mol-1, are calculated with the B3LYP/def2-TZVPP basis set and corrected with the ZPE+solvation contributions calculated at the B3LYP/6-31G(d) level of theory. Energies are relative to the sum of POA and steroidal fragment. Structure 3 was based on a constrained optimization from the geometry identified by Sen and Hackett.18

ACS Paragon Plus Environment

27

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 36

Figure 3. Conformational sampling of pregnenolone (top) and progesterone (bottom) in CYP17A1 from docking and after 10 ns of MD. Heme and protein amino acids are reported in grey sticks, while substrate binding mode A and B are shown in purple and light blue sticks, respectively.

ACS Paragon Plus Environment

28

Page 29 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 4. Accessibilities from the MD simulations for the positions subject of hydroxylase. Cutoff values are based on DFT-calculated geometries reported in Figure 1: d ≤ (1.20 + 2) Å, θ ≥ 125°. The numbers below d + θ are the average percentage with standard deviations of frames satisfying the cutoff criteria in all five simulations. Each system (1-5) was generated from the same starting complex by assigning different random initial velocities for the MD simulations.

ACS Paragon Plus Environment

29

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 36

Figure 5. Accessibilities from the MD simulations for positions undergoing lyase reaction. Cutoff values are based on the DFT calculated geometries reported in Figure 2: d1 ≤ (1.62 + 2) Å, σ ≥ 156°, d2≤ (1.47 + 2) Å, ρ ≥ 112°. Each system (1-5) was generated from the same starting complex by assigning different random initial velocities for the MD simulations. See Figure 2 for the definition of distances d1 and d2, and angles σ and ρ.

ACS Paragon Plus Environment

30

Page 31 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Table 1. B3LYP energies relative to the sum of Cpd I and the substrate for the hydroxylase reaction with Cpd I in the quartet spin state.a

B3LYP H-16α α H-16β β H-17α α H-21 B3LYP-D3 H-16α α H-16β β H-17α α H-21

RC

TS

INT

RC rebound

TS rebound

PROD

TS - RC

-7.2 -4.9 -8.7 -9.2

59.3 65.0 54.2 55.5

7.9 10.8 -39.2 1.7

7.1 16.8 -43.0 -2.7

22.1 22.5 -8.3 24.8

-204.4 -198.9 -198.9 -180.6

66.5 69.8 62.9 64.7

-52.3 -49.1 -55.4 -55.0

10.6 6.18 -10.8 17.7

-39.4 -36.6 33.7 -37.4

-37.2 -35.7 -88.3 -57.3

ND ND -82.0 ND

-240.7 -230.2 -266.5 -210.9

62.9 58.5 44.6 72.8

a

Energies (in kJ·mol-1) determined with the bs2 basis set and corrected for zero-point vibrational energy. Abbreviations: Cpd I (compound I), RC (reactant complex), TS (transition state), INT (intermediate complex), PROD (product), ND (not determined).

ACS Paragon Plus Environment

31

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 36

REFERENCES 1. Miller, W. L.; Auchus, R. J., The Molecular Biology, Biochemistry, and Physiology of Human Steroidogenesis and Its Disorders. Endocr. Rev. 2011, 32, 81-151. 2. Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W., Theoretical Perspective on the Structure and Mechanism of Cytochrome P450 Enzymes. Chem. Rev. 2005, 105, 2279-2328. 3. Groves, J. T., High-Valent Iron in Chemical and Biological Oxidations. J. Inorg. Biochem. 2006, 100, 434-447. 4. Olsen, L.; Rydberg, P.; Rod, T. H.; Ryde, U., Prediction of Activation Energies for Hydrogen Abstraction by Cytochrome P450. J. Med. Chem. 2006, 49, 6489-6499. 5. Rydberg, P.; Jørgensen, F. S.; Olsen, L., Use of Density Functional Theory in Drug Metabolism Studies. Exp. Opin. Drug Met. 2014, 10, 215-227. 6. Yoshimoto, F. K.; Auchus, R. J., The Diverse Chemistry of Cytochrome P450 17A1 (P450c17, CYP17A1). J. Steroid. Biochem. Mol. Biol. 2015, 151, 52-65. 7. Petrunak, E. M.; DeVore, N. M.; Porubsky, P. R.; Scott, E. E., Structures of Human Steroidogenic Cytochrome P450 17A1 with Substrates. J. Biol. Chem. 2014, 289, 3295232964. 8. Yap, T. A.; Carden, C. P.; Attard, G.; de Bono, J. S., Targeting CYP17: Established and Novel Aproaches in Postate Cancer. Curr. Opin. Pharm. 2008, 8, 449-457. 9. Pia, A.; Vignani, F.; Attard, G.; Tucci, M.; Bironzo, P.; Scagliotti, G.; Arlt, W.; Terzolo, M.; Berruti, A., Strategies for Managing ACTH Dependent Mineralocorticoid Excess Induced by Abiraterone. Cancer Treat. Rev. 2013, 39, 966-973. 10. Auchus, R. J.; Miller, W. L., Molecular Modeling of Human P450c17 (17AlphaHydroxylase/17,20-Lyase): Insights into Reaction Mechanisms and Effects of Mutations. Mol. Endocrinol. 1999, 13, 1169-1182. 11. Yoshimoto, F. K.; Gonzalez, E.; Auchus, R. J.; Guengerich, F. P., Mechanism of 17α,20-Lyase and New Hydroxylation Reactions of Human Cytochrome P450 17A1. 18OLabeling and Oxygen Surrogate Evidence for a Role of a Perferryl Oxygen. J. Biol. Chem. 2016, 291, 17143-17164. 12. Gregory, M. C.; Denisov, I. G.; Grinkova, Y. V.; Khatri, Y.; Sligar, S. G., Kinetic Solvent Isotope Effect in Human P450 CYP17A1-Mediated Androgen Formation: Evidence for a Reactive Peroxoanion Intermediate. J. Am. Chem. Soc. 2013, 135, 1624516247. 13. Khatri, Y.; Gregory, M. C.; Grinkova, Y. V.; Denisov, I. G.; Sligar, S. G., Active Site Proton Delivery and the Lyase Activity of Human CYP17A1. Biochem. Biophys. Res. Commun. 2014, 443, 179-184.

ACS Paragon Plus Environment

32

Page 33 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

14. Mak, P. J.; Gregory, M. C.; Denisov, I. G.; Sligar, S. G.; Kincaid, J. R., Unveiling the Crucial Intermediates in Androgen Production. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, 15856-15861. 15. Auchus, R. J.; Lee, T. C.; Miller, W. L., Cytochrome b 5 Augments the 17,20Lyase Activity of Human P450c17 without Direct Electron Transfer. J. Biol. Chem. 1998, 273, 3158-3165. 16. Gregory, M.; Mak, P. J.; Sligar, S. G.; Kincaid, J. R., Differential Hydrogen Bonding in Human CYP17 Dictates Hydroxylation versus Lyase Chemistry. Angew. Chem. Int. Ed. 2013, 52, 5342-5345. 17. Varfaj, F.; Zulkifli, S. N.; Park, H. G.; Challinor, V. L.; De Voss, J. J.; Ortiz de Montellano, P. R., Carbon-Carbon Bond Cleavage in Activation of the Prodrug Nabumetone. Drug Metab. Dispos. 2014, 42, 828-838. 18. Sen, K.; Hackett, J. C., Peroxo−Iron Mediated Deformylation in Sterol 14αDemethylase Catalysis. J. Am. Chem. Soc. 2010, 132, 10293-10305. 19. Bonomo, S.; Hansen, C. H.; Petrunak, E. M.; Scott, E. E.; Styrishave, B.; Jorgensen, F. S.; Olsen, L., Promising Tools in Prostate Cancer Research: Selective NonSteroidal Cytochrome P450 17A1 Inhibitors. Sci. Rep. 2016, 6, Article number 29468. 20. Lee, C.; Yang, W.; Parr, R. G., Development of the Colle-Salvetti SorrelationEnergy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785-789. 21. Becke, A. D., A New Mixing of Hartree–Fock and Local Density‐Functional Theories. J. Chem. Phys. 1993, 98, 1372-1377. 22. Becke, A. D., Density‐Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648-5652. 23. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H., A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. 24. Leth, R.; Rydberg, P.; Jørgensen, F. S.; Olsen, L., Density Functional Theory Study on the Formation of Reactive Benzoquinone Imines by Hydrogen Abstraction. J. Chem. Inf. Model. 2015, 55, 660-666. 25. Rydberg, P.; Sigfridsson, E.; Ryde, U., On the Role of the Axial Ligand in Heme Proteins: a Theoretical Study. J. Biol. Inorg. Chem. 2004, 9, 203-223. 26. Warren J. Henre, L. R., P. von R. Schleyer, John Pople, Ab Initio Molecular Orbital Theory. New York, 1986. 27. Schäfer, A.; Horn, H.; Ahlrichs, R., Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571-2577. 28. Klamt, A.; Schuurmann, G., COSMO: a New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Snergy and Its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799-805.

ACS Paragon Plus Environment

33

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 36

29. Rulíšek, L.; Jensen, K. P.; Lundgren, K.; Ryde, U., The Reaction Mechanism of Iron and Manganese Superoxide Dismutases Studied by Theoretical Calculations. J. Comput. Chem. 2006, 27, 1398-1414. 30. Weigend, F.; Ahlrichs, R., Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. PCCP 2005, 7, 3297-3305. 31. Reed, A. E.; Weinstock, R. B.; Weinhold, F., Natural Population Analysis. J. Chem. Phys. 1985, 83, 735-746. 32. TURBOMOLE, version 6.3; TURBOMOLE GmbH: University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 2011. 33.

Maestro, version 9.8; Schrödinger, LLC: New York, NY, 2014.

34.

Epik, version 2.8; Schrödinger, LLC: New York, NY, 2014.

35. Jones, G.; Willett, P.; Glen, R. C.; Leach, A. R.; Taylor, R., Development and Validation of a Genetic Algorithm for Flexible Docking. J. Mol. Biol. 1997, 267, 727-748. 36. Kirton, S. B.; Murray, C. W.; Verdonk, M. L.; Taylor, R. D., Prediction of binding modes for ligands in the cytochromes P450 and other heme-containing proteins. Proteins: Struct. Funct. Bioinform. 2005, 58, 836-844. 37. Shahrokh, K.; Orendt, A.; Yost, G. S.; Cheatham, T. E., Quantum Mechanically Derived AMBER-compatible Heme Parameters for Various States of the Cytochrome P450 Catalytic Cycle. J. Comput. Chem. 2012, 33, 119-133. 38. Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R., GROMACS: A MessagePassing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43-56. 39. Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C., Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct. Funct. Bioinform. 2006, 65, 712-725. 40. Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Goetz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M. J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; SalomonFerrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12, University of California, San Francisco: 2012. 41. Sousa da Silva, A.; Vranken, W., ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 367. 42. Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M., LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463-1472.

ACS Paragon Plus Environment

34

Page 35 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

43. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R., Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684-3690. 44. Capoferri, L.; Verkade-Vreeker, M. C. A.; Buitenhuis, D.; Commandeur, J. N. M.; Pastor, M.; Vermeulen, N. P. E.; Geerke, D. P., Linear Interaction Energy Based Prediction of Cytochrome P450 1A2 Binding Affinities with Reliability Estimation. PLoS ONE 2015, 10, e0142232. 45. Siegbahn, P. E. M.; Borowski, T., Modeling Enzymatic Reactions Involving Transition Metals. Acc. Chem. Res. 2006, 39, 729-738. 46. Swart, P.; Swart, A. C.; Waterman, M. R.; Estabrook, R. W.; Mason, J. I., Progesterone 16 Alpha-Hydroxylase Activity Is Catalyzed by Human Cytochrome P450 17 Alpha-Hydroxylase. J. Clin. Endocrinol. Metab. 1993, 77, 98-102. 47. Yoshimoto, F. K.; Zhou, Y.; Peng, H.-M.; Stidd, D.; Yoshimoto, J. A.; Sharma, K. K.; Matthew, S.; Auchus, R. J., Minor Activities and Transition State Properties of the Human Steroid Hydroxylases Cytochromes P450c17 and P450c21, from Reactions Observed with Deuterium-Labeled Substrates. Biochemistry 2012, 51, 7064-7077. 48. Lonsdale, R.; Harvey, J. N.; Mulholland, A. J., Inclusion of Dispersion Effects Significantly Improves Accuracy of Calculated Reaction Barriers for Cytochrome P450 Catalyzed Reactions. J. Phys. Chem. Lett. 2010, 1, 3232-3237. 49. Rydberg, P.; Lonsdale, R.; Harvey, J. N.; Mulholland, A. J.; Olsen, L., Trends in Predicted Chemoselectivity of Cytochrome P450 Oxidation: B3LYP Barrier Heights for Epoxidation and Hydroxylation Reactions. J. Mol. Graph. Model. 2014, 52, 30-35. 50. Lonsdale, R.; Harvey, J. N.; Mulholland, A. J., A Practical Guide to Modelling Enzyme-Catalysed Reactions. Chem. Soc. Rev. 2012, 41, 3025-3038. 51. Rydberg, P.; Gloriam, D. E.; Zaretzki, J.; Breneman, C.; Olsen, L., SMARTCyp: A 2D Method for Prediction of Cytochrome P450-Mediated Drug Metabolism. ACS Med. Chem. Lett. 2010, 1, 96-100. 52. Olsen, L.; Oostenbrink, C.; Jorgensen, F. S., Prediction of Cytochrome P450 Mediated Metabolism. Adv. Drug Del. Rev. 2015, 86, 61-71. 53. Poulos, L. T.; Johnson, F. E. Structures of Cytochrome P450 Enzymes. In Cytochrome P450: Structure, Mechanism, and Biochemistry, Ortiz de Montellano, R. P., Ed.; Springer International Publishing: Cham, 2015; Chapter 1, pp 18-21.

ACS Paragon Plus Environment

35

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 36

TOC GRAPHICS

ACS Paragon Plus Environment

36