Theoretical characterization of substrate access/exit channels in the human cytochrome P450 3A4 enzyme: involvement of phenylalanine residues in the gating mechanism.

The human cytochrome P450 3A4 mono-oxygenates approximately 50% of all drugs. Its substrates/products enter/leave the active site by access/exit channels. Here, we perform steered molecular dynamics simulations, pulling the products temazepam and testosterone-6betaOH out of the P450 3A4 enzyme in order to identify the preferred substrate/product pathways and their gating mechanism. We locate six different egress pathways of products from the active site with different exit preferences for the two products and find that there is more than just one access/exit channel in CYP3A4. The so-called solvent channel manifests the largest opening for both tested products, thereby identifying this channel as a putative substrate channel. Most channels consist of one or two pi-stacked phenylalanine residues that serve as gate keepers. The oxidized drug breaks the hydrophobic interactions of the gating residues and forms mainly hydrophobic contacts with the gate. We argue that product exit preferences in P450s are regulated by protein-substrate specificity.

Introduction

Cytochrome P450 3A4 (CYP3A4) is the major drug-metabolizing enzyme in humans, affecting more than 50% of all swallowed drugs.(1) CYP3A4 catalyzes mono-oxygenation reactions such as hydroxylation and epoxidation.(2) The so-formed products are more hydrophilic and hence are more easily excreted from the body. CYP3A4 oxidizes a variety of chemical compounds and drugs with different sizes and shapes.[3,4] The active site is located at the hydrophobic core of the protein and is connected to the surface of the enzyme through access channels (Figure 1).(5) The prosthetic heme group in the active site is covalently bound to a conserved cysteine and is responsible for substrate oxidation. The substrate oxidation process follows a well-characterized catalytic cycle.[6−8] The cycle is initiated when a substrate enters the active site and removes a water molecule coordinated to the ferric heme. The heme’s iron shifts to its high-spin state, thereby facilitating the reduction of the ferric center by NADPH-cytochrome P450 reductase. Following molecular oxygen consumption and subsequent addition of two electrons from the reductase and two protons from the surrounding, the O−O bond is cleaved, and an active species is formed that oxygenates the substrate. Subsequently, the product leaves the active site through the access channels, and the enzyme resumes its resting state. Although the extent of sequence identity of the cytochrome P450 family of enzymes (P450s) is rather low, their structure is conserved(9) with one available fold for P450s in SCOP.(10) The most variable regions in P450s are the substrate recognition sites (SRSs)[11,12] that spatially form the active site and access channels. This variability could account for the different substrate specificity manifested by different P450 isoforms. Wade et al.(13) identified a structurally conserved channel in all bacterial P450s positioned between the FF′ loop, β1, and β4 sheets and termed it 2a (Figure 1); yet, interestingly, using molecular dynamics (MD) and repulsive expulsion molecular dynamics, they also showed that in three different bacterial P450s, the gating mechanism for ligand exit through channel 2a is different and depends on the specific isoform.(14) This difference can modulate protein−substrate specificity in P450s. In contrast to bacterial P450s, in mammalian P450s that are anchored to the endoplasmatic reticulum (ER) membrane, channel 2a is predicted to face the hydrophobic ER membrane, suggesting that highly hydrophobic substrates will reach the active site via this channel. Thus, substrate specificity is influenced not only by interactions at the P450 active site but also by the properties of the access routes to the binding site such as their size, shape, and physicochemical properties.

Figure 1

CYP3A4 structure including substrate channel 2a (shown in blue). The enzyme is represented in cyan cartoons, FF′ loop and F′ helix (connected to G′ helix in cyan) colored red, and β1 sheet colored orange. The heme prosthetic group is represented as white VDW spheres.

