Bio-Computational Evaluation of Compounds of Bacopa Monnieri as a Potential Treatment for Schizophrenia

Schizophrenia is a horrible mental disorder characterized by distorted perceptions of reality. Investigations have not identified a single etiology for schizophrenia, and there are multiple hypotheses based on various aspects of the disease. There is no specific treatment for schizophrenia. Hence, we have tried to investigate the updated information stored in the genetic databases related to genes that could be responsible for schizophrenia and other related neuronal disorders. After implementing combined computational methodology, such as protein-protein interaction analysis led by system biology approach, in silico docking analysis was performed to explore the 3D binding pattern of Bacopa monnieri natural compounds while interacting with STXBP1. The best-identified compound was CID:5319292 based on −10.3 kcal/mol binding energy. Further, selected complexes were dynamically evaluated by MDS methods, and the output reveals that the STXBP1-CID:5281800 complex showed the lowest RMSD value, i.e., between 0.3 and 0.4 nm. Hence, identified compounds could be used to develop and treat neuronal disorders after in vivo/in vitro testing.


Introduction
Schizophrenia is a horrible mental disorder characterized by distorted perceptions of reality. It can produce hallucinations, delusions, and profoundly disordered thought and behavior, impair daily functioning and be catastrophic. Patients with schizophrenia require lifelong care; early intervention may help manage symptoms before severe complications occur and improve the prognosis (https://www.webmd.com/schizophrenia/ mental-health-schizophrenia (accessed on 22 August 2022)).
Approximately 24 million globally, or 1 in 300 persons, have schizophrenia (0.32 percent). Adults had a rate of 1 in 222 (or 0.45%) during this era. It does not occur as frequently as many other mental illnesses. The most common times for onset are in late adolescence and the early twenties, and onset often occurs earlier in men than in women. Investigations have not identified a single etiology for schizophrenia. It is theorized that genetic and environmental factors work together to generate schizophrenia. Psychosocial variables may also influence the onset and progression of schizophrenia. The use of cannabis is associated with an elevated risk for the disease.
There are multiple hypotheses based on various aspects of the disease, some of which are attributable to very well-known mechanisms of treatment therapies. The most widely accepted neurodevelopmental theory of schizophrenia combines environmental factors and causal genes. The dopamine hypothesis is based on the fact that all conventional therapies involve anti-dopaminergic processes, and genes such as DRD2, DRD3, DARPP-32, BDNF, or COMT are linked to dopaminergic system function [1]. The glutamatergic theory has recently resulted in the first successful mGlu2/3 receptor agonistic medication and is supported by significant findings in glutamatergic system-regulating genes (SLC1A6, SLC1A2, GRIN1, GRIN2A, GRIA1, NRG1, ErbB4, DTNBP1, DAAO, G72/30, GRM3) [2].
The pathophysiology of the condition, which is reflected by the participation of genes including GABRA1, GABRP, GABRA6, and Reelin, has consequently been postulated to be modulated by GABA. In addition, multiple genes, including DISC1, RGS4, PRODH, DGCR6, ZDHHC8, DGCR2, Akt, CREB, IL-1B, IL-1RN, IL-10, and IL-1B, have been identified to be involved in the condition. These genes are implicated in immunological, signaling, and networking deficiencies [2].
There is no specific treatment for schizophrenia. Scientists have the privilege to investigate already developed bio-medical databases using advanced computational tools and techniques, which will lead to the discovery of biomarkers and possible treatments for schizophrenia-related disorders. Mejia et al. have explored the potential to reposition 70 drugs investigated in 231 clinical studies as potential candidates for repurposing pharmaceuticals for schizophrenia based on their interactions with the dopaminergic system. During dynamics simulation analysis, flunarizine-D2-like receptors demonstrated good interaction and stability [3]. A recent computational study found the compound ZINC74289318 as a serotonin 5-HT2A and dopamine D2 inhibitor [4]. Human D-amino acid oxidase (h-DAAO) may effectively act on D-serine, which is being investigated as a possible therapeutic target for the treatment of schizophrenia. 1,2,4-triazine derivatives were discovered to have h-DAAO inhibitory characteristics using 3D-QSAR modeling, docking, and dynamics analysis [5].
The lead compound Martynoside (CID:5319292) has been reported for its antioxidant properties [8], while a comparative toxicogenomic database shows literature-based evidence showing that acteoside (CID:5281800) has a therapeutic role in several diseases such as leukemia [9], inflammation [10], skin neoplasms [11], and wounds and injuries [10]. A previous study suggests that acteoside has antioxidant and neuroprotective activity, and herbs containing it are used to enhance memory [12]. Experimental data reported by Chen et al., (2020) proved that after administering echinacoside and acteoside in a rat model, the typical pathological features of osteoporosis and Alzheimer's disease were ameliorated [13]. CID:44559250 (Dehydroapateline) is documented to exhibited antiacetylcholinesterase activity [14].
Unfortunately, a one-size-fits-all natural antipsychotic does not exist. A combination of nutritional recommendations, on the other hand, may help lessen symptoms. However, what helps one person might not always work for the next. Natural compounds' medicinal properties investigation as a complementary and alternative treatment for schizophrenia could also be executed using computational tools.
The majority of molecular targets for schizophrenia are found in the three major components of the Neurotransmitter release cycle (NRC): serotonin, dopamine, and glutamate, according to studies. Therefore, we have extracted and analyzed all 51 gene interactions from NRC. Natural compounds abundant in Bacopa monnieri (B. monnieri) were chosen to target selected biomolecules to assess the interaction properties of these natural compounds as an alternative treatment for schizophrenia. B. monnieri is a plant used as a traditional natural medicine in the Asian sub-continent for centuries. B. monnieri can enhance biochemical molecule activity and protect brain cells, enhancing memory, thinking, and stress reduction could also help treat Alzheimer's disease [15,16].
Our previous investigation, including molecular interactions between human acetylcholinesterase (AChE) and butyrylcholinesterase (BuChE), enzymes with natural compounds from B. monnieri revealed that bacoside X, bacoside A, 3-beta-D-glucosylstigmasterol and daucosterol could be suitable cholinesterase inhibitors and need to be tested against neuro disorders through in vivo/in vitro experimentation [17].
Identified potential biomarker STXBP1 biochemical mechanism includes vesicular transportation and neurotransmitter section. A recent study revealed that STXBP1 mutation leads to amyloid-like fibrils in rat brains [18]. Pathogenic STXBP1 variants cause severe early-onset developmental and epileptic encephalopathy [19]. The patient-based study by Stamberger et al., (2016) found that patients with STXBP1 mutation have severe intellectual disabilities [20]. A comparative study between neurodevelopment disorders and intellectual disability patient groups revealed that the STXBP1-associated group had severe global adaptive impairments, fine motor difficulties, and hyperactivity [21].
In this study, we chose natural compounds from B. monnieri to investigate their inhibitory potency against potential biomarkers of schizophrenia found by pathway analysis and protein-protein interactions (PPI) utilizing a system biology method.

Neurotransmitter Release Cycle Biomolecules Data Mining
We obtained information on the 51 NRC genes corresponding to their UniProt IDs from the Reactome database [22] (Table S1). A computed interaction-based protein-protein network of these 51 gene targets was created after putting UniProt IDs as an input to the STRING (version 11.5) network analyzer [23]. The parameter to the PPI network was set as evidence-based, with the highest confidence level score of 0.9, with 50 interactors in the first and second shells. Generated STRING network was further deeply analyzed by Cytoscape version 3.9.1 [24]. A network analyzer tool was utilized to investigate the topology parameters, e.g., shortest path length, node degree distribution, average clustering coefficient, and average neighborhood connectivity of the PPI network [25].
The network was analyzed using tools integrated into CytoScape, e.g., Molecular COmplex DEtection (MCODE) and GO functional enrichment analysis, ClueGO. Parameters such as the threshold p-value were set as <0.05. The output proteins list generated more confident PPIN after setting 50-50 interactors in the first and second shells. The nodes thus obtained from the network were sorted based on their topological properties, e.g., degree, clustering coefficients, betweenness and bottleneck scores [12].

STXBP1 Preparation as a Receptor Molecule
The 3-Dimensional (3D) crystal structure of Human STXBP1 was not available in Protein Data Bank (PDB). Therefore, it was necessary to develop a 3D model of STXBPI_Human based on the homology modeling approach. SwissModel online workspace [26] was used to generate a 3D model, and further assessment of the model was performed by the Mol-Probity tool integrated into the SwissModel server [27] (Supplementary file S1).

Active Site Prediction by CASTp Server
The CASTp 3.0 server was used to determine the active site of the STXBP1 modeled structure. The CASTp server, in essence, accepts 3D protein structures in PDB format as input for topographic computing [28] ( Figure S1).

Natural Compounds Library Preparation as a Ligand
Structural and chemical information of available natural compounds in B. monnieri was extracted from PubChem online server, and control drug Quetiapine information was retrieved from DrugBank Database (https://go.drugbank.com/drugs/DB01224 (accessed on 24 August 2022)). A total of 123 natural compound details were extracted, and the .sdf library was obtained. The DataWarrior tool was used to predict the compounds' physicochemical properties, drug-likeness, and drugscore, as well as the mutagenic, tumorigenic, reproductive effect, and irritant properties of 123 natural compounds [29].

Energy Minimization
Receptors and natural compounds need to minimize energy to fix missing and nonstandard amino acid residues and atoms for virtual screening/molecular docking analysis. CHARMm force-field parameters were used for minimization, an integrated tool in Discovery Studio Visualizer 2021.  [30], followed by a combination of empirical free energy force field and Lamarckian Genetic Algorithm (LGA) to build receptor and compound complexes. Using AutoDock utilities, polar hydrogen atoms, Kollman united charges, and salvation parameters were added to the selected receptor molecules. Afterward, the ligand molecules were charged by Gasteiger. Grid box to cover the receptor's active site was set to X, Y, and Z coordinate of a grid point with needed values accordingly and for grid center with a default value of grid points spacing 0.375, while other parameters were left at their default values.

Molecular Dynamics Simulation (MDS)
Selected complexes based on docking results were also simulated by the GROMACS tool [31] for 100 ns. The receptor molecules' topology files were built by pdb2gmx tool after implementing the CHARMM27 all-atom force field. In the next step, the control drug and selected natural compounds topology files were created by the SwissParam server [32]. The unit cell triclinic box filled with water was used for the solvation step. The system was stabilized by minimum energy and addition of Na+ and Cl− ions and equilibrated as selected complexes required. It was followed by two-step ensembles NVT (constant number of particles, pressure, and temperature) and NPT (constant number of particles, pressure, and temperature). Both ensembles control temperature and pressure coupling resulting in constancy and system stabilization through complete simulation [33]. GROMACS has many different scripts for complexes MDS analysis and was utilized gmx rms for Root Mean Square Deviation (RMSD) [34], gmx rmsf for Root Mean Square fluctuation (RMSF), gmx gyrate for the calculation of Radius of Gyration (Rg) [35] and gmx H bond for the calculation of numbers of hydrogen-bond (H bond) formed during the interaction.

Results and Discussion
We conducted a network analysis that included the PPI of the 51 genes (Table S1) implicated in the NRC in an attempt to identify the putative target gene or protein that would be responsible for the disease's manifestation in schizophrenia. STRING and CytoScape were utilized to generate and analyze the interaction network.
After analyzing the role of identified proteins from the previously published studies, it was found that SNAP25, STX1A, and VAMP2 are components of the SNARE (soluble Nethylmaleimide-sensitive fusion protein attachment protein receptors) protein family and structurally tail-like small molecules work in an association of other biomolecules. SNAREs and related proteins are essential for vesicle docking, priming, fusion, and neurotransmitter release synchronization [36].
The network was generated by STRING, followed by clustering in three major clusters. Clusters or modules are tightly interconnected nodes in a network that form a dense subnetwork [37]. Selected genes were distributed in the top 3 clusters ( Figure 1A). It was observed that cluster third was most dense where Syntaxin-binding protein 1 (STXBP1) was participating ( Figure 1B).  We established a 3D structure of STXBP1 to explore the structural features and interaction of selected natural compounds. The FASTA format of the STXBP1 receptor sequences was obtained from the UniProt database (UniProt ID: P61764). SWISS-MODEL server 3D model generation quality assessment data revealed the overall sound quality of the STXBP1 model (Figure 2A) (Supplementary file S1). STXBP1 interacts with GTP-binding proteins, most prevalent in the brain and spinal cord and substantially enriched in axons, to regulate synaptic vesicle docking and fusion. It is required for neurotransmission and binds syntaxin, a component of the synaptic vesicle fusion mechanism, in a 1:1 ratio. It can interact with syntaxins 1-3 but not syntaxin 4. Mutation of the STXBP1 gene can lead to STXBP1 syndrome, a rare neurodevelopmental disorder caused by heterozygous variants in the STXBP1 gene and characterized by psychomotor delay, early-onset developmental delay, and epileptic encephalopathy [38].
Due to the established data, we chose STXBP1 as a crucial target protein structurally interacting with other genes in the NRC. Previous research suggests that STXBP1, a component of dysbindin, was co-localized with Munc18-1 at presynaptic terminals in primary cultured rat hippocampus neurons. Overexpression of the mammalian homolog of the unc-18 gene (munc18-1) has also been observed in the brains of schizophrenia patients.
The selection of STXBP1 is that all other top-ranked molecules have thread-like structures and perform the activity only when combined. STXBP1 has a more considerable structure than other top-ranked molecules. Additionally, network analysis revealed that STXBP1 forms a denser network interaction than other genes.
We established a 3D structure of STXBP1 to explore the structural features and interaction of selected natural compounds. The FASTA format of the STXBP1 receptor sequences was obtained from the UniProt database (UniProt ID: P61764). SWISS-MODEL server 3D model generation quality assessment data revealed the overall sound quality of the STXBP1 model ( Figure 2A) (Supplementary file S1). We established a 3D structure of STXBP1 to explore the structural features and interaction of selected natural compounds. The FASTA format of the STXBP1 receptor sequences was obtained from the UniProt database (UniProt ID: P61764). SWISS-MODEL server 3D model generation quality assessment data revealed the overall sound quality of the STXBP1 model (Figure 2A) (Supplementary file S1). In total, 94.55% of peptide residues were in the excellent area of the Ramachandran plot, while only 1.19% were in Ramachandran Outlier and 2.08% in Rotamer Outlier regions, and no bad bonds were found ( Figure 2B). The overall MolProbity score was 1.44. The observed GMQE (Global Model Quality Estimate) was 0.86, and QMEANDisCo global was ±0.05, between the standard cut-off value of 0 to 1. Both combined values can estimate the quality of target-template alignment and the template structure ( Figure 2C,D) [39].
Further, CASTp server probed the active site residues within a radius of 1.4 angstrom (Å). The resulting computation of the total Richard's solvent-accessible area is 685.426 Å 2 and with the volume of 1474.519 Å 3 ( Figure S1). The physicochemical qualities, drug-likeness, and drug score of the 123 discovered natural compounds from B. monnieri and their mutagenic, tumorigenic, reproductive impact, and irritating properties were then evaluated (Table S2).
The molecular docking of these 123 compounds was conducted against the active site residues of the modeled STXBP1. After docking, ten compounds (Table S3) with STXBP1 complexes were shown to have better binding energy (−8.3 to −10.3 kcal/mol) as compared to the control drug Quetiapine (−7.18 kcal/mol) ( Table 2).  Natural compound CID:5319292 showed the best binding affinity with −10.3 kcal/mol and formed ten h-bonds. After visualizing 2D models of the complex, it was observed that amino acid resides GLN576, ILE259, LEU573, TYR254, SER146, SER149, and THR570 were involved in the formation of hydrophobic interaction. ALA150, LEU138, and TYR140 formed the Alkyl/Pi-Alkyl bond, while ARG39, GLU260, and LYS577 formed the Pi-Cation Pi-Anion bond.
Furthermore, 100 ns long molecular dynamics simulation generated multiple output files containing data of RMSD, RMSF, Rg, and the number of H bond that occurred during MDS. After creating plots from these data sets, it was observed that RMSD values ranged from 0.3 to 0.7 nm ( Figure 4A). STXBP1-CID:5281800 complex exhibited the lowest RMSD value, i.e., between 0.3 and 0.4 nm, which was lower than that of the control throughout the simulation; this observation indicates that CID:5281800 has interacted well with STXBP1 with less deviation than other selected compounds. Interestingly, the STXBP1-Quetiapine complex and STXBP1 simulation in water showed almost similar values, approximately 0.45 nm, and a deviation pattern during MDS ( Figure 4A). Furthermore, 100 ns long molecular dynamics simulation generated multiple output files containing data of RMSD, RMSF, Rg, and the number of H bond that occurred during MDS. After creating plots from these data sets, it was observed that RMSD values ranged from 0.3 to 0.7 nm ( Figure 4A). STXBP1-CID:5281800 complex exhibited the lowest RMSD value, i.e., between 0.3 and 0.4 nm, which was lower than that of the control throughout the simulation; this observation indicates that CID:5281800 has interacted well with STXBP1 with less deviation than other selected compounds. Interestingly, the STXBP1-Quetiapine complex and STXBP1 simulation in water showed almost similar values, approximately 0.45 nm, and a deviation pattern during MDS ( Figure 4A). RMSF plot representing per amino acid residues fluctuation with the observed value ranging from 0.25 to 1 nm ( Figure 4B) and remained stable with 0.25 nm except for a few fluctuations. After analysis, the highest fluctuation of the plot was noticed at 500-550 amino acid residues regions, while some other significant fluctuation was observed at 200-230, 260-280, and 300-340 amino acid residues ( Figure 4B). Continuous fluctuation at regular intervals shows that amino acid residues present in these regions are also involved in forming H bonds and different interactions during molecular interaction analysis (Figure 4).
The H bond plot representation 1-6 hydrogen bonds formed between receptor-compound interaction during the 100 ns period ( Figure 4C). STXBP-CID:5281800 interaction formed 6 H bonds. RoG plot critical assessment is the most important for evaluating the compactness and stability of STXBP1 during the whole simulation period due to the presence of Drugs and selected natural compounds. RoG values for all studied molecules ranged between 2.55 and 2.75 nm. It was also shown that compactness was mainly maintained with a value of 2.25 nm. Compared to STXBP1, STXBP1-CID:5281800 and STXBP1-Quetiapine showed less value with better stability except for instability by STXBP1-CID:5319292 complex ( Figure 4D).
Our findings suggest that the identified natural compounds may be effective STXBP inhibitors for the possible treatment of schizophrenia following experimental validation.

Conclusions
STXBP1 has been identified as a potential druggable target for schizophrenia in this study. The binding pattern of Bacopa's natural compounds with STXBP1 was also studied, and it was revealed that only a few compounds (e.g., STXBP1-CID:5281800) have significant binding efficacy with the active site of STXBP1. Additionally, in vivo/in vitro experimental testing is necessary to confirm the pharmacological efficacy of Bacopa monnieri compounds.
Funding: This research received no external funding.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27207050/s1. SWISS-MODEL Homology Modelling Report; Figure S1: Depicting the 3D visualization of predicted active site (in red sphere) key residues are shown by sphere pattern. STXBP1 shown by grey color in ribbon pattern;