Identification of novel inhibitors for TNFα, TNFR1 and TNFα-TNFR1 complex using pharmacophore-based approaches

Background Tumor necrosis factor α (TNFα) is a multifunctional cytokine with a potent pro-inflammatory effect. It is a validated therapeutic target molecule for several disorders related to autoimmunity and inflammation. TNFα–TNF receptor-1 (TNFR1) signaling contributes to the pathological processes of these disorders. The current study is focused on finding novel small molecules that can directly bind to TNFα and/or TNFR1, preventing the interaction between TNFα or TNFR1, and regulating downstream signaling pathways. Methods Cheminformatics pipeline (pharmacophore modeling, virtual screening, molecular docking and in silico ADMET analysis) was used to screen for novel TNFα and TNFR1 inhibitors in the Zinc database. The pharmacophore-based models were generated to screen for the best drug like compounds in the Zinc database. Results The 39, 37 and 45 best hit molecules were mapped with the core pharmacophore features of TNFα, TNFR1, and the TNFα–TNFR1 complex respectively. They were further evaluated by molecular docking, protein–ligand interactions and in silico ADMET studies. The molecular docking analysis revealed the binding energies of TNFα, TNFR1 and the TNFα–TNFR1 complex, the basis of which was used to select the top five best binding energy compounds. Furthermore, in silico ADMET studies clearly revealed that all 15 compounds (ZINC09609430, ZINC49467549, ZINC13113075, ZINC39907639, ZINC25251930, ZINC02968981, ZINC09544246, ZINC58047088, ZINC72021182, ZINC08704414, ZINC05462670, ZINC35681945, ZINC23553920, ZINC05328058, and ZINC17206695) satisfied the Lipinski rule of five and had no toxicity. Conclusions The new selective TNFα, TNFR1 and TNFα–TNFR1 complex inhibitors can serve as anti-inflammatory agents and are promising candidates for further research. Electronic supplementary material The online version of this article (10.1186/s12967-019-1965-5) contains supplementary material, which is available to authorized users.


Background
Tumor necrosis factor α (TNFα) is a cytokine secreted by macrophages in response to septic shock, inflammatory agents and cachexia. TNFα plays a key role in the immune system and cell death (e.g., apoptosis and necrosis) [1]. TNFα is involved in a number of autoimmune diseases, including psoriasis, inflammatory bowel disease, rheumatoid arthritis, systemic sclerosis, systemic lupus erythematosus, multiple sclerosis, diabetes and ankylosing spondylitis [2]. Since TNFα is an important mediator in infections and tumors, a series of biological agents targeted to TNFα has been developed for the treatment of cancer and autoimmunity [3].
TNFα contributes to the pathogenesis of inflammatory, edematous, neovascular, and neurodegenerative diseases of the eye [4]. Injection of TNFα into animal eyes induces breakdown of the blood-retina barrier [5]. Furthermore, increased levels of TNFα and TNF-receptors (TNFRs) have been found in the serum of humans with uveitis. Upregulation of TNFα expression has been shown in keratocytes of patients with rheumatoid corneal ulcerations [6]. Moreover, there is increasing evidence of TNFα Saddala and Huang J Transl Med (2019) 17:215 involvement in the pathogenesis of experimental retinal neovascularization, proliferative vitreoretinopathy, and macular edema [7]. In an in vivo animal model of retinal injury, Berger et al. [8] showed that TNFα played a deleterious role in ischemia-reperfusion injury. Direct neutralization of this cytokine partially preserved retinal function [8]. The diverse characteristics of TNFα were attributed in part to the timing of its expression after injury. Nagineni et al. [9] demonstrated that inflammatory cytokines, including interleukin 1 beta (IL-1β), interferon gamma (IFN-γ) and TNFα, increased the secretion of vascular endothelial growth factor (VEGF)-A and -C by human retinal pigment epithelial (RPE) cells and choroidal fibroblasts, with VEGF being the most important factor for initiating pathological ocular neovascularization [9]. TNFα is crucial for the pathogenesis of diabetic retinopathy in rodents, and its pharmacological blockade leads to the inhibition of retinal cell death [10,11].
A variety of TNFα antagonists, including infliximab, etanercept, adalimumab, certolizumab and glolimumab, were developed for therapeutic applications [12]. However, these biological therapies exhibited inevitable weaknesses, such as risk of infection, high cost, and the requirement for intravenous injections. By contrast, small molecule inhibitors are relatively cheaper and can be taken orally. Therefore, the identification of small molecules that can inhibit TNFα-regulated pathways is a promising research area that has lately received much attention.
Therefore, in the present study we used cheminformatics as part of the pipeline [pharmacophore modeling, virtual screening, molecular docking and in silico ADMET (absorption, distribution, metabolism, excretion and toxicity) analysis] to screen for novel, safe TNFα and TNFR1 inhibitors from the publicly available Zinc database.

