An electron transfer competent structural ensemble of membrane-bound cytochrome P450 1A1 and cytochrome P450 oxidoreductase

Cytochrome P450 (CYP) heme monooxygenases require two electrons for their catalytic cycle. For mammalian microsomal CYPs, key enzymes for xenobiotic metabolism and steroidogenesis and important drug targets and biocatalysts, the electrons are transferred by NADPH-cytochrome P450 oxidoreductase (CPR). No structure of a mammalian CYP–CPR complex has been solved experimentally, hindering understanding of the determinants of electron transfer (ET), which is often rate-limiting for CYP reactions. Here, we investigated the interactions between membrane-bound CYP 1A1, an antitumor drug target, and CPR by a multiresolution computational approach. We find that upon binding to CPR, the CYP 1A1 catalytic domain becomes less embedded in the membrane and reorients, indicating that CPR may affect ligand passage to the CYP active site. Despite the constraints imposed by membrane binding, we identify several arrangements of CPR around CYP 1A1 that are compatible with ET. In the complexes, the interactions of the CPR FMN domain with the proximal side of CYP 1A1 are supplemented by more transient interactions of the CPR NADP domain with the distal side of CYP 1A1. Computed ET rates and pathways agree well with available experimental data and suggest why the CYP–CPR ET rates are low compared to those of soluble bacterial CYPs.


Introduction
CYP enzymes constitute the cardinal xenobiotic-metabolizing enzyme superfamily and they are of considerable interest for the synthesis of novel drugs and drug metabolites, targeted cancer gene therapy, biosensor design and bioremediation. [1][2][3] CYPs are widely used as biocatalysts for regio-and enantioselective C-H hydroxylation reactions. CYPs are heme monooxygenases that activate one molecule of dioxygen and oxidize the substrate by inserting one oxygen atom to form the product and reducing the other oxygen atom to water. Thus, drugs are transformed into more polar metabolites during the process of drug metabolism. 4 In the catalytic cycle, CYP requires two electrons to be transferred sequentially from a redox partner to its heme cofactor (HEME), and the ET steps have been observed to be rate-limiting for the reaction in many cases. 5 For the microsomal CYPs, the redox partner is CPR, a flavoprotein containing three cofactor binding domains: the FMN domain with the flavin mononucleotide (FMN) cofactor, the FAD domain with the flavin adenine dinucleotide (FAD) cofactor, and the NADP domain with the nicotinamide adenine dinucleotide phosphate (NADP) cofactor. Besides these three domains, CPR also contains an amino-terminal membrane binding domain, which serves as a nonspecific membrane anchor, and a flexible tether region that connects this domain and the FMN domains 6 . The membrane binding domain is necessary for efficient ET from CPR to CYP 6 and its interactions with the membrane have been found to depend on cofactor binding and the oxidation state of CPR. 7 Cytochrome b5 (cyt b5) is another natural redox partner for some CYPs and can also allosterically modulate CYP catalysis. [8][9][10][11] All mammalian CYPs have three domains: a globular catalytic HEME-containing domain, an Nterminal transmembrane α -helical domain, and a flexible linker region connecting these two domains. In humans, both the CYP and the CPR are anchored in the endoplasmic reticulum (ER) membrane by a transmembrane helix (TM-helix). The N-terminal region of the CYPs not only anchors the protein in the ER membrane, but also influences their localization to different microdomains of the membrane. 12 Information about the orientation of the globular domain of CYPs, alone or in CYP-CPR/cyt b5 complexes in the membrane, has been gleaned from experiments such as linear dichroism measurements in Nanodiscs, 13,14 atomic force microscopy and antibody mapping. 15 To date, there is no experimentally determined structure of a full CYP-CPR/cyt b5 complex available, although structures of the individual globular domains of a number of mammalian CYPs and CPRs have been determined by crystallography after removal of the N-terminal transmembrane domain to facilitate expression and crystallization. 2 Various modeling and simulation procedures have been developed to predict the orientation of full-length CYPs in phospholipid bilayers. [16][17][18][19][20][21][22] Notably, T.
brucei CYP 51 simulated 20 prior to the determination in 2014 of the crystal structure of yeast CYP 51, 23 the only eukaryotic CYP for which a full-length structure has been determined experimentally, indicated a very similar orientation of the protein with respect to the membrane.
Experiments, such as solid-state NMR and mutagenesis, indicate that the association between a CYP and its CPR/cyt b5 redox partner is driven by electrostatic and hydrophobic interactions with positively charged residues on the proximal side of CYP interacting with negatively charged residues on the CPR FMN domain or cyt b5 24-29 .
We initially took the crystal structure of Bacillus megaterium cytochrome P450 BM3 (P450 BM3 ) 30 , a non-stoichiometric complex of two CYP HEME domains and an FMN domain, as a template, for modeling the CYP 1A1-CPR complex but, despite extensive MD simulation, we did not succeed in obtaining ET-competent complexes. Indeed, considering the P450 BM3 structure as a template for modelling mammalian CYP-CPR interactions has four major drawbacks: (i) The crystal structure of P450 BM3 (PDB ID: 1BVY) is not consistent with measured ET rates due to the high redox center separation distance; (ii) The low sequence identity between P450 BM3 and CYP/CPR of about 30% for the CYP and FMN domains; (iii) Unlike mammalian CYPs, the dimeric form of P450 BM3 is catalytically functional and ET from FMN to HEME is proposed to occur in a trans fashion; 31,32 (iv) Mammalian CYP and CPR are membrane-anchored proteins, whereas P450 BM3 is not and the membrane attachment may affect the binding of CPR to the CYP. Hence, it is preferable to employ a de novo protocol rather than a template-based method for modelling full-length CYP-CPR interactions in a membrane.
Here, we have carried out multiresolution simulations with the aim of achieving a detailed understanding of the structure and dynamics of the membrane-associated full-length CYP 1A1-CPR complex, and consequently of the determinants of the ET process between CYP 1A1 and CPR. CYP 1A1 is an extrahepatic human drug target protein that plays an important role in the bioactivation of procarcinogens and detoxification of carcinogens. 33

