MULTIDISCIPLINARY TECHNOVATION Molecular docking and molecular dynamic simulation of phytochemical derived compounds as potential anti cancer agent against tyrosine kinase

There is a continuous requirement to develop novel, safe, effective and affordable anti-cancer drugs because Cancer is a serious disease at current situation. A huge number of patients die annually due to cancer disease. Phytochemical are the secondary metabolites of medicinal plants and significantly used in conventional cancer research. Bioactive phytochemical is favored as they claim differentially on cancer cell only without altering normal cell. Carcinogenesis is an intricate process and includes multifold signaling procedures. Phytochemical are pleiotropic in nature, function and target these events in multiple manners so they are considered as most appropriate candidate for drug development. The aim of the present research was to find out the anti-cancer activity of the phytochemical constituents through computer aided drug design approach. In this experiment, we have find total 42 natural compounds with anti-cancer activity against the cancer target 1QCF tyrosine kinase. The data set comprising of phytochemical compounds was used for virtual screening and molecular docking in PyRx software. Along with screened compound, hit compound Carnosic acid was further docked to confirm the binding mode and confirmed the effective inhibition of 1QCF and anticancer activity. Molecular dynamic simulation studies were done to confirm the stability of the protein and ligand complex during a simulation. Parameters like RMSD, RMSF, and radius of gyration were experiential to understand the fluctuations. Protein-ligand interaction studies also expose that enough hydrogen and hydrophobic bonds are present to validate our results. Our study suggests that the potential use of Carnosic acid can come out as a potential candidate and in turn prevent cancer.


Introduction
Cancer, also known as a malignant tumor, is a deadliest disease concerning abnormal cell growth with the potential attack or spread to other parts of the body. The characteristics of cancer consists of six biological potentialities to support the development of tumors, which include evading growth suppressors, sustaining proliferative signaling, resisting cell death, enabling replicative immortality, inducing angiogenesis and metastasis [1,2]. Cancer is the major public health burden in both developed and developing countries where the number of cancer patients is in continuous rise. There are more than 100 different known cancers that affect humans and each is classified by the cell type that it is initially affected [3]. According to WHO, in 2012 about 14.1 million new cases of cancer occurred globally and caused about 8.0 million deaths [4]. By 2030, it is expected that there will be 26 million new cancer cases and 17 million cancer deaths per year [5]. At present, although significant efforts, cancer still remains an aggressive killer worldwide. The most general and highly effective methods of cancer treatment are chemotherapy, surgery and radiotherapy [6]. These therapies have numerous limitations and drawbacks. Radiotherapy and chemotherapy have serious side effects as well as complications such as fatigue, pain, diarrhea, nausea, vomiting and hair loss [7]. Millions of species of animals, plants, microorganisms and marine organisms act as eyecatching sources for new therapeutic candidate compounds. Medicinal plants comprise a general substitute for cancer treatment in many countries around the globe. There are more than 2, 70, 000 higher plants accessible on this planet. But only a tiny portion has been investigated phytochemicals. So, it is expected that plants can provide potential bioactive compounds for the development of new 'leads' to fight Int. Res. J. Multidiscip. Technovation, 3(5) (2021) 11-25 | 12 cancer diseases [8]. The key research programs are being undertaken in number of laboratories in various parts of the world for screening plant extracts for anticancer activity. There are several pathways that are involved in the cancer development. A plan to develop anti-cancer drug usually involves blocking one or more components of specific pathway that is involved in the cancer development [9][10][11]. The expenditure of research and development in the pharmaceutical industry has been increasing suddenly and progressively to accomplish any practice in the last decade but the amount of time required in fetching new products to market take around ten to fifteen years. Today we have high computation power with all the necessary information available online and most of them are free. We can fasten the discovery process by using molecular docking and molecular dynamic simulation studies to narrow down our selection for the anticancer compounds at the very initial stage of our process. Even we don't need to have those compounds physically to test their effectiveness and also, we don't need to handle such an experiments. It is being always observed that natural compounds so far have been found to have many of the therapeutic anticancer activity. Computer-Aided-Drug-Discovery approaches have been broadly used in drug investigation to boost the efficiency of the drug discovery and development pipeline, subject to the purpose and systems of interest [12]. Docking is a software practice that fits a molecule into target binding sites. The main objective of molecular docking is to contribute a prediction of the ligand-receptor complex structure by means of computational approaches [13]. Docking can be used to execute virtual screening on enormous libraries of compounds, rank the results based on the scores and put forward structural hypothetical theories on how the ligands inhibit the target receptor, which is vital in lead optimization [14][15][16]. Factors affecting docking are: intramolecular (bond width, bond angle and dihedral angle) and intermolecular (electrostatic, dipolar, Hbonding, hydrophobicity and van der Waals) forces. In this research, the natural phytochemicals are identified using molecular docking for treating cancer. We have searched the literature and found 42 natural compounds, which can come out as a potential therapeutic candidate for cancer. We have used binding affinity and simulation parameters to scrutinize our selected 42 natural ligands.

