Single or Multiple Access Channels to the CYP450s Active Site? An

Letter pubs.acs.org/JPCL

Single or Multiple Access Channels to the CYP450s Active Site? An Answer from Free Energy Simulations of the Human Aromatase Enzyme Alessandra Magistrato,*,†,# Jacopo Sgrignani,‡,# Rolf Krause,§ and Andrea Cavalli*,‡ †

CNR-IOM c/o SISSA, via Bonomea 265, 34136 Trieste, Italy Institute for Research in Biomedicine (IRB), Università della Svizzera Italiana (USI), Via Vincenzo Vela 6, CH-6500 Bellinzona, Switzerland § Institute of Computational Science, Faculty of Informatics, Università della Svizzera Italiana (USI), Via Giuseppe Buffi 13, CH-6900 Lugano, Switzerland ‡

S Supporting Information *

ABSTRACT: Cytochromes P450 (CYP450s), in particular, CYP19A1 and CYP17A1, are key clinical targets of breast and prostate anticancer therapies, critical players in drug metabolism, and their overexpression in tumors is associated with drug resistance. In these enzymes, ligand (substrates, drugs) metabolism occurs in deeply buried active sites accessible only via several grueling channels, whose exact biological role remains unclear. Gaining direct insights on the mechanism by which ligands travel in and out is becoming increasingly important given that channels are involved in the modulation of binding/dissociation kinetics and the specificity of ligands toward a CYP450. This has profound implications for enzymatic efficiency and drug efficacy/ toxicity. Here, by applying free energy methods, for a cumulative simulation time of 20 μs, we provide detailed atomistic characterization and free energy profiles of the entry/exit routes preferentially followed by a substrate (androstenedione) and a last-generation inhibitor (letrozole) to/from the catalytic site of CYP19A1 (the human aromatase (HA) enzyme), a key clinical target against breast cancer, studied here as prototypical CYP450. Despite the remarkably different size/shape/hydrophobicity of the ligands, two channels appear accessible to their entrance, while only one exit route appears to be preferential. Our study shows that the preferential paths may be conserved among different CYP450s. Moreover, our results highlight that, at least in the case of HA, ligand channeling is associated with large enzyme structural rearrangements. A wise choice of the computational method and very long simulations are, thus, required to obtain fully converged quantitative free energy profiles, which might be used to design novel biocatalysts or next-generation cytochrome inhibitors with an in silico tuned Km.

M

connected to the protein surface by several grueling entry/ egress channels.2 A clear understanding of the biological function of these channels as well as of their shape, position, and ligand-channeling energetics has been so far limited, and several controversial aspects are still to be clarified.2−12 Nonetheless, this knowledge is of paramount importance for drug or enzyme design projects aiming at tuning the binding

ammalian cytochromes P450 (CYP450s) are implicated in the metabolism of endogenous hormones and xenobiotic compounds, such as drugs and food additives. Furthermore, some members of this family are targets of breast and prostate cancer therapies, making them particularly relevant to oncology. Mammalian cytochromes are, moreover, involved in the metabolism of about 70% of the clinically used drugs and have been associated with drug resistance in many tumors.1 Many CYP450s are anchored with the N-terminal α-helix to the membrane of the endoplasmic reticulum, and the entire P450 family is characterized by a deeply buried active site © 2017 American Chemical Society

Received: March 23, 2017 Accepted: April 19, 2017 Published: April 19, 2017 2036

DOI: 10.1021/acs.jpclett.7b00697 J. Phys. Chem. Lett. 2017, 8, 2036−2042

Letter

The Journal of Physical Chemistry Letters

channels accessible to ASD and/or LET, followed by US calculations to obtain more accurate estimations of the free energy along the channels. Identif ication of ASD Dissociation Routes. First, to gain a statistically relevant picture of the possible dissociation paths, we performed 100 SMD simulations with a bias on the ASD− Fe centers of geometry distance and a pulling velocity of 0.002 Å/ps. In this way, by applying an isotropic bias, we did not predefine a dissociation direction, but we let the ligand free to choose the more favorable exit route. Then, we analyzed the contacts between ASD and HA during the 100 expulsion trajectories to identify the most probable exit paths. Our analysis revealed that in 72% of the sampled trajectories ASD dissociates trough ChI (tunnel 2a/2b according to the general CYP450 nomenclature of Cojocaru et al.2), while in the remaining 28% of cases ASD dissociates along a second channel ChII (tunnel S of the Cojocaru’s nomenclature (Figure S2 of the Supporting Information)). A full list of residues lining the two channels is reported in Table S1. Consistently with other studies,25 we encountered difficulties in obtaining a convergent free energy profile by using the SMD simulations (see the Supporting Information). Thus, we performed US calculations to obtain quantitatively accurate free energy profiles along the channels.26 For the sake of completeness, together with ChI and ChII, we performed US simulations also for channel ChIII, identified in our previous study on the basis of random expulsion MD (RAMD) simulations,26 that was never observed here in our set of SMD trajectories. The characterization of its free energy profile is, however, useful to assess the reliability of the paths identified on the basis of SMD simulations and compare them with the results of other methodologies such as RAMD.26 Interestingly, the initial part of the free energy profile along ChI is characterized by a steep slope, most likely due to the presence of the Gln225 and Leu228 for which a role as a gatekeeper for the access inside of active site has been suggested (Figure 2).27 The even steeper increase of the free energy profile along ChII is ascribable to a similar gating role involving residues Arg192, Ser224, Leu477, Ser478, and His480. Remarkably, although slightly disfavored energetically, ChII corresponds to the entrance/exit path suggested by Ghosh et al. upon visual inspection of the HA crystal structure.28 Single-

