Pharmacotherapeutics and Molecular Mechanism of Phytochemicals in Alleviating Hormone-Responsive Breast Cancer

Breast cancer (BC) is the leading cause of death among women worldwide devoid of effective treatment. It is therefore important to develop agents that can reverse, reduce, or slow the growth of BC. The use of natural products as chemopreventive agents provides enormous advantages. The aim of the current investigation is to determine the efficacy of the phytochemicals against BC along with the approved drugs to screen the most desirable and effective phytocompound. In the current study, 36 phytochemicals have been evaluated against aromatase to identify the potential candidate drug along with the approved drugs employing the Cdocker module accessible on the Discovery Studio (DS) v4.5 and thereafter analysing the stability of the protein ligand complex using GROningen MAchine for Chemical Simulations v5.0.6 (GROMACS). Additionally, these compounds were assessed for the inhibitory features employing the structure-based pharmacophore (SBP). The Cdocker protocol available with the DS has computed higher dock scores for the phytochemicals complemented by lower binding energies. The top-ranked compounds that have anchored with key residues located at the binding pocket of the protein were subjected to molecular dynamics (MD) simulations employing GROMACS. The resultant findings reveal the stability of the protein backbone and further guide to comprehend on the involvement of key residues Phe134, Val370, and Met374 that mechanistically inhibit BC. Among 36 compounds, curcumin, capsaicin, rosmarinic acid, and 6-shogaol have emerged as promising phytochemicals conferred with the highest Cdocker interaction energy, key residue interactions, stable MD results than reference drugs, and imbibing the key inhibitory features. Taken together, the current study illuminates the use of natural compounds as potential drugs against BC. Additionally, these compounds could also serve as scaffolds in designing and development of new drugs.


Introduction
Cancer is the primary cause of death globally [1], and breast cancer (BC) is the leading cause of cancer mortalities among women [2,3]. Currently available treatments include radiation therapy, chemotherapy, surgery, immunotherapy, and hormone therapy; however, it still lacks effective treatment.
Additionally, the currently available medication is ineffective and induces toxicity thus causing a major hindrance for effective treatment [4]. Adding to these, the acquired resistance that is prone to mutations generated during the cancer treatments and the resistance rendered because of the minor heterogenic subpopulation may enhance the ineffectiveness of the treatment [5]. This warrants the development of more efficient drug formulations with less adverse effects and correspondingly can slow the growth of tumours or reverse the process. In recent years, natural compounds such as plant extracts are being studied for their anticarcinogenic properties. The advantages and importance of natural compounds are greater over synthetic compounds as they are less toxic relative to the concentration of the compounds used and the cellular or the physiological environment. Additionally, they have high selective biological actions [6], easy to extract [6], and their vast abundance. Furthermore, it is reported that over 35% of the cancer cases can be addressed by varying lifestyle and dietary habits [4,7,8], and phytochemicals are potential candidates in suppressing them [9] including BCs [10][11][12][13][14].
Phytochemicals are largely antioxidants in nature at lower concentrations and under favourable cellular conditions that effectively prevent the oxidation of other molecules that have an ability to produce free radicals and thus protect the body. On the contrary, certain phytochemicals tend to show prooxidant activities when used at low pH and high concentrations. These free radicals bear an unpaired electron in their outermost atomic orbital and can either donate or accept an electron from other molecules [15]. Although oxygen is an important element for life, however, under certain conditions, causes transformation into certain chemical compounds called reactive oxygen species (ROS) [16,17] which are highly reactive and unstable [17]. These are further able to cause damage to biologically essential macromolecules such as DNA, carbohydrates, proteins, and lipids [18].
The generated free radicals promote BC [19,20] besides contributing to various diseases [21,22]. However, when an imbalance exists between the genesis of free radicals and their degradation that subsequently leads to oxidative stress (OS) resulting in several "oxidative stress-related diseases" including cancer [23]. Additionally, the treatments offered for cancers such as chemotherapy and radiotherapy enhance the OS condition. The ROS cause damage to genes leading to genetic instability and are involved as intermediaries to certain signals further contributing to cancer progression and angiogenesis. Delineating on the role of free radicals as key players in contributing towards BC, a host of mechanisms have been identified. Free radicals may induce mutations in DNA primarily caused due to oxidation at a frequency of 10 4 lesions/cells/day in humans [24]. This leads to protein alterations considerably impeding its biological activity, thereby leading to genetic instability [25]. Additionally, free radicals increase mitogenic signal intermediaries and are involved in protein remodelling, enhancing proliferation, senescence, cell apoptosis, and autophagy [19].
Plant-derived compounds are ascribed to be the rich sources of antioxidants and can efficiently fight against the deleterious effects of free radicals. Mechanistically, these nonenzymatic antioxidants communicate with the free radicals and thus regulate the deterioration of the biologically important compounds. They exert their activities by several mechanisms [26]. Moreover, the natural antioxidants offer several advantages such as being economical and easy availability complemented by generating low toxic effects when administered at specific physiological doses as reported earlier [27][28][29]. Besides, several reports additionally illuminate the use of phytochemicals in dispelling free radical [30][31][32][33][34].
Phytochemicals have been foremost compounds in possessing antioxidant activities and are embarked as anticancer products. They either exert their action as free radical scavengers or act as quenchers of singlet oxygen in addition to their reducing capabilities employing their bioactive compounds [35]. A majority of the population from Asia, Latin America, and Africa consider the use of phytochemical as drugs [36,37]. Moreover, phytochemicals adapt various mechanisms such as reducing oxidative stress, killing the rapidly dividing cells, inducing programmed cell death, angiogenic hindrance, and targeting the molecular factors that are abnormally expressed. Accordingly, the pleiotropic behaviour of the phytochemicals endorses them as an excellent alternative therapeutics against cancer.
The present study evaluates 36 phytochemicals against BC drug target, aromatase, to probe into the best prospective drug candidate that represents the pharmacophore features and further to determine the mode of their inhibitory mechanism computationally as illustrated in Figure 1.

