Phytoconstituents of Withania somnifera (L.) Dunal (Ashwagandha) unveiled potential cerebroside sulfotransferase inhibitors: insight through virtual screening, molecular dynamics, toxicity, and reverse pharmacophore analysis

Abstract

Cerebroside sulfotransferase (CST) is considered as therapeutic target for substrate reduction therapy (SRT) for metachromatic leukodystrophy (MLD). The present study evaluates the therapeutic potential of 57 phytoconstituents of Withania somnifera against CST. Using binding score cutoff ≤-7.0 kcal/mol, top 10 compounds were screened and after ADME and toxicity-based screening, Withasomidienone, 2,4-methylene-cholesterol, and 2,3-Didehydrosomnifericin were identified as safe and potent drug candidates for CST inhibition. Key substrate binding site residues involved in interaction were LYS82, LYS85, SER89, TYR176, PHE170, PHE177. Four steroidal Lactone-based withanolide backbone of these compounds played a critical role in stabilizing their position in the active site pocket. 100 ns molecular dynamics simulation and subsequent trajectory analysis through structural deviation and compactness, principal components, free energy landscape and correlation matrix confirmed the stability of CST-2,3-Didehydrosomnifericin complex throughout the simulation and therefore is considered as the most potent drug candidate for CST inhibition and Withasomidienone as the second most potent drug candidate. The reverse pharmacophore analysis further confirmed the specificity of these two compounds towards CST as no major cross targets were identified. Thus, identified compounds in this study strongly present their candidature for oral drug and provide route for further development of more specific CST inhibitors.

Introduction

Metachromatic Leukodystrophy (MLD) is a rare, progressive lysosomal storage disorder caused by a deficiency or dysfunctionality of the lysosomal enzyme, aryl sulfatase A (ARSA) as a result of ARSA gene mutation [1,2,3,4,5]. Deficient ARSA activity leads to the accumulation of sulfatide glycolipids, which leads to demyelination of central and peripheral nervous systems, and thus causes severe neurological dysfunctions including motor dysfunction, developmental delay, hyperactivity, lack of communication, coordination, speech problems, behavioral disturbance, and in worst-case scenario it may lead to death [6,7,8,9]. Sulfatides play a critical role in maintaining the equilibrium of the myelin sheath that surrounds the axon and act as insulators for signal transmission for proper communication [10, 11].

Despite the wider clinical spectrum of MLD, the disease is broadly categorized into three clinical forms: late-infantile (< 30 months), juvenile (which is subdivided into early juvenile [30 months–6 years] and late juvenile [7–16 years]), and adult MLD (> 17 years) [5, 12,13,14]. The existing therapies including gene therapy, enzyme replacement therapy, hemopoietic stem cell therapy, chaperone therapy, etc. are ARSA dependent and their focus has been controlling the disease progression at the early stage when the major development of the nervous system occurs [14,15,16]. So far more than 280 mutations have been identified in the ARSA gene, thus the treatment varies from case-to-case basis [3, 7]. The high cost and risk associated with existing treatments restrict larger mass to take advantages from them. Above all, clinical data of available therapies are not inspiring and have shown poor efficacy [3, 7]. In contrast, substrate reduction therapy presents an efficient alternative with a rate-limiting strategy to control the accumulation of substrate of the deficient enzyme [9, 17]. This therapy is centered on the idea of developing an oral drug via developing inhibitors targeting the final enzyme in the biosynthesis of accumulated substrate [7]. In recent time, SRT has been successful in various other lysosomal storage disorders [6, 18,19,20]. In case of Gaucher’s disease, two FDA approved oral drugs are miglustat and eliglustat [20,21,22].

In MLD, substrate reduction therapy is in its nascent stage, although its potential has been well acknowledged [6]. For developing successful substrate reduction therapy, accumulation of sulfatide need to be countered as it gets accumulated with deficiency of ARSA [9, 17]. To control the sulfatide accumulation, development of inhibitor for Cerebroside sulfotransferase (CST), a final enzyme in the biosynthesis of sulfatide, is required to regulate its catalytic activities [7]. CST (CST; EC 2.8.2.11) is a membrane bound protein with 423 amino acid residues, catalyzes the transfer of sulfuryl group from donor co-substrate, PAPS to acceptor substrate, galactosylceramide (GC) for sulfatide synthesis [23]. CST active site comprises two binding sites: one for substrate binding and other for co-substrate binding [7, 24, 25]. To counter this sulfuryl transfer mechanism, the strategy is to develop a competitive inhibitor for the substrate, GC towards substrate binding pocket.

Plants have been used for centuries to treat various health-related complications. Even in the present time, the source of most of the modern drugs can be traced in phytoconstituents, directly or indirectly [7, 26]. Ayurveda, a very popular traditional system of medicine, practiced in India for centuries, shows the efficacy of medicinal herbs in developing therapies for enhancing memory and cognitive functioning. Withania Somnifera (L.) Dunal (Ashwagandha) is one of those potential herbs with tremendous therapeutic and nutraceutical potential. It is a widely known herb for mood and memory enhancement and promotion of cognitive functions [27,28,29,30]. Ashwagandha with its richness of bioactive phytoconstituents can be a potential herb for developing ‘oral drug’ for SRT against MLD as in traditional practice the efficacy of Ashwagandha is well acknowledged [31, 32].

In the era of computational advancement, and the availability of large datasets of natural and synthetic compounds along with numerous software, web tools, and servers for optimization of structure and pharmacokinetic studies, it becomes important to begin the study by utilizing the available cost-effective methodologies for preliminary drug screening [33,34,35,36,37,38]. In this direction, the in silico-based screening of existing phytoconstituents of Ashwagandha available in the database can be a strategy for extracting the most drug-like molecules against CST. For this to be carried out, the first requirement was the availability of a 3D model of CST protein, which our group has successfully developed using a computational modeling approach [7, 25]. The present research focuses on the screening of available phytoconstituents of Ashwagandha from the IMPPAT 2.0 (Indian Medicinal Plants, Phytochemistry, and Therapeutics 2.0) database. The idea is to screen the most suitable drug candidates for further in vitro and in vivo studies.

Materials and methods

The overall work flow of the present study is displayed in Fig. 1.

Fig. 1

The virtual screening workflow where finally two molecules screened

Resources

The computational study for present work was performed in High Performance Super Computing facility, Param Shivay of 837 TFLOPS capacity with Intel(R) Xeon(R) Gold 6148 CPU @ 2.40 GHz and 40 CPU per node. We used softwares including AutoDock MGL Tool 1.5.7, GROMACS 2023, PyMOL, Discovery studio, Chimera 1.17.3, VMD (Visual Molecular Dynamics), origin2024, GraphPad Prism 8.0, Adobe PhotoShop 6.0, etc., and web tools including pkCSM, SwissADME, PreADMET, ProTOX 3.0, FastDRH, etc.

Protein preparation and grid generation

