Int J App Pharm, Vol 15, Issue 1, 2023, 250-255Original Article

MOLECULAR DYNAMICS SIMULATIONS OF THE STK630921 INTERACTIONS TO INTERLEUKIN-17A

FRANSISCUS DEDDY RIANDONO, ENADE PERDANA ISTYASTONO*

Faculty of Pharmacy, Sanata Dharma University, Campus 3 Paingan, Maguwoharjo, Depok, Sleman, Yogyakarta 55282, Indonesia
*Email: enade@usd.ac.id

Received: 16 Sep 2022, Revised and Accepted: 19 Oct 2022


ABSTRACT

Objective: This research aimed to investigate the stability of the STK630921-Interleukin 17A (IL-17A) complex and to predict important residues that interact during molecular dynamics simulations.

Methods: Molecular docking simulations were performed, followed by molecular dynamics (MD) simulations and the free energy of binding calculations using YASARA-Structure. The identification of interacting residues was done using PyPLIF HIPPOS. Molecular docking simulations were performed on the IL-17A binding pocket with the compound 4-[({N-[(4-Oxo-3,4-dihydro-1-phthalazinyl) acetyl] alanyl} amino) methyl] cyclohexane carboxylic acid or known as STK630921. The best-docked pose was selected for the 50 ns MD simulations production run. The MD simulations snapshots were then analyzed to see the stability of IL-17A and for the identification of interacting residues, followed by Molecular Mechanics/Poisson–Boltzmann and surface area (MM/PBSA) analysis for the free energy of binding calculations.

Results: STK630921 is relatively able to stabilize IL-17A. Important interaction residues identified during the MD simulations were: Thr35(A), Pro37(A), Tyr62(A), Pro63(A)(B), Ile66(A)(B), Trp67(A), Ile96(A)(B), Val98(A)(B) and Val117(A)(B).

Conclusion: STK630921 disrupts the interaction of IL-17A to its receptor by binding and stabilizing IL17A.

Keywords: Interleukin 17-A, Molecular docking, Molecular dynamics simulations, MM/PBSA, YASARA-structure, PyPLIF HIPPOS


INTRODUCTION

Interleukin-17A (IL-17A) is a member of the IL-17 family, which consists of six cytokines (IL-17A to IL-17F), of which IL-17A and IL-17F are the main isoforms. All members of the IL-17 family are potent pro-inflammatory cytokines that are primarily secreted by Th17 cells and are also produced by other cells, including NK cells, macrophages, neutrophils, dendritic cells, and mast cells [1]. The pathogenicity of IL-17 has been found in several diseases, including psoriasis [2], rheumatoid arthritis [3], psoriatic arthritis [4], cancer [5], diabetes [6] and end-stage renal disease [7]. Serum levels of several Th17-associated cytokines, including IL-17A and IL-21, were found to be higher in diabetic patients compared to controls [8]. Plasma levels of CD4+, CCR5+, PD-1+, helper T cells, IL-6, and IL-17 in patients with diabetic nephropathy were also found to be greater than in healthy controls [9]. Furthermore, it was found that the blockade of IL-17A can reduce albuminuria and renal injury in diabetic nephropathy [10].

Diabetic nephropathy is a significant microvascular complication of diabetes. Approximately one-third of diabetic patients develop microalbuminuria after 15 y and no less than half develop nephropathy [11, 12]. Current management of diabetic nephropathy relies on optimal control of the renin-angiotensin system, using angiotensin-converting enzyme inhibitors or angiotensin II receptor drugs [13, 14]. The contribution of IL-17A to the pathogenesis of diabetic nephropathy is the background for the need to develop research related to IL-17A as a potential target for diabetic nephropathy therapy. Efforts to find small molecules that can inhibit IL-17A activity have been carried out in cases of intervertebral disc disease (IVD) in the form of compounds 4-[({N-[(4-oxo-3,4-dihydro-1-phthalazinyl)acetyl]alanyl}amino)methyl]cyclo-hexanecarboxylic or known as STK630921 [15]. STK630921, through its binding to the IL-17A receptor site can suppress the expression of COX‐2, IL‐6, MMP‐3, and MMP‐13 in nucleus pulposus (NP) cells.

