The computational study was executed using the high-performance super-computing facility, PARAM Shivay, installed at the Indian Institute of Technology, Banaras Hindu University, Varanasi, India, with a capacity of 837 TFLOPS with Intel(R) Xeon(R) Gold 6148 CPU @2.40 GHz and 40 CPUs per node. This study used key software programs including AutoDock 4.2, GROMACS 2023, Origin 2024, Chimera 1.17.3, PyMOL, Discovery Studio visualizer, VMD, and Bio3D tool with R package and key webtools including PharmMapper, pkCSM, SwissADME, and ProTox 3.0.

For screening of bioactive phytoconstituents, the 3D model of CST developed in our earlier study was used as a receptor. The model comprises 69–336 amino acid residues of the full-length protein with 423 residues and covers the entire catalytic region of the protein with a sulfuryl acceptor substrate (galactocerebroside)-binding site and a sulfuryl donor co-substrate (PAPS)-binding site on a linear horizontal plane (Singh and Singh, 2024a; Singh and Singh, 2024b). The model protein was processed using AutoDock 4.2 by the removal of water molecules and heteroatoms and the addition of hydrogen atoms with proper assignment of atom type and Gasteiger charge to generate a PDBQT file of the protein. A grid box of 90 × 90 × 90 Å and spacing of 0.253 Å around the protein active site was generated considering the key residues LYS82, HIS84, LYS85, HIS141, PHE170, TYR176, PHE177, TYR203, and ARG202, and a grid file (.gpf) was created.

We used a set of compounds from the Indian Medicinal Plants, Phytochemistry, And Therapeutics (IMPPAT 2.0) online database (Vivek-Ananth et al., 2023; Mohanraj et al., 2018). The 3D (.mol2) structures of 81 compounds from C. asiatica Linn (Mandukaparni), 310 from G. glabra Linn (Yastimadhu), 52 from T. cordifolia Miers (Guduchi), and 18 from C. pleuricaulis Chois (Shankhpushpi) were converted to .pdbqt files after energy minimization and assigning proper atom types using AutoDock Raccoon, a virtual screening file preparation tool (Forli et al., 2016).

Following the preparation of the protein, ligands, and grid files, AutoDock Raccoon was further used for the preparation of docking (.dpf) files for each ligand and the subsequent arrangement of all files in a single separate folder for each protein–ligand complex with the generation of a single virtual screening script file (.sh) for performing molecular docking-based virtual screening of all ligands simultaneously under the Linux environment. The parameters used for the generation of .dlg files for docking run were 100 GA run, a population size of 300, maximum number of generations of 27,000, and maximum number of evaluations of 25,000,000. AutoDock applied the Lamarckian genetic algorithm and a gradient-based local search method for protein–ligand interactions. The pKi and ligand efficiency are two major parameters for assessing the binding potential of ligands in the active site. The best-docked conformation was selected, processed using custom Python scripts, and visualized using .pdb visualization tools including PyMOL and Discovery Studio Visualizer to analyze protein–ligand interactions and ligand-binding patterns in the active site.

The selected top hits based on the binding score and number of conformations in the largest cluster were subjected to pharmacokinetic property analysis using pkCSM, SwissADME, and ProTox webservers (Banerjee et al., 2018; Daina et al., 2017; Pires et al., 2015). The canonical Simplified Molecular Input Line Entry System (SMILE) of the selected compounds was used as the entry data in these servers for their pharmacokinetic and drug-likeness property analysis through adsorption, distribution, metabolism, and excretion along with toxicity studies.

All atom molecular dynamic (MD) simulations were performed using the GROMACS 2023.1 software package using the CHARMM27 all-atom additive force field. For the simulation, a dodecahedron simulation box was created with a minimum distance of 1.2 Å from the box edge, and periodic boundary conditions were applied to minimize the edge effect. Each box was solvated with the TIP3P water solvation model, and the charges on the system were then neutralized by the addition of chloride (Cl⁻) ions, and thereafter, the solvated system was energy-minimized using the steepest descent algorithm. The energy-minimized system was then equilibrated to 1,000-ps (1 ns) NVT simulation with a time step of 2 fs at 300 K. Next, the system was equilibrated to 1,000-ps NPT simulations with a time step of 2.0 fs at 300 K. The LINCS algorithm was used to constrain the bond lengths. The final MD simulation run was performed for 100 ns, and trajectories were analyzed using different MD parameters including root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), hydrogen bonds, principal component analysis (PCA), and dynamic cross-correlation matrix (DCCM) analysis of the protein–ligand complexes.

