Cooperative Binding of Aflatoxin B1 by Cytochrome P450 3A4: A

Subscriber access provided by Uppsala universitetsbibliotek

Article

Cooperative Binding of Aflatoxin B1 by Cytochrome P450 3A4: A Computational Study U. Bren, Julian Fuchs, and Chris Oostenbrink Chem. Res. Toxicol., Just Accepted Manuscript • DOI: 10.1021/tx5004062 • Publication Date (Web): 14 Nov 2014 Downloaded from http://pubs.acs.org on November 17, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Chemical Research in Toxicology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

Cooperative Binding of Aflatoxin B1 by Cytochrome P450 3A4: A Computational Study

Urban Bren1,2, Julian E. Fuchs1,3, and Chris Oostenbrink1,*

1

Institute of Molecular Modeling and Simulation, University of Natural Resources and Life Sciences, Muthgasse 18,

AT-1190 Vienna, Austria. 2 3

Laboratory for Molecular Modeling, National Institute of Chemistry, Hajdrihova 19, SI-1001 Ljubljana, Slovenia. Institute of General, Inorganic and Theoretical Chemistry, University of Innsbruck, Innrain 80/82, AT-6020

Innsbruck, Austria

Keywords: aflatoxin

B1 metabolism, cytochrome P450 3A4, positive homotropic cooperativity,

sequential molecular docking, molecular dynamics simulations, free-energy calculations.

*

to whom correspondence should be addressed. E-mail: [email protected]. Telephone: +43-1-

476548302. Fax: +43-1-476548309.

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic

2


Page 2 of 50

Page 3 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Abstract

Aflatoxin B1 (AFB1) – the most potent natural carcinogen known to men – is metabolized by cytochrome P450 3A4 (CYP3A4), either to the genotoxic AFB1 exo-8,9-epoxide or to the detoxified 3α-hydroxy AFB1. The activation of the procarcinogen proceeds in a highly cooperative fashion which differs from common allosteric regulation in the sense that it can be attributed to simultaneous occupancy of a single large and malleable active site by multiple ligand molecules. Unfortunately, unlike in the case of ketoconazole there is currently no experimental structure available of the doublyligated CYP3A4-AFB1 complex. Therefore, we employed a sequential molecular docking protocol to create various possible doubly-ligated complexes and subsequently performed molecular dynamics simulations and free-energy calculations to check for their consistency with the available experimental data on regio- and stereo-selectivity of both AFB1 oxidations as well as with available kinetic data. Only the system in which the first AFB1 molecule was bound in a face-on C8-C9 epoxidation mode and the second AFB1 molecule was bound in a side-on 3α-hydroxylation mode – a result of an unconstrained molecular docking protocol – has successfully fulfilled all the imposed criteria and is therefore proposed as the most likely structure of the doubly-ligated complex of CYP3A4 with AFB1. The empirical Linear Interaction Energy method revealed that shape complementarity through nonpolar dispersion interactions between the two bound AFB1 molecules is the main source of the experimentally observed positive homotropic cooperativity. The reported study represents a nice example of how state-of-the-art molecular modeling techniques can be used to study complicated macromolecular complexes, whose structures have not yet been experimentally determined, and to validate these against the available experimental data. The proposed structure will facilitate future studies on the rational design of successful AFB1 modulators or on human subpopulations characterized by specific CYP3A4 polymorphisms that are especially sensitive to AFB1. 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 50

1. Introduction

Aflatoxin B1 (AFB1) represents the most potent natural carcinogen known to men and is 1000times more carcinogenic than benzo[a]pyrene from cigarette smoke.1 This furofuran mycotoxin is produced by the common molds Aspergillus flavus, Aspergillus parasiticus, and Aspergillus nomius which infest agricultural commodities stored in hot moist conditions and is involved in the etiology of human liver cancer.2 AFB1 constitutes a major public health issue in several areas of the developing world, and represents a major cause of administrative food retraction from the food-store shelves also in the OECD countries.3 Products which are frequently affected include cereals, oilseeds, spices, tree nuts, peanut butter as well as milk, meat and eggs of animals that are fed contaminated feed. AFB1 has also been produced as a biological/chemical weapon, e.g., in 1989 in Iraq.4 To be precise, AFB1 is not carcinogenic per se, but after intake is metabolized by the cytochrome P450 3A4 enzyme (CYP3A4) either to the genotoxic AFB1 exo-8,9-epoxide or to the nontoxic 3α-hydroxy AFB1 (Scheme 1).5 Therefore, AFB1 is referred to as a procarcinogen and its epoxy metabolite is called the ultimate carcinogen. CYP3A4 is a human heme-containing mono-oxygenase enzyme expressed primarily in hepatic tissue and in the small intestine.6 It metabolizes more than 50% of clinically used drugs and is often involved in adverse drug−drug interactions.7,8 CYP3A4 displays atypical kinetic behavior towards AFB1, characterized by a sigmoidal shape of the corresponding titration curves, which is indicative of a positive homotropic cooperativity.9 This requires a participation of at least two ligand molecules, where the binding of the first ligand molecule increases the affinity of CYP3A4 for the binding of the second ligand molecule.10,11 However, the observed CYP3A4 cooperativity differs from common allosteric regulation in the sense that it can be attributed to simultaneous occupancy of a single large and malleable active site by multiple ligand molecules.12

