Article pubs.acs.org/jcim
Ensemble-Based Virtual Screening Led to the Discovery of New Classes of Potent Tyrosinase Inhibitors Joonhyeok Choi,† Kwang-Eun Choi,† Sung Jean Park,‡ Sun Yeou Kim,‡ and Jun-Goo Jee*,† †
Research Institute of Pharmaceutical Sciences, College of Pharmacy, Kyungpook National University, 80 Daehak-ro, Buk-gu, Daegu 702-701, Republic of Korea ‡ College of Pharmacy, Gachon University, Incheon 406-799, Republic of Korea S Supporting Information *
ABSTRACT: In this study, we report new classes of potent tyrosinase inhibitors identified by enhanced structure-based virtual screening prediction; the enzyme and melanin content assays were also confirmed. Tyrosinase, a type-3 copper protein, participates in two distinct reactions, hydroxylation of tyrosine to DOPA and conversion of DOPA to dopaquinone, in melanin biosynthesis. Although numerous inhibitors of this reaction have been reported, there is a lag in the discovery of the new functional moieties. In order to improve the performance of virtual screening, we first produced an ensemble of 10,000 structures using molecular dynamics simulation. Quantum mechanical calculation was used to determine the partial charges of catalytic copper ions based on the met and deoxy states. Second, we selected a structure showing an optimal receiver operating characteristic (ROC) curve with known direct binders and their physicochemically matched decoys. The structure revealed more than 10-fold higher enrichment at 1% of the ROC curve than those observed in X-ray structures. Third, high-throughput virtual screening with DOCK 3.6 was performed using a library consisting of approximately 400,000 small molecules derived from the ZINC database. Fourth, we obtained the top 60 molecules and tested their inhibition of mushroom tyrosinase. The extended assays included 21 analogs of the 21 initial hits to test their inhibition properties. Here, the moieties of tetrazole and triazole were identified as new binding cores interacting with the dicopper catalytic center. All 42 inhibitors showed inhibitory constant, Ki, values ranging from 11.1 nM and 33.4 μM, with a tetrazole compound exhibiting the strongest activity. Among the 42 molecules, five displayed more than 30% reduction in melanin production when treated in B16F10 melanoma cells; cell viability was >90% at 20 μM. Particularly, a thiosemicarbazone-containing compound reduced melanin content by 55%.
■
INTRODUCTION Tyrosinase (EC 1.14.18.1) is an oxidase that catalyzes two distinct reactions: hydroxylation of monophenol and conversion of diphenol to the corresponding quinone. The cognate substrate, tyrosine, is converted into DOPA and further into dopaquinone by the enzyme. A related enzyme, catechol oxidase, can catalyze only the latter reaction. Melanin produced in melanocytes of the skin primarily determines the skin color of humans. The colors of the eyes and hair as well as browning in food are associated with melanin. Since dopaquinone eventually leads to the production of melanin via eumelanin, impaired tyrosinase function causes disorders such as albinism, where very little melanin is synthesized in the body. In contrast, if tyrosinase is overactive, melanin synthesis increases, resulting © XXXX American Chemical Society
in skin disorders. Tyrosinase activity reportedly contributes to cancer and neurodegeneration associated with Parkinson’s disease.1−3 Numerous studies have examined selective tyrosinase inhibitors in order to treat skin disorders or develop skin-lightening cosmetics.4−11 Tyrosinase contains two copper ions, CuA and CuB, at its catalytic center. A copper atom is chelated by three neighboring Nε2 atoms of the histidine residues, classified as a type-3 copper protein. In addition to tyrosinase and catechol oxidase, type-3 copper proteins contain oxygen-carrying hemocyanin. All six histidines that chelate the two coppers are conserved Received: August 5, 2015
A
DOI: 10.1021/acs.jcim.5b00484 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Information and Modeling
There is no universal algorithm that is accurate in all cases, and identifying which is best for specific circumstances must be determined a priori.34−40 Similar difficulties may arise when applying SBVS for identifying inhibitors of tyrosinase. Tyrosinase shows conformational flexibility around the catalytic region,28 suggesting that structural changes may occur upon the binding of inhibitors. Additionally, the oxidation states of catalytic copper ions must be described, which is critical for calculating the interaction energy between the protein and inhibitor using docking software. We report the structure ensemble-based docking in combination with optimization of the structure and software. It can be used to overcome difficulties and to substantially improve VS performance.
among type-3 copper proteins, indicating evolutionary linkage of the proteins. The copper ions shuffle between four possible oxidation states, the oxy, met, deoxy, and deact states and the catalytic cycle is coupled with these states.12 Quantum chemical calculations and spectroscopic studies have characterized the copper states in detail using simplified systems.13 In the oxyform, two additional oxygen atoms form a planar structure with the two coppers, [Cu(II)−O22−−Cu(II)], where oxidation of the substrate is triggered. The met- and deoxy-forms exist as the dicupric ([Cu(II)−Cu(II)]) and dicuprous ([Cu(I)−Cu(I)]) states, respectively. In the met-form, which is assumed to be the resting enzymatic state, one or two hydroxide molecules bridge two Cu(II) atoms. In contrast, the coordination between copper and histidines is thought disrupted in the deact states.12 In vitro enzyme and cell assays have identified numerous natural and synthetic compounds as inhibitors.4−11 The inhibitors exemplify natural compounds containing polyphenols. Arbutin, a glycosylated hydroquinone found in plants, is a representative inhibitor that is used as a skin-lightening agent. In addition to the identification of inhibitors, numerous studies have proposed 3D models of tyrosinase bound to inhibitors through docking simulation, supporting the rationale for new functional moieties.9,14−25 A few complex structures, one between tropolone and mushroom tyrosinase (PDB ID: 2Y9X)26 and another between phenylthiourea and sweet potato catechol oxidase (PDB ID: 1BUG),27 have been reported, where the inhibitors interacted with the proteins in met or deoxy states but not in oxy states. In contrast, in the complex structure involving kojic acid and tyrosinase from Bacillus megaterium (PDB ID: 3NQ1),28 the inhibitor was located 7 Å away from the active site, leaving the understanding of the inhibitory mechanism unclarified. Crystal structures between substrates, tyrosine and DOPA, and tyrosinase from Bacillus megaterium were recently determined after soaking with zinc rather than with copper atoms.29 The structures revealed subtle changes in the orientations for the oxidation of monophenol and catechol, providing a valuable clue for determining the common and distinct features between tyrosinase and catechol oxidase. In this study, we report the use of structure-based virtual screening (SBVS) for inhibitors of tyrosinase. SBVS for identifying bioactive small molecules has advanced with the development of algorithms and computational capacity. Once candidates are narrowed down to several tens or hundreds of molecules, researchers can then experimentally confirm the binding or catalytic activities toward target biomolecules. Despite these improvements, the number of false-positives in SBVS is still considerable. Zhu et al. reported that the percentage of hit molecules among purchased candidate molecules in virtual screening (VS) is approximately 13% (median), based on a literature survey.30 This value was calculated with the results from both ligand-based virtual screening (LBVS) and SBVS. Since LBVS generally outperforms SBVS in hit identification,31,32 the hit rate using only SBVS may be lower. Furthermore, it should be noted that only favorable results are reported in publications. However, SBVS has great potential in identifying the inhibitors of novel scaffolds, accounting for the sustained effort for SBVS improvement.33 Two hurdles for SBVS to overcome for better performance are the conformational changes of receptors upon binding to ligands and the limitations of the docking algorithm. Current docking algorithms require additional optimization to accurately reflect conformational changes upon ligand binding.
■
MATERIALS AND METHODS Generation of Ensemble Structures by Molecular Dynamics Simulation. The structure of mushroom tyrosinase (PDB ID: 2Y9X)26 was used as a template for molecular dynamics (MD) simulation. Gaussian 0341 generated partial charges for Cu(II) and the coordinating atoms of three histidine residues guided by the R.E.D. tools42 with the option of RESP-A1 at the level of HF/6-31G*. The atomic partial charges were converted into the AMBER force field. To generate an ensemble of structures, 100 ns MD simulation was performed using the AMBER package (ver. 12).43 Once the protonation states of ionizable side chains were determined using the PDB2PQR server,44 the coordinate was solvated with TIP3P water molecules in a periodic truncated octahedron box such that its walls were at least 10 Å away from the solute. MD simulation was performed using the PMEMD module with the ff99SB-ILDN force field.45 All bonds involving hydrogen atoms were constrained by SHAKE, permitting a time step of 2 fs. Nonbonding interactions were truncated at 10 Å, and the particle mesh Ewald method was used to calculate electrostatic interactions under periodic boundary conditions. The simulations consisted of five stages: 1000 steps minimization, 50 ps run for heating from 0.1 to 300 K, 50 ps run under a constant pressure of 1 atm and temperature of 300 K, 500 ps run for equilibration, and 105 ns run for production. Positional restraints were applied for the first three stages. Unrestrained equilibration and production runs were carried out at constant temperature (300 K) using a Langevin thermostat with a collision frequency of 2 ps−1 and constant pressure (1 atm) using a Berendsen barostat with a pressure relaxation time of 2 ps. Trajectory at every 10 ps for the later 100 ns of the production run was saved, generating 10,000 data points in total. Selection of Template Structure by Receiver Operating Characteristic Curves with Direct Binders and Decoys. As an engine for SBVS, DOCK 3.646,47 was employed by comparing docked poses of tropolone and phenylthiourea with those by the other software, AutoDock (ver 4.2.6),48 AutoDock Vina (ver 1.1.2),49 and DOCK 6.6.50 The strongest 46 direct binders were selected based on half-maximal inhibitory concentration (IC50) and inhibitory constant (Ki) values from BindingDB.51 Then DUD-E server52 generated their physicochemically matched 2505 decoys, resulting in 2551 molecules in total. Once aligning each structure from MDderived ensemble to the reference one by TM-align,53 SBVSs with the ligands and decoys followed. All the molecules were docked into 10,000 structures. As criteria to choose the best structure for the VS, receiver operating characteristic (ROC) B
DOI: 10.1021/acs.jcim.5b00484 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Information and Modeling
rdkit.org). Once digitizing functional moieties in a molecule with the fingerprints, Tc calculates the portion of common over union elements. The Tc values were between 0 and 1, where 0 and 1 correspond to no and perfect overlap of two chemicals, respectively. The 483 known direct binders from all species of tyrosinase deposited in the BindingDB database (ver. 2015m6)51 were extracted to quantify the similarity between a test molecule and known tyrosinase inhibitors. The ChEMBL database (ver. 20)89,90 was used to search the bioactivities of the inhibitors in this study. In-house written scripts, Automated LIgand Search for PolyPharmacology (ALIS-PP), automated the processes of cheminformatics.
curves were used. All molecules in a docking simulation set were sorted in ascending order in terms of calculated intermolecular energies. The accumulated true-positive and false-positive rates of each molecule from the top were calculated and assigned as y and x coordinates in the ROC curve. The best structure was judged based on three metrics: area under the curve (AUC), logarithmically scaled AUC (LogAUC),46 and enrichment factor at 1% (EF1).54 All procedures were automated using a written in-house protocol, Automated pLatform for Integrative Structure-based DOCKing (ALIS-DOCK), to simultaneously use Linux-cluster consisting of 120 Xeon E5−2670 cores. Structure-Based High-Throughput Virtual Screening. ALIS-DOCK implements the DOCK 3.646 screened ZINC database55 of the ChemBridge virtual library. Only representative molecules in terms of chemical similarities were used. The ZINC database arranges the library into a subset in which no molecule shares similarity higher than 0.9 Tanimoto coefficient (Tc) with other molecules. After filtering out molecules with undesirable physicochemical properties as drug candidates based on Lipinski’s rule of five,56 approximately 400,000 small molecules were chosen. The ZINC-deposited ligand formats (*.ddb) were extracted for DOCK 3.6. The Flexibase format contained a ligand desolvation scoring term.46 Enzyme Activity Assays with inhibitors. All reagents in this study except inhibitor candidates were purchased from Sigma-Aldrich (St. Louis, MO, U.S.A.) or Tokyo Chemical Industry (Tokyo, Japan). Inhibitors identified by HTVS were purchased from ChemBridge (San Diego, CA, U.S.A.). Final reaction mixtures for enzymatic assays contained mushroom tyrosinase, phosphate-buffered saline (PBS), 625 μM substrate (L-tyrosine), and inhibitors. The reaction mixture containing the enzyme and inhibitor was incubated at 30 °C for 30 min, and the change in absorbance at 475 nm after adding substrate was measured over time. All solutions for measuring enzyme activity contained 5% dimethyl sulfoxide (DMSO) to solubilize the organic molecules. After confirming inhibition at single concentration of 50 μM, the activities were quantified by calculating the IC50 with a series of inhibitors at different concentrations. Kinetics assays using substrates at varying concentrations generated Km values, enabling the conversion of IC50 into Ki.57 All fittings in this study were conducted in MATLAB (MathWorks, Natick, MA, U.S.A.). MTT and Melanin Content Assays. To assay cell viability, the strength of violet staining in B16F10 cells stained with 3(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) were examined. Cells at a density of 1 × 104 were treated with MTT for 24 h. The medium was removed from each well, and the cells were washed with PBS. The cells were incubated with 200 μL fresh media and 10 μL MTT (5 mg/mL in PBS). After removing the medium, 200 μL DMSO was added. After 30 min, UV absorption was measured at 570 nm. For the melanin content assay, the amount of melanin released from the cells was measured. B16F10 cells were incubated at a density of 1 × 104 for 24 h. Next, 100 μM of 3-isobutyl-1methylxanthine (IBMX) was combined with inhibitors and incubated for 72 h. UV absorptions were measured at 405 nm with 100 μL of medium using an ELISA reader. Reduced melanin production was expressed as a percentage of untreated controls. Cheminformatics. For cheminformatics studies, the chemical similarity was quantified using Tc with the Morgan circular fingerprints implemented in the RDKit (http://www.
■
RESULTS Selection of Docking Algorithm for High-Throughput Virtual Screening. The apo structure of mushroom tyrosinase contains a molecular oxygen that bridges two copper atoms (PDB ID: 2Y9W).26 The exact position of the oxygen is occupied by an oxygen atom of tropolone in complex (2Y9X). Thus, the area should be available for inhibitors to access this site in SBVS. Therefore, we prepared a structure of tyrosinase that did not contain any molecular oxygen for structure-based docking. We tested which docking algorithm could accurately reproduce the poses found in the X-ray structures by using two inhibitors found in complex structurestropolone and phenylthiourea. According to the crystal structures, the oxygen of tropolone and the sulfur of phenylthiourea are the functional groups responsible for binding to the type-3 copper protein.26,27 Among the academically available software programs AutoDock,48 AutoDock Vina,49 DOCK 3.6,46 and DOCK 6.6,50 only DOCK 3.6 could place both the oxygen and the sulfur at the correct positions observed in the complexes (Figure S1). In contrast, DOCK 6.6 accurately reproduced only the pose of tropolone. These programs required the definition and optimization of the copper parameters for docking. Although the values for the parameters were varied in AutoDock, AutoDock Vina, and DOCK 6.6, the outcomes were still unsatisfactory. The functional groups recognizing the coppers were inappropriately placed. The better performance of DOCK 3.6 is likely associated with its improvement of docking with the metalloenzymes.47 The DUD-E database includes 10 metalloenzymes.52 Of them, six proteins (ada17, comt, def, fpps, hdac2, and hdac8) showed excellent early enrichments of true-positives (LogAUC > 20), again reflecting the reliable performance of DOCK 3.6 in identifying metalloenzyme inhibitors. Therefore, all subsequent dockings were conducted using DOCK 3.6. We also identified inappropriate positions from the viewpoints of metal recognition in the complex models that have been proposed with new inhibitors. For instance, kojic acid,14 arbutin,16 and phenylthiourea analog58 placed their aromatic rings more closely to the catalytic coppers than to the hydroxyl or thione groups in the models. The models were inconsistent with the experimental complex structures. This may be because of the imperfect handling of metal parameters in the software. In contrast, some models that included two oxygen atoms bridging the two coppers in addition to the inhibitors have been proposed.22,59 However, because the surface-exposed bridging oxygen atoms were replaced by the functional groups of ligands in the available complex structures,26,27,29 it will be more reasonable to omit oxygen for docking studies. Selection of Best Structure from Trajectories by Molecular Dynamics Simulation for High-Throughput C
DOI: 10.1021/acs.jcim.5b00484 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Information and Modeling
Figure 1. Receiver operating characteristic (ROC) curves and overlay of selected MD and X-ray structure. (A) ROC curves for initial X-ray structures (2Y9X and 2Y9W, blue) and final structure from molecular dynamics (MD) simulation (red). The values of three metrics, area under the curve (AUC), logarithmically scaled AUC (LogAUC),46 and enrichment factor at 1% (EF1)54 of ROC curve, are shown as well. “Optimized” indicates the selected molecule from MD trajectories. (B) Overlaid structures of the initial X-ray (blue) and the optimized structures (red). The shaded area indicates the substrate binding catalytic core. The root-mean-square deviation of Cα atoms in the two structures was 1.7 Å.
Table 1. Details of the First-Round 21 Hitsa ID
rankb
ZINC ID
MW
LogP
Ki (μM)
MTT (%)
Δmelanin content (%)c
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
5 10 18 19 21 27 28 30 31 34 36 41 43 46 47 49 54 56 57 59 60
ZINC65507416 ZINC67446619 ZINC00318008 ZINC65422151 ZINC04754845 ZINC77471033 ZINC20447630 ZINC15007657 ZINC65501173 ZINC65404460 ZINC04638102 ZINC20255456 ZINC77440228 ZINC01154498 ZINC55108434 ZINC00470926 ZINC00184906 ZINC06596838 ZINC17946040 ZINC72423400 ZINC19781585
347 345 294 377 325 353 250 274 346 329 346 279 340 385 313 302 230 335 255 388 380
3.02 2.69 3.71 3.06 2.46 2.20 2.21 1.63 3.38 2.99 4.05 3.67 1.86 2.74 1.06 1.62 2.41 1.38 3.03 2.05 3.46
3.0 ± 0.8 3.6 ± 0.2 2.6 ± 0.3 940 ± 56 (nM) 3.0 ± 0.2 3.4 ± 0.2 2.4 ± 0.1 29 ± 5 (nM) 3.3 ± 0.1 3.8 ± 0.4 33.4 ± 8.9 4.8 ± 0.5 4.5 ± 0.5 5.6 ± 0.5 908 ± 290 (nM) 8.3 ± 1.2 101 ± 13(nM) 7.5 ± 0.8 705 ± 42 (nM) 5.0 ± 0.5 11.3 ± 1.2
>90 >90 >90 64 ± 6 >90 >90 >90 61 ± 6 50 ± 4 >90 >90 >90 >90 >90 67 ± 8 >90 >90 61 ± 5 >90 >90 84 ± 14
90 >90 >90 >90 >90 >90 >90 >90
55 ± 3*
