Several bioinformatics tools and online resources were used in this virtual screening approach to identify potential activators of ULK1. Molecular docking was performed to predict binding affinities between phytoconstituents and ULK1 protein using InstaDock v1.2 (Mohammad et al., 2021). The docking screening generated output evaluating the interaction strength and orientation of each compound within the ULK1 binding site. Complex protein-ligand complexes were generated for structural visualization and interaction analysis using PyMOL (DeLano, 2002) and Discovery Studio Visualizer (Visualizer, 2005) in a three-dimensional (3D) and two-dimensional (2D) format. Deep-PK (Myung et al., 2024) was used to predict the pharmacokinetic parameters of the screened compounds, assessing their ADMET properties. The PASS server (Lagunin et al., 2000) was used to explore the potential biological activities of the compounds. The phytochemical dataset was retrieved from the IMPPAT 2.0 (Indian Medicinal Plants, Phytochemistry, and Therapeutics Database), a database that contains plant-derived bioactive compounds (Vivek-Ananth et al., 2023). The IMPPAT 2.0 database was chosen due to its comprehensiveness and curation, featuring over 18,000 phytochemicals from Indian medicinal plants. It includes 3D structures, ethnopharmacological relevance, and drug-likeness annotations, making it highly suitable for virtual screening in neurological disorders like ALS.

The RCSB Protein Data Bank (Burley et al., 2022) was used to retrieve the protein structure with PDB ID 6QAS, having a resolution of 1.75 Å (Chaikuad et al., 2019). Co-crystallized ligands and water molecules were removed from the structure using PyMOL software, and the structure was refined to prepare it for molecular docking. Prior to docking, the ULK1 protein structure was energy minimized using Swiss-PDB Viewer to resolve steric clashes and ensure structural optimization. PDB ID 6QAS was selected for its structural completeness, co-crystallized ligand reference, and optimal resolution of 1.75 Å without mutation. Other structures like 4WNO (1.56 Å) were avoided due to introduced mutations that could bias docking/simulation interpretations. The phytochemical library was extracted from the IMPPAT 2.0 database, which contains 3D conformations in PDBQT file format. Compliant with Lipinski’s Rule of Five for drug-likeness and bioavailability, compounds were filtered based on molecular weight, lipophilicity, and the stereochemistry of hydrogen bonds. Molecules were energy-minimized to avoid steric clashes and optimize structural conformation using the Swiss-PDB Viewer (Kaplan and Littlejohn, 2001). The processed compounds were then uploaded into InstaDock for molecular docking screening against ULK1. Ligand Efficiency (LE) was calculated as the docking score divided by the number of heavy (non-hydrogen) atoms.

Molecular docking is a widely used computational technique for identifying compounds that interact with the binding site of a target protein (Hassan et al., 2017). It predicts the binding affinity and ligand efficiency of molecules for use as potential drug candidates. Structure-based virtual screening expedites this process by rapidly analyzing extensive compound libraries to identify suitable candidates based on their docking scores and interaction features with a specific target. InstaDock v1.2 was used for docking screening in a defined search space to identify compounds with potential binding sites on the ULK1 protein (Mohammad et al., 2021). The grid was defined with an X dimension of 28 Å, a Y dimension of 32 Å, and a Z dimension of 22 Å based on the predefined binding site (Liu et al., 2023). The center was positioned at coordinates 7.88 −3.917 and 21.722 to encompass the entire protein surface. A grid spacing of 1 Å was used to ensure precise coverage. The receptor was used as a three-dimensional structure of ULK1, and compounds from the phytochemical library were docked to determine their affinities. The output files were extracted from the docking results, and the compounds were compared based on their binding scores. The best candidates with promising binding affinities were then identified for subsequent evaluations. To validate the docking protocol, the co-crystallized ligand in 6QAS was re-docked into the binding site using InstaDock. The RMSD between the docked pose and the crystal structure was 0.897 Å, indicating good agreement and validating the docking protocol.

