Structural parameter variations with method
Single crystals of commercially available compounds, including 2, 9, 10, 11, 14, 15, and 16, were obtained at room temperature for X-ray diffraction. The crystal data, data collection, and structure refinement details are described in published structure reports50,51,52,53,54,58,59,60,61,62,63. Structural parameters, such as bond lengths and angles, for these compounds, were calculated at the B3LYP level of theory using different basis sets (6-31G(d,p), 6-311G(d,p), and 6–311 + + G(d,p)). A comparison was made between the computational results and the experimental data to assess the method’s precision54,58,59,60,61,62,63,64. This analysis evaluated any discrepancies or variations between the calculated and experimental data50,51,52,53,54,58,59,60,61,62,63,64. The structures of compounds 2, 9, 10, 11, 14, 15, and 16, along with their respective bond lengths and angles obtained from the selected methods, are depicted in Fig. 1 (refer to Figs. S2 to S19 and Tables S1 to S20 in the SI). The summary of the computed and experimental values is presented in Table 1 (See Table S22 in the SI).
Upon comparing the optimized and crystal structures of the compounds, we observed that the position of the carbon six (C6) group attached to the chromone ring differed, resulting in calculated bond lengths in the chromone ring that deviated from the crystal data58,59,60,61,62,63. To investigate the impact of EDG and EWG at the C6 position of the chromone ring, we conducted ab initio calculations. Our findings demonstrated that the bond lengths and angles of the EWG series (Y = F, Cl, or Br) at the C6 position increased from fluorine to bromine, with the bromine atom stabilizing the chromone ring due to its high polarizability. In the gas phase, compounds 2 and 10 exhibited mean deviations (MDs) in bond distances ranging from 0.677 Å to 0.728 Å, while compounds 9 and 11 displayed MDs ranging from 0.812 Å to 0.863 Å. All tested basis sets at the B3LYP level of theory exhibited similar performance in predicting the bond lengths and angles of compounds 2, 9, 10, and 11. It was observed that the inclusion of diffuse functions in the basis set had minimal impact on the geometries of these compounds at the B3LYP/6–311 + + G(d,p) level of theory. Specifically, the B3LYP/6-31G(d,p) method provided good agreement with experimental results64,65,66,67,68,69,70 in estimating the bond lengths and angles of compounds 2, 9, 10, 11, 14, 15, and 16 (Tables 1 and S20 in the SI). Hence, B3LYP/6-31G(d,p) can yield acceptable results while reducing computational time.
Vibrational spectra
The FTIR spectra of 3-formylchromone (1) showed that the characteristic IR bands calculated at B3LYP/6-31G(d,p) and corresponding assignment are presented in Fig. S1 and Table S21 in the SI, which is in good agreement with experimental results71. In general, the bands around 1665–1458 cm−1 at B3LYP/6-31G(d,p) are closer to experimental values71 (around ~ 1652 to ~ 1464 cm−1 in the benzene ring are assigned to the skeletal C–C stretching modes. The calculated splitting was observed in the C‒C stretches from C‒CHO for 3-formylchromone has strong and medium bands, which were exhibited as ~ 1380 to ~ 1329 cm−1, which is closer to experimental values (~ 1360 to ~ 1340 cm−1. In the spectrum of ligand 1, the strong C=O stretching band in chromone appeared in the range 1690–1630 cm−1, and calculated splitting was observed (Fig. S1 and Table S21 in the SI). Similarly, the experimental splitting was observed in the C = O stretches in the same regions. The 1 ligand by both antisymmetric and symmetric vibrations for the C‒H stretching modes was observed at ~ 3005 to ~ 2905 cm−1, which is closer to experimental values (~ 2972 to ~ 2900 cm−1). There is good agreement between the experimental X-ray data71 and theoretical data (Fig. S1 of the SI).
Analysis of Frontier molecular orbitals
Investigating the FMO in the 3-formyl chromone derivatives (1–16) offers vital insights into their stability and reactivity64. The energy gap (Egap) between the HOMO and LUMO known as the HOMO–LUMO energy gap, is a critical factor that determines their chemical reactivity, hardness, softness, chemical potential, and electrophilic index. A narrow Egap indicates softness, signifying high reactivity but low stability, while a wide Egap indicates high stability and low reactivity. The energy levels of the HOMO and LUMO orbitals help in understanding the compounds electron-donating and accepting properties, respectively. Our study presents the findings of FMO analysis, which are summarized in Table 2 and illustrated in Fig. 2 (see Figs. S2 to S19 in the Supplementary Information), providing a comprehensive understanding of the FMOs of the compounds. We have calculated the DFT calculations at 298.15 K. This temperature is close to room temperature and represents conditions relevant to biological systems. Many biological processes, including those targeted by potential drugs, occur at or near body temperature. By performing calculations at 298.15 K, the results become more relevant to understanding how the molecules might behave within a living organism.
(a) HOMO and (b) LUMO with Egap; (c) Maps of electrostatic potential (0.02 electrons bohr−3 surface) (red = electron-rich, blue = electron-deficient) for compound 4. Regions of respective colours indicating electrophilic (blue) and nucleophilic (red) sites and partial nucleophilic sites (yellow regions); (d) NBO charges for compound 4.
The FMO analysis of 3-formyl chromone derivatives (1–16) reveals essential insights into their stability and reactivity. Among the 6-substituted derivatives (9, 10, and 11), the HOMO–LUMO gap decreases in the order F > Cl > Br, ranging from 4.557 eV to 4.542 eV (Table 2). Soft ligands (Y = Br) exhibit a smaller HOMO–LUMO gap due to the greater polarizability of their valence electrons compared to fluorine (Table 2). This enhanced polarizability contributes to the electron cloud distortion in C6-halide substituted 3-formyl chromone derivatives. The soft base (Br), characterized by high polarizability and low electronegativity, and further decreases the HOMO–LUMO gap in C6-X substituted chromone rings. Conversely, the hard base (F) with higher electron density exhibits a larger HOMO–LUMO gap. A smaller gap indicates softness, associated with high reactivity and low stability, while a wider gap implies high stability and low reactivity. The energy levels of the HOMO and LUMO orbitals provide insights into the compounds’ electron acceptor and donor properties, respectively.
For the 6-hydroxy substituted chromone derivatives, the HOMO–LUMO gap of compound 5 (6-hydroxy) is 4.573 eV, lower than that of compound 4 (6-OCH2CH3, 4.505 eV) and compound 3 (6-NH2, 4.233 eV), due to the electronic effects of the electron-rich chromone ring in 3-formyl chromone derivatives. Additionally, the presence of C6 alkyl substituents (–CH3, –CH2CH3, –CH (CH3)2) leads to a decreasing trend in the HOMO–LUMO gap, likely due to steric effects (Table 2). Compounds 2–16 exhibit higher energy gaps than the parent compound 1, indicating changes caused by EDG and EWG at the 6-substituted position. Vitexin and myricetin display lower energy gaps than the other compounds (2–16). The reference drug, dapagliflozin (D), has a higher energy gap than the other compounds. The energy gap values in this series (compounds 1 to 16) decrease in the following order: D (5.328 eV) > 2 (4.601 eV) > 3 (4.592 eV) > 4 (4.585 eV) > 5 (4.573 eV) > 9 (4.557 eV) > 10 (4.543 eV) > 11 (4.542 eV) > 14 (4.538 eV) > 15 (4.502 eV) > 16 (4.493 eV) > 12 (4.432 eV) > 8 (4.233 eV) > vitexin (4.175 eV) > 13 (4.014 eV) > myricetin (3.957 eV) > 1 (3.874 eV).
The results indicate that compounds (2–16) have higher reactivity than the reference drug D, as evidenced by their higher electronegativity (χ), chemical potential (μ), global hardness (η), and electrophilicity index (ω), as well as their lower softness (σ). The electrophilicity index65 is a commonly used parameter for predicting biological activity and identifying reactive sites by measuring the energy lowering resulting from electron transfer between the HOMO and LUMO. Meanwhile, the lower electrophilicity index of compound 8 (3.428 eV) and reference drug D (1.939 eV) provides a basis for further analysis of its potential biological activity through molecular docking with a suitable protein. The higher dipole moment of compounds (2–16) than the reference drug suggested that they may have a better binding affinity. These findings make our compounds promising candidates for further investigation of their biological activity through molecular docking studies.
The DFT results revealed that compound 4 exhibited a strong binding affinity for insulin-degrading enzyme (IDE) with a binding energy of − 8.5 kcal mol−1. Moreover, all docking values obtained for the ligands were superior to the reference drug dapagliflozin (IDE: − 7.9 kcal mol−1), indicating the potential of the ligands for diabetes treatment. Compound 4 also possesses a smaller HOMO–LUMO gap (4.585 eV) compared to Dapagliflozin (5.328 eV). A smaller HOMO–LUMO gap can be associated with greater ease of electron transfer, potentially contributing to stronger binding interactions with the target protein. For instance, compound 4 demonstrated a Homo–Lumo gap of 4.585, electronegativity (χ) of 4.128, and electrophilicity (ω) of 3.716, which are indicative of its favorable molecular characteristics. These properties align well with its strong binding affinity observed in the DFT calculations. Similarly, while compounds like Vitexin and the parent compound exhibited slightly lower binding affinities compared to compound 4, their respective molecular properties, such as Homo–Lumo gap and electronegativity, might provide insights into their interaction mechanisms with the receptor.
Furthermore, comparing the molecular properties of compound 4 with other compounds, such as myricetin and dapagliflozin, offers additional insights. Myricetin, with a binding affinity similar to compound 4, demonstrated comparable molecular characteristics, suggesting potential similarities in their interaction mechanisms with the receptor. Conversely, dapagliflozin, despite its lower binding affinity, exhibited distinctive molecular properties that may contribute to its interaction behavior with IDE.
By correlating the DFT results with the docking outcomes, we aim to provide a more holistic understanding of the ligand-receptor interactions and their underlying molecular mechanisms. This integration will enhance the interpretation of our findings and guide future studies aimed at optimizing ligand design for diabetes treatment.
NBO analysis can provide insights into the bonding and stability of molecules by examining the interactions between atoms and the distribution of atomic charges. The stabilization energies of the studied molecules are affected by intramolecular charge transfer, as shown in Figs. S2 to S19 in the SI. MEP maps in Fig. 2 (Figs. S2 to S19 in the SI) visually represent the electrostatic potential of the molecules and may be used to define areas of electrophilicity and nucleophilicity. The colors used in the MEP maps (blue for electrophilic regions and red for nucleophilic regions) indicate the strength of the electrical potential.
Electron-donating groups (EDGs) and electron-withdrawing groups (EWGs) can significantly influence the biological properties of a compound, particularly for small molecules that interact with biological targets like enzymes or receptors. EDGs tend to increase the electron density around a specific atom, making it more willing to accept protons (H +) and enhance the basicity of the compound, potentially affecting its interaction with acidic residues in biological targets. Depending on the functional group, EDGs can improve the water solubility of a compound by creating hydrophilic interactions with water molecules. This can be crucial for drug delivery and absorption in the body. In some cases, EDGs can hinder the metabolism of a compound by enzymes, potentially extending its duration of action in the body. EWGs pull electron density away from a specific atom, making it more acidic (prone to releasing protons). So that EWGs can influence the interaction with basic residues in biological targets. Depending on the functional group, EWGs can sometimes enhance the ability of a compound to pass through cell membranes, which is important for drug delivery57,72,73,74.
The compound 4 has isopropyl group can slightly increase the electron density around the chromone core, potentially influencing its interaction with biological targets. In some antibiotics, the presence of an EDG like an amine group can enhance its hydrogen bonding with bacterial targets, leading to stronger binding and antibacterial activity75,76. The position of the EDG or EWG relative to the functional group involved in biological activity can significantly impact the effect. The overall balance of EDGs and EWGs within a molecule determines the net effect on its biological properties. Understanding these effects is crucial for medicinal chemistry, where researchers aim to design drugs with optimal potency, selectivity, and pharmacokinetic properties.
Analysis of physicochemical and pharmacokinetic properties
To determine compliance with Lipinski and Veber’s criteria, it is essential to evaluate the physicochemical characteristics of 3-formyl chromone derivatives (1–16). Before a chemical can be administered orally, it must meet the five criteria as per Lipinski’s recommendations: (a) molecular weight (MW) of 500 g/mol or less; (b) octanol–water partition coefficient (log P ≤ 5); (c) no more than five H-bond donors (HBD); no more than ten H-bond acceptors (HBA); and (e) no more than 140 Å2 of topological polar surface area (TPSA)66,67. Additionally, Veber has proposed two more requirements for drug bioavailability: (a) TPSA must be less than or equal to 140 Å2 (as per Lipinski’s standards), and (b) the number of rotatable bonds (nrotb) in the molecules should be less than 10. The compounds (1–16) with the most promising biological activity were analyzed using SwissADME to assess their compliance with the Lipinski and Veber criteria. All compounds (1–16) were found to meet the prescribed limit ranges specified by the Lipinski and Veber criteria (Table 3). Additionally, there was a strong agreement (less than 6) for bioactive compounds regarding molecular weight (< 500 g/mol), MLOGP (4.15), and Log S (ESOL) within the specified limit ranges. Moreover, all compounds (1–16) that obtained a drug-like (bioavailability) score of 1 demonstrated their adherence to the criteria, providing robust theoretical support for developing innovative novel drugs.
The drug-likeness parameters calculated using molinspiration for all 16 compounds are presented in Table 4, which includes their potential as G protein-coupled receptor (GPCR) ligands, ion channel modulators (ICM), kinase inhibitors (KI), nuclear receptor ligands (NRL), protease inhibitors (PI), and enzyme inhibitors (EI). Among these, compound 4 showed a GPCR value of − 0.66, and compared to 0.15 for the reference drug D (dapagliflozin).
In drug development, it is crucial to assess the pharmacokinetic properties, such as absorption, distribution, metabolism, excretion, and toxicity (ADMET), to ensure innovative drugs’ efficient and economical creation. This study utilized SwissADME (http://www.swissadme.ch/index.php) and admetSAR (http://lmmd.ecust.edu.cn/admetsar2/) software to evaluate the ADMET properties of all 16 compounds (1–16). The evaluation involved seven essential ADMET characteristics, as listed in Table 5, including cytochrome P450 enzymes (CYP3A4 and CYP2C19) inhibition, hERG inhibition, plasma protein binding (PPB), blood–brain barrier (BBB) penetration, human intestinal absorption (HIA), and synthetic accessibility (SA) score. It should be noted that drugs that affect the central nervous system (CNS) should have good blood–brain barrier penetration, whereas those that do should not penetrate the BBB. Low absorption is defined as < 0.1, medium absorption as 0.1–2, and high BBB penetration as > 267. The results indicate that all compounds (1–16) meet the ADMET standards for drug-likeness (bioavailability), promising for developing novel drugs.
Based on our study, most compounds showed BBB penetration values between − 0.5 and − 0.7, similar to vitexin (− 0.70) and myricetin (− 0.77). The range of pIC50 predictions for potential risk from hERG activity inhibitors is between − 0.4 and − 0.8268. Our research showed that compound 4 had the lowest hERG pIC50 value of − 0.4744, while all other compounds had values below the reference range. To assess the synthetic accessibility (SA) score of the drug-like compounds (1–16), which ranges from 1 (very simple) to 10 (very difficult), we employed a unique method64. The SA scores of all the compounds (1–16) ranged from 2.4 to 2.6, which was lower than the reference drug, with a high SA score of 4.52. Vitexin and myricetin have higher SA scores than compounds (1–16). Our findings suggest that all compounds (1–16) meet the ADMET standards, indicating drug-likeness (bioavailability) and potential for developing novel medications.
Analysis of pharmacological activities
For a comprehensive investigation of the potential pharmacological effects of the 3-formyl chromone derivatives (1–16), Multilevel Neighborhoods of Atoms (MNA) descriptors were applied using the Prediction of Activity Spectra for Substances (PASS) method. The use of MNA descriptors allowed for a unique and detailed characterization of the chemical structures, which aided in clarifying the compounds’ potential biological functions. The PASS method can simultaneously predict various biological activities, such as mutagenicity, carcinogenicity, teratogenicity, embryotoxicity, and primary and side pharmacological effects. A compound’s biological activity is influenced by its structural and physicochemical characteristics, the biological entity (such as species, gender, age, etc.), and the treatment approach (such as dose, route of administration, etc.). The MNA descriptors are used by PASS to determine the probable activity (Pa) and probable inactivity (Pi) for the anticipated activity spectrum of a drug. These probabilities range from 0.000 to 1.000, where Pa + Pi ≠ 1. PASS predictions can be interpreted in several ways. For example, if Pa > 0.7, there is a high likelihood of finding the activity experimentally. If Pa is between 0.5 and 0.7, the chemical is likely to demonstrate the activity in an experimental setting, but it is likely different from recognized pharmaceutical drugs. The likelihood of finding the activity experimentally is lower if Pa < 0.5.
3-formyl chromone derivatives (1–16), except some of them exhibited Pa values greater than 0.5, as presented in Table 6. The histidine kinase inhibitor values ranged from Pa = 0.614 to Pa = 0.830, nearly identical to those of the vitexin (Pa = 0.819) and myricetin (Pa = 0.892). But the reference drug, dapagliflozin, has a lower Pa (0.439) value as a histidine kinase inhibitor. Compound 4 had a Pa value of 0.669 and a Pi value of 0.010 as a histidine kinase inhibitor, indicating that it was good in promoting antibacterial activity by targeting bacterial histidine kinase, although more potent than dapagliflozin. Dysregulated kinase activity is frequently associated with multiple disorders, including cancer. Kinases are enzymes that are essential for many biological functions, including cell proliferation, differentiation, and signaling. Moreover, all compounds, except for compounds 9 and 13, showed insulysin inhibitor values between Pa = 0.515 to Pa = 0.729, indicating that these compounds could potentially be therapeutic agents for diabetes treatment. Using a clearance mechanism anchored on degradation, IDE controls the amount of circulating insulin in several organs69,70. The reference drug shows no activity as an insulysin inhibitor. HIF1α expression is frequently increased in cancer cells and is linked to metastasis and tumor formation. Compound 4 was found to have a Pa value of 0.738 and a Pi value of 0.016 as a HIF1α expression inhibitor, indicating moderate potency in inhibiting the expression of HIF1α. Additionally, some compounds showed stronger activity as alcohol dehydrogenase inhibitors than the reference drugs.
Our findings revealed that the Pa values of the studied compounds, acting as both insulin inhibitors (anti-diabetic) and alcohol dehydrogenase inhibitors, ranged from 0.5 to 0.7, with only a few exceptions. This indicates a strong potential for these compounds to inhibit both insulin and alcohol dehydrogenase enzymes effectively. Additionally, our compounds exhibited apoptosis agonist activity, which is well-known for inducing cancer cell death71.
In silico molecular docking
Molecular docking investigation and in silico docking was performed on 3-formyl chromone compounds (1–16) and reference standards of vitexin to inhibit targeted proteins, including CAD (3TWO), BHK (3DGE), IDE (6BF8), HIF-α (2WA4), p53 (7EAX), COX (6Y3C), and Mpro of SARS-CoV2 (6LU7), as presented in Table 7. Results revealed that compound 4 had a strong binding affinity for IDE (− 8.5 kcal mol−1). Additionally, all docking values obtained were superior to the reference drug dapagliflozin (IDE: − 7.9 kcal mol−1), indicating their potential for diabetes treatment. Furthermore, compounds 1–13 showed good binding affinity, while compounds 14, 15, and 16, which had additional changes in their chemical structure, exhibited weaker binding affinity. Although Vitexin has a lower binding affinity (− 8.3 kcal mol−1) than compound 4, myricetin has a binding affinity (− 8.5 kcal mol−1) similar to compound 4.
The findings of the molecular docking study indicated that compounds 1–16 exhibited a higher binding affinity for IDE, BHK, and COX, while their binding affinity for human CAD, p53, Mpro, and HIF-α was relatively lower. These results suggest that these compounds have the potential to be used as therapeutic agents for diseases related to the dysregulation of IDE, BHK, and COX. While molecular docking is a valuable tool for drug discovery, it is crucial to consider a compound’s pharmacokinetics, toxicity, and metabolism, in addition to its binding affinity. In vitro and in vivo research is necessary to comprehend how compounds behave in living organisms and evaluate their safety and efficacy. By employing a combination of computational and experimental techniques, researchers can optimize the compounds and create dependable and effective therapeutic drugs.
D (Reference drug: dapagliflozin).
Furthermore, Figs. 3 and 4 present 2D diagrams for compound 4 in IDE and COX, depicting the ligand’s interactions with the protein and how it affects the pathogens’ active components. The bond lengths and residue numbers are indicated in Fig. 3d and Fig. 4d; the ligand–protein interaction of compound 4 through the amino acid residues of IDE is demonstrated, revealing around eight different bonds. Four of these bonds are hydrophobic, including Pi–Pi stacking on residue TYR A: 269, Pi-alkyl on residue LEU A: 156, and alkyl on residue ALA A: 434. Three conventional hydrogen bonds are at residues THR A: 271, ARG A: 432, and THR A: 163. The last bond is Pi-anion on residue GLU A: 160. The protein–ligand interaction of compound 4 in COX, shown in Fig. 4, exhibits five different interactions. Four of them are hydrogen bonds, including a conventional hydrogen bond at residue HIS A: 207, a carbon-hydrogen bond with THR A: 206, and a Pi-donor hydrogen bond with HIS A: 388. One bond is hydrophobic, including amide-Pi stacking with ALA A: 202.
Molecular docking poses: (a) Ligand in protein pocket; (b) Active site; (c) Hydrogen bonding in solid; (d) Ligand–protein interaction for 2D diagram of compound 4 and modeled protein (IDE).
Molecular docking poses: (a) Ligand in protein pocket; (b) Active site; (c) Hydrogen bonding in solid; (d) Ligand–protein interaction for 2D diagram of compound 4 and modeled protein (COX).
Based on the results of molecular docking and interaction analysis, it has been concluded that compound 4 exhibits a stronger binding affinity towards the IDE protein. Additionally, all compounds also show good interactions with the IDE protein. The 2D diagram and bond interactions of the IDE protein with all compounds are shown in SI (Table S2). Thus, the in silico study suggests that compound 4 has promising potential as an insulin-degrading enzyme protein inhibitor. The selection of Insulin Degrading Enzyme (IDE) as the target protein for our in-silico study of potential anti-diabetic agents is well-supported by existing research. IDE plays a critical role in insulin degradation, and its inhibition has been demonstrated as a promising strategy for improving blood sugar control in diabetic patients. This is further evidenced by the existence of previous in vivo and in vitro studies exploring the therapeutic potential of IDE inhibitors17. Our study aimed to leverage computational methods to identify novel candidate molecules that could target IDE and potentially offer new avenues for developing anti-diabetic drugs.
In silico molecular dynamics
The goal of the MD simulations was to examine the stability and interactions of compound 4 with the IDE and COX proteins for a duration of 20 ns. The molecular docking analysis had previously shown that compound 4 had a strong binding affinity of − 8.5 kcal mol−1 with the IDE protein and a favorable binding affinity of − 8.1 kcal mol−1 with COX, indicating its potential as a drug candidate. The MD simulation results were evaluated using RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuation), and Rg values, potential energies, temperature, and hydrogen bonding to gain a comprehensive understanding of the system’s behavior over time.
The RMSD values derived from the MD simulations offer insights into ligand–protein complexes’ stability and conformational changes, as presented in Fig. 5. For the ligand-IDE protein complex (green curve), the RMSD values were in the range of 0.2 to 0.5 nm, while for the ligand-COX protein complex (blue curve), the values were between 0.1 and 0.2 nm, indicating that both complexes are stable and exhibit minimal deviation from their initial positions during the simulation. However, the RMSD values for water and ions in the ligand-IDE protein complex (yellow curve) were approximately 7.7 nm, while those in the ligand-COX protein complex (brown curve) were about 7.3 nm. These observations suggest potential differences in the stability and dynamics of the solvent molecules in the two complexes, which could be due to the size and shape of the binding sites, as well as the specificity and strength of ligand–protein interactions. It is important to note that the RMSD values for water and ions can be affected by various simulation parameters, including force field, simulation time, and binding site definition. Therefore, further analyses, such as solvent density profiles, hydrogen bonding patterns, and residence times, are recommended to elucidate the underlying factors contributing to the differences in RMSD values.
RMSD evolution for merge docked protein–ligand complex between two proteins and modeled ligand 4, (IDE (6BF8), green line and COX (6Y3C), (4, red and black), water and ions in the ligand-IDE protein complex (yellow) and water and ions in the ligand-COX protein complex (brown) during 20 ns MD simulation.
The RMSF is an important measure of the flexibility and mobility of protein–ligand complexes in molecular dynamics simulations. In this study, separate RMSF values were calculated for the protein and ligand in the IDE (6BF8) and COX (6Y3C) protein–ligand complexes with compound 4, and the results are shown in Fig. 6. A lower RMSF value indicates greater rigidity of a particular residue or atom, whereas a higher RMSF value indicates greater flexibility or mobility. The COX protein–ligand complex had an RMSF value of 0.1 nm (red line), indicating that it is relatively rigidly held in place, while the IDE protein–ligand complex had an RMSF value of 0.3 nm (black line), indicating that it is more flexible. A fluctuation observed in the IDE protein–ligand complex at the 15,000 range amino acid atom with a value of around 1.5 nm could be due to a number of factors, such as the inherent flexibility of the amino acid or its interaction with the ligand or solvent molecules. Additional analysis, such as examining the specific interactions of this residue with other parts of the protein or ligand, may provide further insights into the cause of this fluctuation. The differences in RMSF values between the COX and IDE protein–ligand complexes suggest that they have varying degrees of flexibility or mobility, which may have implications for their biological function. The fluctuation observed in the IDE protein–ligand complex at atom 15,000 highlights the importance of careful analysis of simulation results and emphasizes the potential complexity of protein–ligand dynamics.
RMSF evolution for the IDE protein–ligand complex (black) and the COX protein–ligand complex (red) and modeled ligand 4, during 20 ns MD simulation.
The radius of gyration (Rg) is a measure of molecular compactness and is often used to monitor conformational changes in MD simulations. In this study, the Rg values for compound 4 bound to the IDE and COX proteins, were compared over a 20 ns simulation period, as shown in Fig. 7. It was observed that the Rg values for the IDE (protein–ligand) complex exhibited slight fluctuations within the range of 3.4–3.5 nm, while the Rg value for the COX complex remained consistently at 3.1 nm throughout the simulation period. This suggests that the COX complex is more stable and less prone to conformational changes compared to the IDE complex. Additionally, the Rg value for the ligand alone was found to be 0.3 nm in both cases, indicating that the ligand assumes a more compact conformation when bound to the proteins. Furthermore, the Rg value for the water-ion was larger than that of the protein–ligand complex, indicating a more loosely packed structure for water-ions compared to the ligand–protein complex. Specifically, the Rg values for the water-ion were 5.72 nm (red line) and 5.485 nm (black line) in the COX and IDE protein–ligand complexes, respectively, over the simulation period. Understanding the behavior of protein–ligand complexes can be valuable for the rational design of drugs with improved binding affinity and specificity.
Radius of gyration, Rg (nm) versus time (ps) plots of the ligand (Rg), IDE protein–ligand complex (RgX), COX protein–ligand complex (RgY) and water-ions (RgZ,) with modeled ligand 4, during 20 ns MD simulation.
Hydrogen bonds play a crucial role in protein–ligand interactions, providing valuable information about binding strength and specificity. In this study, MD simulations were used to investigate the number of hydrogen bonds formed between compound 4 and the active site of two proteins, IDE and COX, as shown in Fig. 8.
Plots of the number of hydrogen bond versus time (ps) for (a) between protein complex of PDB: 6BF8 and compound 4 and (b) between protein complex of PDB: 6Y3C and compound 4 during 20 ns MD simulation.
The number of hydrogen bonds between the ligand and the active site of the IDE protein fluctuated between 0 and 2 during the 20 ns simulation period (Fig. 8a). This suggests that hydrogen bond formation between the ligand and the IDE protein needs to be consistently maintained and is influenced by many factors. Consequently, the IDE protein and ligand interaction is dynamic rather than static. On the other hand, the number of hydrogen bonds formed between the ligand and the active site of the COX protein ranged from 0 to 1 throughout the 20 ns simulation (Fig. 8b). This indicates that the hydrogen bond formation between the COX protein and the ligand is relatively stable and consistent, implying a more static and less dynamic interaction compared to the IDE protein–ligand interaction.
These findings demonstrate that multiple factors, including conformational changes, ligand movement, and the specific characteristics of the protein, influence the number of hydrogen bonds formed between a protein and a ligand. The observed dynamic nature of hydrogen bond formation in the IDE protein–ligand interaction holds significant implications for drug design and optimization. It emphasizes the importance of considering the dynamic aspects of protein–ligand interactions in rational drug design processes.
In MD simulations, monitoring the temperature of the system is crucial to ensure its stability. In this study, the IDE and COX proteins’ temperature was relatively stable, fluctuating between 298 and 302 Kelvin during the 20,000 ps simulation period (Fig. 9a,c). This indicates that the simulation was well-controlled and the system remained within an appropriate temperature range. However, the system’s potential energy, which reflects the interactions between atoms, showed fluctuations throughout the simulation period. For the IDE protein, the potential energy fluctuated between − 1.680e+06 and − 1.677e+06 kJ mol−1 (Fig. 9b), while for the COX protein, it varied between − 1.325e+06 and − 1.323e+06 kJ mol−1 (Fig. 9d). These fluctuations imply that the interactions between the atoms in the system are dynamic and continuously changing. In summary, the stable temperature observed in the study indicates a well-controlled simulation, but the fluctuation in potential energy highlights the dynamic nature of protein–ligand interactions.
The temperature and potential energy curves over the course of the 20 ns MD simulation for the potential energy graph for the protein–ligand complex of IDE (a) and COX (c); the temperature graph of the protein–ligand complex of IDE (b) and COX (d) with compound 4.
PCA analysis
Conformational principal component analysis (PCA) was performed on the simulated molecular dynamics (MD) trajectories of the IDE protein–ligand complex with compound 4 at 300 K. The aim was to identify the variability, collective motions, and changes in protein conformations observed in subsets of the primary components throughout the MD simulations. The PCA analysis was conducted using the Bio3D program55,56,57. The resulting eigenvalues versus eigenvector plots are displayed in Fig. 10. The first three eigenvectors, namely PC1, PC2, and PC3, were utilized to compare the dominant motions within the smaller trajectory subgroup. Colored dots represent the captured variance by the eigenvectors. Regarding the internal motions observed in the MD trajectory, the protein–ligand complex simulation at 300 K exhibited the most extensive variability in PC1, accounting for 26.13% of the total variance. PC2 showed a lower percentage of variation (16.3%), and PC3 showed 15.48% of the conflict. These three components accounted for 57.9% of the total variance (Fig. 10).
PCA of the IDE complex with compound 4 at 300 K MD trajectory, where red dots represent stable conformations and black dots indicate energetically unstable conformations.
Furthermore, the simulation appeared to have converged, as evidenced by the cosine content value of the eigenvectors, which was computed to be 0.32 based on the MD trajectory. The conformational PCA analysis provided insights into the dominant motions and variability within the IDE protein–ligand complex during the MD simulation at 300 K. These findings contribute to a better understanding of the dynamic behavior and conformational changes occurring in the protein–ligand complex.