4


Page 5 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Consequently, its mechanistic origin is poorly understood and involves a direct interaction between the bound ligand molecules and/or selecting and locking a specific preexisting conformation of CYP3A4 with a higher affinity for the second ligand molecule.13 Interestingly, the degree of cooperativity was found to be practically identical in the case of AFB1 epoxidation and AFB1 3α-hydroxylation.9 Moreover, addition of flavonoids like α-naphthoflavone (Scheme 1) was shown to both stimulate AFB1 epoxidation and inhibit AFB1 3α-hydroxylation while concurrently reducing the degree of positive homotropic cooperativity.9 Finally, aflatoxin B2 (AFB2), which differs from AFB1 only in the presence of a saturated 8,9-bond (Scheme 1), is not oxidized by CYP3A4 and is at the same time capable of inhibiting both AFB1 oxidations.9 In the present study, we addressed the process of cooperative AFB1 activation by CYP3A4 using various computational means. In the past, computational approaches to study the cooperative binding for other Cytochrome P450 enzymes have been reported.14,15 Three plausible doubly-ligated CYP3A4 complexes with AFB1 were initially created by molecular docking techniques. Subsequent molecular dynamics (MD) simulations in conjunction with free-energy calculations were used to judge their correctness. Fast empirical free-energy methods have been shown to be of use in the determination of binding modes and molecular interactions in Cytochrome P450 enzymes before.13,16-18 Trajectories of mono- and doubly-ligated CYP3A4 complexes with AFB1 were carefully compared in terms of structural and thermodynamic properties to identify the docked structure best fitting to the available experimental data. Thereby, it should be clearly stated that due to its large, malleable, and heme-containing active site CYP3A4 represents a challenging system to simulate.19 However, understanding the molecular basis of its cooperative behavior has the potential to facilitate the rational design of successful AFB1 modulators as well as to identify the especially AFB1 endangered human subpopulations characterized by specific CYP3A4 polymorphisms. In a broader context it may also

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 50

help to understand the often-observed CYP3A4-related adverse drug-drug interactions, making it a fundamental question of both scientific and pharmaceutical importance.

2. Computational Methods 2.1. Initial Structure Preparation and Molecular Docking

Since no experimental CYP3A4 structure in complex with AFB1 is currently available, we had to resolve to the use of molecular docking techniques.20 The crystal structure of CYP3A4 in complex with the antifungal drug ketoconazole (PDB code 2V0M)12 served as our starting point, because positive homotropic cooperativity has been unambiguously established by both experimental21 and theoretical13 techniques for this protein-ligand complex. In this structure, the asymmetric unit cell consists of four copies of the enzyme of which the N-terminal transmembrane helix was truncated to allow for crystallisation. Monomer A was selected based on the good resolution and the two bound ketoconazole molecules were removed from it. Coordinates for the missing loops Thr-264 till Arg-268 and Asn-280 till Lys-288 were obtained from monomer C. Missing side-chain atoms of Lys-168, His287, and Asp-497 were manually added. All Lys, Arg, and Cys (with exception of Cys-442 coordinating the heme iron on its proximal side) side-chains were protonated, all Asp and Glu sidechains were deprotonated, the amino and carboxy termini were charged, the His-65 side-chain was charged, the His-30 side-chain was Nε-protonated, and all other His side-chains were Nδ-protonated based on visual inspection of hydrogen bonding networks. This procedure yielded the initial protein structure for the subsequent molecular docking. The observed CYP3A4 cooperativity differs from the common allosteric regulation in the sense that it can be attributed to simultaneous occupancy of a single large and malleable active site by multiple ligand molecules.12 The AFB1 docking was therefore performed in two steps – the first AFB1 6


Page 7 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecule was docked to the prepared apo CYP3A4 structure yielding the mono-ligated complex, which in turn served as starting point for the docking of a second AFB1 molecule. Docking was performed using GOLD (Genetic Optimization for Ligand Docking),22 version 3.3.1, in combination with the Chemscore scoring function23 parametrized for heme-containing proteins.24 The center point for docking was obtained by mirroring the sulfur atom of the proximal Cys-442 over the planar ferricprotoporphyrin IX heme complex. The radius of the docking-grid sphere was set to 1.5 nm to accommodate the entire AFB1 molecule as long as at least one of its atoms is placed less than 0.6 nm away from the heme iron – considered to be the upper limit for a reactive distance in our previous work.25,26 At most, 100000 operations in the genetic algorithm were performed using a population of 100 genes. At least five independent docking simulations were performed to reach statistical significance. 50 poses were stored from every docking simulation, except if the best three docking poses had root-mean-square deviations smaller than 0.15 nm. To prevent an early convergence, a relative pressure of 1.1 was specified, and to account for diversity, the number of niches was set to 2.

