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

Traf2 and Nck-interacting kinase (TNIK) is the downstream molecule of Wnt/β-catenin signal pathway. As the activation kinase of β-catenin/T-cell factor 4 transcription complex, it can fully activate Wnt signalling and promote the growth and invasion of tumor cells. We conducted computer-assisted virtual screening and a series of analyses to find potential inhibitors of TNIK. First, LibDock was used for molecular docking of natural small molecules. Then, ADME (Adsorption, Distribution, Metabolism and Excretion) analysis and toxicity prediction were performed on the top 80 small molecules which have higher scores. Additionally, in order to further determine the affinity and binding mechanism of TNIK-ligands, we analyzed the pharmacophores and used CDOCKER for more accurate molecular docking. Last but not least, molecular, dynamics simulation was used to evaluate the stability of receptor-ligand complexes in natural environment. The results showed that natural small molecules (ZINC000040976869 and ZINC000008214460) had high affinity and low interaction energy with TNIK. They were predicted to have excellent pharmacological properties, such as high plasma protein binding capacity and water solubility, no hepatotoxicity, no blood-brain barrier permeability and tolerant with cytochrome P450 2D6 (CYP2D6). In addition, they have less rodent carcinogenicity, AMES mutagenicity, and developmental toxicity potential. Molecular dynamics simulations showed that the two compounds could achieve the stability of potential energy and Root-Mean-Square Deviation (RMSD) at different time nodes. This study proves that ZINC000040976869 and ZINC000008214460 are ideal lead compounds with inhibition targeting to TNIK. These compounds provide valuable ideas and information for the development of new colorectal cancer targeting drugs.


INTRODUCTION
Colorectal cancer (CRC) is a worldwide disease that seriously harms human health. According to the latest global cancer statistics in 2020, CRC is the third most common cancer in the world, and it has caused more than 900,000 deaths. Its mortality rate has reached the second highest among cancers [1]. Not only that, the incidence of CRC continues to rise along with people's tendency towards red meat diet, reduced physical activity and increased body mass index [2][3][4]. Surgery, as the main treatment, has developed rapidly in recent years, but it is only effective for CRC patients without lymph nodes and distant metastasis [5]. Therapeutic antibodies against epidermal growth factor receptor (EGFR) and vascular endothelial growth factor (VEGF) combined with chemotherapy, as an important treatment for metastatic or recurrent CRC, have been widely used in clinic in recent years, but they have not significantly increased the 5-year survival rate of patients with advanced CRC [6]. Therefore, it is a very hot topic to find related therapeutic targets and to study new targeted drugs for CRC.
The alteration of Wnt signaling pathway plays a crucial role in the occurrence and development of CRC by affecting the stemness, metastasis and tissue repair of tumor cells [7,8]. About 80% of colorectal cancers are thought to be related to the activation of components in Wnt signaling pathway caused by mutation of adenomatous polyposis coli (APC) gene [9]. The mutation of APC gene may account for the degradation failure of β-catenin located downstream of Wnt signaling pathway [10]. The β-catenin protein is transported into the nucleus to further activate Wnt signaling. The accumulation of β-catenin in the nucleus activates T-cell factor 4 (TCF4) and forms a complex with it. TCF4 is a member of TCF/LEF transcription factor family, which is an essential factor leading to CRC [9,11]. Due to the genetic inactivation of APC, only the molecules located downstream of APC are considered as effective targets in Wnt signaling pathway.
Traf2 and NCK-interacting kinase (TNIK) is a member of STE20 serine/threonine protein kinase family, which has an N-terminal activation domain and can specifically activate the c-Jun N-terminal kinase pathways like many germinal center kinases [9,12]. TNIK is the downstream signal protein and nuclear coactivator of Wnt signaling pathway, which is mainly located in the nucleus. In addition, it is a component of the TCF4 and β-catenin transcriptional complex in CRC cells [13,14]. When combined with TNIK, the conserved serine 154 residue on TCF4 protein can be phosphorylated by it [14,15]. This effect is indispensable for the complete activation of Wnt signaling pathway and the growth of CRC cells. Therefore, TNIK is a feasible drug therapy target for rectal cancer caused by abnormal Wnt signaling pathway. NCB-0846 is a TNIK inhibitor, which has been widely studied. It can combine with the ATP activity binding pocket of TNIK and completely inhibit Wnt signal transduction. A study showed that oral administration of NCB-0846 significantly inhibited the growth of tumor transplanted into immunocompromised mice [7]. However, due to its poor pharmacokinetics and low activity, it has not been widely used in clinic [16]. Based on the above reasons, we need to explore safer and more effective TNIK inhibitors for the treatment of CRC.
Natural products and their derivatives have brought new horizons to drug research and development in recent years because of their potential biological functions and unique molecular structure, and have become an important source for the pharmaceutical industry to develop new drugs [17,18]. This study intends to identify potential novel TNIK inhibitors based on the ZINC database using molecular docking, pharmacological analysis, and molecular dynamics simulations [19,20].

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 bloodbrain 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    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.   Table 6.

