Available online on 15.12.2021 at http://jddtonline.info

Journal of Drug Delivery and Therapeutics

Open Access to Pharmaceutical and Medical Research

Copyright  © 2021 The  Author(s): This is an open-access article distributed under the terms of the CC BY-NC 4.0 which permits unrestricted use, distribution, and reproduction in any medium for non-commercial use provided the original author and source are credited

Open Access  Full Text Article                                                                                                                       Research Article 

Study of Molecular Docking, Molecular Dynamic and Toxicity Prediction of Several Quinoline Alkaloid Derivatives as a Bruton Tyrosine Kinase Inhibitor as Anti-Leukemia

Muttaqin, Fauzan Zein^1,2*; Restisari, Ikma Hanifah²; Muhammad, Hubbi Nasrullah³

¹ Faculty of Pharmacy, Bhakti Kencana University, Bandung, West Java, 40614, Indonesia

² Bandung School of Pharmacy, Bandung, West Java, 40614, Indonesia

³ School of Pharmacy, Bandung Institute of Technology, Bandung, West Java, 40132, Indonesia

INTRODUCTION

Cancer is a disease with a high prevalence and is the second-highest disease in the world after heart disease. In Indonesia, cancer is a dreaded disease because it is difficult to cure and causes many deaths. Cancer risk factors such as genetics, carcinogens, and lifestyle behaviors lead to nearly 1.4%, or about 347,792 of all cancer cases ¹.  In the previous research, we found some potential anticancer drug candidate compounds ². Some of the leukemia treatments currently used are surgery and drug treatment. The death rate from leukemia reaches 9.8%, and leukemia recurrence approaches 5.3%. This is due to the lack of available treatments. One of the leukemia treatments is the use of ibrutinib ³. Ibrutinib is the primary choice in chronic lymphocytic leukemia treatment and is one of the potent drugs in leukemia treatment that acts as a Bruton’s Tyrosine Kinase (BTK) inhibitor ^4,5. BTK plays an essential role in the B cell signaling process, especially in the B leukocytes cell proliferase ^6,7. Some studies suggest that other synthetic derivatives can inhibit BTK, namely GDC-0834. GDC-0834 has the same bonding position as ibrutinib (PCI-32765), occupying three inhibition positions on BTK receptor ⁸. Another study suggests that the presence of small molecules, called cinnoline, is known to have the potential to inhibit BTK by occupying one position of inhibitory bonding ⁹. Cinnoline structure has similarities with quinoline alkaloid compound, based on applying the Structure-Activity Relationship (SAR) principle ¹⁰. In this study, we conducted in silico study of some quinoline alkaloid derivatives as Bruton’s Tyrosine Kinase inhibitor by docking and MD for the development of antileukemia drugs.

MATERIALS AND METHODS

Docking, molecular dynamic, and ADMET studies were performed for all designed compounds by Gaussian 09, AutoDock4 and AutoDockTools4, GROMACS, and ADMET PredictorTM installed in a single machine running on a 3.4 GHz Intel (R) Core(TM) i7-6700 processor with 16 GB RAM, 4 GB VGA, and 500 GB hard disk with dual-boot Microsoft Windows 7 Pro 64 bit and Linux Ubuntu as the operating system.

Macromolecule preparation

Crystal structure of Bruton’s Tyrosine Kinase was obtained from the Protein Data Bank with PDB ID 4ZLY ¹¹. From the macromolecule, the small molecule (ligand) and water molecule were removed, and polar hydrogens and Kollman charges were added (Figure 1).resulting

Figure 1: (a) Chemical structure of ibrutinib ¹² (b) interaction of GDC-0834 into the crystal structure of Bruton’s tyrosine (PDB ID 4ZLY) ⁸.

Ligand preparation

The structures of twenty quinoline alkaloid derivative compounds ^13-15 were built using the ChemOffice suite ¹⁶. Geometry optimization and density functional theory (DFT) calculations were performed in Gaussian ¹⁷ using the Becke three-parameter Lee-Yang-Parr (B3LYP) functional at the 6-31G basis sets.

Determination of Physicochemical Properties Parameters

