Identification of Genes Hub Associated with Triple-Negative Breast Cancer and Cannabidiol Analogs Potential Inhibitory Agents: An In-silico Study
GINUMED, Rafael Núñez University Corporation, Cartagena D.T. y C. Colombia.
Pharmacology and Therapeutic group, University of Cartagena, Cartagena D.T. y C. Colombia.
GIBAE Research Group, Faculty of Engineering, University of Cartagena, Cartagena, D.T. y C, Colombia.
*For Correspondence: aalviza@unicartagena.edu.coAbstract
Objective:
Triple-negative breast cancer presents a significant challenge in oncology due to its complex treatment and aggressive nature. This subtype lacks common cancer cell receptors like estrogen, progesterone, and human epidermal growth factor receptor 2 receptors. This study aimed to identify, through bioinformatic analysis, the key genes associated with triple-negative breast cancer. In addition, CBD analogs with potential inhibitory effects on these genes were evaluated through docking and molecular dynamics.
Methods:
Gene expression profiles from the GSE178748 dataset were analyzed, focusing on MDA-MB-231 breast cancer cell lines. Differentially expressed genes were determined through protein-protein interaction networks and subsequently validated. Additionally, the inhibitory effects of cannabidiol analogs on these hub genes were assessed using molecular docking and dynamics.
Results:
Analysis of the hub highlighted RPL7A, NHP2L1, and PSMD11 as significant players in TNBC regulation. Ligand 44409296 showed the best affinity energy with RPL7A, while 166505341 exhibited the highest affinity with NHP2L1 and PSMD11, surpassing CBD. Analyses of RMSD, RMSF, SASA, and Gyration Radius indicated structural stability and interactions of the proteins with ligands over time. MMGBSA calculations showed favorable binding energies for the ligands with the target proteins.
Conclusion:
In conclusion, this study identified key genes, namely RPL7A, NHP2L1, and PSMD11, associated with triple-negative breast cancer and demonstrated promising interactions with cannabidiol analogs, particularly 44409296 and 166505341. These findings suggest potential therapeutic targets and highlight the relevance of further clinical investigations. Additionally, the ligands exhibited favorable ADME properties and low toxicity, underscoring their potential in future drug development for TNBC treatment.
Introduction
Triple-negative breast cancer has emerged as a significant challenge in the field of oncology due to its complex treatment and highly aggressive nature. This cancer subtype is characterized by the absence of typical receptors in cancer cells, such as estrogen and progesterone receptors, and by the overexpression of human epidermal growth factor receptor 2 (HER2) receptors [1,2]. It is estimated that triple-negative breast cancer accounts for approximately 10 to 15% of all diagnosed cases of breast cancer [1]. Furthermore, there is a higher likelihood of this type of cancer developing in women under the age of 40, particularly in those of African descent or those with a mutation in the BRCA1 gene [2].
One of the most studied genetic biomarkers in TNBC is the mutational status of the BRCA1 and BRCA2 genes, which encode proteins involved in DNA repair through homologous recombination. Mutations in these genes are associated with a higher risk of developing TNBC [3]. In addition to genes involved in DNA repair, other genetic biomarkers could be relevant in TNBC, such as those related to cellular signaling pathways, metabolism, angiogenesis, or immune response [4]. Thus, the search for genetic biomarkers and innovative therapies has become a crucial necessity to comprehend the underlying biology of this disease and develop diagnostic strategies for more precise treatments. Genetic and protein biomarkers may play a significant role in the diagnosis, prognosis, treatment selection for TNBC, risk of recurrence, and patient survival [5].
Thus, bioinformatic development has enabled the search and exploration of potential biomarkers from gene repositories identified through their expression profiles using techniques such as Serial Analysis of Gene Expression (SAGE), Single Nucleotide Polymorphism (SNP array), Reverse Transcription Polymerase Chain Reaction (RT-PCR), Massively Parallel Signature Sequencing (MPSS), and high-throughput sequencing (HTS). This has led to the construction of specific interaction pathways and maps, further ensuring new candidates for computational evaluation using tools like molecular docking, pharmacophore studies, and molecular dynamics. These methods have facilitated screenings of multiple potential molecules with inhibitory or adjuvant activity in conventional therapies, minimizing resources and predicting new mechanisms.
Currently, the study of cannabinoid compounds has sparked significant clinical interest in the field of oncology [6]. Most research related to breast cancer and cannabinoids has focused on TNBC (triple-negative breast cancer) models. Among these studies, the phytocannabinoid CBD has received special attention and extensive investigation. It has been demonstrated that CBD can inhibit the proliferation of MDA-MB-231 cells, making it a promising compound in the treatment of this specific form of breast cancer [6,7]. Additionally, the most significant progress has been observed in animal models such as rats inoculated with MDA-MB-231 cells (resistant to DOX), in which cannabidiol in combination with chemotherapeutic agents like doxorubicin was evaluated [8]. Furthermore, it has been demonstrated in other murine models with 6D cell expression that in vivo CBD blocks the tumoral activity of breast cancer cells associated with malignancy by IL-1β [9].
Therefore, this research proposed the search for genes as potential biomarkers associated with triple-negative breast cancer, as well as the study of cannabidiol analogs as promising agents against molecular targets identified through the implementation of computational methodologies for predicting molecular interactions and conformational and structural stability.
Materials and Methods
Access to public data
The GSE178748 expression profile dataset was downloaded from GEO (https://www.ncbi.nlm.nih.gov/gds) and was based on the Agilent-014850 Microarray full human genome microarrays with MDA-MB-231 breast cancer cell lines, by 4x44K G4112F platform (function number version). The MeSH terms used for dataset selection were (Breast cancer, Carcinoma), the study organism was Homo sapiens, the study type was represented by expression profiles by array, and the publication time was limited to one year.
Analysis of gene expression profiles
The differentially expressed genes (DEGs) were analyzed using the interactive online tool GEO2R, which is complemented by the dataset obtained from GEO [10]. The definition of groups was constructed from datasets grouped in MDA-MB231 cancer cell lines regulated by CtBP2 as a control and MDA-MB231 cancer cell lines regulated by CtBP2 silenced with HIPP and P4 inhibitors. The comparative analysis of these three groups was evaluated using a cutoff criterion of P value and adjusted P value <0.05 and fold change (log2 Fold Change) between ≥ 1.5 or ≤ 1.5. The GEO2R analyses obtained were visualized through volcano plots of the DEGs. The intersection analyses between negatively and positively regulated DEGs were illustrated with Venn diagrams, using FunRich tools [11].
Identification and analysis of Hub genes
The DEG profiles identified in the intersection of negatively regulated Vs positively regulated genes were evaluated using the protein-protein interaction network through the online tool STRING (Search Tool for the Retrieval of Interacting Genes). The analysis and visualization of the protein-protein interaction network were performed using Cytoscape software (version 3.9.0) [12], which allowed us to integrate biomolecular interaction networks with DEG data in a unified reading frame. The main genes derived from nodes with the highest scores were obtained by detecting modules with the MCODE algorithm using selection criteria of Cut-off = 0.2; K-core = 2; max Depth = 100. For the identification of Hub genes in the protein-protein interaction network in Cytoscape, the complementary tool cytoHubba was used [13], which employs an algorithm called Maximal Clique Centrality (MCC) [14].
Conversely, were demonstrated the expression patterns of principal hub genes in stages of breast cancer based on Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/index.html), and prediction of survival curve by Kaplan-Meier Plotter (https://kmplot.com/analysis/).
Molecular Docking Simulation
Selection and preparation of ligands and receptors
A search for the cannabidiol structure was performed using the PubChem database [15], which was identified with code CID 644019. In addition, a search for analogs was conducted using a 95% structural similarity based on the Tanimoto Threshold. This was followed by downloading the 3D structure in SDF format of the PubChem database, which revealed 868 molecules. Additionally, 13 molecules were excluded due to incompatibilities or duplication. The structure minimization was carried out using Open Babel algorithms [16], with the MMFF94 force field applied, using conjugate gradient and steepest descent algorithms, including 2000 steps.
Selection and preparation of the receptor
Protein structures associated with the main Hub genes identified in central nodes were selected. To do this, the Uniprot database [17] was used by inserting the corresponding gene name and filtering for reports in humans. The linked codes were: P62424 (RL7A_HUMAN), P55769 (NH2L1_HUMAN), and O00231 (PSD11_HUMAN). The chain related to 60S ribosomal protein L7a (P62424) was obtained by homology from the FASTA sequence using Swiss-MODEL [18]. For the proteins, the crystal structures with the highest resolution were selected, namely chain A of NHP2-like protein 1 (PDB code: 3SIU; 2.63 Å), and chain Q corresponding to 26S proteasome non-ATPase regulatory subunit 11 (PDB code: 5GJQ; 4.35 Å).
Subsequently, each protein system was refined using the online server ModRefiner (https://zhanggroup.org/ModRefiner/). The predicted model was validated using the ERRAT server (Figure 1S). The protein structures were calculated with Kollman atomic charges, and nonpolar hydrogens were merged using AutoDock Tools 1.5.7 [19].
Protein binding pocket prediction
Previously, a predictive search for ligand binding pockets was carried out by PrankWeb (https://prankweb.cz/) for RPL7A and PSMD11, and PockDrug (http://pockdrug.rpbs.univ-paris-diderot.fr/cgi-bin/index.py?page=home) for NHP2L1, using the PDB format files of the generated models and the chains of the downloaded proteins (Figure 2S).
Molecular Docking
The proteins were docked in triplicate using AutoDock-Vina based on PyRx. The grid spacing was set to 1.00 Å, with a center grid space of x = 255.50 Å, y = 235.51, z = 348.78, and offset values of x = 19.03 Å, y = 22.54 Å, z = 20.81 Å for RPL7A; x = 7.258, Å, y = 24.40, z = 9.506 Å and offset values of x = 20.18 Å, y = 16.14 Å, and z = 19.58 Å for NHP2L1; x = 192.99 Å, y = 348.089 Å, z = 293.129 Å and offset values of x = 23.82 Å, y = 21.66 Å, and z = 21.55 Å for PSMD11. Each molecular docking was simulated with a different seed in each run, with an exhaustiveness of 8 and determining 9 conformations based on effectiveness value, free energy binding, and RMSD. The results of the best affinity conformation were shown in pdbqt format and visualized by PyMOL software version 2.3.2 [20], and finally converted to PDB format. The results were expressed as the average of free energy and standard deviation. The 3D ligand-target complexes, interactions, and binding types were visualized using BIOVIA Discovery Studio visualizer version 4.5 [21] and LigPlot+ version 2.2 [22].
Molecular dynamic simulations
The complexes between cannabidiol analogs against RPL7A, NHP2L1, and PSMD11 that showed the highest binding affinity underwent molecular dynamics (MD) simulations. AMBER20 software was employed, running minimization, equilibration, and production processes, along with respective analyses following Alviz-Amador et al. 2019 [23]. Prior, molecules were prepared using the ff14SB force field for proteins and GAFF2 for ligands. The structures were solvated in TIP3P water and subjected to minimization with a combination of 1000 steps of steep descent followed by 1000 steps of conjugate gradient. Subsequently, a 5000-step heating process from 100 to 300 K was carried out with 2 fs intervals. During this heating, a thermostat coupled to constant conditions was employed, and hydrogen bonds were restrained using the SHAKE algorithm with a tolerance of 0.00001. Gradual restraints on peptides and proteins were applied, decreasing from 5 to 0.5 kcal/mol-Å2 in different stages. At each stage, an initial 1000-step steep descent minimization and 500-step conjugate gradient minimization were performed using a 2 fs time step. Then, a 50 ps molecular dynamics simulation was executed at 300 K, maintaining constant pressure and temperature using Berendsen constants of 0.2 ps. The total simulation time is 100 ns. Finally, results were analyzed using AMBER20 CPPTRAJ in terms of root mean square deviation (RMSD), mobility analysis using root mean square fluctuation (RMSF) calculations, solvent-accessible surface area (SASA), and complex compactness through radius of gyration (Rg) measurements. To determine the binding free energy with different ligands, the MM-GBSA method was employed via MMPBSA.py [24], analyzing a 100 ns timeframe for each system.
Prediction of pharmacokinetic properties and toxicity
Pharmacokinetic and drug-likeness prediction for the ligands of higher affinity for each protein was developed by the SwissADME online tool from the Swiss Institute of Bioinformatics (http://www.sib.swiss) [25] and ADMETSAR (http://lmmd.ecust.edu.cn/admetsar2/about/), a web server from East China University of Science and Technology [26]. The canonical SMILES chain in the SwissADME online was included. These parameters as hydrogen bond acceptors and donors, molecular weight, logP, inhibitory cytochrome P450 (CYP450) isoforms, gastrointestinal absorption, binding to P-glycoproteins, blood-brain barrier permeability, plasma protein binding, and Caco2. Likewise, similarity prediction parameters as Lipinki’s rules and bioavailability. Then, in silico toxicity of compounds using the GUSAR-Online server was evaluated; by inserting molecules in format sdf, followed by the predictions of lethal dose 50 (LD50) values for rats oral [26], and finally, showing the classification in rodents based on the OECD Project.
Results
Gene expression analysis
According to the data obtained by GEO2R, a total of 10,251 differentially expressed genes (DEGs) were identified between the cancer cell lines regulated by CtBP2 and those inhibited by HIPP, of which 5,413 genes were upregulated and 4,838 genes were downregulated. Regarding the cancer cell lines regulated by CtBP2 and inhibited by P4, 9,189 DEGs were identified, of which 4,464 were upregulated and 4,725 were downregulated. The volcano plot maps illustrate the significant effect of the difference between the two groups (see Figures 1A and 1B). In the Venn intersection, 81 differentially expressed genes were identified that are expressed in both CtBP2-regulated cell lines and those inhibited with HIPP and P4, as observed in Figure 1C, highlighting their importance in the regulation of triple-negative breast cancer.
Identification and analysis of Hub genes
Protein-protein interaction network
After identifying the most relevant DGEs through the Venn diagram, a protein-protein interaction network was constructed on the STRING platform adjusted to a confidence interval of >0.4. A total of 81 nodes and 39 connections were identified in the network, as shown in Figure 3S. The central nodes generated by the MCODE algorithm identified genes (DHX33, EIF4B, DDX3X), (RPL7A, RRP15, NHP2L1, PSMD11), and (AFF4, CCNT2, and POLR2K) as the most relevant in the interaction network (See Figure 2).
Hub gene identification
For the identification of Hub genes, the protein-protein interaction network module was implemented by running the MCC algorithm of the cytoHubba plugin of the Cytoscape program. The top 5 genes were identified in Table 1; with the genes RPL7A, NHP2L1, and PSMD11 ranking in the top three, which were also identified in the central nodes represented by the MCODE algorithm.
Enrichment of Hub Genes
For the enrichment analysis of the Hub genes, the FunRich software was implemented with a significance threshold of P < 0.05. The obtained results showed that the enriched biological processes pathways were mainly related to protein metabolism and regulation of nucleobase metabolism (33.3%); furthermore, they also included signal translation, cellular communication, energy pathways, and other unknown processes (11.1%). The most relevant molecular functions were linked to transcriptional and translational regulation activities, protease and transferase activities, guanine-nucleotide exchange factor activity, ribosome structural constituent, and RNA binding. The enriched cellular components pathways where the expressed Hub genes can be found were mainly associated with the nucleus (55.6%), nucleolus (44.4%), cytoplasm (44.4%), centrosome (33.3%), ribosomes (22.2%), among others. The most important biological pathways were related to protein (50%), RNA (50%), mRNA (50%) metabolism, gene expression (50%), cap-dependent translation initiation (33.3%), GTP hydrolysis and 60S ribosomal subunit binding (33.3%), 3’-UTR-mediated translation regulation (33.3%), among others. On the other hand, the enrichment pathway of the expression site for the Hub genes shows that these genes are related to their presence in both normal and cancerous tissues. The presence of these genes in breast tissue and breast cancer expression with a percentage of 44.4%; similarly, for other expressions such as ovarian cancer, testicular cancer, cervical cancer, skin cancer, head and neck cancer, lung cancer, and liver cancer, among others. Likewise, an analysis of the Catalogue of Somatic Mutations in Cancer (COSMIC) was performed, which showed the presence of the Hub genes associated with somatic mutations in breast, ovarian, skin, lung, and liver cancer, among others, as shown in Figure 4S.
The expression levels of RPL7A, NHP2L1, and PSMD11 do not vary significantly across tumor stages. However, RPL7A and PSMD11 exhibited prognostic significance according to the Kaplan–Meier plotter, with markedly increased mRNA levels associated with favorable overall survival (OS). In contrast, the levels of NHP2L1 were not correlated with OS. (Figure 5S).
Molecular docking
The molecular docking results indicate that the analog 44409296 had the best affinity energy with a value of -8.10 Kcal/mol for RPL7A, whereas CBD showed values of -6.50, respectively. Similarly, NHP2L1 showed a docking score of -7.37 in interaction with 166505341, while CBD had an energy value of -5.30 Kcal/mol. In the case of the PSMD11 protein, the ligand with the best affinity binding energy was 166505341, with a value of -7.73 Kcal/mol, compared to CBD, which showed values of 6.50 Kcal/mol, respectively (Table 2).
Molecular docking simulation studies reveal that ligand 44409296 exhibited significant binding to RPL7A, whereas 166505341 maintained structural and energetic affinities with NHP2L1 and PSMD11, surpassing those of CBD. Figure 3 illustrates the interactions between the best poses of cannabidiol analogs and proteins associated with the RPL7A, NHP2L1, and PSMD11 genes. For RPL7A, common hydrophobic interactions involving Tyr71, Ile80, and Ile164, along with a single H-bond interaction between I80 and CBD.
Molecular dynamic simulation
In Figure 4, analyses of RMSD, RMSF, SASA, and Gyration Radius of the native RPL7A protein, RPL7A bound to CBD, and RPL7A bound to ligand 4409296, which obtained the lowest interaction energy in molecular docking, Figure 4A, it is observed that the native RPL7A system exhibited the lowest RMSD throughout the trajectory, ranging between 5-10 Å. While RPL7A-CBD exhibits conformational changes exceeding 10 Å between 30 and 60 ns of trajectories, similarly, in RPL7A-Lig, increases in RMSD values exceeding those of the native protein are observed within intervals of 80 to 100 ns. Regarding the RMSF analysis, a similar trend to the RMSD plot is observed in Figure 4B, the native protein exhibited lower fluctuations compared to the RPL7A-CBD and RPL7A-Lig systems, indicating greater mobility in the initial and final portions of the amino acid residues of the protein. In Figure 4C, it is evident that the system with the highest values of Solvent Accessible Surface Area (SASA) is in the range area of 16,000 to 17,000 Ų, and finally, the native RPL7A displays a surface area of 16,000 Ų. Similarly, in Figure 4D, it is shown that RPL7A bound to cannabidiol exhibited the lowest degree of compactness, followed by RPL7A bound to the ligand, and lastly the native RPL7A.
In Figure 5, analyses of RMSD, RMSF, SASA, and Radius of gyration of the native NHP2L1 protein, NHP2L1-CBD, and NHP2L1-166505341, which exhibited the lowest interaction to 18,000 Ų corresponds to RPL7A-CBD. Likewise, the RPL7A-ligand system exhibits surface energy in molecular docking. The RMSD obtained in three systems was between 0.5 to 1.5 Ǻ with major variation for NHP2L1-CBD y NHP2L1-Lig (166505341) in general in all trajectories. NHP2L1-CBD only showed a few fluctuations between residues 10 to 15 aa. The surface area after 100 ns for all three systems ranged between 6500 and 7500 Å2, a consistent aspect across the majority of the simulation time for all systems, with NH2PL1-CBD exhibiting a lower SASA value between 70 to 100 ns of trajectory. As for the gyration radius, the three systems maintained similar behavior on average around 60 to 70 ns, while larger variations were observed for NHP2L1-Lig between 80 to 95 ns.
In Figure 6A, it was observed that the PSMD11-native system exhibited higher variations during the 30 ns trajectory, reaching values above 10 Å. Additionally, there were increments in values between 70 to 100 ns, in comparison to PSMD11-CBD and PSMD11-Lig, which showed lesser changes throughout the simulation time. Regarding RMSF, it was indicated that the native protein established greater fluctuations between residues 120 to 240 amino acids compared to the evaluated ligand systems, as depicted in Figure 6B. For SASA, the highest values were recorded for PSMD11-Native, reaching 26000 Ų after 40 ns, while PSMD11-CBD and PSMD11-Lig confirmed similar behavior as observed in RMSD, as shown in Figure 6C. Figure 6D illustrates the gyration radii behavior, where PSMD11-CBD exhibited overall lower compaction compared to the native protein and Lig at the 100 ns trajectory level.
Free energy calculation
For the stable binding energy analysis, a 100 ns simulation was conducted utilizing 100 complex frames from the trajectory, employing MMGBSA calculations with polar and nonpolar solvation parameters. Table 3 reveals that the average binding energy for RPL7A-44409296 and RPL7A-CBD complex were -40.964 and -4.054, respectively. Similarly, NHP2L1-166505341 and NHP2L1-CBD exhibited binding energies of -30.922 and -23.252. Finally, for the PSMD11 systems, 166505341 and CBD demonstrated binding energies of -35.055 and -28.230, respectively.
Prediction of ADME properties
This computation is anticipated by assessing the likelihood that a molecule exhibits an ideal combination of permeability and bioavailability characteristics. The value of 0.55 indicates adherence to the Lipinski rule, negative BBB permeant, variability in CYP450, and toxicity classified with III and V, respectively Table 4.
| Gene Symbol | Gene Name | Protein Function |
|---|---|---|
| RPL7A | Large ribosomal subunit protein eL8; Ribosomal protein L7a | The ribosome is a large ribonucleoprotein complex responsible for the synthesis of proteins in the cell |
| NHP2L1 | NHP2 Non-Histone Chromosome Protein 2-Like 1 | Involved in pre-mRNA splicing as component of the spliceosome. Binds to the 5'-stem-loop of U4 snRNA and thereby contributes to spliceosome assembly. The protein undergoes a conformational change upon RNA-binding |
| PSMD11 | Proteasome (prosome, macropain) 26S subunit, non-ATPase, 11 | Component of the 26S proteasome, a multiprotein complex involved in the ATP-dependent degradation of ubiquitinated proteins. This complex plays a key role in the maintenance of protein homeostasis by removing misfolded or damaged proteins, which could impair cellular functions, and by removing proteins whose functions are no longer required. Therefore, the proteasome participates in numerous cellular processes, including cell cycle progression, apoptosis, or DNA damage repair. In the complex, PSMD11 is required for proteasome assembly. Plays a key role in increased proteasome activity in embryonic stem cells (ESCs): its high expression in ESCs promotes enhanced assembly of the 26S proteasome, followed by higher proteasome activity |
| RRP15 | Ribosomal RNA processing 15 homolog | Maturation of 5.8S rRNA. Maturation of LSU-rRNA |
| POLR2K | RNA Polymerase II, I and III Subunit K RPABC4 | DNA-dependent RNA polymerase catalyzes the transcription of DNA into RNA using the four ribonucleoside triphosphates as substrates. Common component of RNA polymerases I, II and III which synthesize ribosomal RNA precursors, mRNA precursors of many functional non-coding RNAs and small RNAs such as 5S rRNA and tRNAs, respectively. |
| Ligand | Binding energy (Kcal/mol) | Hydrophobic interaction | Polar interaction |
|---|---|---|---|
| 44409296 | -8,10 ± 0.0 | Val76, Ile80, Asn81, Phe83, Thr84, Val161, Asp160, Asp162, I164 | |
| CBD | -6.50 ± 0.0 | Tyr71, Lys72, Val76, Ile80, Asn81, Thr84, His159, Asp160, Val161, Ile164 | I80 |
| 166505341 | -7.37 ± 0.0 | Leu67, Asn70, Leu71, Tyr80, Glu124, Leu127 | Tyr11, Glu74, Val78 |
| CBD | -5.30 ± 0.0 | Lys9, Leu67, Pro70, Leu71, Tyr80, Glu74, Leu127 | - |
| 166505341 | -7.73 ± 0.05 | Arg132, Phe138, Glu168, Leu171, Ala179, Leu206, Thr209, Ile216 | His78, Asp139, Lys175, Gln213 |
| CBD | -6.50 ± 0.0 | Phe138, Tyr143, His178, Ala179, Thr209, Gln213, Pro317 | Lys175 |
| Protein | Ligand | Average Binding Energy (kcal/mol) |
|---|---|---|
| RPL7A | 44409296 | -40.964 ± 8.67 |
| CBD | -4.054 ± 3.41 | |
| NHP2L1 | 166505341 | -30.922 ± 4.57 |
| CBD | -23.252 ± 7.35 | |
| PSMD11 | 166505341 | -35.055 ± 5,16 |
| CBD | -28.230 ± 3,32 |
| PROPERTIES | 44409296 | 166505341 |
|---|---|---|
| Lipinski #violationsa | 1 | 0 |
| TPSAb | 40.46 | 60.69 |
| GI absorption | High | High |
| BBB permeant | No | No |
| Bioavailability Score | 0.55 | 0.55 |
| Pgp substrate | No | No |
| CYP1A2 inhibitor | No | Yes |
| CYP2D6 inhibitor | Yes | No |
| CYP3A4 inhibitor | Yes | Yes |
| log Kp (cm/s) | -4.02 | -3.91 |
| Caco-2 | 0.8185 | 0.8185 |
| Rat Oral LD50 (mg/kg) | 1409 (III) | 2358 (V) |
Discussion
In the identification of Hub genes linked to TNBC, essential elements such as the RPL7A, NHP2L1, and PSMD11 systems are addressed. Therefore, RPL7A has been established as an RNA-binding protein, known as ribosomal protein, which is an element that can be significantly expressed in multiple cell types. Consequently, altered pathways in ribosomal formation and protein synthesis can substantially induce cancer development. Some research has described the expression of genes like RPL7A in regulating the development of breast cancer cell lines T47D. It has been evidenced that increased ethanol leads to the overexpression of the RLP7A gene and, therefore, greater expressivity in the development of cancer in cell lines, associated with the regulation of translation and oncogenic activation of trk (Modulation of expression of ribosomal protein L7a (rpL7a) by ethanol in human breast cancer cells).
Similarly, Zhang et al. 2022 identified the presence of SH3BGRL as a ribosomal modulator that interferes with the protein translation process, mainly in proteins associated with autophagy. Thus, the interaction of SH3BGRL with a ribosomal subunit RPL7a in different breast cancer cells and 293 T cells was observed in immunofluorescence and co-immunoprecipitation assay [27]. Likewise, studies related to triple-negative breast cancer have identified genes such as RPL27A and RPL15, which were significantly increased in TNBC cancer cells. In vitro studies showed that targeting RPL27A or EIF2 in a cell line like MDA-MB-231 resulted in a significant reduction in cell migration [27]. On the other hand, in different cancer models, such as osteosarcoma, it has been established that ribosomal proteins L7A (RPL7A) showed significantly reduced expression in osteosarcoma tissues compared to regular versus benign tissues. Reduced expression of RPL7A was associated with higher levels of serum alkaline phosphatase and correlated with a critical factor in survival prognosis and potential development of pulmonary metastasis [28].
On the other hand, it has been identified that NHP2L1 participates as a splicing factor in inhibiting the proliferation of TNBC cell lines due to the dysregulation of sister chromatid cohesion (SCC), caused by the retention of intron 1 in sororin. Additionally, the identification of SUN2 is critical in efficient splicing [29]. Similar studies have identified multiple splicing factors, such as NHP2L1 expressed in three breast cancer cell lines (Distinct Splice Variants and Pathway Enrichment in the Cell). Furthermore, it has been observed that NHP2L1 is part of the essential protein set in ribosome biogenesis related to fibrillarin methyltransferase (FBL) of ribosomal RNA (rRNA). These studies have revealed that the overexpression of FBL in MCF7 breast cancer cell lines promotes cell proliferation and resistance to doxorubicin [30,31]. In the same way, rRNA 2’-O-Me enables the alteration of mRNA translation processes encoding oncogenic proteins such as IGF1R or CMYC [32]. On the other hand, studies by Deng et al. 2007 have identified PSMD11 as a gene overexpressed in breast cancer tissue samples, along with the E3 ubiquitin ligase protein (UBE3A). It was determined that the interaction between proteasome proteins and ubiquitin led to an increase in breast cancer [33].
In the link between cannabidiol and TNBC, multiple experimental approaches have been delineated, supporting the molecular potential of the phytocannabinoid against the pathology, as outlined by Surapaneni et al. 2022, they have demonstrated that synthetic cannabidiol (CBD) exerts promising anticancer effects in triple-negative breast cancer (TNBC) cells, especially in three-dimensional (3D) cultures. CBD affects specific genes in MDA-MB-231 cells, reducing fibronectin, vimentin, and integrins in TNBC. It also inhibits autophagy. It has been established that CBD enhances TNBC sensitivity to doxorubicin (DOX), decreasing LOX and integrin-α5, and increasing caspase 9 in MDA-MB-468. These findings suggest that CBD, especially in combination with DOX, could be a promising therapy by inhibiting cell migration, altering genes, and sensitizing cells [34]. Similarly, D’Aloia et al. 2022, found that CBD doses and serum concentrations affect responses in MDA-MB-231 cells. In combination with chemotherapy, CBD protects against cisplatin cytotoxicity under standard conditions, but at threshold concentrations, it induces cell cycle arrest and activates autophagy. Additionally, factors such as IGF-1 and EGF can counteract the antiproliferative effects [35].
Kalevala et al. 2023 established that CBD and THCV significantly enhanced the cytotoxic effectiveness of DOX in 2D and 3D cultures of MDA-MB-231 cells, which were previously resistant to DOX. Transcriptomic and proteomic studies revealed that CBD and THCV, by negatively regulating PD-L1, TGF-β, sp1, NLRP3, and P38-MAPK, while increasing AMPK, induced apoptosis. This action improved the sensitivity of DOX in resistant tumors in BALB/c nude mice. Additionally, the combination of CBD/THCV with DOX reduced H3k4 methylation and H2K5 acetylation, as demonstrated through Western blot assays and RT-PCR [8]. In other studies, CBD was found to decrease cell survival and trigger programmed cell death in both cell lines, as evidenced by viability assays, changes in cell morphology, DNA fragmentation, and apoptosis analysis. CBD-induced apoptosis was linked to the reduction of mTOR, cyclin D1, and the increase in PPARγ protein in nuclei and cytoplasm [36].
Conformationally, the binding of ligand 44409296 to RPL7A revealed the interaction of the 4-methylphenyl aromatic ring with the hydrophobic residues Asn84 and Gly239. These interactions resemble those observed with CBD, although they differ due to the presence of an aliphatic side chain in the 5-pentylbenzene-1,3-diol structural moiety. Additionally, 44409296 exhibited a contrasting pattern with CBD concerning the Phe83 residue, which interacted with the 6-prop-1-en-2-yl portion of the cyclohexene ring. However, while CBD established a polar interaction with Ile80 through the OH group at carbon 3, it did not significantly contribute to the energetic affinity with the RPL7A protein.
In the protein associated with NHP2L1, the analog 166505341 and CBD exhibit similar hydrophobic interactions with residues L67, P70, L71, and L127. For the NHP2L1 systems, it was observed that the analog 166505341 generated a larger coverage area around the protein pocket and a greater number of significant interactions. These interactions were delineated by the potential presence of hydrogen bonding interactions with the residues Tyr11 and Glu74, which formed bonds with the OH groups belonging to the benzene-1,3-diol structural moiety. Additionally, the residue Val78 was identified, establishing a polar interaction with the 3-hydroxyphenyl group of the CBD analog. In contrast to the CBD molecule, it did not exhibit polar interactions but instead displayed hydrophobic interactions, most of which coincided with the higher affinity analog. However, this compound also demonstrated distinct hydrophobic interactions with NHP2L1, involving Cys73 and Asn77 residues.
On the other hand, for the PSMD11 protein, it was identified that the highest affinity analog166505341 maintained a high consistency in common interactions with the associated protein, predominating in apolar interactions such as F138, L136, L171, L172, and L206. The CBD analog exhibited variations in the residues interacting with the protein chain. Notably, for the PSMD11 protein, it was observed that the structural 166505341 displayed a higher number of polar bonds due to the identification of four interactions through hydrogen bridges formed by the residues Gln213 and Asp139 with the diol moiety of the analog. Additionally, other hydrogen bond interactions were formed by Lys175 and His178 with the same OH group present in the 3-hydroxyphenyl structural portion. Similarly, the CBD molecule indicated the interaction of Lys175 with the OH group located as a substituent at C3, forming the diol of CBD. It is also noteworthy that compound 166505341, in comparison to CBD, interacted via the dimethylocta-2,6-dienyl side chain with distinct hydrophobic residues such as Glu168 and Leu206, which appear to be more deeply involved within the PSMD11 pocket.
Furthermore, molecular dynamics for the RPL7A-native system have indicated that its structure remained stable, as it does not have any ligand that would alter its conformation. On the other hand, it can be seen that the RPL7A-ligand of interest system had a slightly higher average RMSD compared to the native system, although quite similar in the first 50 ns. However, towards the end of the trajectory, it reaches values between 10-15 Å, which may be attributed to the ligand of interest having a higher affinity or specificity for the RPL7A receptor, leading to significant changes in its conformation or dynamics. Additionally, it is observed that the RPL7A-Cannabidiol system exhibits a more fluctuating curve compared to the other systems, indicating greater variability in its RMSD. The above could be due to cannabidiol having a lower affinity or specificity for the RPL7A receptor, resulting in weak or transient interactions with the receptor [37].
Overall, RPL7A bound to CBD exhibits the highest residue mobility, followed by RPL7A bound to the ligand, and lastly, the native protein. It can be observed for all systems that there is greater mobility in the residues at the ends, specifically in the region spanning from amino acid 1 to 40 (VVNPLFEKRPKN FGIGQDIQPKRDLTRF VKWPRYIRLQRQ) and from 220 to 240 (SVARIAKLEKAKAKELATKLG). This behavior suggests that when RPL7A is bound to cannabidiol, it demonstrates greater expansibility compared to its binding with the ligand or in its native state. The native state is characterized by a more contracted or less solvent-accessible form. This confirms the trend of higher conformational changes resulting from the binding between RPL7A and cannabidiol. The Rg result has demonstrated that RPL7A bound to CBD displays the highest instability, fluctuation, expansibility, and lower rigidity. This underscores that CBD promotes conformational changes in the RPL7A protein, which could be the atomic-level molecular mechanisms explaining its effect on this protein.
The variations in the RMSD analysis for all three systems remained within the range of 0.5 to 1.5 Å throughout the trajectory. Similarly, the analyses of RMSF, SASA, and Rg did not reveal significant differences in residue fluctuations, expansibility, and compactness of the native NHP2L1 protein, NHP2L1-CBD, and NHP2L1-Lig, respectively. The previous suggests that at the atomic level, neither cannabidiol nor the studied ligand affects this protein system.
In PMSMD11, it is evident that the native protein undergoes substantial conformational changes during the trajectory, indicating relatively low stability, with RMSD values fluctuating between approximately 6 to 11 Å. However, when PSDM11 binds to ligand 166505341 or CBD, the conformational stability is improved, as reflected by RMSD values within the narrower range of 5 to 8 Å (ligand-bound) and 4 to 6 Å (CBD-bound). Equally, the PSDM11-Native exhibits high residue fluctuations, particularly in the region between residues 120 and 240, reaching approximately 8 Å. Conversely, both ligand 166505341-bound and CBD-bound PSDM11 proteins demonstrate reduced mobility in these residues, suggesting enhanced stability upon binding. The SASA analysis indicates that the native PSDM11 protein displays greater solvent-accessible surface area fluctuations between 50 and 90 ns than the proteins bound to CBD and ligand 166505341. These results suggest that drug binding leads to a reduction in solvent accessibility. Lastly, shows that PSDM11-native proteins exhibit substantial variation in compactness throughout the trajectory. However, this variability diminishes when the PSDM11 protein binds to ligand 166505341 and CBD. These findings suggest that CBD and ligand 166505341 have a stabilizing effect on the PSDM11 protein, resulting in reduced conformational fluctuations and increased compactness, potentially influencing its biological function and interactions at the atomic level.
In the evaluated systems, it was identified that the average binding energy was stronger in ligands with higher affinity for each system. Furthermore, in the RPL7A-44409296 system, in comparison to CBD, a tenfold difference in binding was observed. In contrast, systems such as NHP2L1-166505341 and PSMD11-166505341 exhibited affinity energies only approximately 1.3 times compared to the reference analog. This disparity could be attributed to the stable orientation of the ligands and the potential contribution of interactions with amino acid residues in the active pocket of each system. The binding energy obtained in these results validates the affinities observed in the simulations through molecular docking.
In conclusion, the study identified key hub genes, namely RPL7A, NHP2L1, and PSMD11, associated with triple-negative breast cancer. These genes play essential roles in ribosomal formation, protein synthesis, splicing factors, and proteasome function, highlighting their significance in TNBC development and progression. Cannabidiol (CBD) and its analogs, specifically ligands 44409296 and 166505341, showed promising interactions with RPL7A, NHP2L1, and PSMD11. These interactions suggest a potential therapeutic avenue for TNBC treatment, especially when combined with standard chemotherapy drugs like doxorubicin. Molecular dynamics simulations revealed that the binding of ligands to RPL7A, NHP2L1, and PSMD11 induced conformational changes in the proteins. These changes may influence the biological function and interactions of these proteins at the atomic level. Ligand 44409296 demonstrated a notably higher affinity for RPL7A compared to CBD, indicating its potential as a more effective therapeutic agent. Similarly, ligand 166505341 exhibited a strong affinity for both NHP2L1 and PSMD11, suggesting its potential as a promising candidate for TNBC treatment. The stability and binding energies observed in the molecular dynamics simulations validated the results obtained through molecular docking, reinforcing the potential of these ligands as therapeutic options for TNBC. Overall, this study provides valuable insights into the molecular mechanisms underlying TNBC and highlights the potential of CBD analogs, particularly ligands 44409296 and 166505341, as promising candidates for further development in TNBC treatment. Further experimental validations and clinical investigations are warranted to confirm the therapeutic efficacy of these compounds in TNBC therapy.
Acknowledgements
The authors especially thank the University of Cartagena and Rafael Núñez University Corporation
Funding Information
The authors did not receive funding support.
Conflict of interest
The authors report no conflicts of interest.
Data Availability
10.6084/m9.figshare.24587295