Molecular docking procedure and dynamics performed using on a computer with specification Intel Xeon E5 2609 V4 with 16GB Ram, running open suse leap 15.3 operating system.

Eighteen chalcone compounds were 3D modeled with Marvin Sketch software and optimized for geometry with Gaussian 09W (Frisch et al. 2013) using the density functional theory (DFT) method with B3LIP hybrid function and set base 6–31G (d.p). Molecular tethering is performed with Auto Dock Vina 1.1.2 software (Trot and Olson 2009). Receptors in the form of EGFR proteins (PDB ID: 4HJO) are prepared by adding hydrogen atoms to incomplete amino acid residues followed by the addition of charges. Parameters used in molecular tethering using grid box size 20×20×20 with coordinates X = 30.249; Y = 5.739; Z = 56.697 and exhaustiveness 64. Validation of the docking protocol is carried out by re-tetherering erlotinib (AQ4) as a native ligand until an RMSD value of < 2Å is obtained which indicates that the resulting docking protocol has good validity. If a valid docking protocol has been obtained then the protocol is used to anchor the 18 designed chalcone compounds. 20 docking conformations were obtained from tethering and continued with analysis using Discovery Studio 2021.

The three best compounds from molecular docking were studied for their stability on the active side. Molecular dynamics method with GROMACS software - The protein topology was created using pdb2gmx with the CHARMM-jul-2021 force field, The ligand topology was created using the CGenFF web. The system added water molecules and the system charge was neutralized by changing the Cl- anion (Tomanik et al. 2020). Furthermore, the system is minimized and continued with the equilibration of volume and temperature (NVT) and pressure and temperature (NPT) respectively for 100 and 250 ns. Molecular dynamics simulations were run within 35 ns with a system temperature of 300K and a pressure of 1 atm (Al-Anazi et al. 2021). The simulation results are trajectory data (changes in molecular position every time) which will be differentiated using VMD 1.9.2 software. And trajectory plots are visualized using XMGRACE software.

Re-tethering was performed on native ligand, namely erlotinib to validate the docking protocol The docking protocol used was center X, Y, Z = 23.907, 8.951, -0.436 with a grid box size of 20×20×20 Å followed by an exhaustiveness value of 64. The re-tethering results showed an RMSD value of 1.390 Å. The RMSD value of the re-tethering result indicates a value below 2 Å which means that the docking protocol used is valid and feasible for use in re-tethering the designed molecule (Fig. 1).

Figure 1. 

Binding poses of redocked native ligand redocked native ligand shown with green colour and the experimental position native ligand shown by blue colour.

The results of molecular tethering showed that the chalcone compound substituted for the benzyl group in the imidazole skeleton, namely the chalcone (A5, A6, B5, B6, C5, and C6) had a lower bond affinity value than erlotinib. The observed bond affinity difference values of chalcones A5, A6, B5, B6, C5, and C6 are not too far apart compared to erlotinib. The not-so-distant difference in bond affinity is thermodynamically suspected that the designed ligands have the same ease of interacting with erlotinib on the active side of the EGFR protein.

In studying the results of molecular tethering, it is not only seen from the value of bond affinity but also seen from the interactions that occur on the active side. Specific interactions on the active side indicate that the designed compound interacts precisely with the active side and the presence of specific interactions also indicates the compound’s ability to close the interaction of active amino acid residues so that the catalytic process of a protein is inhibited. The active side and interaction that plays an important role in the inhibitory activity of the EGFR protein is the interaction of hydrogen bonds with the amino acid residue Met768 followed by the production of hydrophobic interactions Leu694, Val702, Ala719, Thr766 and Leu820 (Yadav et al. 2014). The results of molecular tethering showed that of the six ligands that had the lowest affinity (chalcones A5, A6, B5, B6, C5, and C6) only chalcones B5, B6, C5, and C6 formed specific interactions with the main catalytic sites in the EGFR protein. Interactions with major amino acid residues in EGFR proteins produced on the B5, B6, C5, and C6 chalcone molecules are thought to provide the same pharmacological effects as native ligands. Visualization results (Fig. 2) of interactions with chalcones that have specific interactions show similar poses. The difference in the position of the substituent halogen group does not give a difference to the bond pose. The similarity of the chalcone bond poses of the B5, B6, C5, and C6 chalcones makes it conjecture that they interact on the active side of the EGFR protein with almost the same conformation.

Figure 2. 

Molecular docking result visualization of chalcones B5, B6, C5, and C6 in EGFR protein active site.

From the description of bond affinity and its interaction chalcone B5, B6, C5, and C6 it is very interesting to further study its stability in dynamic systems using molecular dynamics simulations (Table 2). Molecular dynamics simulations are critical in drug development, as the interaction of compounds developed as drugs with proteins occurs dynamically in body systems where water is involved.

The stability formed between proteins and ligands can be observed through changes in the root mean standard deviation (RMSD) value of protein and ligand complexes. RMSD is the value of the change in position of a structure to its initial position, and is directly related to the stability of the conformation of the structure. Measurement of the RMSD value in ligand atoms is important in assessing the entire complex of proteins and ligands and the smaller the change in RMSD value indicates the better the stability of the conformation (Knapp et al. 2011).

