Pharmacophore-Model-Based Virtual-Screening Approaches Identified Novel Natural Molecular Candidates for Treating Human Neuroblastoma

The mortality of cancer patients with neuroblastoma is increasing due to the limited availability of specific treatment options. Few drug candidates for combating neuroblastoma have been developed, and identifying novel therapeutic candidates against the disease is an urgent issue. It has been found that muc-N protein is amplified in one-third of human neuroblastomas and expressed as an attractive drug target against the disease. The myc-N protein interferes with the bromodomain and extraterminal (BET) family proteins. Pharmacologically inhibition of the protein potently depletes MYCN in neuroblastoma cells. BET inhibitors target MYCN transcription and show therapeutic efficacy against neuroblastoma. Therefore, the study aimed to identify potential inhibitors against the BET family protein, specifically Brd4 (brodamine-containing protein 4), to hinder the activity of neuroblastoma cells. To identify effective molecular candidates against the disease, a structure-based pharmacophore model was created for the binding site of the Brd4 protein. The pharmacophore model generated from the protein Brd4 was validated to screen potential natural active compounds. The compounds identified through the pharmacophore-model-based virtual-screening process were further screened through molecular docking, ADME (absorption, distribution, metabolism, and excretion), toxicity, and molecular dynamics (MD) simulation approach. The pharmacophore-model-based screening process initially identified 136 compounds, further evaluated based on molecular docking, ADME analysis, and toxicity approaches, identifying four compounds with good binding affinity and lower side effects. The stability of the selected compounds was also confirmed by dynamic simulation and molecular mechanics with generalized Born and surface area solvation (MM-GBSA) methods. Finally, the study identified four natural lead compounds, ZINC2509501, ZINC2566088, ZINC1615112, and ZINC4104882, that will potentially inhibit the activity of the desired protein and help to fight against neuroblastoma and related diseases. However, further evaluations through in vitro and in vivo assays are suggested to identify their efficacy against the desired protein and disease.


Introduction
Human cells carry cell cycle-based active enhancers, promoter regions, and a number of super-enhancers. Through analysis of myeloma cell lines, 10,000 active transcription sites, about 8000 enhancers, and 308 super-enhancers have been identified [1]. The bromodomain and extraterminal domain (BET) protein, which maintains gene expression, is one of the epigenetic readers [2]. BET proteins operate as substrates to draw several other proteins to the promoter region and at the enhancer of active genes, particularly super-enhancer genes, and promote transcription by recruiting kinetics as well as interdependency [1,3]. BET series of proteins (Brd2, Brd3, Brd4, and Brdt) are chromatin reader proteins well known for regulating the transcription of genes that participate in oncogenesis and cell proliferation [4]. Different hematologic cancer models show that Bet protein antagonists can produce an anti-proliferative effect partly because they repress the MYCN oncogene by inhibiting the Myc-driven pathway [5]. The MYC family of oncogenes includes n-myc, c-myc (MYC), and l-myc and usually expresses in the central nervous system, usually in the forebrain and hindbrain [6,7]. The structure of the myc-N protein comprises 464 amino acid (AA) residues along with two domains, namely N and C. The N-terminus domain confers the ability to activate the transcription of genes associated with the production of the myc-N family protein, whereas the C-terminal conserved domain is responsible for DNA binding [6]. Overexpression of the MYCN gene causes cancer, specifically neuroblastoma, T-cell leukemia, glioblastoma, and breast adenocarcinoma [8,9]. Through the recruitment and maintenance of the pTEFb complex at gene promoters during mitosis, BRD4 controls various genes vital for cell function. MYCN suppression was varied in a neuroblastoma cell line that expresses myc-N, so the regulation of MYC and the cellular response to BET inhibition differ between hematologic and solid tumor types, such as neuroblastoma [10,11]. BET antagonist's able to block myc-N over expression might occasionally be linked to their anti-cancer effect, although in some cases, these effects are controlled by a distinct group of cancer-related genes [12]. To date, no specific therapeutics have been developed targeting the protein that can be used to fight against the diseases.
Neuroblastoma is the most predominant cancer among children, and the two forms of the disease include localized and metastatic [13]. The survival rate is lower for metastatic tumors, typically less than 40%. Early-stage detection and diagnosis of the tumor play a vital role in the overall survival rate of cancer. If diagnosed early, at less than one year, the tumors can often be easily cured. However, diagnoses after one year have an increased possibility of metastases, with a low survival rate, and require rigorous chemotherapy and other treatment option [14,15]. In most cases, patients resistant to therapies have been observed to have Brd4 amplification along with MYCN [16]. Targeting Brd4 protein is essential in treating several cancers, and few chemical synthesis drugs have been developed that can reduce its overexpression by inhibiting transcription and direct target-based approaches [9]. Retinoic acid works at the mRNA level to downregulate protein expression via cell cycle arrest in neuroblastoma treatment [17]. It has been observed in-vivo studies that several inhibitors can target BET protein in case of acute myeloid leukemia (AML), myeloma, and lymphoma. Proto-oncogenes like Myc, as well as b-cell lymphoma 2 (BCL 2), can suppress transcriptional blood cancer cells and neuroblastoma through the exposure to the BET inhibitors [18]. BET bromodomain inhibitors lead to the differentiation of the T cells, which leads to producing therapeutic effects in autoimmune disease and also suppresses the overexpression of the MYC through transcriptional effects in treating leukemia cells [5,19].
Nowadays, using chemically synthesized compounds as treatment for neuroblastoma has shown long-term side effects, including growth impairment, neurological disorders, the development of leukemia, and hormonal imbalances [20]. Using natural compounds instead of chemically synthesized compounds may help to reduce the side effects arising from the treatment of previously developed drugs.
Computer-aided drug design (CADD), also known as in silico drug design, is currently a popular approach to designing drugs against a specific target. Modern CADD approaches, including pharmacophore modeling, molecular docking, ADMET (absorption, distribution, metabolism, excretion, and toxicity), molecular dynamics (MD) simulation, and MM-GBSA (Molecular mechanics with generalized Born and surface area solvation), can generate hits to lead compounds more rapidly compared to the conventional drug design process from the ZINC database through virtual screening as it contains several libraries of well-characterized [21][22][23]. For example, molecular docking can initially help determine the binding affinity of the ligands with the target protein that shows a desired biological response measure through in vitro and in vivo methods [24]. In conventional drug development, the ADMET parameters can be evaluated by obtaining blood samples from the clinical trial stage or in vivo [25,26]. However, CADD approaches utilize different techniques that can evaluate the ADMET parameters more easily and efficiently. The techniques can determine a broad range of organ and other toxicity properties, indicating the possibility of contradictory effects [27]. Where the MD simulation can determine the stability of the compound to the target protein in different time scales, MM-GBSA can subsequently determine the binding free energy of a drug-like compound to the target protein [16]. CADD approaches facilities to identify drug-like compounds more easily. Therefore, this study aimed to identify potential natural compounds against neuroblastoma disease using pharmacophore modeling, virtual screening, molecular docking, ADMET, MD simulation, and MM-GBSA approaches. The natural lead compounds identified through the modern CADD approaches will be able to reduce the overexpression of the BET protein and, subsequently, the overexpression of MYCN proto-oncogene, resulting in hindering the activity of the neuroblastoma and related disease.

