QSAR Study and Molecular Design of Isoquinolone Derivative JNK1 Inhibitors

TONG Jian-Bo XIAO Xue-Chun

LUO Dinga,b XU Hai-Yina,b WANG Jiea,b

a (College of Chemistry and Chemical Engineering,Shaanxi Uniνersity of Science and Technology, Xi'an 710021, China)

b (Shaanxi Key Laboratory of Chemical Additiνes for Industry, Xi'an 710021, China)

ABSTRACT JNK1 is a drug target for the treatment of type 2 diabetes, and it plays a key mediator role in metabolic disorders. In this paper, holographic quantitative structure-activity relationship (HQSAR)technology and Topomer comparative molecular field analysis (Topomer CoMFA) technology are used to analyze the quantitative structure-activity relationship (QSAR) of 39 isoquinolone derivatives. The cross validation correlation coefficient (q2) is 0.696 (Topomer CoMFA) and 0.826 (HQSAR), and the non-cross validation correlation coefficient (r2) is 0.935 (Topomer CoMFA) and 0.987 (HQSAR). The results showed that the models have good stability and predictive ability. The Topomer search module was applied to define high contribution fragments in the ZINC database, designing 20 new isoquinolone compounds with theoretically high inhibitory activity. The molecular docking was carried out to explore the interaction between the ligand and target JNK1 protein. This study can provide a theoretical basis for the design of new JNK1 inhibitors.

Keywords: Topomer CoMFA, HQSAR, molecular docking, isoquinolone derivatives, molecular design

1 INTRODUCTION

JNK is also known as stress-activated protein kinase[1].The activation of JNK may trigger the phosphorylation of many transcription factors related to cell apoptosis,thereby regulating the signal transduction of various extracellular stimuli from environmental stress or inflammatory cytokines[2]. JNK protein kinase is encoded by three genes, JNK1, JNK2 and JNK3. JNK1 and JNK2 genes are widely expressed throughout the body, while JNK3 has restricted expression, which is only expressed in the brain, heart, and testes[3].

In a study, it is shown that the JNK gene plays a vital role in the metabolic changes and inflammation associated with a high-fat diet lifestyle[4]. High-fat diet can induce type 2 diabetes, which is also commonly known as obesity.The patients with type 2 diabetes caused by obesity are resistant to insulin[5]. When driven by a high-fat diet,JNK1 gene can induce inflammation that leads to diabetes.Therefore, the development of new JNK1 inhibitors has great significance in the treatment of type 2 diabetes[6].

In recent years, computational approaches based on mathematics and chemo-informatics have been widely used to discover novel candidate drugs, such as quantitative structure-activity relationship (QSAR) and molecular docking. Mathematical model of QSAR is used to study the quantitative relationship between the molecular structure and biological activity of compounds. The methods used in this article are as follows: (1) HQSAR[7]is a 2D-QSAR technology. It can determine the relationship between biological activity and structural fragments and realize the ability of regulating molecular arrangement and conformation by converting the chemical structure into corresponding molecular holograms.According to the principle of holographic phase diagram,it can accurately determine the effect of any group or atom in the drug that reflected the drug activity. (2) Topomer CoMFA[8]is the second-generation CoMFA technology. It is the application of topology technology in CoMFA,which overcomes the shortcomings of the first generation of manually superimposing molecules, and greatly accelerates the research speed of drug design. Molecular docking is an effective tool to calculate detailed ligandprotein interactions for further understanding of the binding mode.

In order to design novel JNK1 inhibitors with higher inhibitory activity, a combined QSAR and molecular docking strategy was taken based on a series of isoquinolones analogs. The contour maps and color code maps are provided from QSAR studies to analyze the structural requirements of JNK1 inhibitory. Then, docking studies were carried out for revealing the binding affinity and display of ligand-protein interactions. This work provides a reference for the design of promising JNK1 inhibitors for treating type 2 diabetes.

2 RESEARCH METHODS AND PROCEDURES

2. 1 Data source and preparation

The 39 isoquinolone JNK1 inhibitors were selected from the literature[9]. The molecular structure and pIC50values(-logIC50) of the 39 isoquinolone derivatives are shown in Table S1. According to the molecular structure and activity value range, 9 representative compounds are selected as the test set, and the remaining 30 compounds are used as the training set, so as to ensure the diversity of the compound structure and activity value of the test and training sets. Among them, the training set compounds are used to build the model, and the test set compounds are used to verify the external prediction function of the model.

The three-dimensional structures of the 39 compounds were constructed using the SKETCH MOLECULE drawing module in SYBYL2.0-X. The Tripos force field and Powell energy gradient method were applied to optimize compounds to obtain the lowest energy conformation. The Gasteiger-Huckel charge was added with the maximum number of iterations is 1000. The energy convergence gradient value is set to 0.005 kcal/mol, and the rest parameters are the default values of SYBYL2.0-X[10].

2. 2 Construction of the HQSAR model

The construction of the HQSAR model includes three main steps[11]: Firstly, the molecule under study is cut into fragments of appropriate size. Molecular structure fragments are mainly determined by two parameters:fragment size and fragment distinction; the second step is to encode molecular fragments into sub-holograms; after obtaining molecular holograms, partial least-squares regression was used to establish a quantitative relationship model between molecular holograms and compound properties.

2. 3 Construction of Topomer CoMFA model

The operating process of Topomer CoMFA is roughly as follows[12]: A three-dimensional model of Topomer is built based on the molecular structure fragments of the training set, and the breaking method of the compound structure would be confirmed through the structural characteristics of the compound. Topomer CoMFA uses topological molecular technology to cut the ligand molecule into two or more small fragments, while retaining the common skeleton as much as possible. This method automatically constructs a standard topology 3Dmodel for each fragment and calculates the properties of the electrostatic field and the three-dimensional field of each molecular fragment. The finally obtained parameters are independent variables, and the biological activity value is the dependent variable. The partial least-squares method is used to establish the relationship between compound activity and molecular field characteristics to generate a 3D-QSAR model.

2. 4 Molecular screening

Virtual screening is a method for lead compound discovery in drug development. This study uses the Topomer Search technology in SYBYL2.0-X for molecular virtual screening, which can be used for lead compound optimization and backbone transitions based on ligands.Topomer Search technology can be used in conjunction with ligand-based molecular docking[13], and can also be used as a tool for preliminary screening with uncer- tainly receptor condition. Its basic principle is to use theRgroup in the Topomer CoMFA model as a template, searching for molecular fragments with high similarity in the compound database through Topomer similarity comparison, and selecting fragments with higher activity contribution values through certain criteria. New com- pounds have the theory higher biological activity could be obtained by searching the combination of fragments and basic skeletons.

2. 5 Molecular docking

Molecular docking[14]is to design novel drugs through the characteristics of receptors and the interaction between receptors and drug molecules. It is a theoretical simulation method to predict the binding mode and affinity through geometric, energy matching, and recognition of ligands and macromolecular proteins[15]. In this article, the Surflex-dock module in the SYBYL2.0 software package is selected to perform simulated docking[16], avoiding the tedious process of finding active sites. This process compiles a specific docking method and extracts the ligand bound at the original site to generate a docking protocol.Then the target ligand is docked to the previous specific site. The docking results of the original ligand can also be compared with potential ligand molecules to further screen the more suitable inhibitors.

3 EXPERIMENTAL RESULTS AND DISCUSSION

3. 1 HQSAR results and analysis

HQSAR modeling is performed on 30 training set compounds to establish a JNK1 inhibitor model. The HQSAR model can optimize the molecule by changing the fragment holographic length, fragment size and different fragment distinguishing parameters. In this experiment, the default fragment size is 4～7, all hologram lengths are selected, through combining different fragment discrimination parameters. The best model (A, C and Ch) shows statistical parameters that the cross-validation coefficientq2is 0.764, the non-cross validation coefficientr2is 0.933,and the best holographic length HL is 151. The established models are shown in Table S2.

In Table S2, the molecular fragment lengths are set to be the same, all of which are 4～7, in order to design and build a better model. Different fragment lengths are set under the conditions of selecting the best fragment discrimination parameter A/C/Ch. The experimental results are shown in Table S3, showing the best established HQSAR model gave the fragment length to be 9～10. The cross-validation coefficientq2of statistical parameters is 0.826, the non-cross-validation coefficientr2is 0.987, the best holographic length HL is 151, the standard deviation SEE is 0.080, and the number of principal componentsNis 6.

Compound No. 12 has the highest molecular activity, so it is selected as a template molecule in HQSAR analysis.The HQSAR color code map of the best model shown in Fig. 1 is used to further analyze the influence of each atom on the compound molecule. The green and blue parts in the figure indicate that the fragments have a higher contribution to the activity of the compound, while the white,orange, and red ones suggest no special contribution, so the fragments should be replaced by a better one. This provides a strong basis for the design of new JNK1 inhibitors with higher activity.

Fig. 1. HQSAR model color code map of compound 12

3. 2 Topomer CoMFA results and analysis

The most active molecule No. 12 is used as a template,combined with the analysis results of the HQSAR model,by increasing the contribution value of its white group.The activity of the compounds can be improved, so all the compounds are cut according to the cutting method shown in Fig. 2, and the molecules that are not automatically cut are manually cut. The cut compound is divided intoRa(blue),Rb(red) and common skeleton (green). By using Topomer CoMFA method to analyze the relationship between the structure and activity of the compound JNK1 inhibitor, the results of each parameter are shown in Table 1.

Fig. 2. Template molecule segmentation method and fragment schematic diagram

Table 1. Results of Topomer CoMFA Model

The obtained data show that theq2of the established model is 0.696 andr2is 0.935. It shows that the model constructed by this method is an ideal Topomer CoMFA model, and its statistical results have good predictive ability. At the same time, the estimated standard errorSEEof the predicted value is 0.171,q2stderr is 0.37,r2stderr is 0.17, andFis 207.892, which further shows that the constructed model has high reliability.

The contour map of the Topomer CoMFA model shown in Figs. 3a and 3b is the steric contour map and electrostatic contour map ofRagroup, respectively; Figs. 3c and 3d are the steric contour map and electrostatic contour map ofRbgroup, respectively. In Figs. 3a and 3c, the green part indicates that increasing the volume of the substituent can increase the activity of the compound, and the yellow part indicates that reducing the volume of the substituent can increase the activity of the compound. For example,compared molecule No. 02 (pIC50= 7.000) with No. 12(pIC50= 8.482), in theRagroup position, the No. 12 molecule selected a larger substituent on the No. 3 position of the benzene ring, and its activity was significantly improved. In Figs. 3b and 3c, the red area indicates that the introduction of high electronegative substituents at this position is conducive to the improvement of the activity of the compound, and the blue area suggests that introducing electropositive substituents or low electronegative substituents of the group helps improve the activity, such as No. 21 (pIC50= 5.854) and No. 23 (pIC50= 6.347)molecules. Molecule No. 23 increases the methyl group but removes the carboxyl hydroxyl group, thus improving its activity.

Fig. 3. 3D contour map of Topomer CoMFA model based on template compound 12

3. 3 Comparative analysis of the results of the HQSAR and Topomer CoMFA models

After calculation, the constructed HQSAR model is compared with the data of the Topomer CoMFA model. In the former, based on the results of the training set, the best fragment discrimination parameter A/C/Ch is selected, the fragment length is 9～10, the training set is verified by HQSAR, and the test set is predicted for pIC50. In the Topomer CoMFA model, based on the results constructed from the training set, Topomer CoMFA verification is performed on the training set and the prediction is made on the test set. The comparison of the predicted activity value and the residual value of the Topomer CoMFA and HQSAR models is shown in Table S4.

Fig. 4 is a linear regression correlation diagram between the experimental and predicted values of the training and test sets. It can be seen from the figure that all samples in the linear regression diagram are evenly distributed near the 45° line. As can be obtained from Table S4, the residual value of two models is also around 0, which proves that both have good predictive capabilities.

Fig. 4. (a) Regression analysis diagram of Topomer CoMFA model, (b) Regression analysis diagram of HQSAR model

3. 4 Molecular design results

According to the analysis results of the HQSAR model and the Topomer CoMFA model, the No. 12 template molecule with the highest activity is cut to obtain the corresponding molecular fragments, and then the Topomer similarity and threshold are compared and scored, and the obtained model is used to predict the contribution to the activity. Topomer search is used to perform a conformation search in the ZINC database, andRaandRbare used for searching and screening in the compound database. Among the molecular fragments[17]search process, theRgroup activity contribution value is selected to be greater than theRaandRbstructural fragments, and the value of Topdist is between 170 and 185. After comparison, the selectedRamolecular fragments have no higher contribution to the activity of the compound molecules than theRamolecular fragments of the No. 12 template molecule. Therefore, the new molecules designed continue to use theRamolecules of the No. 12 template molecule fragment. The selectedRbmolecular fragments have a significant increase in the contribution value of the compounds. 20 molecular fragments with better parameters are selected from them to design new molecules, as shown in Table S5.

As shown in Table S5, the predicted activity values(pIC50Pred) of the 20 newly designed molecules based on Topomer CoMFA results are all greater than the experimental activity value of the No. 12 template molecule, and theRbgroup structure of the newly designed molecules conforms to the results of Topomer CoMFA analysis.

3. 5 Molecular docking result analysis

The protein macromolecule used in the molecular docking[18]in this experiment is 4QTD[19]. Before the molecular docking, protein preparation work should be performed, delete all water molecules, other unnecessary ligands and terminal residues, extract the original ligand,and perform pretreatments of the protein such as hydrogenation and charge. Then the eutectic ligand in the protein is extracted from the crystal structure, then the eutectic ligand is re-inserted into the crystal structure through the Surflex-Dock molecular docking protocol, and the original ligand is used as a reference to verify the reliability of the docking. The superimposed conformation of the docking ligand and the original ligand is shown in Fig. 5, where the green ligand stick represents the original reference ligand, and the purple stick shows the respliced ligand. It can be seen from Fig. 5 that there is no obvious difference in the conformation of the two ligands, and the similarity value is 0.788, the rotation trend of the two is basically similar, and theRMSDis 1.484, thus proving the docking method is reasonable[20].

Fig. 5. Superimposed conformation of the docking ligand and the original ligand

Molecular docking analysis was performed in 20 newly designed molecules. The scoring functions including total-score, crash and polar were used as the evaluation criteria. The data obtained are shown in Table 2.

Table 2. Molecular Docking Scoring Function

Total-score[21]is the total scoring function, which ranks the binding affinity of the ligand and the active part of the receptor, and outputs the total score value. The higher the total-score, the higher binding affinity with the target protein for the newly ligand[22]. Crash reveals the inappropriate degree of ligand docked into the receptor.The smaller the absolute value, the more harmonious the docking position of the ligand. Polar is the score of the polarity function. When the binding pocket is on the surface of the molecule, the higher the score is, the better the docking result will be. When the binding site is located inside the molecule, the lower score, the better docking result[23].

Fig. 6. Molecular hydrogen bonding interaction diagram

The docking analysis results are shown in Fig. 6, in which the gray stick represents small molecular ligands,and the green one shows the amino acid residues that form hydrogen bonds. The hydrogen bonds are represented by yellow dashed lines. The hydrogen bond force is the main force for maintaining the protein and compound ligand molecules, which can make the combination more stable between the ligands and amino acid residues.Fig. 6a is the hydrogen bond analysis map of the original ligand extracted from the 4QTD protein crystal after re-docking. It can be seen from Fig. 6a that original ligand was observed to form hydrogen bond with A/MET111 residues in protein crystals, the total-score,crash and polar are 8.0557, -2.9164 and 1.1214, respectively. From these data, it can be seen that docking tool successfully reproduces the natural binding process of the ligand.

Fig. 6b is the docking results of compounds 1～10 with protein 4QTD. It shows that the newly designed compound could form 5 hydrogen bonds with A/GLN37, A/GLY38,A/LEU168 and A/ASN114 residues at the active site of protein. The total-score, crash and polar are 9.1092,-2.9768 and 2.9748, respectively.

Fig. 6c is the docking result of compounds 1～18 at the active site of 4QTD protein. The figure shows that the newly designed compound ligand forms 6 hydrogen bonds with A/GLY38, A/GLN37, A/ASN114 and A/LEU168 residues in the active site of 4QTD protein. The total-score,crash and polar are 8.9930, -4.6780 and 2.8881, respectively.

Fig. 6d is the docking result of compounds 1～4 at the active site of 4QTD protein. The newly designed compound ligand forms 5 hydrogen bonds with A/ARG69,A/GLY38, and A/ASN114 residues at the active site of 4QTD protein. The total-score, crash and polar are 8.6159,-2.4155 and 1.9897, respectively.The analysis results show that the newly designed compound forms strong hydrogen bonds with amino acid residues A/GLN37, A/GLY38 and A/ARG69. These interactions enhance the binding strength between the ligand and the receptor. Therefore, the docking result of the newly designed compound with the protein crystal is reliable and beneficial.

4 CONCLUSION

In the insulin resistance of type 2 diabetes caused by obesity, JNK1 is a drug target for the treatment of type 2 diabetes, and it plays a key mediator role in metabolic disorders. In this article, HQSAR and Topmer CoMFA technologies are used to analyze the relationship between isoquinolone derivative and their anti-JNK1 activity.According to the statistical verification results, a robust predictive ability was finally obtained. Based on the obtained Topomer CoMFA and HQSAR results, 20 new isoquinolone derivative molecules of JNK1 inhibitors were designed through virtual screening. External verification of these new molecules shows that the activity values are higher than that of the original template molecules. Later,the binding mode and target between the ligand and protein receptor were studied by molecular docking technology, and the amino acid residues that formed the hydrogen bond interaction between the ligand and the crystal protein structure were found. It is proved that the newly designed compound has formed a strong hydrogen bond with amino acid residues A/GLN37, A/GLY38,A/ARG69, etc. These interactions enhance the binding strength of the ligand and the receptor. Therefore, the docking result[24]of the newly designed compound with the protein crystal is reliable and beneficial. This study provides a certain idea for the design and development of new JNK1 inhibitor as anti-type 2 diabetes drugs, which is helpful to better understand its inhibitory mechanism and provide a certain theoretical basis for future experimental verification of new compounds.

ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Funds of China (21475081), Innovation Supporting Plan of Shaanxi Province—Innovation Research Team (No. 2018TD-015), and the Graduate Innovation Fund of Shaanxi University of Science and Technology.