Abstract
Molecular dynamics (MD) simulations of 7-ethoxy- and 7-methoxyresorufin bound in the active site of cytochrome P450 (P450) 1A2 wild-type and various mutants were used to predict changes in substrate specificity of the mutants. A total of 26 multiple mutants representing all possible combinations of five key amino acid residues, which are different between P450 1A1 and 1A2, were examined. The resorufin substrates were docked in the active site of each enzyme in the productive binding orientation, and MD simulations were performed on the enzyme-substrate complex. Ensembles collected from MD trajectories were then scored on the basis of geometric parameters relating substrate position with respect to the activated oxoheme cofactor. The results showed a high correlation between the previous experimental data on P450 1A2 wild-type and single mutants with respect to the ratio between 7-ethoxyresorufin-O-deethylase (EROD) and 7-methoxyresorufin-O-demethylase (MROD) activities and the equivalent in silico “E/M scores” (the ratio of hits obtained with 7-ethoxyresorufin to those obtained with 7-methoxyresorufin). Moreover, this correlation served to establish linear regression models used to evaluate E/M scores of multiple P450 1A2 mutants. Seven mutants, all of them incorporating the L382V substitution, were predicted to shift specificity to that of P450 1A1. The predictions were then verified experimentally. The appropriate P450 1A2 multiple mutants were constructed by site-directed mutagenesis, expressed in Escherichia coli, and assayed for EROD and MROD activities. Of six mutants, five demonstrated an increased EROD/MROD ratio, confirming modeling predictions.
Cytochromes P450 constitute a large family of heme-thiolate monooxygenase enzymes widely found in nature (Nelson et al., 1996). These enzymes are capable of catalyzing the oxidation of a wide variety of both xenobiotic and endogenous compounds. However, even closely related isoforms may exhibit different catalytic activities.
Human P450 1A1 and 1A2, the two major isoforms in the P450 1A subfamily, share 72% amino acid sequence identity but display different substrate specificities. P450 1A1 prefers to metabolize benzo-[a]pyrene and other polycyclic aromatic hydrocarbons, whereas P450 1A2 favors the oxidation of heterocyclic and aromatic amines (Kawajiri and Hayashi, 1996; Guengerich, 2005). Likewise, they also exhibit different substrate specificities with resorufin substrates such as 7-ethoxyresorufin and 7-methoxyresorufin (Nerurkar et al., 1993; Burke et al., 1994). Therefore, the P450 1A1/1A2 system provides a good model for exploring the basis for functional differences between individual P450 enzymes.
The experimentally determined structure of a protein can provide valuable insight into its function. Homology modeling is an alternative method for obtaining the structure when the crystal structure is not available (Szklarz et al., 2000). Molecular dynamics (MD) simulations on the enzyme-substrate complex could provide information on whether the substrate is sterically accessible for oxidation, which could be used for predicting the activities of different isoforms. Indeed, this type of approach has been applied to predict the catalytic activities, coupling/uncoupling efficiency, and substrate specificity of several bacterial and eukaryotic P450s and their mutants (Harris and Loew, 1995; Paulsen et al., 1996; Kuhn et al., 2001). Likewise, with a homology model of P450 1A1 (Szklarz and Paulsen, 2000), we have successfully used MD simulations on the enzyme-substrate complex to evaluate the productive binding orientations of benzo[a]pyrene and fatty acids (Ericksen and Szklarz, 2005). This application not only rationalized the metabolite profiles of several P450 1A1 substrates but also indicated that the steric influence could be quantitatively evaluated using the MD-based steric scoring approach.
Docking of resorufin substrates in the active site of the homology model of P450 1A1 suggested that the residue Val-382 may play an important role in binding and O-dealkylation of these substrates (Szklarz and Paulsen, 2002). This suggestion was confirmed experimentally in further studies with P450 1A1 V382L and V382A mutants (Liu et al., 2003). We have also examined the function of five key residues that are different between P450 1A1 and 1A2 and thus may affect enzyme specificity (Liu et al., 2004). Indeed, residue replacement resulted in different levels of specificity changes in reciprocal P450 1A1 and 1A2 single mutants, as indicated by the ratio of 7-methoxyresorufin-O-dealkylation (MROD) to 7-ethoxyresorufin-O-dealkylation (EROD) (EROD/MROD ratio). However, no single reciprocal mutation was capable of entirely conferring the activities of one isoform onto another, and only residue substitutions at position 382 in both P450 1A1 and 1A2 led to the shift in specificity.
The objective of the present study was to investigate whether any multiple mutations in P450 1A2 may convert the resorufin specificity characteristic of P450 1A2 to that of 1A1. A total of 26 multiple mutants representing all possible combinations of five key amino acid residues, which are different between P450 1A1 and 1A2, were evaluated using molecular modeling methods, such as molecular dynamics, and modeling predictions were verified experimentally. This type of approach may help in the determination of P450 specificity.
Materials and Methods
Materials. The DNA sequencing and mutagenic oligonucleotide primers were synthesized at Sigma-Genosys (The Woodlands, Texas). A QuickChange site-directed mutagenesis toolkit and StrataPrep plasmid miniprep kit were purchased from Stratagene (La Jolla, CA). Escherichia coli DH5α competent cells were obtained from Invitrogen (Carlsbad, CA). 7-Methoxyresorufin, 7-ethoxyresorufin, resorufin, NADPH, ampicillin, isopropyl-β-d-thiogalactopyranoside, δ-aminolevulinic acid, CHAPS, dilauroyl-l-3-phosphatidyl choline, and phenylmethanesulfonyl fluoride were from Sigma-Aldrich (St. Louis, MO). Nickel-nitrilotriacetic acid-agarose and a gel extraction kit were purchased from QIAGEN (Valencia, CA). All other chemicals used were of analytical grade and were obtained from standard commercial sources.
Construction of Multiple Mutants of P450 1A2. The clones of P450 1A1 WT, P450 1A2 WT, and 1A2 single mutants (T124S, T223N, V227G, N312L, and L382V) were constructed previously (Liu et al., 2003, 2004). Double to quintuple P450 1A2 multiple mutants were created using pCWori+ plasmids containing the His-tag P450 1A2 WT or its single mutant cDNA as templates (Liu et al., 2004). A QuickChange site-directed mutagenesis kit was used, and mutagenic primers were the same as those used previously for the construction of single P450 1A2 mutants (Liu et al., 2004). A typical protocol described in the QuikChange mutagenesis instruction manual was used with minor modifications. The polymerase chain reaction conditions were adjusted as follows: 1 cycle at 95°C for 1 to 2 min, followed by 20 to 25 cycles of temperature cycling-denaturing at 95°C for 30 s, annealing at 55 to 65°C for 40 to 60 s, and extension at 72°C for an additional minute. The cycling was followed with a final extension at 72°C for 10 min. The polymerase chain reaction product was then treated with restriction enzyme DpnI to digest the parental DNA template. Postdigestion product containing the mutated plasmid with nicked circular strands was transformed into XL 1-B supercompetent cells. The plasmid containing the mutant DNA was purified using a QIAGEN plasmid purification kit. The double mutants were constructed by adding an additional mutation to the single mutants, and the triple mutants were constructed in a similar fashion from the double mutants. The mutations in the plasmids were verified by nucleotide sequence analysis carried out at the Molecular Genetics Instrumentation Facility, University of Georgia (Athens, GA).
Protein Expression and Purification. The His-tag-containing P450 1A1 and 1A2 WT and 1A2 mutants were expressed in E. coli DH5α cells and purified as previously described (Liu et al., 2004), except that during the purification of P450 1A2 WT and mutants, 0.5 M KCl was added to the solubilization buffer. The final purity of the enzymes was assessed by SDS-polyacrylamide gel electrophoresis. Western blotting were performed using anti-human P450 1A1/1A2 (Oxford Biomedical Research, Oxford, MI), and P450 proteins were visualized as described (Kedzie et al., 1991). P450 content was determined by reduced CO/reduced difference spectra (Omura and Sato, 1964), and protein was measured by the method of Lowry et al. (1951) using the Folin phenol reagent.
P450 Activity Assays. 7-Methoxy- and 7-ethoxyresorufin O-dealkylase activities of P450 1A1 WT, 1A2 WT and their mutants were assayed at 37°C by fluorometric detection of resorufin as described previously, using excitation and emission wavelengths of 550 and 585 nm, respectively (Liu et al., 2003, 2004). The reaction mixture contained 50 nM P450, 100 nM P450 reductase, and 10 μM substrate in a 0.1 M potassium phosphate buffer, and the reaction was initiated by the addition of 10 μl of 100 mM NADPH, as described previously (Liu et al., 2003, 2004). Rat cytochrome P450 reductase was expressed in E. coli and purified as described (Liu et al., 2003, 2004). The formation of resorufin was detected as the increase of fluorescence intensity against time. The reaction rate was quantified with resorufin standards. All measurements were performed in triplicate.
Molecular Modeling Methods: General. Molecular modeling was performed on an SGI Octane workstation using InsightII software (Accelrys, San Diego, CA). The homology model of P450 1A2 was constructed previously (Liu et al., 2004). The crystal structure of P450 1A2 (Protein Data Bank code: 2hi4.pdb) was obtained courtesy of Dr. Eric F. Johnson, The Scripps Research Institute (La Jolla, CA) (Sansen et al., 2007). Minimization and MD simulations were performed using the Discover module of InsightII, with the consistent valence force field supplemented with parameters for heme and ferryl oxygen (Paulsen and Ornstein, 1991, 1992). Trajectory data were extracted using the Analysis module of InsightII and graphed with Origin 8 (OriginLab Corporation, Northampton, MA).
MD Simulations with P450 1A2 WT and Mutants. The structures of P450 1A2 mutants were obtained by the replacement of selected amino acids in the homology model of P450 1A2 followed by 500 iterations of energy minimization using the steepest descent gradient, essentially as described earlier (Liu at al., 2003, 2004). 7-Methoxy- and 7-ethoxyresorufins were previously docked into the active site of P450 1A2 using the Affinity module of InsightII (Liu et al., 2004). These enzyme-alkoxyresorufin complexes served as the starting point for the subsequent dynamic docking phase, which refined the position of the substrate in the productive binding orientation.
Dynamic docking involved 5-ps constrained MD simulations at 300 K. During that period, only the substrate and protein residues within 15 Å from the substrate were allowed to move. The nonbond interaction cutoff was set at 15 Å, a residue-based nonbond list cutoff was 16 Å, and the distance-dependent dielectric was used. The distance between the ferryl oxygen and one of the hydrogens to be abstracted was restrained to 2.75 to 3.25 Å with a restoring force (k = 100 kcal·mol-1·Å-1) applied to keep the oxidation site close to the ferryl oxygen. Another distance restraint of 4.10 to 4.35 Å between the carbon atom at the oxidation site and the ferryl oxygen (C···O) was applied concurrently to keep the angle (C—H···O) as close as possible to linear (Kuhn et al., 2001).
In the second part of MD simulations, all of the restraints were abolished and the substrate was allowed to move freely for the subsequent 50 ps. Simulation conditions were exactly as described above.
Sampling and Scoring of Productive Binding Orientations. The coordinates of each enzyme-substrate complex were saved every 250 fs, giving the total of 120 MD frames obtained after the initial 10 ps of dynamic docking. All trajectories were examined with the Analysis module of InsightII. To score the likelihood of hydroxylation, each sampled frame of a given enzyme-substrate complex was evaluated using different sets of geometric criteria.
The productive orientation of the substrate is usually defined by two parameters: 1) the distance, r, between the hydrogen atom at the oxidation site of the substrate, and the ferryl oxygen atom (H···O distance), and 2) the angle, θ, between the ferryl oxygen, hydrogen atom to be abstracted, and the carbon at the oxidation site (C—H···angle). Because the linear transition state is assumed for hydrogen abstraction to occur, ideally, the distance should be close to 3 Å and the angle is postulated to be 180° (Fig. 1A).
In this work, we used four types of criteria to define an effective score or “hit.” Criterion P defined a hit as a frame in which the distance r was equal to or smaller than3Å(r ≤ 3.0 Å), whereas criterion PM combined the condition on distance r with the condition on angle θ, required to be larger than 135° (r ≤ 3.0 Å and θ ≥ 135°). Less stringent criterion Q required the distance r not to exceed 3.5 Å, and criterion QN added the requirement on the angle not to exceed 120° (r ≤ 3.5 Å and θ ≥ 120°). An example of such an evaluation is shown in Fig. 1, C and D, for which criterion QN is used to assess hits, represented by frames located within the gray region.
This procedure was applied to all single and multiple mutants complexed with either 7-ethoxy- or 7-methoxyresorufin. The “E/M score” representing in silico specificity was defined as the ratio of hits obtained with 7-ethoxyresorufin to those obtained with 7-methoxyresorufin.
Prediction of Substrate Specificity of P450 1A2 Multiple Mutants. The E/M scores obtained for single mutants were used to construct a linear regression model that linked these results with experimental specificity expressed as the EROD/MROD ratio defined previously (Liu et al., 2004). The best correlation between the two was seen for criterion P, with R2 equal to 0.9612 (p < 0.001), and for criterion QN, with R2 = 0.9097 (p < 0.001). Therefore, these two criteria were used for prediction of specificity in multiple 1A2 mutants. The specific equation for criterion P is
The specific equation for criterion QN is
where specificity was expressed as the predicted EROD/MROD ratio.
Comparisons with the Crystal Structure of P450 1A2. The homology model of P450 1A2 has been compared with the recently available crystal structure of this enzyme (Sansen et al., 2007). The RMSD between the crystal and the model has been measured by superimposing the backbones of the two structures using different regions as a basis. These included 1) overall structure, 2) the active site region, and 3) the previously mentioned five key residues only. The active site was defined as a region containing residues within 10 Å from docked 7-ethoxyrersorufin, and key residues were T124, T223, V227, N312, and L382.
7-Ethoxyresorufin and 7-methoxyresorufin were placed in the crystal structure of P450 1A2 WT in positions similar to those in the homology model. MD simulations of enzyme-substrate complexes were performed in the same way as with the homology model, and the results were analyzed as described above, by counting hits using different geometric criteria.
Results
Substrate Mobility in the Active Site of P450 1A2 WT and Single Mutants. Our previous studies with P450 1A1 and 1A2 reciprocal active site mutants demonstrated that the replacement of a single residue may affect the specificity of alkoxyresorufin metabolism (Liu et al., 2004). In particular, mutations at position 382 in both P450 1A1 and 1A2 shifted substrate specificity from one enzyme to another. However, none of the single residue substitutions at positions 124, 223, 227, and 312 (location in the 1A2 sequence) completely interconverted P450 1A1 and 1A2 specificities (Liu et al., 2004). Therefore, the goal of our present studies was to identify multiple mutants in which the specificity of one enzyme was fully conferred to another.
In the first part of this project, we used molecular modeling methods, in particular MD, to make predictions concerning the results of multiple residue substitutions in P450 1A2. Initially, P450 1A2 WT and its five single mutants, which were investigated earlier (Liu et al., 2004), were selected as the training set for MD simulations to provide information regarding specificity prediction.
Molecular dynamics simulations were used both as a part of the docking strategy to refine the substrate orientation in the active site and as a means to predict the product profile, similar to our previous work (Ericksen and Szklarz, 2005). Thus, after the energy minimization, the first 5 ps of restrained MD served to optimally place an alkoxyresorufin in a productive binding orientation. Figure 2A shows the resulting orientation of 7-methoxy- and 7-ethoxyresorufin in the active site of P450 1A2, with five key residues (T124, T223, V227, N312, and L382) displayed.
The productive binding orientation has been defined by two geometric parameters, r and θ, as illustrated in Fig. 1A (see also Materials and Methods). During 5 ps of dynamic docking, the distance, r, between the ferryl oxygen of the heme and the hydrogen atom of the substrate to be abstracted was restrained to ∼3.0 Å and remained fairly constant (see Fig. 1B). In contrast, the angle, θ, between the ferryl oxygen, the hydrogen atom to be abstracted, and the carbon at the oxidation site alternated between 120 and 180° (Fig. 1B).
The transition from the restrained state to the free moving state was achieved over 5 ps. For each enzyme-substrate complex, 120 frames from unrestrained MD (productive phase) were then collected, starting at the 10 ps time point. In most cases, after removal of the restraints, the substrate moved rapidly into another, more stable, orientation where it remained throughout the remainder of the simulation. This led to sudden and radical changes in relevant distances and angles, followed by a reasonably steady phase of equilibrium, as illustrated in Fig. 1B.
A detailed analysis of the geometrical parameters was performed for each snapshot of every enzyme-substrate complex to evaluate the tendency of the ligand to remain in the productive binding orientation during the productive dynamic phase. Each enzyme-substrate complex displayed a unique characteristic MD “fingerprint” regarding the distances and angles related to the productive binding orientation. Two examples are shown in Fig. 1, C and D. For instance, for the P450 L382V mutant, both hydrogens at the oxidation site of 7-ethoxyresorufin remained at approximately 3 Å from the ferryl oxygen, whereas all hydrogens of 7-methoxyresorufin were much farther away (Fig. 1D), in agreement with our earlier findings (Liu et al., 2004). Overall, in this mutant, 7-ethoxyresorufin showed a greater preference to occupy a productive binding orientation than 7-methoxyresorufin.
Specificity of 7-Methoxy and 7-Ethoxyresorufin Oxidation by P450 1A2 Single Mutants. In our previous work (Liu et al., 2004), we expressed the specificity of various P450 1A1 and 1A2 mutants in terms of the ratio of 7-ethoxyresorufin O-deethylation to 7-methoxyresorufin O-demethylation (EROD/MROD ratio). This approach allowed us to successfully evaluate alterations in activity upon residue substitution and the degree to which the activity of one enzyme was conferred upon another. Therefore, we also used the EROD/MROD ratio as a measure of enzyme specificity in the present work.
Essentially, the EROD/MROD ratio was expressed in terms of its in silico equivalent, the E/M score, which represented the number of hits for 7-ethoxyresorufin divided by the number of hits for 7-methoxyresorufin. The hits, that is, the number of conformers with the substrate occupying the productive binding orientation, were derived from the analysis of the productive phase of MD simulations using four sets of different criteria with respect to distance r and angle θ. Criteria P (r ≤ 3.0 Å) and PM (r ≤ 3.0 Å and θ ≥ 135°) were more stringent than criteria Q (r ≤ 3.5 Å) and QN (r ≤ 3.5 Å and θ ≥ 120°), and criteria PM and QN included the conditions on both distance and angle (see Materials and Methods).
Table 1 shows the number of hits for the two substrates and the resulting E/M scores for various P450 1A2 single mutants. Depending on the set of criteria used in counting hits, the E/M scores may vary widely. Usually the hit was a positive integer with very few exceptions when no hits were found, such as for the 7-methoxyresorufin-P450 1A2 N312L complex. For most mutants, E/M scores were low and similar to the scores obtained for the WT 1A2, which correlated well with experimental specificity. In contrast, the score obtained for the 1A2 L382V mutant was much higher, as was the corresponding EROD/MROD ratio, consistent with our previous findings (Liu et al., 2004). This general agreement between the in silico scores and experimental results allowed us to construct linear regression models. As shown in Fig. 3, the specificities of single 1A2 mutants as predicted by the models were very close to the experimental ones, confirming good correlations. Therefore, the models were subsequently used to predict EROD/MROD ratios for multiple mutants, as detailed under Materials and Methods.
Specificity of 7-Methoxy and 7-Ethoxyresorufin Oxidation by P450 1A2 Multiple Mutants. To determine which multiple substitutions of five key residues in P450 1A2 could potentially lead to the shift in specificity toward that of P450 1A1 WT, we performed MD simulations with all possible multiple mutants. Figure 4 shows the schematic tree structure of all possible combinations of mutations containing one substitution, namely Leu-382→ Val. The number of mutants with this particular amino acid replacement, varying from double to quintuple, was 15, and, when this procedure was repeated for all five key residues, the total number of all possible mutants, excluding repeated combinations, rose to 26. The MD simulations were performed with all these mutants using exactly the same protocol as that applied in the case of P450 1A2 WT and single mutants (see Materials and Methods). The results were then evaluated using criteria P, PM, Q, and QN, as described for the single mutants.
The E/M scores were used to classify the mutants with respect to the possibility of a significant shift in substrate specificity. The in silico descriptor of P450 1A2 L382V, a mutant known to display a considerable shift in specificity, served as a reference to select likely multiple mutant candidates. As shown in Table 2, 7 of 26 mutants had E/M scores ≥80% of that for the L382V mutant and were thus predicted to display a similarly high shift in specificity. Furthermore, these predictions were upheld with all four geometric criteria. The mutants under consideration were T124S-L382V, T223N-L382V, V227G-L382V, T124S-T223N-L382V, T124S-N312L-L382V, T223N-N312L-L382V, and N312L-L382V. Thus, all the multiple mutants predicted to significantly alter specificity had to include the substitution of Leu-382→ Val, which further confirms the exceptional role of this residue in producing specificity shifts from P450 1A2 to 1A1 and vice versa, with alkoxyresorufin substrates (Liu et al., 2004). Interestingly, the quintuple mutant, T124S-T223N-V227G-N312L-L382V, was not predicted to display altered specificity.
The MD-derived E/M scores were used to calculate the values of EROD/MROD ratios for selected 1A2 multiple mutants on the basis of two linear regression models constructed for the single mutants (see the previous section). The extrapolated ratios then served as a basis for the comparison with the experimental values.
Experimental Validation of Modeling Predictions. By use of modeling, we have significantly reduced the number of P450 1A2 multiple mutants needed for further experimental studies from 26 to 7. Six of 7 multiple mutants have been successfully created by site-directed mutagenesis and were active in biochemical assays after purification. The double mutant V227G-L382V has not been successfully expressed; its expression level was very low, and the protein was unstable.
Overall, mutant enzymes differed in expression levels and activities, as also noted previously (Liu et al., 2004). For example, the P450 content of purified preparations varied from ∼1 to 60 nmol/ml and specific content varied from 0.2 to 60 nmol of P450/mg of protein, with holoenzyme in the range of 1 to 32% (data not shown). Table 3 shows the activities and EROD/MROD ratios for various multiple mutants. Because of considerable variability in the levels of activities, the EROD/MROD ratio is the best discriminator among the mutants, in agreement with our previous findings (Liu et al., 2004). It should be noted that the EROD/MROD ratios for the P450 1A2 WT and the L382V mutant were 0.32 and 2.27, respectively (Table 3), displaying a slight difference compared with the previous values of 0.4 and 2.0 (Table 1) (Liu et al., 2004).
Figure 5 shows a comparison between the predicted and experimental EROD/MROD ratios for the selected P450 1A2 multiple mutants. As mentioned above, the predicted specificities were extrapolated from two linear regression models, which established a correlation between in silico E/M scores and experimental EROD/MROD ratios using P450 1A2 WT and its five single mutants (see Materials and Methods). Five of six multiple 1A2 mutants have demonstrated the expected high shift in specificities from P450 1A2 toward 1A1, with the exception of the double mutant, T124S-L382V (see also Table 3). Moreover, the linear model using criterion P (r ≤ 3.0 Å) seemed to give much closer predictions than that using criterion QN (r ≤ 3.5 Å and θ ≥ 120°). In the case of the latter, the predicted specificity values (EROD/MROD ratios) were all much higher than the experimental ratios (data not shown). The comparison between the experimental and predicted specificity using criterion P for the selected P450 1A2 multiple mutants is shown in Fig. 5. Interestingly, in the case of the double mutant T124S-L382V, both prediction models have yielded much higher specificity (criterion P: 1.8; criterion QN: 3.6) than its experimental value of 0.34. The experimental EROD/MROD ratio for this mutant was thus comparable with that of the P450 1A2 WT (Table 3).
Comparison of the Homology Model of P450 1A2 with the Crystal Structure. At the final stage of this project, the crystal structure of human P450 1A2 was solved, with the resolution of 1.95 Å (Sansen et al., 2007). Thus, it became necessary to compare the homology model of P450 1A2 with the crystal structure to validate our current study.
With respect to the overall structure, the homology model resembles closely its structural template, P450 2C5. This could be the inevitable drawback of the homology modeling methodology itself. The overall RMSD between the model and the P450 1A2 crystal was 4.4 Å when all residues (i.e., residues 34-513) were included in the superimposition. This seemingly large value is not that surprising, as this type of comparison includes both structurally conserved regions (SCRs) and non-SCR segments, especially loops on the protein surface, which are usually created in a random fashion with a lot of variability. Consequently, it would be much more reasonable to compare directly the core regions containing most of the SCRs, which are directly involved in the interactions with substrates; in this case, the RMSD value is reduced to 3.5 Å. Structural comparisons were further performed using additional functional sets as the criteria, including active site residues and the five key residues currently under investigation. The environment formed by only the five key residues and that produced by all active site residues (i.e., residues within 10 Å of 7-ethoxyresorufin docked in the crystal structure) had similar structural divergences of ∼2.5 and 2.6 Å, respectively, much lower than a corresponding value for the global protein structure. This indicates that the active site containing the five key residues is well constructed in the homology model and comparable with the crystal structure, as shown in Fig. 2B. Moreover, in either option of superimposition, i.e., when only key residues or all active site residues were superimposed, all of the five key residues except N312 display reasonably good structural similarity, as indicated by small values of backbone RMSD, which were as follows: for T124, 1.1 and 1.92 Å, respectively; for T223, 2.4 and 2.22 Å; for V227, 1.44 and 1.68 Å; for N312, 3.92 and 4.09 Å; and for L382, 2.11 and 1.53 Å (Fig. 2B). Therefore, N312 may not be properly positioned in the homology model, which might have resulted in artificially high and extremely variable E/M scores observed earlier with the single mutant P450 1A2 N312L (see Table 1).
The dynamic behavior of substrates bound in the active site of the crystal structure has been also investigated using the protocol established for the homology model. As shown in Table 4, regardless of the criteria used, the number of hits for 7-methoxyresorufin was significantly higher than that for 7-ethoxyresorufin, in full agreement with experimental data. However, the E/M scores are similar in both the model and the crystal only in the case of criterion QN (r ≤ 3.5 Å and θ ≥ 120°). This may indicate the inaccuracy of the homology model and/or the specific conditions applied in the computer simulation and data analysis.
Discussion
In the current work, we used molecular modeling methods to predict the changes in specificity of P450 1A2 toward alkoxyresorufin substrates upon introduction of multiple residue substitutions simultaneously in the active site. The objective was 2-fold: to determine which residue combinations could completely interconvert the specificity from P450 1A2 to 1A1 and to develop a fast computational method for prediction of substrate specificity. Among 26 possible multiple mutants, ranging from double to quintuple, the molecular dynamics-based scoring method predicted 7 of them to shift specificity. Five predictions were further confirmed experimentally by site-directed mutagenesis and biochemical assays.
These studies follow our previous work concerning the importance of five active site residues that are different between P450 1A1 and 1A2 (Liu et al., 2004). On the basis of those findings, only mutations at position 382 shifted specificity from one enzyme to another using alkoxyresorufins as substrates. Similarly, our present results indicate that all mutations that likewise affect specificity contain a substitution at this position (Table 2). The only exception was the quintuple mutant, T124S-T223N-V227G-N312L-L382V, which, contrary to expectations, was not predicted to display altered specificity. Moreover, no multiple P450 1A2 mutant exhibited substrate specificity characteristic of P450 1A1. This result further strengthens the idea that the specificity is not simply determined by a linear combination of these five substitutions but that other residues, outside of the active site, must be taken into account. Similar conclusions were reported in earlier studies with single mutants of other related isoforms, such as P450 2B1 and 2B2 (He et al., 1994) and 2B4 and 2B5 (Szklarz et al., 1996). Moreover, in the P450 2B4 and 2B5 system, multiple active site mutations were not sufficient to fully interconvert activities (He et al., 1996), analogous to our studies.
In the case of P450 1A1 and 1A2, a number of residues have been identified by others as important for activity. However, in most cases, the residues studied were different from those we investigated here and in earlier work (Liu et al., 2003, 2004). For example, more recently, Kim and Guengerich (2004) reported that a triple mutant, E163-V193M-K170Q, was five times faster than the wild-type 1A2 in 7-methoxyresorufin O-demethylation. These residues, however, are located far from the active site. Recently, a combinatorial approach was applied by Taly et al. (2007) to discriminate between 1A1 and 1A2 using alkoxyresorufin substrates. However, the residues we identified as important (Liu et al., 2003, 2004) were not confirmed by this study, possibly because of exclusion of this region of the sequence. In a recent study, Lewis et al. (2007) described the construction of a homology model of P450 1A1 based on the structure of P450 1A2 and confirmed the importance of residues 124, 227, and 382 (numbers in 1A2) for binding of 7-ethoxyresorufin. Moreover, Met-121 (Leu-123 in 1A2) was also suggested to have some contributions in binding of this substrate (Lewis et al., 2007).
A comparison of our homology model with the recently crystallized structure of P450 1A2 (Sansen et al., 2007) showed that the active site exhibits reasonably good structural resemblance (Fig. 2B). The placement of the four residues studied (T124, T223, V227, and L382) is very similar; in contrast, the position of N312 appears to be significantly shifted in the model. The latter phenomenon is probably responsible for poor correlations between computational and experimental data concerning the N312L mutant (see Table 1) and might have affected molecular modeling results with other mutants. In the case of P450 2D6, the crystallization of this enzyme helped to explain various mutagenesis data, often confirming the conclusions drawn from previous homology models (Rowland et al., 2006). For example, Asp-301 was implicated in catalytic activity as early as 1993, on the basis of an early model of the active site region of P450 2D6 (Koymans et al., 1993), and Glu-216 was also suggested to play a role in substrate binding (Lewis et al., 1997). In general, most 2D6 models correctly located key residues, with some exceptions, such as residues Phe-481 and Phe-483 (de Groot et al., 1999). On the other hand, the latter discrepancy may reflect protein motion and can be resolved by molecular dynamic simulations (Rowland et al., 2006).
In some respects, the current work is also a continuation of our previous studies on regiospecificity of P450 1A1-mediated oxidations (Ericksen and Szklarz, 2005), in which modeling methodology was successfully applied to rationalize product profiles observed with several 1A1 substrates. This earlier approach has been expanded so that we not only used computational methodology for predictions concerning the specificity of alkoxyresorufin oxidation but also verified them experimentally. As observed with P450 1A1 oxidations, steric effects also seem to be of major importance in the case of alkoxyresorufin oxidations by P450 1A2. In this approach, we neglected substrate electronic considerations, but, nevertheless, electronic factors can also provide contributions. For example, they may account for differences in the number of hits observed for a given enzyme with 7-methoxyresorufin as a substrate versus the number of hits seen with 7-ethoxyresorufin (Table 1). However, the E/M score can serve as a basis for the comparison among different enzymes and can be used as a descriptor to correlate with experimental specificity.
The predictive power of this computationally derived model is affected by the length of the simulation and the number of samples collected for analysis. In our case, after the minimization of the enzyme-substrate complex, we applied restraints to force the oxidation site of the substrate close to the ferryl oxygen and then removed them to evaluate the stability of the productive binding orientation. This approach has been used effectively in our previous work (Ericksen and Szklarz, 2005). Molecular dynamics with distance restraints has also been successfully applied to investigate the metabolism of sirolimus by P450 3A4 (Kuhn et al., 2001). This application, we believe, circumvents the need for long simulations and gives this approach a potential to be routinely used as a rapid method to predict specificity and/or product profile.
An important factor in molecular modeling with the use of homology models is the accuracy of such a model, which depends to a large extent on the template used (Szklarz et al., 2000). In that respect, the 1A2 model used in these studies, on the basis of the structure of P450 2C5, can be expected to differ from the crystal structure of P450 1A2. As discussed above, the active site structure of this homology model was quite similar to that of the crystal. Consequently, new models of highly related enzymes, such as P450 1A1, should now be based on the crystal structure of P450 1A2. As reported recently, a new 1A1 model was superior to the previous ones (Lewis et al., 2007).
However, the question of structural “correctness” is not an easy one. Proteins display structural flexibility, and, in the case of P450s, this is confirmed by the existence of multiple structures for the same isoform that may display various degrees of structural differences, often dependent upon the ligand bound (Poulos, 2003, 2005). Thus, the basic P450 fold is very flexible, which is further supported by molecular dynamics simulations that can “map the route” from one structure to another. For example, in the case of P450 BM3, MD simulations showed that binding of substrate may induce a conformational shift, leading to a structurally distinct P450-ligand complex (Chang and Loew, 1999). A recent review by Hammes-Schiffer and Benkovic (2006) examined the link between protein conformational motion and enzyme catalysis and presented evidence for a network of coupled motions throughout the protein fold that facilitate chemical reactions. Thus, thermal motions of the enzyme, substrate, and cofactor lead to conformational sampling of configurations that provide a favorable environment for catalysis. Therefore, we should take into account P450 flexibility and reexamine the value of MD simulations.
In the current work we successfully combined MD-based computational predictions with their experimental verification, and the majority of predictions concerning changes in substrate specificity in P450 1A2 multiple mutants were confirmed. This technique may thus provide a prototype for rapid evaluation of the catalytic properties of different P450 isoforms, and its application to the crystal structures should further enhance its accuracy.
Acknowledgments
We thank Dr. Spencer S. Ericksen for helpful discussions and advice. Modeling studies were performed at the Computational Chemistry and Molecular Modeling Laboratory, Department of Basic Pharmaceutical Sciences, School of Pharmacy, West Virginia University, Morgantown, West Virginia.
Footnotes
-
This work was supported by National Institutes of Health Grants RR16440 and GM079724.
-
doi:10.1124/dmd.108.022640.
-
ABBREVIATIONS: P450, cytochrome P450; MD, molecular dynamics; MROD, 7-methoxyresorufin-O-dealkylation; EROD, 7-ethoxyresorufin-O-dealkylation; CHAPS, 3-[(3-cholamidopropyl) dimethylammonio]-1-propanesulfonate; WT, wild type; RMSD, root mean square deviation; SRC, structurally conserved region.
-
↵1 Current affiliation: Global Pharmaceutical Research and Development, Drug Metabolism, Abbott Park, Abbott, Illinois.
- Received June 1, 2008.
- Accepted August 13, 2008.
- The American Society for Pharmacology and Experimental Therapeutics