Development of Pharmacophore Model
To identify novel scaffolds as a lead compound against Brd4 protein, structure-based pharmacophore models were generated based on the protein data bank (PDB) ID: 4BJX in complex with the ligand PDB ID: 73B. The inhibitory activity of the ligand 73B with the Brd4 protein has confirmed previously through experimentally inhibitory concentration (IC 50 ): 21 nM, and validated by X-ray diffraction, resolution: 1.59 Å, R-value free: 0.185, indicating that this series of inhibitors could bind the active site of the protein subsequently inhibition of their biological activity. Therefore, the Brd4 protein (4BJX) structure was used to generate the pharmacophore model in this study, and Ligand Scout 4.4 Advance molecular design software was used to create the pharmacophore model ( Figure 1A). Analysis of the complex 3D structure of the protein found six hydrophobic contacts, two hydrophilic interactions, one negative ionizable bond, and fifteen exclusion volumes presented in Figure 1B.

Validation of Pharmacophore Model
Validation is a crucial parameter for an authentic pharmacophore model, as well as determining the model's quality. To validate the pharmacophore model initially, 36 active antagonists of the selected protein were identified from a literature search and the online ChEMBL database (Table S1). Subsequently, the active compounds were submitted to the DUD-E decoy database, and corresponding decoy compounds were retrieved.
The active compounds and corresponding decoy sets were used to validate the pharmacophore model. Initially, the pharmacophore model-based screening process generated a ROC curve based on the discrimination ability of the model from active to inactive. The obtained sensitivity curve (ROC curve) indicated both specificity and sensitivity. The quality of the curve was shown by the area under the curve (AUC) and the enrichment factor (EF) value. AUC values ranging from 0 to 0.5 suggest a likelihood of discrimination; acceptable ranges were 0.51 to 0.7 and 0.71 to 0.8, respectively, indicating good and excellent performance ( Figure 2). Simultaneously, the ER value informs us about the presence of activity from a certain model. Our ROC curve showed good results, with an AUC of 1.0, while ER values of 11.4 to 13.1 showed excellent outcomes with a value greater than 0.9. Analysis of the ROC curve found a total of 36 true positives (TF) and 3 false positives

Validation of Pharmacophore Model
Validation is a crucial parameter for an authentic pharmacophore model, as well as determining the model's quality. To validate the pharmacophore model initially, 36 active antagonists of the selected protein were identified from a literature search and the online ChEMBL database (Table S1). Subsequently, the active compounds were submitted to the DUD-E decoy database, and corresponding decoy compounds were retrieved.
The active compounds and corresponding decoy sets were used to validate the pharmacophore model. Initially, the pharmacophore model-based screening process generated a ROC curve based on the discrimination ability of the model from active to inactive. The obtained sensitivity curve (ROC curve) indicated both specificity and sensitivity. The quality of the curve was shown by the area under the curve (AUC) and the enrichment factor (EF) value. AUC values ranging from 0 to 0.5 suggest a likelihood of discrimination; acceptable ranges were 0.51 to 0.7 and 0.71 to 0.8, respectively, indicating good and excellent performance ( Figure 2). Simultaneously, the ER value informs us about the presence of activity from a certain model. Our ROC curve showed good results, with an AUC of 1.0, while ER values of 11.4 to 13.1 showed excellent outcomes with a value greater than 0.9. Analysis of the ROC curve found a total of 36 true positives (TF) and 3 false positives (FP) from the 472 compounds ( Figure 2). The calculated GH and EF score has also been mentioned in Table S2.