2.2. Molecular Dynamics Simulations

All simulations were performed and analyzed using the GROMOS11 simulation package.27 The GROMOS 45A4 force field was applied.28,29 AFB1 interaction parameters were selected based on analogy with similar groups represented in the GROMOS 45A4 force field and are available in Table S1 in the Supporting Information. The GROMOS force field is a united-atom force field, which means that aliphatic hydrogens are implicitly treated together with the carbon atom to which they are attached. The mono- and doubly-ligated CYP3A4 complexes with AFB1 were energy-minimized and placed into a periodic pre-equilibrated rectangular box of SPC water.30 Minimum solute to solvent and minimum solute to box-wall distances were set to 0.23 and 0.8 nm, respectively. Positions and geometries of 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 50

water molecules were energy-minimized using the steepest-descent algorithm. Subsequently, four of them were selected based on their electrostatic potential and replaced by chloride ions to achieve electroneutrality of the studied complexes. In addition, AFB1 was simulated as a single molecule, solvated in a periodic box of SPC water molecules. The following equilibration protocol was used: Initial velocities were randomly generated from a Maxwell-Boltzmann distribution at 50 K. All solute atoms were restrained to their initial positions through a harmonic potential with force constant of 2.5 × 104 kJ mol-1 nm-2 and the systems were propagated for 20 ps. In each of the five subsequent 20 ps MD simulations, the temperature was raised by 50 K and the harmonic force constant of the positional restraints was reduced by one order of magnitude. Then, positional restraints were removed, rototranslational constraints on all solute atoms were introduced,31 and the systems were simulated for 20 ps at 300 K. Finally, a simulation at a constant pressure of 1 atm was performed for 300 ps. Production runs of 10 ns at 300 K and atmospheric pressure followed. System coordinates and energies were stored every 0.5 ps for further analysis. Temperature and pressure were kept constant using the weak coupling scheme32 with relaxation times of 0.1 and 0.5 ps, respectively. Separate temperature baths were used for solute and solvent to avoid freezing of the protein degrees of freedom and an isothermal compressibility of 4.575 × 10-4 kJ-1 mol nm3 was applied. Bond lengths were constrained using the SHAKE algorithm33 allowing for time steps of 2 fs. The rototranslational constraint was maintained to ensure that the protein never interacts with one of its periodic copies. Nonbonded interactions were calculated using a triple range scheme. Interactions within a short-range cutoff of 0.8 nm were calculated at every step from a pair list that was updated every fifth step. At these points, interactions between 0.8 and 1.4 nm were also calculated and kept constant between updates. A reaction field34 contribution was added to the electrostatic interactions and forces to account for a homogeneous medium outside the long-range cutoff using a relative dielectric constant of 61 as 8


Page 9 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


appropriate for the SPC water model.35

2.3. Free-Energy Calculations

The Linear Interaction Energy (LIE) method of Åqvist and coworkers36 represents a rapid approach to evaluate solvation free energies. It is based on a modified linear response to treat electrostatic interactions and on an empirical term treating the dispersion interactions. The solvation L(S ) free energy ΔGsolv of a solute L in a surrounding solvent S is expressed as:

L(S ) L(S ) ΔGsolv = α VvdW + β VesL(S )

L(S ) where VvdW

(1)

and VesL(S ) represent van der Waals and electrostatic interaction energy between the

solute L and its surrounding solvent S, respectively, averaged over an ensemble of configurations generated by MD simulations. α and β are coefficients of the LIE method, where β has the theoretical L(S ) value of 0.5.36 This approach directly yields contributions to ΔGsolv due to dispersion and electrostatic

interactions. LIE is a variant of the Linear Response Approximation, in which the electrostatic preorganization energy is included in equation (1) and the theoretical background for β = 0.5 is more pronounced. Inclusion of the electrostatic preorganization may be relevant for some enzymes,37,38 but in practice, α and β are often determined empirically, 39 whereby an alternative value of β implicitly includes the electrostatic preorganization.40,41 If the interaction energy V can be decomposed into additive atom-group contributions, the resulting solvation free energy can readily be interpreted in terms of interactions between specific groups of atoms in the system.42,43 Cooperativity of AFB1 binding to CYP3A4 was studied with the help of a dedicated thermodynamic cycle in detail described in Reference 13. Doubly-ligated CYP3A4 (PA1A2) can be obtained from the apo CYP3A4 (P) by the binding of the first AFB1 molecule at its position 1 leading 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 50

to formation of a mono-ligated CYP3A4 (PA1) followed by the binding of the second AFB1 molecule A1/P at its position 2. The free energy of the initial binding ΔGbind can be calculated as a difference between A1(P ) the solvation free energies of AFB1 at the position 1 in the apo CYP3A4, ΔGsolv , and of AFB1 in A(W ) 36,44,45 aqueous solution, ΔGsolv : A1/P A1(P ) A(W ) ΔGbind = ΔGsolv − ΔGsolv

