Development of an MCL-1-related prognostic signature and inhibitors screening for glioblastoma

Introduction: The effect of the conventional treatment methods of glioblastoma (GBM) is poor and the prognosis of patients is poor. The expression of MCL-1 in GBM is significantly increased, which shows a high application value in targeted therapy. In this study, we predicted the prognosis of glioblastoma patients, and therefore constructed MCL-1 related prognostic signature (MPS) and the development of MCL-1 small molecule inhibitors. Methods: In this study, RNA-seq and clinical data of 168 GBM samples were obtained from the TCGA website, and immunological analysis, differential gene expression analysis and functional enrichment analysis were performed. Subsequently, MCL-1-associated prognostic signature (MPS) was constructed and validated by LASSO Cox analysis, and a nomogram was constructed to predict the prognosis of patients. Finally, the 17931 small molecules downloaded from the ZINC15 database were screened by LibDock, ADME, TOPKAT and CDOCKER modules and molecular dynamics simulation in Discovery Studio2019 software, and two safer and more effective small molecule inhibitors were finally selected. Results: Immunological analysis showed immunosuppression in the MCL1_H group, and treatment with immune checkpoint inhibitors had a positive effect. Differential expression gene analysis identified 449 differentially expressed genes. Build and validate MPS using LASSO Cox analysis. Use the TSHR HIST3H2A, ARGE OSMR, ARHGEF25 build risk score, proved that low risk group of patients prognosis is better. Univariate and multivariate analysis proved that risk could be used as an independent predictor of patient prognosis. Construct a nomogram to predict the survival probability of patients at 1,2,3 years. Using a series of computer-aided techniques, two more reasonable lead compounds ZINC000013374322 and ZINC000001090002 were virtually selected. These compounds have potential inhibitory effects on MCL-1 and provide a basis for the design and further development of MCL-1 specific small molecule inhibitors. Discussion: This study analyzed the effect of MCL-1 on the prognosis of glioblastoma patients from the perspective of immunology, constructed a new prognostic model to evaluate the survival rate of patients, and further screened 2 MCL-1 small molecule inhibitors, which provides new ideas for the treatment and prognosis of glioblastoma.