After the structure of the compounds in 2D and 3D was made, these compounds were measured for physicochemical properties parameters, including analysis of RO5 (BM, Log P, H-bond and Rotatable Bond) and determining the similarity of the structure of the compound with overlaying.

Molecular Docking

Each ligand molecule was prepared for docking using AutoDock Tools 1.5.6 ¹⁸. Hydrogen atoms were added, and partial charges of each atomBerberine resulted from the DFT calculations were incorporated. Grid maps were created by centering the grid box at the position of the natural ligand of each macromolecule with a spacing of 1.85 Å and size covering the binding cavity of each target. Lamarckian genetic algorithm and 100 docking runs were used for each simulation.

Molecular Dynamics Simulations

Ten ligands with the best docking score were chosen for further molecular dynamics (MD) study. MD simulations were carried out using the Gromacs 5.1.1 suite ¹⁹. The Amber99sb-ildn force field and the general AMBER force field (GAFF) were used to parameterize the atoms of the macromolecules and ligands, respectively. Energy minimization was carried out on the macromolecules in vacuum using the steepest descent algorithm. Then, the macromolecules were solvated with TIP3P water molecules in an octahedron box. Positive and negative ions were added to the system at a concentration of 0.15 N to neutralize all charges. Energy minimization was again performed on the macromolecule/solvent/ion system to release strains resulted from the solvation procedure; the steepest descent algorithm was used again. Next, the system was carefully heated to 310 K and pressurized to 1 atm using the constant-volume, constant-temperature (NVT) and constant-pressure, constant-temperature (NPT) ensembles. The Berendsen thermostat coupling was used to maintain the system temperature and pressure²⁰. The Particle Mesh Ewald (PME) method ²¹ with a cut-off value of 5.0 Å was used to compute long-range interactions. The system’s stability was evaluated by analyzing the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) of the protein backbones. A production simulation run was carried out on each macromolecule for ten ns. Analysis of the stability of ligand-protein interactions was performed by calculating the RMSD and RMSF values of the atoms at the protein binding sites throughout the simulation.

Toxicity Prediction

Toxicity prediction to ensure the safety of drug candidates was performed using ADMET Predictor 8.1 ²² from SimulationsPlus.inc including endocrine toxicity, maximum recommended therapeutic dose (MRTD), carcinogenicity and genotoxicity, cardiac toxicity, human liver adverse effects, acute rat toxicity, allergenic skin and respiratory sensitization, and reproductive toxicity.

RESULTS

Geometry Optimization

Geometry optimization is a way to obtain structures with minimum or lowest energy from a molecule. The structure with the lowest energy represents the most stable state ²³. In geometry optimization, the priority is how to find the most stable structure with low energy ²⁴. The lowest energy and geometry structure obtained can be used as a comparison in an experiment ²⁵. The geometry optimization was carried out using the Gaussian09 software to achieve the minimum energy and the most suitable/stable position for the molecular test compound for docking at the MCL-1 receptor. The method used in this study was the Density Functional Theory (DFT) method since this method has a high accuracy value, and the basis set used was the 6-31G with the hybrid functional B3LYP.

Bonding Energy (HOMO-LUMO)

The calculation of the HOMO-LUMO energy of the ligands (Table 1) showed the HOMO-LUMO energy values vary in each ligand. Then, the GAP or the difference between HOMO energy and LUMO energy was determined and was used as an indicator of stability.

Table 1: HOMO – LUMO energy and energy gap

Molecular Docking

Validation of Molecular Docking Methods

Previously, the validation of the docking method was conducted by re-docking between natural ligands and the Bruton’s Tyrosine target receptor (PDB ID 4ZLY). The Grid Box settings used were X = 40, Y = 40, and Z = 40. While Grid Spacing was 0.375 with Grid Center X = 20.060; Y = 5.856; and Z = 1.280.

Figure 2: (a) analysis of amino acid residue bonds (b) visualization of redocking results. green color = natural ligand; pink = ligands from re-docking.

Molecular Docking Simulation

The compounds were subsequently docked at the Bruton Tyrosine (PDB ID 4ZLY) target receptor with parameters adjusted to those obtained at the validation stage.   The binding free energy values and inhibition constant values (Ki) were obtained for each compound (Table 2).

