Computational Investigation of Ligand Binding to the Peripheral Site in

Feb 21, 2017 - Human cytochrome P450 3A4 (CYP3A4) is a major drug-metabolizing enzyme responsible for the metabolism of ∼50% of clinically used drug...
0 downloads 7 Views 1MB Size
Subscriber access provided by Iowa State University | Library

Article

Computational Investigation of Ligand Binding to the Peripheral Site in CYP3A4: Conformational Dynamics and Inhibitor Discovery Hanwen Du, Junhao Li, Yingchun Cai, Hongxiao Zhang, Guixia Liu, Yun Tang, and Weihua Li J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00012 • Publication Date (Web): 21 Feb 2017 Downloaded from http://pubs.acs.org on February 22, 2017

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.

Journal of Chemical Information and Modeling 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 37

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

Journal of Chemical Information and Modeling

Computational Investigation of Ligand Binding to the Peripheral Site in CYP3A4: Conformational Dynamics and Inhibitor Discovery

Hanwen Du, Junhao Li, Yingchun Cai, Hongxiao Zhang, Guixia Liu, Yun Tang*, and Weihua Li*

Shanghai Key Laboratory of New Drug Design, School of Pharmacy, East China University of Science and Technology, Shanghai 200237, China

* Corresponding authors, Tel: +86-21-64250811; Fax: +86-21-64251033; E-mail: [email protected] (W. Li); [email protected] (Y. Tang)

1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 Human cytochrome P450 3A4 (CYP3A4) is a major drug-metabolizing enzyme responsible for the metabolism of ~50% of clinically used drugs and is often involved in drug-drug interactions. It exhibits atypical binding and kinetic behavior towards many ligands. Binding of ligands to CYP3A4 is a complex process. Recent studies from both crystallography and biochemistry suggested the existence of a peripheral ligand-binding site at the enzyme surface. However, the stability of ligand bound at this peripheral site and the possibility of discovering new CYP3A4 ligands based on this site remain unclear. In this study, we employed a combination of molecular docking, multi-paralleled molecular dynamics (MD) simulations, virtual screening, and experimental bioassay to investigate these issues. Our results revealed that the binding mode of progesterone (PGS), a substrate of CYP3A4, in the crystal structure was not stable and underwent a significant conformational change. Through Glide docking and MD refinement, it was found that PGS was able to be stably bound at the peripheral site via contacts with Phe215, Phe219, Phe220 and Asp214. On the basis of the refined peripheral site, virtual screening was then performed against the Enamine database. A total of three compounds were finally found to have inhibitory activity against CYP3A4 in both human liver microsome and recombinant human CYP3A4 enzyme assays, one of which showed potent inhibitory activity with IC50 lower than 1 µM, and two of which exhibited moderate inhibitory activity with IC50 values lower than 10 µM. The findings not only presented the dynamic behavior of PGS at the peripheral site, but also demonstrated the first indication of discovering CYP3A4 2

ACS Paragon Plus Environment

Page 2 of 37

Page 3 of 37

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

Journal of Chemical Information and Modeling

inhibitors based on the peripheral site.

3

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 37

Introduction Cytochromes

P450

(CYPs)

are

a

superfamily

of

heme-containing

monooxygenases that play key roles in the metabolism of a large number of xenobiotics such as pollutants, environmental compounds and drugs, as well as many endogenous compounds.1-3 There are 57 CYP isoforms found in humans. Among them, CYP3A4 is the most abundant one in human liver. It is capable of metabolizing ~50% of clinically used drugs and is considered as one of the major drug metabolizing enzymes.1-5 Thus, interactions of ligands with CYP3A4 must always be considered during the discovery and development of new drugs, because CYP3A4 inhibition could lead to xenobiotic-induced toxicity and drug-drug interactions.6, 7 Nevertheless, CYP3A4 inhibition might be useful in some cases. A combination use of HIV protease inhibitor and CYP3A4 inactivator is currently an effective strategy in the treatment of HIV infection. Two potent CYP3A4 inactivators ritonavir and cobicistat8, 9 act as the pharmacoenhancers to enhance the plasma levels of HIV protease inhibitors, which are otherwise quickly metabolized by CYP3A4. Another case is chemical carcinogens like aflatoxins, which are activated by CYP3A4. Inhibition of CYP3A4 activity is able to facilitate the detoxification and elimination of chemical carcinogens.10 Another intriguing aspect of CYP3A4 is that this enzyme exhibits both homotropic and heterotropic cooperativity toward some ligands,11-14 which were manifested by the atypical substrate oxidation kinetic behavior. Several hypotheses have been proposed to explain the mechanism of CYP3A4 cooperativity. One 4

ACS Paragon Plus Environment

Page 5 of 37

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

Journal of Chemical Information and Modeling

hypothesis is that multiple ligand binding sites exist in CYP3A4. Another proposal is the existence of multiple enzyme populations in dynamic equilibrium. More and more evidence support that more than one binding site exists in CYP3A4. Studies from Guengerich group proposed a three-step substrate/inhibitor binding mode based on the kinetic analysis of CYP3A4.

15,16

They suggested that a ligand first binds at a

peripheral site and then moves to the canonical heme active site. Recently, Davydov et al. used fluorescence resonance energy transfer to identify two separated binding sites in CYP3A4.17 They demonstrated the presence of a high affinity binding site at the

enzyme

periphery.

Subsequently,

the same

group

reported

that

the

effector-induced allosteric transitions occurred in CYP3A4, which provided a decisive support to the paradigm of CYP3A4 allosteric regulation.18

Figure 1. The crystal structure of CYP3A4 complexed with PGS (PDB code: 1W0F) (A) The overall structure of 1W0F. (B) Close-up view at the peripheral site with key residues labeled. Evidence that supports multiple ligand binding sites is also from the crystal structures of CYP3A4 complexes. In 2004, Williams et al.19 firstly reported the crystal 5

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

structure of CYP3A4 in complex with the substrate progesterone (PGS) (PDB20 code: 1W0F), in which PGS is bound at the enzyme surface and near the “Phe-cluster”, which is composed of Phe108, Phe213, Phe215, Phe219, Phe220, Phe241, and Phe304 (Figure 1). This peripheral site is remote from the catalytic site above heme. Very recently, Sevrioukova et al.21 reported the crystal structures of CYP3A4 complexed with PGS under two different conditions (PDB codes: 5A1P and 5A1R), in which PGS is also bound at the peripheral site, but located at the interface between symmetry-related CYP3A4 molecules. These two structures are very similar to the one resolved by Williams et al., which strongly proves that PGS bound at the peripheral site is not an artifact. However, in these crystal structures of CYP3A4-PGS complexes, PGS binds to CYP3A4 only via a hydrogen bond formed between PGS and Asp214. It has been pointed out that PGS bound at the peripheral site might be influenced by crystallization, and without the symmetry-related protein contact, the binding affinity of PGS might be weaker.21 Therefore, whether PGS can stay stably at this peripheral site with the binding mode presented in the crystal structure is questionable. In addition, whether the peripheral site serves as the appropriate ligand binding pocket remains an open question. In this study, we used a combination of molecular docking, molecular dynamics (MD) simulations, virtual screening, and experimental bioassay to investigate these questions. Our results indicated that the binding mode of PGS in the crystal structure is not stable and underwent a significant conformational change. Notably, parts of 6