The results of the RMSD analysis of potential chalcones from the molecular tethering results showed that the RMSD graph had a relatively similar pattern in which the RMSD value during the highest simulation was owned by the C6 chalcone with an average value of 0.47 nm, while the average RMSD values during the simulation for the B5 and B6 chalcones were 0.43 and 0.41 nm, respectively. The lowest average RMSD value during the simulation was owned by chalcone B5 with a value of 0.38. The average RMSD value of Chalcone B5, B6, and C5 is relatively close and lower than that of RMSD of C6 chalcone. If observed from the RMSD values, it can be seen that chalcone C5 as the most stable compound forms complexes with ligands. Sequentially During the simulation process, it was observed that before reaching 23 ns RMSD from each chalcone had different values, but the fluctuation patterns were relatively the same, the same relative fluctuation patterns occurred due to the similarity in structure of the chalcones. After 23 ns it was observed that the chalcone complexes B5, B6 and C6 showed a linear and overlapping graph indicating that the ligand protein complex had stabilized, but the RMSD value of the C6 chalcone complex showed fluctuations after 27 ns, in contrast to chalcones B5 and B6 which remained stable after 23 ns. At 25 ns RMSD the C5 chalcone complex showed a decrease in RMSD value but again followed the pattern of chalcones B5 and B6 after 32 ns (Fig. 3).

Figure 3. 

RMSD of protein-ligand system.

The analysis continued by analyzing the root mean square fluctuation (RMSF) value of amino acid residues. Analysis of RMSF values shows the stability of the interactions that occur during the simulation process, which is calculated at each position of the amino acid residue and is characterized by the resulting fluctuations. From all the chalcones it is observed that the RMSF pattern formed is similar (Fig. 4), with the order of fluctuations from the lowest and highest being chalcones B6, B5, C5 and C6 (Fig. 5). The residues around the active side of the EGFR protein are in residues 685–725 and 760–765. Fluctuations of amino acid residues occur due to the influence of ligands that form bonds. The higher the resulting fluctuation value shows high residual flexibility and means that the amino acid residue interaction formed has low stability (Chaudhary and Aparoy 2017). All proposed ligands produced low RMSF values at the main catalytic site of the EGFR protein, namely the amino acid residue interaction of Methionine at position 769 (Met769), indicating that the four proposed ligands have good interaction stability against the main catalytic site of the EGFR protein.

Figure 4. 

RMSF of protein-ligand system.

Figure 5. 

Protein-ligand system radius of gyration.

Analysis of radius of gyration shows the compactness of a protein ligand complex, the higher the radius of gyration value indicates that the reduced compactness of the complex formed (Lobanov et al. 2008). Reduced compactness can be caused by the opening of the active side so that the interaction between proteins and ligands is reduced. The radius of gyration analysis showed that the protein compactness due to the presence of the four ligands was the same during the simulation, but after 27 ns the radius of gyration increased for the C6 chalcone protein complex. The increased radius of gyration supports fluctuations in complex RMSD values observed at 27 ns. From the data it can be said that the C6 chalcone complex is unstable after 27 ns (Fig. 6).

Figure 6. 

Total hydrogen bond formed in protein-ligand system.

Hydrogen bonding during molecular dynamics simulation is related to the stability of the bond between proteins and ligands because hydrogen bonds are the largest contributors to stability when interacting with specific residues on the active side (Fu et al. 2018). Structurally chalcones B5, B6, C5, and C6 have 3 hydrogen bond acceptors that allow hydrogen bonding to occur. The simulation results show that chalcones B5 and B6 have 2–3 hydrogen bonds during the simulation process, while C6 chalcones show fewer hydrogen bonds (Fig. 6). This supports previous data where the C6 chalcone was less stable on the active side. The value of the hydrogen bond formed from the C5 chalcone is less than that of the C6 chalcone will still have a lower RMSD value. The lower RMSD value of the C5 chalcone is due to the fact that the C5 chalcone has more stability contributed by the interaction of hydrophobic and van der Waals forces than the hydrogen bond.

Complex stability energy calculations during simulations can be performed using the molecular mechanics generalized born surface area (MM-GBSA). This technique calculates the average of the energy changes that contribute to the process of complex formation (Wang el al., 2019). The results of the analysis of the energy of the stability of the complex are shown by Table 2. The complex stability energy during the simulation process showed that chalcone B6 is the most stable compound when interacting with the EGFR protein. The results of the MM-GBSA calculation prove and strengthen the molecular dynamics simulation data that the B6 chalcone is a chalcone that is less stable interacting with the EGFR protein. An interesting thing is seen in the C5 chalcone as the compound with the lowest complex RMSD value. The energy during the simulation of the C5 chalcone is higher than that of the halogen-substituted chalcone at the meta position, namely chalcones B5 and B6. Of the total energy-contributing components, the energy that stabilizes C5 chalcones is van deer Waals interactions where the energy due to the interaction of van deer Waals chalcone C5 interactions is the lowest compared to other chalcones. The results of the total energy analysis during the simulation showed that the halogen-substituted chalcones in the meta were higher (chalcones B5 and B6) than the substituted ones (chalcones C5 and C6). The type of halogen substituent to the complex stability energy shows the presence of a chloro group at a higher meta position than bromo (chalcones B5 and C5) while the opposite is true in chalcones whose halogens are substituted in the para position where the chloro group has a lower energy than the bromo group (chalcones B6 and C6).

Reviewing the overall aspects of molecular dynamics parameters, information was obtained that the value of complex RMSD is not so much different, but the presence of a bromo group substituted at the meta position has a higher complex RMSD (chalcone C6). The position of the halogen substituent has no significant difference in the value of RMSF and Radius of Gyration, but the influence of the substituent position is seen in the hydrogen bond formed and the stability energy of the complex. The number of hydrogen bonds formed for chalcones whose halogens are meta-substituted is greater than that of the para-substituted. The stability energy of the complex also indicates the same thing where the stability energy of the complex is lower when the halogen is substituted in the meta position. Judging from these various parameters, the recommended compounds to be studied for activity in vitro are chalcone compounds B5 and B6.