Drug Metabolism and Disposition Fast Forward
First published on January 24, 2007; DOI: 10.1124/dmd.106.014019
0090-9556/07/3504-689-696$20.00
DMD 35:689-696, 2007
Possible Pathway(s) of Metyrapone Egress from the Active Site of Cytochrome P450 3A4: A Molecular Dynamics Simulation
Weihua Li,
Hong Liu,
Xiaomin Luo,
Weiliang Zhu,
Yun Tang,
James R. Halpert, and
Hualiang Jiang
School of Pharmacy, East China University of Science and Technology, Shanghai, China (W.L., Y.T., H.J.); Center for Drug Discovery and Design, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Graduate School of the Chinese Academy of Sciences, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China (W.L., H.L., X.L., W.Z., H.J.); and Department of Pharmacology and Toxicology, University of Texas Medical Branch, Galveston, Texas (J.R.H.)
(Received November 22, 2006;
accepted January 19, 2007)
 |
Abstract
|
|---|
To identify a possible pathway(s) for metyrapone egress from the active site of P450 3A4, a 5-ns conventional molecular dynamics simulation followed by steered molecular dynamics simulations was performed on the complex with metyrapone. The steered molecular dynamics simulations showed that metyrapone egress via channel 1, threading through the B-C loop, only required a relatively small rupture force and small displacement of residues, whereas egress via the third channel, between helix I and helices F' and G', required a relatively large force and perturbation of helices I, B', and C. The conventional dynamics simulation indicated that channel 2, located between the ß1 sheet, B-B' loop, and F'-G' region, is closed because of the movement of residues in the mouth of this channel. The findings suggest that channel 1 can be used for metyrapone egress, whereas both channel 2 and channel 3 have a low probability of serving as an exit channel for metyrapone. In addition, residues F108 and I120 appear to act as two gatekeepers to prevent the inhibitor from leaving the active site. These results are in agreement with previous site-directed mutagenesis experiments.
Cytochrome P450 (P450) 3A4 plays an important role in the detoxification and metabolic activation of a wide variety of xenobiotic and endogenous compounds, and is responsible for metabolizing 50 to 60% of clinically used drugs, including antibiotics, anticoagulants, anticancer drugs, anxiolytics, and immunosuppressants (Guengerich, 1999
; Tang and Stearns, 2001). P450 3A4 often exhibits positive homotropic cooperativity of substrate oxidation as well as heterotropic activation by effectors such as testosterone and 
-naphthoflavone (He et al., 2003). The cooperativity in the reaction mechanism of P450 3A4 has led to speculation that the enzyme either has an adaptive active site to accommodate multiple compounds simultaneously and/or exists as kinetically distinct conformers (Koley et al., 1997; Harlow and Halpert, 1998; Davydov et al., 2003; He et al., 2003).
Recently the crystal structures of P450 3A4 with and without a bound ligand have been determined (Williams et al., 2004; Yano et al., 2004; Ekroos and Sjogren, 2006). The P450 3A4 structure exhibits an overall fold similar to that of other P450s reported previously, and adopts a closed conformation with the active site deeply buried in the center of the fold. Surprisingly, the binding of the inhibitor metyrapone does not induce a significant change in the protein conformation compared with the ligand-free structure (Williams et al., 2004; Scott and Halpert, 2005). Thus, an intriguing problem arises: how does the protein open up to allow metyrapone to enter or exit the active site? Two obvious openings are found on the protein surface in the crystal structure of P450 3A4 bound with metyrapone (Williams et al., 2004). These two openings are separated by the B-B' loop. One opening exists between the ß1 sheet, B-B' loop, and the end of helix F', whereas the other penetrates through the B-C loop. The first opening is similar to that found in some bacterial P450s (Poulos et al., 1987; Cupp-Vickery and Poulos, 1995; Li and Poulos, 1997) and was also identified by molecular dynamics simulations to be the most probable ligand access/exit channel entrance (Ludemann et al., 2000a,b; Winn et al., 2002). The second opening is similar to the open channel described in the crystal structure of P450 51 (Podust et al., 2001) and P450 152A1 (Wade et al., 2004; Cojocaru et al., 2007), and has recently been identified using steered molecular dynamics (SMD) simulation to be a possible channel for substrate egress in P450 2B1 by our group (Li et al., 2005). In addition, the crystal structure of P450 3A4 complexed with progesterone (Williams et al., 2004) shows that the steroid does not bind in the active site above the heme, as expected, but binds at a peripheral site very close to the Phe-cluster formed by F108, F213, F215, F219, F220, F241, and F304. This peripheral binding site was suggested to be involved in initial substrate recognition. It has been speculated that appropriate movement of residues around the Phe-cluster could result in a channel for a ligand to move from the peripheral binding site to the active site (Williams et al., 2004).
View larger version (39K):
[in this window]
[in a new window]
|
FIG. 1. Ribbon schematic representations of a three-dimensional model of P450 3A4 and the chosen channels examined by CMD and SMD. The heme is represented by a stick; metyrapone is represented by a Corey-Pauling-Koltun model. The two-dimensional structure of metyrapone is shown on the right. Major helices and sheets are labeled. During the SMD simulation, the bound metyrapone is pulled away from the active site using a harmonic potential symbolized by an artificial spring that is connected to the center of mass of metyrapone.
|
|
The finding from the crystal structure of P450 3A4 with metyrapone bound has indicated that three putative routes could serve as the exit channel for metyrapone. An understanding of the active site access/egress channel is helpful to explain the broad substrate specificity and regio- and stereospecificities of P450 3A4. However, the most appropriate pathway remains enigmatic. To determine which channel(s) could be used for metyrapone to exit from the active site of P450 3A4 and identify the role of certain residues during the ligand egress, a conventional molecular dynamics (CMD) simulation combined with steered molecular dynamics simulation was performed on the metyrapone complex. Three channels were considered, as shown in Fig. 1. Channel 1 threads through the B-C loop. Channel 2 is located between the ß1 sheet, the B-B' loop, and the end of helix F'. Channel 3 runs through the active site roof formed by the Phe-cluster. A 5-ns conventional molecular dynamics simulation was first performed on P450 3A4 complexed with metyrapone. After that, SMD simulations were conducted to probe the possible exit channel(s). Finally, the important residues associated with metyrapone egress were characterized.
|
Materials and Methods
|
|---|
The starting structure for simulations was taken from the X-ray structure of P450 3A4 complexed with metyrapone (PDB entry 1W0G
[PDB]
) (Williams et al., 2004). In the crystal structure of the P450 3A4 complex, the coordinates for the N-terminal residues 1–24 were missing. Because the N-terminus does not affect ligand binding/unbinding, the N-terminus was not modeled in the present simulations. The amino acid numbering is the same as in the crystal structure (PDB code 1W0G
[PDB]
). The missing coordinates of residues from 262 to 268 were reconstructed on the basis of the coordinates of the corresponding positions in the crystal structure of ligand-free P450 3A4 (PDB entry 1TQN
[PDB]
) (Yano et al., 2004) using SYBYL version 6.7 (Tripos Inc., St. Louis, MO, 1999). The missing coordinates of the other segment (residues 281–288) were reconstructed using the database searching of the Biopolymer module in SYBYL version 6.7. The modified structure was subjected to energy minimization using the steepest descent method up to the gradient tolerance of 0.05 kcal/(mol · Å) to release the steric clashes.
All molecular dynamics simulations were performed with the GROMACS version 3.2 (Lindahl et al., 2001) software package with the GROMOS96 force field (van Gunsteren et al., 1996). The protein was solvated in a rectangular periodic box with water. The simple point charge model (Berendsen et al., 1981) was used to describe the water molecules. Because the system has a total charge of +3, chlorine counterions were added to neutralize the box charge. The resulting system was formed by the protein, the ligand, 29 crystallization water molecules, 11,305 water molecules, and 3 chlorine ions, for a total of 
39,000 atoms. The solvent was relaxed by energy minimization followed by 50 ps of MD at 300 K while restraining proteins' atomic positions with a harmonic potential. The system was then energy-minimized without restraint with the steepest descent method for 100 steps. Afterward, molecular dynamics equilibration was conducted at 300 K, with decreasing harmonic position restraints on the protein atoms over a 15-ps time interval, followed by 5 ns of MD equilibration without restraints. During the MD simulation, the LINCS algorithm (Hess et al., 1997) was used to constrain all bonds. A dielectric constant of 1 and a time step of 2 fs were used. The electrostatic interactions were calculated using the Particle-Mesh Ewald method (Darden et al., 1993; Essmann et al., 1995) with a 1.0-nm cutoff. The temperature was kept constant by coupling solute and solvent separately to a thermal bath at 300 K with a coupling constant 
T = 0.1 ps. Pressure was kept constant by coupling to a pressure bath at 1.0 bar, using a coupling constant P = 0.5 ps.
SMD is an extended molecular dynamics simulation mimicking the principle of atom force microscopy (Binning et al., 1986). This method exerts an external force on the ligand to accelerate its binding or unbinding. SMD has been successfully used to identify the ligand pathways and to explore the elastic properties of several proteins (Grubmuller et al., 1996; Isralewitz et al., 2001; Xu et al., 2003; Li et al., 2005). In SMD simulation, time-dependent forces are applied to a ligand to accelerate the binding or unbinding process. From the accelerated dissociation process of the ligand, a SMD simulation can reveal information about the protein's flexibility and its response to the dissociation of ligand. Analyses of interactions between the ligand and the protein and the relationship between the applied forces and the ligand pathway can yield important information about the structure-function relationships of the protein-ligand complex, the binding and unbinding pathway(s), and possible mechanisms of ligand recognition.
In the subsequent SMD simulations, metyrapone was pulled out from the binding pocket through the putative exit channels (Fig. 1) by an external force. The directions of metyrapone along with the channels were determined using the criterion that the inhibitor passes through the putative entrances with the smallest collision with protein residues. The center of mass of metyrapone was harmonically restrained to a point moving with a constant velocity along the desired directions. The force field parameters for metyrapone were taken from the GROMOS96 library. Atomic charges of metyrapone were calculated using the ChelpG method (Breneman and Wiberg, 1990), an electrostatic potential fitting approach, at the HF/6-31G* level. This ab initio calculation was carried out using the Gaussian98 program (Frisch et al., 1998). During the SMD simulations, the pulling velocity was set to 0.005 nm · ps–1 and the spring force constant was assigned as 500 kJ · mol–1 · nm–2. The pulling forces and the interactions between the protein residues and metyrapone were monitored along the SMD trajectory. During both CMD and SMD simulations, the interaction between the pyridine nitrogen and iron atom was treated as a noncovalent interaction. Thus, only van der Waals and electrostatic interactions were considered between the pyridine nitrogen and the heme group. At each time point the presence of direct hydrogen bonds, water bridges, and hydrophobic interactions between metyrapone and 3A4 was analyzed using the LIGPLOT (Wallace et al., 1995) program.
|
Results
|
|---|
Equilibration of Metyrapone-P450 3A4 Complex. To investigate residue movement within three putative egress channels and obtain an equilibrated starting structure for subsequent SMD simulation, a 5-ns equilibration molecular dynamics simulation was first performed on P450 3A4 complexed with metyrapone. System total energy and heavy atom root-mean-square deviation (RMSD) from the starting structure are two important criteria for the convergence of the free MD simulation. Figure 2 shows the total energy and RMSD of C from the crystal structure during the equilibration. The system total energy sharply decreases during the first 500 ps and reaches a plateau after 1800 ps. The protein atoms do not significantly deviate from the crystal structure and the C RMSD is stabilized to an average value of 0.225 nm. Although the RMSD appears to rise slightly even after 3000 ps, the fluctuation amplitudes obtained are relatively small.
View larger version (19K):
[in this window]
[in a new window]
|
FIG. 2. Total energy and RMSD with respect to simulation time for 5-ns conventional molecular dynamics simulation on the P450 3A4 model.
|
|
To provide a detailed description of the mobility of the protein residues, the backbone RMSD per residue is shown in Fig. 3. In the whole simulation, residues with larger RMSD belong to the surface regions of protein, which are exposed to solvent. For example, residues from 260 to 268 and from 270 to 280 are located in the G-H loop and H-I loop, respectively. Neither loop is clearly identified in the crystal structure of the P450 3A4-metyrapone complex. Figure 3 also shows the root-mean-square fluctuation of the backbone as a function of residue number for the conventional molecular dynamics simulation. The results indicate that large fluctuations of residues occur in loops connecting secondary structure elements. Such residues are located in the exterior regions of protein, and correspond to the protein regions with large RMSD.
View larger version (24K):
[in this window]
[in a new window]
|
FIG. 3. Backbone RMSD (top) and root-mean-square fluctuation (bottom) with respect to residue during CMD simulation.
|
|
Residue Movement in Three Channels during CMD. In the crystal structure of P450 3A4 with metyrapone bound, two openings on the protein surface are found and separated by the B-C loop. One is located between the ß1 sheet, the end of helix F, and the B-C loop, and the other threads through the B-C loop. In addition, a channel was proposed to stretch from the initial substrate recognition site to the active site above the heme. To examine the change of each channel during CMD, the movement of residues at the entrance of each channel to the active site pocket was analyzed. The distances of center of mass between six pairs of residues residing on opposite sides of each channel were monitored during the 5-ns CMD simulation. These six pairs of residues are between P107 and I120 and between F108 and S119 of channel 1, between N79 and P107 and between N79 and T224 of channel 2, as well as between F213 and F241 and between F213 and F304 of channel 3. Figure 4 displays the fluctuation of the distances with respect to simulation time. For channel 1 and channel 3, the distance essentially retains the value in the crystal structure after the system reached equilibrium, whereas the distance in channel 2 changes significantly from an initial 0.8 nm to 0.6 nm. This significant decrease in the distance results in closing the open channel 2. The channel portions of the superimposed structures for the crystal structure and the snapshot at 4000 ps are depicted in Fig. 5A. The backbone of the ß1 sheet region displays the largest difference between the two structures. It undergoes a relatively large shift between the B-B' loop and helices F' to G'. Differences in the locations of helices F through G and the B-B' loop are also evident. Besides the shift of backbone, the orientation of some residues in channel 2 underwent rearrangement during the CMD simulation. The most significant changes in the orientations are observed for R106 in the B-B' loop, D76, N78, and N79 in the ß1 sheet, and F215 and F220 in helices F' to G', as shown in Fig. 5B. Side chains of these residues occupy the space at the entrance of channel 2, resulting in a complete block of this pathway.
View larger version (25K):
[in this window]
[in a new window]
|
FIG. 4. Distances between P107 and I120 (solid line) and between F108 and S119 (dotted line) of channel 1 (A), between N79 and P107 (solid line) and between N79 and T224 (dotted line) of channel 2 (B), and between F213 and F241 (solid line) and between F213 and F304 (dotted line) of channel 3 (C) with respect to time during CMD simulation.
|
|
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 5. Superposition of channel 2 portions between the crystal structure (green) and the snapshot at 4 ns of CMD simulation (yellow). A, secondary structure elements that differ significantly are labeled. B, residues with significant changes in orientation are labeled.
|
|
In contrast to channel 2, the placement of residues in channel 1 and channel 3 exhibits trivial changes, as shown in Fig. 5A. Because channel 2 is completely closed to metyrapone, the subsequent SMD simulations were only performed along channel 1 and channel 3.
Egress of Metyrapone from P450 3A4 via Channel 1. To compare the possibility of metyrapone egress from the active site along channel 1 and channel 3 and to gain insight into the mechanisms of metyrapone egress, we investigated these two channels by use of SMD simulations. In SMD simulations of protein-ligand unbinding, a force is applied to pull a ligand out of a binding site and a rupture force is computed. SMD provides a means to accelerate the unbinding process through the application of external forces that lower the energy barrier and drive the ligand to exit the active site in a nanosecond time scale. The forces of metyrapone with P450 3A4 were monitored throughout the SMD simulations and the interactions between the protein residues and metyrapone were analyzed along the SMD trajectory.
The force profile of extracting metyrapone out from the binding pocket along channel 1 is shown in Fig. 6A. Before the SMD simulation, metyrapone is tightly bound. Analyses of the force profile and SMD trajectory suggest that the unbinding process of metyrapone can be divided into three phases. In phase 1 (from 0 to 250 ps), metyrapone dissociates from the binding pocket and moves toward the entrance of the pocket. In phase 2 (from 250 to 400 ps), metyrapone passes through the entrance of the pocket. After crossing over the entrance of the pocket, the egress process enters the third phase and metyrapone enters the solvent.
View larger version (14K):
[in this window]
[in a new window]
|
FIG. 6. Force profiles for SMD simulations along two different putative channels. A, channel 1; B, channel 3.
|
|
During the first 110 ps of egress along channel 1, a steady increase of applied force was observed. At the beginning of the SMD simulation, the inhibitor-protein complex was stabilized by hydrophobic interactions between the inhibitor and residues R105, S119, R212, A305, and T309. Several water bridges formed by the inhibitor with residues R105, E374, and the A ring propionate of the heme further stabilized the complex, as shown in Fig. 7A. After 110 ps, the water bridges made by the inhibitor to R105 and the A ring propionate of the heme were broken. This event was manifested by a slight decrease of the applied force.
View larger version (24K):
[in this window]
[in a new window]
|
FIG. 7. Snapshots of the relative positions of metyrapone and the protein throughout the simulation along channel 1 (A, 0 ps; B, 240 ps; and C, 440 ps). The residues within 5 Å from metyrapone at the individual times are shown as green sticks. The heme is shown as red sticks. The hydrogen bonds are shown in red dotted lines. Metyrapone atoms are represented as yellow sticks.
|
|
Between 110 and 200 ps, the inhibitor had to overcome the hydrophobic interactions with residues obstructing the exit and water bridge interactions with residues in the binding pocket. At 130 ps, the inhibitor slightly deviated from its equilibrated position and formed hydrophobic interactions with R105, R212, F304, A305, and T309, and two water bridges to R106 and E374. These interactions led to a sharp increase of the force. At 170 ps, the inhibitor mainly formed hydrophobic interactions with R105, S119, F304, and A305. These hydrophobic interactions were so strong that a force up to 550 pN was necessary to break them. The breaking of these interactions occurred between 170 ps and 200 ps. The force needed for unbinding remained nearly constant during this period. After 200 ps, there was a significant decrease of the applied force.
After 200 ps, the unbinding process entered phase 2. At this time, the pyridine ring A of the inhibitor first reached the entrance of the binding pocket and interacted with residues at the entrance of the binding pocket. At 240 ps, the pyridine ring A of the inhibitor formed hydrophobic interactions with F108, R106, P107, and I120, whereas the pyridine ring B of the inhibitor still stayed in the binding pocket and formed hydrophobic interactions with S119 and R105. The oxygen atom of the inhibitor produced water bridges to the oxygen atom of S119 and to a nitrogen atom of R105, as shown in Fig. 7B. At 270 ps, the water bridges were still present, but the inhibitor formed hydrophobic interactions with N104, E122, P107, and I120. This event led to a decrease of the applied force. The following increase in the force, between 270 ps and 310 ps, was due to other hydrophobic interactions between the inhibitor and G109 and V111. At 325 ps, the inhibitor only formed hydrophobic interactions with P107 and K115 and water bridges to N104, R105, and I120. This corresponded to a local minimum in the force profile. The increase of the applied force between 325 and 410 ps was mainly caused by the hydrophobic interactions between the pyridine ring B of the inhibitor and the residues at the entrance of binding pocket. An analysis of snapshots of the simulation shows that after metyrapone left the active site, several solvent water molecules entered the active site to occupy the space left by the inhibitor.
After 410 ps, the inhibitor had passed through the entrance of the channel and moved into the solvent. Analysis of the interactions between the inhibitor and the residues revealed that the applied force between 410 ps and 480 ps was due to the interactions between the inhibitor and K115 and E122, as shown in Fig. 7C. After 480 ps of simulation, the inhibitor was completely pulled out of the protein.
Relatively small and localized displacement and rotations of protein side chains were sufficient to allow the inhibitor to leave the binding pocket. The RMSD per residue is shown in Fig. 8. Apart from the G-H loop and H-I loop, which were not clearly identified in the crystal structure, the large displacement mainly occurred in the region of the B-C loop. However, the displacement is not more than 0.30 nm. The secondary structures of helices B', I, and C were maintained well, as shown in Fig. 9A.
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 8. A comparison of backbone RMSD with respect to residue during SMD simulations along channel 1 (bold line) and channel 3 (dotted line). The shadowed region from left to right represents the B-C loop region, helices F'-G' region, and the helix I region.
|
|
View larger version (65K):
[in this window]
[in a new window]
|
FIG. 9. Secondary structure elements variation with respect to simulation time during SMD along channel 1 (A) and channel 3 (B). Helices I and B' are labeled.
|
|
Egress of Metyrapone via Channel 3. The Phe-cluster formed by F108, F213, F215, F219, F220, F241, and F304 comprises the roof of the inhibitor binding pocket. To exit through channel 3, the inhibitor must pass the Phe-cluster region. Figure 6B shows the force profile of extracting the inhibitor from the binding pocket along the proposed channel. In contrast to unbinding via channel 1, the rupture force for the inhibitor leaving via channel 3 is relatively high, up to 680 pN. The highest force peak is produced at 220 ps. At this time, the inhibitor formed hydrophobic interactions with F108, F220, F213, F304, R212, A305, and T309 and a water bridge to F213. Several of these residues belong to the Phe-cluster, which demonstrated that the inhibitor had moved toward the entrance of the binding pocket. From 200 ps to 450 ps, the inhibitor interacted with residues at the entrance of the binding pocket and had to overcome the strong hydrophobic interactions with the residues of the Phe-cluster. At 320 ps, the inhibitor formed hydrophobic interactions with F108, F213, F220, V240, F241, and F304, a hydrogen bond with the nitrogen atom of F213, and a water bridge to R212. At 405 ps, the B ring of the inhibitor broke through the bottleneck, but the A ring still formed hydrophobic interactions with F213, F219, F220, V240, and P242. After 450 ps, the inhibitor moved into the solvent and only formed several water bridges to residues in the F'-G' region.
The RMSD per residue along channel 3 is shown in Fig. 8. Besides the G-H loop and H-I loop, which are exposed to the solvent, the B-C loop and F'-G' region display large displacement, up to 0.34 nm and 0.41 nm, respectively. In addition, the most conserved helix I has a deviation of 0.22 nm. The secondary structures of helices I, B', and C were not well maintained during the SMD simulation along channel 3, as shown in Fig. 9B.
|
Discussion
|
|---|
Previous evidence from crystal structures and molecular dynamics simulations indicates that many P450s may have one or more channels between the protein exterior and active site (Wade et al., 2004; Cojocaru et al., 2007). However, it is not clear at present whether all P450s use the same channel(s) for ligand passage (Schleinkofer et al., 2005). Crystal structures of P450 3A4 with and without metyrapone bound (Williams et al., 2004; Yano et al., 2004) reveal that the heme could be accessible to solvent via two channels. One channel runs through the B-C loop and is formed by the ß1 sheet, B-B' loop, and the end of helix F'. In addition, it was speculated that appropriate movement of residues around the Phe-cluster could provide a route for a substrate from the initial recognition site to the active site (Williams et al., 2004). The purpose of this study was to determine which channel(s) and residues are involved in metyrapone egress from the active site. To this end, a 5-ns molecular dynamics simulation was performed on P450 3A4 bound with metyrapone, followed by SMD simulations on the possible exit channels.
The CMD simulation results indicate that one of two open channels found in 3A4 crystal structures, located between the ß1 sheet, the B-B' loop, and the end of helix F', is closed because of the movement of residues in this channel. In contrast, residues at the entrance of the channel threading through the B-C loop do not change significantly, and the channel remains open. Although several phenylalanine residues of the Phe-cluster rearrange to block channel 2, the entrance of channel 3 to the active site does not expand. In bacterial P450 101 (Raag et al., 1993) and P450 102 (Li and Poulos, 1997), a channel located between the F-G loop and the ß1 sheet was suggested to be a route for ligand passage. This channel is similar to channel 2 in the present study. Molecular dynamics simulations on both P450 102 (Paulsen and Ornstein, 1995; Chang and Loew, 1999) and P450 101 (Ludemann et al., 2000b) demonstrated that residues at the mouth of this channel underwent a conformational change. As a result, the width of the mouth of this channel was slightly larger than its value in the crystal structures, which readily allows a ligand to pass through. In the current study, however, the width of the mouth of channel 2 is reduced. Several residues rearrange to form a new hydrogen bonding network, resulting in closing of channel 2. The present CMD simulation indicates a low possibility for channel 2 to serve as an exit channel. A similar conclusion has been drawn from a site-directed mutagenesis study on rat P450 2B1 (Scott et al., 2002), in which a number of residues in the F-G region, similar to channel 2, were substituted with only minimal effects on the kinetic parameters for oxidation of three substrates. Comparison of bacterial P450 101 and 102 with mammalian P450 2B1 and 3A4 reveals that the longer length of the polypeptide chain between helices F and G forms two short helices F' and G' in 3A4 and 2B1. The two short helices occupy the channel region between the F-G loop and ß1 sheet in P450 101 and 102. This might be one of reasons the region between the F-G loop and ß1 sheet cannot serve as a route for ligand passage in mammalian P450s.
Inhibitor exit from the active site along channel 1 was explored by the SMD method. The time scale of natural unbinding of metyrapone is not attainable by conventional molecular dynamics simulations, which are limited to a nanosecond time scale. SMD provides a means of accelerating the unbinding process by applying an external force on the ligand. By monitoring the forces applied and the response of the ligand, one can characterize the pathway and clarify the role of certain residues during the egress. SMD simulation along channel 1 indicates that secondary structures of helices I, B', and C are well maintained, and only small backbone motion is required for inhibitor egress. The potential for channel 1 to serve as a ligand channel is also supported by the observations from both crystal structures of P450 51 (Podust et al., 2001) and molecular dynamics simulations on P450 107A1 (Winn et al., 2002), 2C5 (Wade et al., 2004), and 2B1 (Li et al., 2005). The P450 51 structure has a large open channel 1 and a closed channel 2. The suggestion has been made that channel 2 must be closed for channel 1 to open in P450 51 (Podust et al., 2001). This suggestion is in good agreement with our findings. Our SMD simulation results, coupled with previous evidence from both crystal structures and molecular dynamics simulations, demonstrate that ligand can exit from the active site via channel 1.
Analysis of the egress trajectory revealed that F108 and I120 appear to act as gatekeepers to prevent the inhibitor from leaving the active site. Both residues formed hydrophobic interactions with the inhibitor during the egress process. After the inhibitor reached the entrance of the binding pocket, the side chain of F108 underwent a rotation to expand the volume for the inhibitor to pass. To quantify the opening and rotation, the distance between F108 and I120 and the side-chain torsions 
1 (N, C, Cß, C
) of F108 have been monitored, as shown in Fig. 10. In the first 200 ps, the distance remains around 0.8 nm. In the following 120 ps, the distance significantly increases to 1.1 nm in order for metyrapone to pass through. At the same time, the torsion of F108 increases from 55° to a maximal 120°. After the inhibitor crossed over the bottleneck, the aromatic ring of F108 completely flipped up. A similar phenomenon has also been observed in the SMD simulations on P450 101 with camphor (Ludemann et al., 2000b). In that simulation, the two phenylalanine residues F193 and F87 were flipped up for the substrate to exit the active site. The same event was also observed in other enzymes, and the detailed gating mechanism was discussed by McCammon's group (Zhou et al., 1998). Based on these findings, substitution of F108 and I120 by other residues could change the substrate selectivity and regio- and stereospecificities of P450 3A4. This suggestion is in agreement with the mutagenesis results that replacement of F108 and I120 by the larger residue Trp led to an increase in the ratio of 1'-OH/4'-OH midazolam, whereas substitutions with the small residue Ala decreased this ratio (Khan et al., 2002). Another site-directed mutagenesis study indicated that replacement of F108 with Leu resulted in a significant switch of P450 3A4 regioselectivity toward that of P450 3A5 in aflatoxin B1 hydroxylation (Wang et al., 1998).
View larger version (24K):
[in this window]
[in a new window]
|
FIG. 10. Top, the variation of distance between residue F108 and I120 versus simulation time. Bottom, the variation of the dihedral angle 1 of F108 versus simulation time.
|
|
Channel 3 was demonstrated in the molecular dynamics simulations of bacterial P450 101 (Winn et al., 2002). In this previous simulation, although channel 3 was randomly chosen to be a possible exit channel for the substrate camphor, the results indicated a quite low possibility for substrate egress. In addition, this channel could not be observed for product egress. In the current simulations (Fig. 6), the maximum rupture force required for inhibitor egress via channel 3 was 680 pN, which is greater than that of channel 1. Although the rupture forces obtained in the current nanosecond-scale simulations do not reproduce the real rupture force in the natural binding and unbinding process, which usually occur in the milliseconds or longer, the relative values should be meaningful. To compare the protein motions for inhibitor unbinding via these two channels, the RMSDs per residue from the equilibration model are reported in Fig. 8. The two channels exhibit similar RMSD profiles except for three shaded regions indicating residues in the B-C loop, helices F through G region, and helix I, respectively. Among these three regions, inhibitor egress via channel 3 displays a larger displacement than that via channel 1. The results show that egress via channel 3 requires larger protein motion than via channel 1. Analysis of the egress trajectory revealed that egress via channel 3 strongly affected the protein conformation. In particular, secondary structures of helices I and B' appeared to be perturbed, several residues had to be dislodged in order for the inhibitor to pass the roof of the binding pocket, and residues F108, F189, and E308 had to be completely reoriented. It is the rearrangement of these residue side chains that results in the change in secondary structures of helices I and B'. Compared with inhibitor egress via channel 1, the high rupture force and the significant conformational change required for channel 3 indicate a low probability for the latter to serve as an exit pathway for metyrapone.
In conclusion, conventional molecular dynamics simulation and steered molecular dynamics simulations results showed that channel 1 can be used for metyrapone egress, but they do not support a role for channel 2 or channel 3. Residues F108 and I120 play an important role in inhibitor egress. Considering the fact that conformations of P450s have high plasticity upon binding of different ligands, as well as findings from molecular dynamics simulations of other mammalian P450s, the possibility of channel 2 and/or channel 3 serving as a ligand egress pathway cannot be ruled out. This idea is supported by the latest resolved crystal structures of P450 3A4 with the much larger ligands ketoconazole and erythromycin (Ekroos and Sjogren, 2006).
|
Acknowledgments
|
|---|
We thank the molecular dynamics group at Department of Biophysics Chemistry of Groningen University for their kindness in offering us the GROMACS program (www.gromacs.org).
|
Footnotes
|
|---|
The authors gratefully acknowledge financial support from the State Key Program of Basic Research of China (Grant 2002CB512802), National Natural Science Foundation of China (Grants 30672539, 20572117, and 20472094), the Basic Research Project for Talent Research Group from the Science and Technology Commission of Shanghai, the Key Project from the Science and Technology Commission of Shanghai (Grants 02DJ14006 and 03ZR14110), the Key Project for New Drug Research from the Chinese Academy of Sciences, and the 863 Hi-Tech Programme (Grants 2006AA020602 and 2006AA02Z336). These studies were also supported by Grant GM54995 and Center Grant ES06676 from the National Institutes of Health of the United States.
Article, publication date, and citation information can be found at http://dmd.aspetjournals.org.
doi:10.1124/dmd.106.014019.
ABBREVIATIONS: P450, cytochrome P450; SMD, steered molecular dynamics; MD, molecular dynamics; RMSD, root-mean-square deviation; CMD, conventional molecular dynamics; PDB, Protein Data Base.
Address correspondence to: Dr. Hualiang Jiang and Dr. Hong Liu, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, 555 Zu Chong Zhi Road, Zhangjiang Hi-Tech Park, Shanghai 201203, China. E-mail: hliu{at}mail.shcnc.ac.cn
|
References
|
|---|
Berendsen HJC, Postma JPM, van Gunsteren WF, and Hermans J (1981) Interaction models for water in relation to protein hydration, in Intermolecular Forces (Pullman B ed), pp 331–342. D. Reidel Publishing Company, Dordrecht.
Binning G, Quate CF, and Gerber G (1986) Atomic force microscope. Phys Rev Lett 56: 930–933.[CrossRef][Medline]
Breneman CM and Wiberg KB (1990) Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis. J Comp Chem 11: 361–397.
Chang YT and Loew GH (1999) Molecular dynamics simulations of P450 BM3—examination of substrate-induced conformational change. J Biomol Struct Dyn 16: 1189–1203.[Medline]
Cojocaru V, Winn PJ, and Wade RC (2007) The ins and outs of cytochrome P450s. Biochim Biophys Acta 1770: 390–401.[Medline]
Cupp-Vickery JR and Poulos TL (1995) Structure of cytochrome P450eryF involved in erythromycin biosynthesis. Nat Struct Biol 2: 144–153.[CrossRef][Medline]
Darden T, York D, and Pedersen L (1993) Particle mesh Ewald.: an N-log(N) method for Ewald sums in large systems. J Chem Phys 98: 10089–10092.
Davydov DR, Halpert JR, Renaud JP, and Hui Bon Hoa G (2003) Conformational heterogeneity of cytochrome P450 3A4 revealed by high pressure spectroscopy. Biochem Biophys Res Commun 312: 121–130.[CrossRef][Medline]
Ekroos M and Sjogren T (2006) Structural basis for ligand promiscuity in cytochrome P450 3A4. Proc Natl Acad Sci USA 103: 13682–13687.[Abstract/Free Full Text]
Essmann U, Perera L, Berkowitz ML, Darden T, Lee H, and Pedersen LG (1995) A smooth particle mesh Ewald potential. J Chem Phys 103: 8577–8592.
Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Zakrzewski VG, Montgomery JA, Stratmann RE, Burant JC, et al. (1998) Gaussian98, Revision A. 7. Gaussian, Inc, Pittsburgh, PA.
Grubmuller H, Heymann B, and Tavan P (1996) Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force. Science (Wash DC) 271: 997–999.[Abstract]
Guengerich FP (1999) Cytochrome P450 3A4: regulation and role in drug metabolism. Annu Rev Pharmacol Toxicol 39: 1–17.
Harlow GR and Halpert JR (1998) Analysis of human cytochrome P450 3A4 cooperativity: construction and characterization of a site-directed mutant that displays hyperbolic steroid hydroxylation kinetics. Proc Natl Acad Sci USA 95: 6636–6641.[Abstract/Free Full Text]
He YA, Roussel F, and Halpert JR (2003) Analysis of homotropic and heterotropic cooperativity of diazepam oxidation by CYP3A4 using site-directed mutagenesis and kinetic modeling. Arch Biochem Biophys 409: 92–101.[CrossRef][Medline]
Hess B, Bekker H, Fraaije J, and Berendsen HJ (1997) A linear constraint solver for molecular simulations. J Comp Chem 18: 1463–1472.
Isralewitz B, Gao M, and Schulten K (2001) Steered molecular dynamics and mechanical functions of proteins. Curr Opin Struct Biol 11: 224–230.[CrossRef][Medline]
Khan KK, He YQ, Domanski TL, and Halpert JR (2002) Midazolam oxidation by cytochrome P450 3A4 and active-site mutants: an evaluation of multiple binding sites and of the metabolic pathway that leads to enzyme inactivation. Mol Pharmacol 61: 495–506.[Abstract/Free Full Text]
Koley AP, Buters JT, Robinson RC, Markowitz A, and Friedman FK (1997) Differential mechanisms of cytochrome P450 inhibition and activation by alpha-naphthoflavone. J Biol Chem 272: 3149–3152.[Abstract/Free Full Text]
Li H and Poulos TL (1997) The structure of the cytochrome p450BM-3 haem domain complexed with the fatty acid substrate, palmitoleic acid. Nat Struct Biol 4: 140–146.[CrossRef][Medline]
Li W, Liu H, Scott EE, Grater F, Halpert JR, Luo X, Shen J, and Jiang H (2005) Possible pathway(s) of testosterone egress from the active site of cytochrome P450 2B1: a steered molecular dynamics simulation. Drug Metab Dispos 33: 910–919.[Abstract/Free Full Text]
Lindahl E, Hess B, and van der Spoel D (2001) GROMACS 3.0: a package for molecular simulation and trajectory analysis. J Mol Model 7: 306–317.
Ludemann SK, Lounnas V, and Wade RC (2000a) How do substrates enter and products exit the buried active site of cytochrome P450cam? 1. Random expulsion molecular dynamics investigation of ligand access channels and mechanisms. J Mol Biol 303: 797–811.[CrossRef][Medline]
Ludemann SK, Lounnas V, and Wade RC (2000b) How do substrates enter and products exit the buried active site of cytochrome P450cam? 2. Steered molecular dynamics and adiabatic mapping of substrate pathways. J Mol Biol 303: 813–830.[CrossRef][Medline]
Paulsen MD and Ornstein RL (1995) Dramatic differences in the motions of the mouth of open and closed cytochrome P450BM-3 by molecular dynamics simulations. Proteins 21: 237–243.[CrossRef][Medline]
Podust LM, Poulos TL, and Waterman MR (2001) Crystal structure of cytochrome P450 14alpha-sterol demethylase (CYP51) from Mycobacterium tuberculosis in complex with azole inhibitors. Proc Natl Acad Sci USA 98: 3068–3073.[Abstract/Free Full Text]
Poulos TL, Finzel BC, and Howard AJ (1987) High-resolution crystal structure of cytochrome P450cam. J Mol Biol 195: 687–700.[CrossRef][Medline]
Raag R, Li H, Jones BC, and Poulos TL (1993) Inhibitor-induced conformational change in cytochrome P450CAM. Biochemistry 32: 4571–4578.[CrossRef][Medline]
Schleinkofer K, Sudarko, Winn PJ, Ludemann SK, and Wade RC (2005) Do mammalian cytochrome P450s show multiple ligand access pathways and ligand channelling? EMBO Rep 6: 584–589.[Medline]
Scott EE and Halpert JR (2005) Structures of cytochrome P450 3A4. Trends Biochem Sci 30: 5–7.[CrossRef][Medline]
Scott EE, He YQ, and Halpert JR (2002) Substrate routes to the buried active site may vary among cytochromes P450: mutagenesis of the F-G region in P450 2B1. Chem Res Toxicol 15: 1407–1413.[CrossRef][Medline]
Tang W and Stearns RA (2001) Heterotropic cooperativity of cytochrome P450 3A4 and potential drug-drug interactions. Curr Drug Metab 2: 185–198.[CrossRef][Medline]
van Gunsteren WF, Billeter SR, Eising AA, Hunenberger PH, Kuger P, Mark AE, Scott WRP, and Tironi IG (1996) The GROMOS96 Manual and User Guide. Biomos, Zurich, Switzerland.
Wade RC, Winn PJ, Schlichting I, and Sudarko (2004) A survey of active site access channels in cytochromes P450. J Inorg Biochem 98: 1175–1182.[CrossRef][Medline]
Wallace AC, Laskowski RA, and Thornton JM (1995) LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions. Protein Eng 8: 127–134.[Abstract/Free Full Text]
Wang H, Dick R, Yin H, Licad-Coles E, Kroetz DL, Szklarz G, Harlow G, Halpert JR, and Correia MA (1998) Structure-function relationships of human liver cytochromes P450 3A: aflatoxin B1 metabolism as a probe. Biochemistry 37: 12536–12545.[CrossRef][Medline]
Williams PA, Cosme J, Vinkovic DM, Ward A, Angove HC, Day PJ, Vonrhein C, Tickle IJ, and Jhoti H (2004) Crystal structures of human cytochrome P450 3A4 bound to metyrapone and progesterone. Science (Wash DC) 305: 683–686.[Abstract/Free Full Text]
Winn PJ, Ludemann SK, Gauges R, Lounnas V, and Wade RC (2002) Comparison of the dynamics of substrate access channels in three cytochrome P450s reveals different opening mechanisms and a novel functional role for a buried arginine. Proc Natl Acad Sci USA 99: 5361–5366.[Abstract/Free Full Text]
Xu Y, Shen J, Luo X, Silman I, Sussman JL, Chen K, and Jiang H (2003) How does huperzine A enter and leave the binding gorge of acetylcholinesterase? Steered molecular dynamics simulations. J Am Chem Soc 125: 11340–11349.[CrossRef][Medline]
Yano JK, Wester MR, Schoch GA, Griffin KJ, Stout CD, and Johnson EF (2004) The structure of human microsomal cytochrome P450 3A4 determined by X-ray crystallography to 2.05-Å resolution. J Biol Chem 279: 38091–38094.[Abstract/Free Full Text]
Zhou HX, Wlodek ST, and McCammon JA (1998) Conformation gating as a mechanism for enzyme specificity. Proc Natl Acad Sci USA 95: 9280–9283.[Abstract/Free Full Text]
![[Feedback]](/icons/banner/feedbackACT.gif)
![[For Subscribers]](/icons/banner/subscriberACT.gif)
![[Archive]](/icons/banner/archiveACT.gif)
![[Search]](/icons/banner/searchACT.gif)
![[Contents]](/icons/banner/tocACT.gif)
Top
Abstract
Materials and Methods