The three-dimensional model of CST which was used in this study as a receptor, was developed by our group using a computational modelling approach [7, 25]. The model protein spans 69–336 amino acid residues comprises mainly the catalytic domain of the full-length protein sequence of 423 amino acids. The catalytic site of protein comprises two binding sites- one is the substrate (galactosylceramide) binding site and the other is a co-substrate (PAPS) binding site, both sites are opposite to each other in a linear horizontal plane of the protein. The stable model protein at 500 ns of simulation was selected for the screening of inhibitors against CST after the removal of water molecules and other heteroatoms. The refined and optimized protein structure allowed for the appropriate assignment of atom type, which is an important step in molecular docking. After assigning a gastegier charge, PDB file of model CST was converted to a PDBQT file using AutoDock MGL tool. Then the PDBQT file was imported to AutoDock Racoon for further processing [39]. For grid generation substrate binding site was selected, and a Grid box of 90 × 90 × 90 Å was created in x, y, and z coordinates with a spacing of 0.253 Å considering key active site residues of CST namely LYS82, THR83, HIS84, LYS85, ASN113, ASP114, HIS141, PHE170, TYR176, PHE177, GLY178 and HIS212 in the substrate binding site. In AutoDock MGL tool, grid (.gpf) file was created and the .gpf file was then imported into AutoDock-Racoon for the generation of docking files.

Ligand preparation

In the present study, a total of 57 phytochemicals from Withania somnifera were selected as ligand molecules. The 2D (.sdf) and 3D (.mol2) structures of these phytoconstituents were extracted from the online database IMMPAT 2.0 (Indian Medicinal Plants, Phytochemistry And Therapeutics 2.0) [35, 40]. The .mol2 files of these compounds were converted into .pdbqt files using Racoon VS file preparation software.

Molecular docking study

For docking study, AutoDock 4.2.6 software was used. First, in AutoDock Racoon, following the preparation of .pdbqt files of protein, .pdbqt files of ligands and grid file of protein, docking.dpf files were prepared for each ligand. Thereafter, Racoon arranged all these files in separate folders for each ligand along with generation of a single virtual screening script file for molecular docking to run in Linux for all ligands. Then, docking simulation of each ligand with CST model protein was performed. Docking parameters were set to 100 GA run, 300 population size, 27,000 maximum number of generations, and 25,000,000 maximum number of evaluations. Subsequently, the Lamarckian genetic algorithm was applied with a gradient-based local search method. The high throughput molecular docking was carried out in the HPC Linux system with the help of scripts for generating a grid .glg file with other conformational file and output docking .dlg file. From the output docking (.dlg) file the best-docked conformation of the ligands (with the lowest binding energy) was selected, and the corresponding protein–ligand complex was generated using custom Python scripts and pdb-tools (PyMol, Discovery Studio). Next, the protein–ligand docked complexes was analyzed, and then the substrate binding site residues involved in the non-covalent interactions with selected compounds were identified.

Prediction of in silico toxicity study and drug-likeness properties

After molecular docking, the top ten compounds were screened based on binding affinity. The PAINS filter was applied to the compounds to check their specificity and avoid binding to multiple targets. The canonical SMILE of selected compounds was submitted to pKCSM, SwissADMET, and ProTox servers [40,41,42,43]. The portals accessed the structure of compounds and provided the predicted pharmacokinetics and pharmacodynamic data for each compound. The analysis includes blood-brain barrier permeability (BBB), Human intestinal Adsorption (HIA), cytochrome p450 inhibition, toxicity (LD50), hepatoxicity, AMES test, water solubility, etc. along with physiochemical property analysis. All the ten filtered compounds were subjected to carcinogenicity and only those compounds selected which were non-carcinogenic, non-hepatotoxic, and found to be with potential blood-brain barrier permeability. These multiple filtrations were to ensure the druggability of selected compounds with unlikely off-target impact.

Molecular Dynamics (MD) simulations

The Molecular dynamics simulations of protein and protein-ligand complex were carried out to optimize and establish the thermal stability of the protein constructed by homology modeling as well as to study the stability of protein-ligand complexes [44,45,46,47,48,49]. GROMACS 2020 version was used with Charmm 27 all-atom additive force fields. The topology files of the ligands were generated by using SwissParam. The protein-ligand complex was placed at the center of a dodecahedron box with a minimum distance of 1.2 Å from the box edge to apply periodic boundary conditions and minimize the edge effect. TIP3P water model was used to solvate the system, and thereafter, the system was neutralized by the addition of chloride (Cl⁻) ions. Then, the system was energy minimized using the steepest descent algorithm, with energy minimization tolerance set at 100 kJ mol⁻¹ nm⁻¹. The system was then subjected to 1000 ps (1ns) NVT equilibration simulation with 2 fs time step and temperature set at 300 K. Subsequently, the system was subjected to 1000 ps (1ns) NPT simulation with 2.0 fs time step to equilibrate the pressure of the system to 1 bar. The bond lengths were constrained using the LINCS algorithm. Thereafter, a final equilibration simulation was performed for 100 ns with 2.0 fs time steps.

Trajectory analysis

Using the MD simulation dataset, trajectory analysis was done with parameters including root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), inter and intramolecular hydrogen bond (HB), and interaction energy using Coulombic short-range interaction energy (CIE), the short-range Lennard-Jones short range (LJ-SR), total energy required by system, principle component analysis (PCA), free energy landscape (FEL), and Dynamic cross correlation matrix (DCCM) analysis.

RMSD measures the deviation in protein structure in the presence of ligand at t simulation time relative to the structure at t = 0. The stability of the protein in the protein-ligand complex can be determined by its deviation relative to its conformation in a free state or protein-substrate complex state during the course of its simulation. The smaller the deviations, the more stable the protein structure is in the form of the complex [50, 51]. RMSD calculated using protein backbone and ligand for 100 ns simulation to check the stability of both systems. In Eq. 1, x_i^ref is the coordinate of the reference structure, x_i(t) is the coordinate of the protein structure at time t, and N is the number of atoms in the protein.

$$RMSD(t)=\sqrt{\frac1N\sum_i^N\parallel X_i(t)-X_i^{ref}}\parallel^2$$

(1)

RMSF reads the fluctuation in protein residue in protein-ligand complex and also helps to understand the region in protein with different degrees of fluctuation are impacted by ligand binding. RMSF measures the deviation in particular residue throughout the simulation [52]. In RMSF Eq. 2, T is the time over which fluctuation need to be measured. r_i^ref is the reference position of particle i.r_i(t) is the coordinate of particle i at time t.

$$\mathrm{RMSF}=\sqrt{\frac1{\mathrm T}\sum_{\mathrm i=1}^{\mathrm T}{({\mathrm r_i(t)}-{\mathrm r_i^{ref}})}^2}$$

(2)

Rg examined the compactness of the protein and protein-ligand complex and changes in protein in the presence of complex formation [53, 54]. In Rg Eq. 3, r_i displays the coordinate of atom i, r_center is the center of mass, and N is the number of protein atoms.

$$Rg=\sqrt{\frac1N\sum_{i=1}^N\;\left|r(i)-r_{center}\right|^2}$$

(3)

SASA analysis of protein-ligand complex was used to study the stability of protein folding in the presence of ligand binding. SASA describes the area of protein exposed to interact with solvent molecules [55,56,57]. SASA is a determining factor in ensuring protein stability and folding pattern and impact on them in the presence of ligands. A lower value of SASA indicates greater compactness [57].