(2)

A2/PA1 The free energy of the subsequent binding ΔGbind can be calculated as a difference between the A2(PA1) solvation free energies of AFB1 at the position 2 in a mono-ligated CYP3A4, ΔGsolv , and of AFB1 A(W ) in aqueous solution, ΔGsolv . A correction for the altered solvation of AFB1 at position 1 upon the

transition from mono- to doubly-ligated CYP3A4 also needs to be added: A2/PA1 A2(PA1) A(W ) A1(PA2) A1(P ) ΔGbind = ΔGsolv − ΔGsolv + ΔGsolv − ΔGsolv

(3)

Alternatively, doubly-ligated CYP3A4 (PA1A2) can also be obtained from the apo CYP3A4 (P) by the binding of the first AFB1 molecule at the position 2 leading to formation of a mono-ligated CYP3A4 (PA2) followed by the binding of the second AFB1 molecule at the position 1. The A2/P A1/PA2 corresponding free energies ΔGbind and ΔGbind are defined analogously as: A2/P A2(P ) A(W ) ΔGbind = ΔGsolv − ΔGsolv

(4)

A1/PA2 A1(PA2) A(W ) A2(PA1) A2(P ) ΔGbind = ΔGsolv − ΔGsolv + ΔGsolv − ΔGsolv

(5)

A1/P A1/PA2 A2/P A2/PA1 The difference between ΔGbind and ΔGbind or between ΔGbind and ΔGbind offers a

quantification of the cooperative binding.

10


Page 11 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Results and Discussion 3.1. Docked and Simulated Structures

To reduce the large uncertainties related to the subsequent docking of two substrate molecules, a working hypothesis was used that suggests that both AFB1 epoxidation and AFB1 3α-hydroxylation are initiated from a single common doubly-ligated structure of CYP3A4 (PA1A2). This hypothesis can be corroborated by the practically identical degree of cooperativity observed for both competing reactions experimentally.9 Our docking efforts were eased by the fact that AFB1 represents a rather rigid molecule while the preexisting knowledge on the regio- and stereoselectivity of both oxidations allowed us to judge the quality of the obtained structures. In the following, we discuss criteria for the binding modes of AFB1, which we characterize as the epoxidation or the hydroxylation mode of the substrate, according to the product that is most likely to follow from a given pose. The doubly-ligated AFB1 CYP3A4 complex (PA1A2) should fulfill three additional criteria: (i) The AFB1 molecule bound in the C8-C9 epoxidation mode should prevent the 3α-hydroxylation of the second AFB1 molecule, thereby rationalizing the ability of AFB2 (which differs from AFB1 only in the presence of a saturated 8,9-bond) to inhibit both AFB1 oxidations.9 (ii) The AFB1 molecule bound in the 3α-hydroxylation mode could be replaced by the planar molecule of α-naphthoflavone thereby explaining the capability of α-naphthoflavone to both stimulate the AFB1 epoxidation and inhibit the AFB1 3α-hydroxylation while concurrently reducing the degree of positive homotropic cooperativity.9 (iii) The two bound AFB1 molecules should form direct interactions with CYP3A4 amino-acid residues Leu-210, Leu-211, and Phe-304 that were established as indispensable for the positive homotropic cooperativity of both AFB1 oxidations.46 Complex PA1A2_1 depicted in Figure 1A is a result of an unrestrained docking experiment. The first AFB1 molecule is bound in a face-on47 C8-C9 epoxidation mode, while the second AFB1 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 50

molecule adopts a side-on47 3α-hydroxylation mode. The position of the heme iron relative to the AFB1 molecule would lead to the correct stereochemistry for both oxidations and the complex fulfills the three additional criteria set above. Figure 1B depicts complex PA1A2_2 – the first AFB1 molecule is bound in a face-on 3α-hydroxylation pose and the second AFB1 molecule shows a face-on C8-C9 epoxidation pose. The structure is a result of a restrained docking protocol whereby the distances from the 3α-hydrogen of the first AFB1 molecule as well as from the C8 and C9 hydrogens of the second AFB1 molecule to the heme iron were penalized by a harmonic restraint with a force constant of 500 kJ/(mol nm2) if they exceeded 0.3 nm. The complex again exhibits correct stereochemistry for both oxidations and fulfills the three additional criteria set beforehand. Figure 1C shows complex PA1A2_3 – the first AFB1 molecule is bound in a face-on 3α-hydroxylation mode and the second AFB1 molecule adopts a side-on C8-C9 epoxidation mode. The structure is a result of a restrained docking protocol whereby the distances from the 3α-hydrogen of the first AFB1 molecule as well as from the C8 hydrogen of the second AFB1 molecule to the heme iron were penalized by a harmonic restraint with a force constant of 500 kJ/(mol nm2) if they exceeded 0.3 nm. The complex exhibits wrong stereochemistry for the epoxidation mode which in addition fails to prevent the 3α-hydroxylation from taking place. However, this structure is still considered as these issues might be resolved during subsequent MD simulations. During the 300 ps of unrestrained equilibration phase of MD simulation of complex PA1A2_2 significant structural rearrangements occurred. The first AFB1 molecule shifted to a side-on pose exhibiting wrong regioselectivity for the 3α-hydroxylation. The resulting system is depicted in Figure 2A and will be termed PA1A2_2u. Consequently, the equilibration phase of complex PA1A2_2 was repeated using a restrained protocol whereby distances from the C3 atom of the first AFB1 molecule as well as from the C8 atom of the second AFB1 molecule to the heme iron larger than 0.4 nm were penalized by a harmonic restraint with a force constant of 500 kJ/(mol nm2). This procedure resulted in 12