In drug discovery, the pharmacokinetic and toxicity profiles of compounds are crucial in determining their suitability for therapeutic use. Here, Deep-PK (Myung et al., 2024) was used to characterize the ADMET properties of the screened compounds from the docking analysis. Deep-PK thresholds aligned with CNS drug development requirements, focusing on high GI absorption, BBB permeability, and CYP450 inhibition avoidance. The carcinogenic potential of these compounds was evaluated using the CarcinoPred-EL web server (Zhang et al., 2017a). PAINS filtering was applied to avoid Pan-assay interference compounds and reduce the number of false positives, thereby selecting compounds with better pharmacokinetic properties and reduced toxicity risks.

PASS (https://www.way2drug.com/passonline/) is a web-based application that predicts the biological activity spectrum of a compound based on its molecular structure (Lagunin et al., 2000). PASS analysis was performed to explore the potential biological activities of the screened compounds. PASS analysis predicts the active (Pa) or inactive (Pi) probability of a compound for a particular biological function. This suggests that a higher Pa value demonstrates an increasing likelihood of biological activity. PASS analysis was used to determine if there is a possible correlation between the shortlisted compounds and ULK1 activation in ALS, allowing for the prioritization of promising candidates for further assessment. For PASS, predicted activities were accepted if Pa > Pi, indicating functional relevance.

After PASS analysis, molecules with promising results were then subjected to molecular interaction studies to determine their binding to the ULK1 binding pocket. Ligand-protein interactions were visualized from docking output files obtained using InstaDock. Molecular interactions were examined using PyMOL (DeLano, 2002) for the screened molecules, and 3D and 2D diagrams were drafted in Discovery Studio Visualizer (Visualizer, 2005) to assess various interactions, including hydrogen-bonding, hydrophobic interactions, and π-π stacking. The interactions that occurred within the ULK1 binding site between the compounds were the most favorable. To get more detailed information about the stability and dynamic behavior of the selected compounds in the ULK1 complex, we performed MD simulations so that we could explore the binding orientation and structural stabilities of the ligand-protein at the timescale of the simulation. BL-918 was selected as a reference molecule due to its validated ULK1 activation activity and reported efficacy in ALS models, serving as a pharmacological benchmark. The top binding pose was selected based on the lowest docking score, favorable hydrogen bonding, and structural similarity to BL-918 within the ULK1 binding site.

To evaluate the electronic characteristics of the selected compounds, density functional theory (DFT) calculations were conducted using the ORCA software suite (version 6.0.1) (Neese et al., 2020). Molecular structures were constructed in the Avogadro 1.2 package (Hanwell et al., 2012), and input files for calculations by DFT were prepared in XYZ format. Geometry optimizations were carried out using the B3LYP hybrid functional (Becke, 3-parameter, Lee-Yang-Parr) with tight convergence for SCF cycles. Vibrational frequency calculations were performed for each optimized geometry to ensure that each structure represented a true energy minimum and that the thermodynamic parameters (energy and entropy) were obtained. From this study, important electronic descriptors including the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO) and optimized geometries were obtained.

MD simulation is a method for analyzing the time-dependent (dynamical) behavior of atoms and molecules in the context of Newtonian mechanics (Shukla and Tripathi, 2020). Understanding protein-small molecule interactions, including complex stability and dynamic conformational changes, is crucial for drug discovery; therefore, MD simulations are utilized to gain deeper insights. MD simulations of ULK1 and its complexes with selected compounds were conducted using the GROMACS 2020 β suite (Van Der Spoel et al., 2005). The GROMOS 54A7 force field (Schmid et al., 2011) was used in this MD simulation to define atomic interactions, allowing for accurate energy calculations and MD predictions. Individual molecular structures were processed for MD simulation using the PRODRG web server (Zhu, 2019) to obtain force field parameters for the screened compounds. A water model SPC216 (Mark and Nilsson, 2001) was used to solvate the systems. Each system was confirmed to be well-structured through energy minimization using the steepest descent method, which avoided steric clashes and enabled the achievement of a stable starting structure. The system was then equilibrated by increasing the temperature from 0 K to 300 K over a period of 100 ps to stabilize molecular motion following energy minimization. Equilibration was performed in two steps using the NVT (number of particles, volume, temperature) ensemble to simulate 100 ps, ensuring temperature equilibration, followed by the NPT (number of particles, pressure, temperature) ensemble, which is similar to NVT but allows the system’s pressure and density to change. This ensured that the system remained thermodynamically stable under the conditions of the production run. Each system was then run to study the time evolution of the molecules for at 300 ns. Various parameters, including solubility, binding free energy, kinetic stability, and structural integrity of ULK1-ligand complexes, were evaluated. While these observations define the molecular engagement of candidate compounds with ULK1, they provide additional evidence supporting their assessment as potential drug candidates.

