Integrated screening of Indonesian marine natural products as anticancer candidates through ADMET-clustering analysis, molecular docking, and molecular dynamics simulation

2.1. Data collection and ADMET calculation

Initially, 32,000 compounds were found in the CMNPD (https://cmnpd.org/) [12] accessed in May 2024. All compounds were filtered based on the collection site in Indonesia, and as many as 239 natural compounds were found (Supplementary 1). The compounds were then stored in SMILES format and run for ADMET analysis using the web program pkCSM (https://biosig.lab.uq.edu.au/pkcsm/) [13].

2.2. Clustering analysis

In this study, 239 compounds were screened by calculating the similarity to approved EGFR inhibitors (afatinib, osimertinib, and erlotinib). A similarity matrix with dimensions 239 × 3 was generated using a distance formula based on all ADMET variables. To effectively analyze the mixed types of data in our dataset, the Gower distance was calculated, which is suitable for handling different data types [14]. The Gower distance was employed to quantify the similarity between the compounds and the reference drugs, with values approaching 1 indicating greater similarity and values approaching 0 indicating greater dissimilarity. To filter the compounds, a cutoff value of 0.2 was applied, resulting in 67 compounds with similarity to at least one of the reference drugs that is less than 0.2.

Subsequently, hierarchical clustering was performed using the complete linkage method to group the 67 compounds and the three drugs simultaneously based on their ADMET profiles [15]. A circular dendrogram was created to visualize these relationships and to identify the compounds most like the reference drugs. The calculations for Gower distance and hierarchical clustering were conducted using the R software, specifically utilizing the ‘cluster,’ [16] ‘dendextend,’ [17] and ‘circlize’ [18] packages.

2.3. Molecular docking

The compounds sharing properties with EGFR inhibitors (erlotinib, afatinib, and osimertinib) were selected for molecular docking analysis. The proteins used were EGFR (PDB ID: 1M17) and the mutant EGFR^{L858R/T790M/C797S} (PDB ID: 6LUD) for target docking, both of which were freely obtained from RCSB PDB (http://www.rcsb.org/-pdb). The simulation was performed on the active site of the protein, identified based on its native ligand position at the grid points of 42 × 40 × 40 (X, Y, Z) and coordinates of the central grid point of maps at 21.697, 0.303, 52.093 (for EGFR with erlotinib as native ligand). Meanwhile, for the EGFR mutant with osimertinib as native ligand sets at grid points 30 × 24 × 40 (X, Y, Z) and coordinates of the central grid point of the maps at -49.873, 0.027, -18.135. The simulation was performed with AutoDockTools 1.5.7 [19] and visualized with Discovery Studio Visualizer 2016 [20]. The earlier study provides an overview of the entire simulation process [21]. The docking parameters used have been validated by re-docking native ligands with their respective proteins and obtained RMSD values of 1.63 Å for erlotinib and 0.70 Å for osimertinib.

2.4. Molecular dynamics simulation

Molecular dynamic simulations were performed with the YASARA Dynamic (developed by Biosciences GmbH). The simulation parameters were set to physiological temperature and pH (310K, pH 7.4), and counter ions were Na⁺ and Cl^-. The MD simulation was conducted with AMBER14 forcefield for 100,000 ps (100 ns) simulation time. The molecular trajectory was recorded every 25 ps. The number of hydrogen bonds, root mean square deviation (RMSD) Cα, and radius of gyration (RG) were obtained by running the md_analyze macro. Additionally, the md_analyzeres macro was used for the root mean square fluctuation (RMSF), and the md_analyzebindenergy macro was used for the molecular mechanics poisson-boltzmann surface area (MMPBSA).

3.1. ADMET-clustering analysis

The circular phylogenetic tree in Figure 1 illustrates the clustering of 67 compounds based on their ADMET properties in relation to three reference drugs: afatinib, erlotinib, and esimertinib. These compounds are grouped according to similarities in their ADMET profiles, which are crucial for determining their drug-likeness. Notably, compounds such as C178, C100, C162, and C116 cluster closely with afatinib and osimertinib, suggesting they share similar ADMET characteristics. Similarly, compounds C131, C143, C28, C144, C83, C147, C185, C19, and C153 are closely associated with erlotinib. The compounds in these clusters can be considered potential candidates due to their similar ADMET properties to the reference drugs.

Close

Figure 1.

Clustering analysis.

Export to PPT

Conversely, several compounds are significantly distant from afatinib, erlotinib, and osimertinib, indicating dissimilar ADMET profiles. For instance, compounds in the green clusters are far removed from the reference drugs, suggesting they may possess different pharmacokinetic and toxicity characteristics. This substantial divergence implies that these compounds are less likely to be effective alternatives and may not exhibit similar therapeutic benefits.

3.2. Molecular docking

A total of 13 compounds selected based on clustering analysis were bound into EGFR and EGFR mutants through molecular docking. These compounds were C178, C100, C162, C116, C131, C143, C28, C144, C83, C147, C185, C19, and C153. The summary of the simulation has been provided in Table 1. The compounds exhibited a satisfactory interaction with both proteins, as demonstrated by the compound complex’s binding energy value of less than -4 kcal/mol [28].

Table 1 also demonstrates that C162 and C116 continuously exhibited the lowest binding energies to EGFR and EGFR mutants. In contrast to the native ligand (erlotinib), which only has a binding energy of -6.77 kcal/mol, compounds C162 and C116 have binding energies of -10.77 and -9.92 kcal/mol against EGFR, respectively. These compounds are significantly more stable. Additionally, compound C162 demonstrated the lowest binding energy value against the EGFR mutant, followed by compound C116, with values of -9.54 and -9.07 kcal/mol, respectively. Meanwhile, its native ligand (osimertinib) only has a binding energy of -8.04 kcal/mol. In addition to these two compounds, C143 was another one that demonstrated better binding energy compared to three generations of EGFR inhibitors; however, it was shown solely against EGFR and not against EGFR mutants. These complexes, which exhibit higher stability than the approved EGFR inhibitors (osimertinib, erlotinib, and afatinib), were subjected to molecular dynamics simulation analysis.

Therefore, to investigate the similarity of the interactions, in-depth observations of the interactions contributing to complexes stabilization with EGFR in 2D and 3D for compounds C116, C162, C143, and erlotinib were further reviewed. The compounds have promising binding energies, suggesting their potential as EGFR inhibitor candidates, as previously mentioned. According to Figure 2, the amino acids of LEU820, VAL702, LEU694, ALA719, and THR766 appear to be identical residues to stabilize the complex, even in different ways. Nevertheless, no conventional hydrogen bonds were observed between C116/C162 and MET769.

Close

Figure 2.

Visual representation of the 2D and 3D of complex (a) Erlotinib-EGFR, (b) C116-EGFR, (c) C162-EGFR, and (d) C143-EGFR.

Export to PPT

Further observations of interactions in complexes with EGFR mutants for compounds C116 and C162 were compared to osimertinib. These two compounds were chosen because they have better binding energy than approved EGFR inhibitors. Based on the illustration in Figure 3, the three complexes have a lot of comparable amino acids in addition to interacting with the amino acid MET793. These amino acids were SER797, GLY796, LEU718, PRO794, LEU844, LEU792, ALA743, and MET790, stabilizing the complexes in different ways.

Close

Figure 3.

Visual representation of the 2D and 3D of complex (a) Osimertinib-EGFR mutant, (b) C116-EGFR mutant, and (c) C162-EGFR mutant.

Export to PPT

Overall, the molecular docking against wild-type and mutant forms of EGFR indicated that compounds C162 and C116 were the most likely candidates to inhibit both proteins. It can be seen from the binding energy when it forms a complex with the intended protein; this information also demonstrates the compound’s affinity to the target protein. The binding energy of both compounds is consistently the lowest when forming complexes with EGFR and EGFR mutants, even lower than approved EGFR inhibitors. Stabilization of the complex with EGFR generally involves the amino acid MET769 through hydrogen bond interaction. This amino acid is one of the residues that play a role in stabilizing the structure and functional regulation of EGFR [29]. On the other hand, MET793 is the amino acid residue in the EGFR mutant that is most frequently present and forms hydrogen bonds with the compound. This amino acid is critical for stabilizing the structure of the protein, according to Ahmadi et al. [30].

3.3. Molecular dynamics simulation

Initially, the selected compound and native ligand were positioned on the protein binding site, and then the molecular trajectory was observed for 100 ns. The focus of the simulation was to observe several parameters such as energy, total H-Bond, RG, RMSD Cα, RMSF. For the EGFR protein, simulations were run on the complex of compounds C116, C143, and C162, which were then compared with the three generations of EGFR inhibitors (afatinib, erlotinib, and osimertinib), with results illustrated in Figure 4.

Close

Figure 4.

Molecular dynamics simulation of complex selected compound with EGFR (a) Energy, (b) Total H-Bond, (c) RG, (d) RMSD Cα, and (e) RMSF.

Export to PPT

The comparison of the average binding free energy generated by the complexes with EGFR during the 100 ns simulation shows that afatinib has the lowest value, followed by erlotinib, C116, osimertinib, C162, and C143. The comparison of the average binding free energy of these complexes can be seen in Figure 4(a). According to this value, compound C116 has a higher affinity than compounds C162, C143, or osimertinib, but it is still not as high as afatinib and erlotinib. In line with the results, the compound C143 has the least number of hydrogen bonds. Meanwhile, the compound C162 has the highest average of total H-bonds at 9.5 bonds, followed by C116, afatinib, erlotinib, and osimertinib at 6.2, 5, 4.6, and 3.3 bonds, respectively. The difference is clearly visible in Figure 4(b). Compared to molecular docking simulations, the total number of H-bonds in molecular dynamics simulations may appear different. This discrepancy can be attributed to fundamental variations in their approaches, particularly in handling environmental factors and molecular flexibility. The key advantage of molecular dynamics is its ability to detect H-bonds that docking simulations might miss, due to its capacity to simulate flexibility, solvation, and time-dependent interactions.

The RG analysis was conducted to determine whether the compound’s binding to the protein results in structural folding during the simulation [31]. The line chart in Figure 4(c) suggests that the protein complex containing erlotinib exhibits significant fluctuations only in the first 30 ns. Similar results were observed in the complex of compound C116. Meanwhile, the protein complexes with afatinib and osimertinib showed significant fluctuations in the first 10 ns. However, the movement of the protein structure following that was not more stable than the Erlotinib and compound C116 complexes. After the simulation ran for 50 ns, observations on protein complexes containing compounds C143 and C162 revealed stability even though their RG values remained higher.

The RMSD Cα fluctuation in Figure 4(d) demonstrates that the afatinib-EGFR complex only has significant moves in the backbone during the first 10 ns and tends to reach equilibrium afterward. Meanwhile, the protein complex with erlotinib reached equilibrium after 30 ns of simulation. The osimertinib-EGFR complex reached equilibrium after 15 ns; detailed observation revealed that this complex experienced extremely small changes compared to afatinib and erlotinib after equilibrium was reached. In contrast to the three generations of EGFR inhibitors, the compounds C116, C143, and C162 took longer to reach the equilibrium. Compound C162 tends to be stable and has RMSD close to the Osimertinib after 50 ns of simulation. At the same time, compounds C116 and C143 are stable and have RMSD similar to the afatinib complex after 50 ns of simulation.

The fluctuation of amino acid residues in the EGFR protein, as estimated by the RMSF value (Figure 4e), reveals that the complex formed by compounds C116 and C162 with EGFR possesses similarities with the erlotinib-EGFR complex and other EGFR inhibitors. The facts suggest that the amino acids comprising the protein exhibit identical flexibility. Furthermore, the same RMSF value indicates that the binding of compounds C116 and C162 does not affect the stability of the EGFR. The C143-EGFR complex exhibits slightly different behavior, with higher RMSF values per residue, but still has the same pattern.

The EGFR mutant’s structural alterations, interactions, and movement were also observed after docking with compounds C116 and C162. The compounds were selected due to their greater affinity with the protein compared to the native ligand (osimertinib). The activity of both compounds will be compared with that of osimertinib, a third-generation EGFR inhibitor, and the overall findings of the simulation have been displayed in Figure 5.

Close

Figure 5.

Molecular dynamics simulation of complex selected compound with EGFR mutant (a) Energy, (b) Total H-Bond, (c) RG, (d) RMSD Cα, and (e) RMSF.

Export to PPT

Compared to C116 and C162, osimertinib seemed to have a better affinity, according to the simulation results that lasted for 100 ns. The osimertinib complex has the lowest average free binding energy, followed by the C116 and C162 (Figure 5a). However, the compounds C116 and C162 have more average total H-bonds when compared to Osimertinib, refer to Figure 5(b). The average total H-bonds of these compounds are 7.3, 6.9, and 3.8 bonds, respectively.

Osimertinib, C116, and C162 compounds do not significantly alter the protein structure (compactness), as evidenced by the RG analysis in Figure 5(c) as indicated by extremely slight fluctuations in RG values during the simulation. Analogous occurrence observed in the RMSD Cα analysis is depicted in Figure 5(d). The data suggest that the EGFR mutant protein’s backbone is not significantly affected by the compound’s binding. Meanwhile, the amino acid residue fluctuations observed based on the RMSF graph in Figure 5(e) show that the three compounds have relatively similar patterns except at residues 975-100 for osimertinib, which has a higher RMSF value. The higher RMSF in the region indicates higher flexibility affected by osimertinib, resulting in destabilization or alteration of the conformational dynamics of the protein in that region.

The performance of the hit compounds C162 and C116 was extremely encouraging. The compound does not appear to exhibit any notable fluctuation across the different parameters examined. Based on a structural review, the compounds C116 and C162 were derived from butanolide. Butanolide possesses a wide range of potential bioactivities, including pancreatic cancer treatment [32], neuroprotective agents [33], and activity against other cancer cells [34]. Compounds C116 and C162 were Cadiolide A and B, respectively, which were first discovered in Botryllus sp [35]; the compound structures have been shown in Figure 6(a-e). The secondary metabolites exhibit antiviral and antibacterial activity [36,37]. Further investigation into this compound can be based on the perspective of other possible bioactivities based on these computational studies.

Close

Figure 6.

Two-dimensional structures of (a) Afatinib, (b) Erlotinib, (c) Osimertinib, (d) C116, and (e) C162.

Export to PPT