Unveiling the Efficacy of Sesquiterpenes from Marine Sponge Dactylospongia elegans in Inhibiting Dihydrofolate Reductase Using Docking and Molecular Dynamic Studies

Dihydrofolate reductase (DHFR) is a crucial enzyme that maintains the levels of 5,6,7,8-tetrahydrofolate (THF) required for the biological synthesis of the building blocks of DNA, RNA, and proteins. Over-activation of DHFR results in the progression of multiple pathological conditions such as cancer, bacterial infection, and inflammation. Therefore, DHFR inhibition plays a major role in treating these illnesses. Sesquiterpenes of various types are prime metabolites derived from the marine sponge Dactylospongia elegans and have demonstrated antitumor, anti-inflammation, and antibacterial capacities. Here, we investigated the in silico potential inhibitory effects of 87 D. elegans metabolites on DHFR and predicted their ADMET properties. Compounds were prepared computationally for molecular docking into the selected crystal structure of DHFR (PDB: 1KMV). The docking scores of metabolites 34, 28, and 44 were the highest among this series (gscore values of −12.431, −11.502, and −10.62 kcal/mol, respectively), even above the co-crystallized inhibitor SRI-9662 score (−10.432 kcal/mol). The binding affinity and protein stability of these top three scored compounds were further estimated using molecular dynamic simulation. Compounds 34, 28, and 44 revealed high binding affinity to the enzyme and could be possible leads for DHFR inhibitors; however, further in vitro and in vivo investigations are required to validate their potential.


OR PEER REVIEW
It is noteworthy that sesquiterpenes are the prime constituents separated from D. elegans that revealed cytotoxic, antitumor, anti-inflammatory, and antibacterial capacities (Table 1) [29].    [29].
The docking results of the study compounds were listed in Table 2, where they were ranked according to their gscores from highest to lowest in free energy of binding; the more negative scores imply better binding. The top-ranked compounds were 34, 28, and 44, with gscores of −12.431, −11.502, and −10.62 kcal/mol, respectively. The gscores of these compounds exceeded the value of the reference inhibitor SRI-9662, which had a binding  Table 2). The detailed docking parameters and scores are presented in Table S1. It is noteworthy that sesquiterpenes are the prime constituents separated from D. elegans that revealed cytotoxic, antitumor, anti-inflammatory, and antibacterial capacities (Table 1) [29].  [29].

Activity
Compound Name  (1), smenospongine (25) For validation purposes of the docking method, the co-crystallized inhibitor SRI-9662 was redocked inside the active site of the prepared hDHFR, and the original and redocked inhibitors were superimposed. By comparing the two results, the redocked inhibitor produced a nearly identical pose to that of the original crystal structure. The calculated rootmean-square deviation (RMSD) value of the superimposition was 0.2105 Å within the acceptable range ( Figure 6). The docking results of the study compounds were listed in Table 2, where they were ranked according to their gscores from highest to lowest in free energy of binding; the more negative scores imply better binding. The top-ranked compounds were 34, 28, and 44, with gscores of −12.431, −11.502, and −10.62 kcal/mol, respectively. The gscores of these compounds exceeded the value of the reference inhibitor SRI-9662, which had a binding energy of −10.432  The 2D and 3D structures of the co-crystallized inhibitor SRI-9662 revealed that the 4-amino group in the 5-deazapteridine ring interacted with Val115 and Ile7 backbones through hydrogen bonds (H-bond). Moreover, N-1 and N-8 in the ring are also bound through H-bonds to Glu30 with and without a water bridge, respectively. The pyrimidine ring involved in π-π stacking interaction with Phe34, as well as the rest of the molecule, are hydrophobically bound with the nearby hydrophobic residues in the hDHFR active site (Figure 7). The 2D and 3D structures of the co-crystallized inhibitor SRI-9662 revealed that the 4amino group in the 5-deazapteridine ring interacted with Val115 and Ile7 backbones through hydrogen bonds (H-bond). Moreover, N-1 and N-8 in the ring are also bound through Hbonds to Glu30 with and without a water bridge, respectively. The pyrimidine ring involved in π-π stacking interaction with Phe34, as well as the rest of the molecule, are hydrophobically bound with the nearby hydrophobic residues in the hDHFR active site (Figure 7). For compound 34, there was an array of H-bonding interactions between the carbonyl, hydroxyl, and carboxylic groups in the molecule and amino acid residues such as Asn64, Pro66, Gln35, Lys68, and Arg70. Water molecules were involved as well in such interactions  The structure of compound 28 was very similar to that of compound 34; hence, it involved similar binding interactions with the active site residues as well as water molecules ( Figure 9). The structure of compound 28 was very similar to that of compound 34; hence, it involved similar binding interactions with the active site residues as well as water molecules ( Figure 9). The structure of compound 28 was very similar to that of compound 34; hence, it in volved similar binding interactions with the active site residues as well as water molecules ( Figure 9).  The carboxylic acid was replaced by a phenyl ring that formed a π-π interaction with Phe31 in compound 44 ( Figure 10). The carboxylic acid was replaced by a phenyl ring that formed a π-π interaction with Phe31 in compound 44 ( Figure 10).