Introduction
Glioblastoma (GBM) is a brain tumor originating from glial progenitor cells and is the most common primary malignant tumor of the brain, accounting for 81% of malignant brain tumors (Ostrom et al., 2014;Xu et al., 2020;Li et al., 2023). Among them, glioblastoma has the highest and increasing incidence, but no curative treatment is available (Ostrom et al., 2019). The survival time of most patients is much lower than that of patients with other tumors, and the quality of life is very poor (Finch et al., 2021). At present, the conventional treatment methods for GBM include tumor resection, radiotherapy combined with temozolomide (TMZ) and targeted therapy with bevacizumab, etc., but these treatments have more or less obvious limitations (Allahyarzadeh Khiabani et al., 2023;Boongird et al., 2023;Hotchkiss et al., 2023;Jatyan et al., 2023;You et al., 2023). Therefore, how to predict the prognosis of GBM patients more accurately and intervene the factors affecting the prognosis, formulate more reasonable and effective treatment plans, and develop safer and more effective drugs are the key to treating glioblastoma patients.
Myeloid cell ischemia-1 (MCL-1), as a member of the B-cell lymphoma 2 (BCL-2) protein family, is one of the most frequently amplified genes in all human cancers including glioblastoma (Xiang et al., 2018). MCL-1 has three BH domains (BH1, BH2, and BH3), a C-terminal TM domain, and a large N-terminal region (Li S. et al., 2021). The four binding pockets (P1-P4) of MCL-1 interact with hydrophobic residues H1-H4 of only-BH3 protein, respectively, where the P2 and P3 pockets are the locations of "hot spot" residues for protein-protein interaction in MCL-1. This is different from anti-apoptotic proteins such as BCL2 (P4/ P1 and P2) or BCL-XL (P2 and P4) and facilitates the design of specific MCL-1 inhibitors (Denis et al., 2020). Based on the structure of MCL-1, its role in apoptosis is promoting cell survival by interfering in the cascade of the events that cause trigger cell death and MOMP (Sancho et al., 2021). Numerous previous studies have demonstrated that MCL-1 is extremely important for glioblastoma, which can be used as a key target to inhibit the activity of glioblastoma cells. Downregulation of MCL-1 expression significantly induced apoptosis of tumor cells (Premkumar et al., 2013;Gratas et al., 2014;Jane et al., 2016;Juric et al., 2021). Therefore, there is great potential for the development of MCL-1 specific inhibitors for the treatment of glioblastoma, and it is crucial to create inhibitors that are both more efficient and less poisonous for the treatment of glioblastoma.
At present, the main strategies to design inhibitors against MCL-1 are based on the direct binding of BH3-mimetic to MCL-1, thereby releasing proapoptotic proteins and finally activating apoptosis (Li S. et al., 2021). Clinical studies for MCL-1 inhibitors have started for drugs such S64315, AZD5991, AMG 397, AMG 176, MIM1, etc. (Wei et al., 2020). These inhibitors have shown remarkable efficacy in the treatment of non-Hodgkin's lymphoma, multiple myeloma, acute myeloid leukemia, B-cell lymphoma, and other hematological malignancies (Wei et al., 2020;Liu et al., 2021). However, in the treatment of glioblastoma, most MCL-1 small molecule inhibitors are not applicable due to the existence of blood-brain barrier. Among them, MIM1, as an identified BH3-mimetic, has promising biological and biophysical properties such as low molecular weight, ideal solubility, and stability. MIM1 reduced the viability of glioblastoma cells in a dose and time-dependent manner (Respondek et al., 2018). However, for other types of glioblastoma, the role of MIM1 has not been investigated so far.
The aim of this study is to establish an MCL-1-based prognostic model and to screen safe and effective MCL-1 inhibitors. Firstly, the RNA-seq and clinical data of 168 glioblastoma samples were downloaded from the TCGA database. According to the expression level of MCL-1, the samples were divided into MCL-1_L and MCL-1_H groups, and the enrichment levels of 29 immune signals in the two groups were analyzed. Functional enrichment analysis of these differentially expressed genes was performed. The MCL-1-associated prognostic signature was then constructed using lasso cox regression analysis. Finally, a nomogram prediction model was established to estimate the survival rate of glioblastoma patients. In addition, two small molecule inhibitors of MCL-1 were screened by a series of computeraided techniques. With the development of drug research, natural products are playing an increasingly important role in molecular biology and drug exploration, which provide structural patterns for target compounds of new drugs and are an important source of new drugs. MIM1 was used as a reference drug in this study. NP (natural products) database in the ZINC database was virtually screened to explore potential MCL-1 inhibitors. Secondly, the pharmacological and toxicological characteristics of the compounds were analyzed. Molecular docking was then performed to assess the interaction between the selected compounds and MCL-1. The pharmacophore of the compound was also predicted. Therefore, a more suitable smallmolecule MCL-1 inhibitor is required for the treatment of glioblastoma. Natural products are becoming a more significant part of molecular biology and drug discovery as drug research progresses since they offer structural patterns for target molecules of new medications and are a significant source of such pharmaceuticals. In this study, MIM1 served as the reference medication. Virtual screening was done on the natural products database in the ZINC database to look for probable MCL-1 inhibitors. Secondly, the pharmacological and toxicological properties of the compounds were examined. Then, molecular docking was used to evaluate how well the chosen drugs interacted with MCL-1. Additionally anticipated was the compound's pharmacophore. Finally, using a molecular dynamics simulation, we examined the stability of the binding interaction. The study's findings are summarized in a list of potential MCL-1 small molecule inhibitors and their pharmacological characteristics, which can support and assist the research on MCL-1 inhibitors and give further leads for the creation and advancement of glioblastoma therapy medications.
2 Materials and methods 2.1 Immunogenomic analysis, differential gene expression and functional enrichment analysis RNA-seq and clinical data of 168 GBM samples were downloaded from TCGA (Cancer Genome Atlas database) website, and they were divided into 2 groups according to the Frontiers in Pharmacology frontiersin.org expression level of MCL-1: MCL-1_L (n = 84) and MCL-1_H (n = 84). First, 29 immune signal enrichment levels were quantified in all glioblastoma samples, and single-sample gene set enrichment analysis (ssGSEA) score was used in this analysis. ssGSEA scores were used to analyze the activity or enrichment levels of different immune cell functions in each glioblastoma sample. All glioblastoma samples were then evaluated for the level of immune cell infiltration (immunoscore), stromal content (stromal score), and tumor purity. Finally, the expression of HLA genes and immune checkpoint genes in MCL-1_H and MCL-1_L groups was tested by ANOVA. Rstudio and Wilcoxon Rank Sum And Signed Rank Tests were used to analyze the differential expression between MCL-1_H and MCL-1_L. Using |log2 fold change (FC)|>1 and adjusted p values < 0.05 as the cutoff criterion, all genes were analyzed. "limma" package was used for analysis, and 449 differentially expressed genes were obtained. Then "ggpubr" and "ggthemes" packages were used to visualize the expression levels of all genes, and the differentially expressed genes in the MCL1_H and MCL1_L groups were shown in volcano maps. The Metascape website (https://metascape.org) features gene annotation and visualization. The differentially expressed genes were uploaded to this website, and the gene ontology and signal pathway enrichment of these genes were analyzed.

Construction and validation of MCL-1associated prognostic signature (MPS)
To construct an immune prognostic signature, we randomly divided the TCGA_GBM data set into two groups: training set and verification set. LASSO Cox analysis is a widely used highdimensional predictive regressive method. By selecting the optimal penalty parameter lambda and using 10-fold cross-verification, shrinkage and variable identification can be achieved at the same time to prevent overfitting. To establish immune prognostic characteristics, we put DEGs into LASSO Cox regression and use the "glmnet" R package to proceed to analysis. By weighting Cox regression coefficients to estimate the risk score of each patient, MPS was created. Patients are classified as low-risk and high-risk "survivor" R packages obtained based on the best cut-off value of the risk score. The "survival ROC" R package was used to describe the receiver Frontiers in Pharmacology frontiersin.org 04 operating characteristic (ROC) curve to evaluate its sensitivity and specificity. Calculate the area under the curve (AUC) value according to the ROC curve. At the same time, the prognostic prediction ability of the MPS was further verified in the verification set.

Development of the nomogram
We assessed the independent prognostic power of MPS by univariate and multivariate Cox analyses. And based on the result of Cox analyses, we use the "rms" package to develop an innovative nomogram. Calibration charts for the probabilities of observing and predicting 1-year OS were performed to determine accuracy.

Virtual screening based on structure using Libdock, ADMT and TOPKAT
The Discovery Studio 2019 software, from BIOVIA in San Diego, California, United States, provides researchers with easyto-use tools for protein simulation, modification, and precision medicine (Zhong et al., 2021;Zhong et al., 2022). In addition, we used the ZINC database, a free virtual screening database for commercial chemicals, as a ligand database. 3D molecular formats for 17,931 natural, named, and purchasable chemicals were initially obtained from the ZINC15 database. Both the 1.35 Å crystal structure of human MCL-1 (Protein Data Bank identifier: 6UDV) and the 3D structure of the positive reference medication MIM1 were imported into the LibDock working environment. To identify potential MCL-1 inhibitors, the ligand binding pocket domain of MCL-1 was selected as the binding site. The molecular docking between MIM1 and the treated MCL-1 was found to be successful. This site was therefore used as the active site for docking. All downloaded small molecule files were linked to this active site through the libdock module for preliminary virtual screening. All compounds' docking postures were graded based on their LibDock score. The ADME module of Discovery Studio 2019 was used to calculate the absorption, distribution, metabolism and excretion (ADME) levels of the top 30 compounds, and the TOPKAT module was used to analyze the toxicological characteristics of the compounds. Finally, two molecules were selected as candidates.

Molecular docking and pharmacological analysis
Studies on molecular docking of ligands and proteins were conducted using the CDOCKER module of Discovery Studio 2019. Using the CHARMm force field, CDOCKER is a method that generates very accurate molecular docking statistics Yang et al., 2022). During docking, the CDOCKER algorithm is based on a simulated annealing protocol where the receptor remains rigid and the ligand is allowed to bend and dock with protein residues within the binding site to find the most suitable binding mode . The section within a 13 Å radius of the ligand's geometric center is referred as the binding site spot. Frontiers in Pharmacology frontiersin.org Therefore, according to the interaction energy analysis of CDOCKER, the most suitable compound was selected for the next study. In addition, the best binding pose of the selected compounds to the protein was demonstrated using the Schrodinger software.

