Jurnal Computational Studies of Thiourea Derivatives as Anticancer Candidates through Inhibition of Sirtuin-1 (SIRT1)

three thiourea derivatives to the receptor. In addition, pharmacokinetic parameters, toxicity, and selection of Lipinski's Rule of Five were also tested. Molecular docking results showed that compound b ([2-(methylcarbamothioylcarbamoyl)phenyl]benzoa te) had the lowest ∆G value of −9.29 kcal/mol with a KI value of 0.156 µM compared to other thiourea derivatives and was proven by molecular dynamics tests for 30 ns and amino acids that play an active role in the interaction include the residue PheA:297. In terms of pharmacokinetics and toxicity, compound b is better than natural ligands. Compound b is predicted to be used as an anticancer candidate through further research.


Introduction
Cancer is a general term for a large group of diseases that can affect any part of the body. Other terms used are malignant tumor and neoplasm. One of the hallmarks of cancer is the rapid creation of abnormal cells that grow beyond their usual limits. They can then invade adjacent parts of the body and spread to other organs; the latter process is called metastasis. Cancer is the leading cause of death worldwide, accounting for nearly 10 million deaths in 2020 [1]. Cancer arises from transforming normal cells into tumor cells in a multi-stage process that generally progresses from pre-cancerous to malignant tumors. These changes are the result of the interaction between a person's genetic factors and three categories of external agents, including physical carcinogens, such as ultraviolet and ionizing radiation; chemical carcinogens, such as asbestos, components of tobacco smoke, alcohol, aflatoxins (food contaminants), and arsenic (drinking water contaminants); and biological carcinogens, such as infections from certain viruses, bacteria, or parasites [2,3].
The SIRT1 protein is one of the Silent Information Regulator 2 (Sir2) protein members. SIRT1 has received considerable attention as an epigenetic regulator in human physiology [4]. Alteration of sirtuin expression is critical in several diseases, including metabolic syndrome, cardiovascular disease, cancer, and neurodegeneration. SIRT1 protein is an NAD+dependent histone deacetylase, where the cell culture system of SIRT 1 can inhibit apoptosis. Inhibition of SIRT1 or broad-spectrum sirtuin inhibitors may provide therapeutic benefits in certain types of cancer. Many SIRT inhibitory compounds, such as cambinol, sirtinol, and salermide, effectively induce cancer cell apoptosis
Thiourea is a compound containing sulfur and nitrogen, widely used in studies as anticancer, such as nitrosourea, hydroxyurea, and 5-fluorouracil [7]. In another study, it was stated that phenylthiourea derivatives could inhibit melanin production showing strong anticancer activity against melanoma [8], and the compound N-benzoyl-N'-(4-fluorophenyl) thiourea is predicted to inhibit SIRT1 as an anticancer [9,10,11]. In addition, other thiourea derivatives have activity as a central nervous system depressant [12], antimicrobial [13], antitumor, and anti-inflammatory [14]. Based on density functional theory (DFT) studies, molecular docking, and molecular dynamics, it is known that thiourea derivatives provide antineoplastic effectiveness through their inhibition of protein kinases [15].
Several steps were conducted in drug discovery and development, including prediction of the physicochemical properties of compounds by observing their pharmacokinetic properties (ADME), toxicity, and investigating the interaction of compounds (ligands) with receptors (proteins) in silico. The molecular docking method is used to predict the binding orientation of the compound to its target protein in determining the affinity and activity of small molecules [18]. In the molecular docking study, the compound 1-(4-methoxyphenyl)-3-(pyridine-3-ylmethyl)thiourea has a binding affinity value of −7.41 kcal/mol [19]. In another study, the compound 3-(naphthalen-2-yl)-N-(([3-(trifluoromethyl)phenyl]carbamoothioyl)amino)-1,2,4oxadiazole-5-carboxamide had the best potential as an inhibitor SIRT-1 with an IC50 value of 13 IC50 [20].
In previous research, researchers have researched thiourea derivatives as anticancer either in silico or in vitro [21,22,23]. Based on this background, the researchers modified three other thiourea-derived compounds ( Figure 1) and studied their interactions with the SIRT1 protein. The in silico studies investigated pharmacokinetics, molecular docking, and molecular dynamics. Molecular docking is performed using AutodockTools because of the free, open-source software suite used in the virtual screening of small molecules to macromolecular receptors. AutoDock has been applied to various biological systems for proteinligand binding and prediction of binding affinity with good correlation with experimental binding affinity for several protein systems [24]. However, visualization of two-dimensional interactions must utilize other program facilities, for example, discovery studio visualizer.

