Molecular Docking and Molecular Dynamics Simulation using Monascus sp. as a Candidate Cervical Cancer Drug

Cervical cancer is the fourth most common female cancer worldwide and results in over 300000 deaths globally. Given that HPV prophylactic vaccines do not exert a therapeutic effect in individuals previously infected, it is unlikely that HPV-associated cancers will be eradicated in the coming years. Therefore, there is an emerging need for the development of anti-HPV drugs. The purpose of this study is to find out Monascus sp. as cervical anticancer using molecular docking and dynamics methods. The results of docked with AutodockTools were visualized with Pymol, analyzed the effectiveness using the Ramachandran plot. The docking results show that there are 2 pigments that have lower ∆ G than raloxifen in estrogen receptor beta with the lowest ∆ G indicated by the pigment Monascin and Ankaflavin , which is -6.94 kcal/mol with Ki value of 39.49 nM and -6.22 kcal/mol with Ki value of 27.78 nM. The results of molecular dynamics, Ankaflavin and Monascin have good stability at estrogen receptor beta because the outlier area has a value 11.722% and 10.256%. And the amino acid residues in the most preferred area were 68.864% and 70.330%. In addition, Monascopyridine B and Monascuspiloin pigments showed good and stable results at the EGFR receptor because the outlier areas were 14.692% and 10.623%. And the amino acid residues in the most-favored region were 65.403% and 73.260%. Based on the results of this study, we predict that Ankaflavin , Monascin , Monascopyridine B and Monascuspiloin can be used as new cervical anticancer candidates after validated with in vitro and in vivo tests.


Introduction
Cervical cancer is considered as one of the deadliest cancers besides breast cancer [1].The prejudgment that more than 570,000 new cases of cervical cancer in the world.Around of 90% death case cause it happen in developing and middle-income countries [2].The number of cervical cancer sufferers in Indonesia ranks second in Asia after Mongolia, which is 23.4 per 100,000 women on one year.Every year in Indonesia 8000 cases of death occurred from more than 15000 cases of cervical cancer.About 20-25 women die from 40-45 cases of cervical cancer are found in every day [3].
Estrogen through its core receptors, namely estrogen receptor alpha (ERα) and estrogen receptor beta (ERβ), and the GPR membrane receptor, affects physiological processes in various tissues including the female reproductive tract.Therefore, it is not surprising that estrogen is involved in various human diseases including cervical cancer.Estrogen alpha and beta estrogen are present in the cervix, required for dynamic changes in the cervical epithelium.The human estrogen alpha receptor (ERα) is a core receptor that is important in regulating female reproductive development and function.And EGFR has a role in cells as an enhancer of cell proliferation, angiogenesis, and inhibits apoptosis so that it is needed in normal cell cycle conditions.If there is a malfunction of the EGFR, it will cause several conditions such as mutation, amplification, decreased phosphatase, heterodimerization, and overexpression of EGFR (overexpression) so that it triggers the occurrence of tumors or cancer [4].
The importance of this research is based on the benefits of the Monascus sp pigment, which in addition to functioning as a natural dye, also has biological activity, namely as an anticancer.From the pigment groups, both main compounds and derivatives produced during the fermentation process, information about pigment compounds that act as anticancer is still very limited, so it is necessary to carry out developments to determine the selected pigment compounds for further synthesis to produce new compounds that are easier to manufacture.Without going through a long fermentation process and requiring special attention.The new compound produced is expected to be an alternative as one of the candidates for anticancer drugs, especially cervical cancer which is easy to obtain and at affordable prices.
In previous studies, based on in silico predictions, it was stated that monascus pigments have potential anticancer activity at estrogen beta receptors [5].In this research, was means to found a better anticancer activity from Monascus derivatives using molecular docking and molecular dynamics.The docking method is a molecular modeling method to predict the interactions that occur between macromolecules (receptors) and small molecules (ligands) that are bound to form complexes efficiently.In this case, the most stable conformation of the receptor and ligand was sought by analyzing the docking complex which gave the lowest binding affinity.Docking analysis is widely used to predict interactions and the resulting energy [6].Molecular dynamics can study the effect of the presence of an explicit solvent in the system, then it is applied to explore the conformation of protein receptors to improve the drug design process [7].
The main obstacle to cancer is the problem of treatment.Until now drugs to cure cancer is still very limited.Apart from the lack of types and quantities, when available, the price will be very high, so this is one of the factors that makes it very difficult for cancer sufferers to get treatment.So, the purpose of this research is to obtain cervical anticancer drug candidate compounds through molecular docking and dynamics methods.