In silico ADMET Properties
The Maestro's QikProp module in Schrödinger was applied to predict the drug-likeness ADME properties, and toxicity (ADMET) of the metabolites under investigation [30]. Table 3

In silico ADMET Properties
The Maestro's QikProp module in Schrödinger was applied to predict the druglikeness, ADME properties, and toxicity (ADMET) of the metabolites under investigation [30]. Table 3 displayed the ADMET properties and other descriptors. In general, most of the predicted properties for the compounds were within the recommended ranges. However, some compounds fell beyond the recommended ranges of certain descriptors. The logP (QPlogPo/w) and binding to human serum albumin (QPlogKhsa) were high for compounds 66, 67, 68, 69, and 70. Additionally, the toxicity was evaluated in terms of the number of reactive functional groups (#rtvFG) and HERG K + channel inhibition (QPlogHERG). It was observed that none of the compounds exceeded the acceptable range (0-2) of reactive groups. However, 11 compounds (35, 36, 37, 43, 44, 56, 66, 67, 68, 69, and 70) were predicted to inhibit the HERG K + channel. The results suggested that the high lipophilicity of these compounds is a factor that contributes to HERG inhibition and plasma protein binding. Moreover, no CNS activity was predicted for any compound.  Abbreviations: molecular weight (mol_MW), drug-likeness (#Stars), total solvent accessible surface area (SASA), number of hydrogen bond donors and acceptors (donorHB and acceptHB), predicted octanol-water partitioning (QPlogPo/w), estimated binding to human serum albumin (QPlogKhsa), number of the possible metabolites (# metab), predicted blood-brain partitioning (QPlogBB), percentage of human oral absorption, predicted IC 50 for inhibiting HERG-K + channels (QPogHERG), predicted apparent Caco-2 cell permeability in nm/s for gut-blood barrier (QPPCaco), central nervous system activity (CNS), number of reactive functional groups present (#rtvFG), and percent human oral absorption.

MD Simulation
We performed MD simulation for the top three scoring compounds from the docking study (34, 28, and 44), as well as to the co-crystallized inhibitor SRI-9662, using Desmond software in Schrödinger [31,32]. The RMSD of proteins (blue) and ligands (red) are presented at the left and right Y-axes of the plot, respectively.
The results showed that the hDHFR protein and the co-crystallized inhibitor SRI-9662 were stable during the 100 ns of the simulation run time ( Figure 11A). The fluctuations were insignificant and lay within the acceptable range of 1-3 Å (the differences were within 1 and 1.8 Å for the protein and ligand, respectively). This confers a high-binding and stable protein-ligand complex throughout the run.