Methods
This research was conducted in silico by performing physicochemical screening through pharmacokinetics prediction, toxicity, Lipinski's Rule of Five analysis, molecular docking, and molecular dynamics of three thiourea compounds against the SIRT1 receptor.

Equipment and Materials
The tools used in this research include computer hardware with specifications Intel® Core™ i5-6400 CPU @ 3.90 GHz (4CPUs), Nvidia Geforce GTX 970 Gigabyte OC Edition GPU, 8GB DDR4 RAM. The software was Marvin Sketch for an academic license, Molegro Molecular Viewer, MGLTools 1.5.6, Discovery Studio Visualizer, Desmond Package for an academic license, and preADMET.
The materials were 2D structures of three thioureaderived compounds ( Figure 1) drawn using the Marvin Sketch for an academic license program. SIRT1 protein was downloaded through http://www.rscb.org/pdb/ with PDB code 4I5I [25].

Predictive Analysis of Pharmacokinetics, Toxicity and Lipinski's Rules of Five
The screening process of the test compounds was performed to predict the physicochemical properties of the test compounds. The screening was conducted to predict pharmacokinetic properties, toxicity, and compliance with Lipinski's Rules of Five [26,27].

Preparation of Ligands and Receptors
The ligands were drawn in two dimensions using Marvin Sketch software, then protonated to produce a pH 7.4, and then the file was saved in .mrv format. The next step was to determine the conformation of the structure by performing a conformational search with the MMFF94 force field setting. Then the file was saved in the form of mol2 [28].
The 3-dimensional structure of SIRT1 was downloaded from the Protein Data Bank website with the code PDB 4I5I and (6S)-2-chloro-5,6,7,8,9,10hexahydrocyclohepta[b]indole-6-carboxamide as a natural ligand. SIRT1 protein was separated from ligands and water molecules using Molegro Molecular Viewer software; consequently, only the SIRT1 monomeric protein structure was used as the target protein.

Validation of Molecular Docking Method
The method was validated by re-docking the natural ligand to the protein's active site with a grid box size of 40 x 40 x 40 Å and x, y, and z coordinates (42.194; −21,207; and 18,508 Å), and the value obtained was Root Mean Square Deviation. (RMSD) in angstroms (Å) [29]. The validation of the molecular docking method is declared valid if the RMSD value of ≤ 2 Å was obtained [30].

Docking of Test Ligand
Docking of test ligand was carried out for the same target protein on natural ligands using the Autodock program with the Lamarckian GA algorithm and without repetition with each docking yielding 100 conformations. The results seen at this stage were the value of binding free energy (BFE) or ∆G with units of kcal/mol and the value of the inhibition constant (KI) [28].

Docking Result Visualization
The interaction of ligands with receptors in 2D and 3D forms can be determined through docking results visualization. As a result, the position of the ligandreceptor conformation with the lowest BFE was obtained, and a complex file was created in the (.pdb) format for visualization using the Discovery Studio Visualizer software.

Molecular Dynamics Analysis
Molecular dynamics studies were carried out on one of the best compounds resulting from molecular docking and natural ligands for 30 ns using the Desmond software for academic license (DE Shaw Research, New York) program to examine the bond stability of the ligand-4I5I complex. The simulation system used a TIP3P water model and 0.15 M NaCl to mimic a physiological ionic concentration. Energy minimization was performed for 100 ps. The MD simulation was run for 30 ns at 300 K and standard pressure (1.01325 bar) in an orthorhombic box with buffer dimensions of 10 x 10 x 10 Å and an embellished NPT. The energy was recorded at 1.2 ps intervals. Protein-ligand complexes were neutralized by adding Na + or Cl − ions to balance the net charge system of MD simulation. The dynamic Nose-Hoover chain and Martyana-Tobias_klein algorithms were used to maintain the temperature of all MD systems at 300 K and a pressure of 1.01325 bar, respectively [31,32,33].