Page 13 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a structure generally resembling the complex PA1A2_3 – it is depicted in Figure 2B and will be termed PA1A2_2r. The above mentioned restrains were turned off during the subsequent MD production run. MD simulations of five mono-ligated CYP3A4 complexes with AFB1 were also performed: system PA1_1 containing the first AFB1 molecule from the PA1A2_1 complex bound in a face-on C8-C9 epoxidation mode; system PA2_1 containing the second AFB1 molecule from the PA1A2_1 complex bound in a side-on 3α-hydroxylation mode; system PA1_3 containing the first AFB1 molecule from the PA1A2_3 complex bound in a face-on 3α-hydroxylation mode; system PA2_3 containing the second AFB1 molecule from the PA1A2_3 complex bound in a side-on C8-C9 epoxidation mode; and system PA1_2u containing the first AFB1 molecule from the PA1A2_2u complex bound in a side-on pose exhibiting wrong regioselectivity for the 3α-hydroxylation. Systems PA1_2r and PA2_2r were not simulated due to their general similarity with complexes PA1_3 and PA2_3, respectively, nor was system PA2_2u due to its general similarity with complex PA1_1.

3.2. Structural Analysis

The initial crystal structure (PDB code: 2V0M) as well as the final structures of our simulated systems (PA1A2_1, PA1A2_2u, PA1A2_2r, PA1A2_3, PA1_1, PA2_1, PA1_3, PA2_3, PA1_2u), of apo CYP3A4 (P), of mono-ligated (PK1 and PK2) and doubly-ligated (PK1K2) CYP3A4 ketoconazole complexes13 were gathered and the root-mean-square deviation (RMSD) was calculated between backbone atoms of all pairs of structures after a rototranslational fit to remove the translational and rotational motion (Figure S1). During 10 ns MD simulations the RMSD of the CYP3A4 with respect to its initial crystal structure configuration settles in a span from 0.19 nm for system PA2_1 till 0.28 nm for system PA1A2_2u. These RMSD values are larger than the backbone RMSD of our initial structure with respect to the two available human apo CYP3A4 crystal structures (1TQN48 of 0.18 nm and 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 50

1W0E49 of 0.16 nm), or with respect to the four mono-ligated human CYP3A4 crystal structures in complex with progesterone (1W0F49 of 0.17 nm), with metyrapone (1W0G49 of 0.18 nm), with erythromycin A (2J0D12 of 0.18 nm), and with ritonavir (3NXU50 of 0.07 nm). Therefore, the CYP3A4 conformational space explored during the MD simulations seems to be larger than the one covered by the current crystallographic data, although the largest backbone RMSD value observed between the final structures of simulated systems PA1A2_2u and PA1A2_2r of 0.36 nm is comparable to the one observed in MD studies of CYP3A4 ketoconazole complexes (RMSD=0.37 nm). Pairs of CYP3A4 structures with RMSD values smaller than 0.21 nm were considered structural neighbours and the structure with the largest number of neighbours was considered the central member of a first cluster.51 The first cluster consisted of five members: structures PK2, PK1K2, PA1A2_1, PA1_1, and PA2_1. The presence of the second (K2) ketoconazole molecule was found to be responsible for maintaining the CYP3A4 conformation close to its original crystal structure13 and the fact that structures PA1A2_1, PA1_1 and PA2_1 are similar to PK2 and PK1K2 provides additional confidence in the CYP3A4 AFB1 complexes obtained using the unrestrained docking protocol. Similar results were obtained also in the case when clustering was performed on the backbone atoms of the CYP3A4 active-site amino-acid residues only. Root-mean-square fluctuations (RMSF) of the CYP3A4 backbone were calculated over the 10ns MD simulations as well. In accordance with previous results13 on ketoconazole CYP3A4 complexes, both the N- and C-termini as well as the modeled loops Thr-264 till Arg-268 and Asn-280 till Lys-288 exhibit large fluctuations along with the flexible solvent-exposed loops Arg-105 till Gly-109, Arg-161 till Lys-169, Asp-194 till Gln-200, Lys-342 till Thr-346, and Glu-417 till Asp-428. As previous MD simulations revealed the CYP3A4 active-site to be extremely malleable,13,52 the volume of the cavity at the distal side of the ferric-protoporphyrin IX heme complex was monitored using the CASTp server.53 Double-ligated AFB1 CYP3A4 complexes were found to have somewhat larger average active-site 14