Principal component analysis (PCA) is a multivariate statistical data reduction technique, that was applied to examine the movement of protein and conformational subspace of complex to understand the dynamic stable behavior of protein [58,59,60,61,62]. PCA plots revealed different cluster formations. PCA plot was drawn using the MD data obtained by GROMACS. We applied ‘gmx covar’ and ‘gmx analog’ utilities of the GROMACS to get the covariance matrix of the protein backbone (Cα atoms) and its diagonalization in presence of ligand. The covariance matrix produces a set of eigenvalues and corresponding eigenvectors to determine the principal components. These eigenvectors dictate the nature of the transformation of protein over the simulation time and eigenvalue determines the magnitude of transformation, thus the projection of trajectory with these first two or first three eigenvectors determined the collective protein movements [63,64,65]. In Eq. 4 of covariance, N is number of atoms in the protein, and x and y are dimensions.

$$COV_{x,y}=\frac{\Sigma(x_i-\overline x)(y_i-\overline y)}{N-1}$$

(4)

Based on PCA-led reaction coordinates, free energy landscape (FEL) was derived to understand the energy distribution of the protein folding during the simulation. FEL signifies the stability of apoprotein and protein in the form of the complex [66, 67]. The OriginPro 2024 graphical representation tool was used to visualize the free energy landscape in 3D.

DCCM measures the magnitude of correlation coefficient based on degree of atomic fluctuation of the system. Cross correlation was computed with the last 10 ns of the MD simulation trajectory of the protein-ligand complex using the Bio3D package in R studio. DCCMs is majorly characterized with strong relativity along the diagonal spreading out from diagonal and off-diagonal cross-relationships [68].

Molecular mechanics- poisson–boltzmann/generalized born surface area (MM-PB/GBSA)

The protein-ligand simulated complex was taken for calculating binding free energy using MM/PB(GB)SA methods using the FastDRH webserver, a docking post-processing tool to identify near-native binding pose [69]. For this, GAFF2 and ff14SB forcefields with TIP3P water model were used. The binding free energies of the protein–ligand complexes of the top 3 selected compounds were calculated using the PB3, PB4, GB1, GB2, GB5, GB6, GB7, and GB8 procedures for MM/PB(GB)SA calculation. The number followed by PB is a type of Poisson–Boltzmann calculation while GB is a type of generalized Born calculation. The input ligand and protein files were separate files obtained after splitting the simulated conformation file obtained from molecular dynamic simulation using GROMACS 2020. The higher degree of rigidity of the ligand-protein complex is directed by the higher negative MM-PB/GBSA values.

Results

Molecular docking analysis

As a preliminary stage of virtual screening, 57 phytoconstituents of Withania somnifera initially screened using AutoDock 4.6 against CST homology model which was generated in our earlier work [7, 25]. The details of 57 compounds with their respective binding scores are given in Table S1 in supplementary information file. Top 10 hits were selected using binding score cutoff ≤-7.0 kcal/mol and number of conformations cut off ≥ 40 in the largest cluster. The docking score cutoff of top ten compounds fell between the range of -7.15 to -10.25 kcal/mol, which strongly suggested the relatively stronger noncovalent interaction of the top ten hits towards the active site of CST protein among all 57 screened compounds. Table 1 provides details of the binding affinity of the top 10 hits along with their physiochemical properties. The physiochemical properties of these top ten hits supported their qualification on drug likeness criteria. Except Withanolide C and Withanolide S, all other compounds have molecular mass < 500. While other Lipinski rules of five criteria including number of hydrogen bond donor < 5, number of hydrogen bond acceptor < 10, number of rotatable bonds < 10, were followed by all top 10 compounds. These compounds showed remarkable binding potential in the substrate binding site of CST protein. Among 10 compounds, most belong to the Withanolide group of compounds and carry ‘withanolide’ scaffold which is a well-known bioactive scaffold with potential therapeutic actions. These top 10 compounds were subjected to pharmacokinetics studies.

Table 1 Binding affinity estimations and physiochemical properties of the top 10 compounds from molecular docking studies

PAINS and ADMET screening

To ensure the druggability against CST and avoid false positives, the top 10 compounds were subjected to PAINS and ADMET screening. The physicochemical, pharmacokinetic, and toxicity studies ensure the bioavailability and drug-likeness properties of ligands and the suitability of their future investigation. ADMET properties of the top 10 compounds were analyzed using web servers- pkCSM, SWISS-ADME, and ProTOX [41, 43]. These studies aimed to evaluate the bioavailability potential of selected compounds in the bloodstream based on the optimization parameter required for ideal drug candidates, including the ability of the compounds to be absorbed into the bloodstream, their distribution, metabolism, and timely excretion from the body without any form of toxicity. ADMET worked as a rate-limiting step in the chemical safety assessment of the top 10 hits [70]. Among them, four compounds- Withanolide R, Solasodine, Withanolide C, and Withanolide G were found unsafe in toxicity assessment studies. Among these four, Withanolide R had poor excretion potential. One more compound Withanolide A was found to be poor at excretion criteria. Among these five compounds, two (Withanolide C and Withanolide G) were found very poor in blood-brain barrier permeability. Hence, these five compounds were summarily rejected for future studies. Among the remaining five compounds, two compounds- Withanolide P and Withanolide S were found to have low human intestinal absorption (HIA) as well as poor blood-brain barrier permeation, hence they were also rejected. Finally, three compounds- Withasomidienone, 2,4-methylene-cholesterol, and 2,3-Didehydrosomnifericin were found to have good ADMET potential and PAINS pattern, and in ProTOX toxicity scale these three compounds fall into category ‘4’ which suggest no major toxicity; hence they were considered for in-depth docking interaction analysis and further molecular dynamic simulation and trajectory analysis. Table 2 describes the details of pharmacokinetics properties of the top 10 hits. As described in Table 1, all three selected compounds showed a positive result on the Lipinski rule of five under which a number of hydrogen bond donors and acceptors as well as the surface area were taken into account. The molecular weight of the compounds was 438.61 g/mol, 398.68 g/mol, and 488.62 g/mol, respectively. The hydrogen bond donors were found as 1, 1, and 4 while hydrogen bond acceptors were 4, 1, and 7, respectively for Withasomidienone, 2,4-methylene-cholesterol, and 2,3-Didehydrosomnifericin with topological polar surface areas of 63.6 Ų, 20.23 Ų, and 124.29 Ų, respectively.

Table 2 Pharmacological profile of the top 10 ligand molecules that were derived from pKCSM, SwissADME and ProTox webserver

PreMD/docking based receptor-ligand interaction analysis of top 3 compounds

