Molecular Docking and Dynamic Simulation of Erythrina fusca Lour Chemical Compounds Targeting VEGFR-2 Receptor for Anti-Liver Cancer Activity

Abstract


INTRODUCTION
Liver cancer is an aggressive tumor often occurring in the context of chronic liver disease and cirrhosis 1 .This was diagnosed among 905,677 people worldwide in 2020, with reports of 830,180 death cases 2 .Vascular Epidermal Factor Receptor-2 (VEGFR-2), with the PDB code 4ASD, plays a crucial role in promoting angiogenesis by enhancing vascular permeability, endothelial cell survival, proliferation, migration or invasion of surrounding tissue, and capillary formation.VEGFR-2 pathway blocks tyrosine kinase activation 3 .Similarly, sorafenib approved by the Food & Drug Administration (FDA) for the systemic treatment of liver cancer functions as a multi-tyrosine kinase inhibitor 4 .This anti-angiogenic agent inhibits various growth factor tyrosine kinase receptors such as VEGF, PDGF, and c-Kit, thereby blocking the Raf/MEK/ERK signal transduction pathway 5 .However, prolonged use of sorafenib can lead to drug resistance in cancer cells, suggesting the need for new treatment options.Many discoveries and developments of new anticancer agents are selective in inhibiting cancer cell proliferation pathways 6 .
Developing new drugs is costly and timeconsuming, making computer-aided drug design (CADD) crucial for rapid production and expense reduction.CADD also known as an in silico method can efficiently identify various essential anticancer compounds 7 .In this context, molecular docking a structure-based in silico method, permits the virtual screening of approved drugs, natural products, or compounds synthesized under a reasonable timeframe 8 .Molecular dynamics simulation can predict the stability of ligand-receptor complexes identified through molecular docking by producing suitable conformational states 9 .Previous research performing molecular docking and dynamics simulation on fourteen N-aryl methylaniline/chalcone hybrids synthesized against VEGFR-2 showed the inhibition of cancer cell growth 10 .
Erythrina fusca, commonly known as the cangkring plant, has potential for the treatment of liver cancer.Bioactivity exploration performed on E. fusca reported antibacterial, antimalarial, anti-estrogenic, estrogenic, anti-ischemic, reperfusion injury healing, and antiherpetic properties.Previous phytochemical screening of the aerial (above-ground) parts identified the presence of alkaloids, flavonoids, triterpenes, steroids, saponins, lactones, coumarins, reducing sugars, carotenoids, amines, and cardiac glycosides 11 .Therefore, this research aimed to conduct a computational investigation of the interactions between chemical compounds from E. fusca Lour and VEGFR-2 receptor through molecular docking and dynamics simulation.The software used for this investigation included Pymol, AutoDock Tools, Biovia Discovery Studio, and Gromacs.Additionally, the objective was to identify chemical compounds from E. fusca with significant potential and suitable conformation for antiliver cancer activity against VEGFR-2 using the in silico method, which has not been applied in previously published journals.

Molecular Docking Preparation of VEGFR-2 Macromolecule
VEGFR-2 macromolecule with the specific ID 4ASD was downloaded from the PDB website (https://www.rcsb.org/) in (*.pdb) format.This macromolecule was separated from the solvent (H2O) and ligand using PyMOL, and then saved in (*.pdb) format.Subsequently, receptor optimization was performed through hydrogen and charge addition using AutoDock Tools, followed by saving the file in (*.pdbqt) format.

Preparation of native and test ligands
Ligands were downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/), and then optimized by setting the Torsion Tree and adjusting the number of active torsions using Autodock Tools.

Docking Simulation
Redocking was performed with a natural ligand (sorafenib) after separation and optimization.In this context, the grid box was arranged as grid center x = -24.964,y = -0.497,and z = -10.399with dimensions 40x40x40 Å and spacing of 0.375 Å.Subsequently, the grid box configuration was saved in (*.gpf) format, and the grid process was executed using autogrid.exebinary software.The molecular docking stage included selecting receptors and ligands, followed by setting genetic algorithm parameters and configuring docking parameters.The docking document was saved with Lamarckian GA 4.2 in (*.dpf) format, and the docking process was run using autodock.exe.

Analysis and Visualization
Docking results in format (*.dlg) were analyzed using AutoDock Tools, focusing on binding energy (∆Gbind) and inhibition constant (Ki).Furthermore, ligand and amino acid interactions were visualized, and the docked Ligand-Receptor Complex in (*.pdb) format was opened with Biovia Discovery Studio 2019 for 2D and 3D visualization.

Molecular Dynamics Receptor and Ligand Preparation
VEGFR-2 receptor and ligand were prepared for conduction of molecular dynamics simulation.Sorafenib and the best ligand obtained from docking analysis were used, then complexes formed between the protein and each ligand were saved in (*.pdb) format.

