Molecular docking gave clues on rational designing compounds with potential enhanced activities and understanding SAR. In the complex structure of TNF-α with SPD304, SPD304 bound to a pocket in the TNF-α dimer (He et al., 2005; Shen et al., 2014). **EJMC-1** was shown to bind with the same site (He et al., 2005; Shen et al., 2014). Both Glide and AutoDock Vina were used in the docking study. We first tested whether the binding pose of SPD304 can be reproduced. We have tried many times with different parameters, but was unable to get a binding pose that is close to that in the crystal structure (with the minimum RMSD up to 4 Å). We then used AutoDock Vina to dock SPD 304 to the TNF-α dimer and the lowest binding free energy conformation obtained was closed to its crystal conformation with a RMSD of 0.70 Å. Despite of the different binding conformations obtained for SPD304 by using two docking software, the top ranking conformations of **EJMC-1** were almost the same from the docking runs using both Glide and AutoDock Vina. These differences might due to the flexibility of SPD304, which adopted a U shape conformation, and the conformational sampling preference of the docking software. As there are no essential differences in the docking poses of the compounds other than SPD304, we used the Glide docking poses of these compounds to compare to SPD304 in the crystal structure. Compared to **EJMC-1**, **S10** had increased hydrophobic interaction with the Tyr59 residue (**Figures 4A,B**). In addition to the nonpolar interactions with TNF-α as in the case of SPD304, the scaffold of **EJMC-1** and **S10** provide further polar interactions, strengthening the specificity and activity (**Figure 4C**). As **EJMC-1** is smaller than that of SPD304 with unoccupied hydrophobic space in the pocket (**Figure 4C**), several analogs with larger substituted group of sulfonamide of 2-oxo-1,2-dihydrobenzo[cd]indole-6 sulfonamide were designed and docked to this site. Compound **4e,** with larger hydrophobic group size and additional H-bond donor, interacts favorably with TNF-α and might be more potent (**Figure 4D**). Based on the docking analysis, the designed analogs were purchased or synthesized for cell assay.

# Optimization of Compound S10 and Structure-Activity Analysis

As shown in the cell assay, the inhibition activity of **S10** increased about 2-fold than that of **EJMC-1**. The introduction of the naphthalene ring provides stronger hydrophobic interactions. Based on the docking analysis and increased activity of **S10**, we try to: (1) Keep naphthalene ring, changed Nsubstituted groups of dihydrobenzo[cd]indole, (2) Keep N-H of dihydrobenzo[cd]indole, optimize the hydrophobic R group, (3) optimize both N-substituted groups of dihydrobenzo[cd]indole and the hydrophobic R group (**Figure 5**). Seven commercially available analogs of **S10** were purchased for testing (**Table 1**). The SPECS ID of these seven **S10** analogues were listed in Table S3. We further synthesized seven new compounds in three steps from 1,8-Naphthalic anhydride through conventional reactions (**Scheme 1** and **Figure 5**). All compounds passed the PAINS (pan assay interference compounds) remover, which filters out compounds that appear as frequent hitters (promiscuous compounds) in many biochemical high throughput screens (Baell and Holloway, 2010).

All the compounds were tested using the TNF-α induced NF-κB reporter assay. The structures and activities were listed in **Table 1**. For **S10**, methyl or ethyl group substitution on the amide of 2-oxo-1,2-dihydrobenzo[cd]indole-6-sulfonamide had no obviously enhanced inhibition (**S21** and **S22**), and the α or β substitution of the naphthyl group did not affect the inhibition (**S23, S24,** and **S25**). The size of the N-substituted groups of sulfonamide was important for inhibition activity (**EJMC-1, S10,** and **S27**). The flexibility and aromaticity of the N-substituted two-ring group of sulfonamide played dominant role, too rigid or too flexible dramatically reduced the activity (**4g**, **4d**). The fact that N-(5-aminonaphthalen-1-yl) and N-(3-aminonaphthalen-1 yl) group substituted compounds lost functions might be caused by the conformation change due to additional amino group on the naphthalene ring. Introducing heterocycle significantly increases the inhibition activity (**4e and 4f**). The N-(1H-indol-6-yl) substituted sulfonamides (**4e**) were 6-fold more potent than **S10**, even better than SPD304 (**Table 1**, **Figure 3B**). Though **S10** and **4e** had similar size of substitution group on sulfonamide, **4e** shown better inhibition activity than **S10** might due to the additional H-bond that **4e** forms with the backbone carbonyl of Gly121 (**Figure 4C**). Meanwhile, the indolyl group of **4e** was also deeper in the binding pocket than that of naphthyl group on **S10** (**Figure 4D**).

#### CONCLUSION

We have optimized a previously reported TNF-α inhibitor **EJMC-1** using similarity-based VS and rational design. An analog of **EJMC-1**, **S10** was found with 2-fold TNF-α increased inhibition activity. Based on the structures of **EJMC-1**, **S10,** and their interactions with TNF-α, we designed derivatives of 2-oxo-1,2-dihydrobenzo[cd]indole-6-sulfonamide. Several commercially available ones were purchased and seven new compounds were synthesized for SAR study. After two rounds of design, we obtained **4e** with an IC50 of 3.0 ± 0.8µM, which is one of the most potent TNF-α small molecule inhibitors reported so far. Compound **4e** provides a good starting point for developing more potent TNF-α small molecule inhibitors.

#### REFERENCES

# AUTHOR CONTRIBUTIONS

LL and YL designed and guided this study; XD designed the research, performed molecular docking and similarity search, and conducted the chemical synthesis; XZ performed the cell assay; HL and QS participated in the cell assay; BT performed the SPR binding assay; XD, LL, and YL analyzed the data and wrote the manuscript with input from all authors; XD and XZ have equal contribution of this work.

#### ACKNOWLEDGMENTS

This work was supported in part by the Ministry of Science and Technology of China (2016YFA0502303, 2015CB910300) and the National Natural Science Foundation of China (21633001, 21573012).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2018.00098/full#supplementary-material

discovery of a tumor necrosis factor-α antagonist. J. Am. Chem. Soc. 135, 11990–11995. doi: 10.1021/ja405106u

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Deng, Zhang, Tang, Liu, Shen, Liu and Lai. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Insights Into the Bifunctional Aphidicolan-16-ß-ol Synthase Through Rapid Biomolecular Modeling Approaches

#### Max Hirte, Nicolas Meese, Michael Mertz, Monika Fuchs\* and Thomas B. Brück\*

Werner Siemens Chair of Synthetic Biotechnology, Department of Chemistry, Technical University of Munich, Munich, Germany

#### Edited by:

Daniela Schuster, Paracelsus Medizinische Privatuniversität, Salzburg, Austria

#### Reviewed by:

Victor Guallar Guallar, Barcelona Supercomputing Center, Spain Dharmendra Kumar Yadav, Gachon University of Medicine and Science, South Korea Arnout Voet, KU Leuven, Belgium

#### \*Correspondence:

Monika Fuchs<EMAIL_ADDRESS>Thomas B. Brück<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Medicinal and Pharmaceutical Chemistry, a section of the journal Frontiers in Chemistry

Received: 10 January 2018 Accepted: 20 March 2018 Published: 10 April 2018

#### Citation:

Hirte M, Meese N, Mertz M, Fuchs M and Brück TB (2018) Insights Into the Bifunctional Aphidicolan-16-ß-ol Synthase Through Rapid Biomolecular Modeling Approaches. Front. Chem. 6:101. doi: 10.3389/fchem.2018.00101 Diterpene synthases catalyze complex, multi-step C-C coupling reactions thereby converting the universal, aliphatic precursor geranylgeranyl diphosphate into diverse olefinic macrocylces that form the basis for the structural diversity of the diterpene natural product family. Since catalytically relevant crystal structures of diterpene synthases are scarce, homology based biomolecular modeling techniques offer an alternative route to study the enzyme's reaction mechanism. However, precise identification of catalytically relevant amino acids is challenging since these models require careful preparation and refinement techniques prior to substrate docking studies. Targeted amino acid substitutions in this protein class can initiate premature quenching of the carbocation centered reaction cascade. The structural characterization of those alternative cyclization products allows for elucidation of the cyclization reaction cascade and provides a new source for complex macrocyclic synthons. In this study, new insights into structure and function of the fungal, bifunctional Aphidicolan-16-ß-ol synthase were achieved using a simplified biomolecular modeling strategy. The applied refinement methodologies could rapidly generate a reliable protein-ligand complex, which provides for an accurate in silico identification of catalytically relevant amino acids. Guided by our modeling data, ACS mutations lead to the identification of the catalytically relevant ACS amino acid network I626, T657, Y658, A786, F789, and Y923. Moreover, the ACS amino acid substitutions Y658L and D661A resulted in a premature termination of the cyclization reaction cascade en-route from syn-copalyl diphosphate to Aphidicolan-16-ß-ol. Both ACS mutants generated the diterpene macrocycle syn-copalol and a minor, non-hydroxylated labdane related diterpene, respectively. Our biomolecular modeling and mutational studies suggest that the ACS substrate cyclization occurs in a spatially restricted location of the enzyme's active site and that the geranylgeranyl diphosphate derived pyrophosphate moiety remains in the ACS active site thereby directing the cyclization process. Our cumulative data confirm that amino acids constituting the G-loop of diterpene synthases are involved in the open to the closed, catalytically active enzyme conformation. This study demonstrates that a simple and rapid biomolecular modeling procedure can predict catalytically relevant amino acids. The approach reduces computational and experimental screening efforts for diterpene synthase structure-function analyses.

Keywords: homology modeling, aphidicolin, diterpene, diterpene synthase, homology model refinement

# INTRODUCTION

With more than 50,000 different molecules known to date terpenes are the greatest natural occurring product family found in organisms from bacteria to fungi, mammals, and plants. They are all derived from the isoprene units' dimethylallyl diphosphate and isopentenyl diphosphate. Condensation reactions of this molecules lead to the formation of different length phosphorylated linear terpenes, serving as substrate for terpene synthases. This enzyme family carry out highly stereo complex C-C coupling reactions, resulting in structurally complex macrocycles that contribute to the structural and functional diversity of terpenes (Christianson, 2017). Diterpenes are derived from the linear aliphatic precursor geranylgeranyl diphosphate (GGDP) being cyclized by diterpene synthases. More specifically, diterpene synthases are classified into class I and class II enzymes based on the structural presence of the conserved motifs DDXD or DDXXD/E and NSE/DTE, respectively. While class II reactions perform a protonation initiated cyclization reaction to generate phosphorylated bicyclic structures, class I reactions are initiated by hydrolyses of the GGDP pyrophosphate moiety that is coordinated by a Mg2+-triad thereby generating mono- or poly-cyclic structures.

The natural product Aphidicolin, initially isolated from the fungus Cephalosporium aphidicola, is a hydroxylated, tetracyclic diterpenoid that exhibits a broad range of biological activities and applications (Brundret et al., 1972; Dalziel et al., 1973). More specifically, it is a potent inhibitor of the eukaryotic DNA α-polymerase with a commercial application as a cell synchronization agent. The compound is in pharmaceutical development due anti-tumor, anti-viral, and anti-leishmanial activity (Ikegami et al., 1978; Pedrali-Noy et al., 1980; Kayser et al., 2001; Edwards et al., 2013; Starczewska et al., 2016). Recently, other organisms including the fungus Nigrospora sphaerica and the pathogenic fungus Phoma betae have been identified as natural Aphidicolin producers. Current data suggests that Aphidicolin biosynthesis is exclusive to fungal metabolism and that natural sources for Aphidicolin are limited (Starratt and Loschiavo, 1974; Fujii et al., 2011; Lopes and Pupo, 2011). Nevertheless, elucidation of the responsible Aphidicolin biosynthetic gene cluster in P. betae allowed for the identification of a bifunctional diterpene synthase that contains both a functional class I and class II domain (Oikawa et al., 2001). The Aphidicolan-16-ß-ol synthase (ACS) generates the stereo-chemically demanding Aphidicolan-16-ß-ol (AD)—core structure of Aphidicolin—structure via a two-step reaction as depicted in **Figure 1** (Oikawa et al., 2002).