Our study shortlisted three compounds that have high docking score with better in bioavailability and pharmacokinetics and pharmacodynamics properties. The interaction analysis of three drug candidates towards CST was studied in PyMOL and Discovery studio visualizer. The details of interacting residues of CST in complex with three selected compounds are tabulated in Table 3. The docking-based interaction analysis highlighted the key amino acid residues including LYS82, HIS84, LYS85, SER89, HIS141, SER173, TYR176, PHE177, and TYR203, were involved in non-covalent interaction with ligands. Among the selected three compounds, Withasomidienone showed the highest interaction affinity towards CST with the lowest free energy of binding score − 10.25 kcal/mol. In CST-Withasomidienone complex, the compound formed conventional hydrogen bonds with LYS82, LYS85, and SER89 in the polar subsite, while the compound interacted with aromatic residues HIS84, TYR176, and PHE177 in an aromatic subsite of the substrate binding site of CST via Pi-sigma and pi-alkyl bond formation. The Compound, 2,4-Methylene cholesterol with lowest free energy of binding − 9.33 kcal/mol form two hydrogen bonds with LYS82 and SER89, along with maintaining pi-interaction with TYR203, PHE170, HIS141, LYS85 and HIS84 residues in CST active site. With a binding score of -8.67 kcal/mol, 2,3-Didehydrosomnifericin formed four conventional hydrogen bonds with LYS85, LYS82, SER173, and GLY203 and one carbon hydrogen bond with SER89, while TYR176 and PHE177 involved in pi interactions. The docking-based interaction analysis is depicted in Fig. 2. The surface view of these complexes suggested the four aromatic rings of withnolide backbone of all three compounds allowed better and wider occupancy in the substrate binding site which is relatively wider in CST.

Table 3 Non-covalent bond interaction between top three compounds and CST active site residues

Fig. 2

The non-covalent interaction of the top three docked complexes, where, A, B and C indicates the surface view, protein-ligand interaction in 3D view and 2D view of CST complex with A Withasomidienone, B 2,4-methylene cholesterol and C 2,3-Didehydrosomnifericin, respectively, obtained through Discovery Studio 2024, and PyMol software version 2.5.7 package. Hydrogen bonding is represented by green, Pi-sigma by violet, Pi-Alkyl by light pink

Molecular dynamic simulation and trajectory analysis

Molecular dynamic simulation of selected compounds-CST complexes uncovered various dynamic interactions between ligand and receptor and the overall impact on protein structural stability after binding with selected ligands in aqueous environment [71,72,73]. MD simulation was carried out for each complex at the time scale of 100 ns and analyzed their trajectory via RMSD, RMSF, SASA, Rg, hydrogen bond, free energy landscape, and PCA analysis; and compared the data with the result of free CST protein and the complex of protein with its substrate, galactosyl ceramide (GC).

Structural deviations and compactness

The root mean square deviation (RMSD) was calculated to find out the amount of structural deviation in the conformation of CST while complex formation relative to the free CST during the simulation of 100 ns timespan. The RMSD graph of CST-Withasomidienone complex achieved the fastest initial stability after 5 ns while CST-2,4-methylene-cholesterol and CST-2,3-Didehydrosomnifericin attained stability after 45ns and 50 ns respectively. The average deviation was found for CST-Withasomidienone, CST-2,4-methylene-cholesterol, and CST-2,3-Didehydrosomnifericin were 0.61, 0.47, and 0.61 respectively. The average deviation for free CST and CST-GC complexes were 0.49 and 0.76, respectively at the simulation time scale of 100 ns. There had been a slight structural deviation of CST structure under CST-substrate complex formation. RMSD and probability distribution graph of all complexes highlighted the minor deviation of CST structure under complex formation (Fig. 3A). All three complexes showed the minute deviation with Withasomidienone and 2,3-Didehydrosomnifericin were restricted their RMSD between the range of the free CST and CST-GC complex, while the RMSD of CST-2,4-methylene-cholesterol complex was more closure to free CST protein. RMSD plot of all three systems suggested no major shift in the core structure of CST while complex formation with selected ligands, thus structural stability was maintained over the simulation time scale of 100 ns.

Further, root mean square fluctuation (RMSF) analysis determined the residual fluctuation after ligand binding throughout the simulation. The average RMSF of free CST, CST-substrate complex, CST-Withasomidienone, CST-2,4-methylene cholesterol, and CST-2,3-Didehydrosomnifericin were 0.15, 0.17, 0.16, 0.16, and 0.15, respectively Fig. 3B. The RMSF graph and probability distribution function (PDF) graph suggested that the residual fluctuation of selected compounds complexed with CST was similar to the fluctuation in the CST-substrate complex (Fig. 3(B)). This indicates the capacity of selected compounds as competitive inhibitors of CST catalytic activity.

Fig. 3

The representation of structural deviations in CST involves the analysis of the Root Mean Square Deviation (RMSD) and the Root Mean Square Fluctuation (RMSF) over time in CST and its complexes. A The RMSD plot displays the changes in the spatial structure of CST over a specific time scale, highlighting any deviations from its original conformation. B The RMSF distribution displays the magnitude of fluctuations of each residue in the protein structure. Lower panel of A & B represented probability distribution function of RMSD and RMSF

The radius of gyration (Rg) measured the compactness of the CST protein and its overall structural conformation before and after the binding of the ligand. The average radius of gyration of free CST was 1.76 nm, reduced to 1.71 nm after substrate (GC) binding [38]. This slight reduction of Rg in protein-substrate complex suggests that substrate binding increases the compactness of CST structure. A similar trend was observed in CST complex with selected compounds (Fig. 4(A), Upper panel). Rg of CST-Withasomidienone, CST-2,4-methylene-cholesterol, and CST-2,3-Didehydrosomnifericin were 1.72 nm, 1.736 nm, and 1.716 nm, respectively. The probability distribution graph of Rg also suggested no major shift in the structural compactness of CST under different complex formations (Fig. 4(A), Lower panel). The shifting of CST-Withasomidienone and CST-2,3-Didehydrosomnifericin complexes was near to the CST-substrate complex while the shifting of the CST-2,4-methylene-cholesterol complex was between the Rg of the free CST and the CST-GC complex. This indicated no major changes occurred after the binding of selected compounds relative to the protein-substrate binding. Overall, this was a good sign for the development of competitive inhibitors.

Fig. 4

The structural compactness of CST. A Rg plot shows the radial distribution function of CST as a function of time. This plot helps in understanding the changes in the size of CST in the presence of ligand over time and how the ligand affects its structural stability. B The SASA distribution provides information on the exposed surface area of CST and its complexes, giving insight into the protein’s exposure to the surrounding environment. The lower panels of the representation show the PDF (Probability Density Function) of the values, which helps to understand the average distribution of the Rg and SASA values for CST and its complexes

The solvent-accessible surface area (SASA) measured the average surface area of the CST protein exposed to the surrounding solvent and the impact of ligand binding on protein accessibility to solvent via structural changes. The average SASA value of free CST was 180.1 nm² which was reduced to 173. 07 nm² when complexed with the substrate, GC. Reduction in the solvent-accessible area under the protein-substrate complex suggested the compactness of the protein. A similar trend was observed in CST complexed with CST-Withasomidienone, CST-2,4-methylene-cholesterol, and CST-2,3-Didehydrosomnifericin, which possessed solvent-accessible surface area of 178.29 nm², 172.97 nm², and 176.85 nm² respectively (Fig. 4(A), upper panel). SASA of CST-Withasomidienone complex was closer to free CST while SASA of CST-2,4-methylene-cholesterol were aligned to CST-GC complex. The SASA of CST-2,3-Didehydrosomnifericin complex lied between free CST and CST-GC complex. This suggested that presence of Withasomidienone in CST least impacted the CST exposure to the solvent among all three compounds, while presence of 2,4-methylene-cholesterol in CST make similar changes as it was in protein-substrate complex, while 2,3-Didehydrosomnifericin moderately impacted the CST structure as it’s SASA lies between free CST and CST-substrate complex. Further, the probability distribution function of SASA suggested the stability of these complexes while maintaining the structural compactness of CST (Fig. 4(B), Lower panel).