Physicochemical Screening of Test Compounds
Pharmacokinetic predictions have been performed using the PreADMET program with the parameters of Caco-2, human intestinal absorption (HIA), and plasma protein binding (PPB), as shown in Table 2. Caco-2 cells have a role as a drug transport pathway through the intestinal epithelium, a derivative of human colon adenocarcinoma. Based on Table 2, it can be seen that natural and test ligands have a moderate ability to penetrate cell membranes because their values are in the range of 4-70 [34]. A compound must have an excellent ability to penetrate cell membranes to reach the target tissue; hence it can provide good biological activity. HIA parameter relates to the ability of a compound to be absorbed in the body. Drug absorption is the passage of drugs from the application site to the circulatory system. The blood distributes drugs from the body surface or certain parts in the organs to all parts of the body. The results showed that the natural ligands and the test ligands had HIA values in the range of 70-100%, meaning that these compounds could be well absorbed in the body.
The protein can influence the drug's effectiveness in the blood plasma (Plasma Protein Binding). The weaker the drug binding, the more effectively the drug penetrates cell membranes. There are two forms of drugs in the blood: bound and unbound [35]. In Table 2, natural ligands had an excellent distribution ability in the body because they had a PPB of less than 90%, meaning they are weakly bounded. In comparison, the three test ligands had a PPB value of more than 90%, meaning that these ligands are highly bound to plasma proteins. Drugs bound to plasma proteins are inactive, while free or unbound drugs can bind to target receptors and provide biological activity.
Toxicity prediction was conducted to see the level of toxicity of natural and test ligands. The toxicity test used three parameters: the Ames test, carcino mouse, and carcino rat. The Ames test was carried out to assess the potential for mutagenicity that could act as a carcinogen. Ligand compounds as candidates for anticancer drugs in high concentrations are not guaranteed to be safe because they have mutagenic effects.
The carcino mouse and carcino rat parameters are used to determine whether the ligands can have a carcinogenic effect on rats and mice, the results showed that all the tested ligands did not have a carcinogenic effect on mice or rats. However, the test ligands based on predictions can have a carcinogenic effect on mice. As a result, natural and test ligands as candidates for anticancer drugs possess certain limitations. Besides, it is supported by Lipinski's Rule of Five rules as listed in Table 3. The lipophilicity or hydrophobicity of the ligand compound can be determined from the partition coefficient value (log P). A high log P value indicates that the molecule is hydrophobic. Also, highly hydrophobic molecules increase the risk of toxicity. The test ligand has a higher hydrophobic property because the log P value is more than 5. Hydrophobic compounds will be retained longer in the lipid bilayer because the layer has hydrophobic properties. Moreover, these compounds will be more widely distributed in the body; therefore, the selectivity of binding to target receptors is reduced [36].
The molecular weight of a compound can affect the drug distribution process by observing its ability to penetrate biological membranes through a passive diffusion process. Compounds with large molecular sizes have a low ability to penetrate biological membranes. The test ligands and natural ligands had a small molecular weight, making them easily penetrate biological membranes.
The parameters of hydrogen bond donor (less than 5) and hydrogen bond acceptor (less than 10) are related to the biological activity of the drug molecule. The presence of hydrogen bonds can affect the physicochemical properties of a compound. The number of hydrogen bond donors in a compound will bond with the solvent to form hydrogen bonds. Likewise, the hydrogen acceptors affect permeability by interacting favorably with a strong hydrogen-bonding solvent. Therefore, compounds that interact with polar solvents can reduce their ability to permeate the lipid bilayer.
The refractory molar values of all ligands met the requirements in the range of 40-130. This value measures the total polarizability of drug molecules that depend on the refractive index, temperature, and pressure.

Preparation of Ligands and Receptors
At this stage, protonation of the ligand occurs to adjust itself to a blood pH of around 7.4. At the same time, the conformational step was performed to obtain the molecule's position at the lowest or stable energy; thus, it can interact with the receptor's active site. There were ten conformations in the preparation, and one conformation with the lowest energy was continued for the docking process. Receptor preparation was aimed to remove the water content that could block the ligandreceptor interaction [26].

Molecular Docking Validation
In Figure 2, it can be seen that the results of the redocking validation of the SIRT1 protein with a grid box size of 40×40×40 Å and x, y, and z coordinates (42.194; −21,207; and 18,508 Å) against (6S)-2-chloro-5,6,7,8,9,10-hexahydrocyclohepta[b]indole-6carboxamide as a natural ligand obtained an RMSD value of 0.374 Å. The overlay image between the crystallographic structure and the re-docking results can be seen in Figure 2.
The valid receptor parameter (RMSD < 2 Å) showed similarity based on comparing the distance change and the best ligand conformation in the simulation results with the initial ligand distance and conformation.