Preparation of the 3D Receptor Cancer
The Cervical cancer reseptor used in this study were Estrogen alfa, Estrogen beta dan EGFR were downloaded from Protein Data Bank (http://www.rcsb.org/pdb/)with PDB ID 3ERT, 4ZI1 and 3W2S.The receptors preparation was accomplished by removing water molecules and natural ligands, adding polar hydrogen atoms, and calculating the Kollman charge using MGLTools 1.5.6 with AutoDock 4.2.6 [9].

Optimization of the 3-Dimensional Cervical Cancer Receptors and Ligand with the Density Functional Theory (DFT) B3LYP method
The molecular structures of the 3 cervical cancer receptors and five Monascus sp.pigment derivatives compounds from PubChem (https://pubchem.ncbi.nlm.nih.gov/),followed by reoptimization using the best optimization method in Gaussian 09 software.The computational method used is DFT with the B3LYP function and the base set used is 6-311G (Becke, 1993).B3LYP is a functional hybrid proposed by Becke [10].

Ligand Based Drug Likeless Screening (Drug scan)
The physicochemical parameter testing was used the SwissADME software on internet website (http://www.swissadme.ch/)based on the structure of the compound in SMILES format.The physicochemical parameters are based on Lipinski's rule (Lipinski Rule of Five).The parameters used include molecular weight <500 g/mol, lipophilicity <5, hydrogen bond donor <5, hydrogen bond acceptor <10, and refractory molar between 40-130 [11].

Molecular docking Method Validation
Validation method of molecular docking parameters was done using re-docking method.In this process, the value of the Root Mean Square Deviation (RMSD) have to less than 2 Å.The smaller the RMSD value, the better the predicted ligand position because the closer it is to the original conformation [12].

Ligand-Receptor Cervical Cancer Docking Simulation
The optimized Monascus sp pigment compound was docked to the alpha, beta estrogen receptors and EGFR receptors whose native ligands had been removed using the AutodockTools application with the similar docking procedure as validation method.The results of the analysis showed the conformation of compound bonds in proteins with the value of bond energy and hydrogen bonds [8].Protein was rigidity was kept while the ligand was flexible during the docking [9].Based upon the binding affinity the ligands were screened.The docked poses were further analyzed using Pymol for interaction studies.

Molecular docking Data Analysis
The Analysis of molecular docking data result was the bond energy and hydrogen bonds formed.The bond energy was used to indicate the bonds strength between compounds and proteins.The lower the bond energy value, the stronger and more stable the bonds were.The type of hydrogen bond formed is used to analyze the interaction mechanism formed [13].

Molecular Dynamic File Preparation and Simulation
The docked structure from docking analysis is further used for MD simulations for the better understanding of the drug binding.The molecular dynamics simulation was used MOE 2009.10 platform.Before the MD simulation takes place, the energy of the complex system is minimized.The AMBER99 force field and the Generalized Born implicit solvation model (GBIS) were used for this.Initialization time is determined by dynamics of 100 picoseconds (ps) at 300°K.Complex stability was determined by simulation of 5000 ps after being heated to 300°K (cervical cancer patient temperature) for 20,000 ps.The simulation ended with a cooling/annealing step up to 1°K for 10 ps to obtain a complex conformational structure that has the lowest energy.Cooling stage/Annealing aims to find the lowest conformational energy of the molecule.The Ramachandran plot of the most optimal ligand will be described to determine the structural integrity of the protein.Position velocity and acceleration results are saved every 0.5 ps.The output of the simulation results is in the form of a database in mdb format [14].

3
Results and Discussion

Ligand Based Drug Likeless Screening (Drug scan)
One of the vital factors to discover and progress the bioactive ligands as an oral drug is their high oral bioavailability.The significant forecasters of good oral bioavailability are rotatable bonds of molecules which became known under molecular flexibility, a good gastrointestinal absorption, and a low polar surface area (total of acceptor and donor hydrogen bonds).TPSA of a molecule is defined as the surface sum over all polar atoms or molecules, primarily oxygen and nitrogen, also including their attached hydrogen atoms.Additionally, Lipinski et al. reported druglikeness features such as MW, hydrogen bond donor and acceptor values, LogP (partition coefficient) under the 'Rule of Five' term [15].This rule helps to make an easier selection of molecules with better pharmacokinetic properties in the human body for oral formulations.In this step, physicochemical parameters (e.g.CLogP, TPSA, MW, hydrogen bond donor and acceptor, TPSA, and rotatable bonds) of the approved compounds from the previous step were evaluated.Finally, 5 of Monascus sp.pigments were fitted with these properties and revealed a satisfied druglikeness demonstrating perhaps a good permeability through the biological membranes.

Cervical Cancer Ligand-Receptor Docking Simulation
Molecular docking was used AutodockTools 4.2.6 software with flexible ligand position to might the ligand have stable structural adjustments when binding to the receptor, with the same grid box arrangement, with the aim of directing the ligand of the test compound to interact within the receptor region [16].
The purpose of the docking simulation of the test compound is to determine the compound that has the best binding affinity for the target protein receptor.Experimentally binding affinity (∆G) is directly related to Ki according to the equation: ∆G = -RT ln Ki, so that the value of ∆G (Table 1) is able to predict the ability of a compound to inhibit protein [17].
The effectiveness of the molecular docking of the ligand and receptor is realized by the interaction complex of the two molecules which represents the ability of the ligand to bind to the target receptor using the following parameters: 1) Interaction surface; 2) The value of the binding affinity; 3) The value of the inhibition constant (Ki).Surface interactions show the large area of binding of the ligand to the receptor, which is directly proportional to the strength of the binding of the ligand and receptor complex.The amount of free bond energy is the amount of energy needed for the molecular complex to bind, especially on the stability of the binding of the receptor ligand complex [18].The value of the inhibition constant can determine the effectiveness of ligand inhibition on receptor activity.The smaller the value of Ki and the binding affinity, the better the interaction between molecules has the better effectiveness.Energy ∆G ≤ 0 illustrates that the test compound has affinity for the active site of the receptor and is able to inhibit the target protein well and stably.The docking simulation results showed that 2 of the 5 Monascus pigment derivatives had better affinity than Raloxifen for the beta estrogen receptor (4ZI1).The results of the docking show that the 2 pigments have a binding affinity (∆G) lower than raloxifen with the free energy of Monascin pigment -6.94 kcal/mol and Ankaflavin pigment -6.22 kcal/mol.Energy ∆G pigments Ankaflavin and Monascin are more potent when compared to their comparison, Raloxifen, which is a drug provide effective results to prevent and treat cervical cancer with a treatment period in 1 month, and it is a member of the Selective Estrogen Receptor Modulator (SERM) class of drugs.
This phenomenon shows a good signify that Monascin and Ankaflavin have strong interactions and bonds in the active site area of the target macromolecule.After that, all ligandreceptor docking complexes were selected for further study using a simulation of the dynamics of the ligand-receptor interaction.