ACS Paragon Plus Environment

Page 6 of 37

Page 7 of 37

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

Journal of Chemical Information and Modeling

“Phe-cluster” shifted to the heme pocket, which resulted in the surface cleft becoming larger. As a consequence, PGS can move to the inside of the peripheral site and bound comfortably via contacts with Phe215, Phe219, Phe220 and Asp214. On the basis of the refined peripheral binding pocket by MD simulations, virtual screening was performed with aim to discover new CYP3A4 inhibitors. A total of three compounds were finally identified to exhibit inhibitory activity against CYP3A4 with IC50 values lower than 10 µM both in human liver microsome and recombinant human CYP3A4 enzyme assays.

Materials and Methods Protein Preparation The initial structure of CYP3A4 complexed with PGS was taken from the Protein Data Bank (PDB code: 1W0F). The missing coordinates of residues Thr262-Arg268 were reconstructed based on the corresponding positions in the ligand-free CYP3A4 crystal structure (PDB code: 1TQN22) with SWISS-PdbViewer 4.1.23 The missing coordinates of residues Asn281-Lys288 were reconstructed using the Prime module of Schrödinger 2012.24 The fixed structure was then subjected to pretreatment by the Protein Preparation Wizard and MacroModel module, in which hydrogen atoms were added and energy was minimized. The protonation states of ionizable residues were determined by PROPKA 3.0.25

7

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

Molecular Docking Self-docking of PGS into the peripheral site of CYP3A4 was accomplished using Glide 5.8 plugged in the Schrödinger 2012. The center of the grid was placed on the centroid of PGS. Five different sizes of the grid box, namely 10 Å, 12 Å, 15 Å, 18 Å and 20 Å, were taken into consideration in this study. For ligand docking trial, the standard precision (SP) was adopted. The scaling factor for van der Waals (vdW) radii of ligand nonpolar atoms was set to the default value of 0.8. Ligand sampling was set as flexible that allows nitrogen inversions and ring conformations sampling. Thirty poses were generated for each docking run and ranked according to Glide score (G-Score). The root mean squared deviation (RMSD) values between the X-ray ligand and docking poses were computed during the docking. All the output poses were further clustered based on the RMSD values.

MD Simulations MD simulations were executed by Amber12.26 The initial structure of CYP3A4 complexed with PGS for MD simulations was based on the prepared structure of 1W0F. All crystallographic water molecules were retained. The protonation states of the ionizable residues and histidines were determined based on the pKa values calculated by PROPKA25 and the microenvironment. According to the calculation results, His28, His54, His65 and His287 were assigned to be fully protonated at both nitrogen atoms. His30 and His402 were protonated at ε nitrogen and other histidine residues at δ nitrogen atoms. In addition, Asp182 and Glu320 were also protonated 8

ACS Paragon Plus Environment

Page 8 of 37

Page 9 of 37

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

Journal of Chemical Information and Modeling

according to the results. The structural optimization of PGS was conducted at the B3LYP/6-31G* level using Gaussian 03 program [www.gaussian.com]. The restrained electrostatic potential (RESP) fitting procedure27 was used to derive the atomic charge of the optimized ligand. The Amber99SB force field and general Amber force field (GAFF) was used for the protein and ligand, respectively. The system was solvated with TIP3P water28 molecules in a truncated octahedron box with a margin of 10 Å along each dimension and then neutralized by counterions. The force field parameters for the heme group were adopted from our previous study.29 The force constant parameters involving Fe were adopted from the work by Harris et al..30 The detailed procedure for MD simulations was adopted as previously described.31 In brief, all waters and hydrogen atoms were firstly optimized. Energy minimization was then carried out by decreasing the harmonic force restraints to protein heavy atoms. The minimized systems were gradually heated from 0 to 300 K in an NVT ensemble over 40 ps and equilibrated at 300 K for 70 ps. Finally, a 40 ns MD simulation was performed for each system under the NPT ensemble at 300 K and 1 atm. All bonds involving hydrogen atoms were constrained with the SHAKE algorithm.32 The long-range electrostatic interaction was treated using the particle mesh Ewald (PME) method.33 The time step and nonbonding interaction cutoff radius was set as 2 fs and 10 Å, respectively. To test and obtain the stable binding mode of PGS at the peripheral site of CYP3A4, ten independent MD simulations (40 ns for each simulation) were carried 9

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 37

out with different initial seeds.

Clustering analysis The cluster analysis tool Kclust, implemented in the MMTSB package,34 was applied to cluster the MD snapshots based on the RMSD of protein Cα atoms. In this study, the last 4 ns of MD trajectories were considered for clustering. The steps can be summarized as follows. First, the Cα atoms of MD snapshots were aligned with the X-ray structure of the protein. Then, the snapshots were extracted and compared with each other. A snapshot was assigned to a cluster if the RMSD was smaller than a threshold value of 1 Å. Afterwards, the component numbers of each cluster and the representative snapshots which are most close to the average structure were obtained. At last, the largest cluster in each system, in which the representative snapshot belongs to, was chosen to calculate the average structure.

MM-GBSA binding free energy calculation In the present study, the binding free energy between PGS and CYP3A4 was calculated and then decomposed into the contributed residues with the MM-GBSA method,35, 36 according to the following equations:

∆‫ܩ‬௕௜௡ௗ௜௡௚ = ‫ܩ‬௖௢௠௣௟௘௫ − (‫ܩ‬௥௘௖௘௣௧௢௥ + ‫ܩ‬௟௜௚௔௡ௗ )

(1)

∆‫ܩ‬௕௜௡ௗ௜௡௚ = ∆‫ܩ‬ெெ + ∆‫ܩ‬௦௢௟௩ − ܶ∆ܵ

(2)

∆‫ܩ‬ெெ = ∆‫ܩ‬௜௡௧ + ∆‫ܩ‬௘௟௘௖ + ∆‫ܩ‬௩ௗ௪

(3)

∆‫ܩ‬௦௢௟௩ = ∆‫ீܩ‬஻ + ∆‫ܩ‬ௌ஺ = ∆‫ܩ‬ௌ஻ + ߛ(ܵ‫ )ܣܵܣ‬+ ܾ

(4)

where ∆Gbinding represents the binding free energy. ∆GMM is the gas phase molecular 10

ACS Paragon Plus Environment

Page 11 of 37

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

Journal of Chemical Information and Modeling

mechanical energy containing internal (∆Gint), electrostatic (∆Gelec) and van der Waals energy (∆Gvdw). ∆Gsolv denotes the electrostatic solvation energy of the complex, which can be divided into polar and nonpolar components. The nonpolar solvation energy ∆GSA was estimated using the LCPO algorithm37 from the solvent-accessible surface area (SASA, Å2). The electrostatic solvation energy (∆GGB) was calculated by the Generalized Born method. The coefficient surface tension γ and constant b were set to 0.0072 kcal•mol-1•Å2 and 0, respectively.38 T∆S denotes that the entropic contribution upon ligand binding at temperature T. Normal mode analysis was adopted to estimate the entropic contribution using the nmode module of Amber 12. In order to reduce the memory demands, only residues within 12.0 Å around the ligand were included for the entropy contribution calculations. In this study, 400 snapshots extracted from the last 4.0 ns of MD trajectory for the calculations of ∆GMM and ∆Gsolv, and 40 snapshots were used for entropic calculation.