Docking of Test Ligand
In docking of test ligand stage, the parameters used to determine the binding affinity value are analysis of binding free energy (∆G) and constant of inhibition (KI). The stability of the interaction between the ligand and the target receptor was predicted with the bond free energy analysis, which was indicated by a low ∆G value [37]. The results can be seen in Table 1, that compound b was predicted to have an activity value of −9.29 kcal/mol, with the resulting inhibition constant value being directly proportional to 0.156 µM. In contrast, the other test ligands had a higher ∆G value, as revealed by compound a of −9.07 kcal/mol and compound c of −8.68 kcal/mol. However, natural ligands were still lower at −9.74 kcal/mol. The low value of ∆G produced leads to a stable ligand-receptor conformation formed.  On the other hand, the ligand-receptor complex is less stable when the ∆G value is high. The more negative the ∆G value leads to the higher affinity of the proteinligand complex, the stable interaction between the compound and its receptor, and the better the compound activity. The ∆G results will be directly proportional to the KI value, where the KI value describes the ability of the compound to inhibit the enzyme. Figure 3 shows 2D interactions between the ligands and amino acid residues (proteins). The closer the distance between the ligand bonds to the amino acid residues on the receptor, the more the ligand bonds with these residues; thus, the more stable the resulting bond.

Visualization of the Docking Results
All thiourea derivatives have three hydrogen interactions with different amino acid residues. In contrast, the natural ligand (d) has four hydrogen bonds because there is an AspA:348 residue that hydrogen interacts with two different groups in the natural ligand. Hydrogen bonding is a desirable bond in ligand-receptor interactions because it can affect the physicochemical properties of compounds and their biological activities. The hydrophobic bond indicates the level of drug solubility in the cell membrane, which is predicted to be effectively bound to the receptor. The presence of hydrogen and hydrophobic bonds in thiourea derivatives and natural ligands is because these compounds have some of the same groups, including C=O, -NH, and a benzene ring.
The same amino acid residues on the test ligands and target receptors indicate that the test compounds interact stably with the receptor binding site and have competitive activity. The active site area is the site of protein binding to ligand compounds involving amino acid residues by forming a stable conformational structure from ligand-receptor interactions.