In most available CYP3A4 crystal structures, the enzyme exists in a closed state with little difference between the ligand-bound and -unbound enzyme conformations.[15,16] However, in the two latest determined X-ray structures of CYP3A4 complexed with two large inhibitors (ketoconazoles) and a substrate (erythromycin) inside of the active site,(17) the enzyme undergoes a large conformational change. The active site expands by 80% of its volume compared to that of the unbound structures, and there is a large shift in the F′G′ helices region (Figure 1) due to ligand binding. An extreme shift of this region was observed in the mammalian CYP2B4 crystal structure.(18) It is expected that the substrate/product entrance/exit will induce conformational changes that are not observed in most CYP3A4 crystal structures. Previous theoretical studies on bacterial P450s,[13,14,19] mammalian P450s,[20−24] and other related heme-containing enzymes[25,26] were carried out in order to monitor ligand egress from the enzyme’s active site. MD simulations on mammalian P450s resulted in different dominant egress pathways that either expressed different substrate specificity of the P450 isoforms or stemmed from the different ligand probes tested in the simulations.

Here, using MD and steered MD (SMD), we simulated the last step of the catalytic cycle pulling out of the human CYP3A4 active site two distinct products, temazepam (TMZ) and testosterone-6βOH (TST-OH) (Figure 2). We searched for all physically feasible pathways for product egress, their dimensions, and the key residues involved in the gating mechanism and compared the results to previous simulations of other P450 isoforms. We then attempted to predict the preferred product exit routes in CYP3A4, seeking the factors that determine this preference. We expect that the substrate channel characterization, achieved in this work for CYP3A4, will provide further insight into the biophysics of transport to (from) the enzyme.

Figure 2

Chemical structures and carbon numbering systems for the CYP3A4 products temazepam (TMZ) and testosterone−6OH (TST-OH).

Computational Methods

Setting up the System