Molecular dynamic simulation
The best binding conformation of the ligand-MCL-1 complex obtained from the previous molecular docking step was selected for molecular dynamics simulations. Using an explicit periodic boundary solvated water model, the ligand-receptor complex was contained in an orthogonal box and solvated. Sodium chloride was added to the solution with an ionic strength of 0.145 to mimic the physiological environment. The following simulation protocols were used for the system: 500 steps of steepest descent and conjugate gradient minimization; 5 ps-equilibration simulations in a normal pressure ensemble at a temperature of 300 K (slowly driven from an initial temperature of 50 K); and 50 ps-MD simulation (production module) at NPT (normal pressure and temperature) with a time step of 1 fs. Long-range electrostatics calculations were performed using the particle mesh Ewald (PME) technique, and all hydrogencontaining bonds were fixed using an adaptation of the linear constraint solver technique. The Discovery Studio 2019 analysis trajectory procedure was used to construct a trajectory for root- Frontiers in Pharmacology frontiersin.org mean-square deviation (RMSD), potential energy, and structural parameters using the original complex configuration as a reference.

Immunogenomic analysis between MCL-1_H and MCL-1_L
A total of 168 glioblastoma patients were included in this study. According to the expression level of MCL-1, they were divided into 2 groups: MCL-1_H (n = 84) and MCL-1_L (n = 84). Firstly, the expression levels of 29 groups of immune-related genes representing different immune cell types, functions, and pathways in glioblastoma samples were investigated. According to the ssGSEA score, we found that the enrichment levels of immune cells, functions, and pathways in the MCL-1_H and MCL-1_L groups were not significantly different ( Figure 1A). The ESTIMATE results show that ( Figure 1B), The MCL-1_H group had immune score (Kruskal-Wallis test, p < 0.001), stromal score (Kruskal-Wallis test, p < 0.001), and ESTIMATE score (Kruskal-Wallis test, p < 0.05) were higher than those of MCL-1_L group, but tumor purity (Kruskal-Wallis test, p < 0.001) was higher than that of MCL-1_ H group. These results indicated that MCL-1_H contained more immune cells and stromal cells, and MCL-1_L contained more tumor cells.
The expression of HLA genes and immune checkpoint genes was next analyzed in the two groups. In the analysis of 24 HLA genes (ANOVA test, p < 0.05) ( Figure 1C Figure 1D), gene expression levels were higher in the MCL-1_H group than that in the MCL-1_L group, with significantly higher expression of CD274 and TIMP3 in the MCL-1-H group than in the MCL-1_L group. These results suggested that patients in the MCL-1_H group were immunosuppressed, and the use of immune checkpoint inhibitors, especially inhibitors of CD274 and TIMP3, had a positive effect on the treatment of MCL-1_H patients. In addition, when grouped according to MCL1 expression levels, there was no significant difference in prognosis between the two groups (p = 0.59) ( Figure 1E).