Dynamic Molecular Simulation Analysis
Molecular dynamics simulation was conducted on the best-bonded compound (compound b) and natural ligand compounds. Molecular dynamics simulations aimed to discern the stability of the bonded ligandprotein interactions from the docking results. Several factors that influence the results of molecular dynamics simulations include the conformation of compounds, water molecules, ions, cofactors, protonation of compounds, and the conformation of the solvent entropy [31]. Several research results to support the role and importance of molecular dynamics simulation for the refinement of molecular dynamics have been reported. The final structure of the molecular dynamics simulation of compound b showed a stable stereochemical geometry, as shown in the analysis of the Ramachandran plot of Figure 4.
Based on the Ramachandran plot, the percent residues in favored, allowed, and disallowed regions of the b-SIRT1 compound complex can also be observed in Figure 4. From Figure 4, it can be explained that the b-SIRT1 complex is gemlike residue in the favored region was around 84.4%, and there was no residue in the disallowed region. In conclusion, the quality of the b-SIRT1 complex was excellent.  Based on the Ramachandran plot, the percent residues in favored, allowed, and disallowed regions of the b-SIRT1 compound complex can also be observed in Figure 4. From Figure 4, it can be explained that the b-SIRT1 complex is gemlike residue in the favored region was around 84.4%, and there was no residue in the disallowed region. In conclusion, the quality of the b-SIRT1 complex was excellent.
In addition, the stability of the complex interaction of the b-SIRT1 compound can also be seen from the graph of the RMSD during a 30 ns simulation, as shown in Figure 5. Based on the RMSD graph, the complex b-SIRT1 and natural ligand-SIRT1 have a fluctuating RMSD from 0-14 ns, then after 14 ns until the end of the simulation, 30 ns has a reasonably stable RMSD. From the RMSD graph, it can be concluded that for 30 ns of molecular dynamics simulation, the complex with the most stable interaction stability is the b-SIRT1 complex compared to the natural ligand-SIRT1 complex. This result is supported by the average, minimum and maximum values of RMSD during the 30 ns simulation of each SIRTUIN1-ligand complex, as shown in Table 4. Besides RMSD, interaction stability was also evaluated using the RMSF graph, shown in Figure 6. According to the RMSF graph, the SIRT1-ligand can be used to prove the stability of the interaction. It is illustrated that the fluctuations of the residues have fluctuations in the same area. Still, the complex system of b-SIRT1 compounds has lower fluctuations than the natural ligand-SIRT1 if seen from the average RMSF of the compound-complex b-SIRT1 (1.163 Å) and natural ligand-SIRT1 (1.255 Å). The graph proves that the complex system of b-SIRT1 is a more stable interaction than the natural ligand-SITR1 complex system.
The RMSF graphic analysis can also identify amino acid residues that have contact/interaction with the ligand, as shown in Figure 7. Figure 7 shows the RMSF fluctuation and the residues that contact compound b.  Figure 7, it can be explained that the residues that have contact/interaction with compound b are Ala_262, Ser_265, Ile_270, Asp_272, Phe_273, Ile_279, Phe_297, Phe_305, Ile_316, Gln_345, Asn_346, Ile_347, His_363, Ile_411 , Val_412, Phe_413, Phe_414, Gly_440, Ser_441, Ser_442, Leu_443, Lys_444 and Val_445 (25 contact residues) both interacting through hydrogen bonds, hydrophobic bonds, ionic and water bridges, as shown in Figure 7 and Figure 8.  Figure 9 shows that the residues interact with the ligands at each pass. Some protein residues showed more than one specific contact with the ligand, which marked a darker orange color. Overall, six parameters were analyzed to explain the stability of compound b at the SIRT1 receptor during the 30 ns MD simulation, as depicted in Figure 10. Figure 10 shows that the simulation process of compound b RMSD fluctuated. In the early stages, fluctuations were observed from 0 to 17 ns, and after that, constant RMSD was observed during the whole simulation process. Fluctuations for gyration's radius were recorded up to 17 ns, and further, a stable conformation was obtained until the end of the simulation. The simulation radius of 30 ns of gyration ranged from 3.371 to 3.948 for compound b. The SASA plot revealed a fluctuating pattern for 3 ns and stabilized until the simulation was complete. MolSA and PSA plots for compound b up to 17 ns fluctuated but remained stable until the end of the 30 ns simulation. Meanwhile, the plot of the intramolecular H-bond of compound b, from the beginning to the end of the simulation, showed the intramolecular H-bond, which became stable after protein interactions.    Figure 10 shows that the simulation process of compound b RMSD fluctuated. In the early stages, fluctuations were observed from 0 to 17 ns, and after that, constant RMSD was observed during the whole simulation process. Fluctuations for gyration's radius were recorded up to 17 ns, and further, a stable conformation was obtained until the end of the simulation. The simulation radius of 30 ns of gyration ranged from 3.371 to 3.948 for compound b. The SASA plot revealed a fluctuating pattern for 3 ns and stabilized until the simulation was complete. MolSA and PSA plots for compound b up to 17 ns fluctuated but remained stable until the end of the 30 ns simulation. Meanwhile, the plot of the intramolecular H-bond of compound b, from the beginning to the end of the simulation, showed the intramolecular H-bond, which became stable after protein interactions.
3D visualization of compound b complex and natural ligand with SIRT1 after 30 ns simulation is illustrated in Figure 11.  Figure 10, it can be seen that the Phe_297 residue interacts with both natural ligands and with compound b following the results of the research in the reference that Phe_297 is one of the amino acid residues that interact with compound 415 ((6S)-2-chloro-5,6,7,8,9,10-hexahydrocyclohepta[b]indole-6carboxamide) [25].

Conclusion
From the molecular docking results, it was concluded that compound b ([2-(methylcarbamothioylcarbamoyl)phenyl]benzoate) had the lowest ∆G value of −9.29 kcal/mol with a KI value of 0.156 µM compared to the other two thioureas, supported by the result of molecular dynamics simulation for 30 ns. The active site amino acid residue that interacts with the ligand is PheA:297. In terms of pharmacokinetics and toxicity of compound b is better than natural ligands, it can be concluded that compound b is predicted to be conducted further research as an anticancer candidate through inhibition of SIRT1.