A multiscale approach to modelling drug metabolism by membrane-bound cytochrome P450 enzymes

Abstract

Cytochrome P450 enzymes are found in all life forms. P450s play an important role in drug metabolism, and have potential uses as biocatalysts. Human P450s are membrane-bound proteins. However, the interactions between P450s and their membrane environment are not well-understood. To date, all P450 crystal structures have been obtained from engineered proteins, from which the transmembrane helix was absent. A significant number of computational studies have been performed on P450s, but the majority of these have been performed on the solubilised forms of P450s. Here we present a multiscale approach for modelling P450s, spanning from coarse-grained and atomistic molecular dynamics simulations to reaction modelling using hybrid quantum mechanics/molecular mechanics (QM/MM) methods. To our knowledge, this is the first application of such an integrated multiscale approach to modelling of a membrane-bound enzyme. We have applied this protocol to a key human P450 involved in drug metabolism: CYP3A4. A biologically realistic model of CYP3A4, complete with its transmembrane helix and a membrane, has been constructed and characterised. The dynamics of this complex have been studied, and the oxidation of the anticoagulant R-warfarin has been modelled in the active site. Calculations have also been performed on the soluble form of the enzyme in aqueous solution. Important differences are observed between the membrane and solution systems, most notably for the gating residues and channels that control access to the active site. The protocol that we describe here is applicable to other membrane-bound enzymes.

Figure 1. Chemical structure of R-warfarin.

In CYP3A4, R-warfarin undergoes hydroxylation at C10.

Figure 2. Membrane deformation about CYP3A4.

The positions of the phosphate particles of all lipids relative to the position of the protein (shown as a grey backbone trace) averaged over a 1 μs coarse-grained molecular dynamics simulation are displayed as a surface. The surface is coloured according to the position along the normal to the bilayer plane (Z-position) with red corresponding to large distances from the bilayer centre and blue corresponding to the smallest distances from the bilayer centre. Thus, red areas correspond to regions of the bilayer that thicken to accommodate the presence of the protein, whilst blue regions correspond to regions of bilayer thinning. A: Side view of the protein within the bilayer. B: Views from above and below the upper leaflet, highlighting the asymmetric deformation of the bilayer. The membrane is observed to deform asymmetrically, with the membrane thinning (blue) in the region of the A-anchor and thickening (red) in the region of the F′- and G′-helices. The lower leaflet is not shown and the protein is transparent for clarity.

Figure 3. Positioning of CYP3A4 in the lipid bilayer from atomistic MD simulations.

This configuration was generated from unbiased coarse grained simulations, in which a mixed POPC/POPE bilayer was self-assembled around the protein. The transmembrane helix spans the bilayer, with the A-anchor and F′ and G′ helices buried in the hydrophobic region of the membrane. This snapshot corresponds to the position of the protein with R-warfarin bound following conversion to full atomistic resolution, equilibration and simulation. Inset: active site region containing Compound I (Cpd I) and R-warfarin (WARF). The protein is shown in cartoon representation coloured by helix, lipids are shown in transparent grey stick representation, with the phosphate particles shown as orange transparent spheres. R-warfarin and Compound I are both shown in stick representation.

Figure 4. Opening and closing of substrate ingress channels over time for membrane simulations of CYP3A4.

Distance between gating residues (d) [in Å] for simulations of membrane-bound CYP3A4 in the absence (A) and presence (B) of R-warfarin. The channels and gating residues are as described in Table 1. Distances were calculated at 0.02 ns intervals over the final 20 ns of each simulation. Similar patterns of opening/closing of gates are observed between the apo model (left) and the R-warfarin bound model (right).

Figure 5. Opening and closing of substrate ingress channels over time for soluble simulations of CYP3A4.

Distance between gating residues (d) [in Å] for simulations of soluble CYP3A4 in the absence (A) and presence (B) of R-warfarin. The channels and gating residues are as described in Table 1. Distances were calculated at 0.02 ns intervals over the final 20 ns of each simulation. The gating residues are generally more open in the apo simulations, compared to those where R-warfarin is bound. This is in contrast to the membrane-bound simulations (see 4), where the gating residues display similar behaviour for the apo and R-warfarin bound models.

Figure 6. Average Mulliken spin densities (A) and charges (B) for fragments calculated from QM/MM energy minimizations.

Calculated at the B3LYP-D:6-31G/CHARMM27 level of theory. M_A and M_R denote the membrane apo and R-warfarin bound models, respectively. S_A and S_R denote the solvated apo and R-warfarin bound models, respectively. Small differences are observed between the average values in the calculated properties. However, these differences are of similar order to the variation in values between different time points in a given simulation.

Figure 7. Transition state geometries for hydrogen abstraction pathways with lowest barrier.

Left: profile calculated from membrane simulation; right: profile calculated from water simulation. Calculated at B3LYP-D:6-31G/CHARMM27 level of theory. Arg212 shown explicitly and has been proposed to play a role in substrate specificity of CYP3A4. In the lowest-energy profiles generated it is observed to hydrogen bond with the O2 atom of R-warfarin in the membrane-bound simulation but not in the soluble case where Arg212 shows increased flexibility.

Figure 8. Flow diagram for preparation of protein:membrane model for QM/MM calculations.

Blue stages correspond to preparation steps whereas pink stages involve molecular dynamics (MD) simulations and QM/MM calculations. In the initial step (yellow) the structure of the protein of interest is downloaded from the PDB. The following steps describe the unbiased coarse grained (CG) set up of a protein in the presence of a chosen lipid composition. Following these preparatory stages CGMD simulations are run to allow the bilayer to self-assemble in an unbiased manner around the protein. This CGMD protein-membrane simulation box may be converted to atomistic resolution and atomistic MD simulations performed to look at conformational changes. The outcome of these simulations may then be used as the basis of QM/MM calculations or other calculations of interest.