Structure-based virtual screening Virtual

screening

was

performed

against

the

Enamine

database

[http://www.enamine.net], which contains 1,999,274 small molecules. The 3D structures of the small molecules in the database were firstly generated and then were assessed according to Lipinski’s rule of 5.39 The remaining ligands were screened by Glide. The grid of docking was set to the center of mass of PGS in the optimal MD-derived complex (System 10). During the procedure, the residues within 15 Å of the centroid of PGS were defined as the binding pocket. Then, the high throughput 11

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

virtual screening (HTVS), standard precision (SP), and extra precision (XP)40,

Page 12 of 37

41

modes of Glide were successively applied to filter the database. The top compounds were ranked based on the Glide score (G-score). Compounds with scores higher than that of the original ligand PGS were kept for the final MM-GBSA calculation in Prime. Finally, the compounds on the top list were selected for visual analysis.

Biological assays A total of 11 compounds were finally selected for purchasing and bioassay. The purchased chemicals were used without further purification. Compounds 1 and 5 have a purity of 90.1% and 90.6%, compounds 6 and 10 have a purity of 92%, and the other compounds possess a purity of greater than 95%. All compounds were characterized by their 1H NMR and MS spectra, which were provided by the vendor. The 1H NMR spectra were recorded using a 400 MHz Varian INOVA Plus or a 500 Bruker Avance instrument. The LC/MS spectra were measured using Agilent high throughput Liquid Chromatograph Mass Spectrometers with DAD and ELSD.

Human liver microsome inhibition assay The inhibition activity of these compounds against CYP3A4 was performed using a human liver microsome assay. The experiments are performed briefly as follows. In 37°C water bath, 7 different concentrations of test compounds (0, 0.0412, 0.1235, 0.370, 1.111, 3.333, 10.0, 30.0 µM) 42, 43 were respectively incubated for 10 min with the addition of human liver microsome (0.127 mg/mL), NADPH (10 mM) and CYP3A4 probe substrate testosterone (40 µM). The CYP3A4 selective inhibitor 12

ACS Paragon Plus Environment

Page 13 of 37

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

Journal of Chemical Information and Modeling

ketoconazole (3 µM) was used as positive control. All samples including positive controls were conducted in duplicate. In the above system, the reaction was quenched by the addition methanol stop buffer containing corresponding internal standard. After quenching the reaction, the reaction mixture is centrifuged at approximately 4000 rpm for 20 minutes. The activity of CYP3A4 was measured using the substrate testosterone. After centrifugation, the samples were submitted for chromatography– tandem mass spectrometry (LC-MS/MS) method analysis consisting of a Waters Ultra Performance Liquid Chromatography (UPLC) and an Applied BioSystems/MDS SCIEX API 4000 instrument mass spectrometer. The decrease of the metabolite formation in peak area ratio (metabolite/internal standard) compared to that of solvent control was used to calculate the inhibition rates (expressed as % of control activity). SigmaPlot was used to plot the % remaining activity versus the test compound concentrations with non-linear regression analysis. IC50 (test compound concentration which produces 50% inhibition) values were determined using three-parameter logistic equation as shown in the following equation:

y=

௠௔௫

(5)

ೣ ష೓೔೗೗ೞ೗೚೛೐ ଵାቀ಺಴ ቁ ఱబ

where max represents the maximal enzyme activity; x represents the concentrations of test compound and y is the enzyme activity at x. Hillslope is the slope factor and IC50 is the concentration to achieve the half maximal inhibition.

Recombinant human CYP3A4 enzyme inhibition assay 13

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

The compounds with IC50 less than 20 µM by human liver microsome assay were further validated by recombinant human CYP3A4 enzyme inhibition assay. The compounds were dissolved and diluted sequentially in DMSO to ensure that the final DMSO concentration was 10 mM in each sample. The CYP3A4 selective inhibitor ketoconazole (3 µM) was used as positive control. The incubation mixture consisted of 0.25 mM/ml MgCl2 potassium phosphate buffer (PB, pH, 7.4), 4.06 nM/ml recombinant CYP3A4, 40 mM testosterone, and 20 µL of NADPH cofactor in a final volume of 5.0 ml, with various inhibitor concentration of 0.0412–30 mM. Then, the reaction was initiated by addition of NADPH after 10 min of pre-warming at 37°C and was terminated 10 minutes. After incubation, the reaction was terminated by adding 400 µL of cold stop solution into the mixtures. All samples including positive controls were conducted in duplicate. After being centrifuged at 4000 rpm at 4°C for 20 min, the clear 200 µL of supernatant was transferred into 100 µL of HPLC water directly and shake for 10 min for LC-MS/MS analysis.

Results and Discussion Self-docking based on crystal structure In the crystal structure of CYP3A4 complexed with PGS (1W0F), PGS binds at the peripheral site located at the surface of CYP3A4, which comprises Phe213, Asp214, Phe215, Asp217, Phe219, Phe220, Ile238, Val240, and Phe241. PGS forms a hydrogen bond with the backbone oxygen atom of Asp214. Apart from this direct

14

ACS Paragon Plus Environment

Page 14 of 37

Page 15 of 37

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

Journal of Chemical Information and Modeling

hydrogen bond, it appears that no additional interactions can stabilize the binding of PGS. To test whether the binding mode of PGS is reasonable, the self-docking of PGS into the peripheral site was first performed. Five different sizes of docking grid, namely, 10 Å, 12 Å, 15 Å, 18 Å and 20 Å, were considered in this study.

Figure 2. The self-docking of PGS into the peripheral site in the crystal structure of CYP3A4 (PDB code: 1W0F). The crystal structure was represented in gray and the docking poses were shown in other colors. The docking trials revealed that no outputs were generated when the docking grid sizes were set to 10 and 12 Å. The other three sizes can successfully generate the docking results, as shown in Figure 2. The poses can be roughly divided into three categories: 1) the poses that kept the same orientation with the crystallographic ligand; 2) the poses that kept the opposite orientation to the crystallographic ligand; and 3) the poses other than the above two kinds. Although PGS was successfully docked at 15

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 37

the peripheral site, the docking output poses were largely different from the original ligand in the crystal structure. In addition, the docking scores were not satisfactory. With the default value of 0.8 used for the vdW radii scaling factor of ligand nonpolar atoms, the best docking score was only -4.648 and the average score was -3.799, as shown in Table 1. Using the grid size of 15 Å, the self-docking of PGS was also performed with different parameter settings, including scaling factor values of 0.7 and 0.6 as well as docking with the extra precision mode. The similar results were obtained (Table S1 in supporting information).

Table 1. Self-docking results based on the CYP3A4 crystal structure (1W0F) Grid size

BM1a

BM2b

Others

numbers

Average score

Lowest score

numbers

Average score

Lowest score

numbers

15 Å

8

-3.799

-4.648

7

-3.600

-4.262

7

18 Å

13

-3.984

-4.375

14

-4.115

-4.393

2

20 Å

13

-3.731

-4.567

11

-4.029

-4.471

