Computational study on new natural compound inhibitors of Traf2 and Nck-interacting kinase (TNIK)

Virtual screening of potential inhibitors of TNIK

The ligand-binding pocket of TNIK is an essential regulatory site for its activity. NCB-0846 can bind to this pocket region to inhibit the function of TNIK in normal environment. Therefore, we choose this pocket as the reference area. ZINC15 database is a free commercial database provided by Irwin and Shoichet Laboratories in the Department of Pharmaceutical Chemistry at the UCSF. We obtained a total of 14962 natural, named and purchasable small molecules from the ZINC15 database. NCB-0846 was used as a reference compound to evaluate the binding capacity of the candidate small molecules. After the LibDock module calculation of DS4.5, it was found that 2,403 compounds could stably bind to TNIK. Among them, 1246 compounds scored higher than NCB-0846 (LibDock score: 111.272, ranking: 1247). The top 80 compounds are listed in Table 1.

ADME and toxicity prediction

The pharmacological properties of NCB-0846 and 80 candidate compounds were predicted through the ADME module of DS4.5, including aqueous solubility, BBB penetration, CYP2D6 binding, hepatotoxicity, human intestinal absorption, PPB (Table 2). The results showed that all compounds except ZINC000100277550 and ZINC000028538573 were soluble in water (in water at 25°C). Among them, 35 compounds had better water solubility than NCB-0846. In terms of blood-brain barrier permeability, except ZINC000002528486 and NCB-0846 showed Medium permeability, the other compounds were Undefined. CYP2D6 plays an important role in drug metabolism [21]. Except for 13 compounds, all the other compounds including NCB-0846 were predicted to be non-inhibitors of CYP2D6. In addition, we found that NCB-0846 has hepatotoxicity. Among all the candidate compounds, only 34 showed no hepatotoxicity. For human intestinal absorption, NCB-0846 and 5 compounds showed suitable absorption level, while 64 compounds showed poor absorption. Finally, we found that 20 compounds had strong binding ability to plasma protein, while others had weak binding ability.

Next, we conducted a comprehensive test and evaluation of the safety of NCB-0846 and the compounds. The rodent carcinogenicity (based on the U.S. National Toxicology Program (NTP) dataset), AMES mutagenicity and developmental toxicity potential (DTP) properties of the candidate compounds were comprehensively texted and evaluated using the TOPKAT module of DS4.5 (Table 3). The results showed that 7 compounds had AMES mutagenicity and 9 compounds had DTP properties. Based on all the above results, ZINC000040976869 and ZINC000008214460 were identified as ideal potential TNIK inhibitors. They have no inhibitory effect on CYP2D6 activity, no hepatotoxicity, and in addition have strong water solubility and plasma protein binding capacity. What's more, they are safe, with almost no rodent carcinogenicity, AMES mutagenicity and DTP. Therefore, we selected them for follow-up studies (Figure 1).

Analysis of ligand binding and pharmacophore

In order to further analyze the binding mechanism of ligands and TNIK, we used the CDOCKER module of DS4.5 for calculation [22]. The RMSD between the crystal structure and the docking posture was 0.6Å, which indicates that the use of CDOCKER module was reliable. After applying the force field of CHARMm36, we connected the two compounds and NCB-0846 into the 3D structure of TNIK to calculate the interaction energy of CDO CKER. Results as shown in Table 4, the interaction energy of ZINC000040976869 was −57.066 Kcal/mol and the interaction energy of ZINC000008214460 was −58.181 Kcal/mol, which was significantly lower than that of NCB-0846 (−52.062 Kcal/mol). This indicated that both compounds may bind to TNIK more easily and stably compared to NCB-0846. Then, we performed structural analysis of the interactions between TNIK and ligands, including hydrogen bonds and other types of hydrophobic bonds (Figures 2, 3, Tables 5, 6). There were two hydrogen bonds formed by ZINC000040976869 and TNIK (A:CYS108:HN-ZINC000040976869:O24, ZINC000040976869:H80-A:LEU169:O). In addition, five Alkyl interactions were formed in this complex. Similarly, four hydrogen bonds were formed by ZINC000008214460 and TNIK (ZINC000008214460:H76-A:ASN158:OD1, A:GLY34: HA1-ZINC000008214460:O30, ZINC000008214460: H67-A:GLN157:O, ZINC000008214460:H75-A:ASP171:OD2) and four Alkyl interactions were formed. In addition, the calculation results of pharmacodynamic groups of the two compounds were shown in Figure 4. ZINC000040976869 had 29 feature pharmacophores, including hydrogen donors and hydrophobic centers. ZINC000008214460 had 31 feature pharmacophores, including hydrogen bond acceptors, hydrogen donors and hydrophobic centers.

Molecular dynamics simulation

The complex initial conformation of TNIK-compounds was obtained by molecular docking module of CDOCKER. Then, the stability of the complex was simulated by molecular dynamics using the Standard Dynamics Cascade module of DS4.5 under dynamic conditions. Under the action of CHARMm force field, the motion of molecules was displayed dynamically. The potential energy and RMSD results obtained after 1ns simulation are shown in Figure 5. We can see from the figure that the RMSD trajectorie of ZINC000040976869 reaches equilibrium after 300ps, while ZINC000008214460 reaches equilibrium after 800ps. And over time, their potential energy and RMSD tend to be stable. In addition, in order to analyze the volatility of various amino acids in the complex during molecular dynamics simulation, we calculated the RMSF values of all amino acids during simulation. It can be seen from Figure 6 that the TNIK-ZINC000040976869 complex fluctuated greatly around the amino acids Glu12, Gly96 and Arg180, while the TNIK-ZINC000008214460 complex fluctuated greatly around the amino acids Ile13, Gly96, Asn186 and Pro206.