MD Simulation
We performed MD simulation for the top three scoring compounds from the docking study (34, 28, and 44), as well as to the co-crystallized inhibitor SRI-9662, using Desmond software in Schrödinger [31,32]. The RMSD of proteins (blue) and ligands (red) are presented at the left and right Y-axes of the plot, respectively.
The results showed that the hDHFR protein and the co-crystallized inhibitor SRI-9662 were stable during the 100 ns of the simulation run time ( Figure 11A). The fluctuations were insignificant and lay within the acceptable range of 1-3 Å (the differences were within 1 and 1.8 Å for the protein and ligand, respectively). This confers a high-binding and stable protein-ligand complex throughout the run.  In the case of compounds 34, the RMSD plot analysis revealed a stable protein-ligand complex, and the fluctuations of both the protein and ligand RMSD charts were within the acceptable range as well ( Figure 12A). A small jump in the RMSD of the protein between 55 and 75 ns was observed and then resumed to its normal level until the end of the run. This change was within the acceptable range. In the case of compounds 34, the RMSD plot analysis revealed a stable protein-ligand complex, and the fluctuations of both the protein and ligand RMSD charts were within the acceptable range as well ( Figure 12A). A small jump in the RMSD of the protein between 55 and 75 ns was observed and then resumed to its normal level until the end of the run. This change was within the acceptable range. The RMSD of the protein during its interaction with compound 28 was also within range; however, the compound was less stable inside the binding pocket as the RMSD fluctuated slightly beyond the 1-3 Å range ( Figure 13A). A possible explanation would be the absence of the 2-hydroxyethyl group, which was the only difference between the 28 and 34 structures (Figures 8A and 9). The chirality of the OH group and adjacent carbon atom may stabilize the carboxylic acid side chain by restricting its free rotation, leading to a decrease in the number of structural conformations available for binding to the protein active site. The RMSD plot analysis of compound 44 was comparable to 34 in terms of complex stability and RMSD fluctuation range ( Figure 14A). The RMSD of the protein during its interaction with compound 28 was also within range; however, the compound was less stable inside the binding pocket as the RMSD fluctuated slightly beyond the 1-3 Å range ( Figure 13A). A possible explanation would be the absence of the 2-hydroxyethyl group, which was the only difference between the 28 and 34 structures (Figures 8A and 9). The chirality of the OH group and adjacent carbon atom may stabilize the carboxylic acid side chain by restricting its free rotation, leading to a decrease in the number of structural conformations available for binding to the protein active site. The RMSD plot analysis of compound 44 was comparable to 34 in terms of complex stability and RMSD fluctuation range ( Figure 14A).  The secondary structure elements (SSE), alpha helices in orange and beta strands in light blue, of hDHFR (PDB: 1KMV) complexed with each ligand, were monitored during the 100 ns simulation time. The results showed that the total %SSE was maintained for all ligands and that most of the individual SSEs (by residue index) were stable throughout the run. Minor differences were observed in the turn and loop regions between the reference inhibitor SRI-9662 and the three investigated metabolites (Figures 11, 12B, 13B and 14B), as a few parts of these regions were less stable when each of the metabolites interacted with the protein.
Additionally, the MD analysis displayed the specific contact points between the inhibitor and the amino acid residues in the enzyme active site. Figure 15A  The secondary structure elements (SSE), alpha helices in orange and beta strands in light blue, of hDHFR (PDB: 1KMV) complexed with each ligand, were monitored during the 100 ns simulation time. The results showed that the total %SSE was maintained for all The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDHFR-SRI-9662 interactions presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed.
However, the strongest binding interaction was observed with Glu30, which involved multiple H-bonds with the 2-amino, N-1, and N-8 of the 5-deazapteridine ring either directly or through a water bridge ( Figure 15B). The combined effect of these interactions was maintained for about 240% of the run time ( Figure 15A). An H-bond through a water bridge was also formed with Trp24. Additional significant hydrophobic interactions with values near or above 1 were observed with Phe31 and Phe34, as well as, to a lesser extent, with Pro61. Detailed ligand atom interactions that lasted more than 30% of the simulation time with selected amino acid residues are presented in Figure 16B. The The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDHFR-SRI-9662 interactions presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed.
However, the strongest binding interaction was observed with Glu30, which involved multiple H-bonds with the 2-amino, N-1, and N-8 of the 5-deazapteridine ring either directly or through a water bridge ( Figure 15B). The combined effect of these interactions was maintained for about 240% of the run time ( Figure 15A). An H-bond through a water bridge was also formed with Trp24. Additional significant hydrophobic interactions with values near or above 1 were observed with Phe31 and Phe34, as well as, to a lesser extent, with Pro61. Detailed ligand atom interactions that lasted more than 30% of the simulation time with selected amino acid residues are presented in Figure 16B. The dimethoxyphenyl ring formed a π-π stacking interaction with Phe31, and the amino groups and nitrogen atoms of the 5-deazapteridine ring H-bonded with several residues, as detailed above. Moreover, the total specific interactions were presented in the top panel of Figure 16C, while the bottom plot showed the ligand-protein interactions by residues in each trajectory frame. The darker and more continuous the orange color is, the stronger the binding, which was observed with the key residues of Ile7, Trp24, Glu30, Phe31, Phe34, and Val115.
Molecules 2022, 27, x FOR PEER REVIEW 22 of dimethoxyphenyl ring formed a π-π stacking interaction with Phe31, and the amin groups and nitrogen atoms of the 5-deazapteridine ring H-bonded with several residue as detailed above. Moreover, the total specific interactions were presented in the top pan of Figure 16C, while the bottom plot showed the ligand-protein interactions by residu in each trajectory frame. The darker and more continuous the orange color is, the strong the binding, which was observed with the key residues of Ile7, Trp24, Glu30, Phe3 Phe34, and Val115. The coordination that appears at least 30 of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe mo than 100% as some residues might be linked to the same ligand atom more than once; the pi colored arrows represent the hydrogen bond interaction while the % on the arrows represent t coordination % during the simulation time. (C) Timeline of hDRFR-34 interactions presented in The top panel presents the total number of specific interactions the protein made with the liga over the course of the trajectory. The bottom panel presents the residues interacting with the liga in each trajectory frame. The dark orange color indicates more than one specific interaction is ma between some residues and the ligand. As the orange color becomes lighter, the number of conta between the residue and ligand decreases. The white color indicates no interaction is formed.
It was observed that compound 34 was bound to a different set of amino acid re dues in the active site of hDHFR compared to the native inhibitor. The differences cou be related to the polar nature of 34. The major types of interactions between 34 and t active residues were the direct and water-aided H-bonds with Gln35, Asn64, Lys68, an Arg70, where the latter was the strongest (> 200%) and involved in additional ionic inte action with the carboxylate anion in the ligand ( Figure 16A-C). Minor hydrophobic co tacts were formed with Phe31, Phe34, Ile60, and Leu67 with < 50% ( Figure 16A). The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDRFR-34 interactions presented in A. The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed.
It was observed that compound 34 was bound to a different set of amino acid residues in the active site of hDHFR compared to the native inhibitor. The differences could be related to the polar nature of 34. The major types of interactions between 34 and the active residues were the direct and water-aided H-bonds with Gln35, Asn64, Lys68, and Arg70, where the latter was the strongest (> 200%) and involved in additional ionic interaction with the carboxylate anion in the ligand ( Figure 16A-C). Minor hydrophobic contacts were formed with Phe31, Phe34, Ile60, and Leu67 with <50% ( Figure 16A). Similar binding modes were observed with compound 28 and the key residues mentioned above ( Figure 17A-C), owing to the structural similarity between 28 and 34. Similar binding modes were observed with compound 28 and the key residues mentioned above ( Figure 17A-C), owing to the structural similarity between 28 and 34. The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDRFR-28 interactions presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed.
Moreover, it was found that compound 44 bound to amino acid residues of Phe31, Phe34, and Glu30 in the active site ( Figure 18A) similar to the reference inhibitor ( Figure  16A). It might be due to the hydrophobic nature of 44, in which the phenyl ring in the side chain replaced the carboxylic acid in compounds 34 and 28, which formed strong The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDRFR-28 interactions presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed.
Moreover, it was found that compound 44 bound to amino acid residues of Phe31, Phe34, and Glu30 in the active site ( Figure 18A) similar to the reference inhibitor ( Figure 16A). It might be due to the hydrophobic nature of 44, in which the phenyl ring in the side chain replaced the carboxylic acid in compounds 34 and 28, which formed strong hydrogen and ionic interactions with the protein. Compound 44 formed hydrophobic interactions with Phe31 and Phe34 and H-bonds mostly in the presence of water molecules with Gly20, Leu27, Glu30, and Ser59 ( Figure 18A,B). However, the availability of these interactions did not exceed 60% of the simulation time. The inconsistent binding affinity was also demonstrated by the low total number of specific contacts and incomplete orange lines with each residue in Figure 18C, top and bottom panels, respectively.  Figure 18A,B). However, the availability of these interactions did not exceed 60% of the simulation time. The inconsistent binding affinity was also demonstrated by the low total number of specific contacts and incomplete orange lines with each residue in Figure 18C, top and bottom panels, respectively. The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDRFR-44 interactions is presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed. The coordination that appears at least 30% of the designated time period (0.00 to 100.00 nsec) will be presented. It is possible to observe more than 100% as some residues might be linked to the same ligand atom more than once; the pink colored arrows represent the hydrogen bond interaction while the % on the arrows represent the coordination % during the simulation time. (C) Timeline of hDRFR-44 interactions is presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. As the orange color becomes lighter, the number of contacts between the residue and ligand decreases. The white color indicates no interaction is formed.

