Pharmacophoric Investigation of a Natural Product-like Class of Aromatase Inhibitors Using Molecular Modeling

A structural similarity search was conducted using the SciFinder database to build an initial pool of ligands that retain at least 70% structural similarity to 2a _ 1 . The resulting pool of molecules was filtered to remove any that did not retain the scaffold structure of 2a _ 1 such that the position, number, and identity of functional groups was allowed to vary only on the B ring. Molecules containing heavy halogen atoms (Br and I), arsenic, silicon, or phosphate functionalities were also removed from the data set due to either toxicity or lack of sufficiently accurate parameters for simulation. This resulted in a pool of 77 potential ligands. Three-dimensional (3D) geometries of each ligand were built with Avogadro and the built-in geometry optimization function was used to obtain a minimized structure. The ligands were docked to the aromatase active site (PDB code 3s79 ) using AutoDock Vina and the included prepare_receptor4.py and prepare_ligand4.py programs to assign charges and atom types to the protein and ligand. , The docking grid center coordinates used were (85.919, 50.603, 47.323), the box size was (15.654 17.782, 18.848 Å), and docking calculations were conducted with an exhaustiveness of 64 and with the number of generated models capped at 20. The docking scores ranged from −11.1 to −7.9 kcal/mol, and we selected any ligand with a predicted binding affinity less than or equal to the median (−9.5 kcal/mol) for further study with molecular dynamics (MD). This corresponds to 42 ligands (55%) in addition to the scaffold control ( Figure B). We also conducted control simulations of aromatase bound to the natural substrate ASD ( Figure A), and the three FDA-approved aromatase inhibitors (Exemestane, Letrozole, and Anastrozole ). To validate the docking protocol, the natural substrate androstenedione (ASD) was redocked to the receptor and compared to the crystal structure. The two ASD structures have an RMSD of 0.436 Å, as shown in Figure S1 . The AMBER package was used to build and simulate the protein–ligand complexes. The systems were built in tleap using the ff14SB protein parameters and the 3s79 crystal structure, TIP3P water and ion parameters and GAFF for the ligands. , Parameters for deoxygenated heme bound to Cys437 were obtained from Shahrokh et al. The bond from the heme iron to the Cys437 sulfur and the bonds for the Cys437 backbone to the adjacent residues were explicitly called in tleap . Previous work has suggested the importance of the protonation state of Asp309 for substrate entry and exit, postulating that deprotonated Asp309 corresponds to an open entrance cavity to facilitate ligand entry prior to catalysis. Because MD simulations model precatalytic activity, we chose to use deprotonated Asp309 for all simulations. The protein–ligand complex was immersed in a truncated octahedral box of TIP3P water with a buffer of 11 Å, and the system was neutralized with three chloride ions. Energy minimization of the entire system was achieved through seven individual minimizations of 1000 steps of steepest descent and 4000 steps of conjugate gradient with positional harmonic restraints on all non-hydrogen atoms starting at 10.0 kcal/mol*Å 2 , with these gradually relaxed over the seven minimization steps until they have been fully removed. Following energy minimization, the system was heated using the Langevin thermostat and the SHAKE algorithm with a 2 fs time step, with the constant pressure periodic boundary condition and all solute atoms restrained with a force constant of 10.0 kcal/mol*Å 2 . Heating took place first from 10 to 100 K and then from 100 to 310 K over 2 ns. Finally, the solvent was relaxed in a manner similar to minimization over seven individual steps, each 500 ps in length for a total of 3.5 ns and starting from a position restraint force constant of 10.0 kcal/mol*Å 2 that was slowly relaxed until fully removed. Following this system preparation, unrestrained MD was run for 1.0 μs under the NPT ensemble. Each MD simulation used the SHAKE algorithm with a 2 fs time step and a nonbonded energy cutoff of 9 Å. Simulations were visualized using VMD. This protocol was developed based on a standard protocol for explicit solvent systems using AMBER, consistent with the published method in Salomon-Ferrer et al. The procedure is also akin to the tutorials provided online ( https://ambermd.org/tutorials ). Trajectories were processed, and root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), protein–protein and protein–ligand hydrogen bonding, and ligand radius of gyration are calculated with cpptraj in AMBER. CPPTRAJ is freely available in AmberTools for the analysis of MD trajectories. The program is described by Roe and Cheatham, and introductory tutorial is available online ( https://ambermd.org/tutorials ). The MM-GBSA method with MMPBSA.py was used to estimate the average binding free energy over time and to obtain decomposed per-residue free energy. The average binding free energies over the trajectory were used in comparison to those of the scaffold throughout this work. A QSAR model was trained to predict the experimental IC 50 values of the isoflavanone-adjacent compounds. The model was trained on a data set of IC 50 values of competitive aromatase inhibitors obtained from the ChEMBL database and descriptors calculated with rdkit. The data set was curated to remove duplicates and to only include IC 50 values that were obtained under the same experimental conditions, resulting in 1325 ligands total. The model was trained on 85% of the data set (1126 ligands), with the remaining 15% (199 ligands) being excluded from the training process to use for external validation. Of the descriptors calculated by rdkit for the training set, any descriptors that contained NaN for more than 20% of the compounds, descriptors with low variance (0.95) were excluded. Feature selection for training the model from the remaining descriptors was conducted using the LassoCV and RandomForestRegressor modules in sklearn. The union of the top 60 features from each method was selected for training the final model with LightGBM. The final model was used to predict the IC 50 values of the external validation set, and these predicted values were compared to the true values, resulting in a correlation coefficient of 0.72 ( Figure S5 ). The model was then used to predict the IC 50 values of all isoflavanone-adjacent compounds in this study ( Table S1 ). We have prepared a GitHub repository containing all of the docking, system preparation, MD, analysis, and QSAR scripts. We have also included in this repository the P-values for several statistical analyses (see Results ). The repository can be found at https://github.com/aheld7/Supplemental-Materials-Held-2025-Aromatase-Inhibitors-.git .

Molecular dynamics simulations were conducted (see Methods ) of the natural substrate androstenedione (ASD), as well as the three FDA-approved aromatase inhibitors (Exemestane, Letrozole, and Anastrozole) bound to the active site of aromatase. These are compared against the simulation of the scaffold of the structure proposed by Mirzaie et al. Figure shows the binding mode, binding free energy over time, and decomposed per-residue free energy of several residues in the active site for ASD (2A), Letrozole (2B), and scaffold (2C). The binding mode image shows only the top 10 residues from the per-residue decomposition analysis, and each residue is consistently color-coded for easy comparison. Figure S2 shows the same for exemestane and anastrozole. Binding mode, free energy over time with bold running average, and color-coded decomposed per-residue free energy of selected active site residues (A) androstenedione, (B) Letrozole, and (C) isoflavanone-adjacent scaffold. The natural substrate, ASD, relies primarily on nonpolar hydrophobic interactions with Val373, Val370, Trp224, Ile133, Ile305, Ala306, and Phe134. It also forms transient hydrogen bonds with Met374 and Arg115. The binding free energy over time is relatively constant, indicating a stable binding mode with a predicted binding free energy of −42.44 ± 0.03 kcal/mol. Exemestane and Anastrozole maintain similar binding affinities to ASD ( Figure S2 ); however, Letrozole features a diminished binding affinity of −35.80 ± 0.03 kcal/mol. It features many similar interactions to ASD, namely, hydrophobic contacts with Val373, Val370, Trp224, Ile133, and Phe134, while also maintaining a diminished but still strong hydrogen bond with Met374 and hydrophilic contact with Arg115. It introduces an additional strong hydrophilic contact with Thr310 and a hydrophobic contact with Phe221. The binding mode for Letrozole is stable, similar to ASD. The scaffold (2C) mimics several interactions that are formed by ASD, namely, hydrophobic contacts with Val373, Val370, Trp224, Ile133, Ile305, Ala306, and Phe134, and a hydrophilic contact with Arg115, and it also mimics some characteristics of Letrozole with a contact at Thr310, although it lacks any hydrogen bonding or strong interaction with Met374. Although it mimics many contacts and, in some cases, even has contacts stronger than ASD or Letrozole, in general, the strength of the enzyme-ligand contacts is diminished, resulting in a low binding affinity of −28.55 ± 0.03 kcal/mol. Furthermore, the binding affinity over time is unstable, indicating difficulty in finding an energy minimized binding mode at. This is likely due to the small size of the scaffold in comparison to ASD and Letrozole in the binding pocket, allowing the scaffold to experience greater freedom of motion and, therefore, diminished enzyme-ligand interactions. Indeed, the root-mean-square fluctuation (RMSF) of ASD is 0.78 Å, while that of the scaffold is significantly larger at 2.42 Å. The scaffold successfully mimics some enzyme-ligand interactions of both ASD and letrozole; however, given the increased freedom of motion and decreased binding affinity, it is clear that additional functionality is required to fit the ligand to the pocket and strengthen the existing interactions. In the next section, we examine the effects of various functionalities on the B ring of the scaffold ( Figure B). Following the structural similarity search of 2a _ 1 ( Figure B), 77 structures were obtained and docked to the aromatase active site using AutoDock Vina. The ligands with a predicted binding free energy in the top 55% were subjected to a 1.0 μs MD simulation in AMBER, and the binding free energies were calculated with MMPBSA.py (see Methods ). Table S1 contains the SMILES, docking score, MM-GBSA score (if applicable), and ligand RMSF from MD (if applicable) of every ligand studied. Below, we present the effects of functional group identities and positional occupancies on the binding free energy relative to that of the scaffold. Figure A shows the distribution of all functional groups for every ligand in the data set according to identity and positions that they occupy without considering total positional occupancy, while Figure B shows the same thing for only those ligands that qualified for full MD simulation. Figure C shows the frequency of highly ranked (fully simulated) ligands and lowly ranked (docked only) ligands based on position, size, and polarity. Bulky functional groups are defined as any group containing more than one heavy atom, and polar functional groups are defined as exhibiting a dipole moment of any strength. These data are summarized in tabular format in Table S2 . Distribution of all functional groups in the full data set (A); distribution of all functional groups that were ranked highly enough by AutoDock Vina to be fully simulated (B); and frequency of ligands that were ranked highly enough to be fully simulated and frequency of ligands that ranked too low to be fully simulated based on size, polarity, and position, with the total number of ligands given in the table (C). The most explored functionalities in the full data set are fluoro, methyl, chloro, and methoxy groups. This is reflected in the ligands that qualified for full simulation, being mostly represented by fluoro, methyl, and chloro functionalities. The methoxy functionality, however, is less prevalent, being represented rarely alongside ethoxy and amine, which were the only other polar functionalities to qualify for study with MD. Other nonpolar functionalities were even less represented and consist only of ethyl and cyclobutyl functionalities. Regarding positional occupancies, position 1 in the full data set consists almost exclusively of polar groups with the exception being only methyl, and most of the highly ranked ligands exhibited a fluorine at this position. Of all the ligands that have functionality on position 1 within the full data set, very few are bulky. This rarity in the data set of existing molecules indicates that there may be intramolecular steric clashes that preclude bulky groups from occupying this position. The largest frequency of functionalities at this position is polar, and all but one did not qualify for MD analysis. While the data set does generally lack nonpolar groups at position 1, this vast majority of highly ranked ligands with a polar group here may indicate a preference for polar functionality at position 1. In the highly ranked data set, position 2 is similar to position 1 in favoring a larger total percentage of polar groups; however, it does favor nonpolar groups more heavily than position 1 in the full data set. Position 2 also exhibits very few bulky functionalities in the full data set with more being polar than nonpolar. However, while the singular nonpolar functionality qualified for further analysis, none of the polar bulky groups did, which indicates a preference for nonpolar groups at this position if the group is bulky. Regarding small functionalities, there are many more polar groups than nonpolar groups that qualified for further analysis; however, ligands with a small polar or nonpolar group both had a higher percentage of ligands qualify for further analysis than did not. This indicates that there is no clear preference for polarity of small groups at this position, and the larger number of polar groups is likely an artifact of more ligands like this being in the full data set. Regarding position 3, the full data set shows the most diversity of functionalities at this position. The data set consists of many bulky polar groups and a smaller but still significant number of bulky nonpolar groups; however, none of them ranked highly enough to qualify for MD analysis. This indicates that the enzyme pocket is too small to accommodate large groups in an energetically favorable manner when a bulky group is placed at this position. It is difficult to determine a specific steric clash with the active site that is responsible for this low predicted binding affinity, as these ligands lack a consistent binding mode. For smaller groups, however, many ligands with occupancy at this position ranked highly. In the 42 ligands with any functional group on position 3, 28 (67%) were polar, and 14 (33%) were nonpolar. Therefore, similar to position 2, it is unclear if there is a preference for polar or nonpolar groups because any seeming preference could be an artifact of mismatching frequencies in the data set. While fluorine is very prevalent in all other positions, position 4 displays it rarely ( Figure A), instead favoring methyl groups and other nonpolar functionalities as well as chlorine. The highly ranked ligands display mostly a methyl at this position; however, a large percentage of polar functionalities are also found. The majority of ligands exhibiting a small nonpolar group at position 4 ranked highly, and the same is true for every small polar group, again making it difficult to predict a preference for small groups at this position. Regarding bulky groups, many of the available bulky polar groups are highly ranked; however, none of the bulky nonpolar groups are, which may indicate a preference for polar groups should the functionality be bulky. However, the data set lacks a large sample size of bulky groups at position 4, indicating possible intramolecular clashes similar to position 1. Lastly, position 5 seems to equally favor fluoro, methyl, chloro, and methoxy in the full data set, but the highly ranked ligands primarily feature fluoro groups. A similar number of ligands with a bulky polar group on position 5 were ranked highly and ranked low, but all ligands with a bulky nonpolar group were ranked highly, indicating a preference for nonpolar functionalities for larger groups. Regarding small groups, similar numbers of nonpolar groups were ranked highly and ranked lowly, but many more polar groups were ranked highly, indicating a preference for small polar groups. While molecular docking is generally a sufficient method for discriminating ligands with low binding affinity from ligands with high binding affinity, it often lacks accurate binding free energy predictions and struggles to correctly rank the highly performing ligands within themselves. Therefore, to determine more clearly the effects of various functional groups on binding affinity, we will focus on the ligands that underwent a full 1.0 μs MD simulation and use the free energy of binding predicted by the MM-GBSA method in AMBER with MMPBSA.py to view the binding energetics. See Table S1 for a full list of simulated ligands. The simulation of the scaffold structure, as discussed above in Section , bound to aromatase with a binding affinity of −28.55 ± 0.03 kcal/mol ( Figure C). In the following data, ligand binding affinities are taken as a difference from this value (ΔΔ G ) to describe the gain (or loss) of binding affinity relative to the scaffold. Figure A describes the effects of functional group identity on the binding free energy. Note that a ligand is included in the calculation if the functional group in question is present at all in the structure, and did not account for multiple of the same or different functional groups also being present. The wide standard deviations make it difficult to distinguish the effect of the presence of a single functional group, which indicates that the presence of an individual group does not heavily influence the binding free energy. More likely, combinatorial effects of the functional groups determine a ligand’s success. This is supported by two-tailed t -tests, suggesting there is no significant difference between the means of any two combinations of functional groups in Figure A (data not shown). Figure B describes the effect of the total number of functional groups a ligand exhibits on the binding affinity. There are again wide standard deviations here, but a trend of an increase in the total number of groups with an increase in binding affinity is more clearly observed. This seems intuitive, as a greater number of functional groups is likely to increase the frequency of enzyme-ligand interactions and stabilize the ligand in the binding pocket, but it is also clear that in some cases a larger number of functional groups is unable to make up for clashing interactions. Two-tailed t -tests show a statistically significant difference between the mean binding affinity of compounds with 1 FG and 3 FG ( p = 0.015) but not between 2 FG and 3 FG ( p = 0.37), which supports this interpretation. Furthermore, molecules with more than 3 functional groups are rare, only occurring once in the entire data set, and there are zero instances of all 5 positions being occupied ( Figure S3 ). Taken together, these results indicate that the more functional groups does not necessarily equate to greater binding affinity, and the rarity of highly occupied ligands within a data set of existing molecules makes pursuing highly occupied ligands unattractive for eventual experimental studies. Effect of functional group identity (A) and total number of functional groups (B) on binding affinity. It is clear that the effects of polarity, size, and specific positional occupancy are necessary to gain a deeper understanding of ligand binding. Figure shows these effects more clearly. Figure A distinguishes between polar and nonpolar groups at each position, where it is clearly observed that nonpolar functionalities tend to increase the binding affinity. The high nonpolar character of the buried active site of aromatase and the nonpolar character of the substrate ASD also makes this result intuitive. It should be noted that the standard deviations are wide due to the overall binding affinity being dependent on the combinatorial effects of the functional groups bound at different positions. In fact, the only significant difference here is between nonpolar groups on position 2 and polar groups on position 1 ( p = 0.042) and polar groups on position 3 ( p = 0.048). Effect of functional polarity (A), size (B), and polarity and size (C) at each position on the binding free energy relative to the scaffold. Figure B distinguishes the sizes of the functional groups at each position. There is no clear effect of small functional groups at each position, however, it is clear that certain positions cannot accommodate bulky groups as well as others. Bulky groups on position 3 have been discussed previously in Section and result in low-ranking ligands. Bulky groups on positions 1 and 5 seem to have similar effects to their small counterparts, but position 4 shows at best only a mild improvement in binding affinity, while the one ligand with position 2 occupied with a bulky group displays a dramatic improvement. Still, the low sample size of bulky groups in the data set makes it difficult to draw definitive conclusions about bulky groups. Figure C stratifies the data between size and polarity for each position, similar to Figure C but this time against the relative free energy change against the scaffold. Position 1 seems to clearly favor small, nonpolar groups, should it be occupied. Position 2 displays one ligand with a bulky nonpolar group that is very favorable, but small nonpolar groups also seem to provide a reasonable binding affinity at this position relative to polar groups. Position 3 does not display a clear preference for polar or nonpolar groups as the distributions are wide. Position 4 seems to experience an energy penalty when occupied by a bulky polar group, but the effects of small polar and nonpolar groups seem similar. Lastly, position 5 seems to not clearly discriminate based on polarity and/or size. While these plots provide some useful information, there remains wide variation in several energies as well as low sample sizes for certain functional groups, making it difficult to make definitive assertions about which functionalities may be preferable on different positions. These wide standard deviations likely arise from the combinatorial effects of multiple functional groups. It is likely that the ligand binding profile of the scaffold would benefit from a more thorough analysis than can be provided by the general trends presented thus far. In the next section, we examine heat maps that provide a clearer view of these effects. Figure A describes the energy differences relative to the scaffold as a result of different total numbers and identities of functional groups. The number of ligands in each data point is shown, and the standard deviation (if applicable) is shown in parentheses. Effect of functional group identity and the total number of groups on the binding affinity relative to the scaffold (A), and the effect of functional group positional occupancy and polarity (shown respective to the position) (B). The number of ligands in the data point and the standard deviation (if n > 1, given in parentheses) is shown. A darker blue shade indicates a more negative ΔΔ G to −9 kcal/mol. Orange indicates a positive ΔΔ G up to +2 kcal/mol. As the heat map moves to the right, the results of Figure B are presented differently, showing that the average energy difference increases in magnitude as the number of functional groups increases. When the total number of groups is 1, nonpolar groups appear to be favored but the energies are comparatively worse than the ligands with more groups. A notable exception is fluorine, which performs similarly regardless of the total number of groups, indicating that additional highly polar and electron-withdrawing groups do not allow binding affinity to recover. When the total number of groups is 2 or more, a clear favoring of polar or nonpolar groups becomes less obvious, though most high-performing ligands seem to possess at least one nonpolar group while low-performing ligands do not, with the exception being chlorine-containing ligands. Furthermore, ligands with multiple of the same functional groups seem to generally perform better than ligands with diverse functional groups. The entire data set only contained one ligand with four functional groups, but it performed well in the docking and simulation phase. It is possible that novel ligands with four functional groups will also perform well. While Figure A shows more clearly the effect of the total number of groups against their identity, Figure B adds positional occupancy as a parameter against functional group polarity (shown in the x -axis respective to the position on the y -axis). Where positional occupancy features exclusively polar or exclusively nonpolar data points to compare, the nonpolar ligands bind more favorably in every case. Notably, occupancy 2,5 displays the highest number of ligands and the widest range of binding affinities. This higher volume may be partially due to a higher volume of these kinds of ligands in the full data set, however there are several other positional occupancies that have similar frequencies in the data set, and of these the 2,5 ligands saw the greatest percentage be ranked highly in the docking calculation ( Figure S3 ), indicating that this observation is not a result of data set artifacts. As for the wide standard deviation, this is a result of two ligands that each contain at least one fluorine (which has previously been shown to generally reduce binding affinity) that perform similarly poorly against two ligands that are globally nonpolar that score similarly highly (this will be discussed further below). A positional occupancy of 1,4 also corresponds to favorable ligands, indicating that para geometry in general is favorable, though there does appear to be a preference for 2,5 occupancy. Another observation is that ligands containing polar groups tend to perform better when these polar groups occupy higher-numbered positions than lower-numbered positions. As will be discussed in more detail in Section , this is likely not due to a specific interaction made by the polar groups, but rather they allow a more favorable geometry for the benzofuran. Finally, to illustrate the binding affinity of each ligand more clearly, in Figure we show three heat maps, one each for ligands with 1, 2, and 3 total functional groups (the one ligand with four functional groups is included in the last heat map). This allows every ligand to be presented as its own data point, which shows more clearly the effect of the total number, positional occupancy, and identity of functional groups. Heat maps for 1, 2, 3, and 4 total functional groups are shown that relate positional occupancy with functional group identity, allowing for each ligand to be viewed individually. The energy difference relative to the scaffold is shown in each cell, when applicable. A darker blue indicates a more negative ΔΔ G to −9 kcal/mol. Orange indicates a positive ΔΔ G up to +2 kcal/mol. When the total number of functional groups is 1, the binding free energies tend to be overall less favorable compared to those of ligands with more functional groups, consistent with previous findings. The best performing ligand overall for this class used methoxy on position 1, which has polar and nonpolar character, indicating that the success of this functional group at this position may be due to the diversity of interactions it can participate in. The worst-performing ligand in this class uses an amine at position 3. Overall, the best performing ligands tend to use bulky groups that are either completely nonpolar (cyclobutyl) or have high nonpolar character (methoxy and ethoxy), indicating that the reason ligands with one total functional group tend to have low binding affinities is because they lack functionality that can interact with the enzyme and/or stabilize it in the pocket, with the larger functional groups somewhat making up for this deficiency. The largest number of ligands in the data set of simulated ligands contain a total of 2 functional groups, but overall, the binding affinities seem to be comparable to those of 3 total functional groups. As previously noted, the best performing ligands all have groups arranged with para geometry, the most marked example being with chlorine atoms at positions 2 and 5. Chlorine atoms introduce local polarity; however, the para geometry induces a more global nonpolar character. This could indicate that the ligands benefit from a balance between the polar and nonpolar functionality. This is further supported by the high binding affinity exhibited by chlorine and methoxy groups arranged para on positions 2 and 5, but when the methoxy is replaced by a fluorine, the binding affinity decreases markedly regardless of the order of the functional groups. High performance is also observed with exclusively nonpolar groups arranged with para geometry, such as the bulky ethyl groups on 2 and 5 and the small methyl group on 1 and 4. It would appear, however, that the most beneficial positions to be occupied are 2 and 5 because they display the least discrimination between polar and nonpolar groups in high-performing ligands, allowing for greater structural diversity for modifications and possible substitutions to improve bioavailability, safety, etc. Furthermore, the lack of discrimination between types of interactions that the functional groups are capable of indicates that the primary role of the functional groups is to provide the positional stability of the ligand in the pocket to improve interactions near the benzofuran moiety. This is similar to the natural substrate ASD, in which most of the highly favorable interactions are made near the reactive end of the molecule ( Figure A). The worst performing ligand uses a chlorine and a methyl group on positions 2 and 3 respectively; however, in general, the worst performing ligands tend to contain fluorine atoms. This is likely because fluorine atoms are highly polar and the aromatase active site is highly nonpolar. While obviously some polar character is required for a tight binding affinity, the data show that more weakly polar groups, such as chlorine or methoxy, are better suited than fluorine. Lastly, when the total number of functional groups is 3 or 4, there are zero instances of a ligand featuring a bulky group, likely because these would be too strained and were eliminated in the docking phase but also because the data set possessed very few ligands with 3 or more total groups in which at least one was bulky. Here, the trend of ligands arranged with para geometry with at least one nonpolar group being associated with increased binding affinity is maintained, with the highest performing ligand featuring methyl groups at positions 2, 3, and 5. Occupation of these same positions with polar groups results in markedly decreased binding affinity. Another high-performing ligand features chlorine and methyl groups on positions 1 and 4, respectively, with fluorine on position 3. The lowest performing ligands lack para geometry and tend to contain more polar functionalities. In summary, based on these heat maps, there are some characteristics that seem to describe successful ligands. Namely: (1) possessing more than one functional group, (2) displaying para geometry (either 1,4 or 2,5), and (3) containing at least one nonpolar group or globally nonpolar character. When the geometry is 2 and 5, the identity of the functional groups seems not to be a large determining factor in the binding affinity as long as at least one group is nonpolar, which indicates that the role of the functional groups is largely to provide positional stability of the benzofuran. In this section, a detailed analysis of the top 5 performing ligands in the MD analysis is presented. Figure displays the same data as those in Figure , this time for the isoflavanone-adjacent ligands. Each of the best-performing ligands features binding affinities that are similar to or slightly improved from Letrozole ( Figure B). They each contain more than one functional group displaying para geometry, as discussed in the previous section, primarily in the 2,5 positions. Furthermore, two out of the top five ligands contain three total functional groups, even with the small volume of such ligands in the data set of fully simulated ligands. The same is true for ligands containing at least one bulky group. Each of these ligands primarily maintains the same top contributing contacts as ASD and Letrozole. Similar to the scaffold itself ( Figure C), these ligands feature a stronger contact with heme than either ASD or Letrozole ( Figure A,C). Figure shows a heat map of the decomposed per-residue free energy of the active site residues for ASD, letrozole, the scaffold, and each of the five ligands alongside a depiction of the active site, showing the locations of these residues (no ligand is shown). Binding mode, free energy over time with bold running average, and color-coded decomposed per-residue free energy of the top 10 most highly contributing residues to favorable binding of (A) compound 1 , (B) compound 2 , (C) compound 3 , (D) compound 4 , and (E) compound 5 . Decomposed per-residue binding free energies for ASD, letrozole, the isoflavanone-adjacent scaffold, and compounds 1 – 5 . Heme is included, but the color scale is adjusted only to −2.5 kcal/mol (dark blue) to 0.5 kcal/mol (light orange). Figure A shows the binding mode with color-coded residues, binding free energy, and energy decomposition per residue (color-coded to match the binding mode image and the results of other ligands and the controls in Figure ) of compound 1 . It is the only high-performing ligand to display para geometry on positions 1 and 4 rather than 2 and 5. It features strong hydrophobic interactions with Ile133, Phe134, Val370, and Val373, as well as a strong hydrogen bond with Met374. The other interactions it displays are generally medium to weak in comparison to other ligands ( Figure ). The binding free energy is somewhat erratic but still fairly constant in the running average, with an overall average binding free energy of −35.41 ± 0.03 kcal/mol (standard error of the mean), though it slightly increases with time. It is fairly stable in the binding pocket with a ligand RMSF of 1.02 Å; however, the radius of gyration changes, consistent with spikes in RMSD, which indicates changes in the binding mode. The full decomposed free energy, RMSF, RMSD, and ligand radius of gyration for this and the other ligands can be found in Figure S4 . Similar to ASD, where most of the strong interactions arise near the reactive end of the molecule, compound 1 also features mostly favorable interactions with residues near the benzofuran, with the stronger interaction with the heme compensating for the weaker interactions with the other residues. The residues near the B ring (Ile305-Thr310) display the biggest discrepancy compared to ASD. Namely, the strength of the interaction with Ala306 is highly diminished, and compound 1 introduces an unfavorable contact with deprotonated Asp309 via the fluorine group ( Figure ). Protonation of this residue is known to be mechanistically relevant to this enzyme, so it is possible that protonation will introduce a hydrogen bonding opportunity to improve this binding affinity. The effect of the protonation state of Asp309 would be best studied using constant pH MD simulations in a future study. Aside from this clash, there appears not to be any significant interactions formed by the substituents; instead, they serve to stabilize the ligand in the binding site to strengthen the existing interactions in the scaffold (namely, Ile133, Phe134, Trp224, Val373, Met374, and Leu477). The isoflavanone-adjacent ligands also profit from an enhanced π stacking interaction with heme in comparison to ASD and Letrozole via the B ring. Still, the strength of some contacts have been diminished in comparison to the scaffold and ASD, namely, with Arg115, but mostly with the residues near the B ring such as Ile305, Ala306, Ala307, and Thr310, and the clash with Asp309. This indicates that while the functional group geometry does stabilize the ligand to strengthen the existing interactions in the scaffold, the polar functional groups do not contribute much to the binding affinity. Figure B shows the results of compound 2 , which features polar functional groups on positions 2 and 5, though position 5 is occupied by a methoxy group that provides both polar and nonpolar character. In comparison to ASD and Letrozole ( Figure ), it displays generally stronger interactions with Ile133, Phe134, Ala307, Val373, and Met374 and generally weaker interactions with other residues (Arg115, Trp224, Ile305, Ala306, and Leu477). It has a stronger interaction with Met374 than Letrozole but weaker than ASD, and a stronger interaction with Thr310 than ASD but weaker than Letrozole. As with all previous ligands, the primary interactions are with the residues near the benzofuran moiety, and an increased interaction with the heme compared to ASD and Letrozole that compensates for the diminished interactions with other residues. Furthermore, compared with compound 1 , it features improved behavior with the functional groups of the B ring. It lacks any clash with Asp309, in fact, forming a weakly favorable interaction with it (via the chlorine on position 2) that is not present in ASD ( Figure ). Furthermore, while ASD features one strong interaction with Ala306, compound 2 instead features two moderate interactions with Ala306 and Ala307 ( Figure ). Still, compound 2 is only a mild improvement over compound 1 with an average binding affinity of −35.91 ± 0.03 kcal/mol. Furthermore, the ligand RMSF is 1.77 Å ( Table S1 ), indicating an unstable binding mode, which can also be seen by the ligand radius of gyration ( Figure S4 ). Similar to compound 1 , this indicates that while the functional group geometry improves the existing interactions, the polar functional groups still do not improve the binding affinity even without the clash at Thr310. Compound 3 ( Figure C) contains bulky nonpolar functionality on positions 2 and 5 with an average binding affinity of −36.16 ± 0.03 kcal/mol, which is trending down. In comparison to ASD and Letrozole, it features strong interactions with Arg115, Ile133, Phe134, Val370, Met374, and Leu477 near the benzofuran moiety with somewhat diminished interactions at other residues. For residues near the B ring, an improved interaction with Phe221 in comparison to all other isoflavanone-adjacent ligands is observed, as well as strong interactions with Ile305, Ala306, Ala307, and Thr310. This indicates that the presence of nonpolar groups on the B ring with this geometry still stabilizes the ligand to improve the strength of existing interactions in the scaffold, but also provide additional strong interactions with the nearby nonpolar residues and avoids electrostatic clashing with the nearby polar residues. An interesting feature of compound 3 is that the binding mode of the ligand is unique in comparison to that of the other ligands. While other high-performing ligands featured a stable binding mode with a low RMSF, compound 3 took longer to establish the final pose (as shown in the plot of ΔG vs time in Figure C and the ligand radius of gyration in Figure S4 ) and has RMSF of 1.60 Å ( Table S1 and Figure S3 ). The large RMSF indicates that the pocket struggles to accommodate large functionalities. The strength of the interactions with the benzofuran of compound 3 is diminished with most residues in comparison to other isoflavanone-adjacent ligands ( Figure ), including heme, where the π stacking is not as strong, though it does introduce a contact with Phe221 that is not present in other ligands. Along with this, the strong interactions made by the ethyl groups compensate and improve compound 3 in comparison to the other ligands. This indicates that the ligands benefit greatly from nonpolar substituents, but too much bulk causes the geometry of the benzofuran to change, which results in weakened interactions at the benzofuran. Compound 4 ( Figure D) features two chloride groups arranged with para geometry at positions 2 and 5. While it occasionally visits higher energy landscapes, the binding free energy remains relatively stable at −37.24 ± 0.03 kcal/mol, with a very low ligand RMSF of 0.87 Å and stable radius of gyration ( Table S1 and Figure S4 ). At the benzofuran, it forms strong interactions with Ile133, Phe134, Val373, and Met374, with diminished interactions at Arg115, Trp224, Val370, and Leu477. As for the residues near the B ring, compound 4 behaves similarly to compound 1 , with diminished interaction strengths at each residue and a clash introduced with Asp309. It features improved binding affinity in comparison to compound 1 because the interactions at the benzofuran are stronger, likely because the repulsive interactions from chlorine atoms to the pocket residues allow the other end of the molecule to adopt a more favorable geometry. Figure E shows compound 5 , which features three methyl groups at positions 2, 3, and 5. It maintains a relatively stable mode with a ligand RMSF of 1.04 Å and a consistent radius of gyration after approximately 150 ns ( Table S1 and Figure S4 ), but the binding affinity is somewhat erratic (consistent with changes in radius of gyration), averaging at −37.25 ± 0.03 kcal/mol, which is the same as compound 4 . Compound 5 features very strong interactions with Arg115 and Phe134, which are the strongest interactions at these residues for any ligand including ASD. It also features strong interactions with other residues at the benzofuran, namely, Ile133 and Met374, though it is diminished at other residues. The nonpolar functional groups at the B ring do improve the interactions with Ile305-Thr310, similar to compound 3 . Overall, this ligand features strong interactions at every part of the molecule without any clashes but also with a more favorable geometry of the benzofuran, which gives it the strongest binding affinity of any isoflavanone-adjacent ligand studied in this work. To contextualize compounds 1–5 into experimentally meaningful values, we have trained a QSAR model on a data set of experimental IC 50 values from the CHeMBL database and descriptors calculated by rdkit. The model was trained using LightGBM on 85% of the data set, and the remaining 15% was used as an external validation set. The model building is described in more detail in the Methods section. The model was used to predict the IC 50 values of the test set, and these predicted values were compared to the true values, resulting in a correlation coefficient of 0.72 ( Figure S5 ). The model was used to predict the IC 50 of compounds 1 – 5 , which are 1035 nM, 1222 nM, 3249 nM, 631 nM, and 3919 nM respectively. The model was also used to predict the IC 50 values of each isoflavanone-adjacent compound in the study ( Table S1 ).

In this work, we analyzed 77 isoflavanone-adjacent ligands first with docking and then with molecular dynamics to identify potential inhibitors of aromatase using an isoflavanone-like scaffold. The most favorable conditions for binding of ligands with this scaffold are nonpolar functional groups arranged with para geometry at positions 2 and 5. This resulted in several ligands with binding affinity similar to that of Letrozole, an FDA-approved aromatase inhibitor. The para geometry of positions 2 and 5 resulted in an improved binding mode at the benzofuran moiety, with additional improvements to binding affinity provided by nonpolar functional groups on the B ring. The pocket seems better accommodated for small functional groups, but bulky functional groups can also be favorable. This study identified five ligands with strong binding affinities. Given that four of the five ligands have functional groups at positions 2 and 5, it is clear that this geometry is favored for this scaffold. Furthermore, as seen by the polarity of the three strongest binding ligands, the global nonpolar character on the B ring is ideal. This is not surprising considering the nonpolar character of the pocket. Still, the actual identity of the functional groups is not the determining factor for the strength of binding. This is seen in the comparison of compounds 4 and 5 : both are globally nonpolar with the same binding affinity, with compound 4 having polar groups (Cl) and compound 5 having nonpolar groups (methyl). The 2,5 occupancy allows the benzofuran to occupy a favorable geometry, but the pocket struggles to accommodate bulky groups, as seen with compound 3 . Still, while the benzofuran ring geometry seems to be more important for binding than the geometry of the B ring, nonpolar functional groups are more capable of forming favorable local interactions with the enzyme than polar groups. This is evidenced by compounds 3 and 5 having generally strong interactions with the nonpolar residues near the B ring (such as Phe221, Ala306, Ala307, and Val370) while avoiding clashes with Asp309 that the ligands with polar groups do not ( Figure ). Compound 5 also has a methyl group at position 3, which interacts favorably with the nearby Ala306 and Ala307. Interestingly, when a methyl is added to position 3 on compound 4 , the binding affinity decreases to −32.63 kcal/mol, and the ligand RMSF increases to 1.84 Å ( Table S1 ). This indicates that the geometry of the 2,5 occupied B ring differs when the functional groups are polar or nonpolar, and introducing nonpolar groups to the polar ring geometry destabilizes the ligand binding. Of course, an object of great concern with the development of aromatase inhibitors is bioavailability, toxicity, and ease of synthesis. Table provides several metrics provided by the SwissAdme server for Letrozole, the scaffold, and each isoflavanone-adjacent ligand presented in Section . Each ligand features similar properties to Letrozole in terms of solubility, GI absorption, and bioavailability, as well as expected inhibition of other CYP proteins. The isoflavanone-like ligands tend to have higher logP values and only slightly higher synthetic accessibility scores. These ligands are predicted to be orally bioavailable and relatively easy to synthesize. Because these ligands are meant to inhibit aromatase, which is a cytochrome p450, it is expected that they would also target other cytochrome p450s, similar to Letrozole; however, because these ligands are derived from natural products, it is likely that they will exhibit fewer off-target effects with other enzymes. Furthermore, the high degree of importance regarding geometry, rather than the identity, of the functional groups on the B ring (see Results Section and 3.3 ) may increase the area of chemical space that can be explored to tune the synthetic accessibility and toxicity profile while still maintaining a strong binding affinity. Lastly, the isoflavanone-adjacent compounds are in compliance with Lipinski’s rule-of-five. Given the importance of the 2,5 occupancies, the improved interactions of global nonpolar character on the B ring, and the greater ability for the pocket to accommodate small functional groups over bulky groups, we believe that compounds 4 and 5 represent the best lead compounds to focus on for future analyses. Compound 4 is predicted to be more synthetically accessible and bind fewer CYP enzymes than compound 5 ( Table ). Additionally, compound 4 is predicted to have the lowest IC 50 value of the top 5 ligands (631 nM) in the QSAR model (Results, Table S1 ). Still, compound 5 does not clash with Asp309 and has a more consistent binding affinity ( Figure ). In this work, we have presented characteristics of isoflavanone-adjacent ligands that favor or disfavor binding to aromatase and provide examples of ligands with a high binding affinity on par with Letrozole, an FDA-approved aromatase inhibitor. Isoflavanones are a natural product, and as such, a similar scaffold could exhibit a stronger toxicity profile than the available FDA-approved aromatase inhibitors. Future studies should focus on experimental testing of the compounds as aromatase inhibitors. Additionally, while MM-GBSA is a fast and decently accurate method for calculating binding free energies, there are more accurate methods available that might provide a clearer understanding of the effect of functional groups. Lastly, modifications to this scaffold could be explored using computational approaches. The results of this study are expected to be useful in the further development of aromatase inhibitors based on natural product scaffolds.

Endometriosis is a chronic estrogen-dependent inflammatory disease that affects approximately 10% of reproductive-age women and is characterized by the presence of endometrial tissue outside the uterus. , Symptoms are often debilitating and include dysmenorrhea, dyspareunia, and chronic pelvic pain. , Currently there is no curative treatment, and management options are limited. Nonsurgical treatments aim to reduce inflammation and suppress the menstrual cycle, while surgical treatments aim to eliminate lesions or entire pelvic organs, with neither approach offering long-term relief. − Although a complete understanding of the etiology of endometriosis remains elusive, there are several marked differences in gene and protein expression between healthy and diseased tissue, − summarized in the review by Burney et al. A notable finding relevant to the present work is that diseased tissue is self-estrogen producing through a positive feedback cycle, which we summarize in Figure A. Aromatase is the key enzyme for the synthesis of estrogen. Notably, aromatase is undetectable in healthy tissue, but diseased tissue features a marked upregulation of the enzyme. , , The positive feedback cycle is discussed in more detail in the review by Bulun et al. (A) Primary reaction catalyzed by aromatase (red outline) and it is participation in the positive feedback cycle and (B) structure of the flavone scaffold, isoflavanone scaffold, the isoflavanone derivative presented by Bonfield et al., the isoflavanone-adjacent structure presented by Mirzaie et al., and the scaffold of the Mirzaie structure. Endometriosis can be successfully treated with GnRH analogues, but the cumulative recurrence rate is more than 50%, and can be as high as 75%. It is possible that this failure can be attributed to the continued increased production of estrogen in the lesions during treatment, and that interruption of the estrogen production with an aromatase inhibitor may extend the duration of remission. In fact, the use of aromatase inhibitors for this purpose has previously been successful. − The main drawback of aromatase inhibitors (AIs) is the side effects, which commonly include hot flashes, weight gain, insomnia, joint and muscle aches, loss of libido, mood changes, vaginal dryness, and loss of bone density. , − Some side effects of third-generation aromatase inhibitors (exemestane, letrozole, and anastrozole) may be intrinsic to estrogen regulation; however, it is possible that new AIs derived from natural product scaffolds will mitigate negative side effects and have greater synthetic accessibility. Flavonoids are a class of natural products that have been extensively studied as aromatase inhibitors, − most recently in an experimental and computational study that identified submicromolar dual-acting ligands that can inhibit both aromatase and estrogen receptor β Previous studies have focused primarily on two flavonoids: flavone ( Figure B) and flavanone (saturated C2–C3 bond). According to a recent study, natural flavones show no significant activity as aromatase inhibitors, but the structure–activity relationship of synthetic variants of the flavone scaffold have been reviewed for their potential aromatase inhibitory activity. Isoflavanones ( Figure B) are a subgroup of flavonoids that are capable of inhibiting aromatase. Bonfield et al. showed that the nonplanarity of the isoflavanone scaffold can introduce an additional enzyme-ligand interaction compared to the planar isoflavone. Several isoflavanone derivatives with low micromolar IC 50 values were identified, the most promising of which was designated 2a (shown here in Figure B). Of the 26 isoflavanones tested, only one (2a) exhibited a submicromolar potency (0.26 μM). Using 2a as the reference, Mirzaie et al. conducted a structural similarity search and predicted by molecular docking that 2a_1, as well as several other compounds with the same scaffold (shown here in Figure B), displayed improved binding affinity to the aromatase active site. This isoflavanone-adjacent scaffold features a benzofuran and a benzene linked by a ketone. To our knowledge, despite its promise, this unique scaffold has not been explored since. Because the potencies of isoflavanones are generally weak to moderate (even after modification ), and because synthetic variants of flavonoid scaffolds have previously shown promise, , we believe that exploration of these similar scaffolds is worthwhile. It is possible that the natural product adjacent compounds can improve potency without sacrificing a strong safety profile. We have performed a structural similarity search of 2a_1 using the SciFinder database to identify existing compounds with this scaffold ( Figure B), and we used molecular docking (AutoDock ) and molecular dynamics (AMBER ) to characterize the effects of the position and identity of functional groups on the scaffold to the binding affinity. Binding free energy is estimated using the MM-GBSA method, and the five ligands with the lowest binding free energy are analyzed in detail. This methodology and workflow are like other published works but differ in the execution. For example, select benzoxazole analogues were characterized using similar computational methodology following an experimental cytotoxicity screening. We employ an in silico screening of compounds and more extensive computational characterization. Additionally, Ziprasidone was identified as a promising lead compound following a screening of known bioactive molecules with molecular docking and characterization with molecular dynamics and free energy analyses. We employ this workflow focusing on a data set of variants to one scaffold to map the pharmacophore and identify promising compounds for further experimental analysis. Computer-aided drug design (CADD) methods are particularly attractive in the early stages of drug discovery. They are significantly more accessible and cost-effective than experimental techniques, enabling rapid screening and prioritization of candidate compounds without the need for synthesis, acquisition of physical materials, and complex experimental design. , Furthermore, computational approaches offer a deeper understanding than experimental methods are easily capable of; particularly an atomistic understanding of the key protein–ligand interactions. , Computational approaches also allow for the virtual evaluation of compounds that may not be readily commercially or synthetically available. Because compounds with the isoflavanone-adjacent scaffold are not readily available for purchase, computational approaches are especially useful here. We believe that the results of this study provide a strong foundation for further experimental evaluation toward the development of new AIs that are structurally similar to natural products.