Electrostatic steering of the CYP 1A1 and the CPR FMN domains leads to encounter complex formation
The multiresolution computational method used to model and simulate the membrane-associated full-length CYP 1A1-CPR complex is depicted in Fig. 1 Table   S1). The representative encounter complexes and their expected position with respect to the membrane are shown in Fig. S2 and one of these complexes is shown in Fig 2a. The interfacial residues in these complexes are listed in Table S2.
In the encounter complexes, the globular domains of CYP 1A1 and CPR dock specifically due to complementary electrostatic interactions, and the proximal face of CYP 1A1 and the outer FMNbinding loop interact (Fig. 1)

Structural relaxation of encounter complexes in "soluble" MD simulations resulted in three putative ET-competent complexes
All "soluble" MD simulations were performed for the globular domain of apo-CYP 1A1 and the  Table S3). The arrangement of the protein domains in the complexes relaxed during the simulations (Fig. 3). Three of the six encounter complexes were discarded after structural relaxation, either because they were no longer ET-competent (B7 and D3, see Table 1) or because of rearrangements in the hydrogen-bonding of the FMN cofactor (A4, see

Fig. S3).
For the remaining three encounter complexes (C2, C3 and D2), decreases in D CYP-FMN domain and D Fe-N5 were observed during the simulations ( Table 1), indicating the formation of tighter complexes with non-zero ET rates. The interface residues of CYP 1A1 remained similar to those in the initial BD docked complexes (Table S2). Next, these three refined complexes were used to build and simulate complete CYP-CPR complexes in the presence of a phospholipid bilayer.