Role of receptor tyrosine kinases in cancer pathway
Int. Res. J. Multidiscip. Technovation, 3(5) (2021) 11-25 | 13 The RTK family consists numerous subfamilies which contain, among others, epidermal growth factor receptors (EGFRs), fibroblast growth factor receptors (FGFRs), insulin and insulin-like growth factor receptors (IR and IGFR), platelet-derived growth factor receptors (PDGFRs), vascular endothelial growth factor receptors (VEGFRs), hepatocyte growth factor receptors (HGFRs), and proto-oncogene c-KIT [18,19]. RTKs monomers are structured into an extracellular (Nterminal), a trans membrane and a cytoplasmic kinase field [20]. They are activated via ligand-induced dimerization that results in receptor autophosphorylation and tyrosine activation of RTKs' substrates including phospholipase C-γ, mitogenactivated protein kinases and phosphatidylinositol 3kinase. Mutations that affect RTK signaling often lead to cell transformation, which is observed in a wide variety of malignancies. These mutations affect RTKs or components of downstream pathways such as MAP kinase and the PI3K/AKT. This results in increased cell proliferation, survival, invasion and metastasis. Therefore, targeting RTK signaling pathways remains a challenge for scientists and clinicians working in the cancer field. Several small molecule inhibitors and antibodies are being clinically developed to target RTKs, the MAP kinase and PI3K/AKT pathways. This review attempts to highlight the important role played by RTK signaling in carcinogenesis and the therapeutic strategies available, so far, to target these important cellular pathways.

Protein structure preparation
The X-ray diffraction-based crystal structure of tyrosine kinase (1 QCF) with a resolution of 2.00 Å was taken from the protein data bank. The structure was cleaned to ensure maximum quality and reliability [21]. The bound ligands, water molecules were removed and missing atoms and residues were added. Stearic clashes were minimized and hydrogen atoms were added. Formal bond orders were determined, side chains were optimized and fixed, charges added using program implemented in chimera, SWISS PDB viewer, and Chiron minimization and refinement tool [22][23][24].