Page 15 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


volumes (3.0 nm3) than the mono-ligated AFB1 CYP3A4 complexes (2.7 nm3). Upon transition from doubly-ligated to mono-ligated CYP3A4 complexes one AFB1 molecule is on average replaced by 9 water molecules. The water ingress channels54 remain open and the majority of the water molecules experience exchange with bulk solvent during the course of 10 ns MD simulations. The evolution of the CYP3A4 secondary-structure elements55 was followed as well – regardless of the complex they remain stable over simulation time (data not shown).

3.3. Important Interactions, Reactive Distances and Orientations, Conformational Clustering

The presence of hydrogen bonds was monitored over the 10 ns MD simulations using a geometric criterion. A hydrogen bond is considered to be formed if the distance between the hydrogen and acceptor atoms is less than 0.25 nm and the donor-hydrogen-acceptor angle is larger than 135º. The backbone of Arg-106 is hydrogen-bonded to the O1 atom of the AFB1 molecule (Scheme S1 in the Supporting Information) in the PA1_1 complex for 54% of the time and the side-chain of Ser-119 is hydrogen bonded to the O1 atom of the second AFB1 molecule in the PA1A2_3 complex for 25% of the time as well as to the O1 (for 17% of the time) and O11 (for 79% of the time) atoms of the second AFB1 molecule in the PA1A2_2u complex. Both have also been identified as key residues for progesterone binding.56 The side-chain of Arg-105 forms the largest number of hydrogen bonds: with the O11 atom of the AFB1 molecule in the PA1_1 complex for 78% of the time; with the O6A (for 40% of the time) and O7 (for 36% of the time) atoms of the AFB1 molecule in the PA2_3 complex; with the O11 atom of the AFB1 molecule in the PA1_2u complex for 39% of the time; and with the O11 atom of the first AFB1 molecule in the PA1A2_3 complex for 24% of the time. Moreover, hydrogen bonds between AFB1 moieties and various active-site water molecules are also formed. Furthermore, aromatic hydrogens of Phe-304 satisfy the above stated criteria with the O11 atom of the AFB1 molecule in the 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 50

PA2_3 complex for 12% of the time; with the O1 atom of the second AFB1 molecule in the PA1A2_1 complex for 26% of the time; as well as with the O1 atom of the second AFB1 molecule in the PA1A2_2r complex for 12% of the time. This residue has been experimentally established as indispensable for the positive homotropic cooperativity of both AFB1 oxidations46 as well as has been identified by means of MD simulations as influencing the cooperative behavior of diazepam binding to CYP3A4 through stacking interactions.57 The presence of stacking interactions was monitored over the 10 ns MD simulations using a geometric criterion. Ring systems are considered to stack if the distance between the centers of geometry of the rings is less than 0.5 nm and the angle between the planes through the two rings is less than 30 degrees. In the PA1A2_1 complex the furofuran region of the first AFB1 molecule stacks with the CYP3A4 heme group for 63% of the time and the two AFB1 molecules are stacked in an antiparallel orientation for 99% of the time. Similarly, the furofuran region of the second AFB1 molecule of the PA1A2_2u complex also stacks with the CYP3A4 heme group for 93% of the time, but the two AFB1 molecules are stacked in a parallel orientation for 99% of the time. Moreover, the planar region of the first AFB1 molecule stacks with the CYP3A4 heme group for over 90% of the time and the two AFB1 molecules are stacked in a parallel orientation in complexes PA1A2_2r and PA1A2_3. Furthermore, the planar region of the AFB1 molecule stacks with the CYP3A4 heme group for 99% of the time in complex PA1_3 and the furofuran region of the AFB1 molecule stacks with the CYP3A4 heme group for over 80% of the time in complexes PA1_1 and PA1_2u. Finally, occasional stacking interactions of the AFB1 molecule with the Phe-304 residue of CYP3A4 occur in systems PA2_1 and PA2_3. Reactive distances and orientations were also carefully followed throughout the 10 ns MD simulations (Figures 3 and 4). In accordance with our previous work25,26 reactive centers closer than 0.6 nm to the heme iron were considered oxidizable. The stereochemical correctness was determined by 16


Page 17 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