Visualitation Docking Result
Visualization of the docking results was carried out to determine the interaction between the functional group of the ligand and the amino acid residue of the test compound contained in the receptor [19].
The interactions of 5 Monascus pigment derivatives and 3 cervical cancer receptors can be seen in Table 2.There shows the interaction of hydrogen and non-hydrogen bound (hydrophobic) residues of 5 Monascus pigments compared to the standard drug.Some kind of the hydrophobic interactions referred to here there are van der Waals bonds, Pi-Alkyl, Pi-Anions and dipole bonds.Hydrogen bonds are interactions that are formed between hydrogen atoms and electronegative atoms such as flour (F), nitrogen (N), and oxygen (O) [20].Hydrophobic interactions are interactions that avoid the liquid environment and tend to group inside the globular structure of the protein so as to minimize the interaction of nonpolar residues with water [21].
Hydrogen bonding and hydrophobic interactions in the binding between ligandreceptor, because they play an important role in ligand stability to estrogen receptors alpha (3ERT), estrogen beta (4ZI1) and EGFR (3W2S).Hydrogen bonds can stabilize the sequence of amino acid residues in the three-dimensional structure of proteins and determine the dynamics and stability of proteins [22].Besides acting as a determinant of ligand stability to the receptor, hydrophobic bonds play a role in minimizing the interaction of nonpolar residues with water [23].The results of visualization of hydrogen bonds and hydrophobic interactions formed on the test ligands show results that are close to the results of interactions from the comparison, namely Raloxifen.This shows that the compound shows its activity as an inhibitor.This is supported by the test ligands with amino acid residues and hydrogen bonds that are close to Raloxifen which shows similar types of interactions in this case describing similar activities [24].In addition, the more interactions between the test compound and amino acid residues, it is predicted that the interaction will be better [25].

