Molecular docking and dynamics analysis of halogenated imidazole chalcone as anticancer compounds

Cancer is one of the three biggest causes of death in the world. The development of new drugs against this disease is a serious study that must be carried out in order to reduce mortality and extend the life span of sufferers. The aim of this research was to develop a new anti-cancer drug based on halogenated imidazole chalcones have been conducted. A 18 The halogenated imidazole chalcone compound was designed and performed molecular docking. Potential compounds of molecular docking are then carried out molecular dynamics. The results of molecular docking show that the potential chalcones based on bond affinity and specific interactions are chalcones B5, B6, C5, and C6. Analysis of the molecular dynamics result parameters are root mean standard deviation (RMSD) complex, root mean square fluctuation (RMSF), Radius of Gyration, Protein-ligand hydrogen bonding, and complex stability with generalized born surface area (GBSA) method, where the potential chalcone compounds are chalcones B5 and B6.


Introduction
Cancer is a disease caused by abnormal and uncontrolled cell growth that can damage surrounding tissues. Cancer is a worldwide disease where there are 19.3 million new cancer cases followed by 10 million deaths. If not treated properly, it is estimated that there will be 28.4 million cancer cases annually by 2040 (Sung et al. 2021). One of the efforts to control cancer is chemotherapy, but chemotherapy becomes less effective because of mutations that continue to occur in cancer-causing proteins so that the effectiveness of commonly used cancer drugs is greatly reduced. The decrease in the effectiveness of anti-cancer drugs that have been reported includes 35-50% of lung cancer patient deaths after chemotherapy, 50-70% of cases of ovarian cancer that reappear after 1 year of chemotherapy, and 20% recurrence of acute lymphoblastic leukemia after therapy. In addition to the problem of chemotherapy resistance, it also has disadvantages, among others, that it can be toxic to non-cancerous cells due to the less specific nature of the compound (Wang et al. 2019). So that the development of drugs for cancer therapy continues to be evolving.
Chalcones are compounds with 2 aromatic rings connected by a α,β-unsaturated carbonyl skeleton and are compounds with extensive anti-cancer bioactivity (Karthikeyan et al. 2015). Based on the study of its structure and activity, the potential of chalcone as an anti-cancer compound can be increased by modifying the aromatic ring into an imidazole heterocyclic aromatic skeleton as reported by Chouiter et al. (2020) in which imidazole ringed chalcones are toxic to gastric cancer cells (AGS). Imidazole chalcones are also reported to be toxic to lung cancer cells (A549), where halogen groups in those chalcones can significantly increase their toxicity to cancer cells (Oskuei et al. 2021) . This was reinforced by the study of the structure of activity in which the heterocyclic ring of N-atoms and the presence of halogen substituents in chalcones increased the anti-cancer activity of chalcones (Do et al. 2016).
In the development of new drug compounds it is very important to know the interactions that occur between an enzyme and an inhibitor because almost 95% of drug activity arises from these interactions (Santos et al. 2016). One of the methods for determining the formation of interactions between enzymes-inhibitors is computational or in silico modeling (Song et al. 2015). Molecular docking and molecular dynamics are molecular modeling methods that are often used to study the interactions between new compounds and target macromolecules (Naqvi et al. 2019). In the development of anti-cancer drugs, the target protein is epidermal growth factor receptor (EGFR). Excess expression of the EGFR protein is responsible for the formation of cancer cells. Molecular docking studies show that the enone skeleton of the chalcone and the N atom of the imidazole are able to inhibit the active side of the EGFR protein through hydrogen bonds and the results are linear in-vitro activity tests (Al-Anazi et al. 2021;Oskuei et al. 2021). Although molecular docking is able to show linear results between interaction and activity, the use of molecular docking has the disadvantage of not being able to describe the real interaction of a complex system with the presence of pressure and temperature between the enzyme and the substrate involving water. Where in fact when the drug is in the body the drug interacts on the active side not only involves protein but also water and is influenced by physical factors in the form of temperature and pressure (Tripathi and Misra 2017).
Based on this background, 18 halogenated chalcone compounds with imidazole skeletons (Table 1) will be designed which will then be studied for their interaction with EGFR proteins using molecular docking. The results of the molecular docking study were then continued with a molecular dynamics simulation study to study the complex stability of compounds with proteins in water solvents, a system that is close to the actual state in the body. Design compounds that have a low bond affinity and main interaction with EGFR proteins and are stable through molecular dynamics studies and are used as a reference for synthesis and study of their in-vitro activity (Ahmadabadi et al. 2022).

Materials
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.

Molecular tethering
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.

Molecular dynamics
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 CHARMMjul-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 Clanion (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 R 1 = Br R 2 = CH 3 Chalcone B5 R 1 = Cl R 2 = CH 2 -C 6 H 5 Chalcone B6 R 1 = Br R 2 = CH 2 -C 6 H 5 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.

Docking parameter validation
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).

Molecular docking
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.
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.

Evaluation of molecular dynamics results
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 overlap-ping 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).
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 inter- actions 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.
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).
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.

Conclusions
The molecular tethering results of the 18 designed chalcone compounds show that the most potential compounds are chalcones B5, B6, C5, and C6 because they are able to form specific interactions with important residues of met769 and have similar hydrophobic interactions with native ligands. The results of molecular dynamics simulation show that the complex RMSD value is relatively the same where the meta-substituted chalcone has the highest RMSD is the C6 chalcone. The RMSF and radius of gyration values of the four chalcones did not have significant differences but there was an influence of substituent position seen in the hydrogen bonds formed and the stability energy of the complex. The compounds with the most hydrogen bonds and the most stable complex stability energy are those with halogen groups substituted in the meta position.
Thus, the recommended compounds from various molecular dynamics parameter analysis are chalcones B5 and B6.