Differentially expressed genes analysis
R software was used to analyze the information of 168 patients to determine the differentially expressed gene (DEG) data set. To FDR <0.05 and log2 fold change (FC) < 1 or more for the standard, there were 449 genes identified, of which 274 genes downregulated, 175 genes upregulated (Figure 2A). Importing 449 differentially expressed genes into the metascape website, it was found that these genes were mainly enriched in these terms: R-HAS-1474244: Extracellular matrix organization, M5884: NABA CODE MATRISOME, GO:0030198: extracellular matrix organization ( Figures 2B-D). Where each node represents a cluster item and is colored by cluster ID and p-value. MPS were constructed in the training group using lasso cox regression analysis ( Figures 3A, B, E) to calculate the risk score for each sample. The survminer package in R software was used to calculate the optimal cut-off value, and the patients in the training group were divided into 2 groups: high-risk group and low-risk group. The results of the Kaplan-Meier analysis indicated that the patients in the low-risk group had a better prognosis ( Figure 3C). Figure 3H shows the distribution of risk scores [risk score = (−0.112721)*TSHR + (−0.016743)* HIST3H2A+ 0.030476*ARGE+ 0.046739*OSMR +0.005866*ARHGEF25] and gene expression in the training group. We analyzed the MPS of the training group by ROC curve, which showed high accuracy in predicting 1-year survival and 3-year survival (AUC of 1-year survival = 0.741, AUC of 3-year survival = 0.775) ( Figure 3F).
We next validated the prognostic value of MPS using the same formula in the validation set. Similarly, all patients were divided into high-risk and low-risk groups. Figure 3D shows that patients in the high-risk group had a lower survival rate. Figure 3I shows the risk score and gene expression distribution of the validation group. Figure 3G demonstrates that MPS has high accuracy and sensitivity in predicting 1-year survival and 3-year survival (AUC of 1-year survival = 0.71, AUC of 3-year survival = 0.81).
3.5 Virtual screening using Libdock, ADME, and TOPKAT of DS 2019 Based on the above results, MCL-1 proved to be a key target for the treatment of glioblastoma and influencing the prognosis of glioblastoma. Therefore, we used MCL-1 as a target for drug screening for the treatment of glioblastoma. The MIM1-MCL-1 complex binding pocket was an essential regulatory region that was chosen as a significant reference site for screening probable MCL-1 inhibitors. Figure 5A shows the 3D structures of MCL-1 (PDB ID: 6UDV) and the MIM1-MCL-1 compound. In accordance with the LibDock score, 4,854 compounds had a higher LibDock score than MIM1 (106.167). The top 30 compounds were chosen for additional investigation based on the LibDock score, and these 30 compounds are described in Table 1.
Then ADME and TOPKAT modules of DS2019 were used to predict the pharmacological and toxicological properties of the top 30 compounds and MIM1. Based on the analysis in Tables 2, 3, ZINC000013374322 and ZINC000001090002 are not inhibitors of CYP2D6 and have a high intestinal absorption level. What's more, these two compounds are weakly bound plasma proteins. More importantly, these two compounds showed no hepatotoxicity and no Ames mutagenicity. Figure 5B shows that these two compounds and MIM1 have similar six-membered and five-membered annular structures, reactive oxygen and nitrogen species atoms, suggesting that they may play similar roles. The compounds ZINC000013374322 and ZINC000001090002, which are expected to be promising candidates, were then investigated.