(dissociation) rates of ligands (substrates and drugs) to/from their active site to increase substrate specificity, enzyme efficiency, and drug efficacies.13−17 Atomistic simulations have the potential to clarify such mechanisms at the microscopic level, offering a spatial and temporal resolution inaccessible to experimental methods.11,18,19 Here we address some of these controversial biological issues, focusing on the human aromatase (HA) as a prototypical CYP450.20 The main function of HA is to convert androgens to estrogens,21,22 and its deregulated activity is related to multiple pathological conditions, mainly breast cancer, making it one of the major oncological clinical targets.23,24 An integrated computational protocol, based on (i) multiple long steered molecular dynamics (SMD) simulations, (ii) umbrella sampling (US) calculations with a classical force field, and (iii) quantum−classical (QM/MM) simulations, has been used to gain an atomistic view and a quantitative understanding of the energetics of a substrate (androstenedione, ASD, Figure 1) and a nonsteroidal inhibitor (letrozole,11 LET, Figure 1)

Figure 1. Structures of the substrate ASD and the inhibitor LET of the HA enzyme.

traveling along different HA channels. In particular, we performed multiple SMD simulations to identify putative

Figure 2. (A) Positions of the ASD center of mass along channel I (red van deer Waals (vdw) spheres). Residues Gln225 and Leu228 are shown as vdw spheres and colored by atom type. (B) Positions of the ASD center of mass along channel II (green vdw spheres). Residues Arg192, Asp309, and Ser478 are shown as vdw spheres and colored by atom type.28 2037

DOI: 10.1021/acs.jpclett.7b00697 J. Phys. Chem. Lett. 2017, 8, 2036−2042

Letter

The Journal of Physical Chemistry Letters

Figure 3. (A) Free energy profiles obtained for channels I (red), II (green), and III (blue) computed with US simulations. The error in the free energy profiles was estimanted using Monte Carlo bootstrapping analysis (see the Supporting Information for details). (B) Representative trajectories are drawn by the positions of the ASD center of mass. Red, green, and blue vdw spheres are for channels I, II, and III, respectively.