TABLE 2: BINDING FREE ENERGY VALUES, INHIBITION CONSTANT VALUES (Ki) DAN BOND INTERACTION FOR EACH COMPOUND IN THE BRUTON’S TYROSINE RECEPTOR TARGET (PDB ID 4ZLY)

Molecular Dynamic Simulation

Molecular Dynamics is the most recognized simulation method and is considered very close to the original simulation. In various fields of molecular modeling, MD simulations show the flexibility of both ligands and proteins more effectively than other algorithms. The search strategy of new drug compounds by using random search to identify ligands' conformations, followed by MD simulations, is believed to ease researchers in developing new drugs. Molecular Dynamic was carried out using the Gromacs application. After topography preparation, the output file was obtained, which was then analyzed for several parameters. The time used was 10 ns, with the main parameters analyzed were the RMSD and RMSF values. In this test, it can be compared between the RMSD and RMSF values from before and after docking the compounds into the receptor (Table 3).

 Table 3: RMSD and RMSF MD parameter analysis of the test compound

Toxicity Prediction 

Toxicity prediction was carried out using ADMET Predictor ver 8.1 software which included endocrine toxicity, maximum recommended therapeutic dose (MRTD), carcinogenicity and genotoxicity, acute toxicity to mice, cardiac toxicity, side effects on the human liver, skin allergenic and respiratory sensitivity, and toxicity to the reproductive system. The results of toxicity analysis per test compound can be seen as follows:

Table 4: Summary of toxicity predictions analysis

Note: (KG) carcinogenicity, (M) mutagenicity, (E) endocrine, (J) heart, (K) skin, (P) respiratory, (R) reproduction, (T_R) risk of toxicity, (MRTD) maximum recommended therapeutic dose.

DISCUSSION

Parameters of physicochemical properties

From the analysis, it can be seen that the violations carried out were in the Rotatable-bond parameter whereby compound 17 had a value of 8; this would not be too influential because, based on the research, Rotatable-bond values of more than 10 correlates with a decrease in oral bioavailability in mice ²⁶. The mechanistic basis for the rotatable bond parameters was unclear because the number of rotatable bonds did not correlate with the in vivo clearance testing results in mice. However, the in-vitro experimental results allowed it since the average ligand affinity drops 0.5 kcal for each of the two rotatable bonds. Thus, the results of the physicochemical analysis based on the Lipinski Rule of Five found that all the compounds can be used as drug candidates.

Determination of chemical structure similarity

From overlaying process, the actual value obtained was less than 1 Angstrom (<1Å). This showed that the compounds had similarities with natural ligands of receptors. Thus 20 ligands of these compounds were expected to occupy the same active side at the same receptor. Likewise, the biological activity of the compounds was expected to have similar activities as biological activity possessed by natural ligands when occupying the active site of the same receptor.

Bonding Energy (HOMO-LUMO)

The result showed that the highest GAP value was possessed by Tetrahidro Berberine was 0.2032, and the lowest GAP value was possessed by Candidin, which was 0.0876. The large HOMO-LUMO gap indicated high kinetic stability and low chemical reactivity. This is attributed to the reactivity to interact with other compounds and more to the analysis of the formation and breaking of covalent bonds with the test compounds

Molecular Docking

Validation of Molecular Docking Methods

The validation results of the docking method obtained a valid value with the RMSD value that met the requirements of 1.03 Å and was considered valid because the RMSD value was less than 2 Å ²⁷. From the results of docking validation, the binding free energy values and inhibition constant values (Ki) were -11.35 kcal/mol and 0.00483 µM, respectively. The reference value of RMSD and the lowest free energy of the natural ligands that were redocked to the target receptor Bruton’s Tyrosine (PDB ID 4ZLY) were -14.13 kcal/mol and 1.85, respectively. The RSMD value obtained showed that the Docking validation method was considered valid because the RMSD value had fulfilled the requirements of <2Å, and the Docking method parameter could be used for the Docking process of the compounds. Apart from the RMSD values, the following analyzed was amino acid residues that interact with natural ligands (Figure 2a).