Analysis of Molecular Dynamics Simulations
The suitability of amino acids in the comparison of ligand interactions with receptors on molecular docking and molecular dynamics affects in the future the mechanism or activity of a compound on the process of inhibition of a cell, especially cancer which is focused on in this study.The greater the number of the same amino acids from the Monascus sp.It can be assumed that these compounds have good inhibitory abilities.During the simulation there is a change in the interaction between the ligand and the receptor, this can cause the amino acid residues produced in the molecular docking and molecular dynamics to be different, and will experience a reduction in the number of interactions.The increasing match of amino acids before and after performing molecular dynamics parameters indicates that the compound is stable and resistant to thermodynamic changes [26].Based on the analysis and visualization of the data from the in silico study using molecular dynamics on Table 3, it is known that the pigments Ankaflavin and Monascin to the estrogen beta receptor (4ZI1) and the pigments MonascopyridineB and Monascuspiloin to the EGFR receptor (3W2S) have the number of amino acid matches (x) resulting from molecular dynamics and molecular docking with the most means has the lowest conformational value and indicates that the compound is more stable than other compounds.From the results of molecular dynamics, it shows that there is a change in the type and number of amino acid residues in the anchorage area and there are similarities in hydrogen bonds, so that during the simulation the ligand does not come out of the protein which indicates that a stable complex has formed between the ligand and the receptor.This can occur due to conformational changes in the ligands and proteins due to the influence of temperature and the presence of solvents in the simulation system so that the previous atomic distance changes so that finally the ligand groups interact with proteins that form complexes between the drug-receptors.The conformational stability of the protein at a temperature of 312°K which was run for 2000 ps in this simulation indicated that Ankaflavin, Monascin, MonascopyridineB and Monascuspiloin compounds could be predicted to have activity as cervical anticancer candidates.
Determination of protein stability can be analyzed for non-glycine residues in the disallowed region of the resulting plot.It can be seen in the Figure 1, the amino acid residues in the region of the dihedral angle (phi) and the dihedral angle (psi) of the peptide bond.The protein structure is said to have good stability if the non-glycine residue in the disallow region (outlier) is less than 15% and the amino acid residue in the most favored region is more than 50% [27].The dynamic result data showed that all the receptors for the test ligands had good quality and stable protein structure and test compounds (Table 4).The simulation results docking showed that there were 2 pigments that had lower ∆G than raloxifen with estrogen beta receptors, that pigment there are Monascin and Ankaflavin pigments, which is -6.94 kcal/mol with Ki value of 39.49 nM and -6.22 kcal/mol with Ki value of 27.78 nM.The results of molecular dynamics, Ankaflavin and Monascin have good stability at the Estrogen Beta (4ZI1) receptor, because the outlier area has a value of less than 15%, which is 11.722% and 10.256%.And the amino acid residues in the most favored region were greater than 50%, which is 68.864% and 70.330%.
In addition, the pigments MonascopyridineB and Monascuspiloin showed good and stable results at the EGFR (3W2S) receptor because the outlier areas were 14.692% and 10.623%.And the amino acid residues in the most favorite region were 65.403% and 73.260%.Based on the results of this study, we predict that Ankaflavin, Monascin, Monascopyridine B and Monascuspiloin can be used as new cervical anticancer candidates after validated with in vitro and in vivo tests on the four pigments.

Table 1 .
Results of docking test ligands with AutodockTools

Table 2 .
Data Visualization Monascus sp and Raloxifen Pigments on Estrogen Alpha, Estrogen Beta and EGFR Receptors

Table 3 .
Differences in Amino Acid Residues from Docking and Molecular Dynamic Simulations

Table 4 .
Ramachandran Plot Statistics after MD receptors with test ligands a : Most favoured region, b : Generously allowed region c : Disallowed region 4 Conclusions