Hydrogen bonds dynamics

Understanding the dynamics of intramolecular hydrogen bonds were imperative to know the conformational stability of the protein under different protein-ligand complex formations. A dynamic intramolecular hydrogen bond analysis was performed to evaluate the impact of ligand binding on the intramolecular hydrogen bond formation of CST protein. The average intramolecular hydrogen bond formation in free CST was found 175.37, which was increased to 186.13 [25]. This increase in intramolecular hydrogen bond might be due to increase in structural compactness due to the bigger size of the substrate, GC. However, in CST complexed with Withasomidienone, CST-2,4-methylene-cholesterol, and CST-2,3-Didehydrosomnifericin, intramolecular hydrogen bonds were 179.33, 179.61, and 180.20, respectively. This suggested that no significant changes in the number of hydrogen bonds within the CST protein during the formation of complex with selected compounds (Fig. 5I(A, B, C), upper panel). The slight increase in number can be attributed to the increased structural compactness. The probability distribution function of intramolecular hydrogen bonds highlighted the stability of the protein-ligand complexes (Fig. 5 (I) (A, B,C), lower panel).

Fig. 5

(I) The dynamics of intramolecular hydrogen bonding over the simulation time of 100 ns. A, B, and C in upper panel represent intramolecular hydrogen bonding of CST protein in presence of selected ligands Withasomidienone, CST-2,4-methylene-cholesterol, and CST-2,3-Didehydrosomnifericin, respectively and in lower panel represent their respective the probability density function distribution. (II) Intermolecular hydrogen bonds between CST and the selected compounds. A, B and C in upper panel represent the hydrogen bonds between protein and ligand in a protein-ligand complex of Withasomidienone, 2,4-methylene-cholesterol, and 2,3-Didehydrosomnifericin, respectively and in lower panel represents their probability density function

Further, the intermolecular hydrogen bond gives a glimpse of the interaction between CST protein and ligand. The intermolecular hydrogen bond dynamic revealed the range of hydrogen bond formation of selected compounds in CST active site was between 1 and 6. In CST-Withasomidienone, and CST-2,4-methylene-cholesterol complexes, the bond stability was maintained between 1 and 2 with minimal fluctuation (Fig. 5 (II) (A-B)). In CST-2,3-Didehydrosomnifericin complex, 1–3 hydrogen bonds showed consistent stability, while 4–6 bonds showed relatively higher fluctuation (Fig. 5 (II) (C)). This study suggested that these compounds maintained their initial binding position in the active site of CST throughout the simulation, though minor additional hydrogen bond formation occurred during the simulation due to flexibility of the amino acid residues under aqueous environment. These results gave insight into the formation and stability of the hydrogen bonds between CST and the selected compounds, providing a deeper understanding of the molecular interactions between them.

Interaction Energy (IE) analysis

Interaction energy (IE) measures the strength of binding between protein and ligand. To evaluate it, Coulombic short-range interaction energy (CIE) and Lennard-Jones short range (LJ-SR) were calculated along with total energy required by the complexes. The average CIE, LJ-SR and total protein of the free protein were − 731163.8172, 79196.46606 and − 493274.2383 kJ/mol, respectively. The average CIE for CST-Withasomidienone, CST-2,4-methylene cholesterol, CST-2,3-Didehydrosomnifericin were − 726321.2637 kJ/mol, -726461.4137 kJ/mol, and − 726577.7361 kJ/mol, respectively, which were near to the average value of -725797.1006 kJ/mol of CST-GC (Fig. 6(B)). Similarly, the average LJ-SR of CST-Withasomidienone, CST-2,4-methylene-cholesterol, 2,3-Didehydrosomnifericin was found as 76606.02026 kJ/mol, 76618.92037 kJ/mol, and 76618.92037 kJ/mol, respectively, which were closer to the LJ-SR of CST-GC complex (76318.30523 kJ/mol) (Fig. 6(A)). Thus, interaction energy results suggested that ligand binding required relatively more CIE but less LJ-SR interaction energy and thus, the strength of binding of selected compounds in CST active sites were comparable to the substrate binding. Furthermore, shifting of total energy requirement of the protein-ligand system of selected compounds showed similar tilt as was observed in substrate binding mode of CST (Fig. 6(C)). The probability distribution function of interaction energies highlighted the stability of the protein-ligand complexes (Fig. 6(A, B,C), Lower panel).

Fig. 6

Interaction energy (A) Coulombic short-range interaction energy (CIE) (B) short-range Lennard-Jones short range (LJ-SR) (C) total energy required by free energy and different protein-ligand complexes. Lower panel represents probability distribution function of CIE, LJ-SR and total energy

Principal component analysis (PCA), free energy landscape (FEL) and dynamic cross correlation matrix (DCCM) analysis

Principal component analysis (PCA) is a covariance-matrix-based mathematical technique which was used to understand the protein dynamics by capturing the most dominant motion during MD simulation based on eigenvalues of protein functions via coordinated atomic movement [58, 74, 75]. PCA analysis of free CST and its complex with substrate and three selected compounds revealed that the complexes of protein with selected compounds were confined to their minimal subspace. The atomic fluctuation and motion were represented by capturing first three principal components and representing their relationship in two plots- one plot between eigenvectors 1 and 2 (Fig. 7(A)) and a second plot between eigenvectors 1 and 3 (Fig. 7(B)). The free CST eigenvector 1 ranged between − 4 to + 4 while eigenvector 2 ranged between − 3 and 3.5 and eigenvector 3 values from − 3 to + 3. In the CST complexed with substrate, GC, the value of eigenvector ranged from − 8 to 3 nm in eigenvector 1, from − 2.7 to 2.5 for eigenvector 2 while from − 3 to 2 nm for eigenvector 3. The CST-Withasomidienone complex eigenvector 1 was valued from − 1.7 to 6 nm, eigenvector 2 from − 0.7 to 2.5, and eigenvector 3 ranged from − 2.5 to 3.0 nm. In the CST-2,4-methylene-cholesterol complex, the value ranged from − 4 to 2 nm in eigenvector 1, from − 3 to 1.5 for eigenvector 2 while from − 2.8 to 2 nm for eigenvector 3. For CST-2,3-Didehydrosomnifericin complex, the range was from − 3.8 to 5 nm in eigenvector 1, from − 3.5 to 1.5 for eigenvector 2 and from − 2 to 2 nm in eigenvector 3. Among three complexes, CST-2,3-Didehydrosomnifericin complex was found to be relatively better confined with occupying minimal subspace. Overall, the eigenvectors of CST complexed with Withasomidienone, 2,4-methylene cholesterol and 2,3-Didehydrosomnifericin were fall within the range of eigenvector of free CST and CST-GC complex, suggested that the CST maintained its conformational subspace under ligand binding stage.

