Electron Transfer from Cytochrome P450 Reductase to Cytochrome P450: Towards a Structural and Dynamic Understanding

Cytochrome P450 (CYP) enzymes are heme monooxygenases that require two electrons for their catalytic cycle. For the 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). CYP and CPR each attach to the membrane of the endoplasmic reticulum by a single transmembrane helix with their catalytic domains in the cytosol. No structure of a mammalian CYP-CPR complex has been solved experimentally, hindering a detailed understanding of the determinants of electron transfer (ET), which is often rate-limiting for CYP reactions. Here, we investigated the structure and dynamics of the interactions between membrane-bound CYP 1A1, an antitumor drug target, and CPR by a multiresolution computational approach, combining Brownian dynamics, coarse-grained and all-atom molecular dynamics simulations. 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 between the membrane and 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. TOC GRAPHICS


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 regioand 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 (MBD), which serves as a nonspecific membrane anchor, and a flexible tether region that connects the MBD and the FMN domains 6 . The MBD 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 (TM) α-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 TM-helix (Fig. 1). 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 TM 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][25][26][27][28][29] (Fig. 1).
Some molecular simulation studies of CYP-CPR complexes have been carried out on the basis of the crystal structure of Bacillus megaterium cytochrome P450 BM3 (P450BM3), 30 a protein that has the CYP HEME domain and the three reductase domains fused in a single soluble polypeptide chain. Hence, the P450BM3 crystal structure, which contains a non-stoichiometric complex of two CYP HEME domains and an FMN domain, has served as a template for modelling CYP-CPR interactions. 31 We initially took this approach for modeling the CYP 1A1-CPR complex but, despite extensive MD simulation, we did not succeed in obtaining ET-competent complexes.
Indeed, considering the P450BM3 structure as a template for modelling mammalian CYP-CPR interactions has four major drawbacks: (i) The crystal structure of P450BM3 (PDB ID: 1BVY) is not consistent with measured ET rates due to the high redox center separation distance; (ii) The low sequence identity between P450BM3 and CYP/CPR of about 30% for the CYP and FMN domains; (iii) Unlike mammalian CYPs, the dimeric form of P450BM3 is catalytically functional and ET from FMN to HEME is proposed to occur in a trans fashion; 32 (iv) Mammalian CYP and CPR are membrane-anchored proteins, whereas P450BM3 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. Step 1: Brownian dynamics (BD) rigid-body 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: CG and AA 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.
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, see CYP 1A1 is an extrahepatic human drug target protein that plays an important role in the bioactivation of procarcinogens and detoxification of carcinogens. 33,34 Many of its substrates, such as benzo[o]pyrene, a carcinogen found in tobacco smoke, are polycyclic aromatic hydrocarbons.
The standard reference substrate for CYP 1A1 is 7-ethoxyresorufin, which is O-dealkylated by the enzyme. Constitutive expression of CYP 1A1 is low but can be induced in response to environmental contaminants, for example in smokers' epithelial lung cells. Moreover, engineered CYP 1A1 has potential as a biocatalyst. 35,36 Notably, CYP 1A1-CPR fusion proteins have been made to improve ET rates. A fusion between rat CYP 1A1 and rat CPR was found to exhibit four times higher monooxygenase activity than rat CYP 1A1 alone, 35 whereas a fusion between rat CYP 1A1 and yeast CPR had a similar ET rate to P450BM3. 36 Hence, a detailed atomic level understanding of CYP 1A1-CPR interactions and ET kinetics would provide the basis for the design of drugs targeting CYP 1A1, as well as the exploitation of CYP 1A1 as a biocatalyst. We therefore here describe the building and simulation of the full-length CYP 1A1-CPR complex in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer. POPC was chosen because phosphatidylcholine is the main component in the mammalian ER membrane 37 and because it provides a simple representation, which is often used in in vitro studies of CYPs, of the disordered regions of the ER where CYP 1A1 has been observed to localize. 12,38