Principal component analysis (PCA) is a widely used statistical technique in computational biology, as well as chemistry and structural bioinformatics (Yang et al., 2009). It is typically applied to remove high-dimensional complexity in datasets while retaining the most relevant information about the variations of a system. When it comes to results from MD simulations, PCA is applied to summarize atomic fluctuations from MD simulations and extract a few collective motions in the large-scale conformational changes of biomolecules. This method can then determine the dominant modes of motion that inform protein flexibility and rigidity before and after ligand complex formation. PCA was performed on the MD trajectories using the gmx covar and gmx anaeig modules of GROMACS software to investigate the essential dynamics of ULK1 in both unbound and bound states with the selected phytoconstituents. The first few principal components can describe most of the atomic fluctuations, which can be considered to characterize the dynamics and stability of the protein-ligand molecular complex.

Thermodynamic stability and folding dynamics of ULK1 and its complexes were analyzed using free energy landscape (FEL) analysis (Abdelsattar et al., 2021) with the gmx sham module in GROMACS. We subsequently employed the FELs to assign the conformational states of low-Gibbs-energy systems as a function of the projections of their principal components. Our studies reveal the stable and metastable states of ULK1 and how ligand binding reshapes the global energy landscape representation of ULK1. We next evaluated the stability of scaffolding, conformational transitions, and ligand-induced allosteric modulation of free versus ligand-bound ULK1 by analyzing the FELs of both states. FEL analysis of PCA integrated within FEL analysis provided significant insights into ULK1 activation, particularly about physiological motions and energy landscapes of ULK1.

The binding free energies of protein-ligand complexes were evaluated by the Molecular Mechanics Poisson–Boltzmann Surface Area (MM/PBSA) method (Genheden and Ryde, 2015). We used the gmx_MMPBSA program to carry out four of these calculations over selected trajectory fragments that stemmed from GROMACS simulations (Valdés-Tresanco et al., 2021). To do so, we extracted 10-ns segments from each protein-ligand complex trajectory of a 300-ns production run and performed energy evaluations every 0.1 ns. This via MM/PBSA approach allowed for decomposing the various energy components such as van der Waals terms, electrostatic terms, internal molecular energies, the polar and nonpolar solvation, and the electrostatic component of the solvation. ΔG _binding is calculated from the equation: ∆ G b i n d i n g = G c o m p l e x   –   G r e c e p t o r + G l i g a n d where G _complex is the total energy of the protein-ligand system, G _receptor is the energy of the isolated protein, and G _ligand represents the energy of the unbound ligand.

A systematic virtual screening was conducted on a phytochemical dataset comprising ∼18,000 compounds derived from the IMPPAT 2.0 database. The first stage of this screening process consisted of applying the Lipinski rule of five to filter the compounds, which defines drug-likeness in terms of molecular weight, hydrogen bond donors and acceptors, and lipophilicity. This screening resulted in a filtered library of 11,908 phytochemicals. These phytochemicals were docked with the ULK1 protein using InstaDock to determine correlated binding affinities. Docking analysis yielded different conformations of each ligand, ranked by their binding energies and interactions within the ULK1 binding site. The compounds were ranked based on docking scores, and the top 10 compounds were selected, with binding affinities ranging from −11.0 to −10.3 kcal/mol, indicating a promising interaction with ULK1 (Table 1). Importantly, all chosen compounds exhibited considerably stronger binding affinities compared to the reference ULK1 activator, BL-918 (Liu et al., 2023), which has an affinity of −9.4 kcal/mol.