DISCUSSION
As a major cancer, the incidence of CRC has been increasing in recent years, accounting for 10% of cancer-related mortality [1]. This undoubtedly places a heavy health and economic burden on the world every year. Although various treatment methods have made great progress, the cure rate of CRC is still not optimistic. At present, the widely used CRC targeting drugs are mainly targeted at VEGF (bevacizumab) and EGFR (cetuximab and panitumab). But this survival benefit is also limited to a few months [23]. Therefore,   it is of great significance to find new target drugs in the molecular signaling pathway of CRC. It has been demonstrated that 90% of CRC patients have typical Wnt/β-catenin signaling pathway gene mutations [9]. Activation of the Wnt signaling pathway leads to the generation of cancer stem cells [9]. The mutation of APC tumor suppressor gene is more than 80% in CRC and is the earliest step in the carcinogenesis of CRC [24]. Above all, it is perspective to search for molecules downstream of APC in the Wnt/β-catenin signaling pathway as targets for drug therapy.
As a member of serine/threonine protein kinase family, TNIK has been proved to play a cancer-promoting role in CRC, gastric cancer, lung cancer, prostate cancer and other diseases [25][26][27]. As the most downstream molecule of the Wnt/β-catenin signaling pathway, its feasibility as a drug target has been confirmed by many researches [7,14]. In CRC, inactivated APC results in β-catenin accumulation that cannot be degraded. The excess β-catenin will be transported into the nucleus and form a transcription complex with TCF4. As an activating kinase of the β-catenin/TCF4 transcriptional complex, TNIK phosphorylates the TCF4 protein at the conserved serine 154 residue, thereby activating the transcriptional complex. This will fully activate Wnt signaling pathway, thus promoting the growth and invasion of tumor. The study of Gui et al. demonstrated that TNIK gene knockout can block the activation of Wnt signaling and inhibit the growth of tumor cells [28]. At present, several TNIK inhibitors have been developed, and the most representative of which is NCB-0846. However, due to the poor pharmacological characteristics and low activity, it has not been used in clinic so far. As a practical TNIK inhibitor, NCB-0846 can bind to the ATP active pocket of TNIK, making it unable to phosphorylate TCF4, thereby exerting a tumor suppressive effect [7,29]. Consequently, we used NCB-0846 as a reference drug to screen for more ideal TNIK inhibitors.
To explore effective TNIK inhibitors, 14962 natural, named and purchasable small molecules were obtained from the ZINC15 database for screening. Based on the binding stability, pharmacological properties, toxicity and stability of compounds, the ideal TNIK potential inhibitors were determined. LibDock is a preliminary screening method for small molecules, the higher the score, the more stable the receptor-ligand conformation, and the more optimized the energy. According to the results of LibDock, we found that 1246 compounds scored higher than the reference drug NCB-0846. Thus, these compounds may form more stable complexes with TNIK compared to NCB-0846. The top 80 natural compounds were selected based on LibDock score for further analysis.
Next, we analyzed the pharmacological and toxicological properties of the candidate compounds and NCB-0846 using ADME and TOPKAT modules. Among them, ZINC000040976869 and ZINC000008214460 attracted our attention. First of all, these two compounds have strong binding ability to plasma proteins, good water solubility, can be absorbed by the intestinal tract and can not pass through the blood-brain barrier. In addition, they have no inhibitory effect on CYP2D6 and no hepatotoxicity. On the other hand, the rodent carcinogenicity, AMES mutagenicity and developmental toxicity potential of these two compounds are relatively low. Compared with NCB-0846, these two compounds have less hepatotoxicity, better water solubility, do not cross the blood-brain barrier and have lower toxicity. Based on the above results, it shows that ZINC000040976869 and ZINC000008214460 have good prospects in drug development and can be used as ideal lead compounds for further analysis. Although other candidate compounds have not been selected, they can reduce toxicity and improve pharmacological properties by adding or modifying the groups they contain. For this reason, they can still be listed as candidate drugs.
CDOCKER is a technology that produces highprecision molecular docking results in the CHARMm position, which is more precise than LibDock. CDOCKER interaction energy is an index to evaluate the affinity of receptor-ligand. The lower it is, the higher the affinity between ligand and receptor. Our results showed that the CDOCKER interaction energy of ZINC000040976869 and ZINC000008214460 were −57.066 Kcal/mol and −58.181 Kcal/mol respectively, which was lower than that of NCB-0846 (−52.062 Kcal/mol). This indicates that the affinity and stability of the two small molecules we selected were higher than those of NCB-0846. In addition, the protein binding sites of NCB-0846 and these two compounds are the same, and they all have axisymmetric structures. To sum up, compared with NCB-0846, they may have better inhibitory effect and higher safety.
In the calculation of the pharmacophore of the compound, we found that ZINC000040976869 displayed several hydrogen donors and hydrophobic centers, with a total of 29 pharmacophores. Similarly, ZINC000008214460 displayed several hydrogen bond acceptors, hydrogen donors and hydrophobic centers, with a total of 31 pharmacodynamic groups. This means that in future research, the two compounds have the potential to add different functional groups to improve drugs and improve anti-cancer efficacy.
Next, we used molecular dynamics simulations to analyze the two optimal conformations of TNIK-compounds complexes. This can scientifically evaluate its stability in the natural environment. We found that RMSD trajectories and interaction energy of ZINC000040976869 and ZINC000008214460 tended to be stable after 300ps and 800ps, respectively. This indicated that the hydrogen bonds and Alkyl interactions between the two compounds and TNIK played a stable role in the formation of the complexes. Complexes can exist stably in natural environment. In addition, by analyzing the RMSF trajectories of the two compounds, we can see that TNIK-ZINC000040976869 complex fluctuates greatly around amino acids Glu12, Gly96 and Arg180, while TNIK-ZINC000008214460 complex fluctuates greatly around amino acids Ile13, Gly96, Asn186 and Pro206. The fluctuation trends of RMSF trajectories of the two complexes are different, which means that the fluctuations of amino acids of different complexes are different in the simulation process. Therefore, the amino acid positions with large fluctuations in RMSF of these two compounds can be optimized in the subsequent drug modification process.
The most important step in drug design is to identify reasonable lead compounds, followed by continuous modification and improvement to make them clinically applicable. Based on the above studies, the two natural compounds we have identified can provide great help for the development of targeted drugs for CRC. Although our study was rigorously designed, there are some limitations. For example, we may need to conduct animal experiments, molecular biology experiments, etc. to further confirm our results. In the future research, we can also evaluate LD50, ED50, Maximum Tolerated Dosage (MTD) and Aerobic Biodegradability (A.B.) and other indicators.

CONCLUSION
In this study, a series of biological and chemical methods (Virtual Screening, Molecule Docking, ADME, Toxicity Prediction, Molecular Dynamics Simulation) were used under computer assistance to screen and analyze the potential TNIK inhibitors. After comprehensive analysis, two natural fractions, ZINC000040976869 and ZINC000008214460, were identified as ideal drug candidates. They have strong affinity with TNIK and good pharmacological properties and safety, which makes them promising to be further developed into new CRC targeted drugs.

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, threedimensional 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.

AUTHOR CONTRIBUTIONS
This study was completed with a teamwork. Every author has made substantial contributions to the study. Lushun Ma has come up with the conception and was responsible for the creation of new software used in the work. Additionally, Rui Li, Chunxiang Liu and Shuxian Chen has done the design of the work and drafted the work. Furthermore, analysis of the data was done by Zhiwei Yao and Yong Liu. As for interpretation of the data, Bo Wang and Heng Wang have contributed a lot to this part. And Daqing Sun has substantively revised it. In addition, Lushun Ma and Zhiwei Yao responded and revised the manuscript according to the revised opinions.