Molecular Docking Simulations

The compounds were subsequently docked at the Bruton Tyrosine (PDB ID 4ZLY) target receptor with parameters adjusted to those obtained at the validation stage.   The binding free energy values and inhibition constant values (Ki) were obtained for each compound (Table 2). The interactions formed between natural ligands and receptors were the bonds of amino acid residues, namely MET477, GLU475, LEU408, and THR474. Although the interaction of the compounds had a slightly different bond pattern with amino acid residues, it showed overall interaction patterns similar to natural ligands. From this docking stage, ten compounds were selected for the next MD stage. The selection of these ten compounds was based on the same interaction pattern with natural ligands and had a Ki value smaller than natural ligands, so it was suspected that these compounds could stabilize and have better biological activity among other test compounds. The ten compounds were S1, S2, S3, S4, S5, S6, S7, S8, S9, and S10.

Molecular Dynamic Simulations

Molecular Dynamics is the most recognized simulation method and is considered very close to the original simulation. In various fields of molecular modeling, MD simulations show the flexibility of both ligands and proteins more effectively than other algorithms. The search strategy of new drug compounds by using random search to identify ligands' conformations, followed by MD simulations, is believed to ease researchers in developing new drugs. Molecular Dynamic was carried out using the Gromacs application. After topography preparation, the output file was obtained, which was then analyzed for several parameters. The time used was 10 ns, with the main parameters analyzed were the RMSD and RMSF values. RMSD analysis describes changes or shifts of a molecular conformation. RMSD shows the stability of a molecule based on time. If the RMSD receptor value is higher than the RMSD of the compounds, it is estimated that the bond of the compounds is more stable. It can be concluded that the compounds S1, S4, S5, S6, S8, S9, and S10 are more stable than natural ligands with the RMSD value range is 0.075-0.5) (Table 3).

RMSF
 Root Mean Square Fluctuation (RMSF) analysis describes a change or shift from each atom/receptor amino acid residue that will produce flexibility from the receptor. The RMSF curve of the compounds was lower than the RMSF receptor curve. Thus, it can be stated that the compounds are better. To observe the fluctuations in the movement, the interaction of each atom was analyzed and then visualized according to the atomic analysis, where it was obtained the interaction movement per atom and visualized per residue that interacted just with the receptor graph and the receptor with the docking compound. It produced several compounds proven at the atomic level can stabilize receptors, including S1, S2, S5, S6, S8, S9, and S10, with fluctuations values in the range of 0.05-0.1.

Toxicity Prediction 

Toxicity prediction was carried out using ADMET Predictor version 8.1 software which included endocrine toxicity, maximum recommended therapeutic dose (MRTD), carcinogenicity and genotoxicity, acute toxicity to mice, cardiac toxicity, side effects on the human liver, skin allergenic and respiratory sensitivity, and toxicity to the reproductive system. The results of toxicity analysis per test compound can be seen as follows:

From the analysis of the results in Table 4, it can be seen that each compound has a different risk of toxicity, which can help researchers to develop new drugs by modifying existing structures to reduce the possible side effects. This data can be used as a reference if the effects have not been detrimental to the host in the long term, then the compounds can still potentially be developed into a new Leukemia drug.

Figure 3: Docking Visualization, RMSD, and RMSF of S8.

After analyzing the prediction of physicochemical properties, molecular docking, molecular dynamics, and toxicity prediction, it can be concluded that the compound, Yukositrin (S8), is a potent candidate compound to be developed to the next stage to become a Bruton’s Tyrosine Kinase inhibitor (Figure 3).

CONCLUSION

1. The results of molecular docking of 20 test compounds revealed that some compounds could be used as candidates for Bruton Tyrosine Kinase Inhibitors, including Candidine (S1), N-Methylaphinine (S2), Vebilocine (S4), Citracridone (S5), Yukocitrine (S8) compounds, Acrifoline (S13), Glycofolinine (S14), Roxiamine A (S16), Evoxine (S17) and Tripthantin (S20) which were predicted to have a higher potential than other compounds in inhibiting B cells signaling (MR).
2. The results of molecular dynamics of 10 compounds obtained seven compounds that have lower RMSD and RMSF values , namely Candidine, Vebilocine, Citracridone, Yukocitrine, Acrifoline, Roxiamine A, and Triptanthin. Thus it can be stated that the compounds have the best stability.
3. The compound, Yukositrin (S8), was a potential candidate to be developed to the next stage to become a Bruton’s Tyrosine Kinase inhibitor.