Preparation of target proteins
We took the crystal structures of the target proteins: TNFα (2AZ5) with resolution 2.1 Å [13] and TNFR1 (1EXT) with resolution 1.85 Å, from the PDB (https :// www.rcsb.org/). The TNFα-TNFR1 complex protein was downloaded from a public web site (http://www.cblig and. org/downl oads/TNF_TNFR1 .pdb). We removed all the hetero atoms and crystal water molecules from the target proteins and minimized the energy.

Active site prediction
The active sites of TNFα, TNFR1 and the TNFα-TNFR1 complex were predicted using the CASTp (Computed Atlas of Surface Topography of proteins) server (http:// sts.bioe.uic.edu/castp /index .html?2pk9). CASTp measures and identifies pockets and pocket mouth openings, in addition to the cavities. We uploaded the target proteins as input to predict the ligand binding sites. The CASTp server predicted the key amino acids for binding interactions to the inhibitors [14].

Target-ligand pharmacophore generation
The structure-based pharmacophore technique can be used to advance the drug development process. For pharmacophore modeling, we selected three PDB (protein data bank) structures, i.e. TNFα (2AZ5), TNFR1 (1EXT), and the TNFα-TNFR1 complex (modelled protein) and their inhibitors (default inhibitor: 307), physcion-8-glucoside (ZINC33832439), Erythrosine B (ZINC08214556). The structure of the TNFα-307 complex was used as the fundamental of the pharmacophore modelling while TNFR1 with the physcion-8-glucoside and TNFα-TNFR1 complex with Erythrosine B complexes were used for pharmacophore design. ZINCPharmer (http://zincp harme r.csb.pitt.edu) is an online interface for searching the purchasable compounds of the Zinc database using the Pharmer pharmacophore search technology. ZINCPharmer can automatically extract a set of pharmacophore features from the molecular structure. Each feature comprises the feature type (hydrophobic, hydrogen bond donor/acceptor, positive/negative ion or aromatic), a position, and a search radius [15]. We provided both a receptor and bound-ligand structures, and ZINCPharmer is automatically identified an interaction pharmacophore. All possible pharmacophore features on the ligand were computed; however, only those that were within a distance cutoff of complimentary features on the receptor are enabled. We modified the pharmacophore features for TNFα, TNFR1 and the TNF-TNFR1 complexes and set parameters such as the hydrogen bond acceptors (HBA)/donors (HBD) are within 4 Å, charged and aromatic features are within 5 Å of the receptor.

Pharmacophore based virtual screening
The modelled pharmacophores features were used as query features for screening for small molecules against the Zinc purchasable compound database. The TNFα has three aromatic spheres, two hydrophobic spheres, and one HBA spheres; TNFR1 has two aromatic spheres, one hydrophobic sphere, one HBD and one HBA spheres; TNF-TNFR1 has two aromatic spheres, one hydrophobic sphere, one HBD and one HBA spheres. Each pharmacophore model feature searched on Zinc purchasable compounds to get the hits based on matched features. The TNFα pharmacophore model obtained 39 hits, while TNFR1 had 37 and the TNF-TNFR1 complex had 45 hits. These hits were used for the docking studies.

Docking simulation
Molecular docking studies were carried out with the AutoDock 4.2 in PyRx Virtual Screening Tool, which was used to generate the docking key files. Experiential free energy utility and Lamarckian genetic algorithm (LGA) with the following settings: a maximum of 2,500,000 energy evaluations, a preliminary population of 150 randomly placed individuals, a maximum of 27,000 generations, a transmutation rate of 0.02, and a crossover velocity of 0.8, along with an exclusiveness rate (numeral of top individuals to endure to the next generation) of one were designed for docking. The supposed Solis and Wets law was useful to a maximum of 300 iterations for each look for confined search. Default principles were thought to be designed for all the parameters not previously mentioned.

Chemical analysis of drug-likeness
The drug-likeness properties were analyzed using Mol-Soft Drug-Likeness explorer (http://www.molso ft.com/ mprop /) and the FAF-Drugs4 server (http://fafdr ugs4. mti.univ-paris -dider ot.fr/). Drug-likeness was indicated by the Lipinski "Rule of 5'' [16]. Drug likeness can be described as a complex balance of various molecular properties and structural features, that determine whether a molecule is a drug or a non-drug. In accordance with the Lipinski "Rule of 5'' a compound has a lot of possible elected membrane permeability and merely captivated through the body, if its relative molecular mass is a less than 500, its lipophilicity, expressed as an amount referred to as LogP is a less than five, the number of groups within the compound that may give hydrogen atoms to hydrogen bonds is a less than five, and the number of group that may settle for hydrogen atoms to make hydrogen bonds is a less than 10 [17].

Prediction of physicochemical descriptors and ADMET parameters
The physicochemical profiles of lead compounds can increase the quality of clinical candidates [18]. The individual consideration of ADMET behaviors in the early stages of drug discovery have decreased the fraction of global pharmacokinetics related to failures in later phases of development. ADMET parameters of the best 15 compounds were predicted by SwissADME tools [19]. Swis-sADME predicts BBB (blood brain barrier) penetration and GI (gastro intestine) absorption by BOILED-Egg method [20].

Cheminformatics pipeline
TNFα is a cell signaling protein (cytokine) involved in systemic inflammation and is one of the cytokines that comprise the acute phase reaction. TNFα is produced chiefly by activated macrophages, although it can be produced by many other cell types. It is associated with a variety of important physiological processes and pathological conditions [21]. To control the adverse effects of TNFα, the current efforts have focused on blocking TNFα binding to its receptor. The overview flow chart of the cheminformatics pipeline for the present study is shown in Fig. 1.

Target proteins preparation
The crystal structures of TNFα with a small molecule inhibitor (307) (2AZ5) [13] with resolution 2.1 Å, R-value free 0.278, and R-value work 0.220 and an extracellular domain of the 55 kda TNFR1 (1EXT) with resolution 1.85 Å, and R-value free 0.243, R-value work 0.203 were downloaded from the protein data bank (PDB) (https :// www.rcsb.org/) [22]. The TNF α-TNFR1 protein complex was downloaded from Chen et al. [23] public web site (http://www.cblig and.org/downl oads/TNF_TNFR1 .pdb). All hetero atoms and crystal water molecules were removed from target proteins. We performed energy minimization by using the AutoDock Vina tool [24] with the following energy minimization parameters: a UFF (Universal force field) force field, a steepest descent optimization algorithm, 2000 steps for run, 1 step for update and an energy difference of less than 0.1. After energy minimization, target proteins were used for further analysis.

Target-ligand pharmacophore generation
The molecular binding process relies on several properties and features of the amino acids presenting in the active site [25]. According to the IUPAC (International Union of Pure and Applied Chemistry) definition, a pharmacophore is the ensemble of steric and electronic features that are necessary to ensure the optimal supramolecular interactions with a specific biological target structure and to trigger (or to block) its biological response. A pharmacophore model comprises several features organized in a specific 3D pattern. Each feature is typically represented as a sphere. Such pharmacophore features are typically used as queries to screen small molecule libraries of compounds [26]. The TNFα pharmacophore contains six pharmacophoric features with a default ligand 307 (

Pharmacophore based virtual screening
The pharmacophore design models were used for screening molecules against the Zinc database (https ://zinc. docki ng.org/). The pharmacophore features acted as query parameters along with hit reduction, hit screening and subset selection, which we set. For examples, Max Hits per Conformation was set to 10, Max Hits per Molecular was set to 1, Max Total Hits was set to maximum, Max RMSD (root mean square deviation) was set to 2, molecular weight (MW) was set to 300-500, rotatable bonds was set to 0-10. All the hit reduction and hit screening parameters were applied against the purchasable Zinc subset. Accordingly, the success of a virtual screening performance can be quantified by the enrichment factor (EF) and hit rate (HR) when the percentage of active compounds in the screening database is known. The TNFα pharmacophore model obtained 39 hits, TNFR1 had 37 and TNF-TNFR1 complex had 45 hits. The pharmacophore hits were mapped with the pharmacophore models. All the hit molecules were submitted to molecular docking studies.

Target proteins and small compounds docking simulation
Protein-ligand docking is the generally used docking algorithm. It predicts the site of a ligand when it is bound to its protein [27]. Mostly docking algorithms can make many possible structures; thus, the means to score each structure is also required, to identify those of greatest interest. Docking was performed using the AutoDock in PyRx Virtual Screening tool [28,29]. The hit molecules were docked into the active site of TNFα, TNFR1 and the TNFα-TNFR1 complex. Based on the binding conformation AutoDock generated binding energies for all molecules.  (Fig. 4). Ile58, Leu120, Gly121, and Tyr151 are key residues for interacting with small molecules to inhibit TNFα trimer formation. The TNFR1 inhibitors, ZINC02968981 molecule ON, ON, ON, OC and OC groups interacted with five hydrogen bonds to the Lys75, Gln82, Gln82, Arg104 and Tyr106 amino acids CO, CO, CN, CN, and CO functional groups with − 10.1 kcal/mol −1 binding energy respectively. The ZINC02968981 molecule ON, ON, ON, OC and OC groups were interacted with five hydrogen bonds to the Lys75, Gln82, Gln82, Arg104 and Tyr106 amino acids CO, CO, CN, CN, and CO functional groups with − 10.1 kcal/mol −1 binding energy respectively. The ZINC09544246 molecule OC, NH, NH, OC, OC, NH and NH groups interacted with seven hydrogen bonds to the Glu56, Glu56, Ser57, Ser59, Cys70, Cys73 and Ser74 amino acid CN, CO, CO, CO, CN, CO and CO functional groups with − 9.8 kcal/mol −1 binding energy respectively. The ZINC58047088 molecule NH, and NH groups interacted with two hydrogen bonds to the Ser74 and Asn110 amino acid CO and CO functional groups with − 9.5 kcal/mol −1 binding energy respectively. The ZINC72021182 molecule OH group interacted with one hydrogen bond to the Arg104 amino acids CN functional group with − 9.3 kcal/mol −1 binding energy respectively. The ZINC08704414 molecule NH, NH, OC, OC, and NH groups interacted with four hydrogen bonds to the Ser74, Lys75, Agr77 and Asn110 amino acid CO, CO, CN, and CO functional groups with − 9.1 kcal/ mol −1 binding energy respectively. The ZINC09544246 (Query) molecule OC, OH, OH, OH, OH, OH, OH and OH groups interacted with eight hydrogen bonds to the Arg77, Cys96, Cys96, Arg104, Arg104, Arg104, Tyr106 and Asn110 amino acid CN, CN, CO, CN, CN, CN, CO and CO functional groups with − 7.6 kcal/mol −1 binding energy respectively (Fig. 5).
The TNFα-TNFR1 complex inhibitors, ZINC05462670 molecule OH, OH, OC, and OH groups interacted with four hydrogen bonds to the Ser74, Thr94, Glu109 and Asn110 amino acid CO, CO, CO, and CO functional groups with − 10.0 kcal/mol −1 binding energy respectively. The ZINC35681945 molecule NH, ON, ON, NH, OC and NH groups interacted with six hydrogen bonds to the Lys75, Gln82, Gln82, Asp93, Arg104 and Asn110 amino acid CO, CO, CN, CO, CN and CO functional groups with − 9.7 kcal/mol −1 binding energy. The ZINC23553920 molecule OC, NH, and NH groups interacted with three hydrogen bonds to the Ser74, Thr94, and Thr94 amino acid CO, CO, and CO functional groups with − 9.5 kcal/mol −1 binding energy respectively. The ZINC05328058 molecule OH, OC, OH, OH and OH groups interacted with five hydrogen bonds to the Ser74, Gln82, Pro90, Asn92, and Cys96 amino acid CO, CO, CO, CN and CN functional groups with − 8.9 kcal/mol −1 binding energy respectively. The ZINC17206695 molecule NH, NH, NH, ON, NO and OC groups interacted with six hydrogen bonds to the Ser74, Ser74, Asn93, Asn110, Asn110 and Ser147 amino acid CO, CO, CO, CN, CN and OC functional groups with − 8.5 kcal/mol −1 binding energy respectively. The ZINC08214556 (Query) molecule OC, and OC groups interacted with six hydrogen bonds to the Thr94 and Asn110 amino acid CO, and CN functional groups with − 7.2 kcal/mol −1 binding energy respectively (Fig. 6). The functional key residues play a vital role in TNFR1 and TNFα-TNFR1 complex   [1,2,4] triazolo [1,5-      formation and activation of the TNFα signaling pathway in proinflammation.

Prediction of physicochemical descriptors and ADMET parameters
We analyzed physicochemical descriptors and ADMET parameters by FAF-Drugs4 and SwissADME analysis to find the solubility and permeability of the 15 ligand molecules in order to use them for experimental assays and to reach their site of action in an accurate drug ability. The molecular complexity of the fifteen ligands could be measured by the number of rings and aromatic rings, the fraction of carbons that were sp3 hybridized (Fsp3), or the number of stereocenter properties and ADMET properties, which were all computed by FAF-Drugs4 (Additional file 1: Figs. S1 and S2). The TNFα, TNF1 and TNFα-TNF1 complex of best compounds (ZINC09609430, ZINC02968981 and ZINC05462670) FAF-Drugs4 results are presented in Fig. 7. The in silico ligand toxicity and biological property predictions are faster and more reliable approaches to take before further exploring experimental authentications such as in vitro and in vivo tests. All fifteen inhibitors were screened with the SwissADME server. The results revealed that all 15 ligands were safe and passed the lipophilicity, water solubility, pharmacokinetics, drug likeness and medicinal chemistry properties. All fifteen TNFα, TNF1 and TNFα-TNF1 complex inhibitor molecules obeyed the Lipinski rule of five and ADMET properties with biologically possible activity (Additional file 1: Tables S1, S2, S3). Therefore, these TNFα, TNF1 and TNFα-TNF1 complex inhibitors are most appropriate for additional drug discovery approaches to drug discovery.

Discussion
The present study screened for novel small inhibitors that can specifically inhibit TNFα-TNFR interaction and downstream signaling. We identified the key amino acid residues involved in the interactions of the TNFα and TNFR1 proteins. We screened the lead compound hits for TNFα, TNFR1 and the TNFα-TNFR1 complex from the Zinc library using structure-based pharmacophore modeling, virtual screening, and molecular docking along with in silico ADMET analysis. The identified novel small inhibitors can potentially be utilized for antiinflammatory agents to treat relevant disorders. The pharmacophore features are the key elements to screen for the best, potent small molecules binding to target proteins from publicly available databases. Pharmacophore-based approaches were widely used in virtual screening, de novo design and other applications such as lead optimization and multitarget drug design [30]. For TNFα, TNFR1 and the TNFα-TNFR1 complex, we used six pharmacophore features with the default ligand (307), five pharmacophore features with a Physcion-8-glucoside ligand, and five pharmacophore features with an Erythrosine B ligand respectively. Based on these pharmacophore features, we identified 39, 37 and 45 best hits from the Zinc database. Molecular docking results revealed that the aforementioned hits exactly docked into the active site of TNFα and TNFR1. Protein-ligand interactions suggested that the functional groups (residues) mimic the binding of hits and fit well into the active domain of TNFα, TNFR1 and the TNFα-TNFR1 complex. In particular, the Ile58, Leu120, Gly121, Tyr515, Glu56, Ser57, Ser59, Cys70, Cys73, Ser74, Lys75, Arg77, Gln82, Cys96, Arg104, Tyr106, and Asn110 residues are critical for the inhibitory interaction between TNFα, TNFR1 and the TNFα-TNFR1 complex. These key residues are located in the TNFR1 binding site of the TNFα protein. The in silico ADMET results revealed that all the top five of TNFα, TNFR1 and the TNFα-TNFR1 complex's inhibitors are virtually safe and active.
TNF is a cytokine protein expressed by activated monocytes/macrophages (including central nervous system [CNS] microglia), activated NK (Natural killer) and T cells, and by a diverse array of non-immune cells such as endothelial cells and fibroblasts [31]. TNFα is produced in two forms, soluble TNFα (sTNFα) and membranebound TNFα (tmTNFα). The soluble form of TNFα is created from the tmTNFα extracellular domain by the matrix metalloproteinase TNFα converting enzyme (TACE) [32]. Membrane-bound TNFα is able to serve as a ligand binding to TNFR or as a receptor mediating the transfer of external signals back to the cell which has exprimed it on its surface [33]. Both cytokine forms, i.e. soluble and membrane bound are active as homotrimers with a characteristic cone-shape. The five best TNFα inhibitors interact with the TNFα homodimer and inhibit the active form of homotrimers. Oanh et al. [34] reported that the triterpene saponins had a good binding affinity       with protein TNFα and were docked to the pore at the top of the bell or cone shaped TNFα trimer. Mehreen et al. [35] also reported that the novel small molecules interacted with TNFα trimer. The docking results are in agreement with the findings from the literatures. Our in silico method identified that TNFα inhibitors may disrupt the trimer formation of TNFα. The trimer form of TNF binds to TNFR1, activates the downstream signaling, and predominantly promotes inflammation and tissue degeneration [36]. The five best TNFR1 inhibitors interacted with the TNFα binding site of TNFR1 and inhibited the TNF/TNFR1 signaling pathway. Chen et al. [23] also reported that small molecules that directly bind to TNFα or TNFR1, inhibit the interaction between TNFα and TNFR1, and/or regulate related signaling pathways. Fischer et al. [37] reported the sTNF/ TNFR1 signaling as a new therapeutic target pathway. Recent researchers have focused primarily on identifying small molecules that directly bind to TNFα or TNFR1 [38], inhibit the binding of the TNFα and TNFR1 [39] and regulate related signal pathways [40]. The TNFR1 docking results are consistent with the results by other investigators. Identifying potential inhibitors of TNFα, TNFR1 and TNFα-TNFR1 complex and its analogues is thus an attractive strategy for treating inflammatory diseases, such as in central nervous system (i.e. brain and retina). Using our established cheminformatics pipeline, we identified 15 inhibitors. These novel inhibitors are worthy of further assessment for safety and efficacy in vitro and in vivo.

Conclusion
In the present work, we established a pharmacophore model to recognize vitally assorted lead hits for TNFα, TNFR1 and the TNFα-TNFR1 complex. The recognized hit compounds were utilized to create novel, strong inhibitors for the targets, and further assessed by docking and in silico ADMET studies. Fifteen lead compounds satisfied all the criteria and serve as novel, structurally diverse inhibitors for TNFα, TNFR1 and the TNFα-TNFR1 complex.

Additional file
Additional file 1: Fig. S1. FAF-Drugs4 ADME results for the TNF-α best ligand molecules and their respective properties such as: 2D structure of each ligand atoms, physicochemical filter positioning, compound complexity, oral property space, oral absorption estimation and Pfizer 3/75 rule positioning. Fig. S2. FAF-Drugs4 ADME results for the TNFR1 best ligand molecules and their respective properties such as: 2D structure of each ligand atoms, physicochemical filter positioning, compound complexity, oral property space, oral absorption estimation and Pfizer 3/75 rule positioning. Fig. S3. FAF-Drugs4 ADME results for the TNF-α-TNFR1 complex best ligand molecules and their respective properties such as: 2D structure of each ligand atoms, physicochemical filter positioning, compound complexity, oral property space, oral absorption estimation and Pfizer 3/75 rule positioning. Table S1. TNF-α and its inhibitors to compute physicochemical descriptors as well as to predict ADME parameters, pharmacokinetic properties, druglike nature and medicinal chemistry friendliness properties predicted by SwissADME tool. Table S2. TNFR1 and its inhibitors to compute physicochemical descriptors as well as to predict ADME parameters, pharmacokinetic properties, druglike nature and medicinal chemistry friendliness properties predicted by SwissADME tool. Table S3. TNF-α -TNFR1 complex and its inhibitors to compute physicochemical descriptors as well as to predict ADME parameters, pharmacokinetic properties, druglike nature and medicinal chemistry friendliness properties predicted by SwissADME tool.