monitoring the sign of the improper dihedral angle defined by three ring atoms of AFB1 and the heme iron. For epoxidation, this was the improper dihedral angle defined by atoms C2A-C3A-C3-Fe and for hydroxylation we used C9-O7-C8A-Fe (See scheme 1 in the supporting information for atom names). More elaborate ways to study the reactivity would involve some means of quantum mechanical calculations,19,58,59 but for the current purpose of determining plausible binding poses, a simple geometric criterion suffices. In Figures 3 and 4, distributions of the distances between the sites of metabolism and the heme iron are drawn in solid lines when the orientation corresponds to correct stereochemical product formation, in dotted lines when the orientation would lead to the incorrect stereo-product and in dashed lines when both orientations are observed. The first face-on AFB1 molecule of the PA1A2_1 complex has the C8 and C9 atoms within the reactive distance for epoxidation for 99% of the time and solely exhibits the correct exo-stereochemistry. The second sideon AFB1 molecule of the PA1A2_1 complex has the C3 atom within the reactive distance for hydroxylation for 5% of the time and preferentially samples the correct α-stereochemistry. The absence of the second AFB1 molecule in the monoligated PA1_1 complex has virtually no effect on the reactive distances and stereochemistry of the remaining face-on AFB1 molecule in the epoxidation mode, while the absence of the first AFB1 molecule in the monoligated PA2_1 complex facilitates the approach of the remaining side-on AFB1 molecule in the hydroxylation mode towards the heme iron and brings its C3 atom within the reactive distance for 51% of the time. This is in line with the requirement that the binding pose leading to epoxidation should be able to inhibit hydroxylation, offering an explanation of the experimental finding that AFB2 inhibits both reactions. The doubly-ligated complexes PA1A2_2r and PA1A2_3 show a very similar behavior: The first face-on AFB1 molecule places the C3 atom within the reactive distance for hydroxylation for 99% of the time and solely exhibits the correct αstereochemistry. The second side-on AFB1 molecule places its C8 and C9 atoms within the reactive distance for epoxidation for 15 and 6% of the time for complexes PA1A2_2r and PA1A2_3, 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 50

respectively, and exclusively displays the incorrect endo-stereochemistry. The absence of the second AFB1 molecule in the monoligated PA1_3 complex has virtually no effect on the reactive distances and stereochemistry of the remaining face-on AFB1 molecule in the hydroxylation mode, while the absence of the first AFB1 molecule in the monoligated PA2_3 complex facilitates the approach of the remaining side-on AFB1 molecule in the epoxidation mode towards the heme iron and brings its C8 and C9 atoms within the reactive distance for 99% of the time also occasionally providing the correct exostereochemistry for the reaction. Finally, the first side-on AFB1 molecule of the complex PA1A2_2u has all of its reactive centers more than 0.6 nm away from the heme iron during the entire 10 ns MD simulation, while the second face-on AFB1 molecule has the C8 and C9 atoms within the reactive distance for epoxidation for 99% of the time and solely exhibits the correct exo-stereochemistry. The absence of the second AFB1 molecule in the monoligated PA1_2u facilitates the approach of the remaining AFB1 molecule towards the heme iron placing it in a face-on mode with its C8 and C9 atoms within the reactive distance for epoxidation for 45% of the time and with the correct exostereochemistry. This pose greatly resembles the complex PA1_1 – a fact that can also be confirmed by the significant decrease of the RMSD of the AFB1 atoms in the PA1_2u system with respect to their final structure in the PA1_1 system from 0.38 to 0.27 nm during the course of the 10 ns MD simulation after a rototranslational fit based on the backbone atoms of CYP3A4. All in all, complex PA1A2_1 represents the only system that agrees with the experimentally determined regio- and stereo-selectivity of both oxidations. 2000 AFB1 structures were collected for each of the thirteen (A1_PA1A2_1, A2_PA1A2_1, A1_PA1A2_3, A2_PA1A2_3, A1_PA1A2_2r, A2_PA1A2_2r, A1_PA1A2_2u, A2_PA1A2_2u, PA1_1, PA2_1, PA1_3, PA2_3, and PA1_2u) CYP3A4-bound AFB1 poses at 5 ps intervals of the MD simulations and the RMSD of the AFB1 molecules was calculated between all pairs of structures after a rototransalational fit based on the backbone atoms of CYP3A4. Pairs of AFB1 structures with RMSD 18


Page 19 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


values smaller than 0.25 nm were considered structural neighbors and the structure with the largest number of neighbors was considered the central member of a first cluster.51 After removing all structures belonging to this first cluster the procedure was iteratively repeated until no structure was left and 23 clusters were formed. The population of the seven largest clusters with overall more than 1000 structures is presented in Table 1 and the corresponding central member structures are depicted in Figure 5. The first face-on AFB1 molecule of the PA1A2_1 complex can primarily be found in Cluster 2 and so is the monoligated analogue PA1_1. Both the second side-on AFB1 molecule of the PA1A2_1 complex as well as its monoligated analogue PA2_1 form Cluster 4. However, the absence of the first AFB1 molecule in the monoligated PA2_1 state results in an increased conformational freedom of the second AFB1 molecule indicated by an almost identical partition between Clusters 1 and 4. Structural resemblance between complexes PA1A2_3 and PA1A2_2r is reflected by the first face-on AFB1 molecule being almost exclusively found in Cluster 1 and by the second side-on AFB1 molecule being primarily found in Cluster 3. The monoligated reference state of the first face-on AFB1 molecule PA1_3 also primarily belongs to Cluster 1, while the absence of this first AFB1 molecule in the monoligated PA2_3 reference state results in an increased conformational freedom of the second AFB1 molecule indicated by an almost identical partition between Clusters 5 and 7. The first side-on AFB1 molecule of the PA1A2_2u complex can primarily be found in Cluster 6, while its monoligated reference state PA1_2u in line with the increased mobility due to the absence of the second AFB1 molecule crosses from Cluster 6 to Cluster 2. This fact was already described by the significant decrease of the RMSD of the AFB1 atoms in the PA1_2u system with respect to their final structure in the PA1_1 system from 0.38 to 0.27 nm during the course of the 10 ns MD simulation. The second face-on AFB1 molecule of the PA1A2_2u complex is primarily observed in Cluster 5, therefore system PA2_3 is selected as its most appropriate monoligated state. All in all, the AFB1 molecules bound in a face-on fashion largely maintain their original position in the CYP3A4 binding-site upon transition to 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 50