Validation of Pharmacophore Model
Validation is a crucial parameter for an authentic pharmacophore model, as well as determining the model's quality. To validate the pharmacophore model initially, 36 active antagonists of the selected protein were identified from a literature search and the online ChEMBL database (Table S1). Subsequently, the active compounds were submitted to the DUD-E decoy database, and corresponding decoy compounds were retrieved.
The active compounds and corresponding decoy sets were used to validate the pharmacophore model. Initially, the pharmacophore model-based screening process generated a ROC curve based on the discrimination ability of the model from active to inactive. The obtained sensitivity curve (ROC curve) indicated both specificity and sensitivity. The quality of the curve was shown by the area under the curve (AUC) and the enrichment factor (EF) value. AUC values ranging from 0 to 0.5 suggest a likelihood of discrimination; acceptable ranges were 0.51 to 0.7 and 0.71 to 0.8, respectively, indicating good and excellent performance ( Figure 2). Simultaneously, the ER value informs us about the presence of activity from a certain model. Our ROC curve showed good results, with an AUC of 1.0, while ER values of 11.4 to 13.1 showed excellent outcomes with a value greater than 0.9. Analysis of the ROC curve found a total of 36 true positives (TF) and 3 false positives (FP) from the 472 compounds ( Figure 2). The calculated GH and EF score has also been mentioned in Table S2. del eter for an authentic pharmacophore model, as well as o validate the pharmacophore model initially, 36 active were identified from a literature search and the online sequently, the active compounds were submitted to the sponding decoy compounds were retrieved. rresponding decoy sets were used to validate the pharharmacophore model-based screening process genercrimination ability of the model from active to inactive. C curve) indicated both specificity and sensitivity. The y the area under the curve (AUC) and the enrichment ing from 0 to 0.5 suggest a likelihood of discrimination; and 0.71 to 0.8, respectively, indicating good and excelltaneously, the ER value informs us about the presence ur ROC curve showed good results, with an AUC of showed excellent outcomes with a value greater than nd a total of 36 true positives (TF) and 3 false positives ure 2). The calculated GH and EF score has also been ) indicating H-bond acceptors, and one red star (

Validation of Pharmacophore Model
Validation is a crucial parameter for an authentic pharmacophore model, as w determining the model's quality. To validate the pharmacophore model initially, 36 antagonists of the selected protein were identified from a literature search and the ChEMBL database (Table S1). Subsequently, the active compounds were submitted DUD-E decoy database, and corresponding decoy compounds were retrieved.
The active compounds and corresponding decoy sets were used to validate the macophore model. Initially, the pharmacophore model-based screening process ated a ROC curve based on the discrimination ability of the model from active to in The obtained sensitivity curve (ROC curve) indicated both specificity and sensitivi quality of the curve was shown by the area under the curve (AUC) and the enric factor (EF) value. AUC values ranging from 0 to 0.5 suggest a likelihood of discrimi acceptable ranges were 0.51 to 0.7 and 0.71 to 0.8, respectively, indicating good and lent performance ( Figure 2). Simultaneously, the ER value informs us about the pr of activity from a certain model. Our ROC curve showed good results, with an A 1.0, while ER values of 11.4 to 13.1 showed excellent outcomes with a value greate 0.9. Analysis of the ROC curve found a total of 36 true positives (TF) and 3 false po (FP) from the 472 compounds ( Figure 2). The calculated GH and EF score has als mentioned in Table S2.
Curr. Issues Mol. Biol. 2022, 2, FOR PEER REVIEW 5 Figure 2. Representing a receiver operating characteristic (ROC) curve obtained during the pharmacophore model validation process. The total number of active and decoy sets was determined using the ChEMBL and DUD-E decoy databases, respectively.

Dataset Identification for Pharmacophore Screening
Dataset generation is essential for identifying lead compounds. The ZINC is the most commercially available database, containing 230 million purchasable compounds and 1.37 million compounds (including natural and synthetic compounds) ready-to-dock for achieving binding affinity and providing the 3D structures for clinical development [28]. We utilized the natural product database and the natural derivatives database, which con- Figure 2. Representing a receiver operating characteristic (ROC) curve obtained during the pharmacophore model validation process. The total number of active and decoy sets was determined using the ChEMBL and DUD-E decoy databases, respectively.

Dataset Identification for Pharmacophore Screening
Dataset generation is essential for identifying lead compounds. The ZINC is the most commercially available database, containing 230 million purchasable compounds and 1.37 million compounds (including natural and synthetic compounds) ready-to-dock for achieving binding affinity and providing the 3D structures for clinical development [28]. We utilized the natural product database and the natural derivatives database, which contain millions of currently marketed drugs, natural products, and semisynthetic medical products. The database was initially filtered based on molecular weight <500 kDa and retrieved for the further screening process.

Pharmacophore Model-Based Virtual Screening
The validated pharmacophore model obtained from the LigandScout 4.4.8 advanced tools was used to screen the ZINC natural compounds library. Initially, the validated model was submitted in the ZINCPharmer pharmacophore search tools for screening the purchasable natural subset of the ZINC database. During the screening process, the molecular weight limit set was 500 kDA, and the RMSD was ≤1. One hydrophobic feature with the inclusion-exclusion volume was omitted before screening the ZINC natural compound library. The pharmacophore model-based virtual screening process identified 136 hits with fit scores ranging between 85.2 to 96.7 that were saved for further screening.

Binding-Site Detection and Ligand-Protein Interaction Prediction
The desired Brd4 protein (PDB ID: 4BJX) had two associated ligands and separate attachment sites for interacting with the target ligand. The bond formation with the various amino acid (AA) residues of the protein was demonstrated by analyzing the protein-ligand complex structure using Discovery Studio software and a CASTp server ( Figure 3A). The numbers of active sites interacting with the ligand are shown ( Figure 3B), and the number of active sites determined by CASTp is provided in Figure S1A. The program was used to count the number of hydrophobic and hydrogen-bond-acceptor contacts between amino acids of the protein. The hydrogen-bond-acceptor contacts are represented by red arrows ( Figure S1C). The Ramachandran plot was used to determine the beta-sheet, right-handed helix, and left-handed helix in the protein structure ( Figure S1B).
Curr. Issues Mol. Biol. 2022, 2, FOR PEER REVIEW 6 with the inclusion-exclusion volume was omitted before screening the ZINC natural compound library. The pharmacophore model-based virtual screening process identified 136 hits with fit scores ranging between 85.2 to 96.7 that were saved for further screening.

Binding-Site Detection and Ligand-Protein Interaction Prediction
The desired Brd4 protein (PDB ID: 4BJX) had two associated ligands and separate attachment sites for interacting with the target ligand. The bond formation with the various amino acid (AA) residues of the protein was demonstrated by analyzing the proteinligand complex structure using Discovery Studio software and a CASTp server ( Figure  3A). The numbers of active sites interacting with the ligand are shown ( Figure 3B), and the number of active sites determined by CASTp is provided in Figure S1A. The program was used to count the number of hydrophobic and hydrogen-bond-acceptor contacts between amino acids of the protein. The hydrogen-bond-acceptor contacts are represented by red arrows ( Figure S1C). The Ramachandran plot was used to determine the beta-sheet, right-handed helix, and left-handed helix in the protein structure ( Figure S1B).

Molecular Docking
Docking is a technique used in drug development to assess the binding affinity of a protein and its ligand [29]. The method aids in predicting the behavior of proteins (enzymes) and their interactions with small molecules (ligands) to the binding sites of target

Molecular Docking
Docking is a technique used in drug development to assess the binding affinity of a protein and its ligand [29]. The method aids in predicting the behavior of proteins (enzymes) and their interactions with small molecules (ligands) to the binding sites of target proteins. Therefore, this study used a molecular docking simulation to determine the ligand-receptor binding process. The 136 hits screened through the structure-based virtual screen process were docked to the desired protein. Initially, the protein was prepared using the Biovia Discovery Studio tools by adding a force field, namely CHARMm. Later the binding site residues identified from the complex structure analysis process were used to make a receptor grid, where the grid dimension was X = 30.1404, Y = 25.00, and Z = 25.4020. The hits were docked with the intended protein using AutoDock Vina, and the four compounds with the highest binding energy (ZINC4104882, ZINC2509501, ZINC2566088, and ZINC1615112) were chosen for further analysis ( Table 1). The four compounds chosen in the study were also validated through the re-docking process (Table S3). For docking validation, the same grid box position was used, and the docking was visualized, as shown in Figure S2.

Pharmacophore Features Analysis
Prior to the development of medicine, lead development and screening were important aspects of the pharmaceutical company [14]. For the virtual scanning of the database library, many methods for identifying pharmacophore features are required. To optimize a biologically active compound, compounds with similar properties should be identified. Using the rule of five, we analyzed the drug-likeness and non-drug properties of the top four higher-binding-energy compounds, ZINC4104882, ZINC2509501, ZINC2566088, and ZINC1615112. The pharmacophore features generated from the selected compounds had better or similar properties with respect to the attached antagonist For ZINC1615112 (binding affinity, −7.2 kcal/mol), there were five pi-alkyl bonds (PRO82, VAL87, ILE146, LEU94, and VAL87), two carbon-hydrogen bonds (PRO86 and PRO82), two pi-pi T-shaped bonds (TRP81), and five pi-alkyl bonds (PHE83, THR97, THR139, LEU92, and ILE146).

Pharmacophore Features Analysis
Prior to the development of medicine, lead development and screening were important aspects of the pharmaceutical company [14]. For the virtual scanning of the database library, many methods for identifying pharmacophore features are required. To optimize a biologically active compound, compounds with similar properties should be identified. Using the rule of five, we analyzed the drug-likeness and non-drug properties of the top four higher-binding-energy compounds, ZINC4104882, ZINC2509501, ZINC2566088, and ZINC1615112. The pharmacophore features generated from the selected compounds had better or similar properties with respect to the attached antagonist (ZINC95504909). Therefore, the selected compounds should have shown better efficacy against the Brd4 protein and, subsequently, the MYCN gene. The features obtained for the four selected compounds are shown in Figure 6.

MM/GBSA Analysis
The MM-GBSA techniques were employed to determine the binding free energy of the selected four compounds in complex with the targeted protein. The desired ligand confirmation obtained through the docking process was used to estimate the binding free energy [22,30]. The selected four compounds ZINC2509501, ZINC2566088, ZINC1615112, and ZINC4104882 in complex with the targeted protein show ∆G binding (∆GB) energies value −19.91 kcal/mol, −16.69 kcal/mol, −17.98 kcal/mol, −15.34 kcal/mol (Table 2). Additionally, analysis of the physical and chemical properties found an invaluable achievement of ∆G Bind Covalent (∆GBCo), ∆G Bind Hbond (∆GBH), ∆G Bind Lipo (∆GBL), ∆G Bind Packing (∆GBP), ∆G Bind Solv GB (∆GBS), and ∆G Bind vdW(∆GBV) energies for the selected four compounds listed in Table 2. Therefore, it can be said that the compounds would be able to make a prolog binding affinity with the desired Brd4 protein.
(ZINC95504909). Therefore, the selected compounds should have shown better efficacy against the Brd4 protein and, subsequently, the MYCN gene. The features obtained for the four selected compounds are shown in Figure 6.

MM/GBSA Analysis
The MM-GBSA techniques were employed to determine the binding free energy of the selected four compounds in complex with the targeted protein. The desired ligand confirmation obtained through the docking process was used to estimate the binding free energy [22,30]. The selected four compounds ZINC2509501, ZINC2566088, ZINC1615112, and ZINC4104882 in complex with the targeted protein show ∆G binding (∆GB) energies value −19.91 kcal/mol, −16.69 kcal/mol, −17.98 kcal/mol, −15.34 kcal/mol (Table 2). Additionally, analysis of the physical and chemical properties found an invaluable achievement of ∆G Bind Covalent (∆GBCo), ∆G Bind Hbond (∆GBH), ∆G Bind Lipo (∆GBL), ∆G Bind Packing (∆GBP), ∆G Bind Solv GB (∆GBS), and ∆G Bind vdW(∆GBV) energies for the selected four compounds listed in Table 2. Therefore, it can be said that the compounds would be able to make a prolog binding affinity with the desired Brd4 protein.   Absorption, distribution, metabolism, and excretion (ADME), combinedly known as the pharmacokinetics properties of a lead compound, play an important role in the likelihood of success of a drug. These parameters determine the kinetic behaviors of a drug inside the body from administration to excretion through sweating, urine, or feces [14]. The volume of the drug distributed to the tissue and target site is important for its bioavailability and for reducing the side effects and toxic effects. The determination of the efficacy and toxic effect of a drug before designing it for a clinical trial can help in evaluating the possible unwanted effects in human or animal models. Thus, the evaluation of the drug candidate by using the online server Swiss ADME may help to reduce the cost by minimizing drug rejection. We evaluated the ADME properties, such as the lipophilicity, water-solubility, drug-likeness, and medicinal chemistry of our four selected compounds. The ADME properties found for all of the selected compounds in the optimized range and listed in Table 3.

Toxicity Determination
Toxicity analysis is a popular way to select an appropriate drug candidate through computer-based drug design due to its accuracy, rapidity, and accessibility for synthetic and natural compounds [12]. The TEST tool and ProTox-II are free tools that can be used to determine the toxicity of compounds. The cytotoxicity, mutagenicity, carcinogenicity, hepatotoxicity, and LD 50 (mg/kg) were calculated using the tools listed in Table 4. According to ProTox II, ZINC2509501 belongs to class 6, meaning nontoxic. On the other hand, ZINC4104882, ZINC2566088, and ZINC1615112 belong to class 5, meaning they may be toxic if 5000 mg/kg is swallowed. Table 4. Analysis of several aspects of toxicity (organ toxicity, toxicity class, Tox21-nuclear-receptor signaling pathways, Tox21-stress-response pathway, fathead minnow LC 50 (96-h), developmental toxicity, water solubility, oral rat LD 50 , and bioaccumulation factor) of the four selected compounds.

Protein-Ligand Interaction Analysis through Dynamic Simulation
Molecular dynamics simulations help atomic-scale investigations into the movement of atoms and molecules. It also provides the number of bond types and their interaction with the respective amino acids within a specific time [15]. The complex structure file obtained from the molecular docking was analyzed using the Schrödinger Desmond simulation tool.

RMSD Value Analysis
In the dynamic simulation, the number of displaced atoms was generated concerning the time-generated RMSD (root mean square deviation) value. First, the protein Cα, backbone, side chains, and heavy atoms were determined, and protein-fitting ligands were also measured at the observed duration (100 ns).

Analysis of Protein RMSD
Most of the compounds were stable based on the reference protein variability (1-3 Å) and their interactions with the ligand; values above 3 Å indicate large conformational changes in the protein, implying that the system is unstable. Analysis of the selected proteins (ZINC2509501, ZINC2566088, ZINC1615112, and ZINC4104882) and control ligand (ZINC95504909) showed that most were stable, although the ZINC2566088 compound had some fluctuation at 80 ns (4.59 Å). The other compounds were more stable (Figure 7) than the control compound (ZINC95504909).

Analysis of Ligand RMSD
The data gathered from the interactions between proteins and ligands showed that all the ligand interactions were stable. For the compounds, ZINC2566088 was unstable at 70 ns and stable from the beginning of 72 ns, but the compound did not exceed the RMSD by more than 2 Å. The selected control ligand (ZINC95504909) showed instability with respect to the protein at 30, 54.2, and 62.9 ns, and from 75.9 to 80.4 ns, it showed maximum fluctuations and, again, reached an equilibrium state. The natural compound ZINC2509501 was unstable from 25 to 39.5 ns; after that, it reached a stable state (Figure 8. From 0 to 35 ns, compound ZINC2509501 showed no fluctuations, but from 35 to 40 ns showed instability and again came to the stable expression. The stability of the ligand with respect to the protein and its binding pocket was determined by ligand RMSD analysis.

Analysis of Protein RMSD
Most of the compounds were stable based on the reference protein variability (1-3 Å) and their interactions with the ligand; values above 3 Å indicate large conformational changes in the protein, implying that the system is unstable. Analysis of the selected proteins (ZINC2509501, ZINC2566088, ZINC1615112, and ZINC4104882) and control ligand (ZINC95504909) showed that most were stable, although the ZINC2566088 compound had some fluctuation at 80 ns (4.59 Å). The other compounds were more stable (Figure 7) than the control compound (ZINC95504909).

Analysis of Ligand RMSD
The data gathered from the interactions between proteins and ligands showed that all the ligand interactions were stable. For the compounds, ZINC2566088 was unstable at 70 ns and stable from the beginning of 72 ns, but the compound did not exceed the RMSD by more than 2 Å. The selected control ligand (ZINC95504909) showed instability with respect to the protein at 30, 54.2, and 62.9 ns, and from 75.9 to 80.4 ns, it showed maximum fluctuations and, again, reached an equilibrium state. The natural compound ZINC2509501 was unstable from 25 to 39.5 ns; after that, it reached a stable state ( Figure  8. From 0 to 35 ns, compound ZINC2509501 showed no fluctuations, but from 35 to 40 ns showed instability and again came to the stable expression. The stability of the ligand with respect to the protein and its binding pocket was determined by ligand RMSD analysis.

RMSF Analysis
The root mean square fluctuation (RMSF) was used to evaluate the conformational changes between the compound and the selected ligands. All the compounds had minimal RMSF values, but most showed low instability at position 50 (LEU92). The maximum RMSF in the N-terminal was observed at a maximum of 4.52 Å. All the compounds had lower RMSF values, increasing the possibility of discovering new compounds targeting

RMSF Analysis
The root mean square fluctuation (RMSF) was used to evaluate the conformational changes between the compound and the selected ligands. All the compounds had minimal RMSF values, but most showed low instability at position 50 (LEU92). The maximum RMSF in the N-terminal was observed at a maximum of 4.52 Å. All the compounds had lower RMSF values, increasing the possibility of discovering new compounds targeting the Brd4 protein (Figure 9). Figure 8. RMSD value determined from the protein-fitting ligand interaction. Different colors of the selected compounds, i.e., ZINC4104882 (orange), ZINC2509501 (gray), ZINC2566088 (gold), ZINC1615112 (blue), and control ligand ZINC95504909 (light blue) indicated the numbers of ligands and their expression with respect to the 100 ns simulation time.

RMSF Analysis
The root mean square fluctuation (RMSF) was used to evaluate the conformational changes between the compound and the selected ligands. All the compounds had minimal RMSF values, but most showed low instability at position 50 (LEU92). The maximum RMSF in the N-terminal was observed at a maximum of 4.52 Å. All the compounds had lower RMSF values, increasing the possibility of discovering new compounds targeting the Brd4 protein (Figure 9).

Protein-Ligand Interaction Analysis
Several bonds were involved in the protein-ligand interactions, such as hydrogen bonds, hydrophobic interactions, and water bridges. These bonds are important when considering a compound as a drug molecule. In our experiment, all the compounds showed tangible contact with the targeted protein (Brd4) molecule. In ZINC2509501, the amino acids ALA89 and THR103 did not show contact with the protein ( Figure 10B).
All the selected ligands were compressed during the simulation time except the control ligand ZINC95504909, where the structural conformation changed highly ( Figure 11A). The number of intra-molecular hydrogen bonds was higher in ligand ZINC1615112 and increased the binding properties to the protein ( Figure 11B).

Protein-Ligand Interaction Analysis
Several bonds were involved in the protein-ligand interactions, such as hydrogen bonds, hydrophobic interactions, and water bridges. These bonds are important when considering a compound as a drug molecule. In our experiment, all the compounds showed tangible contact with the targeted protein (Brd4) molecule. In ZINC2509501, the amino acids ALA89 and THR103 did not show contact with the protein ( Figure 10B). All the selected ligands were compressed during the simulation time except the control ligand ZINC95504909, where the structural conformation changed highly ( Figure  11A). The number of intra-molecular hydrogen bonds was higher in ligand ZINC1615112 and increased the binding properties to the protein ( Figure 11B).  considering a compound as a drug molecule. In our experiment, all the compounds showed tangible contact with the targeted protein (Brd4) molecule. In ZINC2509501, the amino acids ALA89 and THR103 did not show contact with the protein (Figure 10B). All the selected ligands were compressed during the simulation time except the control ligand ZINC95504909, where the structural conformation changed highly ( Figure  11A). The number of intra-molecular hydrogen bonds was higher in ligand ZINC1615112 and increased the binding properties to the protein ( Figure 11B).

Discussion
Every year, more people die from cancer than other chronic diseases, such as HIV, malaria, and tuberculosis combined. In 2020, the death rate was close to 55%, compared to nearly 25% in 1960 [31]. Several genes are typically involved in carcinoma development through the deregulation of cellular growth, nucleic acid repair, chromosome instability, cell interaction and communication, angiogenesis, and changes in the apoptotic pathway [32,33].
Several factors are responsible for the delayed development of cancer among neuroblastoma patients, such as early diagnosis and molecular and genetic factors. Neuroblastoma has also been observed among children, with a mortality rate of 13% [34,35]. Nowadays, BET proteins are used as therapeutic targets because of research linking pharmacologic suppression to patterning in Brd4 localization on chromatin and/or modifications in overall genetic expression [36,37]. Brd4 is a crucial part of many illness-associated gene regulation networks. This could explain why BET proteins have recently become experimental targets in a remarkable variety of disease models, including cancer [18,37]. Neuroblastoma inhibition by targeting BET protein overexpression is a novel research finding, as previous studies revealed that several scientists tried to find a BET protein inhibitor using a combination of two drugs, some of which used recombinant protein inhibitors using a combination of two drugs with some of them using recombinant protein [38,39]. Current neuroblastoma treatment includes surgery, chemotherapy, and radiotherapy. The treatments improve the condition but are not patient-friendly due to their side effects. The MYCN oncogene has great potential to induce the development of more aggressive neuroblastomas. Most neuroblastoma cases initially respond to treatment, but the treatment ultimately fails due to drug resistance [40]. BET inhibitors are becoming resistant in cancer therapy due to increased Wnt signaling catenin-mediated MYC in AML and activation of the RAS pathway in lymphoma [1]. The proto-oncogenes MYC was transcriptionally suppressed by Brd4-NUT oncoprotein, released by the presence of BET inhibitors [41]. BET inhibitors' anti-cancer efficacy can occasionally be linked to their effects on MYC transcription. Still, in other cases, these effects are controlled by a distinct group of cancer-related genes [42].
The study indicated the possibility of developing a new class of compounds that would be more suitable for cancer patients due to the lower side effects but needs further experimental and animal studies. We first selected a suitable protein candidate from the protein data bank, and the currently available antagonists were identified from a literature search. The number of fit compounds with pharmacophore features similar to those of the attached ligand was retrieved from the online database [22]. The obtained compounds were validated by downloading all the structures of the 36 antagonists with the specific numbers of decoy sets from the Ligand Scout 4.4 software [43,44]. The ROC curve was satisfactory, and the validated models were further used to generate fit compounds. All the compounds were transferred to the PyRx software to obtain the binding affinity for the fit compounds by molecular docking [30]. The top four compounds with higher binding affinity and lower side effects, as determined by the software, were selected for further experiments. All the selected compounds were evaluated for their pharmacokinetic properties, such as absorption, distribution, metabolism, and excretion [22]. All the compounds had better absorption capability and drug-likeness properties, as there was no violation of the rule of five; they also had better bioavailability scores. The compounds with good ADME profiles were used for further toxicity analysis [45].
The molecular dynamic simulation provided us with data regarding the stability of the binding between proteins and ligands for 100 ns. The developed trajectory file was further studied to find the RMSD and RMSF values of both the proteins and ligands as well as the protein structure predictions, also the numbers of hydrophobic, hydrophilic, H-bond, ionic, and water bridges [12]. This identification of potential lead natural compounds would be able to inhibit the Brd4 protein in neuroblastoma treatment by silencing the expression of MYCN. This research focused on identifying potentially natural compounds that would be more patient-friendly, due to having fewer toxic effects and side effects, for treating MYCN overexpressing neuroblastoma patients by targeting the Brd4 protein. The crystallographic structure of the Brd4 protein was obtained from the PDB database (PDB ID: 4BJX) and filtered depending on the protein resolution (1 to 2.0 Å). The structurebased preparation protein was selected by the attachment of ligands. Meanwhile, the currently available antagonists were identified by a literature search and from the ChEMBL database (https://www.ebi.ac.uk/chembl/; accessd on 10 May 2022). In total, 36 antagonists were identified (Table S1) from online databases and rigorous literature searches. The ligands were screened and examined to determine their toxicity (through the ProTox-II server; (https://tox-new.charite.de/protox_II/; accessd on 9 May 2022) and inhibitory activity in neuroblastoma. All the antagonists underwent docking, and the best one was selected for structure-based drug design. The selected proteins (PDB ID: 4BJX) were analyzed using the LigandScout 4.4.8 Advance software for the number of hydrogen-bond donors or acceptors, hydrophilicity, lipophilic characteristics, and ionizable charges. A total of five hydrophobic and five hydrogen-bond acceptors were determined. For better results, we followed the Lipinski rule of five and added one positive feature from the protein for better model generation.

Model Validation
Validation was essential when evaluating the quality of our model through proteinligand interactions; several known active inhibitors (Table S1) were used to distinguish between the selected active and inactive compounds. Compounds with larger molecular weights and higher IC 50 values were not selected [18,25]. The identified active compounds (Table S1) were evaluated using the DUD-E decoy set (obtained from the DUD-E decoys database) for better measurement between the active and inactive compounds. The selected DUD-E decoy and our compounds were converted to an "idb" file using LigandScout before preparing an ROC curve. The quality of the curve was evaluated based on the true and false positives, enrichment factor, and GH score. The accuracy of the model was examined considering important parameters such as the active hits (AH), decoy compounds (DC), early enrichment factor (EF), overall compounds in the dataset (D), the total number of hits retrieved (TH), and goodness of hit score (GH).

Dataset Preparation and Virtual Screening
Identifying structurally similar novel and active molecules is important for structurebased drug discovery based on the prepared pharmacophore model. Freely accessible ZINC databases containing 230 million compounds were screened for the further identification of structurally similar compounds. The 2D and 3D structures and physical and chemical properties can therein be obtained for a compound [46,47]. The compounds most similar to our selected ligand features were selected. Several libraries were available in ZINC Pharmer (http://zincpharmer.csb.pitt.edu/; accessd on 11 May 2022) for dataset generation, and we used the ZINC natural compound library for searching compounds by matching the maximum features. The pharmacophore features of the selected compound were analyzed and prepared for screening depending on adding one feature and omitting exclusion volumes. A total of 136 hits were retrieved (Supplementary file renamed Hits), with the maximum query features considered for further validation. The compounds were further analyzed by ligand scout 4.4 software to generate the pharmacophore fit score.

Protein and Ligand Preparation
The macromolecular structure of the protein obtained from the PDB was prepared for molecular docking. The protein's crystal structure (PDB ID: 4BJX) attached to the 73B ligand was used as a positive control, and the attached ligand was separated using the Discovery Studio Visualizer (version 16.1.0). Adding any necessary bond and deletion of the water molecules was not part of the structural refinement process. The desired protein structure was obtained and analyzed, determining the free R value (0.185), resolution (1.59 A), and observed R value (0.158). A few bonds were missing in the currently selected protein, and we used Discovery Studio to attach a new bond by adding the force field (CHARMm). All the ligands obtained through screening the ZINC database were prepared and optimized using the OPLS_2005 force field Schrodinger Suit. The ionization pH range was set to 7.1 ± 2.

Grid Preparation and Active-Site Identification
The prepared structure of the selected protein was analyzed using UniProtKB and PrankWeb (https://prankweb.cz/; accessd on 11 May 2022) to identify the active site of the protein. Furthermore, CASTp (CASTp 3.0: Computed Atlas of Surface Topography of proteins (uic.edu)) was used to identify the number of active pockets ( Figure S1). The binding affinity of the protein and ligand depends on the presence of hydrogen bonds, lipophilic or hydrophilic interactions, and ionizable charges. Additionally, we used the PrankWeb (http://prankweb.cz/; accessd on 11 May 2022) database for active-site identification. Using the PyRx tool, the receptor grid was prepared through selecting the active site of the protein.

Binding Affinity Determination by Docking
Identified hits of the compound's 2D structure were downloaded from the NCBI databank and the PubChem compounds. All the compounds were transferred to the PyRx software, and docking was carried out using the AutoDock Vina tools. PyRx is a popular tool in drug design for selecting a drug active against different animal diseases to identify potential drug candidates. Based on the binding affinity and RMSD value, the compounds were further transferred to the BIOVA Discovery Studio Visualizer Tool 16.1.0 for analyzing the interaction [22,44].

MM/GBSA Data Calculation
The molecular mechanics generalized Born surface area (MM/GBSA) is a popular method usually used to calculate the free binding energy of ligands [48]. We used the Prime MM/GBSA package from Schrödinger (released in March 2020) to evaluate the free binding energy.

Pharmacokinetic and Pharmacodynamic Evaluation
The most important parameters to consider for treating disease are a drug's absorption, distribution, metabolism, and excretion after administration inside the body by any of the available routes (oral, nasal, intravenous, subcutaneous, and rectal) [49,50]. The metabolism of a drug by different enzymes and microbiota may affect its release and cause it to produce toxic effects or change its efficacy [51]. Computer-aided drug design has emerged as a new pathway for choosing the best candidate for drug development before proceeding to the clinical trial phase [22]. A drug's bioavailability is affected by gender, age, disease condition, lipophilicity, hydrophobicity, microbiota, body enzymes, and administration pathway [52]. We used the freely accessible database Swiss ADME (http://www.swissadme.ch/index. php; accessd on 13 May 2022) to identify the solubility pattern, enzymatic degradation, and absorption for each selected compound.

Toxicity Analysis
Measuring a drug's toxicity is important in establishing a drug candidate for drug discovery and development. Computer-based drug design has made it easy to determine toxicity profiles (hepatic failure, carcinogenicity, immune response, and membrane potential pathway) both quantitatively and qualitatively [25,53]. Toxicity measurement tools (Toxicity Estimation Software Tool, TEST version 4.2.1) are usually used to determine the toxic effects of chemicals based on their molecular structures. Our toxicity study (Table 4) determined the fathead minnow LC 50 (96 h), 48-h daphnia magna LC 50 , developmental toxicity, oral rat LD 50 , bioaccumulation factor, and water solubility (at 25 • C). The hepatotoxicity, carcinogenicity, mutagenicity, immunogenic response, and several toxicological pathways were determined using the ProTox-II server (https://tox-new.charite.de/protox_II/; accessd on 13 May 2022).

Preparation of the Protein and Ligand for Simulation
To justify our ligand's binding to the protein after it was obtained from the docking studies, we performed a 100 ns dynamic simulation. It was necessary to determine the complex's stability and predict every atom-binding movement of the ligand and protein molecules over a specific period. Our dynamic simulation was processed using the commercial software Schrödinger Release 2020-3 (Academic version) with the Linux command line. The protein-ligand interaction was first solved with the water model, and an orthorhombic box shape boundary was obtained. The complex atom buffer box calculation method was used by combining Na + and Cl − to create a 0.15 M sodium chloride solution.
The simulation was carried out with an ensemble temperature of 300 K at a pressure of 1.01325 bar, with a recording at 50 ps. The simulation was run using the OPLS-2005 force field.

Trajectory File Analysis
All the trajectory files from the simulation were analyzed using the Simulation Interaction Diagram (SID) of the Desmond module in the Schrödinger package. The simulation trajectory file provided information regarding the stability of the protein-ligand interaction complex based on the RMSD and RMSF values and ligand-protein complex. The radius of gyration among the selected ligand was evaluated to measure the extendedness of the ligand, and the intra-molecular hydrogen bond was determined.

Conclusions
The study discovered four natural compounds ZINC2509501, ZINC2566088, ZINC1615112, and ZINC4104882-which can attach to the Brd4 protein and have the possibility to block the activity of this protein overexpression and subsequently suppress the MYCN overexpression expression. All the substances had greater protein-binding affinities, ranging between −7.6 to −8 kcal/mol. The selected compounds show substantially improved pharmacophore features compared to the compound (PDB:73B) attached to the Brd4 (PDB ID: 4BJX) protein. The ADME profile of the selected four natural compounds has a good value that maintains all the rules of the Lipinski rule of five. Additionally, a computer-based toxicity study revealed that all the compounds were less harmful than the currently available medications for treating neuroblastoma patients. The MD simulations revealed the ligand's stability in its interaction with the protein complex for 100 ns. The computational study suggested four compounds that have the potential to combat neuroblastoma and may aid in reducing neuroblastoma-related deaths in children.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cimb44100329/s1, Figure S1. Structural description and the probable active site for the ligand molecule. The LigandScout software was used to determine the interaction with the ligand. Figure S2: Binding process validation of the docking results. Table S1: List of known antagonists for the target protein, with the PubChem CIDs, molecular weights, and IC 50 . The binding affinity of each antagonist is provided with the docking score. Table S2: Enrichment factor and goodness of hit score of the model are included in the validation of the pharmacophore model utilizing the GH scoring method. Table S3: The docking of all selected compounds is provided with the standard deviation.