5

a

Binding mode 1:the orientation similar to that of PGS in the crystal structure;

b

Binding mode 2: the orientation opposite to that of PGS in the crystal structure. In terms of the docking poses and scores, we concluded that the binding modes

of PGS in the crystal structure (1W0F) cannot be reproduced by Glide docking. Therefore, MD simulation was used to check whether PGS was stably bound at the peripheral site as well as the stability of the residues that constitute the peripheral site.

MD Simulations In order to examine the stability of PGS bound at the peripheral site, MD simulations were performed on the CYP3A4-PGS complex. Ten MD runs were 16

ACS Paragon Plus Environment

Page 17 of 37

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

Journal of Chemical Information and Modeling

carried out, each of which lasted 40 ns. The RMSD values of protein backbone atoms and PGS relative to the initial structure were calculated to examine the stability of the systems during the MD simulations. The RMSD variations with respect to simulation time were monitored as shown in Support Information Figure S1. From Figure S1, it is apparent that the RMSD values of protein backbone atoms of all ten systems converged to ~2 Å. By contrast, the RMSD values of the ligand exhibited significant changes relative to the initial structure. The ligand in the systems 6 and 7 completely dissociated from the binding site, which were manifested by the RMSD values of over 10 Å. Therefore, these two systems were excluded in the subsequent analysis. In the remaining eight systems, the ligand was still bound in the binding site, but the RMSD values had ~5-10 Å deviation from the initial structure. This implies that the ligand underwent significant changes in its binding mode. The hydrogen bond between PGS and Asp214 found in the crystal structure was lost in the remaining systems after MD simulations, except system 10. This also meant that the binding mode of PGS in the crystal structure was not stable. This is not surprising because in the crystal structure PGS binds at the peripheral site only via a hydrogen bond formed between the acetyl oxygen of PGS and the backbone N atom of Asp214. Sevrioukova et al. also pointed out that PGS bound at this peripheral site was stabilized by the symmetry-related protein and PGS molecules and without the symmetry contacts, the interactions of PGS with CYP3A4 might become weaker.21 The steroid rings are totally exposed to the bulky solvent and no strong interactions stabilize the binding of PGS. To identify the conformational changes of PGS after MD simulations, the 17

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

average structures of eight systems were superimposed to the crystal structure, as shown in Figure 3. In all the systems, PGS had very large displacement from the binding mode in the crystal structure. The binding modes of PGS can be roughly classified into three groups. Systems 4, 9 and 10 were classified into Group 1, in which PGS retains a similar orientation with the original binding mode. Systems 2, 3, 5 and 8 constitute Group 2, in which PGS is perpendicular to the original binding mode. Group 3 has only system 1, in which PGS has 180° flip from the one in the crystal structure.

Figure 3. The equilibrated binding modes of PGS after 40 ns of MD simulations. The original binding mode of PGS in the crystal structure was shown in green stick. The binding modes after MD refinement were represented with colorful sticks. Groups 1-3 were shown in yellow, cyan, and rose, respectively. Panels A-H corresponded to systems 1-10, while systems 6 and 7 were excluded. In Group 1, PGS had 4 to 6 Å deviation from its original binding conformation. In this mode, residues Phe108, Phe215, Phe219 and Phe220 underwent conformational changes and pointed to the heme pocket to different levels, which made the space of the peripheral site become larger so that allowed PGS to bind 18

ACS Paragon Plus Environment

Page 18 of 37

Page 19 of 37

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

Journal of Chemical Information and Modeling

tightly. Moreover, Phe219 formed π-π stacking with PGS in the system 10, and the surrounding hydrophobic residues in enzyme periphery made PGS easily fit in. In Group 2, PGS rotated ~90°, which resulted in the methyl ketone of PGS turning over to the other side. The residues Phe108, Phe215 and Phe220 also deflected towards the heme site in systems 4 and 5. Phe219 kept stacking with the steroid rings of PGS, which enhanced the stability of PGS. In system 3, PGS had a flip of ~180° relative to its original mode. Apart from Phe108 and Phe215 stretched to the heme pocket, Phe219 and Phe220 deflected and formed hydrophobic interactions with the ligand. Although the final equilibrated binding modes of PGS differed in various systems, PGS in most of the systems can stay at the peripheral site, which signified the stability of this binding site.

Self-docking based on MD snapshot structures From the above analysis of MD simulations, three groups of binding modes of PGS were found to be stable during MD simulations. In order to select an appropriate model for virtual screening, the representative MD snapshot structures from eight systems were chosen for optimization and used for self-docking of PGS into the peripheral site with Glide. In contrast to the self-docking based on the crystal structure, the self-docking of PGS can well reproduce the original binding modes in the majority of MD generated systems, as shown in Support Information Figure S2. The RMSD values ranged from 0.38 to 2.79 Å. The docking scores of PGS were presented in Table 2. Except for system 4, the docking scores of seven systems were higher than 19

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

that of the crystal structure. Among them, the best result was from system 10, which has the lowest docking score of -7.989. From the docking scores, it is obvious that the binding of PGS in these MD refined structures was more stable than that in the crystal structure. The binding free energies of these eight systems were also calculated (Supporting Information, Table S2). From the results, PGS in system 10 exhibited the lowest binding free energy with CYP3A4. From the self-docking results based on the MD generated snapshot structures, we found that MD simulation was a useful approach that refines the binding of PGS with CYP3A4 at the peripheral site. In fact, at present, MD simulation has become one of widely used protein conformational sampling methods and has been used in many systems including P450s.21, 44, 45 According to the MD simulations, self-docking and the binding free energy based on the MD snapshot structures, system 10 exhibited the best result. Therefore, this model was chosen for analyzing the interactions between PGS and CYP3A4 in detail.

Table 2. Glide docking results based on MD-derived structuresa

a

Systems system 1 system 2 system 3 system 4 system 5 system 8

Lowest docking score -7.084 -4.265 -6.584 -4.778 -5.951 -7.600

RMSD 0.450 0.684 0.552 2.789 0.798 1.042

system 9

-5.171

0.478

system 10 1W0F

-7.989 -4.648

0.382 3.985

The results were obtained with the docking grid size of 15 Å. 20

ACS Paragon Plus Environment

Page 20 of 37

Page 21 of 37

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

Journal of Chemical Information and Modeling

Interactions between PGS and the residues at the peripheral site The crystal structure of CYP3A4 in complex with PGS shows that a hydrogen bond is formed between the acetyl oxygen atom of PGS and the backbone N atom of Asp214. Analysis of the MD trajectories revealed that the hydrogen bond disappeared in systems 3-9 after the systems reached equilibrium, due to the altered binding mode of PGS. By contrast, in system 10, the distance between the two atoms was ~3.2 Å during the last 20 ns of MD simulation, indicating that the hydrogen bond maintained after the system reached equilibrium, which also have contributed to its binding free energy.

Figure 4. (A) Comparison of the MD refined structure (system 10) with the crystal structure. The key residues were shown in lines. PGS in the crystal structure and MD average structure were shown in cyan and purple sticks, respectively. (B) The variation of the distance between the carbonyl oxygen of PGS and the amino nitrogen of Asp214 with respect to simulation time. Superposition of the MD average structure of system 10 with the crystal structure revealed that both PGS and the residues around the peripheral site underwent 21

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