This research aimed to study the stability of the IL-17A complex resulting from molecular docking with STK630921 using molecular dynamics (MD) simulation methods and Molecular Mechanics/Poisson–Boltzmann and surface area (MM/PBSA) calculations, followed by identification of its important interaction residues.

MATERIALS AND METHODS

Materials

The instrument used in this research is a personal computer with a specification of an Intel® CoreTM i5-10400 CPU @ 2.90GHz and 8 GB of RAM. The operating systems used are Windows 11 Pro-64-bit and Ubuntu 20.04 focal (on the Windows Subsystem for Linux). The software used is YASARA-Structure version 21.12.19 [16] and PyPLIF HIPPOS [17]. The material used is the crystal structure of human IL-17A obtained from the RCSB Protein Data Bank (https://www.rcsb.org/) namely: 5HI5 [18]. The ligand used is compound STK630921 [15]. Co-crystal ligand of the 5HI5 structure, i.e. (4S,20R)-7-chloro-N-methyl-4-{[(1-methyl-1H-pyrazol-5-yl)carbonyl]amino}-3,18-dioxo-2,19-diazatetracyclo [20.2.2.1~6,10~.1~11,15~]-octacosa-1(24),6(28), 7,9,11(27),12,14,22,25-nonaene-20-carboxamide or known as compound 63Q [18] is used as the control ligand.

Methods

The study started by a computer-aided visual analysis to identify binding pockets in the 5HI5. Subsequently, the IL-17A inhibitor molecule i.e., STK630921 was molecularly docked to the suspected 5HI5 binding pocket from the earlier visual analysis result.

The docking of the STK630921 ligand on the 5HI5 structure was conducted using the YASARA-Structure with the default macro command accessed from http://www.yasara.org/dock_run.mcr [19]. The docking result with the best score was then used for MD simulation studies.

MD simulations were performed on the 5HI5 structure with docked STK630921 and on the 5HI5 structure with co-crystal ligand 63Q. MD simulations were also carried out on the 5HI5 structure without ligand. YASARA-Structure is used to run MD simulations with modified macros from http://www.yasara.org/md_run.mcr [20]. The duration and save interval parameters were changed to run simulation durations of up to 50 ns, and simulation snapshot intervals were saved every 10 ps. The simulation cell was set as a periodic boundary with a simulated temperature of 310 K and density of 0.993 g/cm3.

The MD simulation results were analyzed using macros accessed from http://www.yasara.org/md_analyze.mcr by adding the parameter ligandsel='obj 2'. Free energy of binding analysis using the MM/PBSA was performed on each snapshot using a macro accessed from http://www.yasara.org/md_analyzebindenergy.mcr [21] with the method parameter set as Poisson-Boltzmann at a simulated temperature of 310 K.

Identification of interaction residues between the ligand and 5HI5 was performed using the PyPLIF HIPPOS software for the last 5 ns of the MD simulation. The pdb2plif. sh and md2plif. sh scripts [22] were used to produce outputs in the form of the names of amino acids and the types of interactions that are formed.

RESULTS AND DISCUSSION

The fluctuations in the root mean square deviation (RMSD) of the 5HI5 backbone atoms (RMSDBb), RMSD of the ligand movement (RMSDLm), and the Radius of Gyration (RoG) during the 50 ns simulation were used as indicators of the stabilization activity of the ligands. Results from the calculation of free energy of binding were used as a source of information on the binding affinity formed during the simulation.

The RMSDBb values provided information on the average distance between the backbone atoms in each snapshot compared to the initial position. RMSDBb of STK630921 for the last 30 ns is stable (fig. 1). Although the absolute value of RMSDBb from the complex with STK630921 was higher than complex without ligand, the comparison of ΔRMSDBb every 5 ns block of time showed that complex with STK630921 had a lower range value than complex without ligand and complex with the co-crystal ligand 63Q. In the last 5 blocks of time, ΔRMSDBb of the complex with STK630921 ligand is lower than RMSDBb of the complex without ligand as well as the complex with co-crystal ligand 63Q (fig. 2).

Fig. 1: The graphs of the RMSDBb of the 5HI5 without ligand (orange), the 5HI5-STK63021 complex (purple), and the 5HI5-63Q complex (cyan) vs simulation time

Fig. 2: The graphs of the ΔRMSDBb of the 5HI5 without ligand (green), the 5HI5-STK63021 complex (yellow), and the 5HI5-63Q complex (blue) vs simulation block of time

At the interval of 45 ns-50 ns, the value of ΔRMSD of complex with STK630921 was 1.369 Å, this was lower than ΔRMSD of complex without ligand, which is 1.872 Å, and complex with the 63Q which was 2.166 Å. Liu et al. [23, 24] found that the stability of the complex was achieved if the ΔRMSD value in the last 5 ns in a 10 ns simulation is less than 2 Å. It can be concluded that the STK630921 ligand does not destabilize the 5HI5 backbone atoms. On the contrary, the presence of this ligand can increase its stability.

Information about the movement of the ligands in its binding pocket during the MD simulation is objectively represented by the value of RMSDLm (fig. 3).

Fig. 3: The graphs of the RMSDLm of the 5HI5-STK63021 complex (yellow) and the 5HI5-63Q complex (blue) vs simulation time

STK630921 has less movement than the co-crystal ligand 63Q. The average movement distance of STK630921 during the simulation is 1.753 Å, which was smaller than the average movement of the 63Q, which was 4.455 Å. The range of movement of the STK630921 ligand was 2,991 Å, which is lower than the range of movement of the 63Q ligand, which was 6.648 Å. This shows that the movement of the STK630921 ligand was more stable than the movement of the 63Q ligand.

RoG provides information on the compactness of the structure during the simulation. The RoG value of the complex with the STK630921 ligand during the simulation had a lower mean compared to the other two complexes (fig. 4). This indicates that the STK630921 ligand can maintain the compactness of its complex structure.

The average free energy of binding of the complex with STK630921 was lower than the average value of the complex with co-crystal ligand 63Q (fig. 5). This indicated that the binding affinity of the complex with co-crystal ligand 63Q was better than that of the complex with STK630921. The complex with STK630921 in absolute terms had lower free energy of binding value than the complex with 63Q but referring to fluctuations in the free energy of binding per 5 ns block of time, it was found that in the last 6 blocks of time from 20 ns to 50 ns, the complex with STK630921 had lower fluctuations than the complex with 63Q (fig. 6). This indicates that the STK630921 ligand can maintain the stability of its free energy of binding in the 5HI5 structure.

Protein-ligand interaction analysis was performed by utilizing the ability of PyPLIF HIPPOS to identify the presence or absence of interaction (fig. 7 and fig. 8). PyPLIF HIPPOS was originally designed to identify the outputs of molecular docking simulations, namely AutoDock Vina [25] and PLANTS [26]. In its development, this tool has also been successfully used to identify the output of the MD simulation [27]. Therefore, it can be used to cover the weakness of molecular docking simulations that treat proteins and ligands as rigid entities without considering external forces such as temperature and pressure according to the real nature of biomolecules in the body that are always dynamic. Identification was carried out in the last 5 ns because at that interval, a stable complex had occurred, as indicated by the RMSDBb value. The summary of results and the visualization (table 1 and fig. 9) show that the type of interaction formed is mostly hydrophobic. Besides that, there are also some aromatic interactions which occur in the region between monomers A and B of the 5HI5 structure.

This study provides new information that the STK630921 compound is not only active at the IL-17A receptor site as previously found by Suyama et al. [15] but is also active at the IL-17A molecular site. This study addressed some of the questions raised by Lavoz et al. [28] regarding the discovery of new drugs targeting IL-17A for diabetic nephropathy. The design of new drugs can employ the molecular determinant template of interaction residues that have been successfully identified in this study.

Table 1: Interacting residues identified by PyPLIF HIPPOS

Interacting residue Interacting type Interaction percentage for 5HI5-63Q complex Interaction percentage for 5HI5-STK630921 complex
Thr35 A hydrophobic 0.02% 0.17%
Pro37 A hydrophobic 2.25% 4.50%
Tyr62 A aromatic (edge to face) 6.61% 0.04%
Pro63 A hydrophobic 2.23% 0.39%
Pro63 B hydrophobic 2.23% 0.39%
Ile66 A hydrophobic 6.74% 7.80%
Ile66 B hydrophobic 6.74% 7.46%
Trp67 A hydrophobic 5.67% 19.32%
Ile96 A hydrophobic 1.32% 9.68%
Ile96 B hydrophobic 1.32% 9.73%
Val98 A hydrophobic 0.08% 0.13%
Val98 B hydrophobic 0.08% 0.13%
Val117 A hydrophobic 0.99% 9.13%
Val117 B hydrophobic 0.99% 9.30%

Note: The letters A and B in the interacting residue column indicate monomers A and B of the 5HI5 structure, respectively.

Fig. 4: The graphs of the RoG of the 5HI5 without ligand (green), the 5HI5-STK63021 complex (yellow), and the 5HI5-63Q complex (blue) vs simulation time

Fig. 5: The graphs of the free energy of binding of the 5HI5-STK63021 complex (yellow), and the 5HI5-63Q complex (blue) vs simulation time

Fig. 6: The graphs of the free energy of binding fluctuations of the 5HI5-STK63021 complex (yellow), and the 5HI5-63Q complex (blue) vs. simulation block of time

Fig. 7: Interacting residues heatmap of the 5HI5-STK630921 complex during the last 5 ns simulation. Red indicates the presence of interaction, and green indicates the absence of interaction

Leu97A Hydrophobic

Leu97 B Hydrophobic

Leu26 A Hydrophobic

Ile66 B Hydrophobic

Ile66 A Hydrophobic

Tyr62 B Aromatic (edge to face)

Tyr62 A Aromatic (edge to face)

Leu99 A Hydrophobic

Leu99 B Hydrophobic

Trp67 A Hydrophobic

Leu112 A Hydrophobic

Leu112 B Hydrophobic

Trp67 B Hydrophobic

Pro37 A Hydrophobic

Pro63 A Hydrophobic

Pro63 B Hydrophobic

Ile96 A Hydrophobic

Ile96 B Hydrophobic

Val117 A Hydrophobic

Val117 B Hydrophobic

Met23 A Hydrophobic

Pro37 B Hydrophobic

Tyr62 B Aromatic (face to face)

Tyr62 A Aromatic (face to face)

Thr33 A Hydrophobic

Asn34 A Hydrophobic

Tyr62 B Hydrophobic

Tyr62 A Hydrophobic

Val98 B Hydrophobic

Val65 A Hydrophobic

Val98 A Hydrophobic

Val65 B Hydrophobic

Glu95 A Hydrophobic

Glu95 B Hydrophobic

Asn34 A H-Bond (residue as the donor)

Thr35 A Hydrophobic

Fig. 8: Interacting residues heatmap of the 5HI5-63Q complex during the last 5 ns simulation. Red indicates the presence of interaction, and green indicates the absence of interaction

Fig. 9: Close view of interacting residues (in surface representation) at the binding pocket of IL-17A. Red and blue sticks are STK630921 and 63Q pose at the final MD simulation snapshot

CONCLUSION

MD simulations for 50 ns showed that the STK630921 compound was able to stabilize IL-17A. The important interaction residues that were successfully identified were: Thr35 (A), Pro37 (A), Tyr62 (A), Pro63 (A)(B), Ile66 (A)(B), Trp67 (A), Ile96 (A)(B), Val98 (A)(B), and Val117(A)(B).

FUNDING

This research was funded by the Directorate of Research, Technology and Community Services, the Directorate General of Higher Education, Research, and Technology, the Indonesian Ministry of Education, Culture, Research, and Technology (No. 157/E5/PG.02.00. PT/2022).

AUTHORS CONTRIBUTIONS

E. P. I. conceptualized the research project. F. D. R conducted all the computational simulations and initiated the original draft of the manuscript. E. P. I. reviewed and edited the manuscript. All authors have given approval for the final version of the manuscript.

CONFLICT OF INTERESTS

All authors have none to declare

REFERENCES

  1. Ma J, Li YJ, Chen X, Kwan T, Chadban SJ, Wu H. Interleukin 17A promotes diabetic kidney injury. Sci Rep. 2019;9(1):2264. doi: 10.1038/s41598-019-38811-4, PMID 30783187.

  2. Blauvelt A, Chiricozzi A. The immunologic role of IL-17 in psoriasis and psoriatic arthritis pathogenesis. Clin Rev Allergy Immunol. 2018;55(3):379-90. doi: 10.1007/s12016-018-8702-3, PMID 30109481.

  3. Robert M, Miossec P. IL-17 in rheumatoid arthritis and precision medicine: from synovitis expression to circulating bioactive levels. Front Med (Lausanne). 2018;5:364. doi: 10.3389/fmed.2018.00364. PMID 30693283.

  4. Syngle A, Verma I, Kaur S, Syngle T. Interleukin-17 inhibition with secukinumab improves sudomotor dysfunction in psoriatic arthritis. Int J Pharm Pharm Sci. 2018 Mar 1;10(3):167-9. doi: 10.22159/IJPPS.2018V10I3.22660.

  5. Hyun YS, Han DS, Lee AR, Eun CS, Youn J, Kim HY. Role of IL-17A in the development of colitis-associated cancer. Carcinogenesis. 2012 Apr;33(4):931-6. doi: 10.1093/carcin/bgs106, PMID 22354874.

  6. Kuriya G, Uchida T, Akazawa S, Kobayashi M, Nakamura K, Satoh T. Double deficiency in IL-17 and IFN-γ signaling significantly suppresses the development of diabetes in the NOD mouse. Diabetologia. 2013 Aug;56(8):1773-80. doi: 10.1007/s00125-013-2935-8, PMID 23699989.

  7. Coto E, Gomez J, Suarez B, Tranche S, Diaz Corte C, Ortiz A. Association between the IL17RA rs4819554 polymorphism and reduced renal filtration rate in the Spanish RENASTUR cohort. Hum Immunol. 2015 Mar;76(2-3):75-8. doi: 10.1016/j.humimm.2015.01.027. PMID 25636567.

  8. Baharlou R, Ahmadi Vasmehjani A, Davami MH, Faraji F, Atashzar MR, Karimipour F. Elevated levels of T-helper 17-associated cytokines in diabetes type i patients: indicators for following the course of disease. Immunol Invest. 2016;45(7):641-51. doi: 10.1080/08820139.2016.1197243, PMID 27611173.

  9. Zhang N, Tai J, Qu Z, Zhang Z, Zhao S, He J. Increased CD4+CXCR5+T follicular helper cells in diabetic nephropathy. Autoimmunity. 2016 Sep;49(6):405-13. doi: 10.1080/08916934.2016.1196677, PMID 27477820.

  10. Lavoz C, Matus YS, Orejudo M, Carpio JD, Droguett A, Egido J. Interleukin-17A blockade reduces albuminuria and kidney injury in an accelerated model of diabetic nephropathy. Kidney Int. 2019;95(6):1418-32. doi: 10.1016/j.kint.2018.12.031. PMID 30982673.

  11. Gheith O, Farouk N, Nampoory N, Halim MA, Al-Otaibi T. Diabetic kidney disease: worldwide difference of prevalence and risk factors. J Nephropharmacol. 2016;5(1):49-56. PMID 28197499.

  12. Cheng HT, Xu X, Lim PS, Hung KY. Worldwide epidemiology of the diabetes-related end-stage renal disease, 2000-2015. Diabetes Care. 2021;44(1):89-97. doi: 10.2337/dc20-1913, PMID 33203706.

  13. Coates PT, Wong G. Current controversies in nephrology-how to cross-match for transplantation? Kidney Int. 2020;97(4):662-3. doi: 10.1016/j.kint.2020.02.002.

  14. Dhodi JB, Mestry SN, Juvekar AR. Diabetic nephropathy-genesis, prevention and treatment. Int J Pharm Pharm Sci. 2014 Sep 1:42-7.

  15. Suyama K, Sakai D, Hirayama N, Nakamura Y, Matsushita E, Terayama H. Effects of interleukin-17A in nucleus pulposus cells and its small-molecule inhibitors for intervertebral disc disease. J Cell Mol Med. 2018;22(11):5539-51. doi: 10.1111/jcmm.13828, PMID 30207057.

  16. Krieger E, Vriend G. New ways to boost molecular dynamics simulations. J Comput Chem. 2015 May 15;36(13):996-1007. doi: 10.1002/jcc.23899, PMID 25824339.

  17. Istyastono EP, Radifar M, Yuniarti N, Prasasty VD, Mungkasi S. PyPLIF HIPPOS: A molecular interaction fingerprinting tool for docking results of AutoDock Vina and PLANTS. J Chem Inf Model. 2020 Aug 24;60(8):3697-702. doi: 10.1021/acs.jcim.0c00305. PMID 32687350.

  18. Liu S, Dakin LA, Xing L, Withka JM, Sahasrabudhe VPV, Li W. Binding site elucidation and structure-guided design of macrocyclic IL-17A antagonists. Sci Rep. 2016;6:30859. doi: 10.1038/srep30859, PMID 27527709.

  19. Krieger E. YASARA MACRO-docking a ligand to a receptor. Available from: http://www.yasara.org/dock_run.mcr. [Last accessed on 15 Mar 2022]

  20. Krieger E. YASARA MACRO-running an accurate molecular dynamics simulation in water with slow, normal or fast speed. Available from: http://www.yasara.org/md_run.mcr. [Last accessed on 15 Mar 2022]

  21. Krieger E. YASARA MACRO-analyzing the ligand binding energy during a molecular dynamics simulation. Available from: http://www.yasara.org/md_analyzebindenergy.mcr. [Last accessed on 15 Mar 2022]

  22. Istyastono EP, Octa Riswanto FDO. Molecular dynamics simulations of the caffeic acid interactions to dipeptidyl peptidase IV. Int J App Pharm. 2022 Jul 7:274-8. doi: 10.22159/ijap.2022v14i4.44631.

  23. Liu K, Kokubo H. Exploring the stability of ligand binding modes to proteins by molecular dynamics simulations: a cross-docking study. J Chem Inf Model. 2017 Oct 23;57(10):2514-22. doi: 10.1021/acs.jcim.7b00412. PMID 28902511.

  24. Liu K, Watanabe E, Kokubo H. Exploring the stability of ligand binding modes to proteins by molecular dynamics simulations. J Comput Aided Mol Des. 2017;31(2):201-11. doi: 10.1007/s10822-016-0005-2, PMID 28074360.

  25. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. 2010 Jan 1;31(2):455-61. doi: 10.1002/JCC.21334, PMID 19499576.

  26. Korb O, Stutzle T, Exner TE. Empirical scoring functions for advanced Protein-Ligand docking with PLANTS. J Chem Inf Model. 2009;49(1):84-96. doi: 10.1021/ci800298z, PMID 19125657.

  27. Perdana Istyastono E, Gani MR. Identification of interactions of ABT-341 to dipeptidyl peptidase IV during molecular dynamics simulations. J Farmasi Galenika (Galenika Journal of Pharmacy) (e-Journal). 2021 Oct 1;7(2):91–8. doi: 10.22487/J24428744.2021.

  28. Lavoz C, Rayego Mateos S, Orejudo M, Opazo RL, Marchant V, Marquez Exposito L, Tejera Munoz A, Navarro Gonz JF. Could IL-17A be a novel therapeutic target in diabetic nephropathy. J Clin Med. 2020 Jan;9(1):272. doi: 10.3390/jcm9010272.