The CYP3A4 structure was taken from the PDB(27) (code 1TQN(15)). Missing residues of 1TQN were modeled with InsightII (http://www.accelrys.com/products/). pKa values were calculated using the MCCE(28) software with DelPhi(29) as a Poisson−Boltzmann solver with a dielectric constant of 4. Hydrogen atoms were added with CHARMM.(30) The structure was solvated in a cubic box of ∼21000 TIP3P(31) water molecules, in addition to crystallographic waters in order to apply periodic boundary conditions. Chloride ions were added to neutralize the system. The structures of TMZ and TST-OH were derived from their substrate structures, diazepam and testosterone, taken from the Cambridge structural database (http://www.ccdc.cam.ac.uk/products/csd/, codes DIZPAM10 and TESTOM01, respectively). We substituted the oxidized C−H groups in diazepam and testosterone with C−OH groups. We took into account the flexibility of the TMZ phenyl ring by generating several conformations using the systematic search algorithm of the Cerius2 software. Each TMZ conformer and TST-OH were geometrically docked to the CYP3A4 active site using the PatchDock(32) software, taking into account shape complementarity. The TMZ flexibility was addressed by including different orientations of its rotatable phenyl ring in the docking procedure. Docking solutions with a distance less than 3 Å between the product hydroxyl and the heme’s iron were considered for further steps. The remaining solutions were minimized by the CHARMM software, and the lowest CHARMM energy docking solution was used as a starting point for further MD simulations.

Force Field Parameters

Heme and protein parameters were taken from CHARMM22(33) and CHARMM27(34) force fields, respectively. Fe and S related missing parameters were added from Autenrieth et al.(35) and Bathelt et al.(36) The heme’s atoms’ partial charges were calculated by a QM(B3LYP)/MM optimization at the DFT level of theory, with a sextet multiplicity as the ground state.(37) The basis sets used were an effective core potential (ECP) and the LACVP for the iron and 6-31G for the rest of the heme and C442 atoms. The heme Mulliken charges were added to the force field. The QM/MM calculation used the ChemShell package.(38) The force field parameters for TMZ were taken from a previous diazepam force field(39) with the addition of the hydroxyl moiety. TST-6OH and -16OH parameters were derived from a previously published CHARMM force field for cholesterol.(40) TMZ and testosterone products’ atomic partial charges were modified to account for the addition of the hydroxyl group by a QM/MM optimization as done with the heme.

MD Procedure

For each run, the system was initially minimized with CHARMM, with a decreasing harmonic force constraint on the enzyme. Further procedures were performed using the NAMD code.(41) The minimized system was gradually heated up to 310 K and equilibrated for 50 ps. The equilibrated system was then run for 500 ps of MD simulation followed by a 1.6 ns of steered MD (SMD(42)) pulling each product out of the active site. We applied a canonical NVT ensemble dynamics at 310 K with a 1 fs time step. The Shake algorithm was applied to fix all bond lengths involving hydrogen atoms. A final nonbonded interaction cutoff radius was 12 Å, with a switching function starting from 10 Å. Nonbonded electrostatic interactions were calculated with Particle Mesh Ewald summation.(43) Snapshots were saved each 500 steps. The SMD procedure used a pulling spring force constant of 7 kcal/mol/Å2 with a constant velocity of 0.000025 Å/fs. The chosen force constant and pulling velocity were three times larger and twice slower than those applied in previous SMD simulations of the human CYP3A4.(22) The pulling force was applied after the equilibration phase, with the spring connected between the center of mass of the product and a “dummy” atom in the general direction of the channels exit. Three channels (2b, 2e, and the solvent channel, termed S; see Figure 3) are observed to be partially open in the PDB file 1TQN, which we used as input.(15) Therefore, the directions of the SMD pulling forces in the case of the above channels were determined according to the visible channel exit. For those channels that evaded detection, 18 different vectors were tested that are in the general direction of the roof and sides of the active site, and three additional channels were detected (see ). All of the SMD pulling force directions and the SMD starting point structures are available in the Supporting Information (SI).

Figure 3

CYP3A4 secondary structures related to the identified substrate channels. (A) Secondary structure elements in CYP3A4 that define the different channels. FF′ loop, F′ and G′ helices, G helix, I helix, BC loop, and β1 and β4 sheets are colored blue, silver, purple, yellow, black, orange, and red, respectively. The rest of the enzyme is colored cyan and represented as cartoons. The heme prosthetic group is represented by white VDW spheres. (B) Substrate channels detected during SMD. Channel surfaces were calculated by MolAxis. Channels 2a, 2b, 2c, 2e, 3, and S are colored red, blue, yellow, green, orange, and purple, respectively. F′G′ helices are colored silver, and the BC loop is in black. The rest of the enzyme is colored cyan and represented as cartoons.

In addition to the SMD simulations in the presence of a product, we performed a substrate-free 1.6 ns of classical MD simulation on CYP3A4.

Results

Global CYP3A4 Stability

We performed different SMD runs to find the minimum constant velocity and force constant values that ensure low protein fluctuations but permit the escape of the product from the active site. We ended up with 12 different SMD simulations (6 distinct channels for each product) and another SMD run of TST-160H escaping through channel 2c (see ). A velocity of 0.000025 Å/fs was found as a minimum cutoff. As judged from the total energy, pressure, rmsd, and RMSF profiles (Figures S1−S4, Supporting Information), the entire system and protein structure appeared stable throughout the simulations, even when further applying an external pulling force. The average Cα rmsd’s along all of the CYP3A4-unbound, TMZ-bound, and TST-OH-bound simulations were 1.2, 1.2−1.3, and 1.8 Å, respectively, with not more than ±0.1 Å fluctuations (Figure S3, Supporting Information). Major local fluctuations observed in the unbound dynamics were in residues 255−270 and 280−290, which comprise a modeled loop, missing in several PDB files of CYP3A4.(16) In the bound SMD simulations, fluctuations occur mainly in the TST-OH simulations and are larger than TMZs. Major RMSF peaks detected upon addition of the products are at residues 100−125 and 200−250 (Figure S4a−S4d, Supporting Information), which belong to the BC loop and FG helices’ secondary structure elements, respectively (Figure 3A). These secondary structures are known to be the most flexible in P450s.(44)

Substrate Channels’ Structure and Dimensions

The substrate exit pathways are named according to Cojocaru et al.’s(5) nomenclature and include substrate channels 2a, 2b, 2c, 2e, and 3 and the solvent channel (S) (Figure 3). The channels are defined according to the secondary structures that form their entrance/exit (Figure 3A). Channel 2a is located between the FF′ loop and the β1 and β4 sheets. Channel 2b is positioned between the F′ and G′ helices and the β1 sheet and the BC loop. Channel 2c is between the GG′ helices and the BC loop. Channel 2e is surrounded by the BC loop. Channel 3 passes through the F′G′ helices, and channel S is above the I helix and between the FF′ loop and the β4 sheet.

The channel lengths were measured, starting 3 Å above the heme Fe up to the bottleneck residues that form the gates. The channels lengths are 19, 18, 17, 12, 17, and 13 Å for 2a, 2b, 2c, 2e, 3, and S channels, respectively. Thus, channels 2e and S are the shortest paths from the active site to the enzyme surface. The channels’ bottleneck opening radii were measured during the simulations using the MolAxis tool that identifies and calculates channel surfaces within macromolecules by computational geometry techniques.(45)Figure S5 (Supporting Information) depicts the bottleneck radius of all of the channels along the simulations, and Table 1 summarizes the channels’ maximum opening radii. We noticed that for both products, the S channel manifested the largest opening, and substrate channel 2c had the smallest bottleneck radius.

Table 1

Maximum Bottleneck Values for Each Substrate Channel Observed for the Egress of TMZ and TST-OH in Angstroms

channel	TMZ	TST-OH
2a	3.70	4.14
2b	3.65	3.96
2c	3.09	3.91
2e	3.87	4.20
3	3.31	4.12
S	3.98	4.47

Channel−Product Interactions

Table 2 summarizes the major residues which interact with the products during the simulations (Tables S1−S10, Supporting Information, give secondary structure details and electrostatic interactions). These residues are derived from contact analysis, taking into account interactions of less than 5 Å between the product and the enzyme. The majority of the nonbonded interactions between the products and the channel-forming residues are hydrophobic aromatic and hydrophobic aliphatic. In contrast to observations on the ligand egress in the bacterial P450cam,[46,47] for the CYP3A4, we observe no salt bridge perturbations during product exit through any of the channels. The two salt bridges located at channels 2a and 2b between D76-R106 and D76-R372 are very stable during the exit of both products. Polar and some hydrophobic residues serve as hydrogen-bond donors and acceptors that H bond usually the hydroxyl and carbonyl groups of the products (Table 2, superscript). The active site key residues that hold the product by hydrogen bonds are S119 and R212. S119 was experimentally shown to alter testosterone hydroxylation selectivity.(48) Mutations of S119 to hydrophobic residues switched the 6β to 2β hydroxylation. Thus, the loss of the S119 directed hydrogen bond to the substrate by site-specific mutations changes the regioselectivity of testosterone oxidation. Product release from the CYP3A4 active site is initiated by hydrogen-bond breakage between the product and R212 and S119 (Tables S7−S10, Supporting Information) and detachment from hydrophobic residues, I301, F304, A305, F213, A370, and L482. Overcoming these barriers leads to product movement toward the gate keeper residues. The products pass the gates of channels 2a, 2b, 2c, 2e, 3, and S after 900, 800, 1150, 700, 900, and 750 ps, respectively.

Table 2

Major Channel Forming Residues of Substrate Channels 2a, 2b, 2c, 2e, 3, and S in the CYP3A4 Enzyme Obtained by Close Contact Analysisa with the TMZ and TST-OH Products

channel	channel residues
2a	L47b	Y53b	F57	D76	R105	R106	F108	S119c
	R212c	F213b,c	F215b	L216	L221	T224	R372	L482b
2b	F57	D76b	Q79b,c	R105	R106	F108b	S119c	R212b
	F213	F215	F220	T224b	P227b	A370b	R372	L482
2c	R105	F108	M114b,c	I118	S119b,c	I120	R212c	F213c
	F215	V240	F241c	T246b	I300c	F304c
2e	R105b,c	F108	G111	K115	I118	S119b,c	I120b	E122b
	R212c	F213	F215	V240	A370
3	F108b	S119b,c	I120	K209c	R212c	F213b	F215	F219
	F220	V240	F241	P242	R243	V245c	F304b
S	S119c	K173c	L211	R212c	F213c	E308b	T309	S311b
	S312b,c	I369b	A370	M371	L482b	L483c	Q484c

Close contacts were defined by product atoms having a distance of less than 5 Å from enzyme residues.

The residue serves as a H-bond acceptor.

The residue serves as a H-bond donor.

Substrate Channels’ Gating Mechanism

The invested total SMD work (WSMD) during each product egress was measured every 1 ps according to eq 1 where WSMD is the mechanical work done along a SMD trajectory, n is the number of MD steps of a trajectory, F⃗(t) is the SMD force invested at time t, D⃗(t) and D⃗(t − 1) are the coordinates of the center of mass of the product at times t and t − 1, respectively, and θ is the angle between the force and the product displacement vectors. The calculated total WSMD values are summarized in Table 3, and the entire SMD work profiles for all simulations are presented in Figure S6 (Supporting Information). Most of these profiles have a similar behavior. At the beginning of the simulations, WSMD increases, reaching a global maximum that represents the potential barrier of product release via the channel gate. Following the product exit, there is usually a fast decline in the work to almost zero. For TMZ (TST-OH), the largest total WSMD observed was while exiting from channel S (2c), and the lowest was from channel 2e (2e and 3). In order to test the robustness of our results, we performed three different SMD runs for each CYP3A4-detected channel with either TMZ or TST-OH. The TMZ (TST-OH) exit from the CYP3A4 active site through channels 2a, 2b, and 3 (2c, 2e, and S) was repeated three times. The results presented similar exit pathways, Cα rmsd and SMD work profiles (Figures S7−S9, Supporting Information). In addition, the calculated total SMD work was similar for each triplicate, with a relative standard deviation of not more than 3% (Table 3).

Table 3

Total SMD Work (Calculated According to Equation 1) Done along the TMZ and TST-OH Trajectories while Exiting through Substrate Channels 2a, 2b, 2c, 2e, 3 and S (Units Are in kcal/mol)

channel	TMZ	TST-OH
2a	265.0 ± 5.5a	283.2
2b	235.1 ± 7.5a	268.7
2c	292.1	445.4 ± 14a
2e	203.6	244.6 ± 6.9a
3	274.8 ± 8.9a	235.1
S	337.1	276.2 ± 5.8a

Average of three runs.

In order to define the gating mechanism of product egress from the substrate channels, we focused on MD snapshots 250 ps prior to and following the time frame where the maximum WSMD was observed. From visual inspection, rmsd, and distance profiles, we identified the gate keepers of each channel for both products. Figure S10 (Supporting Information) presents the distance along MD between the two major gating residues of each channel found for both products. Table 4 summarizes the gating residues of each channel for both products.

Table 4

Gating Residues of Substrate Channels 2a, 2b, 2c, 2e, 3, and S Identified in the SMD Simulations of TMZ and TST-OH

channel	TMZ	TST-OH
2a	F57−F215	F57−F215
2b	F108−F220	F108−F220
2c	F108−F241	F108−F241
2e	F108−I120	R105−S119
3	F213−F241	F213−V240
S	R212−L482	R212−L482

We observe that for channels 2a, 2b, 2c, and 3, the gates consist of hydrophobic residues, mainly phenylalanine residues. These residues parallel- or T-shape π-stack with each other, and the product disrupts this interaction when exiting the enzyme. This disruption is shown in the sequence of frames in Figures 4 and 5. In channel 2e, the release of TST-OH resulted in a breakage of a hydrogen bond between the carbonyl group of S119 and one of the guanidinium moiety hydrogen atoms of R105. However, for TMZ, the exit through channel 2e is coupled to breakage of aliphatic hydrophobic interactions between F108 and I120. This observation is experimentally supported by F108W and I120W mutations(49) that showed higher 1OH over 4OH midazolam (structurally related to TMZ) hydroxylation compared to the wild-type enzyme. Mutations to tryptophan are believed to close channel 2e and impose a switching of the midazolam entrance to another channel that may lead to a different active site orientation and product formation. In channel S, we observe interruption of aliphatic interactions between R212 and L482, while the product escapes the enzyme, as portrayed in Figure 6.

Figure 4

TST-OH release from substrate channel 2a. TST-OH is colored red, heme is in green, and gating residues are in blue. The enzyme is represented by cyan cartoons. The figures have the same orientation. The time frames are after (A) 0, (B) 450, and (C) 715 ps from the start of the simulation.

Figure 5

TMZ release from substrate channel 2b. TMZ is colored red, heme is in green, and gating residues are in blue. The enzyme is represented by cyan cartoons. The figures have the same orientation. The time frames are after (A) 10, (B) 400, and (C) 900 ps from the start of the simulation.

Figure 6

TST-OH release from substrate channel S. TST-OH is colored red, heme is in green, and gating residues are in blue. The enzyme is represented by cyan cartoons. The figures have the same orientation. The time frames are after (A) 0, (B) 460, and (C) 960 ps from the start of the simulation.

We therefore observe three types of channel gating mechanisms, (i) breaking F−F π-stacking; (ii) breaking aliphatic contacts, as seen at the S channel between R212 and L482; and (iii) breaking the H bond, as observed in the exit of TST-OH through channel 2e between R105 and S119. The predominant mechanism is the π-stacking.

Discussion

Dominant Channel in CYP3A4

Of all known P450s, the isoforms CYP2C9, CYP2D6, and CYP3A4 account for the vast majority of drug metabolism.(3) In the CYP2C9 structure bound to flurbiprofen,(50) the solvent channel was observed to be relatively open. It was suggested that for the human CYP2D6, channel S, defined earlier by Haines et al.(51) for the bacterial CYPBM3, is the main substrate access channel.(52) Here, TST-OH leaving the active site through channel S resulted in the largest opening compared to the other channels with a similar total WSMD to most detected channels (Table 3). Hence, the solvent channel in CYP3A4 is likely to serve as a substrate/product exit channel as in the human P450 isoforms.

CYP3A4 is tethered to the ER membrane by its N-terminus and contacts the membrane through other parts of the enzyme, including the F′G′ loop in the vicinity of substrate channels 2a and 2b.(53) Therefore, channels 2a and 2b are predicted to lie close to or dip into the ER membrane. TST-OH contains three hydrophilic positions, and it is less likely that it will exit through those channels.

It was shown by mutation and kinetic studies that cholesterol, the precursor of testosterone (TST), enters the microsomal CYP7A1 through the membrane and interacts with the F′G′ loop residues.(54) Williams et al.(53) suggested that lipophilic substrates and products could enter and exit through the lipid bilayer. As the exact orientation of channels 2a and 2b with respect to the membrane is unknown, the exit of TST-OH through channels 2a and 2b cannot be ruled out.

In terms of the invested SMD work, TST-OH was found to exit through channel 2c with the largest difficulty compared to any other channel and with a thermodynamically favorable exiting through channels 2e and 3 (Table 3). We therefore conclude that channel 2c is less likely to serve as a gateway for testosterone major product in CYP3A4, while channels 2e and 3 are predicted to form a gateway for TST-OH. In addition, channel S was found to be highly flexible. It faces the cytosol and contains a charged gating residue which could lead the product out of the CYP3A4 active site.

In contrast to TST-OH, although the egress of TMZ resulted in a large opening of the S channel, it needed the largest WSMD compared to the rest of the channels (Table 3). This is accounted for by a high potential barrier derived from electrostatic interaction between the positive R212 and the oxygen of the hydroxyl group of TMZ (Figure S11, Supporting Information). We observe the lowest work profiles when TMZ exits from the CYP3A4 active site via channel 2e, suggesting that this channel is preferred for the exit of the hydroxylated product of diazepam.

Thus, there is no one specific access/exit route for substrates/products in CYP3A4.

Dominant Channel in P450s

Schleinkofer et al.(21) carried out MD simulations on CYP2C5, which is closely related to CYP2C9, pulling out of the active site progesterone, 21-hydroxyprogesterone (PROG-21OH), and benzenesulphonamide. The 2c channel was the predominant exit route. Li et al.(20) tested by MD simulations the egress of TST from the CYP2B1 active site. They found that the 2c channel is one of the preferred exit pathways for TST, in contrast to our results of CYP3A4. The differences in these observations may arise from different substrate specificities between P450 subfamilies.

We hypothesize that the product outlet pathway is determined in part by the post-oxidation product orientation inside of the active site. To test this hypothesis, we pulled TST-16OH, positioned in a 16 hydroxylation orientation, out of CYP3A4 from channel 2c. The TST-16OH exit resulted in a much faster and easier exit than that of TST-6OH, with a maximum pulling force of about half of that observed for TST-6OH (Figure S12, Supporting Information). In addition, the mechanical work measured along the TST-16OH escape from channel 2c was 189.4 kcal/mol. This easy escape from channel 2c is in accord with the previous simulations of the exit of TST-16αOH(20) and PROG-21OH(21) from CYP2B1 and CYP2C5, respectively. Both TST-16OH, which is hydroxylated mainly at the TST 16 position(55) in CYP2B1, and PROG-21OH are positioned inside of the active site in the general direction of channels 2c and 2e exits. However, the major hydroxylation product of TST by CYP3A4 is the TST-6βOH, with minor hydroxylation at position 16 (Figure 2),[56,57] and docking it accordingly resulted in a different orientation than that of TST-16OH and in a favorable orientation for exiting via channel 3 (Figure S13, Supporting Information).

We therefore conclude that the oxidized position in large products that have less rotational and translational freedom inside of the active site can determine the outlet pathway. Furthermore, there is no single dominant exit pathway for all P450s, and different products of the same substrate (e.g., TST) leave the active site through channels based on the substrate oxidation specificity.

Conservation of the Gating Mechanism

We find that the gating residues are bulky and hydrophobic in most CYP3A4 substrate channels (Table 4). Schleinkofer et al.(21) pointed to V106 and K241 in CYP2C5 as the gate keepers of channel 2c. In contrast, in the channel 2c of CYP3A4, the gating residues are F108 and F241, despite the fact that the channel faces the hydrated cytosol. Li et al.(20) also observed that channel 2c in CYP2B1 is gated by two phenylalanine residues. They also pointed to a phenylalanine residue as the gate keeper of channel 2e in CYP2B1. Winn et al.(14) showed that CYPcam channel 2a is gated by F87 and F193, similar to our predicted gating mechanism. However, they observed different gating mechanisms for other bacterial P450s that include a positive arginine.

From all of the above and considering that the channels are part of the variable SRSs regions, the gating mechanisms are different for corresponding channels in different P450 isoforms. However, the phenylalanine gating mechanism exists in many P450s, which implies a nonspecific substrate channel.

Conclusions

We performed SMD simulations on CYP3A4 pulling out of its active site two products resulting from the hydroxylation of diazepam and testosterone. We identified several putative access channels connecting the active site to bulk solvent. We observed no single dominant channel in CYP3A4. For both temazepam and testosteron-6βOH, we observed a thermodynamic preference for exiting through channel 2e. In addition, testosteron-6βOH showed a low mechanical work profile while exiting through channel 3, which faces the cytosol as does channel 2e. The solvent channel was found to be the most flexible in CYP3A4 and is predicted to be a substrate/product channel as proposed earlier for other P450s.(52) Moreover, we found that the product shape and orientation after being oxidized determined the preferred channel from which it would exit. In the majority of the channels, the gate keepers are π-stacked phenylalanine residues; to exit, the product has to break this interaction. This gating mechanism is not conserved in P450s but exists in many isoforms and in different channels.

The versatility of the egress mechanism in CYP3A4 makes perfect sense, considering the fact that this isoform is promiscuous and oxidizes a great variety of substrates. This variety requires, in turn, a good number of egress pathways. We believe that these data will assist in identifying some structural determinants related to substrate specificity and in understanding drug metabolism and toxicity in CYP3A4.