ChI and ChII are essentially identical (ΔΔGbind# = 1.0 ± 0.6 kcal/mol for ChI and ChII, respectively). In line with these considerations, the fact that also ChII is viable is consistent with experimental28 (ChII is detectable in the X-ray structure; mutations of residues forming this channel influence both Km and Vmax) and computational evidence (ChII is identified by a simple geometrical analysis of the protein cavities performed by us).26 Finally, ChIII, never observed in the SMD simulations, has a dissociation ΔGdiss# (18.9 ± 0.4 kcal/mol, Figure 3) markedly higher than that of ChI and ChII and therefore does not appear to be a viable conduit. RAMD simulations reported so far in the literature were carried out with the aim to observe ligand dissociation in a very short simulation time.8,26,33−35 This can provide a useful qualitative picture of possible CYP450s channels. 26 However this result suggests that RAMD simulations must be performed with milder simulations conditions (i.e., longer simulation time, less strong bias, etc.8) or different computational methods are required to obtain quantitativelly more reliable results. Identif ication of LET Dissociation Routes. Third-generation nonsteroidal HA inhibitors currently in the therapy (anastrozole, LET and vorozole) are characterized by the common presence of a triazole moiety that, at least in the case of anastrozole, has been proven to interact with the heme iron atom.36 Their exact inhibition mechanism is, however, still debated as noncompetitive, or mixed modes have been recently put forward.37 As a representative of this class of inhibitors, we selected LET. This inhibitor has a different chemical structure and shape compared to the HA natural substrates such as ASD (Figure 1). In fact, it is remarkably bulkier (181.5 Å3 for LET and 160 Å3 for ASD) and more lipophilic than ASD (predicted log Po/w: 1.56 vs 3.00, respectively). These properties make it an interesting probe to assess the functional role of various HA channels. Indeed, one of the hypotheses is that channeling preferences depends on ligands size and hydrophobicity.38 Because of the direct coordination bond between LET and the HEME iron atom, to gain insights on the entire dissociation process, the protocol used for ASD was complemented with QM/MM simulations.39,40 Namely, we estimated the contribution to the free energy profile given by the breakage of the coordination bond and to the regeneration of the enzyme

point mutations of the residues lining both channels (Lys119, Gly121, Met374, and Ile395 for ChI and Ile474, Asp476, His475, Ser478, and His480 for ChII)22 have been shown to reduce the Km and Vmax values, strengthening the idea that both paths are involved in ligand in/out channeling.22 The free energy profiles (Figure 3) obtained from US simulations revealed interesting differences: (i) the relative free energy barrier for substrate dissociation (ΔGdiss#) along the three channels is in qualitative agreement with the dissociation frequency along the different paths observed with SMD simulations (ΔGdiss# = 8.5 ± 0.2, 10.6 ± 0.3, and 18.9 ± 0.4 kcal/mol for ChI, ChII, and ChIII, respectively); (ii) the main transition states (TSs) are distant from the active site as they are located near the entrance/exit of the channel, consistently with other studies.29−31 In fact, in ChI, the main TSIASD lies at a distance of ∼27−28 Å from the catalytic site, and it is lined by Ser90, Gly117, and Asp232. In ChII, the main TSIIASD is at a distance of about ∼27 Å and ASD interacts with Val215, Glu483, and Thr484, while in ChIII, the TSIIIASD is at about ∼20 Å, with ASD in the vicinity of Ile255, Leu251, and Leu301. This suggests that a change in solvation free energy, difficult to evaluate computationally, could partially contribute to the main free energy barrier. The frequency of the dissociation events along the distinct channels obtained by SMDs, together with results of mutagenesis experiments,22 suggests that, despite the difference between the main free energy barrier encountered along ChI and ChII (ΔΔGdiss# ≈ 2 kcal/mol) is small, ChI may be the preferential route through which ASD exits from the catalytic site. We are, of course, aware that the absolute values of the binding free energy barriers, as well as the binding free energies, may not be accurate because of finite size/boundary effects associated of the simulation box.32 However, here, we compare the relative free energies along different binding paths traveled by the same ligand within HA. Thus, finite size effects most likely affect the two free energy profiles by the same amount, allowing a ranking the route preferences of a substrate. However, since the absolute ΔGbind# is not meaningful, we report only the difference of the free energy barriers for ASD (or LET) along the two distinct channels. This comparison shows that the free energy costs to enter the active site along 2038

DOI: 10.1021/acs.jpclett.7b00697 J. Phys. Chem. Lett. 2017, 8, 2036−2042

Letter

The Journal of Physical Chemistry Letters

Figure 4. (A) Free energy profiles as a function of the Fe−LET distance (i.e., the Fe−LET distance) computed with US for ChILET (red) and ChIILET (green). Error bars have been computed with Monte Carlo bootstrapping analysis (more details are available in the Supporting Information). The orange line represents the energy profile for the dissociation of the FE−N1@LET bond computed by QM/MM simulations. (B) Representative dissociation routes are shown by the positions of the LET center of mass. Red and green vdw spheres are for ChILET and ChIILET, respectively. (C) Focus on the β8−β9 loop rearrangement observed during US calculations. The loop position in the last frame of the US simulation is depicted in green, while the X-ray structure 3EQM is shown in red.

residue is part of the β8−β9 loop. In ChIILET, the main TS (TSIILET) is located, as for ASD, at the entrance of the protein channel at a distance of 27 Å from the active site near Thr484, Gln218, Glu483, and Asp222. It is worth mentioning that the convergence of the estimated free energy profile along ChILET was reached only after about 100 ns of simulation for each point and that the final profile substantially differs from that obtained after 50 ns (Figure S10A of the Supporting Information), which is the time scale often used in US simulation studies. This behavior, which is different from what we observed for ASD (Figure S10A of the Supporting Information) clearly suggests that large and bulky ligands require significant structural rearrangements when traveling inside of the CYP450s channels. A comparison of the conformation of the protein at the end of the US simulations (after 105 ns US MD calculations at every point along the channel), with the HA crystal structure of HA (pdb code 3EQM)28 did indeed show a marked movement of the loop β8−β9, resulting in a widening of ChILET (Figure 4C). Remarkably, this result highlights that enzyme flexibility makes the channels capable of fully adapting even to large ligands. In this respect, it is clear that even 100 ns long SMD simulations, due to the relatively fast pulling velocity, hamper proper relaxation of the channels. Therefore, the relative ranking of substrate preferences toward the different channels may be qualitatively incorrect. Are the HA Channels Conserved? Despite common structural features shared among all of the over 6800 mammalian CYP450s identified to date,44 their sequence similarity can be rather low (