significant conformational changes. Notably, the side chains of Phe215, Phe219 and Phe220 shifted to the heme pocket, which resulted in the surface cleft becoming larger. As a consequence, PGS can move to the inside of the peripheral site and bind with a comfortable mode, as shown in Figure 4. Part of the “Phe-cluster” maintained the hydrophobic interactions with PGS. Especially, Phe219 formed π-π stacking with the hydrophobic part of PGS. At the same time, Asp214 remained the hydrogen bond with the carbonyl oxygen of PGS. Previous studies have shown that some residues that constitute the “Phe-cluster” were reported to be involved in CYP3A4 activity.19, 46-48 For instance, Phe304 was reported to have a dual impact for substrate recognition and cooperative effects.46 Residues such as Phe108, Phe215, Phe220, Phe231 and Phe241 were suggested to act as the “gatekeeper” in the substrates access/egress channels.49, 50 In this study, we found that the “gatekeeper” residues Phe215 and Phe220 altered their conformations and shifted away from the PGS, which tended to open the channel so as to let ligands enter the heme pocket from the enzyme surface. Recently, Nicholas et al. found that the “Phe-cluster” fell into a medium to high exchange regime between the helices F′-G′ and F-G region by the H/D exchange mass spectrometry, which confirmed the enzyme periphery has a moderate degree of solvent accessibility and thus could provide an intermediate site for hydrophobic ligands access to the active site.51 Our results were consistent with their findings. To further identify the interactions between PGS and the peripheral site residues, the binding free energy of PGS with CYP3A4 was calculated with the MM-GBSA 22

ACS Paragon Plus Environment

Page 22 of 37

Page 23 of 37

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

Journal of Chemical Information and Modeling

method. A detailed analysis of the binding free energy showed that the Van der Waals energy was the major driving force for the substrate–enzyme binding. The binding free energy was then decomposed on the individual residue (the detailed data of each key residue energy contributions were provided in Support Information Figure S3). It is obvious that the “Phe-cluster” residues made the large contribution to the binding of PGS. The largest contribution to the binding free energy was from Phe219, which formed π-π stacking with PGS via the hydrophobic ring. Phe220 also has a large contribution to binding free energy of PGS. The contribution of Asp214 is mainly from the electrostatic interaction due to its hydrogen bonding interaction with PGS. This residue has been shown by mutagenesis studies that affected the homotropic and heterotropic cooperativity in CYP3A4.52 Our results coupled with previous studies demonstrated that the “Phe-cluster” has an important effect on the recognition of CYP3A4 with ligands.46-49

Virtual screening based on the peripheral site Ligand binding with CYP3A4 is a complex process together with accommodating multiple molecules simultaneously. In 2004, Williams et al. proposed that the peripheral site could be the effector/substrate recognition region,19 which agreed with structure and mutagenesis data that mutations of Leu210, Leu211 and Asp214 into bulkier residues affected cooperativity in the substrate binding and metabolism.52, 53 The Guengerich group suggested that substrate/inhibitor binds with CYP3A4 via a multiple-steps process.15,

16

They proposed that the substrate or

23

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 24 of 37

inhibitor molecule firstly binds at a peripheral site and subsequently moves toward the active site above the heme, which results in a low to high spin iron shift, then causes a conformational change in the enzyme active site. Recently, a study by Davydov et al. demonstrated the existence of a high affinity binding site at the enzyme periphery with

fluorescence

resonance

energy

transfer,

which

strongly

signified

a

conformational transition triggered by the allosteric ligand binding site.17 Sevrioukova et al. also confirmed that PGS bound at the previously identified peripheral site through co-crystallized with CYP3A4 under two different conditions.21 Our docking and MD simulations results revealed that PGS can stably bind at the peripheral site of CYP3A4, although it displayed multiple binding modes at this site. By combing our results and other groups’ findings,15-18, 21 we proposed that the peripheral site might be a suitable ligand binding site and could be used for ligand design of CYP3A4. To test this hypothesis and discover new CYP3A4 inhibitors, virtual screening was performed against the Enamine database based on this peripheral site. The prepared ligands were docked into the peripheral site based on the model of system 10. The workflow is shown in Figure 5. The high-throughput virtual screening (HTVS) mode was firstly applied to get rid of the molecules that have obvious collision with protein residues. The top 20% of compounds ranked by G-score were retained for subsequent screening. Next, the screening was executed with the standard precision (SP) and Extra precision (XP) modes. A total of 2809 compounds were retained according to the G-score values. Afterwards, Prime MM-GBSA calculation was conducted for these compounds. The top 200 compounds on the MM-GBSA 24

ACS Paragon Plus Environment

Page 25 of 37

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

Journal of Chemical Information and Modeling

scoring list were stored for visual analysis. Based on the compounds’ availability, interaction modes, chemical similarity, and the Pan Assay Interference Compounds (PAINS) check,54 11 compounds were finally selected and tested for the inhibition activity against CYP3A4 with the human liver microsome and recombinant human CYP3A4 enzyme inhibition assays.

Figure 5. The workflow of virtual screening used in this study.

Biological assays The bioassays results of 11 compounds in SMILES format are presented in Table 3. Among the 11 tested compounds, 8 compounds had IC50 values less than 20 µM in human liver microsome assay. It is generally thought that compounds with IC50 < 1 25

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 26 of 37

µM were considered as the strong inhibitors for CYP enzymes, compounds with 1 µM < IC50 < 10 µM as moderate inhibitors, and those with IC50 > 10 µM as weak inhibitors.55 According to this general evaluation criterion for CYP inhibition, compound 10 showed potent inhibitory activity with the IC50 value of 0.142 µM and compounds 2 and 6 exhibited moderate inhibitory activity with the IC50 values between 1.99-8.44 µM. Compounds 3, 4, 5, 7, and 11 had the weakest inhibitory effects. In order to further validate the inhibitory activities of these compounds against CYP3A4, the recombinant human CYP3A4 enzyme inhibition assay was conducted. The results showed compound 10 showed potent inhibitory activity with the IC50 value of 0.0618 µM towards recombinant human CYP3A4 enzyme, and compounds 2, 4, 6, 7, and 11 exhibited moderate inhibitory activities with IC50 values between 2.62-7.81 µM. Compounds 3 and 5 displayed weaker activities with IC50 values of 11.6 and 10.8 µM, respectively. Overall, three compounds displayed moderate to potent inhibitory activities against CYP3A4 in both human liver microsome and recombinant human CYP3A4 enzyme inhibition assays. These results further demonstrated the possibility of the peripheral site serving as a binding site for discovering new inhibitors against CYP 3A4.

Table 3. Activities of CYP3A4 inhibitors discovered via human liver microsome and recombinant human CYP3A4 enzyme assay Cpds

1

Structure

MM-GBSA

IC50

IC50

∆Gbindinga

valuesb

valuesc

-135.6

>72.30

-

SMILES

O=C1NC(N2C1CN(C(C3= CC=C(C4=CC=C(Cl)C=C4 )S3)=O)CC2)=O

26

ACS Paragon Plus Environment

Page 27 of 37

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

Journal of Chemical Information and Modeling

2