Fig. 7

Representation of the two-dimensional projections of the protein conformational changes throughout the simulation trajectory using (A) the first two eigenvectors, (B) the first three eigenvectors. Free energy landscape created from 100 ns MD simulation trajectories of (C) CST-Withasomidienone, (D) CST-2,4-methylene-cholesterol, and, (E) CST-2,3-Didehydrosomnifericin complexes, and energy required by CST molecules in each complex over time. The color bar denotes the relative free energy value. Dynamic cross correlation map for protein complexed with (F)Withasomidienone, (G) 2,4-methylene-cholesterol, and, (H) 2,3-Didehydrosomnifericin, using Bio3D package in R studio

Thereafter, based on PCA, free energy landscape (FEL) analysis was performed to monitor the distinct conformation of binding as well as identifying the most dominate internal mode of motion via the principal components. Gibbs free energy expresses the global energy minima state. The FEL values ranged from 0 kJ/mol to 20.40 kJ/mol, 15.1 kJ/mol, and 13.30 kJ/mol for CST-Withasomidienone, CST-2,4-methylene-cholesterol, CST-2,3-Didehydrosomnifericin, respectively (Fig. 7 (C, D,E)). The most energetically favored region closer to the native state is represented in dark blue while the relatively unfavorable region is represented in yellow. For free CST, a free energy landscape with a global energy minimum of 16.60 kJ/mol was reported as the stable conformation, though substrate binding slightly shifted the free energy landscape with a global energy minimum of 17.10 kJ/mol [25]. As per the FEL graph, 2,3-Didehydrosomnifericin showed the most stable complex among all three compounds.

Furthermore, in order to explore the impact of selected top hits, dynamic cross correlation matrix (DCCM) analysis was performed to understand the conformational motion of amino acid residues of the binding pocket of CST protein in the presence of ligand. For DCCM analysis, last 10 ns of simulated trajectories were considered. The 2D correlation matrix ranged from − 1.0 to + 1.0 represents different degree of correlation between residues with different color from white to dark blue, where 0 to + 1 shows positive correlation while 0 to -1 for negative or anticorrelation (Fig. 7(F, G,H)). In CST-2,3-Didehydrosomnifericin complex, with relatively more areas represented with dark blue color indicate higher residue correlation followed by CST-Withasomidienone complexes. Positively, both complexes showed better correlation than protein-substrate complex. CST-2,4-methylene cholesterol showed least correlation. Conclusively, the correlated residues in the complexes supports the stability of 2,3-Didehydrosomnifericin and Withasomidienone in the protein active site, with CST-2,3-Didehydrosomnifericin complex found to be the best performing.

Post MD interaction analysis

Following the molecular dynamic simulation of the selected docked complexes and their trajectory analysis, it was essential to understand the positioning of compounds in the active site pocket of the protein. Thus, post-MD interaction analysis was performed to validate the docking result along with analyzing the degree of stability of protein-ligand complexes in dynamic environment. Under 100 ns time scale, the protein-ligand complexes were simulated, and allowed to pass through conformational changes in aqueous environment. The compound Withasomidienone maintained its original docked positioning in the active site of the protein by forming hydrogen bond with LYS82 and LYS85. Additionally, HIS84 and PHE177 maintained its interaction with aromatic ring of withanolide backbone of the compound via Pi-sigma bond (Fig. 8(A)). In CST-2,3-Didehydrosomnifericin complex, the compound maintained its initial docked position in the active site through its withanolide backbone, however, the aromatic fragment linked by methylene group oriented due to flexible rotatable bonds, thus this part of the compound was slightly shifted away from LYS85 but maintained hydrogen bond formation with LYS82. In this complex also HIS84 involved in bond formation with rings of withanolide backbone of compound (Fig. 8(C)). In post MD complex of CST- 2,4-methylene cholesterol, the compound maintained its position in the active site mainly because of the four consequent aromatic rings in similar fashion as it was in docked conformation. Like docked complex, HIS84 and HIS141 maintained their interaction via Pi-sigma bond with withanolide backbone of the compound. However, its aliphatic chain oriented outward and away from its earlier position of nearby SER173 and this caused loss of interactions with SER173, TYR202 and PHE170 which was observed in docked conformation. The Post MD interaction suggest CST- 2,4-methylene cholesterol as least stable complex among three (Fig. 8(B)). The details of the MD interaction study including residues in contact, interaction type, and distances of bonds of protein-ligand complexes of three selected compounds are given in Table 4.

Fig. 8

Post MD interaction Analysis of selected compounds in the active site of CST with surface view, protein-ligand interaction and 2D representation of interaction of (A) CST- Withasomidienone complex, (B) CST-2,4-methylene-cholesterol complex, and (C) CST-2,3-Didehydrosomnifericin complex after 100 ns of simulation. Conventional hydrogen bonds, carbon hydrogen bond, Pi-sigma bonds, Pi-alkyl bonds, and Vander Waal interactions are represented in dark green, light green, violet, light pink, and green colors, respectively

Table 4 Post MD interaction analysis of simulated protein-ligand complexes

MM-PB(GB)SA led binding energy analysis and per residue decomposition of top 3 hits

PostMD protein-ligand interaction was further evaluated based on binding energy derived from MM/PB(GB)SA analysis to estimate the free energy of binding of selected compounds towards CST protein. For this, FastDRH was used which provided eight different MM/PB(GB)SA procedures (PB3, PB4, GB1, GB2, GB5, GB6, GB7, GB8) to calculate binding free energy of protein-ligand interaction. The numbers followed by PB and GB are the type of Poisson-Boltzmann calculation and generalized Born calculation, respectively. The MM-PB(GB)SA showed better affinity of compounds than binding affinity obtained through molecular docking. In line of MD trajectory analysis, consistent trend was observed in MM-PB(GB)SA study which showed the highest binding affinity of 2,3-Didehydrosomnifericin in the active site pocket of CST protein among three test compounds, while Withasomidienone and 2,4-methylene-cholesterol showed mixed affinities. Table 5 provides the details based on different binding free energy parameters. Further, per residue decomposition analysis of CST-2,3-Didehydrosomnifericin complex showed the dominance of TYR176 (-2.08 kcal/mol), PHE177 (-1.94 kcal/mol), and ALA211 (-2.18 kcal/mol)). While in case of CST-Withasomidienone complex also, HIS84 (-2.12 kcal/mol), TYR176 (-2.36 kcal/mol), and PHE177 (-1.46 kcal/mol), were found to be the dominant residues in terms of energy decomposition. In CST-2,4-methylene-cholesterol complex, TYR176 (-2.14 kcal/mol) was found to be the dominant energy decomposed residue while the residues PHE177 (-1.36 kcal/mol) showed relatively lesser binding affinity. Thus, the per residue decomposition analysis highlighted the role of HIS84, TYR176, and PHE177 in determining the binding affinity of compounds towards CST.

Table 5 Calculated Energy parameters for the post-MD CST complexes with selected compounds

Reverse pharmacophore analysis for cross target identification of top 3 selected compounds