ACKNOWLEDGEMENT

The researchers would like to thank Sekolah Tinggi Farmasi Bandung (Bandung School of Pharmacy) for funding this research.

REFERENCES

1. Data and Information Center, the Ministry of Health Republic of Indonesia, 2015. Available at: http://www.pusdatin.kemkes.go.id/folder/view/01/structure-publikasi-pusdatin-profil-kesehatan.html [accessed 02 July 2017]

2. Muttaqin FZ, Fakih TM, and Muhammad HN: Molecular docking, molecular dynamics, and in silico toxicity prediction studies of coumarin, n-oxalylglycine, organoselenium, organosulfur, and pyridine derivatives as histone lysine demethylase inhibitors, Asian Journal of Pharmaceutical and Clinical Research 2017; 10(12):212-215. https://doi.org/10.22159/ajpcr.2017.v10i12.19348

3. Furman RR, Cheng S, Lu P, Setty M, Perez AR, Guo A, Racchumi J, Xu G, Wu H, Ma J, Steggerda SM, Coleman M, Leslie C, and Wang YL: Ibrutinib resistance in chronic lymphocytic leukemia. N Engl J Med. 2014 Jun 12; 370(24):2352-4. Epub 2014 May 28. https://doi.org/10.1056/NEJMc1402716

4. Brown JR, Barrientos JC, Barr PM, Flinn IW, Burger JA, Tran A, Clow F, James DF, Graef T, Friedberg JW, Rai K, and O'Brien S: The Bruton tyrosine kinase inhibitor ibrutinib with chemoimmunotherapy in patients with chronic lymphocytic leukemia. Blood. 2015 May 7; 125(19):2915-22. Epub 2015 Mar 9. https://doi.org/10.1182/blood-2014-09-585869

5. Ivanescu AM, Oprea M, Turbatu A, Colita A, and Lupu AR: Ibrutinib, a novel agent in relapsed or refractory chronic lymphocytic leukemia. Maedica (Buchar). 2014 Jun; 9(2):217-8.

6. Maas A and Hendriks RW: Role of Bruton's tyrosine kinase in B cell development. Dev Immunol. 2001; 8(3-4):171-81. https://doi.org/10.1155/2001/28962

7. Wiestner A: Targeting B-Cell receptor signaling for anticancer therapy: the Bruton's tyrosine kinase inhibitor ibrutinib induces impressive responses in B-cell malignancies. J Clin Oncol. 2013 Jan 1 https://doi.org/10.1200/JCO.2012.44.4281

31(1):128-30. DOI: 10.1200/JCO.2012.44.4281. Epub 2012 Oct 8. https://doi.org/10.1200/JCO.2012.44.4281

8. Young RM and Staudt LM: Ibrutinib treatment of CLL: the cancer fights back. Cancer Cell. 2014; 26(1):11-3. https://doi.org/10.1016/j.ccr.2014.06.023

9. Reiff SD, Mantel R, Smith LL, Greene JT, Muhowski EM, Fabian CA, Goettl VM, Tran M, Harrington BK, Rogers KA, Awan FT, Maddocks K, Andritsos L, Lehman AM, Sampath D, Lapalombella R, Eathiraj S, Abbadessa G, Schwartz B, Johnson AJ, Byrd JC, and Woyach JA: The BTK Inhibitor ARQ 531 Targets Ibrutinib-Resistant CLL and Richter Transformation. Cancer Discov. 2018 Oct; 8(10):1300-1315. Epub 2018 Aug 9. https://doi.org/10.1158/2159-8290.CD-17-1409

10. Thomas G: Fundamentals of Medicinal Chemistry. John Wiley & Sons, 2004