O=C(NC1(CSC2=C(C(N3C CN(CC4=CC=C(Cl)S4)CC

-100.0

8.44±0.67

4.25±0.50

-97.6

19.20±3.25

11.60±2.03

-96.1

13.40±1.43

5.68±0.45

-95.6

15.40±1.20

10.80±1.23

-91.6

1.99±0.22

2.62±0.19

-91.2

14.90±6.56

5.28±0.71

-90.2

53.0±8.64

-

-88.9

28.80±6.37

24.70±3.01

-87.9

0.142±0.015

0.0618±0.004

-86.5

13.40±1.04

7.81±1.09

0.0017±0.001

0.0018±0.001

3)=O)C=CC=C2)C)NC1=O 3

O=C(NC1(CSC2=C(C(N3C CC(C(N4CCCCCC4)=O)C C3)=O)C=CC=C2)C)NC1= O O=C(N1)NC(C)(CSC2=C(

4

C(NC3CCN(CC4=CC=C(C l)C=C4)CC3)=O)C=CC=C 2)C1=O

5

O=C(N1)NC(C1C2CCCCC (N3CCN(CC4=C(C=CC=C 5)C5=CC=C4)CC3)=O)CS 2(=O)=O

6

ClC1=CC=C(C(N2CCCC( C(NCC3=CC(OC)=C(OCC 4=CC=CC=C4)C=C3)=O)C 2)=O)C=C1

7

O=S(C1=CC=CC=C1)(N2 CCN(C(CNC3=CC=CC(C( N4)=NC(CC)=CC4=O)=C3 )=O)CC2)=O CC1=CC=CC=C1C2=NN(

8

CCO)C(NC(NC3CCC(O)C C3)=O)=C2

9

O=C(C1=O)N(CC(C2=CC( CCC(N3)=O)=C3C=C2)=O )C(N1CCC4=CCCCC4)=O

10

ClC1=C(C2=CSC(NC(CN C3=CC(O)=CC=C3)=O)=C 2C(OCC)=O)C=CC=C1

11

O=C1C(CCC(NCC2=CC= C(CN3C(C=CC=C4)=C4C C3)C=C2)=O)C(C)NC(N1) =O Ketoconazole a

: The unit of MM-GBSA ∆Gbinding was kcal/mol; : IC50 values (µM) of compounds by human liver microsome assay;

b

27

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

c

: IC50 values (µM) of compounds by recombinant human CYP3A4 assay. The binding modes of three compounds with IC50 less than 10 µM at the

peripheral site of CYP3A4 are shown in Figure 6. All the three hits formed π-π stacking interactions with Phe219 and kept the hydrogen bonding interaction with Leu211. In addition, compounds 2 kept the hydrogen bonding interactions with both Asp214 and Lys209. Compounds 6 and 10 formed hydrogen bonding interactions with Asp214. These interactions contribute to better docking G-scores and MM-GBSA scores of these three compounds. It is worth noting that these binding modes were obtained by docking and warrant further validation by x-ray crystallography of CYP3A4 in complex with these ligands. However, the determination of the crystal structures of CYP3A4 complexes is laborious and costly and remains a challenge, especially for ligands bound at the peripheral site.21 The purpose of this study is to explore the stability of PGS at the peripheral site and test the hypothesis of discovering the new CYP3A4 inhibitors binding to the peripheral site. Our results have proven that it is possible to find the CYP3A4 inhibitors based on the peripheral site. The determination of crystal structures of CYP3A4 in complex with these inhibitors has not been conducted. Nevertheless, the compounds that we discovered provide the substance basis for future efforts by the experts in this field.

28

ACS Paragon Plus Environment

Page 28 of 37

Page 29 of 37

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

Journal of Chemical Information and Modeling

Figure 6. The binding mode of three hits. Panels A-C corresponded to compounds 2, 6, and 10. Although the inhibitory activities of these compounds were not as potent as the known CYP3A4 inhibitors such as ketoconazole, the three compounds with new scaffold may serve as effectors that could trigger a functionally important conformational transition to affect the catalytic metabolism of CYP3A4. In addition, these compounds could provide lead structures for modification in the discovery and development of drugs targeting CYP3A4.

4. Conclusions In the present study, we used a strategy that combined molecular docking with MD simulations of multiple paralleled systems to verify the stability of CYP3A4 peripheral site bound with PGS and discover new inhibitors of CYP3A4 from a commercial database based on the optimized peripheral site. Through the MD refinement and clustering analysis, it was found that though it is solvent-exposed, PGS can bind stably in enzyme periphery by contacts with several residues including Phe215, Phe219, Phe220 and Asp214. In addition, the stability of PGS bound in the 29

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

peripheral site was verified by the MD-derived self-docking. Based on these findings, the CYP3A4 peripheral site-based virtual screening model was firstly built and used for screening against the Enamine database. Three compounds were identified to have inhibitory activity towards CYP3A4 and one of them showed potent inhibitory activity with IC50 lower than 1 µM. To our knowledge, our study for the first time demonstrated that it is possible to discover CYP3A4 inhibitors based on the peripheral site other than the canonical site above heme. The discovered compounds could serve as the starting point for future drug design targeting CYP3A4.

30

ACS Paragon Plus Environment

Page 30 of 37

Page 31 of 37

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

Journal of Chemical Information and Modeling

Acknowledgements This work was supported by the National Natural Science Foundation of China (Grants 81373328 and 81373329) and the National Key Research and Development Program (Grant 2016YFA0502304). WuXi AppTec (Shanghai) Co., Ltd is gratefully acknowledged for human liver microsome and recombinant human CYP3A4 enzyme inhibition assays.

31

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

ASSOCIATED CONTENT Supporting Information Available: [The RMSD variations of ten systems (Figure S1), the alignments of CYP3A4 crystal structure and average conformation in eight systems (Figure S2). Energy contributions of key residues to the binding energy (Figure S3), the self-docking results of 1W0F (Table S1), and the binding free energy of 3A4 complexed with PGS (Table S2).] This material is available free of charge via the Internet at http://pubs.acs.org.

32

ACS Paragon Plus Environment

Page 32 of 37

Page 33 of 37

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

Journal of Chemical Information and Modeling

References 1.

Nebert, D. W.; Wikvall, K.; Miller, W. L. Human Cytochromes P450 in Health and Disease. Phil.

Trans. R. Soc. B 2013, 368, 20120431. 2.

Kirchmair, J.; Williamson, M. J.; Tyzack, J. D.; Tan, L.; Bond, P. J.; Bender, A.; Glen, R. C.

Computational Prediction of Metabolism: Sites, Products, SAR, P450 Enzyme Dynamics, and Mechanisms. J. Chem. Inf. Model. 2012, 52, 617-648. 3.

Ji, L.; Faponle, A. S.; Quesne, M. G.; Sainna, M. A.; Zhang, J.; Franke, A.; Kumar, D.; van Eldik,

R.; Liu, W.; de Visser. S. P. Drug Metabolism by Cytochrome P450 Enzymes: What Distinguishes the Pathways Leading to Substrate Hydroxylation Over Desaturation? Chem. Eur. J. 2015, 21, 9083-9092. 4.

Sevrioukova, I. F.; Poulos, T. L. Understanding the Mechanism of Cytochrome P450 3A4: Recent