The pharmacokinetic and toxicological properties of small molecules are key considerations in the drug discovery and development process (Pagan et al., 2019). ADMET is an analysis of potential therapeutics that helps estimate their drug-likeness, bioavailability, and safety (Ferreira and Andricopulo, 2019). The computational methods can also be successfully used to predict these properties, thereby saving time and cost in experimental studies. The pharmacokinetic parameters of the selected 10 compounds were further screened with Deep-PK (Myung et al., 2024). From the screened ten compounds, AMDET analysis identified four compounds, Withametelin B, Candidine, Jervine, and Delavinone, with promising pharmacokinetic profiles (Table 2). The compelling finding from the study was that all four phytoconstituents exhibit permeability across the blood-brain barrier (BBB), a key physicochemical property for the discovery of drug molecules targeting ALS.

PASS analysis was performed to explore the potential biological activities of the screened hits, Withametelin B, Candidine, Jervine, and Delavinone. Evaluating the biological activities of small molecules is a crucial step in drug discovery, as it predicts whether compounds will exhibit the intended pharmacological attribute. The PASS analysis can be used to extract two important values: Pa, the likelihood that the compound is active for a particular function, and Pi, the possibility of inactivity (Lagunin et al., 2000). PASS analysis indicated that two of the four molecules screened, Candidine and Delavinone, have relevant biological activities. These compounds exhibited a range of Pa values of 0.944–0.143, indicating considerable biological potential relevance in ALS and kinase-associated pathways (Table 3).

Interaction analysis was conducted to investigate the binding mechanism of Candidine and Delavinone with ULK1. The docked confirmation files of both compounds were extracted and analyzed for their interaction with key residues in the ULK1 structure. Since the binding site is essential for the functional activity of ULK1, the binding behavior of Candidine and Delavinone at this site was of particular interest. The structural studies revealed strong binding complementarity between both compounds and the binding site pocket of ULK1, providing a favorable binding orientation (Figure 1). It was demonstrated that Candidine and Delavinone fit snugly within the binding site cavity, indicating potential activation of ULK1 (Figure 1A). A detailed investigation into molecular interactions revealed that both compounds exhibit persistent interactions with corresponding amino acid residues, including Arg18, Leu21, Val29, Cys47, Ile48, Ser86, and Tyr89 (Figure 1B). These interactions closely resembled those of known ULK1 activators, endowing these compounds with the potential to act as modulators of ULK1 (Ouyang et al., 2018). Structural analysis further revealed that Candidine and Delavinone bound inside the deep binding pocket of ULK1, suggesting that these compounds may act as activating ligands and mediate a modification in ULK1 that promotes its activation (Figure 1C). These results demonstrate the ability of Candidine and Delavinone as potential ULK1-targeting molecules, which merit further therapeutic evaluation.

The crystal structure of ULK1 (PDB ID: 6QAS) includes a well-defined ATP-binding pocket and key residues such as Arg18, Leu21, Cys47, Ser86, and Tyr89 that contribute to ligand recognition. These residues were found to be critical in ligand interactions with Candidine and Delavinone, suggesting the importance of the kinase domain in stabilizing the ligand-protein complex. The binding modes of Candidine and Delavinone were compared with those of the reference ULK1 activator BL-918 using Discovery Studio Visualizer to perform a detailed interaction analysis. We generated 2D interaction diagrams to illustrate all potential molecular interactions, providing insight into the key residues involved in binding (Figure 2). The docking analysis revealed that both Candidine and Delavinone formed multiple interactions within the binding site of ULK1, thus exhibiting strong binding affinity towards the target region (Figures 2A,B). The compounds were found to have similar key binding interactions with BL-918, as observed in the plots, suggesting a similar mode of action (Figure 2C).