System Setup
Ligand-receptor complex format (*pdb) files were input into CHARMM-GUI (https://www.charmmgui.org/)using Solution Builder on Input Generator.During this process, topology generation, water box creation, ion addition, and force field selection were conducted.

System Minimization
System minimization was performed using the GROMACS application to prevent inappropriate collisions in the simulation system.Additionally, this helped loosen structural complexes by reducing potential energy to meet the requirements for simulation

Equilibration
Solvent and ions in the simulation system were equilibrated to achieve thermostat and barostat conditions at 310 o K temperature and 1 atm pressure.

Production
The ligand-protein complex was subjected to molecular dynamics simulation to obtain an atomic trajectory.

Analysis of Molecular Dynamics Results
Production results were analyzed for parameters including RMSD, RMSF (Root Mean Square Fluctuation), Radius of Gyration (Rg), and hydrogen bonds, then atomic trajectories were visualized using Grace software.

RESULT AND DISCUSSION
The Lipinski Ro5 analysis was conducted to predict the physicochemical properties of the native and test ligands.Lipinski's Ro5 is a rule for designing drugs expected to penetrate membranes when administered orally.This includes molecular weight (BM) of no more than 500 g/mol, a partition coefficient (logP) value < 5, hydrogen bond donors < 5, and hydrogen bond acceptors <10.Based on the rule, the positive control Sorafenib and 25 of the 36 test ligands met the four specified criteria.
In this context, physicochemical criteria for acceptable drug absorption and permeability constitute the first step in oral drug bioavailability 12 .
The binding or active site of the receptor was determined 13 using PyMOL by inputting the GDP file of the original ligand, followed by typing 'center of mass' on the PyMOL command line to obtain x, y, and z values (x = -24.964;y = -0.497;z = -10.399).Subsequently, the grid box was arranged as grid center x = -24.964,y = -0.497,and z = -10.399with dimensions of 40x40x40 Å and spacing of 0.375 Å.This box was used as a reference for the molecular docking process of compounds to be tested.Additionally, the redocking results were analyzed with PyMOL, which produced an RMSD value of 1.136 Å, suggesting the applied docking procedure to be valid due to the generation of an RMSD < 2 Å.A greater RMSD value often leads to a larger deviation, while a smaller RMSD value makes the ligand pose resulting from docking to be closer to the crystallographic ligand pose. 14.

Figure 1. RMSD Result
Docking results stored in (*.dlg) format were analyzed using AutoDock Tools for parameters including binding energy and inhibition constants (Ki), then ligand and amino acid interactions were visualized.The binding energy value and Ki are parameters of the conformational stability of the investigated compounds in terms of affinity towards VEGFR-2 receptor 15 .The binding energy value was obtained through the interaction of the ligand with a receptor, which showed the suitability between the shape and size of the ligand and the binding site on receptor 16 .Moreover, molecular docking results of the native ligand Sorafenib showed a binding energy of -11.51 kcal/mol and a Ki value of 3.67 nM.The test ligand Isobavachalcone (ISB) had the smallest binding energy of -11.45 kcal/mol, which was close to the value obtained for Sorafenib, along with Ki of 4.05 nM.Negative binding energy leads to a lower Ki value due to a directly proportional relationship between both parameters 17 , as presented in the results shown in Table 1.  2. These results presented the interactions of amino acid residues with ligands, showing potential contact between ligands and the protein that initiated inhibitory activity.The binding area or site is the point where ligands attach to a particular part of the protein, thereby influencing conformation and function.Additionally, the site shows amino acid residues forming interactions, such as hydrogen bonds, hydrophobic bonds, and electrostatic bonds, with the macromolecule and ligands 18 .Table 2 shows the interactions between the macromolecule and sorafenib, as well as the top 10 test ligands.
The native ligand Sorafenib as the positive control formed hydrogen bonds with VAL 899, CYS 919, ASP 1046, and GLU 885, as well as a carbonhydrogen bond at GLU 917.The 31 test ligands formed hydrogen bonds with amino acid residues partly identical to those bonded by the positive control.ISB also formed hydrogen bonds at residues CYS 919, ASP 1046, GLU 885, and GLU 917.Based on this perspective, more hydrogen bond formation increases the potential of ligand activity 19 .
Sorafenib formed hydrophobic bonds with amino acid residue LEU 1035, which was similar to the interactions shown by the test ligand ISB.Moreover, ISB formed other identical hydrophobic bonds at residues ILE 892, LEU 1019, LEU 840, LEU 1019, ILE 1044, PHE 918, LEU 889, VAL 848, LYS 868, VAL 916, CYS 1045, ALA 866, and CYS 919.These hydrophobic bonds play an essential role in drug design due to the ability to enhance binding affinity between a drug and the target 20 .The molecular dynamics simulation process included selecting the best test ligand (ISB) based on molecular docking results of E. fusca compounds.Docking results in (*.dlg) format were analyzed using AutoDock Tools, generating a protein-ligand complex saved in (*.pdb) format.This complex file was then used for molecular dynamics simulation, starting with system preparation on the CHARMM-GUI website.At this simulation stage, topology generation, water box creation, ion addition, and force field selection were performed.The topology was generated automatically from the input complex file, followed by the automatic formation of a rectangular, cubic crystal-type water box.Subsequently, the system was solvated by adding 13,737 water molecules for the Sorafenib complex and 13,702 water molecules for the ISB complex.
The addition of 0.15 M NaCl ions (Na+ = 39, Cl-= 38) was performed to simulate biological systems, which generally contain several salts.To model conditions closely resembling those in biological systems, molecular docking simulation is conducted under electroneutral conditions 21 , and total charges are neutralized by adding ions before minimization 22 .In this research, the CHARMM36 (C36) force field and an input generator were selected, then hydrogen mass repartitioning was applied and simulation was run using GROMACS at a temperature setting of 310 o K.The C36 force field allows for easy simulation of complex heterogeneous systems in various packages, including GROMACS 23 .This force field was optimized using the standard Transferable Intermolecular Potential Three-Point (TIP3P) solvent model for proteins 24 .TIP3P is part of the CHARMM simulation package frequently used to simulate biological systems 25 .
Minimization was performed for 5,000 steps, with termination conditions set to a maximum force of <10 kJ/mol.The results showed that sorafenib had a potential energy value of -7.2752488e+05 kJ/mol, with a maximum force of 9.0506854e+02 kJ/mol on 4913 atoms and a normal force of 1.7884336e+01 kJ/mol.Meanwhile, ISB minimization generated a potential energy of -7.1923694e+05 kJ/mol, with a maximum force of 9.5177234e+02 kJ/mol on 4913 atoms and a normal force of 1.5471794e+01 kJ/mol.The systems were equilibrated for 125,000 steps and 125 picoseconds (ps) using the Parrinello-Rahman barostat method.Generally, this method maintains constant pressure and is recommended for NVT/NPT production simulation 26 .
The molecular dynamics production stage was run for 5,000,000 steps and a time of 20 nanoseconds (ns).The production simulation of Sorafenib required 18 hours 46 minutes, while the duration for ISB was 18 hours 33 minutes.Simulation results were analyzed for parameters, including RMSD, RMSF, Rg, and hydrogen bonds, then atomic trajectories were visualized using xmgrace 27 .These analyses determined the stability of the interactions between the ligands and receptors over space and time 28 .
RMSD analysis determines the magnitude of binding pose changes over time, assessing the stability of ligand-protein interactions 21 .The graphical representation of RMSD values shows variations that occur during simulation 27 Initially, the backbone RMSD simulation presented an increasing RMSD value, showing the protein structure opening and the process of ligand searching for suitable binding sites.Stabilizing RMSD values suggest the presence of interactions between protein residues and the ligand, which facilitates protein structure maintenance 29 .The protein-BAX and protein-ISB complexes had RMSD values ranging from 0.175-0.225nm, which were less than 3Å (0.3 nm), thereby suggesting good system (complex) stability 30 .Meanwhile, RMSF analysis determines fluctuations in ligand interactions with amino acids during simulation.Bonded amino acid residues are presented in Table 2 showing the results of VEGFR-2 macromolecule interactions with ligands 28 .RMSF for binding residues ranged from 0.0394-0.0730nm in the positive control complex and 0.0392-0.1026nm in the ISB complex, with values less than 3Å (0.3 nm) suggesting good system (complex) stability 30 .The radius of Gyration (Rg) is defined as the distribution of protein atoms around the axis of a particular protein.When a ligand binds, a conformational change occurs, leading to alteration in Rg, which aids in predicting drug-protein compactness and binding patterns 27 .During the 20000 ps simulation, Rg values ranged between 1.95-2.0nm, with a graphical representation showing stable compactness for the two ligand-protein complexes formed in this research.Hydrogen bonds play an important role in the formation of receptor-ligand complexes 31 and are often found among relatively stable compounds in molecular dynamics simulation 32 .From the results, it was observed that lower hydrogen bonding correlated with increased RMSD deviations 27 .Meanwhile, the graph shows total hydrogen bonds formed during the 20,000 ps simulation, which presents stable interactions in both complexes.

CONCLUSION
In conclusion, molecular docking results of E. fusca chemical compounds showed that the best test ligand (ISB) had a binding energy of -11.45 kcal/mol, compared to the native ligand (Sorafenib) with a value of -11.51 kcal/mol.These chemical compounds interact with amino acid residues of the VEGFR-2 receptor through hydrogen, hydrophobic, and electrostatic bonds.Some of the residues that interacted with the investigated chemical compounds were similar to those bonded by the positive control.The amino acids included in the hydrogen bonds were CYS 919, ASP 1046, GLU 885, and GLU 917.Moreover, the stability of the ligandprotein bonds was confirmed through molecular dynamics simulation using ISB and Sorafenib.Molecular dynamics simulation produced RMSD, RMSF, Rg, and hydrogen bond parameters.ISB and Sorafenib complexes met the RMSD and RMSF value requirements of less than 3Å (0.3 nm).Additionally, the Rg value was stable until the completion of the simulation, and the movements of the ligand-protein complexes were identical.The results of the hydrogen bond interactions between the two complexes were relatively stable per unit of time.Therefore, these hydrogen bonds played an essential role in receptorligand complexes formation.

Table 1 .
Binding energy results and Ki against VEGFR-2 macromolecule

Table 2 .
The results of VEGFR-2 macromolecule interactions with ligands