Docking software and ligand library

Discovery Studio 4.5 (DS4.5, Accelrys, Inc.) is a software for molecular modeling and environmental simulation in the field of life sciences. It has various functions including characterization of protein, three-dimensional molecular construction, drug design, molecular docking, database screening and so on. Through this method, many lead compounds and drug candidates have been identified and improved. Our candidate small molecules were obtained from natural products (NP) database in ZINC15 database. ZINC15 is a free database of commercial compounds provided by Irwin and Shoichet laboratories in the Department of Pharmaceutical Chemistry of UCSF (University of California, San Francisco, USA). First of all, we used the LibDock module of DS4.5 to screen the candidate compounds. Furthermore, ADME (absorption, distribution, metabolism, excretion) and TOPKAT (toxicity prediction by computer assisted technology) modules of the software were used to analyze the pharmacological characteristics of compounds. Then, the CDOCKER module was used to dock the compounds with TNIK more accurately. Finally, in order to analyze the stability of TNIK-ligand complexes, molecular dynamics simulation was carried out.

Virtual screening using LibDock

Firstly, we obtained the crystal structure of TNIK in complex with NCB-0846 from the protein database (PDB) (PBD ID: 5d7a, 2.9 Å). LibDock is a rigid docking program in DS4.5. It used a grid located at the binding position and uses polar and apolar probes to calculate the hotspots of the complex. Then hot spots were used to coordinate ligands for favorable binding. The 3D chemical structure of TNIK was shown in Figure 7. The structure was prepared by removing crystal water, ligands and other heteroatoms, then adding hydrogen, protonation, ionization and energy minimization. Furthermore, Smart Minimiser algorithm and CHARMm force field (Cambridge, MA, USA) were used to minimize energy [30]. In order to screen TNIK inhibitors more accurately, the binding pocket region of TNIK and NCB-0846 was selected as the docking site. 9.3 Å was set as the spherical docking site diameter based on the PDB site records. Next, all prepared small molecules were butted to the defined active sites for virtual screening using LibDock. Finally, the compounds with different docked poses were ranked according to LibDock score.

ADME and toxicity prediction

The ADME (Absorption, Distribution, Metabolism and Excretion) module of DS4.5 was used to evaluate the water solubility, blood-brain barrier (BBB) permeability, cytochrome P450 2D6 (CYP2D6) inhibition, hepatotoxicity, human-intestinal absorption level and plasma protein binding (PPB) level of the selected compounds. Additionally, the TOPKAT (Toxicity Prediction by Komputer Assisted Technology) module was employed to assess and analyze the toxicity of all potential compounds. This module can accurately analyze and verify the toxicity and environmental effects of compounds based on their 2D molecular structures. The molecular toxicities we have mainly analyzed including rodent carcinogenicity, AMES mutagenicity, developmental toxicity potential (DTP). We comprehensively analyzed the pharmacological properties predicted by ADME and the results of toxicity analysis to select potential candidate TNIK inhibitors.

More precise molecular docking and pharmacophore prediction

The CDOCKER module of DS4.5 was used to conduct molecular docking research based on CHARMm36 force field. During the docking process, the receptor is held rigid, while the ligands are allowed to flex, so that higher precision docking results can be obtained. Initially, we obtained the crystal structure of TNIK from PBD. Then, DS4.5 was employed to prepare and process protein crystals. Considering that fixed water molecules may affect the formation of receptor-ligand complexes during rigid and semi-flexible docking, we removed crystalline water molecules from protein before this process [31, 32]. After that, hydrogen atoms were added to protein. In order to make the docking results more reliable, we removed NCB-0846 from the protein structure, and then re-docked into the crystal structure of TNIK. The binding site of TNIK was defined as a binding sphere with a diameter of 9.3Å centered on the NCB-0846 binding region. CHARMm36 force field is applied to the docking process of receptor and ligands. During the docking process, the ligand gradually recognized and bound with the residues in the receptor binding sphere. The CHARMm energy (interaction energy plus ligand strain) based on each complex posture and the interaction energy representing ligand binding affinity were calculated after docking process. The ligand poses with the lowest interaction energy and appropriate docking direction was selected for follow-up research. In addition, the 3D-QSAR pharmacophore generation module was used to display the pharmacophore of the compound. Up to 255 conformations were generated per molecule to represent a small molecule. However, only conformations with the energy below 10 kcal/mol were preserved.

Molecular dynamics simulation

To further investigate the binding process between TNIK and the two candidate compounds, molecular dynamics simulations were performed on the two optimal ligand-receptor complexes. We employed the Solvation module of DS4.5 to put the ligand-receptor complex into an orthogonal box and solvated it with the explicit periodic boundary solvent-water model. In order to simulate the action in physiological environment, chlorides with ionic strength of 0.145 were added to the system. After that, we applied CHARMm force field to the system and minimized the energy (500 steps for steep descent and 500 steps for conjugate gradient). The temperature of the system was slowly driven from 50K to 300K, and the equilibrium simulation time was 20 ps. Molecular dynamics simulations (production module) were run for 1 ns with time step of 2 ps. This process was carried out with NTP (Atmospheric Pressure and Temperature) system at a constant temperature of 300K. The particle mesh Ewald (PME) algorithm was used for long-range electrostatic calculations, and the Linear Constraint Solver (LINCS) algorithm was used to identify all hydrogen-containing bonds. In the next process, we used DS4.5 to analyze molecular dynamics trajectories, including Root-Mean-Square Deviation (RMSD), Root-Mean-Square Fluctuation (RMSF), potential energy and structural characteristics.