Advances and Remaining Problems. Dalton Trans. 2013, 42, 3116-3126. 5.

Guengerich, F. P. Cytochrome P450 Enzymes in the Generation of Commercial Products. Nat. Rev.

Drug Discov. 2002, 1, 359-366. 6.

Denisov, I. G.; Grinkova, Y. V.; Baylon, J. L.; Tajkhorshid, E.; Sligar, S. G. Mechanism of

Drug-Drug Interactions Mediated by Human Cytochrome P450 CYP3A4 Monomer. Biochemistry 2015, 54, 2227-2239. 7.

Wang R. W.; Newton D. J.; Liu, N.; Atkins, W. M.; Lu, A. Y. Human Cytochrome P-450 3A4: In

Vitro Drug-Drug Interaction Patterns are Substrate-Dependent. Drug Metab. Dispos. 2000, 28, 360-366. 8.

Kempf, D. J.; Marsh, K. C.; Kumar, G.; Rodrigues, A. D.; Denissen, J. F.; McDonald, E.; Kukulka,

M. J.; Hsu, A.; Granneman, G. R.; Baroldi, P. A.; Sun, E.; Pizzuti, D.; Plattner, J. J.; Norbeck, D. W.; Leonard, J. M. Pharmacokinetic Enhancement of Inhibitors of the Human Immunodeficiency Virus Protease by Coadministration with Ritonavir. Antimicrob. Agents Chemother. 1997, 41, 654-660. 9.

Xu, L.; Liu, H.; Murray, B. P.; Callebaut, C.; Lee, M. S.; Hong, A.; Strickley, R. G.; Tsai, L. K.;

Stray, K. M.; Wang, Y.; Rhodes, G. R.; Desai, M. C. Cobicistat (GS-9350): A Potent and Selective Inhibitor of Human CYP3A as a Novel Pharmacoenhancer. ACS Med. Chem. Lett. 2010, 1, 209-213. 10. Zhou, S.; Chan, S. Y.; Goh, B. C.; Chan, E.; Duan, W.; Huang, M.; McLeod, H. L. Mechanism-Based Inhibition of Cytochrome P450 3A4 by Therapeutic Drugs. Clin. Pharmacokinet. 2005, 44, 279-304. 11. He, Y. A.; Roussel, F.; Halpert, J. R. Analysis of Homotropic and Heterotropic Cooperativity of Diazepam Oxidation by CYP3A4 Using Site-Directed Mutagenesis and Kinetic Modeling. Arch. Biochem. Biophys. 2003, 409, 92-101. 12. Atkins, W. M.; Wang, R. W.; Lu, A. Y. Allosteric Behavior in Cytochrome P450-Dependent in Vitro Drug-Drug Interactions: A Prospective Based on Conformational Dynamics. Chem. Res. Toxicol. 2001, 14, 338-347. 13. Davydov, D. R.; Halpert, J. R.; Renaud, J. P.; Hoa, G. H. B. Conformational Heterogeneity of Cytochrome P450 3A4 Revealed by High Pressure Spectroscopy. Biochem. Biophys. Res. Comm. 2003, 312, 121-130. 14. Müller, C. S.; Knehans, T.; Davydov, D. R.; Bounds, P. L.; von Mandach, U.; Halpert, J. R.; Caflisch, A.; Koppenol, W. H. Concurrent Cooperativity and Substrate Inhibition in the Epoxidation of Carbamazepine by Cytochrome P450 3A4 Active Site Mutants Inspired by Molecular Dynamics Simulations. Biochemistry 2015, 54, 711-721. 15. Isin, E. M.; Guengerich, F. P. Kinetics and Thermodynamics of Ligand Binding by Cytochrome 33

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

P450 3A4. J. Biol. Chem. 2006, 281, 9127-9136. 16. Isin, E. M.; Guengerich, F. P. Multiple Sequential Steps Involved in the Binding of Inhibitors to Cytochrome P450 3A4. J. Biol. Chem. 2007, 282, 6863-6874. 17. Davydov, D. R.; Rumfeldt, J. A.; Sineva, E. V.; Fernando, H.; Davydova, N. Y.; Halpert, J. R. Peripheral Ligand-Binding Site in Cytochrome P450 3A4 Located with Fluorescence Resonance Energy Transfer (FRET). J. Biol. Chem. 2012, 287, 6797-6809. 18. Sineva, E. V.; Rumfeldt, J. A.; Halpert, J. R.; Davydov, D. R. A Large-scale Allosteric Transition in Cytochrome P450 3A4 Revealed by Luminescence Resonance Energy Transfer (LRET). PloS one 2013, 8, e83898. 19. Williams, P. A.; Cosme, J.; Vinkovic, D. M.; Ward, A.; Angove, H. C.; Day, P. J.; Vonrhein, C.; Tickle, I. J.; Jhoti, H. Crystal Structures of Human Cytochrome P450 3A4 Bound to Metyrapone and Progesterone. Science 2004, 305, 683-686. 20. Protein Data Bank. http://www.rcsb.org/pdb/home/home.do (accessed Feb 18, 2017) 21. Sevrioukova, I. F.; Poulos, T. L. Anion-Dependent Stimulation of CYP3A4 Monooxygenase. Biochemistry 2015, 54, 4083-4096. 22. Yano, J. K.; Wester, M. R.; Schoch, G. A.; Griffin, K. J.; Stout, C. D.; Johnson, E. F. The Structure of Human Microsomal Cytochrome P450 3A4 Determined by X-ray Crystallography to 2.05-Å Resolution. J. Biol. Chem. 2004, 279, 38091-38094. 23. Guex, N.; Peitsch, M. C. SWISS-MODEL and the Swiss-Pdb Viewer: An Environment for Comparative Protein Modeling. Electrophoresis 1997, 18, 2714-2423. 24. Schrödinger Suite; Schrödinger, LLC, New York, NY, 2012. 25.

Li, H.; Robertson, A. D.; Jensen, J. H. Very Fast Empirical Prediction and Rationalization of

Protein pKa Values. Proteins: Struct., Funct., Bioinf. 2005, 61, 704-721. 26. Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham III, T. E.; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.; Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.; Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; York, D. M.; Kollman, P. A. AMBER12; University of California: San Francisco, 2012. 27. Gilson, M. K.; Sharp, K. A.; Honig, B. H. Calculating the Electrostatic Potential of Molecules in Solution: Method and Error Assessment. J. Comput. Chem. 1988, 9, 327-335. 28. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. 29. Li, W.; Ode, H.; Hoshino, T.; Liu, H.; Tang, Y.; Jiang, H. Reduced Catalytic Activity of P450 2A6 Mutants with Coumarin: A Computational Investigation. J. Chem. Theory Comput. 2009, 5, 1411-1420. 30. Harris, D. L.; Park, J. Y.; Gruenke, L.; Waskell, L. Theoretical Study of the Ligand-CYP2B4 Complexes: Effect of Structure on Binding Free Energies and Heme Spin State. Proteins: Struct., Funct., Bioinf. 2004, 55, 895-914. 31. Cai, J.; Li, J.; Zhang, J.; Ding, S.; Liu, G.; Li, W.; Tang, Y. Computational Insights into Inhibitory Mechanism of Azole Compounds against Human Aromatase. RSC Adv. 2015, 5, 90871-90880. 32. Rychaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular dynamics of n-Alkanes. J. Chem. Phys. 1977, 23, 327-341. 33. Essmann, U.; Perera, L.; Berkowitz, M. L; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle 34

