{"paper_id":"325a6ce5-ef70-49d2-9e57-02f2cf3b90fc","body_text":"Human aromatase, also known as human CYP19A1,\nis a member of the\ncytochrome P450s (CYPs), a superfamily of monooxygenases that contain\nheme as a prosthetic group and play essential roles in drug metabolism\nand homeostasis. 1  Aromatase is responsible\nfor the last step of the endogenous biosynthesis of estrogens in both\nmen and women, as it catalyzes the aromatization of androgens (androstenedione/testosterone)\nto estrogens (estrone/17β-estradiol). 2  Estrogens are a group of steroid hormones that, in addition to their\nphysiological function in the development and regulation of the female\nreproductive system and secondary sex characteristics, also play a\ncrucial role in the development and progression of breast cancer (the\nmost common cancer among women worldwide) by promoting cell division\nand growth in breast tissue. Therefore, aromatase has been targeted\nfor the treatment of ER+ (estrogen receptor positive) breast cancer,\nespecially in postmenopausal women. 3  In\naddition to breast cancer, other hormone-related conditions, such\nas short stature in boys, 4  infertility\nin men, 5  and women, 6  endometriosis, 7  leiomyomatosis, 8  and Klinefelter’s syndrome 9  are increasingly being treated by administering aromatase\ninhibitors. Several generations of aromatase inhibitors have been\nmarketed and widely employed in clinical practice, with the third\ngeneration of aromatase inhibitors ( Figure  1 ) being nowadays the most commonly used.\nCommonly\nprescribed third-generation aromatase inhibitors.\nAromatase inhibitors are commonly divided into\ntwo types: nonsteroidal\naromatase inhibitors (NSAIs; type 1), such as anastrozole and letrozole,\nand steroidal aromatase inhibitors (SAIs; type 2), such as exemestane. 10  While pharmacological therapies with AIs are\nuseful adjuvants in treating the conditions mentioned above, there\nis increasing evidence that systemic suppression of aromatase can\nlead to adverse health effects, including musculoskeletal symptoms, 11  cognitive dysfunction, and other neurological\nsymptoms. 12  In light of these challenges,\nit would be highly beneficial to identify new inhibitors with significantly\nlower permeation of the blood-brain barrier (BBB). 13  In addition to the adverse effects that have been previously\ndiscussed, another problematic issue that arises during aromatase\ninhibitor administration is drug resistance. This issue can be divided\ninto two categories: innate resistance and acquired drug resistance.\nThe latter can occur as a result of continued treatment with aromatase\ninhibitors in breast cancer patients. 14  Consequently, there is a need for novel aromatase inhibitors, designed\nto mitigate side effects or drug resistance and tailored to individual\npatient profiles. Furthermore, these inhibitors should demonstrate\nthe potential for use in conjunction with various adjunctive therapies,\naiming to improve tolerability and reduce the incidence of side effects\ncompared to current treatment protocols.\nComputer-aided drug\ndesign (CADD) strategies, specifically tailored\nto the discovery of CYP19A1 inhibitors, have been extensively utilized\nto streamline drug development. High-quality structures of CYP19A1\nin complex with substrates or inhibitors have been resolved down to\n2.75 Å resolution. 15 − 17  These structures provide crucial\ninsights into key binding residues and their roles in catalysis. Molecular\ndocking techniques have proven invaluable in elucidating structure–activity\nrelationships (SARs) for potential CYP19A1 inhibitors, as demonstrated\nin studies of sulfonamide derivatives, where docking was used to predict\nSARs as CYP19A1 inhibitors. 18  Furthermore,\ndocking has been utilized to expedite the identification of aromatase\ninhibitors, facilitating the ranking of compounds based on their predicted\nbinding affinities and interactions with the active site. 19  Additionally, molecular dynamics (MD) simulations\nhave been employed not only for binding energy calculations 20  but also to analyze the SAR and stability of\nCYP19A1 inhibitor-protein complexes, offering insights into how dynamic\nprotein–ligand interactions influence inhibitor potency. 21\nIn addition to structure-based approaches,\nmachine learning models\nare particularly useful as ligand-based approaches, as they permit\nthe prediction of aromatase inhibitors from a range of scaffolds,\nincluding nonsteroidal and steroidal inhibitors. Although some research\nhas investigated the potential of machine learning in predicting CYP19A1\ninhibitors, 22 − 24  incorporating such models in structure-based virtual\nscreening pipelines remains an underexplored area.\nThe development\nof synthetic biology has provided powerful tools\nfor the expression and application of heterologous proteins. 25  Heterologous proteins generally have the characteristics\nof low cost, good specificity, suitability for analysis in complex\nenvironments, and simple operation; based on these advantages, heterologous\nproteins have broad application prospects in clinical diagnosis, food\nand drug analysis, and research on biomaterial. 26  The CYP19A1 inhibitors screening and evaluation method\nbased on heterologous expression of human CYP19A1 in yeast was adopted\nin this study to validate candidate inhibitors predicted by our machine\nlearning model. 27  This method allows the\nscreening of a large number of candidate chemicals faster and more\neconomically than other methods. 28 − 30\nIn our work, we\ntrained machine learning models on data from the\nChEMBL database 31  and PubChem 32  BioAssay to predict the probability of small\nmolecules to inhibit CYP19A1. These models were applied for virtual\nscreening in combination with additional filtering strategies, including\ndocking and visual inspection. We developed classifiers optimized\nfor early enrichment, making them more suitable for virtual screening\nthan previous classification models for CYP19A1 inhibition. We utilized\nan enzymic assay based on heterologous protein to evaluate the activity\nof the proposed CYP19A1 inhibitors identified through  in silico  modeling in our study. The screening process yielded several active\ncompounds with previously unknown scaffolds, among which the most\npromising hit was then subjected to further investigation through\nstructure-based modeling and MD simulations to understand its binding\ninteractions better.\n\nTo build predictive models for aromatase (CYP19A1) inhibition,\nreliable sources of CYP19A1 inhibition activity data are essential.\nThe ChEMBL 33  database is a curated resource\nthat compiles bioactivity data from scientific literature and assays,\nmaking it a key source of activity data in drug discovery. PubChem 32  BioAssay is a comprehensive repository of bioactivity\nresults from high-throughput screening and other experimental studies.\nWe compiled data on the inhibition of human CYP19A1 by drug-like small\nmolecules from the ChEMBL database and the PubChem BioAssay database,\nadhering to the protocol described in the  Methods  section. For the data set extracted from the ChEMBL database, compounds\nwith pChEMBL values ≥6 were labeled as active compounds, and\ncompounds with pChEMBL values <6 were labeled inactive. For data\nextracted from the PubChem BioAssay database, compounds with a PUBCHEM_ACTIVITY_SCORE\nof 0 were labeled as inactive, compounds with scores between 40 and\n100 were labeled as active, and compounds with a score outside this\nrange were excluded from the study. The final data sets were composed\nof several hundred compounds, as summarized in  Table  1 .\nThe Venn diagrams in  Figure  2  show minimal overlap between\nthe curated data sets derived\nfrom the ChEMBL database and PubChem BioAssay database. Only five\naromatase inhibitors and seven inactive compounds are shared between\nthe two data sets. The small overlap underscores the substantial differences\nbetween the data sets, which could contribute to distinct model preferences\nand performance when trained on each data set.\nVenn diagrams illustrating\nthe overlaps between the curated (a)\nactive and (b) inactive data sets derived from the ChEMBL database\n(turquoise) and PubChem Bioassay (pink).\nTo better understand the chemical space covered\nby the data sets,\nwe performed Principal Component Analyses (PCAs) based on 13 physicochemical\nproperties (see the  Methods  section) of the\ncompounds and molecular similarity comparisons. The PCA plots shown\nin  Figure  3 a,b put\nthe two data sets and the molecular diversity of approved drugs 36  into perspective. From these plots, we conclude\nthat the ChEMBL-derived data set is more focused on the chemical space\nmost densely populated by the approved drugs, whereas the PubChem\nBioAssay-derived data set spreads more broadly across the chemical\nspace of approved drugs. For both data sets, an accumulation of bioactive\ncompounds in the area most populated by the approved drugs is apparent.\nThe pairwise Tanimoto similarities of compounds in the two data sets,\ncalculated using ECFP4 fingerprints, are compared in  Figure  3 c,d. The plots indicate that\nthe active and inactive compounds from the ChEMBL and PubChem data\nsets are of moderate similarity. We further analyzed the intragroup\nsimilarity of active and inactive compounds within the ChEMBL and\nPubChem BioAssay-derived data sets. Although the PCA did not reveal\nclear distinctions in chemical features between active and inactive\ncompounds, particularly within the drug-like range, Tanimoto similarity\ncomparisons ( Figure  3 e,f) indicate lower similarity between the two classes, suggesting\nthat machine learning models could effectively classify the compounds.\nThe bioactive compounds in the ChEMBL-derived data set include 313\nMurcko scaffolds, 37  while the inactive\ncompounds encompass 294 distinct scaffolds. The PubChem data set represents\n146 scaffolds among the active compounds and 1,166 scaffolds within\nthe inactive set. The most prominent Murcko scaffolds’ distribution\nwithin the ChEMBL and PubChem Bioassay-derived data sets are illustrated\nin  Figure  3 g,h. In\nthe ChEMBL-derived data set, 26 Murcko scaffolds are represented by\nat least 10 compounds, while in the PubChem data set, 11 Murcko scaffolds\nmeet this criterion ( Table S1 ). The two\nbar charts further demonstrate that compounds with the same Murcko\nscaffold often have different activity classifications across data\nsets, especially for the ChEMBL data set. This allows the model to\nlearn the importance of various derivatives beyond the primary scaffolds.\nChemical\nspace analysis of the ChEMBL and PubChem BioAssay-derived\ndata sets. (a, b) PCA scatter plots comparing the chemical space of\nactive and inactive compounds in the (a) ChEMBL derived data set and\n(b) PubChem BioAssay data set with approved drugs from the DrugBank\ndata set as a background. (c–f) Proportion of active or inactive\ncompounds in the ChEMBL and PubChem BioAssay data sets at a given\nminimum similarity of nearest neighbors in the other group, based\non the Tanimoto coefficient calculated from ECFP4 (2048 bits). (g)\nNumber of molecules in the ChEMBL data set represented by the most\nprominent Murcko scaffolds. (h) Number of molecules in the PubChem\ndata set represented by the most prominent Murcko scaffolds. The structures\nand further information for the scaffolds listed in (g, h) are listed\nin  Table S1 . (i) PCA scatter plot comparing\nthe chemical space of the experimentally tested compounds in this\nstudy with approved drugs from the DrugBank data set as a background.\nTo\nbuild predictive models for aromatase inhibition, we trained four\nrandom forest (RF) classifiers selectively using the two data sets\n(ChEMBL and PubChem BioAssay) separately, with two types of molecular\ndescriptors (208 physicochemical descriptors) and ECFP4 fingerprints\n(2048 bits), following the steps of cross-validation (CV), hyperparameter\noptimization, and model validation.\nBefore model building, the\ndata sets were each split into a training set and a test set in ratio\nof 4:1. We performed 5-fold CV on the training set with similarity-based\nsplitting and oversampling (see the  Methods  section:  Model development and validation ) with\ndifferent hyperparameters, and recorded the best Boltzmann-enhanced\ndiscrimination of receiver operating characteristic (BEDROC) scores 38 —a metric used to evaluate the effectiveness\nof virtual screening methods in ranking active compounds early in\na ranked list—of the average performance in 5-fold CV models\nand the corresponding hyperparameter (“max_features”).\nAs shown in  Table  2 , the best model we obtained from the two data sets and descriptor\ncombinations obtained BEDROC scores of 0.35 to 1.00. The high BEDROC\nscores of ChEMBL-based models indicate superior early enrichment.\nThe PubChem-based models showed lower BEDROC values, likely due to\nthe broader chemical diversity ( Figure  3 a,b) and data imbalance within the PubChem test set.\nMaximum fraction of features considered\nper split.\nThe best hyperparameters identified through CV were\nused to train\nRF classifiers on the training sets, with oversampling applied to\naddress the class imbalance. The RF classifiers were evaluated on\nthe test set using receiver operating characteristic (ROC) curves\n( Figure  4 a,b), with\nperformance measured as the area under the ROC curve (ROC-AUC) and\nthe BEDROC score ( Table  3 ). The ROC-AUC values obtained by the four RF classifiers were between\n0.69 and 0.86, whereas their BEDROC values ranged from 0.29 to 1.00.\nThe BEDROC values revealed a discrepancy: while the ChEMBL-derived\nmodels performed well with strong early enrichment, the PubChem-derived\nmodels exhibited a lack of early enrichment ability. This might be\ndue to artifacts arising from data imbalance and the smaller test\nset size of the PubChem data set (33 active data points), as BEDROC\nplaces significant weight on the top-ranked data points; in this case,\nthe few top-ranked molecules dominate the BEDROC value.\nReceiver operating\ncharacteristic (ROC) comparison of (a) ChEMBL-based\nmodels evaluated on ChEMBL test set; (b) PubChem-based models evaluated\non PubChem test set; (c) ChEMBL-based final models evaluated on PubChem\ndata set; and (d) PubChem-based final models evaluated on ChEMBL data\nset.\nTo maximize the use of the available data resources\nfor the virtual\nscreening campaign, the test sets were combined with the training\nsets to create a full-size data set for training the final models.\nThis approach ensures that all available data contribute to model\ndevelopment, enhancing the predictive power and robustness of the\nresulting models. We also cross-validated the two data sets with the\ncorresponding models ( Figure  4 c,d and  Table  3 ). The ROC-AUC values were lower in this testing scenario, ranging\nfrom 0.60 to 0.64. The drop in performance compared to the test with\ndata of the same origin is expected and consistent with common observations\nand reports. One important reason is that the two data sets cover\ndistinct areas of the (drug-like) chemical spaces. Another major factor\nis differences in data annotation and data heterogeneity in general.\nCertainly, some—although not all—of these discrepancies\nwere mitigated by following a binary classification rather than a\nregression approach.\nInterestingly, the PubChem_ECFP4 model\nachieved a high BEDROC score\nwhen tested on the ChEMBL data set, which contrasts with its performance\nwhen trained on 80% of the PubChem data set and tested on the remaining\n20%. This result highlights the model’s capacity to generalize\nacross data sets despite its relatively low performance during internal\ncross-validation with the PubChem data set. This result led to the\nhypothesis that the active and inactive rate of data in the ChEMBL\ndata set is more balanced than the PubChem test set, making it more\nsuitable for validation with BEDROC. In addition, the PubChem_ECFP4\nmodel may show bias against novel scaffolds from a different chemical\nspace, leading to distinct performance when tested across data sets.\nTo identify potential aromatase inhibitors, we employed the four\ntrained models in a virtual screening campaign, prioritizing compounds\nfor experimental validation based on their predicted activity, with\nthe workflow shown in  Figure  5 . For this purpose, we used the MolPort database ( www.molport.com ), which contains\napproximately 4.6 million molecules. After applying a washing and\ncanonicalization process (as detailed in the  Methods  section), a total of 4,310,620 molecules were retained for further\nanalysis. The database was preprocessed by filtering for compounds\nthat obey the rule of five (ro5) 39  and\nthe presence of substructures common to CYP inhibitors (imidazole,\ntriazole, or pyridine), resulting in 1,434,904 molecules for screening.\nThese molecules were then evaluated using the RF classifiers, which\npredicted the probability of a compound being a CYP19A1 inhibitor.\nMolecules ranked in the top 10,000 by each model and with a probability\ngreater than 0.7 (20,444 molecules) were kept for further analysis.\nTo ensure novelty, compounds with a similarity score greater than\n0.7 to known active compounds in the ChEMBL and PubChem bioassay data\nsets were filtered out, narrowing the selection to 20,317 molecules.\nFurther refinement involved the removal of undesirable substructures\nand properties, such as Pan-Assay Interference Compounds (PAINS), 40  long alkanes and molecules with more than eight\nrotatable bonds, which reduced the pool to 11,527 molecules, followed\nby additional substructure filtering to avoid steric hindrance for\nCYP inhibition moieties, which yielded 6182 molecules. These molecules\nwere subjected to similarity-based clustering followed by visual inspection\nto ensure the selection of chemically diverse hits, ultimately reducing\nthe number to 1503. A final round of docking and visual inspection\nconsidering steric hindrance, shape complementary, and feasibility\nto reach the heme iron in the binding site led to the final selection\nof 10 promising compounds for experimental testing. The structures\nof the selected compounds are listed in  Table  4 , together with the structure of the most\nsimilar compound in the respective training set. The chemical space\nof the tested molecules was analyzed with the same dimensionality\nreduction using PCA as previously described, as shown in  Figure  3 i. The data points\nlocalized within the drug-like range, and also in the similar range\nof active data points of both data sets.\nSchematic workflow of\nthe virtual screening campaign, including\nthe number of compounds remaining after each step.\nWith\nadditional information on the\nTanimoto similarity compared with the inhibitor candidate, and activity\nlabel of the compound in the data set.\nAccording\nto a SciFinder search.\nThe measured activity reported by literature was listed. 19 , 41 − 44  N.A.: not available.\nTwo yeast strains were employed for\nscreening and evaluating candidate compounds against the human CYP19A1\nprotein. One strain, engineered to express human CYP19A1 and cytochrome\nP450 reductase (CPR), served as an enzyme system (“enzyme bag”)\ncapable of converting testosterone to β-estradiol. In the other\nstrain, the β-estradiol then acts as an inducer, displacing\nan Hsp90 chaperone complex from a chimeric transcription factor that\npasses through the nuclear pore and binds to short lex operators.\nThis binding enhances the expression of the reporter protein yEGFP\n(yeast-enhanced green fluorescent protein), as illustrated in  Figure  6 a. The corresponding\nfluorescence signals can then be monitored by flow cytometry. When\ninhibitors are added together with the aromatase substrate testosterone,\nthe enzymatic activity of human CYP19A1 decreases, leading to reduced\nβ-estradiol production. In turn, a lower fluorescence signal\nis observed.\nIn vitro testing of CYP19A1 inhibitor candidates. (a)\nIllustration\nof the enzymatic assay based on heterologous expression of human CYP19A1\nused for testing candidate CYP19A1 inhibitors. (b) Screening of compounds\nfor human CYP19A1 inhibition by the 10 selected compounds and the\nreference compound letrozole at inhibitor concentrations of 10 μM\n(left panel) and 1 μM (right panel). (c) IC 50  curves\nfor  8 ,  9 ,  9a , and  9b . Data are presented as mean ± SEM ( n  = 3).\nTen compounds identified through virtual screening\nwere tested\nin this aromatase inhibitor screening assay to evaluate their relative\npotency compared to letrozole. 45  Initially,\nthese inhibitors were applied at a final concentration of 10 μM.\nCompounds  1 ,  2 ,  4 ,  7 ,  8 , and  9  demonstrated significant inhibition\nof CYP19A1 activity at 10 μM and were further tested at 1 μM.\nAt this concentration, only  8  and  9  significantly\ninhibited CYP19A1 activity ( Figure  6 b). Among these two, the imidazole and coumarin ring\ncontaining  9  ( N -(1-(4-(1 H -imidazol-1-yl)phenyl)ethyl)-6-chloro-2-oxo-2 H -chromene-3-carboxamide)\nwas identified as the most effective inhibitor of CYP19A1, with an\ninhibition rate of 13 ± 2%, as compared to 17 ± 1% for letrozole\n( Table S2 ). As  9  was received\nas a racemate and exhibited potent CYP19A1 inhibition, two analogs\nof  9  ( 9a  and  9b ), lacking the\nmethyl group at the stereocenter, were identified from the Enamine\ndatabase. In addition,  9a  lacks a Cl atom and  9b  lacks a keto group in comparison to  9 . These analogs,\nalong with  8  and  9 , were also tested for\ntheir IC 50  values ( Table  5  and  Figure  6 c). Again,  9  proved to be the most potent inhibitor,\nwith an IC 50  value of 271 nM.\nSeveral factors were considered to avoid assay artifacts. 46  The biosensor gene circuit used in this study\nis orthogonal to the yeast chassis, 47  ensuring\nno interference between the biosensor and the yeast host, thus preventing\nfalse positives. Furthermore, the concentrations of β-estradiol\nused in the tests were demonstrated to be noncytotoxic, which is crucial\nin avoiding both false-positive and false-negative results. 48  To further confirm the absence of interference\nbetween the inhibitors and the biosensor, a control test was conducted,\nwhere a 10 μM inhibitor sample was incubated with the biosensor\nin the presence of 30 nM testosterone and 30 nM β-estradiol.\nThe results showed no significant difference from the control group\n(without the inhibitor), as shown in  Figure S1 .\nThe most\nactive inhibitor of CYP19A1 in this study, an imidazole and coumarin\nring-containing molecule,  9 , possesses a single stereocenter\nand was purchasable only as a racemic mixture. To further understand\nthe stability of the interactions and the binding mode in a dynamic\nsystem, both  (R)- 9  and  (S)- 9  were docked to the CYP19A1 binding site ( Figure  7 a,b), and these bound\nconformations served as starting points for MD simulations. Both isomers\nconsistently inhibited the heme iron during our MD simulations of\n100 ns in five replicas: in 99.4% of the trajectories of  (S)- 9  and 100% for  (R)- 9 .\nThe hydrophobic contacts between both isomers’ imidazole-connected\nbenzene ring and the nearby lipophilic residues (i.e., Ile133, Trp224,\nThr310 and Val310) in the binding site of CYP19A1 were observed in\nmore than 99% of the trajectories. The hydrophobic contacts between\nboth isomers’ methyl groups at the stereo centers and Leu477,\nTrp224 and Phe134 of CYP19A1 were also stable during the simulations\n(100% of trajectories for  (S)- 9  and\n99.8% of trajectories for  (R)- 9 ). We\nalso compared  9a  and  9b  with the same MD\nsimulation settings. The absence of the methyl group did not cause\nany decrease in heme inhibition, stressing the importance of the overall\nshape complementary of the protein with the molecule rather than the\ndifference of a single methyl group at the stereocenter. Interestingly,  (S) - 9  exhibited an additional interaction between\nits carbonyl oxygen and Ser478 in 19.9% of the trajectories, which\nwas absent in  (R) - 9 . We visualized the\ninteraction between  9  and CYP19A1 using our in-house\ndynamic pharmacophore (dynophore) tool, 49 , 50  which summarizes\nthe interactions during MD simulations in each trajectory. Interestingly,\nwe found two plausible binding modes for  (R)- 9  and one for  (S)- 9 , as shown\nin  Figure  7 c,d. The\nbinding modes 1 and 2 of  (R)- 9  mainly\ndiffer in the positioning of the interactions contributed by the coumarin\nmoiety. The interaction patterns of the two binding modes are similar—π–π\ninteractions between the 2-pyrones and Phe221 or His480, π–π\ninteractions between the coumarin benzene ring and His480, and hydrophobic\ncontacts between the coumarin benzene ring and Phe221. In addition,\na halogen bond interaction between the chlorine and Arg192 was observed\nduring 4% time of the MD simulation while the molecule is arranged\nin binding mode one, and hydrophobic contacts with the chlorine group\nwith Phe221 in binding mode 2. The two binding conformations majorly\ndiffer due to protein flexibility but not different interactions.\nFor  (S)- 9 , the coumarin moiety was stabilized\nthrough aromatic interactions with Phe221, and hydrophobic contacts\nwith Phe221 and Val313.\n(a, b) Docking poses of  (R) -9  and  (S) - 9  shown in\nthe CYP19A1 binding site. (c,\nd) Dynamic pharmacophore (dynophore) interactions shown as interaction\nclouds showing the spatial extent of  (R) -9  and  (S) - 9  binding to the CYP19A1 active\nsite during MD simulations. Iron binding is shown as light blue clouds,\nhydrophobic contacts are shown as yellow cloud, hydrogen bond donors\n(HBD) are shown as green clouds, hydrogen bond donors (HBA) are shown\nas red clouds, aromatic interactions are shown as dark blue clouds,\nand halogen bond interactions are shown as purple clouds.\nDuring our analysis of the binding mode, we hypothesized\nthat either\nAsp309 or Ser478 could potentially be targeted by  9  as\na covalent binder. For Asp309 to serve as a viable target, it must\nbe protonated. There are reported cases in the literature where aspartic\nacid residues have been successfully targeted by covalent binders,\nhighlighting the feasibility of such an approach in drug design. 51  In our MD simulation setup, Asp309 was automatically\ncalculated to be protonated, based on the microenvironment of the\nbinding site. The crystal structures also suggest that protonation\nof Asp309 is necessary to stabilize a water molecule involved in steroidal\ntransformation. 16 , 52  With respect to 9, its aromatic\n2-pyrone ring enables electron delocalization toward the α,β-unsaturated\ncarbonyl structure. When linked via the amide bond, this structural\nfeature suggests that the compound could function as a weak covalent\nbinder, potentially acting as a Michael acceptor. Ring-opening reactions\nof 2-pyrones have been observed in nature. 53  After the hypothetical reaction ( Figure S2 ), the resulting carboxylic acid carboxylate could form an ionic\nbond with His480, while the remaining benzene ring could maintain\nπ–π interactions. However, no significant differences\nwere detected when we tested this hypothesis by comparing the activity\nof compound  9  with and without preincubation during culturing.\nOur inhibitor occupies the binding pocket that exemestane (PDB code:  3S7S ) 16  faces to establish the covalent bond, as compared in  Figure S3 . The residues potentially targeted\nby the warhead of exemestane are Asp309 or Ser478.  9  does\nnot seem to target these residues as seen in the assay. Despite being\nclose to Asp309 or Ser478,  9  is far away from heme iron,\nmaking it unlikely to undergo the same inhibition reaction as exemestane.\n\nThe conformational flexibility and complex\ncoordination chemistry\nof CYP19A1 pose formidable challenges to structure-based methods such\nas docking. In view of these challenges, we explored the capacity\nof ML to learn the existing measured data on CYP19A1 inhibition and\nidentify novel, promising CYP19A1 blockers from libraries of purchasable\ncompounds. By training the ML models on two distinct data sets (derived\nfrom the ChEMBL database and PubChem BioAssay), our ML approach could\ncover a vast chemical space. Using BEDROC as the primary performance\nmetric, we obtained models yielding sound early enrichment, which\nis key to the success of screening campaigns.\nHarnessing the\nML models, established compound filters, and clustering,\nwe reduced a library of commercially available compounds from 4.6\nmillion to 1500 compounds. From these 1500 compounds, we selected\nten compounds for experimental validation, of which seven showed activity\nin an enzymic assay.\nML accelerates the early filtering stage\nin virtual screening campaigns,\nwhich were traditionally performed using high-throughput docking for\nother targets but are less feasible for CYPs, and docking, combined\nwith visual inspection, complements this process by incorporating\nstructure-based knowledge into the final selection for the success\nof our screening process.\nWith an IC 50  value of 271\n± 51 nM, we identified\ncompound  9  as the most potent inhibitor of CYP19A1 among\nthe ten selected compounds. Interestingly, the compound, containing\nan imidazole moiety and a coumarin ring connected as a different scaffold,\nis structurally clearly distinct from any known NSAIs. We also derived\nthe likely binding mode of  9  interacting with CYP19A1\nby utilizing structure-based modeling approaches and MD simulations.\nWe hope these insights and results will support the further development\nof much-needed, innovative inhibitors of CYP19A1.\n\nTwo separate data\nsources were utilized to investigate compound activity against CYP19A1\n( Homo sapiens ). The first source was\nthe ChEMBL database 33  (Target ID: CHEMBL1987,\naccessed March 2023). From this database, compounds with reported\nIC 50  values and corresponding pChEMBL values were selected.\nCompounds with a pChEMBL value ≥6 were categorized as active,\nwhile those with lower values were categorized as inactive. The second\nsource was PubChem 32  BioAssay 743139 (Tox21),\nwhich summarizes data from BioAssay 743083 (aromatase antagonist mode\nassay) and BioAssay 7431084 (cell viability counter assay), both targeting\nCYP19A1. Compounds were labeled as inactive if their PUBCHEM_ACTIVITY_SCORE\nwas 0, active if the score ranged between 40 and 100, and those with\nscores between 5 and 30 were excluded from the study. The data sets\nfrom the ChEMBL database and PubChem BioAssay database were analyzed\nseparately and were not combined for model training.\nAll compounds\nfrom both data sources were standardized using the ChEMBL Structure\nPipeline. 54  Additional standardization\nsteps were applied, including tautomer canonicalization, removal of\nstereoisomer information, and elimination of salt components from\nSMILES strings according to the rules described in Hit Dexter. 55  Compounds containing elements other than the\ndefined common organic compounds elements (i.e., H, B, C, N, O, F,\nSi, P, S, Cl, Se, Br, and I) or with molecular weights outside the\nrange of 250–900 Da were excluded. These preprocessing steps\nwere carried out using the “csp_wash” method from the\n“MoleculePreprocessorExtended” class in the RingSystems\nlibrary. 56  Additionally, duplicates within\neach data set were removed based on canonical SMILES. If compounds\nin the same data set had identical labels, only one was retained,\nwhile compounds with conflicting labels were removed. When checking\nthe ChEMBL and PubChem data sets for inconsistently classified compounds,\nwe identified a single compound, nordihydroguaiaretic acid, which\nwere labeled differently in the two data sets. In the ChEMBL data\nset, the activity data was curated from the summary of a virtual screening\nstudy, reporting an IC 50  value of 11 nM. However, upon\ntracing back to the original study cited by this virtual screening\nwork in 1993, the compound was reported with an activity of 11 μM. 57  To ensure consistency in data treatment across\nthe entire data set, we retained the data point as it appeared in\nthe ChEMBL data set, as it would not significantly impact the overall\nperformance of the model.\nDrugBank 36  (download date: March 14th,\n2023) was preprocessed using the same workflow as the ChEMBL and PubChem\ndata sets and subsequently used for chemical space comparison. The\nprepared data sets were analyzed via principal component analysis\n(PCA), which was based on 13 physicochemical properties of the compounds\ncalculated with RDKit (version 2021.03.2): 58  number of nitrogen atoms, number of oxygen atoms, number of chiral\ncenters, molecular weight, number of heavy atoms, number of hydrogen\nbond acceptors, number of hydrogen bond donors, log  P , topological polar surface area, number of aromatic atoms,\nsum formal charge, number of rings, and the fraction of sp3 hybridized\ncarbon atoms.\nThe molecular structures\nwere featurized using two types of molecular descriptors: RDKit 2D\ndescriptors (208 physicochemical descriptors) calculated with “MolecularDescriptorCalculator”\nfunction and ECFP4 fingerprints (2048 bits) calculated with “GetMorganFingerprintAsBitVect”\nfunction from RDKit. 58  The ChEMBL and PubChem\nBioAssay data sets were each randomly split into a training set and\nan external test set using the “train_test_split” method\nfrom scikit-learn (version 1.3.2), 59  with\na 4:1 train-test ratio, shuffling enabled, and the random_state set\nto 42.\nIn total, four Random Forest (RF) classifiers were trained\nusing the ChEMBL and PubChem BioAssay data sets, with either RDKit\n2D descriptors or ECFP4 fingerprints. The RF classifiers were developed\nusing the scikit-learn library, with the following fixed hyperparameters:\n“n_estimators” was set to 1000 to build a robust ensemble\nof decision trees, “min_samples_split” was set to 2\nto allow more flexible splits with minimal samples, and the “random_state”\nwas fixed at 42 to ensure reproducibility. In addition, “max_features”\nwas optimized during cross-validation exploring values of None, “sqrt”,\n0.2, 0.4, and 0.8 to control the number of features considered for\neach split. The optimal parameters were selected based on maximizing\nthe Boltzmann-enhanced discrimination of receiver operating characteristic\n(BEDROC). 38\nTo ensure the model’s\nability to predict novel structures,\ncompounds with a Tanimoto similarity (calculated by “BulkTanimotoSimilarity”\nfunction from RDKit) greater than 0.3 were kept separate in the training\nand testing splits during cross-validation: the training sets, as\ndescribed in the  Data Set Preparation and Analysis  section, were clustered based on the similarity using Butina clustering\nalgorithms, 60  allowing lowest Tanimoto\nscore of 0.3. The resulting clusters were shuffled and distributed\ninto five groups, with the principle of keeping the total number of\ncompounds in each group as similar as possible. Followed by clustering,\nto address class imbalance, we employed Synthetic Minority Oversampling\nTechnique (SMOTE) or Synthetic Minority Oversampling Technique for\nNominal and Continuous data (SMOTENC) from the Imbalanced-learn library 61  to the training data within each fold during\ncross-validation (CV) with the steps of: 1. Split the data into training\nand test sets for the current fold. 2. Apply SMOTE to the training\nset to balance the class distribution. 3. Train the classifier using\nthe oversampled training set. 4. Predict probabilities for the test\nset using the trained classifier. 5. Record predictions for performance\nevaluation.\nThe optimum hyperparameters were used to train the\nRF classifiers\non the whole training set, with same oversampling method applied for\nboth training models on the training data for external testing, and\nmodels on the complete data for screening. The performance was evaluated\non the external test set using metrics such as the area under the\nreceiver operating characteristic curve (AUC-ROC) and BEDROC. After\nevaluation, each test set was combined with the corresponding training\nset to develop the final model, which was then applied for virtual\nscreening.\nAll stock compounds from MolPort\n( www.molport.com ) were retrieved\n(downloaded date: March 10th 2021) and used for virtual screening.\nThe same standardization steps used for the ChEMBL and PubChem BioAssay\ndata sets were applied to the MolPort compounds. The database was\nalso preprocessed by filtering for compounds that obey the rule of\nfive (ro5) and the presence of the common substructure of CYP inhibitors\n(imidazole, triazole, or pyridine). ECFP4 and RDKit 2D descriptors\nwere calculated for all the remaining compounds in the database. These\ncalculated descriptors were then used as input for the trained RF\nclassifiers, and the predicted probabilities of the compound being\nactive were assigned from each model.\nDuring virtual screening,\nseveral filters and criteria were applied to refine the compound selection.\nCompounds were retained if their predicted probabilities exceeded\n0.7 and were ranked among the top 10,000 by any of the four models.\nTo ensure novelty, only compounds with similarity to any known active\ncompounds below 0.7 were included. Compounds matching any PAINS pattern\nwere excluded, as were compounds with specific undesirable features\nsuch as a positive charge, long alkane moieties, bromine or iodine\ngroups, nitro groups, or basic amine groups. Additionally, compounds\nwith fewer than eight rotatable bonds were retained to prioritize\nthose with favorable conformational flexibility.\nThe remaining\ncompounds were clustered with a maximum Tanimoto\nsimilarity of 0.6 based on ECFP4 fingerprints (2048 bits), and compounds\nof interest from different clusters were manually selected for further\nanalysis. The selected compounds were docked to assess potential steric\nhindrance and their ability to reach the heme iron in the binding\nsite.\nThe atomistic\nstructure of human CYP19A1 was obtained from the Protein Data Bank 62  (PDB entry:  5jkv ( 15 )). Co-crystallized\nwater molecules, androstenedione, pentaethylene glycol, and phosphate\nions were removed, while the protein structure, including the heme\nprosthetic group, was retained. The structure was prepared using MOE\nv.2020.0901 (Molecular Operating Environment; Chemical Computing Group\nULC, Montreal, QC, Canada), utilizing the integrated Structure Preparation\ntools. The structure was optimized and protonated using Protonate3D 63  with the OPLS-AA 64  force field.\nSelected small molecules were docked\ninto the binding pocket of the prepared CYP19A1 structure using GOLD\nv.5.8.1 (Genetic Optimization for Ligand Docking; CCDC Software, Cambridge,\nU.K.). 65  The search efficiency was set\nto 100%, generating 10 docking poses per molecule. The docking center\nwas defined by the coordinates of the heme iron, with a surrounding\nsphere of 10 Å radius. The scoring function used was “goldscore_p450_cds”,\nwhile all other settings remained at their default values. The docking\nresults were visualized using LigandScout v.4.4.3, 66 , 67  and the conformations were further minimized with the MMFF94 force\nfield. 68\nIn order to investigate the interaction\ndynamics of  9  and its analogs to CYP19A1, all-atom MD\nsimulations were carried out with the most plausible binding hypotheses\nas starting conformation. The prepared complex was loaded into Maestro\nv. 13.1.137 (Schrödinger Release 2022–1: Maestro, Schrödinger,\nLLC, New York, NY) for structure optimization with the implemented\nfunctionality “Protein Preparation Wizard”. The heme\niron atom type was manually corrected to Fe 3+ . The simulation\nenvironment was prepared with “System Builder”. The\nprotein–ligand structures were solvated in cubic water boxes\nwith padding of 10 Å filled with TIP3P water model. 69  Automatically calculated numbers of chloride\nor sodium ions were added to each system for neutralization, and another\n0.15 M NaCl was added to mimic the physiological environment. The\ngenerated systems were simulated with the Desmond simulation engine\nv. 6.9 70  on water-cooled NVIDIA RTX 2080\nTi graphical processing units (GPU) for 100 ns in five replicas for\neach system, using the OPLS-AA force field. The simulation temperature\nand pressure were kept at their default values of 300 K and 1.01325\nbar, respectively. The system was relaxed and equilibrated following\nthe standard seven-step protocol. System coordinates were recorded\nevery 100 ps. The coordinate and trajectory files were wrapped and\naligned with VMD v. 1.9.3. 71  For the full\nsimulation trajectory, the frequencies of inhibitor-enzyme interactions\nwere derived and analyzed using the dynamic three-dimensional (3D)\npharmacophores analysis method (“dynophores”). 49 , 50\nThe inhibitor candidates ( 1 – 10 ) were purchased from MolPort (Riga,\nLatvia). Compound  1  was obtained from BIONET—Key\nOrganics Ltd., compounds  2 – 4  were\nobtained from Life Chemicals Inc., and compounds  5 – 10  were obtained from Vitas M Chemical Limited. Compounds  9a  and  9b  were purchased from Enamine. Letrozole\nwas obtained from Enamine Ltd. (Riga, Latvia).\nCompounds  1 – 10, 9a, 9b  and letrozole were >95%\npure\nbased on HPLC analysis, as listed in the  Supporting Information . The chromatographic analysis was performed using\nan Agilent 1290 Infinity HPLC system, consisting of a binary pump,\nautosampler, and column compartment, coupled with an Agilent 1260\nDAD VL+ Detector and an Agilent 6130B Single Quadrupole MS. Separation\nwas achieved on an Agilent Poroshell C18 column (100 × 2.1 mm 2 , 2.7 μm particle size) maintained at 30 °C. The\nmobile phase comprised Solvent A (Water with 0.1% Formic Acid) and\nSolvent B (Acetonitrile with 0.1% Formic Acid), delivered at a flow\nrate of 0.400 mL/min with the following gradient: 5% B from 0.0 to\n1.0 min, linearly increased to 95% B from 1.0 to 8.0 min, held at\n95% B until 10.0 min, then returned to 5% B by 10.5 min, with a total\nruntime of 11.00 min. The injection volume was 0.5 μL. The DAD\nmonitored signals at 254, 210, and 220 nm (with the latter used for\nanalysis), at a scan rate of 40 Hz. The MS operated in both positive\nand negative scan modes, with an  m / z  range of 50–700.\nSynthetic defined complete (SDC) liquid\nmedium was used as the\nculture medium of  Saccharomyces cerevisiae . The medium was composed of 79.2 mg histidine (H108260 Aladdin Beijing,\nCN), 396 mg leucine (L104898 Aladdin Beijing, CN), 79.2 mg tryptophan\n(T103480 Aladdin Beijing, CN), 79.2 mg uracil (U102087 Aladdin Beijing,\nCN), 2 g amino acid mix (0.5 g adenine, and 2.0 of each of the following\namino acids: adenine, alanine, arginine, asparagine, aspartic acid,\ncysteine, glutamine, glutamic acid, glycine, isoleucine, lysine, methionine,\nproline, serine, threonine, tyrosine, valine, phenylalanine—Aladdin\nBeijing, CN), 6.7 g Yeast Nitrogen Base ( Q30009  Thermo Fisher Beijing,\nCN) and 20 g glucose (G116307 Aladdin Beijing, CN). These components\nwere dissolved with 1 L double distilled water and autoclaved to prepare\nthe SDC liquid medium. Steroids β-estradiol (E2758) and testosterone\n(T5411) were purchased from Merck (Beijing, CN).\nThe biosensor strains, initially\nobtained from synthetic selective medium plates (diluted to an OD 600  of 0.04), were precultured in 3 mL of SDC liquid medium\nat 30 °C and 240 rpm for 14 h. Then, 20 μL of the precultured\nbiosensor cell solution (diluted to an OD 600  of 5.00) was\nadded to 1980 μL of SDC medium containing β-estradiol,\ntestosterone, and an inhibitor. This 2 mL solution, containing the\ninducer and biosensor cells, was incubated at 30 °C and 240 rpm\nfor 24 h. After incubation, 5 μL of this cell suspension was\ndiluted into 195 μL of double-distilled water (1:40 dilution)\nbefore analysis using a BD FACSVerse flow cytometer (laser 488 nm,\nFITC filter 527/32 nm). A total of 10,000 cell events were collected\nper measurement. Performance Quality Control (PQC) was conducted monthly,\nwith fluorescent beads (BD FACSuite CS&T Research Beads 650621)\nused to adjust the FITC voltage, ensuring reliable and reproducible\nFACS results. The relative difference between bead peaks was required\nto be less than 5% when comparing each measurement to the initial\nbead measurement. The  flowCore  R-Bioconductor package\nwas used for the analysis of the FACS data, with each mean value calculated\nfrom flow cytometry measurements of three independent experiments.\nThe  S. cerevisiae  strain byMM1712, which expresses both\nhuman cytochrome P450 reductase (CPR) and human CYP19A1, was incubated\nin 10 mL of Synthetic Defined Complete (SDC) liquid medium at 30 °C\nand 240 rpm for 14 h. A 14-h precultivation period ensured that the\naromatase catalytic capacity remained consistent across all independent\nexperiments and the inhibitor cannot act on P450 expression levels.\nThe initial OD 600  was diluted to 0.04, and then the yeast\nculture’s OD 600  was measured and adjusted to 2.00\nwith fresh SDC medium. One mL of yeast solution (OD 600  =\n2.00) was poured into the wells of a 48-well deep-well cell culture\nplate (Axygen P-5 ML-48-C). Subsequently, 1 mL of SDC medium containing\n120 nM testosterone and either 10 μM or 1 μM of the inhibitor\nwas added to each well containing the yeast solution. The mixture\nwas incubated at 30 °C and 240 rpm for 2 h. After incubation,\n1.5 mL of the mixture from each well was transferred to 2 mL tubes\nand centrifuged at 16,000 rpm to pellet the yeast cells. Then, 1.2\nmL of the supernatant was moved to 1.5 mL tubes and labeled according\nto the inhibitor concentration. For each sample, 1 mL of the supernatant\nwas transferred to a different well of a 48-well deep-well plate.\nNext, 980 μL of SDC medium and 20 μL of the precultured\nbyMM1984 biosensor cell solution (OD 600  adjusted to 5.00\nafter 14 h of preculture) were added to each well. Finally, flow cytometry\nwas performed as previously described. The entire assay was repeated\nin three independent experiments to calculate the mean values.\nAll repeated experimental results\nare presented as mean ± SD or mean ± SEM. Statistical significance\nwas determined using a two-tailed  t  test. Differences\nwere considered significant if  P  < 0.05. IC 50  curve was obtained by fitting the raw data to the empirical\nHill function, where y is fluorescence (A.U.), and x is proportional\nto the inhibitor concentration (log 10 , nM). Statistical\nanalysis was done with GraphPad Prism 8.01 (GraphPad Software Inc.,\nLa Jolla, CA). All repeated experimental results are presented as\nmean ± SD or mean ± SEM. Statistical significance was determined\nusing a two-tailed  t  test, with differences considered\nsignificant if  P  < 0.05. The IC 50  curve\nwas generated by fitting the raw data to the empirical Hill equation,\nwhere  y  represents fluorescence (A.U.), and  x  is the log10-transformed inhibitor concentration (nM).\nStatistical analyses were performed using GraphPad\nPrism 8.01 (GraphPad Software Inc., La Jolla, CA).","source_license":"CC-BY-4.0","license_restricted":false}