Ligand binding analysis
ZINC000013374322 and ZINC000001090002 were accurately connected to the binding pocket of MCL-1, and the mechanism of ligand binding was examined by using CDOCKER module. As shown in Table 4, the interaction energies of ZINC000001577210 (−40.7616 kcal/mol) and ZINC000001090002 (−43.3771 kcal/mol) are both much lower than MIM1 (−35.0968 kcal/mol), indicating that they may have higher binding affinity to MCL-1.
By applying DS2019 and Schrodinger software, we thoroughly analyzed the ligand conformation in the MCL-1 binding pocket and the protein-ligand complex interaction ( Figures 6A-C). In these figures, the binding pattern of the two molecules to the MCL-1 binding pocket is visually shown. As shown in Figure 6D, there is a significant overlap between the two molecules and MIM1 in the binding pocket posture. According to Figure 6; Table 5, these two molecules and MIM1 are essentially identical in the way they bind and interact with MCL-1, proving that they have similar inhibitory effects on MCL-1. At the same time, it can be inferred that the two amino acid residues of PHE270 and MET250 play a crucial role in the functional domain of MCL-1.
Through the precise analysis of DS 2019, we showed detailed information on the interaction between ligand and protein, including bond length, bond type, bond atoms, and so on (Tables 5, 6). The results showed that ZINC000013374322 and ZINC000001090002 and MIM1 formed 2, 2, 1 pair of hydrogen bonds with MCL-1, respectively. In addition, ZINC000013374322 and ZINC000001090002, and MIM1 formed 12, 5, 8 pairs of hydrophobic interaction with MCL-1, respectively. Among them, although the hydrophobic interaction formed by ZINC000013374322 and MCL-1 are few, the bond lengths are small. So ZINC000013374322 forms a more stable hydrophobic interaction bonds with MCL-1. In conclusion, these results indicate that ZINC000013374322 and ZINC000001090002 may have better binding affinity to MCL-1 than MIM1, indicating that these two compounds have broad application prospects.