As a final screening these three compounds were subjected to cross-target validation using reverse pharmacophore analysis using Swiss Target Prediction. Using probability cut off ≥ 0.15 in Swiss target prediction, no potential cross target was identified for both 2,3-Didehydrosomnifericin and Withasomidienone, while 11 potential targets were identified for 2,4-methylene-cholesterol, details are provided in Table 6. The network of cross target proteins is shown in Fig. 9. Most of proteins belonged to neural cell, thus may strongly pose side effect while targeting CST with for 2,4-methylene-cholesterol.

Table 6 Reverse pharmacophore mapping for the top compounds to identify potential targets

Fig. 9

Connection between probable 11 cross targets of 2,4-methylene-cholesterol

Comparative validation of potency of 2 top performing compounds with standard sulfotransferase inhibitors

To understand the potency and efficacy of 2,3-Didehydrosomnifericin and Withasomidienone, we performed a comparative molecular docking led interaction analysis with 11 standard sulfotransferase inhibitors including Mefenamic acid, Monobutyl phthalate, Mono-cyclohexyl phthalate, Pentachlorophenol, Clomiphene, Danazol, Triclosan, Curcumin, Quercetin, Bisphenol A, and Paraben against CST. Based on binding score ≤-7.0 kcal/mol and number of conformations in largest cluster ≥ 40, Monobutyl phthalate, Mono-cyclohexyl phthalate, Danazol, and Mefenamic acid were considered for molecular dynamic simulation to study the relative stability of these compounds inside the binding pocket of CST. In preMD interaction study, mono-cyclohexyl phthalate was found to be highest in binding affinity towards CST binding pocket residues with binding score of -8.39 kcal/mol. In supplementary file, Table S2 provides the complete details of molecular docking study of these standard sulfotransferase inhibitors against CST. The docking based binding study indicated that these compounds were relatively smaller in size and the enough room were available in the binding pocket for their fluctuation (Fig. 10 (A-D)). Among four, Monobutyl phthalate and Mono-cyclohexyl phthalate were too small and occupying the corner of the left polar site in the binding pocket. To understand further the stability of these top four standard sulfotransferase inhibitors, MD simulation was performed. These compounds were largely found unstable in the active site pocket because of their relatively smaller size, as can be seen in RMSD and RMSF graphs in Fig. 10 (E & F). These comparative studies were imperative to assure the specificity and potency of the selected CST inhibitors in this study.

Fig. 10

Protein-ligand interaction analysis of CST with top scoring standard sulfotransferase inhibitors: (A) CST-Danzol (B) CST- Mefenamic acid, (C) CST- Monobutyl phthalate (D) CST-Mono-cyclohexyl phthalate complexes; MD based trajectory analysis (E) RMSD and (F) RMSF graphs

Discussion

Metachromatic leukodystrophy is a single gene neurodegenerative disorder occur due to deficiency of lysosomal enzyme, aryl sulfatase A (ARSA). Existing therapies revolve around this deficient enzyme. In contrast to these therapies substrate reduction therapy focusses on the development of inhibitor against precursor enzyme, cerebroside sulfotransferase to counter the accumulation of substrate, sulfatides. Though SRT is in its nascent stage in MLD but presents a strong alternative to the existing challenges with providing possibilities to develop oral drugs for this fatal genetic disorder. For developing substrate reduction therapy in MLD, the target protein cerebroside sulfotransferase was studied thoroughly by our group [7, 25], and we developed a computational 3D model of the CST protein with a thorough study of its binding sites and catalytic mechanism of action [25]. CST protein contains two binding sites (one for substrate binding and the other for co-substrate binding), present at the same plane to facilitate the smooth sulfotransferase activity where the co-substrate, PAPS acts as a donor of sulfuryl group and substrate, galactosylceramide in the substrate site acts as sulfuryl group acceptor. The model of CST protein spans between 69 and 336 amino acids of the full-length CST. In the present study, a newly generated 3D model of CST was used to screen bioactive compounds from Ashwagandha, one of the widely known herbs for stimulating neurological functioning. The major reason for selection of this herb was its wider acceptability in Ayurvedic form of medication against metachromatic leukodystrophy though its scientific validation is yet to be done. Thus, in this study, through screening of available phytoconstituents of Ashwagandha, we have attempted to decipher potency of Ashwagandha in developing drug candidates for developing SRT for MLD.

A total of 57 phytoconstituents of Ashwagandha showed overall good interaction with CST protein, however the strength of interaction and occupancy in the active site of the protein was varied. The binding score and stability of the compound with the largest number of conformations in a cluster with the lowest free energy of binding were used as parameters for the initial screening of compounds after molecular docking led virtual screening of 57 phytoconstituents against CST. Total 10 best compounds were selected by keeping the cutoff of lowest free energy of binding ≤-7.0 kcal/mol and number of conformations in the largest cluster cut off ≥ 40. Keeping the high cutoff for screening was an attempt to remove all possible false positives based on the interaction potential. The next stage of screening was ADMET and drug likeness-based screening which delivered 3 best compounds out of the top 10 compounds obtained from molecular docking-based screening. These three compounds were Withasomidienone, 2,4-methylene-cholesterol, and 2,3-Didehydro-somnifericin. ADMET screening showed that these compounds have good adsorption, distribution, excretion and metabolism capacity. Since in this study we are targeting neurodegenerative disease, it was essential to select compounds with good blood-brain barrier permeability. Withasomidienone, 2,4-methylene-cholesterol, and 2,3-Didehydro-somnifericin stands positively on blood-brain barrier permeability parameter. On the toxicity scale also, the selected three compounds showed no major toxicity. All three compounds belonged to withanolide groups of steroidal compounds which typically contained a four-ring ‘withanolide’ scaffold that give inherent bioactivity to these compounds.

The in-depth docking interaction analysis of three selected compounds confirmed the tight positioning of compounds in the substrate binding site of the protein. Key interactions in the active site revealed that two binding subsites were crucial for inhibitor binding mediated by two histidine residues HIS84 and HIS141. One binding subsite at left end is dominated by LYS82, LYS85, and SER89, while another binding subsite at right end is dominated by aromatic residues- TYR203, PHE170, and PHE177. Polar residues HIS84 and HIS141 stacked opposite to each other interacted with two lower rings of withanolide backbone in the compound via Pi-interaction. The surface view clearly showed the horizontal occupancy of these compounds. Withasomidienone and 2,3-Didehydrosomnifericin followed similar orientations with the positioning of their withanolide backbone closer to the aromatic site and nonpolar interactions including pi-alkyl, pi-pi and pi-sigma, etc. took place between the aromatic site and the withanolide backbone of the compound. Throughout the MD simulation compounds maintained their positioning in the active site of protein broadly in a similar fashion as it was in the docked conformations of selected compounds. In postMD interactions also polar site was dominated by LYS82, LYS85, and SER89 while the aromatic site was dominated by PHE170, TYR176, and PHE177. Vander Waal interaction was an additional interaction in the aromatic site which helped to hold the molecule in that site. Further, the superimposition of docked and post-MD docked protein-ligand complex after 100 ns molecular dynamic simulation are represented in Supplementary Figure S1 (A, B, C). The structure of Withasomidienone and 2,3-Didehydrosomnifericin were better aligned to their pre-MD docked conformations as shown in Supplementary Figure S1 (A&C). While in the CST-2,4-methylene-cholesterol complex, the compound though bound in a similar orientation but position is sifted upwards, leads its aliphatic part to protrude outwardly (Supplementary Figure S1(B)). In contrast to other two, in 2,4-methylene-cholesterol, the positioning of its withanolide backbone at the polar site was mainly due to the absence of any polar fragment on the other end of the compound and the presence of one oxygen atom on the withanolide backbone oriented the compound in the polar site. Thus, among three, major divergence was observed between pre- and post-MD interactions in CST-2,4-methylene-cholesterol complex. Without drastic changes, the conformational adjustments of 2,3-Didehydrosomnifericin and Withasomidienone showed the behavior of the compound in the active site with maintaining the integrity of protein with negligible structural shift. The detailed movement of 2,3-Didehydrosomnifericin and Withasomidienone in the binding pocket at different time interval are depicted in Figs. 11 and 12, respectively.

