Molecular docking and MD simulation approach to identify potential phytochemical lead molecule against triple negative breast cancer

Gene expression profiling of TNBC microarray datasets

A thorough literature mining effort encompassing all eligible studies on gene expression in TNBC was conducted. The search involved querying the Gene Expression Omnibus (GEO) datasets. Gene expression profiling was performed using GEO2R to identify significantly upregulated genes. Figure 1 presents an overview of the methodology.

Figure 1. Overview of methodology.

During the literature mining process, a microarray dataset was obtained from the NCBI GEO repository using the accession number GSE45498 annotated in the GPL16299 platform. This dataset encompasses 40 samples from healthy normal tissues, 160 from individuals with cancer, and 54 from metastatic cases. NGS datasets were obtained from the NCBI GEO repository using accession number GSE214101 annotated in the GPL20301 platform. This dataset included 24 samples derived from the MDA-MB-231 and MDA-MB-436 cell lines. Gene expression profiling values underwent log (base2) transformation and percentage shift normalization was applied. To assess the differences in gene expression between normal and diseased samples, the fold change for each gene was individually calculated. A threshold of 1.25-fold change was used to categorize genes as being upregulated. Gene expression profiling followed the protocol reported previously.¹⁰

Study of protein–protein interactions

The selected genes were subjected to the Bisogenet plug-in of Cytoscape to identify protein-protein interactions of all genes differentially regulated in TNBC. STRING is an open-source bioinformatics platform integrated in Cytoscape, designed for the study of both predicted and known protein-protein interactions. This database gathers, evaluates, and integrates information on protein-protein interactions from all publicly available sources. Additionally, it augments these data with computational predictions.¹¹ These interactions encompass both indirect (functional) and direct (physical) associations.¹² The genes were uploaded and a string network was built. Molecular Complex Detection (MCODE) detects Protein-Protein Interactions subnetworks and highly interconnected clusters within the PPI network.¹³ PPI networks were broken down into top-ranked dense cliques (sub-clusters) using the MCODE plugin. The top-ranked dense clique was selected for further analysis.

Building a library of phytochemicals with anti-breast cancer activity

Phytochemicals are naturally occurring biologically active chemical compounds found in plants that serve as medicinal ingredients and nutrients, offering health benefits to humans.¹⁴ Many natural products and their analogs have been identified as potent anticancer agents and the anticancer properties of various plants and phytochemicals.¹⁵ Phytochemicals were identified through a systematic literature search indicating anti-breast cancer activity were selected, and their 3D structures in SDF format were retrieved from PubChem database. Subsequently, phytochemicals that did not conform to Lipinski’s rule of five were excluded, and the remaining compounds were subjected to further analyses.

Virtual screening

Understanding the fundamental principles governing how protein receptors recognize, interact, and form associations with molecular substrates and inhibitors is crucial for drug discovery. PyRx v0.8 software¹⁶ with an inbuilt AutoDock Vina 1.2.5¹⁷ for molecular docking was used to scan phytochemicals conforming to Lipinski’s rule of 5. AutoDock Vina uses a semi-empirical free-energy force field to predict the binding free energies of small molecules to macromolecular targets.