Ligand Selection and Preparation.
For the current investigation, 36 phytochemical antioxidants have been identified from different literature reports [9,17,[38][39][40]. Their 2D structures were drawn on BIOVIA Accelrys (http://accelrys. com/products/collaborative-science/biovia-draw/) and were imported onto DS to subsequently generate their 3D structures ( Figure 2). The ligands were preprocessed by rectifying their bond angles and bond orders and were subsequently minimized on DS employing CHARMm force field.

Protein Selection and Preparation.
For the current investigation, the validated drug target, aromatase, was imported from protein data bank (PDB) with the PDB code 3EQM. The protein was prepared by removing all the water molecules and adding hydrogen atoms. The orientation of the histidine residues was placed in accordance with the crystal structure. The active site was defined at 10 Å around the cocrystal ligand (hereinafter it is referred to as cocrystal) and the key residues were determined as Arg115, Ala306, Asp309, Val370, Leu372, Met374, and Leu477. The prepared protein and the ligands were forwarded to molecular docking exploration to evaluate the binding affinities between the protein and ligands and further to deduce the ideal binding mode.

Molecular Docking Studies and Binding Energy
Calculations. Molecular docking is one of the increasingly popular tools in the modern day drug discovery that imparts knowledge on small molecules that accommodate well in the proteins active site [41]. Generally, two steps govern the success of the docking protocol; the location of the ligand and its orientation within the binding pocket and further assessing the binding interactions between the key residues and the chosen ligands. To understand the potency of the phytochemicals, two approved drugs exemestane and letrozole were chosen and were labeled as reference. The Cdocker programme implemented on the DS was recruited and the results were analyzed based upon the -Cdocker interaction energies. The higher the -Cdocker interaction energies, the greater the binding affinity between the protein and the ligand. Cdocker is a grid-based docking approach that utilizes CHARMm in which the protein is rigid while the ligands were allowed to move. Mechanistically, Cdocker utilizes high temperature molecular dynamics to arbitrarily generate ligand conformations and is subsequently redirected towards the binding site to facilitate the formation of poses using the simulated annealing method [42]. In order to ensure the accuracy of the docking protocol, the cocrystal was docked into the proteins active site. The cocrystal and the docked pose have resulted in an acceptable root mean square deviation (RMSD) of 0.7 Å (Supplementary Figure 1A) that assures the docking parameters. Therefore, the same parameters were chosen for the investigation. The selected 38 ligands were docked into the proteins active site. Each ligand was allowed to generate 50 conformations and were subsequently clustered to deduce the best binding pose. Additionally, the binding energies were computed which renders information on the approximate binding energy between the protein and the ligand. This investigation was accomplished with calculate binding energies protocol available with the DS and was conducted employing the equation energy binding = energ y complex − energy ligand − energy receptor . Based upon the highest dock scores and interactions with the crucial active site residues of the best binding conformations, they were probed for possessing the inhibitory features that are predominantly essential for repressing the target enzyme exploiting the structure-based pharmacophore approach.
2.4. Structure-Based Pharmacophore Generation. In order to probe into the inhibitory chemical features of the chosen small phytochemical molecules, the structure-based pharmacophore approach (SPB) was employed. SBP is regarded as one of the finest methods adapted in the field of drug discovery that exploits the features of the protein ligand interaction. SBP was generated considering the protein aromatase bearing the protein data bank (PDB) code: 3EQM with androstenedione as cocrystal. The residues in close proximity to the cocrystal were considered during the generation of the pharmacophore structure. Accordingly, the clean protein protocol was enabled to check for any gaps existing in the protein. Subsequently, the receptor ligand pharmacophore generation module embedded with the DS was chosen and the parameter for maximum pharmacophore was opted as  10 with minimum and maximum features as 4 and 6, respectively. Additionally, the maximum charge distance was opted as 5.6 Å with an interfeature distance of 2.0 Å, while retaining the maximum hydrogen bond distance and maximum hydrophobic distance as default.

Validation of the Generated Pharmacophore by
Güner-Henry Method. Pharmacophore validation is one of the major criteria in assessing the robustness of the generated pharmacophore in identifying the active compounds from the inactive compounds. For its accomplishment, the decoy set method of validation was carried out formulating an external database (D) of 1229 compounds with 15 active compounds (A) in it. The ligand pharmacophore mapping module accessible with the DS with rigid fitting method was initiated. The results were evaluated based upon the enrichment factor (EF) and the goodness of fit (GF) score as described earlier [43] employing the formula where EF refers to the enrichment factor, Ha refers to the  2.6. Retrieving the Compounds with the Pharmacophore Features. The 36 phytochemicals that have displayed the highest dock scores than the approved drugs were subjected to pharmacophore mapping using the ligand pharmacophore mapping module embedded with the DS. Logically, it is assumed that the small molecules that obey to all the features of the pharmacophore might be imbibed with the key inhibitory features recruiting the rigid fitting method.

Molecular Dynamics Simulation
Studies. The best fitted compounds resulted from the docking studies along with the reference compounds were forwarded to molecular dynamics (MD) simulations to assess the stability and further affirm the docking results employing GROningen MAchine for Chemical Simulations v5.0.6 (GROMACS). The simulation run was performed for 10 ns using CHARMm ff [44]. The ligand topologies were obtained employing SwissParam [45][46][47], and a dodecahedral water box consisting of (transferable intermolecular potential 3P) TIP3P water model with a thickness of 1 nm was subsequently generated and neutralized with counter ions. The system was minimized through 1,000 steps using the steepest descent algorithm to expel bad contacts. Following this, the equilibration process was conducted using constant number, volume, and temperature (NVT) and constant number, pressure, and temperature (NPT) refraining the protein backbone and allowing the solvent molecules and the counter ions to flex. The NVT was executed for 1 ns at 300 K with V-rescale thermostat employed to maintain constant temperature. The NPT was performed for 1 ns at 1 bar using Parrinello-Rahman barostat [48]. The bonds of heavy atoms were restrained recruiting the LINCS algorithm [49]. The long-range electrostatic interactions were measured using particle Mesh Ewald (PME) [50] method. The short-range interactions were computed selecting a cut-off value of 12 Å. The MD was conducted under periodic boundary conditions to escape edge effects with a time step of 2 fs saving the coordinate data for every 1 ps. The results were evaluated employing visual molecular dynamics (VMD) [51] and DS.

Molecular Docking Studies and Binding Energy
Calculations. Molecular docking studies disclosed that all the selected natural antioxidants have rendered higher dock scores than both the reference compounds. The reference compounds have generated a dock score of 16.86 and 20.03, respectively, for exemestane and letrozole. Conversely, phytochemicals have rendered dock scores that exist between 18.61 and 55.14 demonstrating their therapeutic usability. Among them, only one compound, ascorbyl palmitate, generated a score of 18.61, which is above exemestane and below letrozole (Table 1). Furthermore, curcumin, capsaicin, rosmarinic acid, and 6-shogaol have demonstrated highest -Cdocker interaction energies, representing the quintessential binding modes and additionally have anchored at the active site with the key residues. Moreover, the binding energies reinforce the results generated by the docking. Most of the phytochemicals have computed lower binding energies than the reference compounds correspondingly implying greater affinity towards the target (Table 2). However, the phytochemical safrole and p-cymene have generated a different result; safrole has displayed a greater binding energy than exemestane and p-cymene has rendered a higher binding energy than both references (Table 2). Additionally, to assess the inhibitory features of the phytochemicals in accordance with the pharmacophore, these compounds were promoted to map with the SBP.
3.2. Structure-Based Pharmacophore Generation. The structure-based pharmacophore has generated only one pharmacophore consisting of three features that are complimentary to the key residues. One hydrogen bond acceptor feature was formed between the ligand and the residue Met374. The two hydrophobic features were noticed between the critical residues Trp224 and Val370, respectively ( Figure 3).

Güner-Henry Method of Pharmacophore Validation.
Decoy set of validation was fundamentally executed to determine the ability of the pharmacophore in retrieving the active compounds from a given database. Accordingly, the generated pharmacophore has mapped to 20 compounds (Ht) with 15 active compounds (Ha). Correspondingly, the EF and the GF values were computed as 61.45 and 0.77, respectively. The GF score critically evaluates the quality of the pharmacophore and can lie between 0 and 1 representing the null and ideal model. Since the current pharmacophore has generated a GF score of 0.77 is near to the value 1, this model can be considered as a good model and can be utilized for further studies (Table 3).

Retrieving the Compounds with the Pharmacophore
Features. According to Ehrlich (1909) a pharmacophore is defined as "a molecular framework that carries (phoros) the essential features responsible for a drug's (pharmacon) biological activity" [52]. In another definition according to IUPAC, a pharmacophore model is "an ensemble of steric and electronic features that is necessary to ensure the optimal supramolecular interactions with a specific biological target and to trigger (or block) its biological response" [53]. Therefore, a SBP was generated and subsequently, the phytochemicals were mapped against it.
The ligand pharmacophore mapping has rendered that only 19 compounds exhibited the features imbibed by the pharmacophore including the approved drugs (Table 1). From these compounds, the top four compounds with highest dock score exhibiting the interactions with the key residues were scrupulously examined further and were upgraded to molecular dynamics simulations to assess their behaviour at the atomic level and interpret the molecular mechanism of inhibition.

Molecular Dynamics Simulation
Studies. MD studies elaborate on the dynamic behaviour of the ligand at the proteins active site thereby additionally affirming the molecular docking results. MD run for 10 ns was initiated with the best dock poses for top-ranked phytochemicals (4) along with the reference compounds (2). A total of six systems were subjected to MD simulations and the results were read as RMSD and potential energy profiles. The RMSD of the four antioxidants and the reference compounds were computed to be below 0.2 nm (Figure 4 and Supplementary Figure 1C). The average RMSD of the reference compounds was calculated to be 0.13 nm and 0.12 nm for exemestane and letrozole, respectively. The antioxidants on the other hand have rendered 0.13 nm, 0.13 nm, 0.12 nm, and 0.12 nm correspondingly for curcumin, capsaicin, rosmarinic acid, and 6-shogaol. Introspecting the RMSD plots, it was revealed that marginal variations were noticed during the initial MD  Figure 1D). Accordingly, the representative structures were extracted from last 4 ns and were correspondingly superimposed. Subsequently, it was discovered that the four phytochemicals have occupied the binding site in the similar fashion as that of the reference compounds ( Figure 5). Elaborating on the molecular interactions, it was observed that the reference compounds have anchored to the protein with two key residues, Arg115 and Met374, respectively, conferred with a relatively acceptable bond length (Table 4). Furthermore, exemestane has interacted with the protein through π-alkyl bonds involving the residues such as Ile133, Val370, and Val373. Other residues have participated in the van der Waals interactions including Phe134, Ile305, Ala306, Asp309, Thr310, Leu372, and Leu477 which help in proper positioning of the compound in the protein active site ( Table 4, Figure 6(a), and Supplementary Figure 2A). Letrozole has displayed two hydrogen bonds with important residues Arg115 and Met374 contributed by an acceptable bond length ( Table 4). The benzene ring A of letrozole has participated in anchoring to the protein with residue Val370 as π-alkyl bond. The pentane ring has formed the π-alkyl bond with Ile133 residue. Phe134, Phe221, Trp224, Asp309, Thr310, Val369, Leu372, and Val373 have formed the van der Waals interactions. The benzene ring B has additionally participated in π-cation interaction with Arg115 enabling the accommodation of the ligand in its binding pocket (Table 4, Figure 6(b), and Supplementary Figure 2B).
Rosmarinic acid displayed higher number of hydrogen bond interactions involving the catalytic residues Arg115, Glu302, Met374, and Leu477 with an acceptable bond length (Table 4). Several residues aid in appropriate residing of the ligand conferred by van der Waals interaction such as Phe134, Phe221, Trp224, Ile305, Ala306, Asp309, Thr310, Val369, Leu372, Val373, and Ser478, while only one π-alkyl bond was resulted through the benzene ring A of the ligand and Val370 residue (Table 4, Figure 6(e), and Supplementary Figure 3C). 6-Shogaol exerted its inhibitory activity by interacting with two key residues, Arg115 and Met374 correspondingly. The benzene ring has participated in anchoring to Ile133 and Val370 residues by π-alkyl bond. Furthermore, residues Phe134, Trp224, Ile305, Ala306, Asp309, Val373, Leu372, Leu477, and Ser478 have secured the ligand firmly at the binding site of the protein (Table 4, Figure 6(f), and Supplementary Figure 3D). From the above results, it can be noted that the phytochemicals have apparently interacted with a higher number of residues upon comparison with reference compounds. Therefore, the above findings guide us logically to infer that the prospective drug candidates could be deemed as alternative therapeutics against BC.

Discussion
BC has been described as one of the major impediments caused to the normal well-being and is one of the undefeated diseases. Consequently, new strategies to combat this disease have been on a high demand. Nature has been a rich source of medicines for various diseases [54][55][56][57], besides demonstrating high antioxidant activities. For the current investigation, 36 phytochemical antioxidants have been assessed computationally to discover the highly effective compounds against BC protein, aromatase. Molecular docking results have disclosed that all the 36 phytochemicals have generated higher Cdocker interaction energies (Table 1) and lower binding energy than the reference compounds (Table 2). These phytochemicals have secured a similar binding mode as was seen with the reference and the cocrystal (Supplementary Figure 1B). Structurally, all the selected 36 phytochemicals have exhibited diversity. The relatively small compounds had more freedom to be accommodated within the active site and adapted a linear conformation, while the larger compounds have adapted an inverted "C" (Ͻ) conformation to be accommodated firmly within the proteins active site. Expounding on the molecular interactions has imparted information on the key   Table 1). An interesting observation was noticed with the phytochemical p-cymene that was devoid of oxygen and this could probably be the reason for showing no interaction with Met374 residue. This drives us to perceive that the effectiveness of the ligand in inhibiting the aromatase lies if it bears an oxygen atom. Additionally, it is worth discussing regarding the interactions with residues Phe134 and Val370.
The residue Val370 has interacted with all the phytochemicals represented by pi-alkyl bond. However, the interaction Val370 was missing with the ascorbyl palmitate. Phe134 on the other hand displayed interactions with all the phytochemicals by van der Waals interaction. However, this interaction was missing with phytochemicals carnosol, flavone, myristicin, carvacrol, and p-cymene. Nevertheless, these compounds have interacted with Phe134 showing pi-alkyl bond. Additionally, the phytochemical ascorbyl palmitate lacked the interactions with these two residues and further rendered a lower -Cdocker interaction energy of 18.61. This guides us to predict that the three residues Phe134, Val370, and Met374 are crucial in showing enhanced ability to interact with the protein.
Interesting results were obtained upon performing the search for compounds with key inhibitory features for all the 36 compounds via the SBP. It was observed that only 19 compounds have aligned to the pharmacophore, although all the 36 compounds have demonstrated a higher dock score. This finding bespeaks the therapeutic ability of 19 compounds against the enzyme aromatase.
The study further focuses on the compounds that have generated highest dock scores displaying interactions with key residues. The top scored phytochemicals were curcumin, capsaicin, rosmarinic acid, and 6-shogaol. Rosmarinic acid has rendered four hydrogen bonds and curcumin displayed three hydrogen bonds while capsaicin and 6-shogaol have represented two hydrogen bonds each. Curcumin was the only phytochemical whose benzene ring B has interacted with Agr115 through π-cation bond similar to the reference compound letrozole. This makes curcumin a unique compound that also exhibited the highest dock score. Additionally, it was found that several amino acids have favoured proper positioning of the ligands at the protein active site firmly (Table 4 and Figure S3). Complementing on the results, it was observed that the prospective drugs have accommodated in the similar fashion as the reference drugs and were shown to interact with the key residues. Moreover, the MD results have disclosed that the RSMD and the potential energy profiles of all the systems was relatively stable after 2500 ps (Supplementary Figures 1C  and 1D). These identifications guide us to secure the four antioxidants as alternative therapeutics against BC.
The target for the present study, aromatase, is a product of CYP19A1 gene that catalyzes the conversion of androstenedione to estrogen and testosterone to estradiol  [58]. The overexpression of this enzyme promotes hormone-responsive BC. Moreover, inhibition of aromatase expression can prompt selectivity in blocking the production at tumour site. There are several aromatase inhibitors and are broadly grouped into type I (steroidal) and type II (nonsteroidal) inhibitors [59]. In order to develop and discover potential chemotherapeutic compounds, we examined the role of phytochemicals as an alternative. For the efficient execution of the current research and to garner sufficient knowledge on the phytochemical therapeutic ability and their mechanism of action against aromatase, two reference compounds exemestane (type I) and letrozole (type II) inhibitors were chosen. The computational results led us to draw important information involved in the inhibitory mechanism. The interactions of the ligands with the residues Phe134 and Val370 were determined to be imperative in inducing the inhibitory mechanism.
The sulphur-containing residue, Met374, was significant in demonstrating strong hydrogen bond interaction with all the ligands and the reference compounds (Table 4 and Supplementary Table 1). These findings drive us to infer the importance of the residues Phe134, Val370, and Met374 in contributing to the inhibitory activity as illustrated in Figure 7. Furthermore, these three residues are located as a "triad" at the active site facilitating to lock the ligand firmly from three sides and such a triad clamp has been demonstrated by a majority of the chosen phytochemicals, reflecting their substantial responsibility in alleviating the BC.