monoligated complexes, whereas the AFB1 molecules bound in a side-on fashion mostly experience significant structural rearrangements.

3.4. Thermodynamic Analysis

The average electrostatic and van der Waals interaction energies of AFB1 in water as well as in the mono- and doubly-ligated complexes are reported in Table 2 along with the corresponding configurational entropies calculated by the Schlitter60 approach. Nonpolar dispersion interactions clearly represent the main driving force of AFB1 binding to CYP3A4. Therefore, shape complementarity taken together with the large active-site size, its flexibility and malleability has the potential to explain the promiscuous behavior of CYP3A4.61 These van der Waals interactions become even more favorable upon transition from mono- to doubly-ligated complexes pointing to the stacking interaction between the two bound AFB1 molecules as the likely source of the experimentally observed positive homotropic cooperativity. On the other hand, electrostatic interactions were found to disfavor AFB1 binding to all of the studied complexes. Although several hydrogen bonds were detected with various active-site residues and water molecules, which facilitate a proper orientation of the AFB1 molecules in the binding pocket, their overall number seems to decrease upon transfer from the aqueous solution into the CYP3A4 active-site. The AFB1 binding process is entropically disfavored from the standpoint of the ligand as indicated by the decrease in the conformational entropy. The configurational entropy provides insight into the conformational restriction due to the protein environment, but on the other hand does not allow for any meaningful comparison with experimentally determined binding entropy, because it lacks contributions from the protein and the displaced water molecules.62,63 Electrostatic and van der Waals interaction energies between AFB1 molecules and their 20


Page 21 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


hydrated protein environment distribute in a normal fashion and therefore we applied the LIE method for the calculation of the corresponding binding free energies (Figure S2 in the Supporting Information). The only notable exception represents the first AFB1 molecule of the system PA1A2_2r, where the dissociation of the CYP3A4 active-site hydrogen bond connecting the side-chain of Arg-76 with the backbone oxygen of Ser-90 leads to a decrease in the van der Waals interaction energy and a mutually compensating increase in the electrostatic interaction energy. The LIE model with empirical coefficients α and β of 0.375 and 0.152, respectively, that gave a good agreement with the experiment for ketoconazole binding to CYP3A413 was therefore applied to calculate the corresponding AFB1 binding free energies using Equations 2-5. This choice of α and β parameters is quite similar to the optimal set of parameters in simulations of other cytochromes P450 as well.65 The obtained results collected in Table 3 show a large extent of positive homotropic cooperativity as computed from the A1/P differences in affinity in the singly complexed and in the doubly complexed form, e.g. between ΔGbind A1/PA2 and ΔGbind . The cooperativity ranges from 17 (system PA1A2_2r) to 28 kJ/mol (systems PA1A2_3

and PA1A2_2u) – basically twice the amount of positive homotropic cooperativity established in the case of ketoconazole binding to CYP3A4 using the identical LIE model (10 kJ/mol).13 Experimental reaction kinetics of AFB1 epoxidation and AFB1 3α-hydroxylation by CYP3A4 can best be described using the Hill equation

v=

Vc n K n + cn

(6)

which relates the reaction rate v to the given AFB1 concentration c via a set of kinetic parameters of maximal velocity V (0.9±0.1 and 8.1±1.2 min-1 for AFB1 epoxidation and 3α-hydroxylation, respectively), degree of cooperativity n (2.3±0.3 and 2.1±0.1 for AFB1 epoxidation and 3αhydroxylation, respectively), and the half-saturation concentration K (27±2 and 43±7 µM for AFB1 epoxidation and 3α-hydroxylation, respectively).9 The latter can be used to calculate the approximate 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 50

app binding free energy ΔGbind via app ΔGbind = RT ln K

(7)

which will be lower than the binding free energy of the first AFB1 molecule and larger than the binding free energy of the second AFB1 molecule to CYP3A466 and where R is the universal gas constant and T is the absolute temperature. Consequently, AFB1 epoxidation should fulfill the following criteria epo ΔGbind

A

/P

hydro ΔGbind

A

epo < -26.4 kJ/mol < ΔGbind

/PAepo

A

/PAhydro

hydro and AFB1 3α-hydroxylation ΔGbind

A

/P

< -25.2 kJ/mol

Cooperative Binding of Aflatoxin B1 by Cytochrome P450 3A4: A

Recommend Documents