Electrostatic steering of the CYP 1A1 and the CPR FMN domains leads to encounter complex formation
The BD docking method allows generation of diffusional encounter complexes that satisfy biochemical constraints. It has been successfully applied to predict protein-protein complexes, 39 including those between cytochrome P450cam and putidaredoxin. 40 The six representative encounter complexes of the globular domains of CYP 1A1 and CPR generated by BD rigid-body docking simulations and selected for subsequent MD simulations are listed in Table 1   , cofactors are shown in pink stick representation, and the bilayer is shown in cyan lines with orange spheres representing the phosphorous atoms. The angles defining the arrangement of the proteins are shown: θ is the angle between the CYP C-helix (green) and the FMN domain α1-helix (residues 91-105; red); the FMN domain α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 zaxis 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).
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). The FMN domain binding mode differs from that in the crystal structure of the P450BM3 CYP-FMN domain complex. 30 In the encounter complexes, both the α1and α3-helices of the FMN domain contribute to the interface whereas, in P450BM3, it is mainly the α1-helix of the FMN domain. The distance between the redox centers, DFe-N5, in the six representative encounter complexes ranges from 15.4-20.2 Å and is shorter than the distance (23.6 Å) observed in the crystal structure of P450BM3 even though the distance between the centers of mass of the globular domain of CYP 1A1 and the FMN domain, DCYP-FMN domain, is higher in the representative diffusional encounter complexes (34.6-37.9 Å) than in the P450BM3 crystal structure (31.7 Å). ( Table 1). The relative orientation of the FMN domain and the CYP globular domain can be characterized 41 by the angle θ between the α1-helix of the FMN domain and the CYP Chelix (see Fig. 2a). Angle θ ranges from 54° to 99° in these 6 complexes, whereas in the P450BM3 crystal structure, it is 95°. Nevertheless, the binding face of CYP 1A1 in all these complexes is similar to that identified in mutagenesis studies of CYP 1A2-CPR 42 and CYP 2B4-CPR binding. 27 Next, we relaxed the BD rigid-body docked complexes by performing MD simulations.

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 hydrogenbonding of the FMN cofactor (A4, see Fig. S3).
For the remaining three encounter complexes (C2, C3 and D2), decreases in DCYP-FMN domain and DFe-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).

(a) Weakening of CYP 1A1-membrane interactions and reorientation of CYP globular domain
During the simulations, the axial distance between the centers of mass of the globular domain of CYP 1A1 and the membrane, DCYP-mem, increased by 3-8 Å (Table S4) (Table S4), indicating a weakening of the peripheral CYP interactions with the membrane.
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 hemetilt 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). Thus, the interactions with CPR, resulting in an increase of the area of the FMN domain interface of up to ~700 Å 2 , weakened the CYP-membrane interactions, and facilitated reorganization of the CYP with respect to the membrane.  Table S4). This domain motion, as well as the initial binding mode of the encounter complexes, brought the FAD and NADP domains nearer to the globular domain of CYP 1A1 in the C2 and D2 systems (Fig. 4). For the D2' simulation, the decrease in the FAD-NADP interdomain distance resulted in a more compact CPR structure, similar to that of the semi-open conformation in the crystal structure (PDB ID 3ES9 chain A) (Fig. S5). In this structure, the center-to-center distance between the FMN and FAD cofactors was reduced to 44.

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 P450BM3 (31.7 Å), possibly because the membrane prevents the two domains from coming as close as in the membrane-free P450BM3.
No contacts between the TM-helices of CYP 1A1 and CPR were observed, which is consistent with the observations of Sundermann and Oostenbrink 32 in their 10 ns simulation of a membrane-