Conclusion
To develop and design a drug with least side effects has been one of the challenging avenues in identifying potential BC inhibitors. It can therefore be concluded that the use of plant-derived compounds rich in antioxidant potential that are characterized by excellent anticarcinogenic activities is an ideal choice to combat BC. Our findings demonstrate the superior quality of the phytochemicals over the known drugs. Additionally, based upon the obtained results, we propose the use of compounds curcumin, capsaicin, rosmarinic acid, and 6-shogaol that can represent as valuable compounds to fight against the BC. Additionally, they can serve as novel, fundamental scaffolds paving way for designing new drugs.

MD:
Molecular dynamics DS: Discovery Studio OS: Oxidative stress ROS: Reactive oxygen species 2D: Two dimensional 3D: Three dimensional PDB: Protein data bank RMSD: Root mean square deviation GROMACS: GROningen MAchine for Chemical Simulations NVT: Constant number, volume, and temperature NPT: Constant number, pressure, and temperature PME: Particle mesh Ewald LINCS: LINear constraint solver VMD: Visual molecular dynamics.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
We wish to confirm that there are no known conflicts of interest associated with this publication, and there has been no significant financial support for this work that could have influenced its outcome.  Figure 7: Proposed mechanism of phytochemical inhibition forming a triad.

Supplementary Materials
Supplementary Figure