The PharmMapper server was used to identify a wide range of targets using an innovative reverse pharmacophore mapping approach (Wang et al., 2017; Liu et al., 2010). The best mapping poses of the submitted molecule (.mol2) were aligned against all target proteins available in PharmTargetDB. The algorithm used to perform the reverse pharmacophore matching protocol comprises a sequential combination of triangle hashing (TriHash) and genetic algorithm (GA) optimization (Liu et al., 2010). Based on the calculated highest fit-score (cutoff >5.0) between the small compound and the pharmacophore models, the probable protein targets were ranked. It was imperative to check the disease-causing potential of the targets.

Virtual screening is a computational approach used to screen libraries of compounds available in databases against the target protein to identify a potential drug candidate for a targeted disease (Singh et al., 2018; Singh et al., 2021; Zare et al., 2024; Pirolli et al., 2023). In the CST-led research to develop substrate reduction therapy for metachromatic leukodystrophy, the state-of-the-art in silico approach could be the most promising and directional method for future studies. Toward this approach, the development of a homology model of the CST protein using various computational algorithms is an important step (Singh and Singh, 2024a). In the present study, this 3D model was utilized for screening specific and potent phytoconstituents of four major herbs of Medhya Rasayana. Out of a total of 461 bioactive compounds, 81 compounds from C. asiatica Linn., 310 from G. glabra Linn., 52 from T. cordifolia Miers, and 18 from C. pleuricaulis Chois were screened against CST, which are detailed in Supplementary Tables S1–S4. In the initial level of screening, using the lowest free energy of a binding cutoff of ≤ −7.5 kcal/mol and the number of conformation cutoff of ≥75 in the largest conformation cluster, the top 15 compounds were selected, i.e., 3 from C. asiatica and 12 from G. glabra. Compounds belonging to T. cordifolia and P. Chois were found to be less potent and less specific toward CST and, hence, were not considered for further study. The docking score of these top 15 compounds falls between −7.57 and −10.32 kcal/mol. With a ligand efficiency of ≤0.19 per non-hydrogen atom, the top 15 compounds were considered promising for drug-likeness property analysis. Apart from α- and β-carotene, all 13 compounds strictly followed the Lipinski rule of 5 parameters in terms of molecular mass (<500), number of hydrogen bond donors (<5), hydrogen bond acceptors (<10), number of rotatable bonds, and topological surface area (TPSA). Despite failing at the Lipinski rule of five parameters, α- and β-carotene were considered for further studies because of their wider therapeutic potential in neurodegenerative diseases, and these molecules have already been tested in the human body (Banerjee et al., 2023; Abrego-Guandique et al., 2023; Hira et al., 2019; Narisawa et al., 1996). Additionally, these compounds showed good lipophilicity with a log p-value > 0, which is imperative for compounds to cross the lipophilic membrane and, thus, strengthen the movement of these molecules to the target site. Therefore, all 15 hits were considered for absorption, distribution, metabolism, and excretion (ADME) and toxicity analysis to eliminate false positives from this list, which is provided in Table 1.