Fig. 11

The surface view of the CST-2,3-Didehydrosomnifericin docked complex in molecular dynamics simulation. The snapshot was taken from 0, 20, 40, 60, 80 and 100 ns, respectively, for the complex

Fig. 12

The surface view of the CST-Withasomidienone docked complex in molecular dynamics simulation. The snapshot was taken from 0, 20, 40, 60, 80 and 100 ns, respectively, for the complex

Furthermore, the trajectory analysis of simulated protein-ligand complexes revealed key factors that confirmed the stability and integrity of the protein in the presence of the ligand in its active site pocket. RMSD of selected compounds did not exceed 1.0 and fell between the range of RMSD of free CST protein and CST-substrate complex, which specified their overall integrity and stability. RMSF analysis revealed the conformational shift of the residues of protein in the presence of the selected compounds in a similar fashion as it was observed in the CST-substrate complex. Furthermore, the radius of gyration and the solvent-accessible surface area analysis validated the relatively better structural compactness in the presence of 2,3-Didehydrosomnifericin and Withasomidienone among the selected inhibitors. Thus, RMSD, RMSF, Rg and SASA analysis were instrumental in determining the stability and mechanistic aspects of the protein-ligand interactions in all protein-ligands systems [76]. Further, the impacts of ligand binding on the overall protein structure were assessed by non-covalent interaction under dynamic simulation environment. Large molecule adjust itself in the active site pocket by pushing the molecules of protein and this leads to change the overall non-covalent interactions [77]. Hydrogen bonds are considered as the most important non-covalent interaction parameter as they play a critical role in the formation of the helix and β strand [78]. In this study, the average intramolecular hydrogen bond of free CST protein was 175.37 which was increased to 186.13 when complexed with substrate, GC [25]. The increasing trend was also observed in CST complex with Withasomidienone, 2,4-methylene-cholesterol and 2,3-Didehydrosomnifericin. Since these potential inhibitors were relatively smaller than the substrate, despite showing good binding affinity with the active site residues, they showed relatively lesser impact on the overall structural integrity of the protein and it was a positive sign for becoming a competitive inhibitor with minimal impact on the core of the protein structure. Further, the protein-ligand interaction was studied through intermolecular hydrogen bond formation, interaction energy calculation and free energy landscape analysis. The interaction energy analysis through Coulombic short-range interaction energy (CIE) and Lennard-Jones short range (LJ-SR) showed that the energy requirement of protein-ligand system of selected compounds, and were aligned more towards the energy requirement of the protein-substrate complex than free protein. Through PCA and free energy landscape analysis, deviation was least observed in CST-2,3-Didehydrosomnifericin complex than other two. In dynamic cross correlation matrix, CST-2,3-Didehydrosomnifericin complex was found to be the best correlated complex, followed by CST-Withasomidienone complex, while CST-2,4-methylene-cholesterol complex was found to be the least correlated complex. Thus, overall trajectory analysis presented 2,3-Didehydrosomnifericin as the best potent performer throughout different analysis, followed by Withasomidienone.

Further, reverse pharmacophore analysis was performed to rule out possibility of cross contamination with other potential targets. This study also confirmed the potency of 2,3-Didehydrosomnifericin and Withasomidienone towards CST. While 2,4-methylene-cholesterol showed its possibility to interact 11 biological targets and that makes it a weak molecule to compete for CST inhibitor category. Therefore, 2,4-methylene-cholesterol was eliminated for further study as potential inhibitor against CST. The potency of the final two compounds (2,3-Didehydrosomnifericin and Withasomidienone) were further confirmed by comparative docking of CST with standard sulfotransferase inhibitors. Thus, the combinatorial high throughput computational studies including molecular docking, toxicity studies, molecular dynamic simulation, reverse pharmacophore screening and comparative docking approaches strongly favored the candidature of 2,3-Didehydrosomnifericin as a top performer, while Withasomidienone also presents strong candidature for a potential inhibitor of cerebroside sulfotransferase for developing substrate reduction therapy for metachromatic leukodystrophy. These compounds need further validation through in vitro inhibition studies as well as preclinical in vivo studies.

Conclusion

In the present study, virtual screening of 57 phytoconstituents of Ashwagandha retrieved from IMMPAT database was carried out against CST, a drug target for substrate reduction therapy for metachromatic leukodystrophy. The initial screening based on binding score cut off of ≤-7.0 kcal/mol delivered top 10 hits which were subjected to pharmacokinetic studies. The pharmacological estimation and toxicity analysis confirmed the better adsorption, distribution, metabolism and excretion profile of Withasomidienone, 2,4-Methylene-cholesterol and 2,3-Didehydrosomnifericin with no major toxicity possibilities. Thereafter, interaction analysis of the selected docked complexes showed their higher occupancy in the active site of CST for better interaction with protein. In substrate binding site, LYS82, LYS85 and SER89 interacted with drug candidate through polar interaction while TYR176, PHE170 and PHE177 were involved in hydrophobic or nonpolar interaction. HIS84 and HIS141 guarded the ligand in the middle from both side in the binding pocket. Among three selected compounds, Withasomidienone and 2,3-Didehydrosomnifericin interacted in similar mode where aromatic withanolide backbone interacted with aromatic region while their opposite polar part interacted with polar region of the active site. Whereas, 2,4-methylene-cholesterol interacted just opposite as it has absence of polar part in the other side. These interaction patterns were maintained in the post MD simulated complex of Withasomidienone, and 2,3-Didehydrosomnifericin, while 2,4-Methylene-cholesterol showed relatively lesser stability in the binding pocket of CST. Further, the trajectory analysis of the simulated complex confirmed the overall stability of 2,3-Didehydrosomnifericin, and Withasomidienone in the active site without disturbing the structural integrity of the CST protein. Furthermore, reverse pharmacophore screening and comparative binding with available standard sulfotransferase inhibitors confirmed the potency of 2,3-Didehydrosomnifericin and Withasomidienone towards CST. The present study is an attempt to screen the most suitable drug candidates from Ashwagandha for CST inhibition. To further validate the potency of these two drug candidates, it is imperative to confirm the result via in vitro enzymatic studies and in vivo diseased animal model study. Thus, the study will give direction to future researchers in the selection of precise target molecule for drug development against metachromatic Leukodystrophy with saving the random use of molecules with the expense of huge resources.