Pharmacophore analysis and molecular dynamic simulation
We performed pharmacophore analysis on these three molecules ( Figure 7A). 9 and 5 hydrogen bond acceptors are  Molecular dynamics simulations were performed in simulated natural environments to evaluate the stability of ZINC000013374322-MCL-1 and ZINC000001090002-MCL-1complexes. As shown in Figures 7B, C, the potential energy and RMSD of each compound became stable over time, and the trajectory of the complex basically reached equilibrium after 15 ps. Molecular dynamics simulation results show that the interaction between these compounds and MCL-1 is beneficial to the stability of the complex. In conclusion, ZINC000013374322-MCL-1 and ZINC000001090002-MCL-1

Discussion
Glioblastoma is a common brain tumor with a high degree of malignancy. Most patients have a poor prognosis and a very short survival time. The current conventional treatment methods for brain tumors have no significant improvement in the survival time and quality of life of glioblastoma. The occurrence and development of most tumors, including glioblastoma, are closely related to cell apoptosis. At present, a large number of studies have proved that the induction of apoptosis through a variety of ways can inhibit the progression of glioblastoma, which is the key way of drug treatment for glioblastoma. MCL-1 is a widely studied and potent anti-apoptotic protein that regulates cells by various mechanisms, including interactions with cell cycle regulators to affect cell division, acting as a molecular switch for double-strand break (DSB) DNA repair, regulation of autophagy and mitophagy through BH3-like protein interactions, etc. (Widden and Placzek, 2021). Downregulating the expression and function of MCL-1 in tumor cells can effectively promote the apoptosis of tumor cells. Among them, the widely studied MCL-1 inhibitor is BH3-mimic, which has good effects on blood system tumors and multiple myeloma and has made great progress in combination with established therapies (Caenepeel et al., 2018). Among them, the widely studied MCL-1 inhibitors are BH3-mimics, including AZD5991, S63845, MIM1, etc. (Tron et al., 2018;Mallick et al., 2019;Moujalled et al., 2019;Wei et al., 2020), which have good effects on hematological tumors and multiple myeloma and have made great progress in combination with established therapies (Caenepeel et al., 2018). In recent years, targeted therapy has been widely used in the treatment of a variety of tumors and has achieved good results. By matching with genetic testing, targeted agents have significantly improved patient outcomes. However, GBM shows strong drug resistance, and the use of targeted drugs is severely limited [21] due to the permeability of the blood-brain barrier that limits drug delivery, low mutation burden, and suppression of the immune microenvironment (Yuan et al., 2022). However, the relatively mature small molecule inhibitors of MCL-1 that have been studied and developed are largely not used in the treatment of glioblastomas, probably due to difficulties in crossing the blood-brain barrier. Therefore, it is essential to study the role of MCL-1 in the prognosis prediction of glioblastoma patients and to develop safer and more effective MCL-1 inhibitors for the treatment of glioblastoma. The present study aimed to establish an MCL-1-based prognostic model and to screen safe and effective MCL-1 inhibitors. Firstly, GBM samples were divided into two groups according to the expression level of MCL-1, and the immune signal enrichment level and differential gene expression were analyzed. Functional enrichment analysis of differential genes showed that these genes were in the extracellular matrix organization, and the expression level of MCL1_H group is significantly higher than that of MCL1_L group. NABA CORE MATRISOME and other aspects were enriched. Next, the prognostic MPS model was constructed and validated, which had high accuracy and sensitivity in predicting 1-year and 3-year survival rates of patients. By LASSO cox regression analysis, 5 genes were found to be independent prognostic factors: TSHR, HIST3H2A, AREG, OSMR, ARHGEF25. Glioblastoma highly expresses TSHR and TSH in the tumor microenvironment promotes its proliferation, invasion and immune evasion, which limits the T cell killing of glioblastoma. Treatment targeting intracranial TSH may reverse the immunosuppressive state of glioblastoma (Wu et al., 2022). miR-516a-5p downregulates the expression of HIST3H2A, thereby reversing the anti-proliferation effect induced by miR-516a-5p in NSCLC cells. miR-516a-5p may inhibit the proliferation of NSCLC cells by targeting HIST3H2A (Ye et al., 2019). AREG, one of the seven ligands that bind and activate EGFR, can promote the differentiation of T cells into Tregs in the tumor microenvironment, and targeting AREG in the tumor microenvironment may inhibit tumor invasion and immunosuppression (Coniglio and Segall, 2021). OSMR is a cell surface receptor that is a key component of EGFRvIII-STAT3 signaling, which forms a feedforward signaling mechanism with these molecules to drive glioblastoma genesis and progression Frontiers in Pharmacology frontiersin.org (Jahani-Asl et al., 2016). ARHGEF25 promotes tumor cell migration, and serum-induced ARHGEF25 activation plays a key role in chemotactic migration by restricting lamellipodia formation to the direction of cell movement and keeping it at the leading edge (Hayashi et al., 2013). Next, we chose the crystal structure of MCL-1 (PDB ID: 6UDV) and used MIM1 as the positive reference drug for the entire study. We first used the LibDock module of DS2019 to analyze the LibDock scores of compounds downloaded from the ZINC15 database. The 30 compounds with the highest LibDock scores were selected to analyze their pharmacological and toxicological properties using ADME and TOPKAT modules. In this step, several pharmacological properties of ZINC000013374322 and ZINC000001090002 were found to be superior to MIM1 and less toxic, so these two compounds were selected for further analysis. Next, the binding modes of the compounds and proteins were precisely analyzed using the CDOCKER module, and the results showed that the complex formed by these two compounds with MCL-1 had lower interaction energy, proving that their binding was more stable. Two amino acid residues, PHE270 and MET250 of MCL-1, were found to interact with two compounds and MIM1, indicating that these two amino acid residues were the key sites for inhibiting MCL-1 protein. This is the latest discovery in our study. At present, the specific functions of these two amino acid residues in MCL-1 have not been mentioned in other existing studies. We consider PHE270 and MET250 as the key amino acid residues in the MCL-1 pocket, and through modification and modification of them, the small molecule inhibitor can be more stably bound to the corresponding domain of MCL-1. In the development of more MCL-1 inhibitors, PHE270 and MET250 can be used as effective binding sites to select more reasonable inhibitors. We then analyzed the pharmacophores of the compounds and found that the compounds had similar characteristic pharmacophores. Finally, molecular dynamics simulations were performed for both compounds. In the simulated natural environment, the complex trajectories reached equilibrium after 30 ps, and the potential energy and RMSD of each complex tend to stabilize with time. Based on these results, ZINC000013374322 and ZINC000001090002 can be further modified to make the ligand binding to MCL-1 more reasonable and stable. Furthermore, few studies have been conducted on these two compounds, ZINC000013374322 (Aurantiamide Acetate) and ZINC000001090002 (Bisdionin B), and studies that have shown no effect on cancer, especially in the treatment of glioblastoma by inhibiting MCL-1 function. However, our study did demonstrate that they can effectively inhibit the MCL-1 function, which provides more prodrugs for the treatment of glioblastoma. Through further modification and improvement, these two compounds show excellent development prospects as MCL-1 small molecule inhibitors.
We have to admit that despite accurate measurements and virtual calculations in this study, there are still some limitations. Since there are many changes in the metabolism and transformation of drugs in vivo, corresponding experiments will be carried out in the future to verify other safety indexes of these two compounds, such as IC50, AB (aerobic biodegradability) and MTD (maximum tolerated dose), etc., to continuously optimize the structure of compounds and develop more reasonable drugs.

Conclusion
This study demonstrated that MCL-1 was a key factor affecting the prognosis of glioblastoma patients, and inhibition of MCL-1 can improve the prognosis of glioblastoma patients. We used lasso cox regression analysis to construct an MCL-1 related prognostic evaluation model and prognostic-related nomogram to predict the survival rate of glioblastoma patients. In addition, we found that TSHR, HIST3H2A, AREG, OSMR, and ARHGEF25 are novel independent factors affecting the diagnosis, treatment, and prognosis of GBM. Based on the role of MCL-1 in glioblastoma development, we used a series of computer-aided techniques to screen safer and more effective MCL-1 small molecule inhibitors from the ZINC15 database. ZINC000013374322 and ZINC000001090002 are safe and ideal drug candidates. In addition, this study also provided 30 candidate drugs and their pharmacological properties, which provided a new idea for the development and study of MCL-1 inhibitors.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.