Ligand archive research
A deep literature survey was performed to discover the natural compounds having said anticancer properties. A total of 42 compounds were compounds were identified (table 1) and their 3D structures were generated by Galaxy 3D generator tool of online server Mo inspiration Cheminformatics (https://www.molinspiration.com/) (Figure 2 & 3) and 3D structure of target protein tyrosine kinase (PDB ID: 1QCF). The 3D structures were visualized by Discovery Studio Visualizer. The compounds were imported in to the PyRx (V 8.0) and energy minimization was prepared using Open Babel (Version 2.3.1) [23] module of the same software. The binding energies between receptor and the ligands are attained in terms of Kcal/mol. Energy minimization was done via the Universal force field (UFF) using the conjugate gradient algorithm. A total number of 200 steps were set and the number of steps to update was set to 1. The minimization was set to stop at an energy difference of less than 0.1 Kcal/mol [25].

Drug screening using PyRx (V 8.0)
All the 42 phytochemicals are subjected to Drug Screening process against 1QCF receptor using PyRx (V8.0) a standalone virtual screening software

Molecular Docking studies
Molecular docking is a method used to analyze the position and the inhibition interaction between the protein and the small molecules. Molecular docking was executed with PyRx (V 8.0), which is an extension of the python molecular viewer. A Lamarckian genetic algorithm was used to perform the automated molecular docking of the protein with each ligand. The torsion bonds and side chains were kept to rotate freely, while the protein structure was kept rigid. As a reprocessing step, the PDB format of macromolecule and SDF format of small molecules are converted to PDBQT format. Gustier charges were computed, and all the charges of non-polar hydrogens were assigned [26].  The docking of small molecules to the macromolecule was focused on the specific binding site. The total number of rotatable bonds of the ligand is calculated. The grid was defined to the binding site of the protein structure with the configurations of x/y/z coordinates was set to size x= 59.60, size y=76.61 and z=69.28 centre of the grid box was set to centre x= 12.53, centre y=29.84, centre z=32.10 in X, Y, Z dimensions, in which the grid was covered to the binding site of macromolecule and the grid was spaced at 0.375 Å.

ADME and Toxicity predictions
The selected 42 phytochemical compounds were further studied for their Adsorption, distribution, excretion, metabolism, and toxicity profile using SWISS ADME [27] and data warrior tools [28,29]. The predicted properties considered were blood-brain barrier penetration properties, Human intestinal absorption, inhibition to cytochrome P450 enzyme, and bioavailability. Compounds showing satisfactory properties were further studied for their toxicity profile using data warrior tools. Toxicity profiles included were mutagenicity, tumorigenicity, irritability, reproducibility, Ames toxicity, and carcinogens.

Molecular dynamic simulation
Docked protein and ligand complexes were subjected to molecular dynamics simulation using NAMD software [30]. The success of MD simulation depends on the selection of the initial protein and ligand structures. Initially, the structure was checked for inconsistencies. Out of 42 selected compounds from the docking results, we have selected the three final compounds having a PubChem number 65126, 5350 and 969516. The docked complexes were studied for their stability during the simulation. The root means square deviation, root mean square fluctuation, and radius of gyration was studied for protein backbone residue and ligand within the binding site of the simulated system [31][32][33]. The stabilities of the complexes were examined by monitoring their root mean square deviation (RMSD) during 50, 00,000 steps for a 10 ns simulation. MD simulations were performed using the CHARMM36 force field [34]. Visual molecular dynamics (VMD) was used to generate PSF files for complex. All complex was solvated in cubic water boxes containing transferable intermolecular potential with 3 points (TIP3P) water molecules. The box size was chosen to match the molecular dimensions so that there was a distance of 5˚A between the protein surface and the edges of the periodic box. A 5˚A cut off distance was used to calculate short-range non bonded interactions. The particle mesh Ewald (PME) method was used to calculate long-range electrostatic interactions. The SHAKE method was used to constrain all bonds involving hydrogen atoms. A conjugated gradient system was used for energy minimization, with all parameters set to default. The system first performed 10000 steps of conjugated gradient with energy minimization. We used Langevin dynamics with pressure control so our system was not an NVT ensemble. The Nose-Hoover method was used to Int. Res. J. Multidiscip. Technovation, 3(5) (2021) 11-25 | 18 maintain a constant temperature. The time step of each simulation was set to 2 fs [30,35,36]. Visualizations and data analysis were performed with VMD software [37].

Virtual screening and docking results
Virtual screening helps us to screen the biological molecules with good binding affinity. In this study, we have used PyRx 8.0 tool to screen out the molecules. A total of 42 natural ligands were selected and were docked to the target protein. The docked compounds were examined in the Auto dock tool and binding free energy was calculated [38,39]. We have selected a total of 42 compounds on their binding affinity ranging from -10.4 to -4.0 kcal/mol ( Table 2).

ADMET analysis
We have selected 42 compounds and the same compounds were studied for their ADMET properties. The properties like Human intestinal absorption, irritability, reproductive effect, inhibition to cytochrome p450 enzyme, and several others were predicted. It was clear from the results that compound number 3220, 3314, 5350,7794 ,65030 ,65036  ,65126,73201 ,442793, 558221 ,637563, 637566,  969516, 1549026, 1550607, 3083634, 3707243,  5281650 ,5281654,5281794 ,5318556, 5321010,  6442492, 6917864, 9793905, 13033058, 17174331, 3220 have a high intestinal absorption value. From the Table 2, it was also observed that we had several compounds like 73201, 442793, 3707243, 5281654, 5281794, 6442492, 135352544 have high intestinal absorption values but at the same time they are also showing their inhibitory properties towards the Cytochrome P450 enzymes, so such compounds were removed from the further selections. Similarly, from the selected compounds with high intestinal absorption values and negative inhibitory actions to cytochrome P450 enzymes, compounds were also studied for their mutagenic, tumorigenic, irritability, and reproductive effect. Table no 2 indicates that compounds number 5350, 969516 have either of the said effects, so these compounds were also removed from the study. Finally, we selected one ligand namely Carnosic acid CID 65126 for further analysis. The ligand selected for further study was having either hydrogen bonds or presents a hydrophobic interaction with the protein Table 3 and 4.

Protein-ligand interaction
The hydrogen bond and hydrophobic interactions of protein-ligand complexes were analyzed by LigPlot+ (v 1.4.5)38 and Protein-ligand interaction profiler. "LigPlot+" is a graphical system that generates multiple two-dimensional (2D) diagrams of ligandprotein interactions from docked complexes. PLIP is complementary to another state of the art tools like a SWISS dock, galaxy site, or ProBis and thus it can be used to study the protein-ligand complex. The server allows comprehensive detection and visualization of protein and ligand complexes along with interaction pattern (Figure 4a-b).

Molecular Dynamic Simulation studies
We assessed the residue RMSD to study the residue behaviour of the protein during the simulations. In general, a residue's RMSD value was considered to represent the local flexibility of a protein and ligand complex. It reflected the mobility of an atom during the MD simulation trajectory. Therefore, a higher residue RMSD value indicated higher mobility; conversely, a lower residue RMSD value indicates lower mobility. To investigate the fluctuations in the ligand-binding energy as well as the motions of the amino acid residues within the complex during the simulation, the root means square fluctuation (RMSF) of the complex was also monitored. Besides, the compactness of each complex was determined by carefully examining how folded or unfolded the protein-ligand complex was by calculating the radius of gyration [33]. Based on the docking analysis 42 compounds were selected for further ADMET investigation and it leads us to select the final compound Carnosic acid (CID 65126) to consider the structural stability of each protein-ligand complex by molecular dynamic simulation. The stability of complex (1QCF-Carnosic acid) was monitored using root mean square deviation (RMSD) during 10 ns simulation studies ( Figure 5 & 6).
The values presented in (Table 5) for proteinligand complex studied for its stabilities during 10 ns simulation. From the values, it is clear that the range of RMSD obtained for the complex follows the acceptance range between 1 to 3.5 (Å). It is also observed from the graphs that the complex was also equilibrated as the average RMSD values are stabilized at the end of the 10 ns simulation. This fixed range of RMSD was indicating the interaction between bound ligand and flexible loop region, as it reduces the flexibility of the protein-ligand complex (Figure 7).

Protein-Ligand complex Mean RMSD (Å) Min RMSD (Å) Max RMSD (Å)
1QCF-Carnosic Acid 1.9627 0.4300 3.7910    The root means square fluctuations (RMSF) were assessed and plotted to equate the flexibility of each residue in the-ligand-protein complexes. The RMSF of the protein-ligand complex denoted the minimized fluctuation for all the complexes. The RMSF did not deviate much during the simulation period of 10 ns and the average RMSF values were kept constant for all the complexes. The radius of gyration was also monitored during the 10-ns MD simulation for each protein-ligand complex to ascertain whether the complex was stably folded or unfolded. If the radius of gyration remained relatively constant, the complex was considered to be stably folded, otherwise, it was considered to be unfolded.
In this study, the radius of gyration value obtained is listed in Table 6. The values obtained for the Carnosic acid showed a relatively constant radius gyration during the simulation. So, we can conclude that the complex is formed relatively stable folded polypeptide structures during the 10 ns MD simulation

Conclusion
The globe is moving on the way to use natural products due to their low cost and trustworthiness over side effects resulted from existing drugs. Researchers are raising their efforts for the development of phytopharmaceuticals against severe metabolic syndromes including cancer. Bioactive phytochemicals/formulations are potential leads for the development of safer anticancer drugs. Several plants and their constitutive phytochemicals have been screened for this purpose but only a very few have reached up to the clinical level. Anticancer phytochemicals described in this article must be further researched in clinical trials for their effectiveness and toxicological documentation. They must be developed as drug gable forms with sufficient bioavailability. Moreover, we know that a traditional herbal preparation has greater medicinal effect than the same phytochemical/molecule taken in a pure form. So therapeutic intervention based upon the combination of anticancer molecules may give potent and effective therapeutic results. The treatment of chronic disease like cancer with phytochemicals is critical in the research studies. Us in silicon approach on phytochemicals against cancer target 1QCF is carried out using virtual screening, molecular docking and ADMET methods. Virtual screening of three compounds showed the binding affinity towards target 1QCF. The compounds were screened with least binding affinity and compound Carnosic acid was selected as hits. The molecular docking of the hit showed the binding mode and interaction energy. H-bond pattern was analyzed and confirmed the inhibition of cancer target 1QCF to show the anti-cancer activity of phytochemicals. Our result based on in silico studies, concluded that Carnosic acid inhibit 1QCF and possessed better anticancer activity against 1QCF. Further studies on lead Carnosic acid can be done in experimental studies to confirm the inhibition and may be used in the treatment of cancer.