After the administration of the drug in the body, it goes through absorption, distribution, metabolism, and excretion processes. In this process, the drug interacts with desirable or undesirable targets and causes a pharmacological impact (Zhang and Tang, 2018). Thus, the bioavailability of a drug depends on the rate of absorption, metabolism, and transportation of the drug to the target site (Zhang and Tang, 2018; Stielow et al., 2023). Human intestinal absorption and lipophilicity of compounds are other important parameters to be considered to ensure the easy diffusion of compounds. The cytochrome p450-based drug response to body metabolism exhibits genetic variability and plays a critical role in the detoxification of drugs and homeostasis (Zhao et al., 2021). In this study, CYP2C9, CYP2D6, and CYP3A4 inhibitors were used to analyze the drug response at the metabolism stage. Renal OCT2 is a renal transporter that plays an important role in determining the renal clearance or deposition of drugs (Zou et al., 2021; Wright, 2019). Except few, most of the selected compounds showed good intestinal absorption capacity, along with body metabolism and renal excretion potential. Since MLD is a brain disorder, the major deciding parameters for the potential drug candidates were BBB permeability, which screened four compounds—IMPHY012226 (18alpha-glycyrrhetinic acid), IMPHY012473 (lupeol), IMPHY011609 (alpha-carotene), and IMPHY011707 (beta-carotene)—with considerable blood–brain barrier permeability. In order to select a potential lead molecule, in silico toxicity study becomes another crucial criterion due to its accuracy, accessibility, and rapidity in preclinical-level screening and providing a potential lead scaffold for further optimization (Raies and Bajic, 2016; Hemmerich and Ecker, 2020; Yang et al., 2018). pkCSM and ProTox 3.0 are freely available in silico toxicity study servers that use the SMILE of each compound to evaluate various toxicity parameters including Ames mutagenicity, hepatotoxicity, and cytotoxicity (Banerjee et al., 2018; Pires et al., 2015). These four phytoconstituents belonged to class “4” of the ProTox-predicted toxicity class based on their LD₅₀ range between 560 and 2,000 mg/kg, suggesting that the selected four compounds are mainly nontoxic. Table 2 provides the details of ADME and toxicity analysis of the top 15 hits, and Figure 1 shows the structural details of the best 4 compounds screened through ADME and toxicity analysis.

Protein–ligand interaction analysis was imperative to obtain an insight into the binding pattern of compounds in the substrate-binding site of the CST protein. The substrate-binding site comprises a left-end polar site dominated by Lys85, Phe170, Lys82, and Ser88; a middle part flanked on both sides by His84 and His141; and an aromatic site at the right end dominated by Tyr203 and Phe170. With the lowest free binding energy of −10.32 kcal/mol and 98 similar conformations in the largest cluster, 18alpha-glycyrrhetinic acid was found to be the most potential ligand with the highest binding affinity among the 4 in the binding pocket of the CST protein. The compound consists of five steroidal rings with a carboxylic group attached at one end and a hydroxyl group at the other end. These ends determined the orientation of the compound toward the polar site based on the degree of polarity of the carboxylic group, which interacted with Lys85 by hydrogen bonding. One oxygen atom was attached by a double bond at the third ring in the middle of the compound positioned at the middle of the active site and was flanked on both sides by His84 and His141 and formed a hydrogen bond with His84. At the right side of the binding pocket, methyl groups located on the rightmost ring of the compound interacted with the residues Tyr203 and Phe177 via pi–sigma and pi–alkyl bonds, respectively (Figure 2I). Both α-carotene and β-carotene complexes shared similarity in the interaction pattern in the CST-binding pocket. α-Carotene and β-carotene are isomers, differing in the position of the double bond on the ring at one end oriented toward the polar site in the binding pocket. Lys82 and Lys85 played a critical role in the interaction with the isomeric fragment of these compounds. In α-carotene, the methyl group at the C2 position was oriented toward Lys85 and, thus, interacted with it, whereas in β-carotene, the methyl group at the C2 position was oppositely oriented and away from Lys85. In the CST–β-carotene complex, the two methyl groups at the C6 position were oriented toward Lys85 and interacted via the pi–alkyl interaction, whereas in the CST–α-carotene complex, the methyl groups at the C6 carbon were positioned away from Lys85 and interacted with Lys82. The aliphatic chain between the two aromatic rings interacted in a nearly similar fashion in both complexes with key residues Phe177, Tyr203, phe170, and His84, while Arg202 interacted with the aromatic ring on the right side of the binding pocket (Figures 2II, III). In contrast to the other three compounds that broadly occupied the active site core, in the CST–lupeol complex, the ligand predominantly occupied the right side of the binding pocket, which might be due to its small size and the missing of a polar carboxylic group at one end as it was in 18alpha-glycyrrhetinic acid. The lack of occupancy of the polar site by lupeol is due to its five-carbon nonpolar ring that interacted strongly at the aromatic site of the binding pocket (Figure 2IV), while the six-membered polar ring in 18alpha-glycyrrhetinic acid stretched the compound toward the polar site and, thus, was strongly accommodated in the active site (Figure 2I). The details of the interaction types are given in Table 3. To further understand the protein–ligand interaction under a dynamic environment, all four compounds were considered for molecular dynamic simulations.