The human Androgen Receptor (PDB ID: 1E3G) was sourced from the RCSB Protein Data Bank. Initially, the protein structure underwent a curation process to remove any crystallographic water molecules and heteroatoms that might interfere with docking simulations. Subsequently, energy minimization was performed using UCSF Chimera vs 1.54 (https://www.cgl.ucsf.edu/chimera/) to optimize the geometry of the protein. The steepest descent algorithm was applied for 100 steps, which is a common approach to relieve steric clashes and achieve a more stable conformation. Partial charges were then assigned to the protein using the AMBER ff14SB force field, which is well known for accurately modeling protein dynamics and interactions. The co-crystallized ligand metribolone (R18) was used as the control, and the ligands were docked at its active site.

ADMET - ProTox II

The development of high-quality in silico ADMET models will enable compound efficacy and druggability features to be optimized concurrently, thereby improving the overall quality of drug candidates.¹⁸ ProTox-II was used to experimentally validate the chemical toxicity and their combination. It uses machine learning models, the most common features, pharmacophore-based, fragment propensities, and chemical similarity to forecast different toxicity endpoints.¹⁹ Based on the virtual screening results, the top ten phytochemical compounds were chosen for ADMET analysis.

Induced fit docking

Induced fit docking was carried out using Schrodinger vs. 2020.3, which takes into account the flexibility of both the protein receptor and ligand, allowing for conformational changes to occur upon binding. The energy-minimized ligands were saved in PDB format for compatibility with the Schrodinger software, and the partial charges of the ligands were assigned, such as Gasteiger charges, which estimate the distribution of charges on the molecule based on its structure. Similarly, the protein charges may also be assigned using OPLS_2005 force fields to accurately capture its electrostatic properties. The grid box is a crucial parameter in docking simulations, as it defines the search space where the ligand can orient itself around the protein receptor. The dimensions of the grid box are typically specified in terms of the number of grid points along each axis (nx, ny, nz) and the grid spacing (Å) around the binding cavity residues LEU701, LEU707, MET742, MET745, ARG752, MET780, MET787, ALA748, LEU880, LEU873, PHE876, MET895, ILE899, THR877, GLN774, PHE764, LEU746, GLY708, GLN711, TRP741, ASN705. The dimensions were set to (58, 64, and 52 Å), providing a sufficient volume to explore potential binding modes of the ligand within the protein’s active site with a charge cutoff polarity set for a charge cutoff of 0.25 Å.

Molecular dynamics simulation

Molecular dynamics (MD) simulations were conducted for the docked complex of the human Androgen Receptor with the best-docked molecule, employing Schrodinger Desmond 2020.1.²⁰ The OPLS-2005 force field,²¹ along with an explicit solvent model using SPC water molecules,²² were employed in this system. The simulation was performed in a periodic boundary solvation box with dimensions of 10 × 10 × 10 Å. To neutralize the charge, Na+ ions were added, and a 0.15 M NaCl solution was added to mimic the physiological environment. The initial equilibration was carried out using an NVT ensemble for 10 ns to allow the system to relax over the protein-ligand complexes. Subsequently, a short run of equilibration and minimization was conducted using an NPT ensemble for 12 ns. The NPT ensemble utilized the Nose-Hoover chain coupling scheme²³ with a temperature set at 37 °C, relaxation time of 1.0 ps, and pressure maintained at 1 bar in all simulations. A time step of 2 fs was used.

Pressure control was achieved using the Martyna-Tuckerman-Klein chain coupling scheme²⁴ with a relaxation time of 2 ps. The long-range electrostatic interactions were calculated using the particle mesh Ewald method,²⁵ and the Coulomb interaction radius was fixed at 9 Å. A RESPA integrator with a time step of 2 fs was used for each trajectory to calculate the bonded forces. The final production run was extended for 100 ns for the Human Androgen Receptor with the best-docked molecule complex. To track the stability of the MD simulations, a variety of parameters were computed, including the number of hydrogen bonds, radius of gyration (Rg), root-mean-square fluctuation (RMSF), and root-mean-square deviation (RMSD).

Binding free energy analysis

Molecular Mechanics Generalized Born Surface Area (MM-GBSA) approaches are less computationally intensive than biochemical free energy methods and more precise than most molecular docking scoring systems. This method is useful for predicting the binding free energy in molecular systems. MM-GBSA is a useful technique for comprehending the impact of mutations on large biomolecular systems.²⁶ Biomolecular research has been utilized in investigations of protein folding, protein-ligand binding, protein-protein interactions etc.²⁷

The MM-GBSA approach was used to determine the binding free energies of the ligand-protein complexes. The MM-GBSA binding free energy was computed using the Python script thermal mmgbsa.py in the simulation trajectory with the VSGB solvation model and OPLS5 force field over the last 50 frames with a 1 step sampling size. The binding free energy of MM-GBSA (kcal/mol) was calculated using the additivity principle, wherein the differences in free energies, GBSA solvation energies, and surface area energies of ligand-protein complexes compared to their respective total energies of them individually were calculated.

Differentially expressed genes (DEGs) analysis

Gene expression in TNBC and normal microarray datasets was compared to assess the underlying molecular pathways driving TNBC, and further network analysis was performed. Boolean operators and relevant filters were used to filter the microarray datasets using the GEO2R. The Benjamini-Hochberg-Yekutieli approach was used to adjust the P-value for the DEGs, and only the top 10% of the upregulated genes (P-value < 0.05) were selected. Tables 1 and 2 display the list of elevated genes with LogFC > 1.25 and P-value < 0.05 in dataset GSE45498 and GSE214101, respectively.

The STRING tool was used to identify potential connections between DEGs in different tissues.¹² To build PPI networks, active interaction sources such as databases, co-expression, text mining, experiments, and species restricted to “Homo sapiens” were used, along with an interaction score greater than 0.4. The PPI network was displayed using Cytoscape v3.6.1 software as depicted in Figure 2.

Figure 2. Protein–protein interaction network where Androgen receptor (AR) is the central hub gene.

The MCode plugin was employed to identify the highly linked regions inside the PPI network, while the CentiScape plugin was utilized to calculate the network topology parameters. Using degree and betweenness as the primary parameters, hub genes were identified. A complete set of algorithms, called CentiScape, was used to analyze the centrality of the network nodes. It can calculate multiple centralities for weighted, directed, and undirected networks.²⁸ The human Androgen Receptor was determined to be an appropriate hub gene in the protein-protein interaction network consisting of DEG genes.

Virtual screening of phytochemical library

The human Androgen Receptor (hAR), covering the C-terminal amino acid residues (1E3G) with the co-crystallized ligand metribolone (R18), consists of 250 amino acid residues arranged in a three-layered α-helical sandwich structure. The ligand-binding pocket is located within the hydrophobic cavity formed by helices. A total of 1358 compounds were initially identified through systematic literature search, and their structures were retrieved from the PubChem database. Of these, only 543 compounds met the criteria outlined by Lipinski’s rule of five. These 543 compounds were then selected for the initial virtual screening against human Androgen Receptor using PyRx, and their binding affinities were tabulated²⁹ (refer to extended data Table S1). The top 50 ranked compounds were subjected to ADMET analysis on the ProTox II server. Only the top 10 ranked compounds that showed favorable binding affinity towards hAR based on their docking interaction and ideal ADMET properties were chosen for further analysis. The initial docking results and ADMET properties are shown in extended data.

Induced fit docking and the molecular interactions

Molecular interaction studies of the binding cavity of the human Androgen Receptor and molecules are listed in extended data This was compared with the co-crystallized ligand associated with hAR protein R18 and analyzed by Schrodinger-induced fit docking. The ligand 2-hydroxynaringenin demonstrated high affinity for flexible residues within the binding pocket of the Human Androgen receptor protein. The calculated free energy of binding (ΔG) was determined to be -8.59 kcal/mol, indicating a strong binding interaction. While couple of other molecules 8-Prenyldaidzein and 5-Hydroxy-7-acetoxy-8-methoxyflavone also exhibited significant binding with hAR having ΔG = -8.54 kcal/mol and -8.26 kcal/mol, respectively. The highest affinity with a low negative binding energy was observed for 2-hydroxynaringenin, where the ligand formed conventional hydrogen bonds with Leu704, Asn705, Gln711, Met745, Arg752, and Thr877. Leu707, Met780, Leu873, and Phe876 were found to be involved in pi-alkyl and alkyl interactions with the 2-Hydroxynaringenin ligand. The binding energies of 2-Hydroxynaringenin and protein-ligand interactions are displayed in Figure 4 and the binding energies of other molecules are depicted in extended data.

Figure 3. Role of Androgen receptor (Source modified from Ref. 30).

Figure 4. Induced fit docking pose of the ligand (A) 2-Hydroxynaringenin and co-crystallized (B) R18 molecules with hAR (PDB ID: 1E3G) displaying the ribbon shaped 3D protein and ligand interaction, 3D image of binding cavity residues and 2D interaction profile of bidning cavity residues with the respective ligands.

Molecular dynamics simulation studies

Molecular dynamics simulation (MD) investigations were performed to ascertain the convergence and stability of 1E3G-Apo (no ligand hAR protein), 1E3G+R18 (R18 co-crystallized ligand) and 1E3G+2-Hydroxynaringenin complexes. When comparing the root mean square deviation (RMSD) measurements, the 100 ns simulation showed a stable conformation. The Apo protein’s Cα-backbone’s RMSD showed a 3.0 Å divergence ( Figure 5A). While 1E3G+R18 and 1E3G+2-Hydroxynaringenin both showed 2.9 Å, the overall RMSD is shown to be 2.9 Å ( Figure 5A).

The root mean square fluctuations (RMSF) plot of the 1E3G+2-Hydroxynaringenin complex protein revealed notable variations at residues 60–70, 110–120, and 180–185, which may have been caused by the residues’ increased flexibility. The rest of the residues fluctuated less during the course of the 100 ns simulation ( Figure 5B). Radius of gyration (Rg) in this study, 1E3G Cα-backbone bound to Apo protein displayed increment of Rg values indicating lesser compactness while stable Rg was observed from 20.2 to 20.3 Å in 1E3G+R18 ( Figure 5C). The number of hydrogen bonds was significantly different between 1E3G+2-Hydroxynaringenin, throughout the simulation time of 100 ns ( Figure 5D). The average number of hydrogen bonds observed in 1E3G+2-Hydroxynaringenin was two on average in MD simulation studies ( Figure 5D, red color).

Figure 5. MD simulation analysis of 100 ns trajectories of (A) Cα backbone RMSD of 1E3G+2-Hydroxynaringenin (red), RMSD of 1E3GApo (black), and 1E3G+R18 (blue) (B) RMSF of Cα backbone RMSD of 1E3G+2-Hydroxynaringenin (red), RMSD of 1E3GApo (black), and 1E3G+R18 (blue) (C) Radius of gyration (Rg) of Cα backbone of Cα backbone RMSD of 1E3G+2-Hydroxynaringenin (red), RMSD of 1E3GApo (black), and 1E3G+R18 (blue) (D) Formation of hydrogen bonds in 1E3G+2-hydroxynaringenin (red) and R18 (black).

Mechanics generalized born surface area (MM-GBSA) calculations

Utilizing the MD simulation trajectory, the binding free energy along with other contributing energy in form of MM-GBSA were determined for hAR+2-hydroxynaringenin. The results ( Table 4) suggested that the maximum contribution to ΔG_bind in the stability of the simulated complexes were due to ΔG_bindCoulomb, ΔG_bindvdW and ΔG_bindLipo, while, ΔG_bindCovalent and ΔG_bindSolvGB contributed to the instability of the corresponding complexes. The hAR+2-hydroxynaringenin complex showed significantly higher binding free energies. These results supported the potential of 2-hydroxynaringenin, showed the efficiency in binding to the selected protein and the ability to form stable protein-ligand complexes.

Limitations of the study

The current study identified the human Androgen Receptor as a potential candidate drug target to combat TNBC and recognized 2-hydroxynaringenin as a potential lead molecule. The in vitro and in vivo efficacies of 2-hydroxynaringenin require further investigation. Safety, pharmacokinetics, and pharmacodynamics tests need to be performed to further develop hydroxynaringenin for clinical use.