ACS Paragon Plus Environment

Page 34 of 37

Page 35 of 37

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

Journal of Chemical Information and Modeling

Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593. 34. Feig, M.; Karanicolas, J.; Brooks, C. L. MMTSB Tool Set: Enhanced Sampling and Multiscale Modeling Methods for Applications in Structural Biology. J. Mol. Graphics Modell. 2004, 22, 377-395. 35. Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Donini, O. Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889-897. 36. Kuhn, B.; Kollman, P. A. Binding of a Diverse Set of Ligands to Avidin and Streptavidin: An Accurate Quantitative Prediction of Their Relative Affinities by a Combination of Molecular Mechanics and Continuum Solvent Models. J. Mol. Model. 2000, 43, 3786-3791. 37. Weiser, J.; Shenkin, P. S.; Still, W. C. Approximate Atomic Surfaces from Linear Combinations of Pairwise Overlaps (LCPO). J. Comput. Chem. 1999, 20, 217-230. 38. Li, J.; Du, H.; Wu, Z.; Su, H.; Liu, G.; Tang, Y.; Li, W. Interactions of Omeprazole-Based Analogues with Cytochrome P450 2C19: A Computational Study. Mol. BioSyst. 2016, 12, 1913-1921. 39. Lipinski, C. A.; Lombardo, F.; Dominy, B. W.; Feeney, P. J. Experimental and Computational Approaches to Estimate Solubility and Permeability in Drug Discovery and Development Settings. Adv. Drug Deliv. Rev. 2012, 64, 4-17. 40. Friesner, R. A.; Murphy, R. B.; Repasky, M. P.; Frye, L. L.; Greenwood, J. R.; Halgren, T. A.; Sanschagrin, P. C.; Mainz, D. T. Extra Precision Glide: Docking and Scoring Incorporating a Model of Hydrophobic Enclosure for Protein-Ligand Complexes. J. Med. Chem. 2006, 49, 6177-6196. 41. Shen, J.; Tan, C.; Zhang, Y.; Li, X.; Li, W.; Huang, J.; Shen, X.; Tang, Y. Discovery of Potent Ligands for Estrogen Receptor β by Structure-Based Virtual Screening. J. Med. Chem. 2010, 53, 5361-5365. 42. Midde, N. M.; Rahman, M. A.; Rathi, C.; Li, J.; Meibohm, B.; Li, W.; Kumar, S. Effect of Ethanol on the Metabolic Characteristics of HIV-1 Integrase Inhibitor Elvitegravir and Elvitegravir/Cobicistat with CYP3A: An Analysis Using a Newly Developed LC-MS/MS Method. PloS one 2016, 11, e0149225. 43. Ren, S.; Zeng, J.; Mei, Y.; Zhang, J. Z.; Yan, S. F.; Fei, J.; Chen, L. Discovery and Characterization of Novel, Potent, and Selective Cytochrome P450 2J2 Inhibitors. Drug Metab. Dispos. 2013, 41, 60-71. 44. Baylon, J. L.; Lenov, I. L.; Sligar, S. G.; Tajkhorshid, E. Characterizing the Membrane-Bound State of Cytochrome P450 3A4: Structure, Depth of Insertion, and Orientation. J. Am. Chem. Soc. 2013, 135, 8542-8551. 45. Lonsdale, R.; Rouse, S. L.; Sansom, M. S.; Mulholland, A. J. A Multiscale Approach to Modelling Drug Metabolism by Membrane-Bound Cytochrome P450 Enzymes. PLoS Comput. Biol. 2014, 10, e1003714. 46. Domanski, T. L.; He, Y. A.; Harlow, G. R.; Halpert, J. R. Dual Role of Human CytochromeP450 3A4 Residue Phe-304 in Substrate Specificity and Cooperativity. J. Pharm. Exp. Ther. 2000, 293, 585-591. 47. Domanski, T. L.; He, Y. A.; Khan, K. K.; Roussel, F.; Wang, Q.; Halpert, J. R. Phenylalanine and Tryptophan Scanning Mutagenesis of CYP3A4 Substrate Recognition Site Residues and Effect on Substrate Oxidation and Cooperativity. Biochemistry 2001, 40, 10150-10160. 48. Domanski, T. L.; Liu, J.; Harlow, G. R.; Halpert, J. R. Analysis of Four Residues within Substrate Recognition Site 4 of Human Cytochrome P450 3A4: Role in Steroid Hydroxylase Activity and α-Naphthoflavone Stimulation. Arch. Biochem. Biophys. 1998, 350, 223-232. 35

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

49. Fishelovitch, D.; Shaik, S.; Wolfson, H. J.; Nussinov, R. Theoretical Characterization of Substrate Access/Exit Channels in the Human Cytochrome P450 3A4 Enzyme: Involvement of Phenylalanine Residues in the Gating Mechanism. J. Phys. Chem. B 2009, 113, 13018-13025. 50. Krishnamoorthy, N.; Gajendrarao, P.; Thangapandian, S.; Lee, Y.; Lee, K. W. Probing Possible Egress Channels for Multiple Ligands in Human CYP3A4: A Molecular Modeling Study. J. Mol. Model. 2010, 16, 607-614. 51. Treuheit, N. A.; Redhair, M.; Kwon, H.; McClary, W. D.; Guttman, M.; Sumida, J. P.; Atkins, W. M. Membrane Interactions, Ligand-Dependent Dynamics, and Stability of Cytochrome P4503A4 in Lipid Nanodiscs. Biochemistry 2016, 55, 1058-1069. 52. Harlow, G. R.; Halpert, J. R. Analysis of Human Cytochrome P450 3A4 Cooperativity: Construction and Characterization of a Site-Directed Mutant that Displays Hyperbolic Steroid Hydroxylation Kinetics. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 6636-6641. 53. Harlow, G. R.; Halpert, J. R. Alanine-Scanning Mutagenesis of a Putative Substrate Recognition Site in Human Cytochrome P450 3A4 Role of Residues 210 and 211 in Flavonoid Activation and Substrate Specificity. J. Biol. Chem. 1997, 272, 5396-5402. 54.

Baell, J. B.; Holloway, G. A. New Substructure Filters for Removal of Pan Assay Interference

Compounds (PAINS) from Screening Libraries and For Their Exclusion in Bioassays. J. Med. Chem. 2010, 53, 2719-2740. 55.

Nettleton, D. O.; Einolf, H. J. Assessment of Cytochrome P450 Enzyme Inhibition and

Inactivation in Drug Discovery and Development. Curr. Top. Med. Chem. 2011, 11(4): 382-403.

36

ACS Paragon Plus Environment

Page 36 of 37

Page 37 of 37

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

Journal of Chemical Information and Modeling

Graphical Table of Contents

The stability of progesterone at the peripheral site of CYP3A4 was investigated by multiple-paralleled MD simulations. On the basis of the MD refined peripheral site, the new CYP3A4 ligands were successfully discovered.

37

ACS Paragon Plus Environment