CYP-membrane interactions weaken in the presence of CPR and the CPR structure adopts a more compact form
All three models of membrane-bound CYP-CPR complexes were simulated for about 500 ns. In addition, a second simulation (D2'), with identical initial coordinates but with different initial assigned velocities, was run for the D2 system. The structures of the individual protein domains were well maintained as shown in C α -RMSD plots ( Fig. S4) but there were rearrangements of the domains within CPR and of the proteins with respect to the membrane (Table S4).   In the presence of CPR, the CYP 1A1 globular domain underwent an orientational rearrangement with respect to the CPR-free membrane-bound state as indicated by the changes in α , β and heme-tilt angles 16 (Fig. 2b) along the simulations (Fig. 4, Table S4). Despite the almost identical starting orientations of CYP 1A1 in the membrane, the simulations converge to quite different orientations of the CYP globular domain due to the different binding modes of CPR to CYP 1A1. In the presence of CPR, the heme tilt angle dropped from ~71° to around 40° for the C2 and D2 systems (Table S4)

CYP 1A1 and CPR form stable complexes in the membrane
In all the simulations, the interdomain distance D CYP-FMN domain was stable at 33-36 Å. It remained longer than the interdomain distance in the crystal structure of P450 BM3 (31.7 Å), possibly because the membrane prevents the two domains from coming as close as in the membrane-free P450 BM3 . Nevertheless, in all the CYP 1A1-CPR complexes, D Fe-N5 was within the range (14-16 Å) for biological ET (Table S4) and converged stably in all the simulations ( Fig. 4).
No contacts between the TM-helices of CYP 1A1 and CPR were observed, which is consistent with the observations of Sundermann and Oostenbrink 46 in their 10 ns simulation of a membrane-bound CYP 2D6-CPR complex. However, the computed interface contact area of the globular domain of CYP 1A1 and the CPR FMN domain in the model complexes was ~2000 Å 2 (Table S4), much higher than the interface area of about 240 Å 2 observed in the simulation of a membrane-bound CYP 2D6-CPR complex, 46 and greater that the interface area of many transient protein-protein complexes, which is usually is less than 1500 Å 2 . 47,48 This large interface area is consistent with the measured binding affinity of CPR to rat CYP 1A2 (which has 68% identity to human CYP 1A1) of 47 nM. 45 The CPR NADP domain interactions with the distal side of CYP create an additional contact interface. The high contact area of the CYP 1A1-CPR complexes indicates the formation of a rather strongly bound transient complex between the proteins.

The interface between CYP 1A1 and the CPR FMN domain mostly has similar residues to other CYP-CPR complexes but shows differences to the P450 BM3 CYP-FMN domain interface
The  (Table S2 and Fig. 5f). In all cases, a hydrogen-bonding interaction between the FMN phosphate group and the side-chain of Q139 in the C-helix of CYP 1A1 was also present. NMR spectra of a membrane-anchored CYP 2B4-CPR FMN domain complex 24 showed the same region of the FMN domain at the interface.
Site-directed mutagenesis and structural studies have been performed to identify CYP residues that interact with the redox protein for CYP 2B4 (CPR), 27,49 CYP 17A1 (cyt b5), 50,51 CYP 1A2 (CPR) 45 and P450 BM3 (FMN domain), 30 revealing the importance of positively charged residues for redox protein binding. Even though the residues at the FMN domain binding interface of CYP 1A1 do not align completely with the interfacial residues of the other CYPs in a multiple sequence alignment (Fig. S6), the secondary structure regions at the interface are similar, mostly involving the B, C and L-helices, and the CC'-loop, HI-loop and the loop near the HEME. The alignment revealed that the positively charged CYP residues at the FMN domain interface were most conserved in the C3 simulation, with the conservation being highest for CYP 1A2 (80%).
However, there are some differences, e.g. V267 and L270 in the H-helix of CYP 2B4 mediate CPR interactions 27 whereas neither the corresponding aligned residues in CYP 1A1, S284 and E287, nor the H-helix make any contacts with CPR in our simulations. An NMR study of the binding of the CPR FMN domain to CYP 17A1 51 indicated that distal residues were affected by FMN binding along with proximal residues close to those identified in our CYP 1A1 complexes ( Fig. S6). In CYP 19A1, the B-helix (K108), C-helix (S153 and G154), K-helix (K352), K''-Lloop (K420, Y424 and R425) and L-helix (Y441) were reported to be involved in interactions with FMN domain of CPR in a modelled CYP 19A1-CPR structure. 52 In a model of the membrane-bound CYP2D6-CPR complex, 46 the B-helix (R88), C-helix (R129, R133 and R140), CC'-loop (K146, K147 and S148) and loop near heme moiety (K429, E431 and R440) were found to interact with the FMN domain of CPR. In our simulations, we found that a similar, positively charged region of CYP 1A1 was involved interactions with the FMN domain of CPR (Table S2).

Transient interactions are formed between CYP 1A1 and the NADP domain of CPR
In the initial model for the D2 system, the NADP domain of CPR was relatively far from the membrane and from the globular domain of CYP 1A1 (Fig. 2a). During the simulations, the NADP domain came closer to the CYP 1A1 globular domain and formed interactions (Fig. 5cd (Fig. 5e). In the D2' simulation, the NADP domain moved towards the CYP globular domain after 200 ns and residues 664, 667 and 668 formed contacts (with ≥ 50% occupancy) with the EF-loop, the C'helix, the G-helix and the GH-loop (Fig. 5e) During the simulations of the C2 complex, the FAD and NADP domains of CPR moved towards the phospholipid bilayer and started to interact with the head groups (Fig. 5a). Otyepka and colleagues made a similar observation in their simulations of a full-length CPR-membrane system. 56

Simulations reveal two ET pathways and ET rates consistent with experiment
CPR transfers two electrons sequentially to CYPs during the catalytic monooxygenation cycle.
The second ET step is rate-determining for most bacterial CYPs, 58  As regards the first path, experiments showed that the I458V point mutation in CYP 1A1 enhanced N-demethylase activity by about two-fold due to a decreased K m value, indicating slightly tighter binding of the substrate, possibly due to realignment of the heme in the mutant. 65 The I458P mutation had a smaller and opposite effect to the I458V mutation. 65  To investigate the relation between ET rate and the structure of the complex, 2D histograms of the computed binding free energies between the CYP globular domain and the FMN domain in the complexes (Fig. 5i), D Fe-N5 and the θ angle ( Fig. S8a-b) were plotted against log(k ET ). The maximum population of structures has log(k ET ) between 0.1 to 0.7. The calculated ET rates ( Table 1) agree well with the measured value 45 of k ET for the rat CYP 1A2-CPR complex expressed in yeast of 5.9 s -1 . Fig. 5i shows that the highly populated clusters in the trajectories have intermediate binding free energies, and this is consistent with evidence from surface plasmon resonance measurements that the association of CYP and CPR in ET-compatible orientations is governed by entropic rather than enthalpic contributions. 70 The histogram plot of D Fe-N5 vs log(k ET ) (Fig S8a) shows that, consistent with the distances observed experimentally in several other redox proteins, 71 D Fe-N5 ranges from about 13 to 17.5 Å, and that the experimental k ET value is reproduced with this range of D Fe-N5 . Notably, the histogram plot (Fig. S8b) of θ vs log(k ET ) has a region from 82° to 90° where no structures were sampled during the simulations as this arrangement of the CYP 1A1 C-helix and FMN domain α 1-helix is disfavored as discussed above.

CPR affects CYP ligand tunnels due to CYP reorientation and CYP-CPR interactions
We analyzed the effects of CPR binding on the opening of ligand tunnels between the buried CYP active site and the protein surface in the conventional MD simulations described above (Fig. 7) and in RAMD simulations of ligand egress from the CYP active site (Table 2).
In the conventional MD simulations, two ligand tunnels in CYP 1A1 with diameters sufficient for passage of a water molecule were identified (Fig 7a-d) were due to BC-loop motion resulting in hydrogen-bond formation closing tunnel 2ac (Fig. 7e) and movement of the side-chain of F224 away from the ligand to interact with F319 in the Ihelix and close the solvent tunnel (Fig. 7f). A direct influence of CPR on these motions could not be identified by contact map analysis of the trajectories.
The RAMD simulations, in which 7-ethoxyresorufin travelled from the active site to the protein surface, revealed five ligand egress routes in the CYP-CPR complexes (S, 2ac, 2c, 3s and 3) and three routes (S, 2c and 3) in CPR-free CYP 1A1 ( Table 2). Although the transient interactions between F224 and F319 were perturbed in the RAMD ligand egress simulations and thus egress via the solvent tunnel was observed in all systems, the solvent tunnel tended to close and either the 2c or the 3s tunnel tended to open in the CPR-bound systems. Tunnel 2ac was rarely observed as an egress route, even though it was detected in the conventional MD simulation.
Egress via tunnel 3s, between the F-and G-helices, was observed in two of the CPR-bound systems and occurred when the contacts of N221 with A250 and D253 (occupancy >90% in the CPR-free system) were broken, allowing the F-and G-helices to move apart and permit ligand passage. In summary, in these simulations, the binding of CPR to CYP 1A1 alters the distribution of egress pathways. An alteration in ligand tunnels due to CPR binding has also recently been observed in conventional MD simulations of the CYP 19A1-CPR complex in a membrane. 52 Our results indicate, however, that more extensive studies of substrate access and product release would be necessary to confirm the effect of CPR.    75 . Thus, the present simulations give detailed atomistic insights into structural and functional aspects of CYP 1A1-mediated drug and carcinogen metabolism. Moreover, they provide insights into how CYP-redox partner interactions can vary, and provide a basis for further work to explore the interactions of different CYPs with CPR.

Methods
The computational workflow for building systems for simulations of full-length CYP-CPR complexes in a membrane bilayer is illustrated in Fig. 1 and described briefly below; further details are given in Supplementary Methods and Fig. S1.

Protein structures
Modeling of CYP 1A1 was based on the crystal structure of the globular domain of human CYP 1A1 (PDB ID: 4I8V; 2.6 Å resolution), in complex with the inhibitor α-naphthoflavone, 76 following refinement and rebuilding using the PDB-REDO server 77 1.75 Å resolution). 83 Missing residues were modelled using VMD 84 and Modeller.

Brownian dynamics simulations
BD rigid-body docking simulations of the four structures of the globular domain of CYP 1A1 and the N-terminally truncated CPR were carried out using the SDA 7 software. 39,85-87 All cofactors were retained but the ligand was removed from the CYP. Docking was performed subject to electrostatic interaction, electrostatic desolvation, and non-polar desolvation forces. [85][86][87] 240,000 independent trajectories were run starting with CPR at a random relative orientation and position at a center-to-center distance of 450 Å from the CYP, and they were all terminated when the inter-protein center-to-center separation reached 600 Å. ET in redox proteins typically occurs when the edges of the redox centers lie within 14-15 Å. 71 Hence, the coordinates of diffusional encounter complexes were recorded subject to two distance constraints: (i) the cofactor-cofactor distance from FMN:N to HEME:Fe, D Fe-N5 < 20 Å, and (ii) the center-to-center distance between the CYP globular domain and the CPR FMN domain, D CYP-FMN domain < 60 Å, to ensure close approach of the domains.
For each of the four CYP structures, the energetically ranked top 5,000 docked encounter complexes were clustered into 10 clusters ranked by cluster size using a hierarchical method. 39 Thus, 40 representative encounter complexes were generated and named according to the chain identifier in the crystal structure of CYP 1A1 and the docking cluster number of CPR. E.g., C3 for chain C of CYP 1A1 and the third largest docked cluster of CPR. These complexes were further filtered according to interaction energy, and to ensure that the clusters were structurally distinct and did not have orientations of CPR inconsistent with membrane anchoring if the CYP was assumed to be positioned in the membrane as observed in MD simulations of CYP 1A1 in a POPC bilayer. As a result, six structures of diffusional encounter complexes were selected for further studies.

All-atom molecular dynamics simulations in aqueous solution
To

All-atom molecular dynamics simulations of membrane-bound full-length CYP-CPR systems
All the simulations with phospholipid bilayers, "membrane simulations", were carried out for complete structures of full-length CYP and CPR. The initial configurations of the complexes of the globular domain of CYP and the FMN domain of CPR were obtained from the last frames of the soluble simulations.

(a) Preparation of full-length CYP structures
We previously developed a protocol for inserting the full-length structure of a CYP into a membrane bilayer without predetermining the orientation or insertion depth. 16

(c) Simulation setup and trajectory generation
Each CYP-CPR complex in a bilayer of POPC phospholipids was immersed in a periodic box of water molecules with 150 mM ionic strength maintained by adding Na + and Cl − ions. Production simulations were performed for protein-phospholipid bilayer systems based on three encounter complexes (with one system simulated in duplicate) using the NAMD v2.10 software.   The formation of the CYP-CPR complex is necessary for the transfer of electrons to the CYP active site in the CYP catalytic cycle, as indicated in the schematic cycle.

Post-analysis of MD trajectories
Step 1: Brownian dynamics (BD) rigidbody docking of CYP and CPR globular domains. Molecular electrostatic isopotential contours at ±1 kT/e show a highly positive (blue) patch on the proximal face of CYP and a highly negative (red) patch on the CPR that interact complementarily in the docked complexes.
Step 2: Coarse-grained (CG) and all-atom molecular dynamics (MD) simulation of CYP (blue cartoon representation) in a phospholipid bilayer (cyan with orange spheres representing phosphorous atoms).
Step 3: Relaxation of the BD docked complexes by MD simulation in aqueous solution followed by superimposition on the CYP in the bilayer (with CPR shown in cartoon representation colored by domain (FMN: red, FAD: green, NAD: pink)).
Step 4: Atomic detail MD simulation of the CYP-CPR complexes in a phospholipid bilayer.   3-helix (150-158; cyan) and N-terminal residue of FMN domain (blue sphere) are indicated; the heme tilt angle is the angle between the heme plane and the z-axis perpendicular to the membrane plane; α , β and γ are the angles between the z-axis and the vectors v1, along the I-helix (cyan), v2, orthogonal to v1 and connecting the C-helix (red) and the F-helix (magenta); and v3 along the TM-helix (residues 7-26; green).     1-helix and the CYP C-helix. In contrast, in P450 BM3 , Dubey et al. 44 found that after MD simulation, the α 1-helix reoriented and E494 lost its hydrogen-bond to H100 and approached K97 of P450 BM3 , located at the beginning of the C-helix, close enough to form a salt-bridge. This reorganization resulted in a change from a perpendicular to a parallel arrangement of the C-and α 1-helices. CYP 1A1, however, lacks a corresponding positively charged residue at the N-terminus of the C-helix: P129, A132 and R135 of the C-helix of CYP 1A1 structurally align with K94, K97 and H100 of P450BM3, respectively. Color scheme: P450 BM3 FMN and CYP domains: white; CPR FMN and CYP 1A1 globular domains : salmon and light blue, respectively; CYP Chelices: green; FMN domain α 1 helices of P450 BM3 and CPR: pink and red, respectively; FMN cofactors: pink, with P450 BM3 and CPR cofactors in licorice and ball-and-stick representation, respectively; HEME: magenta carbons; 7-ethoxyresorufin: orange; selected residues are shown in yellow carbon representation with hydrogen bonds indicated by black dashed lines. ) (e-f) Changes in the active site of CYP 1A1 during simulations of the CYP-CPR complexes in the membrane bilayer resulting in closure of ligand tunnels. Structures are shown before (cyan) and after (green) simulation. (e) Formation of a hydrogen bond with 95% occupancy between the side chains of N117 and E256 in the C3 simulation blocking tunnel 2ac. The occupancy of this interaction during the other simulations was less than 30%. (f) Movement of F224 and F319 in the I-helix in the D2 simulation to make van der Waals contact with an occupancy of 56%, thus blocking the route for ligand passage through the solvent tunnel. In the other simulations, the interaction between F224 and F319 was absent. Instead, there were parallel π -π interactions between the ligand and F224.