Initially, GGDP is rearranged in the class II active site cleft by protonation to the bicyclic syn-copalyl diphosphate (syn-CDP). Subsequently, syn-CDP is elaborated to AD in the class I active site (Adams and Bu'Lock, 1975; Oikawa et al., 2002). As depicted in **Figure 2** the cyclization mechanism in the class I active site, initiated by the hydrolysis of the pyrophosphate group, results

in 8-ß-pimaradienyl carbocation formation. A subsequent attack of the vinyl group, bridging the C ring, directly undergoes a Wagner-Meerwein rearrangement and results in the formation of the aphidicolenyl carbocation. Eventually, this cation is quenched by water thereby generating AD.

Terpene cyclization mechanisms are conventionally elucidated by radio labeling of protons and carbons (Dickschat, 2017). This substrate specific labeling provides for identification of unusual hydride shifts and rearrangements. Alternatively, the enzyme's cyclization mechanisms can be probed by altering amino acids, trying to terminate the reaction cascade at a specific transition state (Morrone et al., 2008; Janke et al., 2014; Schrepfer et al., 2016; Jia et al., 2017). Therefore, random mutagenesis can be performed but the screening effort for this methodology is elaborate without an efficient high throughput screening options (Lauchli et al., 2013). Biomolecular modeling allows for the rational identification and in silico modulation of amino acid networks that are involved in complex reaction cascades (Pemberton et al., 2015; Schrepfer et al., 2016; Christianson, 2017; Escorcia et al., 2018). This methodology provides for a knowledge based approach of enzyme mutagenesis and screening. Nevertheless, a particular challenge for this strategy is based on the missing structural information for most terpene synthases. However, as their structural elements and domains are highly conserved (Christianson, 2017), homology

**Abbreviations:** ACS, Aphidicolan-16-β-ol synthase; AD, Aphidicolan-16-β-ol; GGDP, geranylgeranyl diphosphate; LRS, labdane related diterpene synthase; syn-CDP, syn-copalyl diphosphate.

modeling is a potential route to identify catalytically relevant amino acids despite the low primary sequence identities in this enzyme family (Xu and Li, 2003). Unfortunately, most available crystal structures of terpene synthase are deposited in the open apo-enzyme configuration that is catalytically inactive. This open enzyme conformation presents an additional obstacle when catalytically relevant amino acids have to be identified in silico. At present, only two diterpene synthase structures have been reported in the closed, catalytically active form (Liu et al., 2014; Serrano-Posada et al., 2015). Therefore, automated homology modeling approaches will almost always result in catalytically non-relevant open enzyme configuration. Moreover, while prediction tools can place large cofactors (i.e., FAD, NADH, Heme) correctly in the apo-protein framework, ligand-metal interactions are difficult to predict because of the multiple coordination geometries and the lack of sufficiently accurate force field parameters (Khandelwal et al., 2005). Hence, structure function predictions that depend on the interplay between the amino acids of the protein framework with small metal ions cannot be conducted solely by application of automated software tools. In this context, a rational combination of structural information by superposition and extraction of cofactors is performed to prepare the protein structure for docking studies. Nevertheless, this approach often neglects reliable positioning of the cofactor coordinating amino acids. Additionally, falsely predicted positioning of amino acid side chains in the active site cleft can lead to invalid interpretation of a homology model based protein-ligand complex. To improve this situation, this study elucidated rapid and simple methodologies to refine diterpene homology models for docking studies thereby allowing for reliable structure-function predictions. In this context, an ACS class I homology model of the α-domain was predicted from the primary sequence. Subsequently, these models were compared to catalytically relevant closed terpene synthases structures. The location of metals was refined and fitted against specifically selected structural templates and multiple docking studies were carried out and validated. Our in silico results were experimentally evaluated by ACS mutagenesis studies. This lead to an identification of essential amino acid residue sidechains that are necessary for retaining the enzymes activity. Additionally, we detected amino acid substitutions that abort the catalytic reaction cascade en- route from syn-CDP to AD. Structural analyses and elucidation of these compounds

revealed the formation of syn-copalol and a labdane related, non-hydroxylated diterpene by the ACS mutants Y658L and D661A. Our approach of a protein homology model based structure function analysis can be easily adapted for other terpene synthases. This methodology allows for rapid and simple analysis of the catalytically relevant amino acid network that help studying complex reaction cascades and developing new biocatalysts.

### MATERIALS AND METHODS

#### Materials and Chemicals

All genes used were synthesized by Life technologies GmbH and the codon usage was optimized for E. coli if not stated otherwise. Primers were obtained from Eurofins Genomics GmbH. Strains and plasmids were obtained from Merck KGaA. All chemicals used were obtained at highest purity from Roth chemicals or Applichem GmbH. Enzymes were purchased from Thermo Fisher Scientific.

#### Software and Web-Tools

RaptorX was applied for homology modeling studies (http:// raptorx.uchicago.edu; Källberg et al., 2012). The initial predicted structure was analyzed and further modified in the environment of UCSF Chimera software package (Pettersen et al., 2004; http:// www.cgl.ucsf.edu/chimera). Comparative modeling by spatial restraints was performed by MODELLER (Eswar et al., 2006), and all substrate docking studies performed by AutoDock Vina (Trott and Olson, 2010; http://vina.scripps.edu). Chemical structures were drawn by PerkinElmer ChemBioDraw Ultra (http:// www.cambridgesoft.com). For ligand preparation the Avogadro (Hanwell et al., 2012; https://avogadro.cc/) software package was used. A syn-CDP toppar stream file was generated by CHARMM General Force Field program version 1.0.0 for use with CGenFF version 3.0.1 (https://cgenff.paramchem.org; Vanommeslaeghe et al., 2010, 2012; Vanommeslaeghe and MacKerell, 2012). Two ns molecular dynamic studies of the docked ACS model B in a water sphere have been performed under CHARMM general force field by NAMD (Phillips et al., 2005; http://www.ks.uiuc. edu/Research/namd/). NAMD was developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana-Champaign. For high resolution pictures the protein was prepared by Visual Molecular Dynamics (http://www.ks. uiuc.edu/Research/vmd/; Humphrey et al., 1996) and rendered by Tachyon implemented in the VMD software package (Stone, 1998).

#### Docking

Ligand structures were downloaded from https://pubchem.ncbi. nlm.nih.gov/ available and geometrically optimized by 500 steps of steepest descent under MMFF94 force field parameters included in Avogadro. Protein structures were prepared by Dock Prep, which is part of the Chimera software environment. The AMBER force field (AMBERff14SB) was applied to the receptor while Gasteiger charges were added to the ligand and co-factors. As recently reported, docking can be improved by assigning partial charges to metal ions (Hu and Shelver, 2003). In this context, Mg-ion charges were set to +1. Syn-CDP charge was set to −3. Docking was performed by AutoDock Vina using standard parameters. Docking poses were chosen based on a structural comparison to the pyrophosphate group that is co-crystallized in pdb 5A0J (see Figure S1B). The chosen pose was furthermore validated by re-dock approaches. Therefore, the predicted syn-CDP pose was de novo geometrically optimized by 500 steps of steepest descent under MMFF94 force field parameters included in Avogadro software environment prior to docking repetition (see Figure S1A).

### Model Generation

An initial homology model of the ACS α-domain was predicted by RaptorX starting from the amino acid 565. A model based on the pdb crystal structure 5A0J, referring to a labdane related diterpene synthase, was manually selected for further structure function analyses. In order to prepare the model for docking studies, the coordinating Mg2+-ion triad and water molecules were implemented in the structure by different methods. Model A was generated by structural alignment to 5A0J. Cofactor positions were transferred from the structure template to Model A without any further adjustment prior to docking studies. Model B was created by MODELLER implemented in the Chimera software environment using the 5A0J as template structure. In this model hetero atoms and water molecules in the structure environment were computationally implemented. The pyrophosphate group was removed prior to docking with syn-CDP. Model C was prepared analogously to Model B but prior to refinement by MODELLER, syn-CDP was docked into the template structure 5A0J.

#### Model Validation

The protein ligand complex of Model B was validated by molecular dynamics studies. Therefore, syn-CDP was initially extracted from Model B and parameterized by CHARMM General Force Field program version 1.0.0 for use with CGenFF version 3.0.1. VMD was used to parameterize the protein and for merging ligand and protein. Subsequently, a water sphere was added around the protein-ligand complex. Two nanoseconds of molecular dynamic studies under CHARMM General Force Field was applied to the protein complex by NAMD. The calculated rmsd of the generated frames was plotted over time (Figure S2). A constant rmsd value was chosen as the criteria for an equilibrated protein-ligand complex. The last frame obtained was compared to the initial model B (Figure S3).

#### Plasmids for Diterpene Production

For all cloning procedures E. coli HMS 174 (DE3) was used. Clones were cultivated at 37◦C in Luria-Bertani (LB) medium. Chloramphenicol (34 µg/L) and Kanamycin (50 µg/L) were added as required. For efficient production of the diterpene AD, E. coli's internal 1-deoxy-xylulose-5 phosphate pathway flux was increased by overexpression of deoxy-xylulose 5 phosphate synthase (dxs: GenBank: YP001461602.1), isopentenyldiphosphate delta isomerase (idi: GenBank: AAC32208.1), and further extended by expressing geranylgeranyl diphosphate synthase (crtE: GenBank: KPA04564.1) and Aphidicolan-16-ß-ol synthase (acs: GenBank: AB049075.1). Therefore, dxs and acs were amplified from original sources by PCR. Polycistronic operons (**Table 1**) were constructed by BioBrick cloning standard (Shetty et al., 2008).

Site directed mutations of acs were generated by PCR. Forward primers were designed exhibiting the respective mutation at the 5′ end while the corresponding reverse primers were phosphorylated at 5′ end (Table S1). PCR products were ligated by T4 Ligase prior to transformation. All amino acid exchanges were confirmed by sequencing.

### Production of Diterpenes

All diterpene production experiments were performed in E. coli BL 21 (DE3). To investigate the product outcome of ACS mutants, pACYC acs plasmids were co-transformed with pAX dic. Cultivation was performed in minimal media supplemented with 6 g/L yeast extract and 30 g/L glycerol at 25◦C. After 60 h the culture was extracted with a mixture of hexane, ethanol and ethyl acetate (1:1:1) (v/v/v) for 1 h. The extract was centrifuged at 10,000 g for 2 min. The upper, organic phase was directly analyzed for diterpene products via GC-MS.

### Diterpene Analytics

GC-MS analyses of diterpenes was performed by a Trace GC Ultra with DSQII (Thermo Fisher Scientific). Therefore, 1 µL sample was loaded (Split 1/10) by TriPlus AS onto a SGE BPX5 column (30 m, I.D 0.25 mm, Film 0.25µm). The initial column temperature was set to 160◦C and maintained for 5 min before a temperature gradient at 8◦C/min up to 320◦C was applied. The final temperature was kept for additional 3 min. MS data were recorded at 70 eV (EI) and m/z (rel. intensity in %) as total ion current (TIC). The recorded m/z range was in between 50 to 650.

NMR spectra were recorded in CDCl3 with an Avance III 500 MHz (Bruker) at 300 K. 1H NMR chemical shifts are given in ppm relative to CDCl3 (δ = 7.26 ppm). The 2D experiments

(HSQC) were performed using standard Bruker pulse sequences and parameters.

# RESULTS AND DISCUSSION

#### Homology Model Refinement

The steady increase in published protein crystal structures provides for an accelerated improvement of computational homology prediction. Especially due to the high structurally conservation of the terpene synthase enzyme families, biomolecular tools can predict structures solely based on the amino acid sequence. In this context, structure prediction of the bifunctional ACS was performed to analyze the highly complex conversion of GGDP via syn-CDP to the tetracyclic AD which is the core structure of the cytostatic compound Aphidicolin. ACS belongs to the diterpene synthase family and we identified three highly structurally conserved domains. The initial conversion from the universal diterpene precursor GGDP to syn-CDP occurs in class II active site, located between the ACS ß- and γ-domain. The subsequent syn-CDP cyclization to AD is then conducted in the class I active site that is positioned in the middle of an α-helical bundle forming the ACS α-domain. Notably, the fungal ACS is structurally highly similar to the previously crystallized plant diterpene synthases Abietadiene (pdb: 3S9V) and Taxadiene synthase (pdb: 3P5R), respectively. Homology prediction based on the full ACS sequence took those two structures into account, but for both, crystals could only be achieved in N-terminal truncated forms. Furthermore, these crystal structures have only been solved in an open conformation that is catalytically inactive. In order to circumvent the consideration of these catalytically inactive templates, only the ACS α-domain sequence was used for homology prediction. A model based on the labdane related diterpene synthases (LRS) (pdb: 5A0J), which is provided in a catalytically active holo-complex (Serrano-Posada et al., 2015), was selected for ACS homology refinement. The structural superposition of Abietadiene (pdb: 3S9V), LRS (pdb: 5A0J), and the ACS model, as depicted in **Figure 3**, explicitly demonstrates that there is a better fit between the ACS model and the LRS crystal structure. While the structural fit between LRS and the ACS model is visually well apparent, we have not calculated an rmsd value qualifier as structural domains that do not constitute the active site region are highly variable.

Co-crystallized cofactors (Mg2+-ion triad) and waters, both provided in the LRS structure, are also involved in the ACS reaction en-route from syn-CDP to AD. Therefore, we differentially adapted both, the positions of the Mg2+ ions and waters into the ACS models that resulted in the generation of three ACS models (A–C). Model A was prepared by adaptation of cofactor positions from the template structure LRS after structural alignment. Initial evaluation of this model indicated that this un-refined modeling method results in clashes of cofactors positions with amino acids side chains. Generally, in homology prediction the active site's cavity is not reserved for the substrate or cofactors specifically. Therefore, we presume that amino acid sidechains occupy this free space due to applied energy minimization optimizations. This is demonstrated in

our docking studies of model A, where ACS amino acid Y658 is preventing syn-CDP to completely access the active site cavity. With the MODELLER package, which is based on comparative protein structure modeling by spatial restraints, a protein structure can be refined based on a template structure. Additionally, hetero-atoms and water molecules can be included directly in the model refinement. This refinement methodology applied to our initial model structure lead to the generation of ACS model B. This model B computationally included the three Mg2+-ions, a pyrophosphate group (conventionally derived by Mg2+ based hydrolysis of the phosphorylated substrate [syn-CDP] substrate) and water molecules directly as they are all present in the LRS template structure. Model B provides reliable positioning of the conserved amino acids that constitute the class I diterpene synthase signature DDXXD/E and NSE/DTE motifs in relation to the adapted Mg2+- ions, water and pyrophosphate moieties, respectively. Subsequently, we removed the pyrophosphate group from the model B structure to enable docking with the native syn-CDP substrate. Our docking data indicated that in Model B syn-CDP can completely access the active site's cavity. A specific syn-CDP conformation was selected pointing toward the ACS G-Helix, as this flexible helix is proposed to be involved in terpene cyclization reactions (Yoshikuni et al., 2006; Baer et al., 2014; Jia et al., 2017). This docking pose was validated by multiple re-docking approaches (Figure S1A). Additionally, we validated the pose while the position of the pyrophosphate moiety was compared to the pyrophosphate group co-crystallized in LRS (Figure S1B). Finally, a third approach for structure-function analyses was performed by docking syn-CDP into LRS prior to ACS refinement with MODELLER. Again, a syn-CDP conformation was chosen with close proximity toward the G-Helix. On the basis of this LRS holo-protein complex, an ACS holo-complex model C was generated. This method provided for a protein model that was refined around the substrate and cofactors. This methodology also provided for a precise specification of amino acids involved in the AD cyclization reaction. For all three models amino acids located within a five ´Å vicinity to the docked substrate syn-CDP (thereby neglecting the pyrophosphate moiety) were analyzed by mutational studies to elucidate their catalytic relevance (see **Figure 4**, Table S2).

## Mutational Validation of Catalytically Relevant ACS Amino Acids

Due to their stereo-chemical diversity, natural diterpene scaffolds are attractive research leads. The enormous stereo-chemical demand of diterpene macrocycles renders them difficult to access via total chemical synthesis approaches. Therefore, biosynthetic routes to generate these complex structures are currently an intense research focus (Dickschat, 2016; Bian et al., 2017; Jones, 2017). The ability to access new diterpene macrocycles via selective alteration of amino acids in diterpene synthases provides for a highly varied accessible chemical space. For the class I cyclooctat-9-en-7-ol synthase, which naturally generates a tricyclic fusicoccin type diterpene, amino acid mutations in the vicinity of the active site lead to intermittent abortion of the reaction cascade. Hence, alternative macrocyclic structures, such as the bicyclic dolabellane and the monocyclic cembrane, could be generated thereby elucidating the reaction cascade (Görner et al., 2013; Janke et al., 2014). In this study, insights into the class I reaction of the ACS were achieved by mutational studies. In that respect, we intended to quench the reaction from syn-CDP to AD at previously proposed transitional states (Adams and Bu'Lock, 1975; Oikawa et al., 2002). Based on the proposed ACS transitional states we presume that syn-labdatriene and syn-copalol (termination product of the syn-copalyl carbocation), stereoisomers of syn-pimaradiene (termination products of the pimaradienyl carbocation), or aphidicolene and stemodene (termination products of the aphidicolenyl carbocation) are potential abortion products (see **Figure 5**).

For an intermittent abortion of the reaction cascade from syn-CDP to AD, we have selected amino acids within a range of five ´Å to the docked ligands as prime targets for mutagenesis (see **Figure 4**). Preliminary studies revealed that sidechain substitutions encompassing amino acids exchanges that inherently change physico-chemical properties frequently resulted in inactive enzyme variants (Janke et al., 2014; Schrepfer et al., 2016). In this context, we focused on changing the size of the respective amino acid sidechain thereby trying to preserve physico-chemical characteristics. Alternatively, we chose amino acid side chain substitutions that would replace polar groups with similar size amino acids (Table S2).

ACS syn-CDP docking results pointed toward a strong interaction between the decalin core and surrounding hydrophobic sidechains. However, as the decalin structure of syn-CDP remains untouched in further cyclization steps most of the implemented mutations near this particular moiety resulted in inactive (I626A, Y923L, F789L) or wildtype activity variants (F629L, Y658F, C831G, C831T, T920G, Y923F). Based on our modeling results, we also identified specific amino acids located in the ACS G-helix that in other studies have been proposed to be of catalytic relevance (Baer et al., 2014;

Jia et al., 2017). While mutational changes in the G-Helix of Kaurene synthase like diterpene synthases resulted in alternative product profiles (Jia et al., 2017), our analogous approaches with ACS only provided inactive (A786L, F789L) or wildtype active (A786G, F789Y) variants. Nevertheless, our results support previous findings that propose the G-Helix as an essential flexible motif which is involved in the catalytically relevant structural change from the open to the closed enzyme configuration (Baer et al., 2014).

Only the substitution of ACS Y658L and D661A provided for a varied product outcome. In addition to amino acids that constitute the DXXDD/E and NSE/DTE signature motifs that are responsible for Mg2+-ion coordination, our combined in silico and experimental study identified only a few amino acids (see **Figure 6**, colored in pink) capable to terminate activity. Our successful mutations (D661A, Y658L) indicated that the unusual cyclization from syn-CDP to AD proceeds in a spatially restricted area of the active site's cleft. Additionally, our data suggests that the pyrophosphate group remains in the active site and coordinates the reaction cascade. This is in accordance to the recently postulated Taxadiene synthase reaction mechanism (Schrepfer et al., 2016).

#### ACS Mutants D661A and Y658L

GC-MS analyses of the ACS mutants Y658L and D661A revealed that this mutations lead to the formation of two unknown diterpene products (see **Figure 7**). In contrast to the native AD, which had a GC retention time of 17.67 min, these new diterpenes had a retention time of 12.79 and 13.46 min, respectively. The latter product with a retention time of 13.46 min, showed a total mass of 290 m/z. Comparison of the MS spectral data suggests that this was a hydroxylated diterpene with a similar structure to syn-copalol (Hoshino et al., 2011). Subsequently, this compound was isolated and structural characterized by NMR (Figures S4, S5). The results are in accordance to previous spectral data for syn-copalol (Yee and Coates, 1992). One plausible explanation for syn-copalol formation is the quenching of the syn-copalyl carbocation intermediate by water in the active site of the enzyme. The other diterpene product with a retention time of 12.79 min had a total mass of 272 m/z indicating that this structure was not-hydroxylated. While we expected the formation of syn-labda-8(17),12E,14-triene, comparison with published MS-spectra revealed significant differences (Morrone et al., 2011). Unfortunately, due to the low amounts produced and purification issues for this highly hydrophobic compound, we could not conduct NMR analysis. However, we presume that this compound is also originated from the syn-copalyl carbocation and that a labdane related diterpene with high structural similarity to syn-labda-8(17),12E,14-triene was generated by the ACS mutants. The newly generated diterpenes are of great interest as copalol derivatives display various biological activities analogous to aphidicolin (Hanson, 2015).

The structural changes (D661A and Y658L) still allowed syn-CDP binding in the active site with subsequent hydrolyses of the pyrophosphate group. The syn-copalyl carbocation was then quenched either by water (release of syn-copalol) or an amino acid side chain (release of non-hydroxylated diterpene). Furthermore, as we did not find other substitution that stopped cyclization at the proposed transitional states and as we could not even detect changes in the byproduct formation of the active mutants, we presume that the ACS cyclization occurs in a spatially restricted area and that the pyrophosphate group remains in the active site, which is in accordance to recent reports (Schrepfer et al., 2016).

Former diterpene centered production processes were limited by low target compound yields. However, optimization of recombinant diterpene production hosts has extensively progressed to provide gram per liter yields (Ajikumar et al., 2010; Schalk et al., 2012). Today, access to novel diterpene lead structures is limited by the effective identification of relevant enzyme systems from large scale genome sequencing projects. Therefore, rational alteration of known terpene synthase product profiles by using a combination of in silico prediction and knowledge based mutagenesis studies can allow for a more rapid and targeted expansion of the desired chemical space.

# CONCLUSION

A model of ACS synthase was computed that required the application of various methods for model refinement to improve the quality of in silico structure function analysis. A model of the

acid substitutions Y658L and D661A in the vicinity of the ACS active site lead to formation of the alternative cyclization products syn-copalol and a minor labdane related diterpene. Formation of these products were delineated by quenching of

the syn-copalyl carbocation en-route to AD. Additional mutants leading to inactive enzyme variants (A786L, F789L) provided insights into catalytically relevant amino acid residues within the G-Helix. The cumulative in-silico and experimental data suggests that amino acids constituting the G-loop motif of class I terpene cyclases are involved in the transformation of the open to the closed, catalytically active enzyme conformation. Moreover, as we only obtained a limited number of alternative cyclization products in our mutational screens, we presume that AD formation occurs in a rather confined location of the ACS active site. With respect to our biomolecular modeling approaches, we demonstrated that application of simple and rapid computational methodologies can be employed for prediction and structure function analyses of class I diterpene

# AUTHOR CONTRIBUTIONS

TB and MF supervised this study. MH initiated this study and performed virtual modeling and docking studies. NM and

#### REFERENCES

synthases.

MM conducted mutagenesis experiments and screening under supervision of MH and MF. Data was analyzed by MH, MF, NM, MM, and TB. All figures were created by MH. All authors verified the data, contributed to the manuscript, and approved the final version.

#### ACKNOWLEDGMENTS

MH, MF, and TB would like to acknowledge the financial support of the German ministry for Education and Research (BMBF) with the grant number 031A305A. TB gratefully acknowledges funding by the Werner Siemens foundation for establishing the field of Synthetic Biotechnology at the Technical University of Munich (TUM).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2018.00101/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Hirte, Meese, Mertz, Fuchs and Brück. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# How to Achieve Better Results Using PASS-Based Virtual Screening: Case Study for Kinase Inhibitors

Pavel V. Pogodin1 , Alexey A. Lagunin1,2, Anastasia V. Rudik 1 , Dmitry A. Filimonov 1 , Dmitry S. Druzhilovskiy 1 , Mark C. Nicklaus 3 and Vladimir V. Poroikov 1 \*

1 Department of Bioinformatics, Institute of Biomedical Chemistry, Moscow, Russia, 2 Department of Bioinformatics, Medical-Biological Department, Pirogov Russian National Research Medical University, Moscow, Russia, 3 Computer-Aided Drug Design Group, Chemical Biology Laboratory, Center for Cancer Research, National Cancer Institute, NIH, NCI-Frederick, Frederick, MD, United States

#### Edited by:

Daniela Schuster, Paracelsus Medizinische Privatuniversität, Salzburg, Austria

#### Reviewed by:

Victor Kuz'Min, National Academy of Sciences of Ukraine (NAN Ukraine), Ukraine Alexandre Varnek, Université de Strasbourg, France

\*Correspondence: Vladimir V. Poroikov<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Medicinal and Pharmaceutical Chemistry, a section of the journal Frontiers in Chemistry

Received: 29 January 2018 Accepted: 09 April 2018 Published: 26 April 2018

#### Citation:

Pogodin PV, Lagunin AA, Rudik AV, Filimonov DA, Druzhilovskiy DS, Nicklaus MC and Poroikov VV (2018) How to Achieve Better Results Using PASS-Based Virtual Screening: Case Study for Kinase Inhibitors. Front. Chem. 6:133. doi: 10.3389/fchem.2018.00133 Discovery of new pharmaceutical substances is currently boosted by the possibility of utilization of the Synthetically Accessible Virtual Inventory (SAVI) library, which includes about 283 million molecules, each annotated with a proposed synthetic one-step route from commercially available starting materials. The SAVI database is well-suited for ligand-based methods of virtual screening to select molecules for experimental testing. In this study, we compare the performance of three approaches for the analysis of structure-activity relationships that differ in their criteria for selecting of "active" and "inactive" compounds included in the training sets. PASS (Prediction of Activity Spectra for Substances), which is based on a modified Naïve Bayes algorithm, was applied since it had been shown to be robust and to provide good predictions of many biological activities based on just the structural formula of a compound even if the information in the training set is incomplete. We used different subsets of kinase inhibitors for this case study because many data are currently available on this important class of drug-like molecules. Based on the subsets of kinase inhibitors extracted from the ChEMBL 20 database we performed the PASS training, and then applied the model to ChEMBL 23 compounds not yet present in ChEMBL 20 to identify novel kinase inhibitors. As one may expect, the best prediction accuracy was obtained if only the experimentally confirmed active and inactive compounds for distinct kinases in the training procedure were used. However, for some kinases, reasonable results were obtained even if we used merged training sets, in which we designated as inactives the compounds not tested against the particular kinase. Thus, depending on the availability of data for a particular biological activity, one may choose the first or the second approach for creating ligand-based computational tools to achieve the best possible results in virtual screening.

Keywords: ChEMBL, bioactivity data, kinase inhibitors, SAR, PASS, virtual screening, classification, SAVI

# INTRODUCTION

Discovery of novel pharmaceutical agents with improved safety and efficacy is the perpetual task of medicinal chemistry (Pammolli et al., 2011). In addition to the traditional methods of chemical synthesis and pharmacological studies of various drug-like substances, in recent years substantial attention has been paid to the analysis of the general chemical-biological space (Lipinski and Hopkins, 2004; Baell and Holloway, 2010; Bon and Waldmann, 2010; López-Vallejo et al., 2012; Deng et al., 2013; Medina-Franco et al., 2013; Buonfiglio et al., 2015; Rodriguez-Esteban, 2016; Horvath et al., 2017). Such approaches significantly increase the diversity of the studied chemical libraries as well as the chances to identify the pharmaceutical agents interacting with multiple molecular targets and causing additive or synergistic desired pharmacological action (Sidorov et al., 2015; Lauria et al., 2016).

Nowadays, available chemical libraries can be divided into four categories: (1) databases containing information about structure and properties of publicly disclosed chemical compounds, e.g., PubChem (Li et al., 2010; Wang Y. et al., 2014) and ChEMBL (Bento et al., 2014); (2) databases containing information about structure of commercially available chemical samples, e.g., ZINC (Sterling and Irwin, 2015); (3) databases of virtually generated structures comprehensively covering the particular chemical space, e.g., GDB-17 (Ruddigkeit et al., 2012); (4) databases of virtually generated, synthetically accessible, structures with data on starting materials and proposed synthetic routes, e.g., SAVI (Synthetically Accessible Virtual Inventory) (Pevzner et al., 2017). Although GDB-17 is one of the largest1 currently known libraries of chemical structures containing 166.4 billion possible molecules up to 17 atoms of C, N, O, S, and halogen, SAVI looks more attractive for utilization in drug discovery because of the synthesability of its molecules. Furthermore, it was shown (Pevzner et al., 2017) that the overlap between the 93 million structures from PubChem with the 238 million SAVI database is only about 0.03%. Thus, SAVI represents a significant previously unexploited reservoir of novel structures, presumably useful for drug discovery.

To reveal the hidden pharmacological potential of the synthesizable molecules from SAVI, computer-aided virtual screening could be applied (Jorgensen, 2004; Nettles et al., 2006; Bajorath, 2014; Fujita and Winkler, 2016; Lee et al., 2016). Although structure-based methods are widely used now, ligand-based methods have important advantages (Leelananda and Lindert, 2016). In several case studies, machine learning approaches were shown to surpass the performance of both chemical similarity assessment and reverse docking (Anusevicius et al., 2015; Druzhilovskiy et al., 2016; Murtazalieva et al., 2017).

Thus, it is reasonable to analyze the probable biological activity of SAVI molecules using our computer program PASS that recently received high marks: "One of the earliest and most widely used examples of data-mining target elucidation is the continuously curated and expanded Prediction of Activity Spectra for Substances (PASS) software, which was assimilated from the bioactivites of more than 270,000 compound-ligand pairs" (Mervin et al., 2015). The PASS development started more than 25 years ago (Poroikov et al., 1993; Filimonov et al., 1995), and during this time its performance has continuously and significantly improved. PASS in its 2017 version predicts over 7,000 kinds of biological activity with an average accuracy of 94% based on the analysis of structure-activity relationships for more than 1 million known biologically active compounds.

Initially, in the PASS training set a molecule is designated as "active" if reliable information about some biological activity is found in a authoritative source (publication in a peerreviewed journal, record in curated database, etc.); otherwise, it is designated as "(conditionally) inactive." This would seem to be a reasonable approach as it has been found that if the same set of chemical compounds is studied against the same molecular target in the three different assays, only 35% of active compounds completely coincided (Lipinski and Hopkins, 2004).

Since no one chemical compound has been tested for all known biological activities, this may appear to be the incorrect designation in some cases. However, it has been shown that PASS provides reasonable estimates of structure-activity relationships despite the incompleteness of information in the training set on both chemical structures and biological activities, due to the robustness of the Naïve Bayes approach in general (Rish, 2001; Rennie et al., 2003) and the MNA descriptors and the biological activity representation used in PASS in particular (Poroikov et al., 2000).

Quantitative data on structure and activity of many chemical compounds freely available from ChEMBL and PubChem databases allow one to consider alternative approaches for creating training sets that may improve the performance of machine learning methods. Such possibilities were recently considered in several studies (Heikamp and Bajorath, 2013; Smusz et al., 2013; Kurczab et al., 2014; Afzal et al., 2015; Mervin et al., 2015).

In this work we evaluated the PASS performance in virtual screening for kinase inhibitors with training performed using three approaches, which differ with respect to what compounds were selected as inactives: (1) only experimentally validated ("true") inactives; (2) combining true and conditionally inactives; (3) only conditionally inactives. The first and second approaches have the drawback that they require enough data on true inactives.

These training strategies are both related to the multi-label classification (Tsoumakas et al., 2010; Cherman et al., 2011; Afzal et al., 2015) and positive unlabeled learning (Kilic and Tan, 2012), because one and the same classifying object may simultaneously belong to several categories [have multiple labels, i.e., inhibit more than one kinase (Martin et al., 2011) in our case study] and the problem of inactives' selection may be solved using more than one method. In contrary to various approaches of inactives' selection described by the authors (Kilic and Tan, 2012), we used only straightforward approaches, since in chemoinformatics we

1The Danish biopharmaceutical company Nuevolution announced that it had created a library of 40 trillion unique molecules (C&EN, 2017, 95: 28–33); however, the web site (https://nuevolution.com/technology) states that the company enables DNA encoded synthesis of billions of chemically diverse drug-like small molecule compounds.

are forced to deal with extremely sparse data about ligand-protein interactions and, thus, introduction of data about target-to-target relations during the training may lead to strong overfitting.

The kinases were chosen for this study because of the strong family ties among kinases that manifest themselves through common structural features and predispose kinase inhibitors to polypharmacological action (Knight et al., 2010; Gani et al., 2015; Sidorov et al., 2015). Thus, the aforementioned differences in the training set formation may lead to visible changes in the virtual screening performance. Although this class of protein targets has a privileged place in contemporary drug discovery and there are thus many compounds that have been assayed against several or even numerous kinases (Fedorov et al., 2007; Gao et al., 2013; Christmann-Franck et al., 2016; Elkins et al., 2016), multitarget action is found only for a small and diverse subset of the whole chemical-biological space (Jasial et al., 2016).

Therefore kinases and their inhibitors represent an interesting and challenging case that provides useful insights into the influence of the multitarget action of chemical compounds on the success of virtual screening studies (Merget et al., 2017). Moreover, since the multitarget action is by definition an attribute of thoroughly studied compounds, such as FDAapproved drugs (Law et al., 2014), whereas most known compounds are not thoroughly studied, our results may be extrapolated to the target classes (Barelier et al., 2015; Munoz, 2017) less extensively studied compared to kinases, to help achieve better results in virtual screening of a huge chemical library.

## MATERIALS AND METHODS

# Brief Description of PASS

PASS (Filimonov et al., 2014) is a computer program for analysis of structure-activity relationships (SAR) that allows users to perform ligand-based virtual screening for ligands of multiple targets and/or compounds with desired biological activities (Abdou et al., 2017; James and Ramanathan, 2017; Stasevych et al., 2017; Yildirim et al., 2017). Structures of chemical compounds are represented in PASS as a set of 2D atom-centric substructural descriptors called MNA (Multilevel Neighborhoods of Atoms). It was previously shown that MNA descriptors are suitable for implementation in a wide range of qualitative (classification) SAR studies and reflect structural features important for ligand–target interactions (Filimonov et al., 1999). PASS predicts biological activity profiles for chemical compounds in standardized representation: uncharged, single-component, containing at least three carbon atoms, with molecular mass not exceeding 1,250 Da. The majority of druglike molecules fulfill these conditions and clipping of the nondrug-like compounds allows us to avoid dealing with non-specific and atypical biological activities. The mathematical approach of PASS is based on a naïve Bayes classifier and its particular realization in PASS has been previously described in detail elsewhere (Filimonov et al., 2014).

The result of PASS prediction is a list of probable biological activities arranged in descending order of Pa-Pi values, where Pa is the probability of belonging to the class of "actives," while Pi is the probability of belonging to the class of "inactives". By default, Pa-Pi > 0 is considered as the cutoff for discrimination between "active" and "inactive" molecules. The result of PASS-based virtual screening for a chemical library is the list of molecules predicted as "actives"; and these could be recommended for biological testing.

#### Training and Test Datasets Data Acquisition

Every dataset used in this study was formed based on the data contained in the ChEMBL database. We chose ChEMBL because this is one of the largest freely available sources of experimental bioactivity data, its data are well-organized and documented, they are easy to acquire (via graphical web interface or API), and easy to manipulate by setting-up a local version of the database. We used the list of protein kinases and their IDs that is available via the ChEMBL web interface by browsing targets by assigned protein classes to select the subset of targets for this case study.

The training set of chemical structures and activities of chemical compounds tested for inhibition of protein kinases was extracted from the 20th version of the ChEMBL database. The ChEMBL SQL-format file dump (dump itself and instructions are available from here: ftp://ftp.ebi.ac.uk/pub/databases/chembl/ ChEMBLdb/releases/chembl\_20/)was handled in MySQL, SQL queries and PHP scripts were used to manipulate the data and write them to external SD files. Basic validation and comparison of the virtual screening performance were executed using 5-fold cross-validation.

The external test sets contained data from the up-todate 23rd version of ChEMBL on structures and activities not present in ChEMBL 20 (ftp://ftp.ebi.ac.uk/pub/databases/ chembl/ChEMBLdb/releases/chembl\_23/). ChEMBL 23 contains 1 154 583 new data on activities, among which we searched for those related to the targets involved in our study using the following procedure:

#### Data Preparation

It is known that some noise and various contradictions are stored in, and migrate from one source of bioactivity data to another, along with correct records (Kramer and Lewis, 2012; Kalliokoski et al., 2013; Tiikkainen et al., 2013; Papadatos et al., 2015). Thus, it is necessary to filter the data before using them in order to eliminate incorrect data and records that are inconsistent with the goal of the virtual screening study (Fourches et al., 2016). To achieve this goal, we used the procedures described in our previous work (Pogodin et al., 2015) with slight differences, designed to reflect the peculiarities of the targets selected for this study.

#### **Training data preparation**

First, chemical structures were filtered to eliminate incorrect molecular representations and to provide PASS with unambiguous (in the given feature space of MNA descriptors) examples for training and validation. We used an in-house command-line utility (SDF-check) to check structures for PASS compatibility and remove unsuitable ones. In addition to this, we identified structures having different ChEMBL IDs, but the same sets of MNA-descriptors, i.e., equivalent structures. We treated such structures as a single one. Thus, data on their activities were joined together, and all structures except first one encountered were deleted from the set.

After the filtering and preparation of the structures, data on bioactivities were processed to remove unreliable and inconsistent data points. In this study we used the following endpoints: Ki , Kd, IC50, Potency—assessed as concentration of compound that induces the given response; Activity, Inhibition and Residual Activity—assessed as response of the kinase, induced by the given concentration of the compound. In addition to duplicates and incomplete records, the following data were excluded:

Measurements assessed as response of the kinase (Activity, Inhibition, Residual Activity), induced by the given concentration of a chemical compound were transformed to Inhibition for convenience. The problem with the "Activity" records is their ambiguity. Such records may mean both Inhibition and Residual Activity. We clarified the meaning based on the content of the assay description field. Residual Activity and Inhibition are unambiguously connected (Residual Activity = 100 − Inhibition) and it was easier for us to deal with only one (Inhibition) type of measurement.

Records on the bioactivities were filtered semi-automatically, utilizing the content of the "Description" field from the "Assays" table. Distinct "Description" fields were reviewed and, in the cases of detection of ambiguous data, analogous records were found using suitable set of words or regular expressions. Thus, identified suspicious entries were inspected using the original publications and deleted, if the suspicions were confirmed.

To improve the validation reliability, we included in the study only those kinases that had at least 100 actives and 100 inactives (determined at the concentration 1µM). These limitations also help with the creation of accurate classifiers, which may be used for their primary purpose: to search for novel kinase inhibitors. Attempts to balance sets in terms of actives to inactives ratio were not conducted, not in the least because the assessment of the difference in the quality of classifiers built on the training data with a different ratio of actives to inactives was of interest, since two of the studied approaches for the training set creation may be considered as a method to fight skewed training data distribution (Rennie et al., 2003).

After the filtering of the bioactivities, different measurements of the inhibitory activity were used to create overall qualitative assessments for each compound designating it as active or inactive against the particular kinase. As it was mentioned earlier, we had different types of data on activities in our set for some compounds. Within these types (percentage of kinase inhibition and compound concentration producing response), median values were calculated in case a given kinase-ligand pair had multiple assessments. If concentrations of compound were available and it was less than or equal to 1µM, we designated it as active against the particular kinase. In cases where data on concentration of compound were absent we designated it as active if inhibition of the particular kinase produced by this compound was greater than or equal to 50%. Otherwise the compound was designated as inactive.

Initially we extracted from ChEMBL 458 863 records on kinase inhibition. After the completion of the all procedures described above we were left with 173 275 data points on kinase inhibitors evaluated relative to the cut-off value of 1µM (62 309 on true actives and 110 966 on true inactives at given cut-off). These data characterize interactions of 55 162 compounds with one or more of 152 human protein kinases selected for this study. These kinases represent all major families of human kinases. Our data cover a significant portion of the human kinome and allow one to search for inhibitors for all kinase families (**Figure 1**).

#### **External test data preparation**

Preparation of the data for external test set was performed in the same way as for the training set data, except for the following differences:

In total, we were able to identify 81 563 new activities against the kinases involved in this study in the 23rd version of ChEMBL. After filtering, 35 317 activities describing the action of 23 004 compounds against kinases remained.

#### **Training set formation approaches**

Filtered training set data on kinase inhibitors were stored in the local MySQL database and used to create three different training sets described below and presented in **Figure 2**. In addition, each training set was divided into the five non-overlapping and equivalent subsets for subsequent stratified 5-fold crossvalidation (5-f CV).

#### **Individual sets (I-sets)**

The tested compounds for each kinase were sorted from the most active to the most inactive and, in this order, they were written to the five SD files: the first compound in the rank was placed into the first subset, the second compound into the second subset, the fifth compound into the fifth subset, the sixth then again into the first subset and so on; until each compound was placed into the each corresponding subset. The subsets were created in this way to be equivalent in terms of the total number of compounds and similar to each other in the degree of inhibitory activity of the placed compounds.

#### **Merged actives and inactives set (MAI-set)**

Then, we merged the first, second etc. subsets for each of the 152 kinases. If identical compounds were found in different subsets, only the structural formula was retained with all its kinase inhibiting activity data. As a result, we obtained 5 combined MAI-subsets, which were equivalent to the I-subsets because these subsets contained the same active compounds.

#### **Merged actives set (MA-set)**

This set was created in the same manner as MAI-set, but the true inactives were excluded.

#### Quality Metrics

We used the following metrics to evaluate the results of our ligand-based virtual screening of kinase inhibitors:

$$\text{SENSITIVITY(RECALL)} = \text{TP/(TP} + \text{FN)}\tag{1}$$

$$\text{SPECIFICITY} = \text{TN/(TN} + \text{FP)}\tag{2}$$

$$\text{Recent use presented in section }\mathsf{M}\text{-graded \"act\" sites and \"rates\" (\mathsf{A}\mathsf{A}\mathsf{A}\mathsf{L}\mathsf{-set})\text{ and \"Mover\" (\mathsf{A}\mathsf{C}\mathsf{A}\mathsf{L}\mathsf{-set})\text{.}}$$

$$\text{BALANCEED ACCCY} = \underbrace{\frac{1}{2}}\_{\text{---}} \text{\*} \underbrace{\text{TP}}\_{\text{---}} + \underbrace{\text{TN}}\_{\text{---}} \text{)} \qquad \text{(3)}$$

$$\text{PRECISION} = \text{TP/(TP} + \text{FP)}\tag{4}$$

$$\text{F2ECCISION} \* \text{RECALL}\tag{5}$$

$$F1 = 2 \* \frac{\textbf{\underset{PECCION}} + \textbf{\underset{RECALL}}}{\textbf{PRECSION} + \textbf{RECALL}} \quad \text{(5)}$$

$$\text{OCAUC} = P \left( \text{Rank}\_{\text{active}} < \text{Rank}\_{\text{inactive}} \right)$$

$$ROCAUC = P\left(Rank\_{active\_i} < Rank\_{inactive\_i}\right)$$

$$\text{in Uniform distribution} \tag{6}$$

$$BEDROC = P\left(Rank\_{active\_i} < Rank\_{inactive\_i}\right)$$

$$\text{in exponential Probability Density}$$

Function (PDF) with parameter α, IF α ∗Ra << 1

(7)

Metrics (1–6) are appropriate for the evaluation of the performance of the classification procedure, which determines the upper limits of the virtual screening quality under condition where every compound predicted as active is screened experimentally.

Boltzmann-Enhanced Discrimination of Receiver Operating Characteristic (BEDROC) (Truchon and Bayly, 2007) (Equation 7) represents the adaptation of ROC AUC metric to conditions under which detection of maximal number of TPs in a certain top fraction of the set is more important than general recognition. Thus, it is designed to evaluate the early detection rate, i.e., to assess the quality of virtual screening under the limitation that it is possible to evaluate experimentally only small fraction of top rated compounds from the whole library. Parameter α in the BEDROC AUC is inversely related to the size of the top fraction that will contribute to 80% of the score value while the other 20% will come from the assessment of the remaining part of the set. Values of α that were used in this study, and the corresponding top fractions of the sets, are given in **Table 1**.

TABLE 1 | Values of BEDROC parameter α and corresponding top fractions of sets.

# Performance Assessments

#### Stratified 5-Fold Cross-Validation

The training data had been divided into the five subsets in such a way that the average numbers of actives and inactives were approximately equal in all subsets (Refaeilzadeh et al., 2009). Four subsets from each set were used for the training, while one subset was used as the external test set. This procedure was repeated five times; each time a different subset was used as the external test set. The main differences from the standard 5-fold CV were:

The overall scheme for performance evaluation is given in **Figure 3**.

Such validation procedure provides reliable quality assessments for classifiers, since every compound in the test sets had experimental test results against a particular kinase. Besides, such an approach provides the conditions for comparison that are close to those observed in real research projects when one tries to find novel activity for a compound already included in the training set with some other activities. Such situations occur in drug repurposing projects or in in silico toxicological studies (Wang Y. J. et al., 2014).

The results of the predictions were assessed using the metrics described in the Materials and Methods section. Unfortunately, at least one of them, BEDROC, may suffer from saturation. To avoid this, the ration of actives to inactives for a set (Ra in Formula 7) must be low enough to fulfill the condition given in Formula 7.

The condition of low fraction of actives in the set seems acceptable and reasonable in the context of high throughput screening, which typically provides a number of hits below 5% (Murray and Wigglesworth, 2017). However, the data on kinase inhibitors from our set do not fulfill this condition. Thus, the saturation effect on BEDROC was expected to affect the results of our study. To avoid BEDROC saturation, we implemented the procedure of random sampling with replacement as realized in R package mlr (Bischl et al., 2016) applied to the prediction results. We undersampled the portions of actives and oversampled the portions of inactives for each kinase. Factors to under- and oversample actives and inactives were chosen in such a way that numbers of actives and inactives in the resampled set became equal to approximately 60 and 60 000, respectively (Formulae 8, 9). Thus, we maintained the same actives rate in the resampled sets, which was chosen to be approximately 0.001. This rate is low enough to calculate BEDROC values for each α level selected for this study without the risk of saturation.

The resampling procedure was repeated 5 000 times for each type of sets and each kinase to achieve statistical significance in the subsequent assessment of differences between the results. BEDROC values were calculated on the resampled data using the R package enrichVS (http://cran. r-project.org/web/packages/enrichvs/index.html) for each resampled set. ROC AUC was also calculated using the R package pROC (Robin et al., 2011). To increase the speed of obtaining resampling results, we performed calculations in parallel mode using R package "parallel" (https://stat.ethz.ch/ R-manual/R-devel/library/parallel/doc/parallel.pdf). Values of the classification quality metrics achieved in cross-validation and training set composition could be found in Supplementary Table 1.

#### Virtual Screening of the External Test Set

Prepared data from 23rd version of ChEMBL was used for forming the test sets according to the procedure used for preparation of the training I-sets. During the external validation (Chen et al., 2012) with these sets we calculated BEDROC values for the resampled prediction results. Values of the classification quality metrics achieved in external validation and training set composition could be found in Supplementary Table 2.

### Comparison of the Results Obtained Using Different Training Approaches

The Tukey honest significant difference (HSD) test was used along with the analysis of variance to compare the quality of the created PASS classifiers based on the different types of training sets. These quality parameters include BEDROC for the resampled results; sensitivity, specificity, balanced accuracy, precision, F1 score and ROC AUC for the original results. The analysis was performed at a P-value < 0.05 using the functions "aov" and "TukeyHSD" from the R standard library. This provides the ranked lists for three PASS classifiers, which allows one to evaluate their performance.

# RESULTS

### Stratified 5-Fold Cross-Validation

All classification metrics values averaged over all kinases except the sensitivity values were slightly higher for the results achieved by classifiers trained on I-sets. Statistical analysis indicates that results obtained using the I-sets differ significantly from those obtained with the MA and MAI sets (**Figure 4**). The results of classifiers trained on the MA- and MAI-sets do not differ at the given level of significance from each other.

We used the resampled results to calculate values of BEDROC at different degrees of early recognition of TP (via varying values of α). These values were grouped according to the types of sets used for the training, and then averaged over the kinases in a manner similar to the way the original results were obtained. Statistical analysis of these data shows that classifiers trained on I-sets significantly outperform classifiers trained on MAI-sets and those, in turn, outperform classifiers trained on MA-sets (**Figure 5**) for any α value used in the study.

Also, using the resampled results, we were able not only to compare different approaches for the training by averaging values of the selected metrics across kinases, but to select the most adequate approach for each kinase individually. This was because

after the resampling procedure repeated 5,000 times, we had enough data points to estimate the statistical significance. Such estimation was performed as follows: at the level of the P-value chosen earlier, less than 0.05, we found that for most of the kinases the best approach for training is to use I-sets; nonetheless, for some kinases it is better to use MA- or MAI-sets (**Figure 6**) according to our evaluation. In total, we depicted 13 kinases for which the classifiers trained using MA- or MAI-sets performed better in early recognition of TP at at least three levels of α.

# Virtual Screening of External Test Set

Since we did not impose any limitations on the number of actives and inactives in our external test set, we were not able to calculate values for all the metrics for each kinase. We excluded such

kinases before averaging the values of the classification metrics across the different training approaches, thus only results for 128 kinases were compared.

The main conclusions of the comparison of Specificity, Balanced Accuracy, and AUC values are similar to those obtained using 5-f CV: The training approach I provided significantly better results than those introducing conditionally inactives (MA and MAI). No significant difference for the other metrics was found (**Figure 7**).

To compare the earliness of actives detection achieved using different training approaches, we resampled results of the inhibitory activity prediction for each kinase and calculate BEDROC values. In this part of the study only results related to kinases having at least 20 actives and 20 inactives in the external test set were included. This restriction was imposed to exclude the influence of extreme cases, where only few actives and inactives exist. Despite the introduced restrictions, we were forced to change the resampling protocol in some cases; if the kinase had less than 60 actives, we used an oversampling procedure instead of undersampling to make sure we had 60 actives.

The main result of the comparison of BEDROC values was concordant to those obtained using 5-f CV: at each value of the criterion α, training using I-sets led to the better results than training performed using MA- or MAI-set, while MAI-sets outperformed MA-sets (**Figure 8**).

# Correlations Between the Values of Metrics and Actives to Inactives Ratio in the Sets

shape of the points corresponds to the type of the training set.

We also analyzed the behavior of the employed accuracy metrics for different actives/inactives ratios, to be sure that they give an unbiased picture.

Values of Precision and F1-score were found to show correlations with the actives to inactives ratio in the test sets. Thus, we conclude that sets' imbalance affects Precision and F1 score values, while the other metrics are significantly more robust (see Supplementary Figure 1), especially AUC and Balanced Accuracy.

# Applicability Domain Estimation

of α. Results were obtained using the external test set.

To estimate the applicability domain, we calculated the values of the classification quality metrics for those cases where compounds had a certain number of new MNA-descriptors not found in the training set. In this case we merged the results over all kinases to obtain sufficient numbers of data points.

We showed that in the case of the results achieved using I-sets for training, the performance of the classifiers decreases linearly with increasing number of new MNA descriptors. In contrast to this, for the results achieved using MA- and MAI-sets for training, we were unable to find a strong dependence between the number of new MNA descriptors and the performance of the classifiers. Still, these results should be treated with caution, since the percentage of data points involved in this assessment decreases drastically with increasing number of new MNA descriptors, especially for the classifiers built using MAI- and MA-training sets (see **Figure 9**).

In the case of the classifiers built using I-sets for training we can judge that the applicability domain includes those compounds which have 4 or fewer new MNA descriptors, since the average balanced accuracy and AUC exceeded 0.7.

### DISCUSSION

In contrast to the many contemporary studies in the field of the virtual screening, in this work no decoys (Irwin, 2008) were used to assess the enrichment achieved in virtual screening of large datasets. Instead, validation and subsequent comparison of the different training approaches were performed using only experimentally tested compounds, both actives and inactives. Today, due to the constant growth of available computational resources and amount of bioactivity data, it is possible to do this using 5-f CV and true external test sets. Moreover, since negative influence of the conditionally inactive compounds involved in training was shown, this makes us wonder: if conditionally inactives can do harm during training, are decoys good for testing? The exact answer is not known yet, but the risk of reaching wrong conclusions may be mitigated by using resampling-based approaches in parallel with, or instead of, decoys.

Our study represents a quantitative assessment of the tradeoff between the initial requirements on the training data and the quality of PASS-based virtual screening. We have shown that the most efficient training approach for the ligand-based virtual screening system is to use the true actives and inactives for each target. This approach outperformed those where conditionally inactive compounds were introduced, in both classification quality and earliness of the detection. Moreover, in this case we observe a strong dependence of the performance depending on the number of new descriptors in the structures of the test compounds.

According to the analysis of the data from our training set, the higher the number of kinases for which compounds are tested, the more activities are found. Thus, using MA and MAI sets for training, some unknown actives could be treated as conditionally inactives (**Figure 10**). This may shed some light onto the problem of promiscuity of kinase inhibitors, which are often discussed as polypharmacological drugs. However, analysis of the content of bioactivity databases such as ChEMBL has shown that the average degree of promiscuity of such compounds is not so high (Hu et al., 2014). According to our results there is no contradiction between these points of view: kinase inhibitors tend to show promiscuity, but at the moment most of them have been studied against only a rather limited number of kinases.

Nevertheless, using MA and MAI approaches, it is possible to achieve good virtual screening results too, despite the softer requirements on the amount and quality of the training data. These approaches may be implemented in cases when only few active compounds are known, even in the absence of inactives, which helps expand the druggable target space and find new modes of action for existing molecular targets.

From this perspective it is surprising that we also found 13 kinases for which virtual screening may be performed more efficiently using training approaches introducing conditionally inactive compounds. This means that using machine learning it is easier to distinguish between inhibitors of these kinases and compounds tested against other kinases, than between their inhibitors and inactives at the given concentration cutoff. This fact can possibly be explained by the systematical shift in compounds selection for testing against these kinases. Also, it may indicate the importance of small structural changes in related targets leading to larger changes in inhibitor potency, since these 13 kinases are diverse, they belong to different families represented in our set and, in the case of other members of their families, introduction of the conditionally inactive compounds leads to the observed negative consequences. Thus, we show that virtual screening performance may benefit from the introduction of conditionally inactive compounds if these

compounds are unfamiliar to the main target. Unfortunately, this knowledge is risky to apply to achieve better results in ligandbased virtual screening, since our knowledge on target-target relations mediated by common ligands are generally based on sparse training sets.

We obtained rather good results of both external (quasi prospective) and cross-validation. However, in case of data on kinase inhibitors extracted from ChEMBL, one initially deals with the pre-selected compounds studied in the appropriate biological activity area, which provides good predictivity, particularly using the approach based on individual sub-sets.

Big libraries like SAVI contain diverse and previously not investigated chemical structures, including compounds other than those possessing known ligand-related target signatures (Sidorov et al., 2015). To achieve the best predictivity for such library, it seems reasonable to make pre-selection with the standard PASS approach using conditionally inactive compounds. As we already mentioned above, PASS provides satisfactory results of prediction despite the incompleteness of data in the training set (Poroikov et al., 2000). Moreover, in this work, we showed that classifiers created using the merged training sets did not exhibit the significant dependence between the prediction quality and the number of new MNA descriptors contained in the predicted chemical structures.

Consequently, we propose two-steps procedure to analyze the big and diverse chemical libraries. At the first step, pre-selection is performed using the general classifier that took into account the conditionally inactives. At the second step, one may more thoroughly discriminate between the active hits and putatively inactive structures using the specific classifier that is based only on the real actives and inactives.

# CONCLUSIONS

In this study, we compared the performance of three approaches for the analysis of structure-activity relationships that differ in their criteria for selecting "active" and "inactive" compounds for the training sets. We used the program PASS to build classifiers based on different subsets of kinase inhibitors extracted from ChEMBL 20 (for training and 5f-CV) and ChEMBL 23 (for external, quasi-prospective validation). The highest classification and early recognition quality was obtained by using individual

training sets for each kinase containing only experimental data. Nevertheless, other training strategies can provide acceptable results even in the absence of data on known inactives, which is often the case with the novel targets (Russ and Lampel, 2005; Nguyen et al., 2017). We assessed the applicability domain of our classifiers: while classifiers trained using individual sets expose strong dependence of the prediction quality on the predicted compounds' novelty, training strategies employing merged sets are much less sensitive to the novelty of predicted compounds.

Taken together these findings allow us to suggest that one can benefit most from using combinations of different training strategies when exploring huge chemical libraries containing diverse structures of unexplored chemical compounds.

### REFERENCES

Abdou, W. M., Shaddy, A. A., and Kamel, A. A. (2017). Structure-based design and synthesis of acyclic and substituted heterocyclic phosphonates linearly linked to

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

This work is supported by the Russian Foundation for Basic Research grant No. 17-54-30015-NIH\_a.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2018.00133/full#supplementary-material

Afzal, A. M., Mussa, H. Y., Turner, R. E., Bender, A., and Glen, R. C. (2015). A multi-label approach to target prediction taking ligand

thiazolobenzimidazoles as potent hydrophilic antineoplastic agents. Chem. Pap. 71, 1961–1973. doi: 10.1007/s11696-017-0190-z

promiscuity into account. J. Cheminform. 7, 1–14. doi: 10.1186/s13321-015-0 071-9

and comparison with the other descriptors. J. Chem. Inf. Comput. Sci. 39, 666–670. doi: 10.1021/ci980335o

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Pogodin, Lagunin, Rudik, Filimonov, Druzhilovskiy, Nicklaus and Poroikov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Discovery of Natural Products as Novel and Potent FXR Antagonists by Virtual Screening

Yanyan Diao1 , Jing Jiang1 , Shoude Zhang1 , Shiliang Li 1 , Lei Shan2 , Jin Huang1 , Weidong Zhang2 and Honglin Li 1 \*

1 Shanghai Key Laboratory of New Drug Design, State Key Laboratory of Bioreactor Engineering, School of Pharmacy, East China University of Science and Technology, Shanghai, China, 2 Department of Phytochemistry, School of Pharmacy, Second Military Medical University, Shanghai, China

#### Edited by:

Simone Brogi, University of Siena, Italy

#### Reviewed by:

Marco Tutone, Università degli Studi di Palermo, Italy Denis Fourches, North Carolina State University, United States Michal Brylinski, Louisiana State University, United States Francesco Ortuso, Università degli Studi Magna Græcia di Catanzaro, Italy

> \*Correspondence: Honglin Li<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Medicinal and Pharmaceutical Chemistry, a section of the journal Frontiers in Chemistry

Received: 08 January 2018 Accepted: 12 April 2018 Published: 30 April 2018

#### Citation:

Diao Y, Jiang J, Zhang S, Li S, Shan L, Huang J, Zhang W and Li H (2018) Discovery of Natural Products as Novel and Potent FXR Antagonists by Virtual Screening. Front. Chem. 6:140. doi: 10.3389/fchem.2018.00140 Farnesoid X receptor (FXR) is a member of nuclear receptor family involved in multiple physiological processes through regulating specific target genes. The critical role of FXR as a transcriptional regulator makes it a promising target for diverse diseases, especially those related to metabolic disorders such as diabetes and cholestasis. However, the underlying activation mechanism of FXR is still a blur owing to the absence of proper FXR modulators. To identify potential FXR modulators, an in-house natural product database (NPD) containing over 4,000 compounds was screened by structure-based virtual screening strategy and subsequent hit-based similarity searching method. After the yeast two-hybrid (Y2H) assay, six natural products were identified as FXR antagonists which blocked the CDCA-induced SRC-1 association. The IC50 values of compounds 2a, a diterpene bearing polycyclic skeleton, and 3a, named daphneone with chain scaffold, are as low as 1.29 and 1.79µM, respectively. Compared to the control compound guggulsterone (IC50 = 6.47µM), compounds 2a and 3a displayed 5- and 3-fold higher antagonistic activities against FXR, respectively. Remarkably, the two representative compounds shared low topological similarities with other reported FXR antagonists. According to the putative binding poses, the molecular basis of these antagonists against FXR was also elucidated in this report.

Keywords: FXR, antagonist, virtual screening, molecular docking, similarity searching, natural product

### INTRODUCTION

The farnesoid X receptor (FXR, NR1H4), a member of the metabolic nuclear receptor superfamily, regulates the expressions and activities of a broad spectrum of genes. Since 1995 when FXR was isolated from a rat cDNA library for the first time (Forman et al., 1995), the studies on the physiological functions of FXR have been appealing and challenging. FXR is conserved from teleost fish to human beings (Maglich et al., 2003) and is abundantly expressed in liver, intestine, and kidney. As the endogenous receptor of bile acids, FXR can be activated by chenodeoxycholic acid (CDCA), lithocholic acid (LCA), deoxycholic acid (DCA), and many other bile acids (Makishima et al., 1999). In addition, FXR is also reported to exert regulatory roles in lipoprotein and glucose homeostasis, fatty acid and triglyceride synthesis, liver regeneration, and bacterial growth in the intestine (Lee et al., 2006; Wang et al., 2008). All these accumulating data make FXR a promising pharmaceutical target for multiple diseases, especially those related to metabolic disorders such as diabetes and cholestasis (Schaap et al., 2014; Gonzalez et al., 2016; Yuan and Li, 2016; Filho et al., 2017).

As a typical nuclear receptor, FXR shares common structural characteristics with other members of this superfamily, which comprises a highly conserved DNA-binding domain (DBD), a moderately conserved ligand-binding domain (LBD), and a ligand-dependent transcriptional activation domain (AF-2) (Pellicciari et al., 2005).Upon the binding of proper ligand to the LBD, FXR will undergo a conformational change, which is critical to determine whether a coactivator or a corepressor binds efficiently to the AF-2 motif. If activated by appropriate agonists, the recruitment of coactivators (such as SRC-1, DRIP, and PRMT) to FXR occurs, which further up- or down-regulates the expressions of certain target genes. While for antagonists, the association of FXR with activators will be hindered (Lew et al., 2004). Although it is widely accepted that FXR participates in many biological processes, owing to the diversity and complexity of target genes involved in the FXR signaling pathways (Zhang and Edwards, 2008), the physiological functions of FXR haven't been clearly defined. Therefore, it is still an essential step to identify potential FXR modulators, which may contribute to the elucidation of physiological effects of FXR and provide novel opportunities for the treatment of metabolic diseases by targeting FXR.

Apart from the natural bile acid ligands with steroidal skeleton, over 700 structurally diverse FXR modulators have also been identified (Gaulton et al., 2017), most of which function as agonists (Carotti et al., 2014). Obeticholic acid (6αethyl-chenodeoxycholic acid, 6-ECDCA), a semi-synthetic bile acid analog with highly potent FXR agonistic activity (EC50 = 0.099µM) (Pellicciari et al., 2002), is the first FDA-approved drug that is used for treating primary biliary cholangitis (PBC) (Nevens et al., 2016). In contrast, the development of FXR antagonists, which are also useful chemical tools to unravel the physiological roles and relative clinical significance of FXR (Li et al., 2013), does not seem to be satisfactory due to the scanty number of potent FXR antagonists that have been reported so far. Guggulsterone, a natural product extracted from the resin of the guggul tree, is the most described FXR antagonist, with the ability of blocking the agonist-induced coactivator recruitment and decreasing the hepatic cholesterol in wild-type mice (Urizar et al., 2002). However, the researches on guggulsterone are still confined to preclinical and academic studies because of the complexity of its mechanism of action (Fiorucci et al., 2010; Yamada and Sugimoto, 2016). Although other natural or synthetic FXR antagonists have also been developed (**Figure 1**; Wu et al., 2002; Dussault et al., 2003; Nam et al., 2007; Choi et al., 2011; Huang et al., 2012; Xu et al., 2015), further pharmaceutically relevant activities were rarely reported. Herein, six natural products were identified as antagonists from an inhouse natural product database (NPD) through virtual screening strategy and subsequent biological experiment validation. In the yeast two-hybrid (Y2H) assay (Fields and Sternglanz, 1994; Lin and Lai, 2017), these compounds could abolish CDCAinduced FXR activation at micromolar level. We hope the natural products revealed in this study will offer novel scaffolds for uncovering new FXR regulatory mechanism and provide insights into potential development for further discovery of FXR modulators.

# MATERIALS AND METHODS

#### Structure-Based Virtual Screening Protein Preparation

The crystal structures of FXR-LBD in complex with 6-ECDCA (Mi et al., 2003) (PDB code 1OSV, a dimer and chain B was used) and fexaramine (Downes et al., 2003) (PDB code 1OSH, monomer) were obtained from the Protein Data Bank. The synthetically modified bile acid ligand 6-ECDCA was derived from the structure of CDCA, but showed almost 100-fold more potent FXR agonistic activity than CDCA and did not activate other nuclear receptors. Fexaramine is a synthetic nonsteroidal FXR agonist, which was identified by optimization of a benzopyrane-based combinatorial derived library. The coactivators and all the water molecules were removed. Hydrogen atoms and charges were added during a brief relaxation performed using the "Protein Preparation Wizard" workflow in Maestro 10.1. After the hydrogen bond network was optimized, the crystal structure was minimized until the root-mean-square deviation (RMSD) between the minimized structure and the starting structure reached 0.3 Å with OPLS\_2005 force field.

#### Glide Docking

The grid-enclosing box was placed on the centroid of the crystallographic ligand in the optimized protein structure and defined to enclose residues located within 15.0 Å of the binding pocket. A scaling factor of 0.8 was set to van der Waals (VDW) radii of those receptor atoms with partial atomic charges of less than 0.15 to soften the nonpolar parts of the receptor. After addition of hydrogen atoms and ionization at a pH range of 5.0- 9.0, the three-dimensional structures of compounds in the NPD were generated with Ligprep v3.3 module. Standard precision (SP) and extra precision (XP) approaches of Glide (Friesner et al., 2004; Halgren et al., 2004) were respectively adopted to dock the molecules into the binding site with the default parameters, and only the top one pose for each molecule were retained. After parallel Glide SP scorings using two different protein structures (PDB codes 1OSV and 1OSH), the top 500 docking poses were reserved for each docking calculation and subjected to XP calculation with a more precise scoring function, and the top 200 docking poses were retained, respectively, for further visual observation.

# Hit-Based Similarity Searching

The similarity searching process was accomplished in Pipeline Pilot v7.5, and two of the most potent FXR antagonists **2a** and **3a** were used as query molecules, respectively. The Tanimoto coefficient (Tc) of similarity between the query molecule and the target molecule was calculated using SciTegic functional connectivity fingerprints of radius 4 (FCFP\_4) (Bender et al., 2009). The minimum Tc was set to a low value of 0.3, to maximize the number of obtained analogs.

# SRC-1 Recruitment Assay

#### Materials

The restriction and modification enzymes in this work were obtained from New England Biolabs (Beijing, China). P-nitrophenyl α-D-galactopyranoside (PNP-α-Gal), yeast nitrogen base without amino acids, agar, PEG3350, dimethyl sulfoxide (DMSO), lithium acetate, and glucose were all purchased from Sigma (Shanghai, China). The yeast expression plasmids pGADT7 and pGBKT7 were from Clontech (Palo Alto, CA), and CDCA was from Merck. The dropout supplement free from leucine and tryptophan (-Leu/-Trp DO supplement) was bought from Takara, and Salmon Sperm DNA was obtained from invitrogen. The yeast strain AH109 was purchased from Clontech (Palo Alto, CA).

#### Plasmid Construction

Based on the genome sequences of FXRα (GenBank accession no. NC 000012.10), human FXRα-LBD (200-473 AA) was sub-cloned into vector pGBKT-7 using NdeI and BamHI restrict enzyme sites. The primers used for PCR amplification were listed as follows: FXRα-LBD (sense) 5 ′ -ATCATATG-GAAATTCAGTGTAAATCTAAGCG-3′ ,

(anti-sense) 5′ -ATGGATCCTCACTGCA-CGTCCCA-3′ . The combination plasmid pGADT7-SRC-1 was prepared as described previously (Lin et al., 2008), by amplifying with the following primers: (sense) 5′ -CAGAATTC-CATAACAATGACAGACTTTCA-3′ and (anti-sense) 5 ′ -AAGGATCCCACCTTTA- CATCATCCAGGCT-3′ .

#### Y2H System Construction

We constructed the Y2H for FXR by yeast co-transformation with pGBKT7-FXR LBD (BD) and pGADT7-SRC-1 (AD) according to the lithium acetate method. Briefly, 500 ng of BD and AD were added to 50 µL of the yeast competent cells and mixed with 36 µL of lithium acetate, 240 µL of 50% PEG3350, and 50 ng single-strain DNA at 30◦C for 30 min, followed by heat-shock (250 rpm) at 42◦C for 30 min. The mixture was subsequently spread on a drop-out-agar plate without leucine and tryptophan (6.7 g/L yeast nitrogen base without amino acids, 1.54 g/L -Leu/-Trp DO supplement, 20 g/L glucose, 20 g/L agar). The plates were incubated at 30◦C for 48 h for yeast growth and the PCR method was used to confirm the successful transformation.

#### Y2H Assay

We performed Y2H assay to determine the agonistic or antagonistic activities of the compounds. Yeast transformations were incubated with either a control vehicle (DMSO) or the indicated compounds for 24 h in an hFXR agonist testing, and in antagonist assays treated with tested compounds plus 10µM CDCA. The quantitative α-galactosidase activity assays were carried out by using PNP-α-Gal as the substrate according to the Clontech manual. Each experiment was repeated three times independently.

The α-galactosidase activity was calculated according to the following formula:

$$\begin{aligned} \alpha &\quad \text{-galactoside activity} \text{[milliumits/(mL} \times \text{cell)]}\\ &= \frac{OD\_{410} \times V\_f \times 1000}{(\varepsilon \times b) \times t \times V\_i \times OD\_{600}} \end{aligned}$$

where t is the elapsed time of incubation (min), Vf is the final volume of assay (200 µL), Vi is the volume of culture medium supernatant added (16 µL), OD600 is the optical density of overnight culture, and ε×b is the p-nitrophenol molar absorptivity at 410 nm×the light path (cm) = 10.5 mL/µmol (Yeast Protocols Handbook PT3024-1, Clontech).

The agonistic activation and inhibition rates (%) were calculated as follows:

$$\text{againistic activation} = \frac{GA\_{treated}}{GA\_{DMSO}}$$

$$\text{inhibition rate(\%)} = \frac{GA\_{CDCA} - GA\_{treated}}{GA\_{CDCA} - GA\_{DMSO}}$$

where GA indicates α-galactosidase activity.

#### Chemistry

The NPD is our in-house collection of over 4,000 natural products isolated from about 100 plants and their structures were established by extensive spectroscopic. The purities of all compounds were checked by using NMR and HPLC (purities ≥ 95%).

#### RESULTS AND DISCUSSIONS

#### Protein Structure Selection and Redocking Validation

The binding of proper modulators to LBD is the molecular basis for FXR activation that triggers the conformational change of FXR, the subsequent coactivators association, and the final target genes regulation. Previous studies have shown that FXR can be activated by structurally diverse agonists, and an array of crystal structures of FXR-LBD complexed with agonists have been solved. The agonist-binding pocket is positioned in the interior of the LBD, and agonists with different skeletons display enormously distinct binding features (Maloney et al., 2000; Soisson et al., 2008; Flatt et al., 2009; Jin et al., 2013). Computational studies, along with crystallographic experiments (Xu et al., 2015), support the notion that the agonist-binding pocket can also be occupied by antagonists (Meyer et al., 2005). Given the variability of the binding pockets occupied by distinct modulators, it is necessary to use different FXR crystallographic models in the virtual screening process, in order to maximize the diversity of hit compounds. At the time when we initiated this study, there were in total 27 crystal structures of FXR in complex with different modulators in the Protein Data Bank. After deleting the structures bound with structurally analogous ligands and those without published biological data, 10 unique protein structures were reserved for further analysis.

Protein flexibility has vital influence on ligand recognition, and even subtle protein conformational changes can significantly affect the results of docking simulations (Kitchen et al., 2004). Nevertheless, the receptor is usually held rigid for most docking procedures, including Glide used in this study, to speed up virtual screening of large databases. To compensate for the limitations of rigid protein conformation, as well as to simply the computational simulations, we decided to select two receptors representing dramatically dissimilar conformations to the other reported crystal structures to perform two independent docking calculations. Both the binding site similarity and ligand similarity profiles were taken into consideration for docking model selection. On the one hand, pairwise binding pocket similarities among the 10 structures were calculated using our in-house program PocketShape, which is designed for computational evaluation of the binding site similarity based on pocket shape and property and could be accessed through the webserver SiteMapper (http://lilab.ecust.edu.cn/). Residues within 5Å distance around the ligand were extracted as the binding site for each structure. Typically, two binding sites with a score value over than 0.8 are considered similar. On the other hand, pairwise molecular similarities among the 10 crystallographic ligands were calculated with SciTegic FCFP\_4 fingerprints in Pipeline Pilot v7.5, and a Tc cutoff value of 0.6 was set to define similarity.

The binding pocket and ligand similarity values of the 10 unique crystal structures were plotted in Figure S1. Meanwhile, for each model, the average pocket and ligand similarity values were calculated, respectively, to evaluate its uniqueness. Four crystal structures 1OSH, 1OSV, 3OLF, and 4WVD were considered, as all of them have top five minimum average values ranked by either pocket or ligand similarity calculations. Although both the pocket and ligand similarity values of 4WVD scatter in a low range, its ligand ivermectin is a large macrocyclic lactone (Jin et al., 2013), making the ligand binding pocket expand to a volume of 1081 Å3 (Dundas et al., 2006), which is likely to cause artificial enrichment of molecules with large sizes and high molecular weights in docking procedure (Verdonk et al., 2004). Therefore, we did not choose 4WVD. The Tc values of the FXR agonist 6-ECDCA in 1OSV to other ligands are relatively low with predominant scattering below 0.2, and its binding pocket also reveals uniqueness with the average similarity value of 0.62. Moreover, 6-ECDCA is the only FDAapproved FXR modulator so far. Therefore, the 1OSV model was first selected for docking simulation (Mi et al., 2003). After thoroughly analyzing the binding interactions, the biological activities against FXR, as well as the X-ray crystal parameters, the second model of 1OSH was preferentially reserved, the ligand of which has a lower Tc value to 6-ECDCA than that of 3OLF (0.13 vs. 0.15).

Alignment of the 6-ECDCA and fexaramine binding sites exhibited substantial differences in both shape and surrounded residues (Figure S2A). And the tremendous dissimilarities of the two ligands in terms of topological structures and induced binding conformations render them extend to non-overlapping space. To evaluate the docking accuracy of Glide, the two cocrystallized ligands were redocked into the active pocket using Glide XP scoring for each structure. Superimposition of the best redocked poses and the experimental structures (Figures S2B,C) gave the RMSD values of 0.53 Å for 6-ECDCA and 0.65 Å for fexaramine, indicating the robustness of Glide in accurately reproducing the bioactive conformations of the ligands for our two docking models. Cross-docking calculations were also carried out where each crystal structure ligand was docked into the binding pocket of the other. The ligand 6-ECDCA could be docked into 1OSH by Glide XP mode, but the predicted binding pose deviated greatly from the crystallized bioactive conformation (Figure S2D). In agreement with the enormous differences of the two binding sites, fexaramine could not be accommodated by the smaller 6-ECDCA-binding site of 1OSV, hence no proposed docking pose was obtained.

#### Virtual Screening

To search for potent FXR modulators, structure-based virtual screening strategy, an effective method to identify novel ligands

based on predicted binding poses and docking scores, was executed using Glide v6.6 (Maestro v10.1, Schrödinger Inc.). And we speculated that if the docking pose of a certain compound to the agonist-binding pocket was computationally favorable, it could be an effective FXR modulator, either an agonist or an antagonist. The screened natural products database is a collection of over 4000 natural products isolated from about 100 plants, the structures of which have been validated by our researchers. After hierarchical virtual screenings independently implemented with Glide (**Figure 2**) by using two crystal structures of the receptor, a total of 400 top-ranking compounds were retrieved as candidates from the NPD. These candidates were then subjected to visual inspection to remove those that are likely to be nonbinders. With consideration of key interactions observed from the crystal structures, such as predominant VDW complementarity and some critical hydrogen bonds with polar residues, each docked pose of these candidates was carefully checked to delete the inappropriate compounds. Meanwhile, a specific focus was put on the sizes of the candidates, and those compounds with relatively large groups protruding out of the binding pocket were not considered. Additionally, the compounds with the same scaffold were reserved with a maximum of three to maximize the structural diversity. A total of 30 candidates were finally selected for further bioactivity assay. In the coactivatorrecruitment assay based on the Y2H system, none of the 30 compounds could enhance the association of SRC-1 to FXR-LBD, thus no agonist was found. Intriguingly, four compounds (**1a, 2a, 3a, 3c**) strongly inhibited the CDCA-induced SRC-1 recruitment with the inhibition rate higher than 50% in the concentration of 25µM, displaying apparent FXR antagonistic profiles. The IC50 values of the four compounds were determined (**Table 1,** Figure S4), and guggulsterone was used as the reference compound.

Starting from the potent FXR antagonists **2a** and **3a** as hit compounds, the in-house NPD was re-screened using similarity searching method, to obtain more potent derivatives as well as to establish underlying structure-activity relationships (SAR). A prior minimum Tc value of 0.6 was first set to retrieve analogous compounds. Unfortunately, only several derivatives were obtained for each hit compound, probably owing to the heterogeneity of the in-house NPD. We, thus, chose a relatively low threshold of 0.3 to maximize the number of obtained analogs. Subsequently, the analogs were also manually checked, and only compounds possessing the same skeleton to the hit compound and with proper sizes were selected. After further in vitro Y2H assay, two compounds **3d** and **3e** with the Tc values of 0.4 and 0.41, respectively, to compound **3a**, were found to display moderate antagonistic activities against FXR. As illustrated in **Table 1**, compounds **3d** and **3e** share the phenol moiety with compound **3a**, and their IC50 values were 14.1 and 19.3µM, respectively.

To evaluate the performance of the virtual screening strategy adopted in this study, the rankings and docking scores of the newly identified natural products were retrospectively examined. The distribution of docking scores of the top ranking candidates reserved from Glide calculations were presented in Figure S3. The natural compounds bearing different scaffolds, including **1c**, **2a**, **2b**, **3b,** and **3c**, could be ranked in the top 500 candidates during the Glide SP docking process using both crystal structures 1OSV and 1OSH, whereas the Glide XP results are totally different. As shown in **Table 2**, compounds containing similar chemical skeleton to the crystallographic ligand tend to score higher. Moreover, none of the identified FXR antagonists could be simultaneously ranked in the top 200 candidates by the two Glide XP calculations, confirming the rationality of using two distinct crystal structures for structure-based virtual screening. Despite of the confirmed moderate antagonistic effects against FXR, compounds **3d** and **3e** with smaller sizes were excluded after the initial Glide SP screening, which may be ascribed to the recognized bias of structure-based virtual screening method toward high molecular weight compounds (Pan et al., 2003). The results also demonstrate that 2D molecular similarity search method is a powerful and complementary approach to structurebased virtual screening, which could retrieve biologically active compounds that are regarded as false-negatives by docking simulations.

TABLE 1 | Chemical structures and activities of FXR antagonists and their analogs reported in this studya .

(Continued)

#### TABLE 1 | Continued

aData shown are the average values of triplicate measurements determined by Y2H assays. This system employs the interaction between hFXR-LBD and the coactivator SRC-1. bAttempts to determine IC50 values were made if the inhibition rate at 25µM was larger than 50%.

aCompounds that were ruled out by structure-based virtual screening process but recovered using similarity searching method are labeled with \*.

#### Structural Novelty Assessment

The six natural products were first reported to show FXR antagonistic activity. To evaluate their structural novelty with respect to known FXR antagonists, the pairwise Tc values of chemical similarity were calculated based on the FCFP\_4 fingerprints. The Tc value between the similar compounds **3a** and **3c** is 0.67, a Tc cutoff value of 0.6 was therefore set to define similarity. 15 structurally diverse FXR antagonists including seven natural products (compounds **4**–**21**, Table S1) were collected from literatures, among which the maximum Tc value is 0.46. As shown in **Figure 3**, all the Tc values of the six newly identified hits to the known 15 FXR antagonists were below 0.4, and the maximum Tc values of compounds **1a**, **2a,** and **3a** were 0.38, 0.25, and 0.31 (Table S2), respectively. Accordingly, the six natural products could be considered to be structurally novel as FXR antagonists.

#### Analysis of Predicted Binding Poses

From the structural point of view, the six antagonists can be simply categorized into two classes: terpenes possessing polycyclic skeletons and phenols with chain scaffolds. In order to better delineate SAR, their inactive analogs were also displayed and analyzed in this study.

FIGURE 3 | Heatmap presentation of topological similarities of the six natural products to the 15 reported FXR antagonists.

#### Terpenes

Compounds **1a** and **2a** belong to triterpenes and diterpenes, respectively (Yang et al., 2010), and both of them are isolated from Abies georgei which grows exclusively in China. In a previous study, compound **1a** was reported to have moderate agonistic effect against estrogen receptor (ER). The polycyclic ring skeletons of the two compounds are much similar to that of the bile acids, especially for the tetracyclic triterpene compound **1a**. However, the Tc values between compounds **1a** and **2a** and 6-ECDCA are as low as 0.36 and 0.28, respectively. Intriguingly, despite the opposite activities against FXR, the proposed binding poses of the two compounds closely resemble that of 6-ECDCA when interacting with FXR (**Figure 4**). Similar to 6-ECDCA, compound **1a** adopts cis-orientation in the A/B rings linkage, which is considered to be a unique feature for bile acids. The ring skeleton fits the binding pocket well through favorable VDW

contacts and hydrophobic effects with adjacent residues. The 3α-OH group extends to the space between helix 7 and helix 10/11 and putatively participates in hydrogen bond interactions with residues Tyr358 and His444. Additionally, the terminal carboxyl group could interact with residue Arg328, located at the entrance of the binding pocket, through salt bridge or hydrogen bond interactions. All the interactions described above are beneficial to the binding of compound **1a** to FXR. However, the hydrogen bond formed between the 7α-OH of 6-ECDCA and Tyr366 is absent for compound **1a** due to its structural variation on ring B. Under the physical-shape discrimination mechanism employed by FXR (Mi et al., 2003), bile acids without a 7α-OH group, such as LCA and DCA, showed extremely weak affinity with FXR, which may also cause the moderate antagonistic activity of compound **1a** against FXR.

For the analogs where the 3α-OH groups are replaced by carbonyl groups (compounds **1b** and **1c**), no detectable antagonistic activity was found in the coactivator recruitment assay. The planarity of the double bond may restrict the carbonyl oxygen atom to a position distant from residues Tyr358 and His444, and the absence of corresponding hydrogen bond interactions presumably results in the loss of agonistic effects against FXR.

Compared with compound **1a**, compound **2a** has a relatively smaller volume, but displayed 10-fold stronger antagonistic activity. Apparently, compound **2a** doesn't fit the canonical mechanism that nuclear receptors' antagonists are usually voluminous than agonists (Meyer et al., 2005). In the proposed binding pose with FXR, compound **2a** reveals the same amphipathic properties as the bile acid ligands. The oxygen atom of the carboxyl group putatively forms hydrogen bond interactions with residues Ser329 and Tyr366. At the other end of compound **2a**, hydrogen bond interactions could also form between the hydroxyl group and Arg328. The hydroxyl group seems to be essential for the antagonistic activity of compound **2a**, as its analog with an acetyl group (compound **2b**) exhibited no observed activity against FXR. Owing to the methyl group located at the10α-position, the carboxyl-substituted hydrocarbon ring A protrudes from the benzene ring panel, making VDW contacts and hydrophobic effects with Met325 on helix 5. Besides, the sequential ring structure could interact with loop H1-H2 (Met262), helix 3 (His291 and Met287), and helix 6 (Leu345) by favorable hydrophobic and VDW interactions.

Because of the relatively smaller volume, compound **2a** is not able to extend to the pocket that is occupied by rings A and B of 6-ECDCA, hence no direct interactions with helix 10/11 and helix 12 were observed. Previous studies have suggested that the π-cation interaction between His444 on helix 10/11 and Trp466 on helix 12 plays a critical role for the active conformation of helix 12 induced by endogenous bile acids (Mi et al., 2003; Pellicciari et al., 2005). Steroid agonists with 3α-OH group could facilitate the π-cation interaction by providing appropriate disposition of His444 through the steric restriction of hydrogen bonds formed between 3α-OH and residues His444 and Tyr358. Consequently, without the ability of establishing the triad of Tyr358, His444, and Trp466, compound **2a** couldn't secure helix 12 in the active conformation, thus preventing the recruitment of coactivators.

#### Phenols

Compounds **3a** and **3c**, named daphneone and daphnenone respectively (Zhang et al., 2006), are constituents of Daphne odora Thunb. var. marginata, an ornamental plant whose growth is restricted to the south of China. The two compounds, together with their simple structures and small volume, are extraordinarily peculiar to present antagonistic effects against the CDCAinduced SRC-1 association with FXR. Moreover, it is difficult to find common structural features between the two compounds and known agonists or antagonists. Accordingly, we turned to the initial docking poses to probe the structural basis of the FXR antagonistic profiles of this chemical series (**Figure 5**).

The two compounds were selected from the virtual screening process using the crystal structure of FXR-fexaramine complex. Compound **3a** is sandwiched in the cleft enclosed by helices 5, 7, and 10/11, partially overlapping with the fexaraminebinding pocket. The hydroxyl group points toward helix 5 and TABLE 3 | In silico predicted properties of the six FXR antagonists.

The recommended ranges by QikProp are as follows:

aMolecular weight, 130.0–725.0.

bNumber of hydrogen bond donors, 0.0–6.0.

cNumber of hydrogen bond acceptors, 2.0–20.0.

dNumber of non-trivial rotatable bonds, 0–15.

ePredicted aqueous solubility, −6.5–0.5.

fPredicted octanol/water partition coefficient, −2.0–6.5.

gPredicted apparent Caco-2 cell permeability in nm/sec, <25 poor, >500 great.

hNumber of violations of Lipinski's rule of five, maximum is 4.

iNumber of violations of Jorgensen's rule of three, maximum is 3.

may form hydrogen bonds with residues Ser336 and His298. On the other terminal, the benzene ring moiety extends to the aromatic residues-rich groove formed by Phe288 (helix 7), Trp458 (helix 10/11), and Phe465 (loop H10/11-H12), and probably contacts with these residues by advantageous πstacking interactions. The linker between the two benzene rings fits the pocket by hydrophobic interactions and VDW contacts with surrounding residues such as Leu291, Met294, Ile356, and Ile361.

Previous molecular dynamics simulation studies assumed that the intrinsically unstable loop H10/11-12 controlled the flexibility of helix 12 (Costantino et al., 2005). Through offset face-to-face π-stacking interaction between the benzene ring and Phe465, both compounds **3a** and **3c** could contact with the loop directly, which may interfere with the conformation of the loop and push helix 12 away from its active conformation. For compound **3c** which has a double bond in the linker region, the antagonistic effect in the SRC-1-recruinment assay is slightly weaker. Presumably, the relatively flexible hydrocarbon linker is more suitable for the binding process with FXR-LBD, hence compound **3a** displayed stronger antagonistic activity. When a hydroxyl group was introduced into the hydrocarbon linker, the antagonistic activity of compound **3b** markedly decreased, displaying an IC50 value higher than 25µM. Whereas the much smaller compounds **3d** and **3e**, which share the phenolic moiety with **3a**, exhibited moderate antagonistic effects against FXR. Compounds **3d** and **3e** were docked to FXR by Glide XP mode using the crystal structure 1OSH. As illustrated in Figure S5, the two small molecules occupy merely a fraction of the fexeramine-binding pocket. Hydrophobic effects and shape complementarities presumably dominate the interactions with FXR, as no hydrogen bond was detected in their proposed binding poses.

Notably, compounds **3a** and **3c** have been previously reported to show cytotoxic activities against a variety of human tumor cell lines, including K562, A549, MCF-7, LOVO, HepG2, and A375-S2, with the IC50 values ranging from 3.12 to 51.0µM (Zhang et al., 2006; Wang et al., 2012). Besides, the agonistic profiles of compounds **3d** and **3e** against ER have also been described in a previous study (Cao et al., 2013). The phenol FXR antagonists identified in this study are relatively small, especially for compounds **3d** and **3e**, which probably have effects on other targets in living cells. Further thorough investigations are ongoing to better elucidate the exact mechanisms of action of the newly identified natural FXR antagonists and their implications regarding in vivo pharmacological effects.

#### Druglikeness Evaluation

To assess the drug-like profiles of the six natural products, an in silico prediction of ADME properties was performed using QikProp v4.3 module integrated into Maestro 10.1, and 6- ECDCA was used as the reference compound (**Table 3**). All physically significant descriptors and pharmaceutically relevant properties of the natural FXR antagonists, except for compound **1a**, fall into the recommended ranges of 95% of known drugs, suggesting remarkable potential of druglikeness. The QPlogPo/w and QPlogS values of compound **1a** exceed the limits of either Lipinski's rule of five or Jorgensen's rule of three, therefore the aqueous/lipid solubility should be taken into consideration if further structural optimization was carried out based on the tetracyclic triterpene compound **1a**.

#### CONCLUSION

In summary, we have established a small NPD containing over 4,000 compounds that were previously isolated from about 100 medicinal plants. From the database, six FXR antagonists were identified by strategic virtual screening method, which validated the feasibility of virtual screening to explore the potential targets of natural products. Although procured on the basis of known agonist-binding pocket, two of the most potent compounds **2a** and **3a** could antagonize the CDCA-induced SRC-1 recruitment to FXR-LBD with the IC50 values of 1.29µM and 1.79µM, respectively. The predicted docking mode of the diterpene **2a** against FXR exhibited partially similar binding interactions to those of the crystallographic ligand 6-ECDCA bound to FXR, whereas the daphneone **3a** showed noncanonical proposed binding mode, which may directly contact with the intrinsically unstable loop H10/11-12 by π-stacking interactions with the aromatic residue Phe465. Moreover, as assessed by QikProp, most of the natural FXR antagonists displayed comparable drug-like properties to that of 95% of known drugs. We hope our discovery will provide promising chemical scaffolds for further hit-to-lead optimization and for the study of FXR-related biological mechanisms.

### AUTHOR CONTRIBUTIONS

HL, WZ, and LS conceived the study. YD and SL performed molecular simulations. YD analysed the data and drafted the manuscript. JJ, SZ, and JH carried out experimental studies. HL

#### REFERENCES

reviewed and revised the manuscript. All authors have read and approved the final manuscript.

#### FUNDING

The research is supported in part by the National Key Research and Development Program (Grant 2016YFA0502304), Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under Grant No. U1501501 and the Fundamental Research Funds for the Central Universities. SL is also sponsored by Shanghai Sailing Program (No. 18YF1405100).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2018.00140/full#supplementary-material

orally active agonist of the farnesoid X receptor (FXR). J. Med. Chem. 52, 904–907. doi: 10.1021/jm8014124

receptor FXR. Mol. Cell 11, 1093–1100. doi: 10.1016/S1097-2765(03)00112-6 Nam, S. J., Ko, H., Ju, M. K., Hwang, H., Chin, J. W., Ham, J., et al. (2007). Scalarane sesterterpenes from a marine sponge of the genus spongia and their FXR antagonistic activity. J. Nat. Prod. 70, 1691–1695. doi: 10.1021/np070024k Nevens, F., Andreone, P., Mazzella, G., Strasser, S. I., Bowlus, C., Invernizzi, P.,

et al. (2003). Structural basis for bile acid binding and activation of the nuclear

agonist and induces GLUT4 expression in CHO-K1 cells. J. Steroid Biochem.

Maglich, J. M., Caravella, J. A., Lambert, M. H., Willson, T. M., Moore, J. T., and Ramamurthy, L. (2003). The first completed genome sequence from a teleost fish (Fugu rubripes) adds significant diversity to the nuclear receptor superfamily. Nucleic Acids Res. 31, 4051–4058. doi: 10.1093/nar/gkg444 Makishima, M., Okamoto, A. Y., Repa, J. J., Tu, H., Learned, R. M., Luk, A., et al. (1999). Identification of a nuclear receptor for bile acids. Science 284,

Maloney, P. R., Parks, D. J., Haffner, C. D., Fivush, A. M., Chandra, G., Plunket, K. D., et al. (2000). Identification of a chemical tool for the orphan nuclear receptor

Meyer, U., Costantino, G., Macchiarulo, A., and Pellicciari, R. (2005). Is antagonism of E/Z-guggulsterone at the farnesoid X receptor mediated by a noncanonical binding site? A molecular modeling study. J. Med. Chem. 48,

FXR. J. Med. Chem. 43, 2971–2974. doi: 10.1021/jm0002127

6948–6955. doi: 10.1021/jm0505056

110, 150–156. doi: 10.1016/j.jsbmb.2008.03.028

	- Wang, L. B., Dong, N. W., Wu, Z. H., and Wu, L. J. (2012). Two new compounds with cytotoxic activity on the human melanoma A375-S2 cells from Daphne giraldii callus cells. J. Asian Nat. Prod. Res. 14, 1020–1026. doi: 10.1080/10286020.2012.701206
	- Wang, Y. D., Chen, W. D., Moore, D. D., and Huang, W. (2008). FXR: a metabolic regulator and cell protector. Cell Res. 18, 1087–1095. doi: 10.1038/cr.2008.289