11. Smith CR, Dougan DR, Komandla M, Kanouni T, Knight B, Lawson JD, Sabat M, Taylor ER, Vu P, and Wyrick C: Fragment-Based Discovery of a Small Molecule Inhibitor of Bruton's Tyrosine Kinase. J. Med. Chem., 2015; 58(14):5437-5444 https://doi.org/10.1021/acs.jmedchem.5b00734

12. de Claro RA, McGinn KM, Verdun N, Lee SL, Chiu HJ, Saber H, Brower ME, Chang CJ, Pfuma E, Habtemariam B, Bullock J, Wang Y, Nie L, Chen XH, Lu DR, Al-Hakim A, Kane RC, Kaminskas E, Justice R, Farrell AT, and Pazdur R: FDA Approval: Ibrutinib for Patients with Previously Treated Mantle Cell Lymphoma and Previously Treated Chronic Lymphocytic Leukemia. Clin Cancer Res. 2015 Aug 15; 21(16):3586-90. https://doi.org/10.1158/1078-0432.CCR-14-2225

13. Michael JP: Quinoline, quinazoline and acridone alkaloids. Nat. Prod. Rep., 2000; 17:603-620. https://doi.org/10.1039/a904850b

14. Kayser O, Kiderlen AF, and Croft SL: Natural products as antiparasitic drugs. Parasitol Res. 2003 Jun; 90 Suppl 2:S55-62. Epub 2003 Feb 20. https://doi.org/10.1186/1475-2832-2-S1-S55

15. Rueping M, Antonchick AP, and Theissmann T: A highly enantioselective Brønsted acid catalyzed cascade reaction: organocatalytic transfer hydrogenation of quinolines and their application in the synthesis of alkaloids. Angew Chem Int Ed Engl. 2006 May 26; 45(22):3683-6. https://doi.org/10.1002/anie.200600191

16. Kerwin SM: ChemBioOffice Ultra 2010 Suite. J. Am. Chem. Soc., 2010; 132(7):2466-2467. DOI: 10.1021/ja1005306 https://doi.org/10.1021/ja1005306

17. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, et al: Gaussian 09, Revision A.02. Gaussian, Inc., Wallingford CT, 2016.

18. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, and Olson AJ: AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J Comput Chem. 2009 Dec; 30(16):2785-91. https://doi.org/10.1002/jcc.21256

19. Páll S., Abraham M.J., Kutzner C., Hess B., and Lindahl E. (2015) Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. In: Markidis S., Laure E. (eds) Solving Software Challenges for Exascale. EASC 2014. Lecture Notes in Computer Science, vol 8759. Springer, Cham https://doi.org/10.1007/978-3-319-15976-8_1

20. Berendsen HJC., Postma JPM, van Gunsteren WF, DiNola A, and Haak JR: Molecular-Dynamics with Coupling to an External Bath. Journal of Chemical Physics. 984; 81 (8):3684-3690. https://doi.org/10.1063/1.448118

21. Di Pierro M, Elber R, and Leimkuhler B: A Stochastic Algorithm for the Isobaric-Isothermal Ensemble with Ewald Summations for All Long-Range Forces. J Chem Theory Comput. 2015 Dec 8; 11(12):5624-37. https://doi.org/10.1021/acs.jctc.5b00648

22. Simulation plus., Inc.: ADMET PredictorTM. 42505 10th Street West, 2016.

23. Leach AR: Molecular Modelling: Principles and Applications. Prentice-Hall, Second Edition 2001.

24. Schlegel HB: Geometry optimization. Advanced Review 2011; 1(5):790-809. https://doi.org/10.1002/wcms.34

25. Cramer CJ: Essentials of Computational Chemistry: Theories and Models. John Wiley and Sons Second Edition 2004.

26. Lipinski CA: Drug-like properties and the causes of poor solubility and poor permeability. J Pharmacol Toxicol Methods. 2000 Jul-Aug; 44(1):235-49. https://doi.org/10.1016/S1056-8719(00)00107-6

27. Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, and Olson AJ: Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J Comp Chem. 1998; 19:1639-1662. https://doi.org/10.1002/(SICI)1096-987X(19981115)19:14<1639::AID-JCC10>3.0.CO;2-B