MD simulation is a computational approach used to analyze and optimize the overall stability of the protein–ligand complexes under atomistic simulation conditions with a dynamic aqueous environment (Guterres and Im, 2020; Kurniawan and Ishida, 2022; Khan et al., 2016). MD simulation provides a cumulative idea about the movement of every atom or atom in the protein over the simulation time span to study important biological processes, including the impact of ligand binding on the overall protein dynamics and the way the macromolecule responds at the atomic level with the binding or unbinding of the ligand (Alrouji et al., 2023; Hollingsworth and Dror, 2018). In this study, to understand the in-depth dynamic behavior of the protein–ligand interaction, MD simulation was carried out for a time span of 100 ns to evaluate the strength and stability of the protein–ligand complexes through trajectory analysis with various parameters including RMSD, RMSF, Rg, SASA, and hydrogen bonding. We also carried out principal component analysis, free energy landscape analysis, and dynamic cross-correlation analysis to understand the dominant motions responsible for the binding pattern and stability of the protein–ligand complex.

Post-simulation interaction analysis was an attempt to understand the binding pattern of the selected compounds under a dynamic aqueous environment, where the protein was free to take its most stable state. Of the four hits, 18alpha-glycyrrhetinic acid successfully retained its interaction with key residues, as well as its orientation and positioning in the binding pocket of the protein. In the CST–18alpha-glycyrrhetinic acid complex, as in the docked conformation, C=O of the carboxylic group formed a hydrogen bond with Lys85 in the polar site. In the middle of the binding pocket, His84 and His141 flank the compound on both sides, and hydrogen bond formation took place with His84, while Tyr203 at the aromatic side interacted with the compound through two pi–alkyl bonds (Figure 9A). Additionally, Ser173, as in the docked conformation, retained van der Waals interaction with the C=O group of the third aromatic ring of the compound. Thus, the CST–18alpha-glycyrrhetinic acid complex showed a strong protein–ligand interaction. During simulation, β-carotene also largely maintained a similar orientation in the binding pocket, with a slight change in the interaction pattern. Lys85 and Tyr176 in the left polar region, His84 and Phe177 in the middle, and Tyr203, Phe170, and Arg202 at the right end of the binding pocket determined the positioning of β-carotene in the binding pocket (Figure 9C). In contrast to the similarity in the docked conformation of CST–α-carotene and CST–β-carotene complexes, α-carotene was found to be the least stable in the binding pocket during the simulation. The isomeric aromatic ring of α-carotene shifted outward from the polar region away from Lys85, which was a key interacting residue in the docked conformation (Figure 9B). The simulation-led structural shift is shown in pre- and post-MD aligned complexes (Figure 10B). As in the docked conformation, the CST–lupeol complex maintained its position in the binding pocket, occupying mostly the right side of the binding pocket. Due to its relatively small size, lupeol left major room for ligand flexibility, which might not be suitable for competitive inhibition (Figure 9D). The distance of most of the pi-interactions was relatively short in both the CST–18alpha-glycyrrhetinic acid and CST–β-carotene-simulated complexes than their respective docked complexes, suggesting the strongest binding potential of 18alpha-glycyrrhetinic acid and the second strongest binding potential of β-carotene in the protein-binding pocket in the dynamic environment (Tables 3 and 4). Therefore, based on the interaction pattern analysis and alignment of pre- and post-MD complexes, 18alpha-glycyrrhetinic acid and β-carotene were found to be suitable candidates for further reverse pharmacophore analysis.

Reverse pharmacophore mapping for potential cross-target identification of the CST bioactive compounds (18alpha-glycyrrhetinic acid and β-carotene) was performed through PharmMapper, which compared the pharmacophore of the active compounds with the pharmacophore models of 300 proteins deposited in the database (Wang et al., 2017; Liu et al., 2010). Based on their fitness score cutoff ≥5.0, the key cross-targets were identified. Table 5 provides the details of the cross-target interactions of 18alpha-glycyrrhetinic acid and β-carotene. The two potential targets of 18alpha-glycyrrhetinic acid were corticosteroid 11-beta-dehydrogenase isozyme 1 and estradiol 17-beta-dehydrogenase 1, and positively, none of these are associated with the brain, while in the case of β-carotene, potential cross-targets are transthyretin, cellular retinoic acid-binding protein 2, and retinol-binding protein 4 (RBP4), among which transthyretin and RBP4 are associated with the brain.