DFT is a powerful tool for developing novel organic corrosion modulators in modern drug discovery (Obot et al., 2015). We performed DFT calculations for the optimized geometries and electronic properties of Candidine, Delavinone, and BL-918. These evaluations were made based on the HOMO and LUMO, as a measure of the electronic reactivity of the compounds. HOMO is the highest energy orbital accessible to the system to donate an electron, while LUMO is the lowest energy receiving orbital to accept an electron in the process. The energy gap of HOMO and LUMO (ΔE H-L) gives information about the stability and reactivity behaviour of the compound. A narrow energy gap indicates a high reactivity, because a low energy is required for the transition of electrons, while a higher gap, less reactivity and more stability. It can be seen from Figure 3 that Candidine had higher reactivity with a HOMO-LUMO energy gap of 2.5570 eV. The Delavinone HOMO–LUMO gap was 4.9119 eV, while 4.2471 eV for the HOMO–LUMO gap of BL-918. The HOMO and LUMO distributions and their energy level are shown in Figure 3. The orbital maps show the electron profile, and blue and pink correspond to the positive and negative lobes of the molecular orbitals.

MD simulations were performed to elucidate the thermodynamic properties and binding behavior of ULK1 in its free and complex forms with Candidine and Delavinone. The simulations were performed for 300 ns, providing detailed insights into the conformational stability and flexibility of ULK1 after its binding to the ligands. The stability of ULK1 in the presence of the Candidine and Delavinone was further examined using the various structural parameters described in the following subsections.

PCA is a widely used method for essential dynamics analysis, which identifies the conformational movements of proteins and protein-ligand systems, reducing high-dimensional data into principal components that capture the dominant motion patterns of the biomolecule (Yang et al., 2009). We performed PCA to analyze the collective motions of the ULK1 protein and its complexes with Candidine and Delavinone during simulations (Figure 8). The PCA analysis revealed that the ULK1 protein and its ligand-bound forms exhibited distinct projections along two principal eigenvectors (EVs) calculated from the Cα atomic displacements (Figure 8A). The enrichment analysis revealed that ULK1-Delavinone exhibited a more compact projection than ULK1, ULK1-Candidine, and ULK1-BL-918. The ligand-bound complexes showed stable conformational motions, with minor deviations from the conformation of the native state, even with diverse normative trajectories (Figure 8B). The PCA results showed clustering of the ligand-bound complexes compared to free ULK1, indicating ligand-induced restriction of conformational mobility. Such clustering reflects stable and restricted motions, consistent with binding-induced stabilization.

FEL analysis is another valuable technique for evaluating the impact of ligand binding on protein folding and conformational stability (Abdelsattar et al., 2021). FEL provides insights into the folding dynamics by mapping the energy states of a protein through simulations, revealing the most stable regions and the extent of reshaping in response to ligand binding. To assess the interaction of Candidine and Delavinone with elongating ULK1, FELs were generated from the PCA trajectories to visualize the global minima and folding mechanism of ULK1. Energy links to ligand-bound and unbound conformations of ULK1 that are presented in deep blue on the contour plots represent the states that are energetically preferred and nearest to the respective native state (Figure 9). ULK1 presented a single global minimum, stable in 1-2 basins (Figure 9A). In complexes with Candidine and Delavinone, the expansion of energy basins was observed in the FELs, where ULK1 adopted diverse stable conformations with two to three global minima (Figures 9B,C). At the same time, ULK1-BL-918 also adopted diverse stable conformations with a single global minimum (Figure 9D). The FEL structural analysis revealed distinct conformational preferences of ULK1 in its unbound and ligand-bound states (Figure 9, lower panels). Structural snapshots extracted from the global minima highlight the most stable conformations, offering mechanistic insights into how each ligand modulates ULK1 stability and flexibility.

To evaluate the binding affinities of the studied ligands for ULK1, we carried out MM/PBSA calculations for the protein-ligand complexes using the production trajectories. These studies dissected the binding energy (ΔG _binding) into contributions from the gas-phase (van der Waals (vdW) and electrostatic) and solvation (polar and non-polar) effects. The results show that Candidine and Delavinone exhibited better binding energy as compared to BL-918 (Table 6). Concretely, ULK1-Delavinone reached the lowest binding free energy (19.99 ± 3.73 kJ/mol), with the next being ULK1-BL-918 (−18.26 ± 2.71 kJ/mol) and ULK1-Candidine (−6.93 ± 2.79 kJ/mol). The negative values indicate a strong interaction and stability of ULK1 with the selected phytochemicals during the binding process. Overall, the data signify the possibility of Candidine and Delavinone as potential ULK1 binders are worth deeper validation.