CYP-CPR complexes but shows differences to the P450BM3 CYP-FMN domain interface
The CYP-FMN domain interactions involve charge pairing of the positively charged proximal side of CYP 1A1 with the negatively charged surface of the FMN domain surrounding the cofactor.
These electrostatic interactions are complemented by van der Waals and hydrophobic interactions.
The interface residues in CYP 1A1 are located on the B, C and L-helices and the JK-loop and the loop structure near the HEME (Table S2 and Fig. 5e). The interface residues of the CPR FMN domain surround the cofactor and involve the α1 and α3-helices, Lα′ and Lβ1-5 loops and β1 and β2-strands (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,45 CYP 17A1 (cyt b5), 46,47 CYP 1A2 (CPR) 42 and P450BM3 (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) along with proximal residues close to those identified in our CYP 1A1 complexes (Fig. S6).   41 found that after MD simulation, the α1-helix reoriented and E494 lost its hydrogen-bond to H100 and approached K97 of P450BM3, 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 Cand α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: P450BM3 FMN and CYP domains: white; CPR FMN and CYP 1A1 globular domains : salmon and light blue, respectively; CYP C-helices: green; FMN domain α1 helices of P450BM3 and CPR: pink and red, respectively; FMN cofactors: pink, with P450BM3 and CPR cofactors in licorice and ball-andstick representation, respectively; HEME: magenta carbons; 7-ethoxyresorufin: orange; selected residues are shown in yellow carbon representation with hydrogen bonds indicated by black dashed lines.

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. 5c-d) occupancy) with the J-helix and the JK-loop of CYP 1A1 (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. 51  systems.

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, 53 whereas for eukaryotic CYPs, different rate-determining steps -formation of the membrane-bound CYP-CPR complex 54 , product release 36 , or the first 55 or second ET 56 -have been reported. The first ET to mammalian CYPs from CPR occurs at ~2-10 s −1 57,58 and the second ET is usually within this range or slower. 42,57,59 Irrespective of the different modes of interaction between CPR and CYP in the three plausible complexes (C3, D2 and D2'), only two major ET pathways were identified. These are FMN-I458-C457-HEME (56%) and FMN-K456-C457-HEME (35%), see Fig. 5g-h. Only the latter pathway was present in the D2' simulation, whereas the former pathway was the predominant ET route in the other two simulations (C3: 77% and D2: 96%). Evidence supporting both these paths, which involve the backbone of L-helix residues and the heme-coordinating cysteine, comes from experiments and simulations.
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 Km value, indicating slightly tighter binding of the substrate, possibly due to realignment of the heme in the mutant. 60 The I458P mutation had a smaller and opposite effect to the I458V mutation. 60  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), DFe-N5 and the θ angle ( Fig. S8a-b) were plotted against log(kET). The maximum population of structures has log(kET) between 0.1 to 0.7. The calculated ET rates ( Table   1) agree well with the measured value 42 of kET 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. 64 The histogram plot of DFe-N5 vs log(kET) (Fig S8a) shows that, consistent with the distances observed experimentally in several other redox proteins, 65 DFe-N5 ranges from about 13 to 17.5 Å, and that the experimental kET value is reproduced with this range of DFe-N5. Notably, the histogram plot (Fig. S8b) of θ vs log(kET) 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) tunnel was present and, in the D2 simulation, only tunnel 2ac was present. These tunnel closures 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 I-helix 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. (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.
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. 67 Our results indicate, however, that more extensive studies of substrate access and product release would be necessary to confirm the effect of CPR.  The importance of identifying ET pathways is highlighted by the fact that mutation of I458, a key residue in the dominant computed ET pathway, has been shown to affect prodrug activation. 60 The most common non-synonymous polymorphisms of CYP 1A1 have the mutations T461N and/or I462V, resulting in high-risk variants for non-small cell lung cancer in non-smokers. 68

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, 71 following refinement and rebuilding using the PDB-REDO server 72 78 Missing residues were modelled using VMD 79 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,80-82 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. [80][81][82] 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 Å. 65 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, DFe-N5 < 20 Å, and (ii) the center-to-center distance between the CYP globular domain and the CPR FMN domain, DCYP-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 account for protein flexibility, all-atom MD simulations in 150 mM aqueous solution (referred to as "soluble simulations") were carried out for all 6 representative structures of CYP 1A1 ̶ CPR encounter complexes obtained from BD simulations. For computational efficiency, only the globular domain of CYP 1A1 and the FMN domain of CPR were simulated in the soluble simulations. After equilibration, all-atom MD production simulations were performed for about 50 to 140 ns using the NAMD v2.10 software. 83

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 Hence, we superimposed the CYP domain of the CYPglobular-domain ̶ CPRFMN-domain complexes onto a modelled structure of the membrane-bound full-length CYP 1A1 in a POPC membrane. 77

(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.

Post-analysis of MD trajectories
CYP globular domain-FMN domain interaction free energies were computed using the MMGBSA 86,87 protocol implemented in the AMBER 14 package for snapshots collected at 100 ps intervals after completion of 350 ns simulation (i.e. for the last 150-200 ns of the trajectories).
Substrate access tunnels were computed using CAVER 3.0 88 for the last 500 snapshots of each trajectory collected at 100 ps intervals. Random acceleration molecular dynamics (RAMD) 89 simulations, as implemented in NAMD 2.10, were performed to study the ligand egress pathways from the CYP 1A1 active site in the presence and absence of CPR. ET pathways and rates were computed using Beratan's model with the VMD pathway plugin module. 90

AUTHOR CONTRIBUTIONS
GM, PPN and RCW designed the research; GM performed computations and GM, PPN and RCW analyzed the computational data; GM and RCW wrote the manuscript.