Ligand and Protein Preparation
The docking study was performed with the Schrodinger program (Schrödinger Release 2022-3: Schrödinger, LLC, New York, NY, USA, 2021). The crystal structure of hDHFR complexed with the SRI-9662 inhibitor was downloaded from the protein data bank (PDB; ID: 1KMV) [33]. The protein was prepared using the "Protein Preparation Wizard" tool in Maestro software, where the missing hydrogens were added to the residues, the metal ionization state was corrected, and the water molecules > 5 Å from protein residues were deleted. The protein was then refined by predicting the pKa of the ionizable residues using PROPKA and water molecules >3 Å (not involved in the water bridge), which were removed [34]. Finally, the protein was minimized by applying the OPLS4 force field. All ligands (sesquiterpene metabolites in addition to the co-crystallized inhibitor SRI-9662) were prepared before docking using the "LigPrep" tool, where their 2D structures were converted to 3D and energy-minimized using the OPLS3 force field [35]. The hydrogens were added, and all possible ionization states and tautomeric forms were created at a pH of 7.0 ± 0.2 by Epik; a desalt option was also chosen. The H-bonds were optimized by predicting the pKa of ionizable groups using PROPKA [34].

Grid Generation and Molecular Docking
Before docking, a grid box was generated around the active site of hDHFR (PDB: 1KMV) containing the co-crystallized inhibitor SRI-9662, aided by Glide's "Receptor-Grid-Generation" tool in the Schrödinger suite [36]. The box was built around the co-crystallized ligand by selecting the "centroid of workspace ligand" function. The length of the box in each of the X, Y, and Z dimensions was set by default at 10 Å. All compounds under investigation were docked inside the grid box once with standard precision (SP) and three times with extra precision (XP) protocols (for validation purposes), and all other parameters were set to default [37]. The non-polar atoms were set for the VdW radii scaling factor at 1.0, and the partial charge cutoff was 0.25. Finally, the "Ligand Docking" tool was implemented for docking [38]. To further validate the docking method, the co-crystallized inhibitor was re-docked inside the grid box and evaluated. The docking results were assessed in terms of the gscore (ranks different compounds), emodel (ranks different conformers), and XP gscore. Glide uses emodel scoring to select the best poses of the docked compounds; then, it ranks the best poses based on the given gscores. The XP gscore ranks the poses generated by the XP Glide mode. The XP Glide takes into consideration the major driving forces and structural motifs that contributed to protein-ligand binding affinity, as described in the following equations [37]: XP Glide Score = E coul + E vdW + E bind + E penalty E bind = E hyd_enclosure + E hb_nn_motif + E hb_cc_motif + E PI + E hb_pair + E phobic_pair E penalty = E desolv + E ligand_strain where: E is energy (calculated for each of the following descriptors); E coul is Coulomb energy, E vdW is van der Waal, E bind is the energy that favors binding, E penalty is the penalty that disfavors binding, E hyd_enclosure is hydrophobic enclosure, E hb_nn_motif is special neutralneutral hydrgen-bond motifs, E hb_cc_motif is special charged-charged hydrogen-bond motifs, E PI is pi-cation interactions, E hb_pair is hydrogen bond pair, E phobic_pair is lipophilic pair, and E desolv is desolvation energy [37].
The Molecular Mechanics Generalized Born Surface Area (MMGBSA) method was used to calculate the binding free energy (∆G) of the protein-ligand complexes in a solvent based on the following equation: ∆G binding = G complex − (G protein + G ligand ) where G complex is the free energy of the protein-ligand complex, G protein and G ligand are the free energies of unbound protein and ligand in the solvent, respectively [39,40].

ADME Properties
The ADME properties (absorption, distribution, metabolism, and excretion) and toxicity of compounds under investigation were predicted using the QikProp module of the Schrodinger suite [30]. This module is useful in predicting the physicochemical properties and other descriptors to facilitate the drug discovery and development process by identifying and eliminating non-drug-like compounds from entering the clinical stage and failing. The predicted descriptors were molecular weight (mol_MW), drug-likeness (#Stars), total solvent accessible surface area (SASA), number of hydrogen bond donors and acceptors (donorHB and acceptHB), predicted octanol-water partitioning (QPlogPo/w), estimated binding to human serum albumin (QPlogKhsa), number of the possible metabolites (# metab), predicted blood-brain partitioning (QPlogBB), percentage of human oral absorption, predicted IC 50 for inhibiting HERG-K + channels (QPogHERG), predicted apparent Caco-2 cell permeability in nm/s for gut-blood barrier (QPPCaco), central nervous system activity (CNS), number of reactive functional groups present (#rtvFG), and percent of human oral absorption. The predicted values were then compared to the recommended range derived from values determined for 95% of known drugs.

MD Simulation
MD simulations were performed using Desmond in the Schrödinger package [31,32]. The software created a simulated environment that resembles the dynamic nature of the molecular system under physiological conditions to evaluate the virtual stability of interactions between the protein and ligand [41]. The RMSD plots assess the stability of protein-ligand complexes by calculating the deviation of the protein and ligand atoms inside the binding pocket at the end of the simulation time of 100 ns and comparing the results to their initial positions before the simulation at 0 ns [42]. The hDHFR protein was first complexed with the desired ligand in the docking experiment. The proteinligand complex was then tuned through the "System-Builder" tool to generate the solvated system for simulation. The solvent model was set to TIP3P, the selected box shape was orthorhombic, and the box dimensions were 10 Å. Sodium ions (Na + ) were added to neutralize the system. The simulation parameters were set up in the "Molecular Dynamic" tool, where the protein-ligand complex was evaluated at pH 7.0 ± 0.2 over a simulation time of 100 ns. The ensemble class was set as NPT to maintain a constant temperature and pressure of 300 K and 1.01325 bar, respectively, throughout the run. The generated results were analyzed at the end of the MD simulation.

Conclusion
Our findings suggested that 34, 28, and 44 possessed potent interactions with the active site of hDHFR compared to SRI-9662 (reference standard), as demonstrated by the docking studies. The MD simulation revealed detailed information about protein-ligand complex stability and specific binding contacts between each ligand and the hDHFR binding site. Additionally, the ADMET prediction demonstrated that all physicochemical parameters and ADMET properties are within the satisfactory range described for human treatment for the majority of sesquiterpene metabolites. These metabolites could be possible leads for DHFRI candidates; however, more in vitro, in vivo, and mechanistic investigation, as well as developing semi-and synthetic derivatives to enhance their DHFRI effectiveness, should be the focus of future research.