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

Abstract

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.

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α 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.

Methods

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.cbligand.org/downloads/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://zincpharmer.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 MolSoft Drug-Likeness explorer (http://www.molsoft.com/mprop/) and the FAF-Drugs4 server (http://fafdrugs4.mti.univ-paris-diderot.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]. SwissADME predicts BBB (blood brain barrier) penetration and GI (gastro intestine) absorption by BOILED-Egg method [20]. It classified compounds as targets of p-glycoprotein (p-gp) efflux, inhibitors of cytochrome P450 enzymes CYP2C9, CYP2C19, CYP2D6 and CYP3A4 and substrates for metabolism by CYP2D6 and CYP3A4. It has predicted drug likeness by Lipinski, Ghose, Veber, Egan, Muegge methods and medicinal chemistry parameters by the Pan-Assay Interference Compounds (PAINS), Brenk methods and other parameters.

Results

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.

Fig. 1

The detailed work flow of the present study. The cheminformatics part of the pipeline part of the pipeline indicates pharmacophore modeling, virtual screening, molecular docking and in silico ADMET analysis

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.cbligand.org/downloads/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.

Active site prediction

The active amino acid sites of TNFα, TNFR1 and the TNFα–TNFR1 complex were predicted using the CASTp server. CASTp identified pockets, pocket mouth openings and the cavities of TNFα, TNFR1 and the TNFα–TNFR1 complex. TNFα had a pocket ID of 2, an area of 104.507 and a volume of 35.048. TNFR1 had a pocket ID-2, an area of 24.807 and a volume of 26.187. The TNFα–TNFR1 complex had a pocket ID of 3, an area of 521.964 and a volume of 26.187. The CASTp server predicted the binding site amino acids of TNFα (Val91, Asn92, Leu93, Phe124 of chain-A, His15, Val17, Ala18, Pro20, Arg32 Ala33, Asn34, Ala35, Phe144, Glu146, Ser147, Gly148, Gln149 and Val150 of chain-B) [14], TNFR1 (Ser74, Lys75, Arg77, Asn110 and Leu111 of chain-A), and the TNFα–TNFR1 complex (His15, Val17, Ala18, Pro20, Arg32, Ala33, Asn34, Ala35, Phe144, Ala145, Glu146, Ser147, Gly148, Gln149, Val150, and Tyr151 of chain-A, Thr77, His78, Thr79, Ser81, Pro-90, Val91, Asn92, Leu93, Ser95, Ile97, Phe124, Glu135, Ile136, and Asn137 of chain-B, Phe60, Leu71, Ser72, Cys73, Ser74, Lys75, Arg77, Gln82, Cys96 and Leu111 of chain-R for binding interactions to the inhibitors). The active sites of the three target proteins are presented in Fig. 2. These active site amino acids play vital roles in the pathological consequences of the dysregulated TNF-TNFR1 signaling pathway.

Fig. 2

The active sites of TNF-α, TNFR1 and the TNF-α–TNFR1 complex. TNFα has a pocket ID of 2, an area of 104.507 and a volume of 35.048. The TNFR1 has a pocket ID of 2, an area of 24.807 and a volume of 26.187. The TNFα–TNFR1 complex has a pocket ID of 3, an area of 521.964 and a volume of 26.187

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 (Fig. 3a): three aromatic spheres (pink), two hydrophobic spheres (green) and one hydrogen bond acceptor (HBA) (orange) taken into consideration. The hydrogen acceptor (orange) has a 0.50 radius, along with x = − 17.85, y = 78.31, z = 34.03, θ = 50.782 and φ = 70.610. The aromatic sphere1 (pink) has a 1.10 radius, along with x = − 18.80, y = 76.42, z = 35.78, θ = 0.0 and φ = 0.0. The aromatic sphere2 (pink) has a 1.10 radius, along with x = − 16.78, y = 73.23, z = 34.14, θ = 0.0 and φ = 0.0. The aromatic sphere3 (pink) has a 1.10 radius, along with x = − 21.71, y = 71.83, z = 34.02, θ = 0.0 and φ = 0.0. The hydrophobic sphere1 (green) has a 1.00 radius, along with x = − 17.18, y = 75.75, z = 38.08, θ = 0.0 and φ = 0.0. The hydrophobic sphere1 (green) has a 1.00 radius, along with x = − 15.71, y = 71.21, z = 36.46, θ = 0.0 and φ = 0.0. The TNFR1 pharmacophore contains five pharmacophoric features: Physcion-8-glucoside (ZINC33832439), one hydrophobic (green), two aromatic (pink), one hydrogen bond acceptor (HBA, yellow), and hydrogen bond donor (HBD, white) features taken into consideration (Fig. 3b). The HBD sphere (white) has a 0.50 radius, along with x = 0.25, y = 36.09, z = − 10.12, θ = 138.098 and φ = − 78.856. The HBA sphere (yellow) has a 0.50 radius, along with x = 1.05, y = 32.22, z = − 10.32, θ = 137.203 and φ = − 102.266. The hydrophobic sphere (green) has a 1.00 radius, along with x = 4.71, y = 31.01, z = − 5.13, θ = 0.0 and φ = 0.0. The aromatic sphere1 (yellow-green sphere) has a 1.10 radius, along with x = 2.76, y = 31.26, z = − 7.24, θ = 0.0 and φ = 0.0. The other aromatic sphere2 (green) has a 1.10 radius, along with x = − 1.67, y = 30.08, z = − 9.17, θ = 0.0 and φ = 0.0. The TNFα-TNFR1 complex pharmacophore contains five pharmacophoric features with Erythrosine B (ZINC08214556), one hydrophobic (green), two aromatic (pink), one HBA (yellow), and HBD (white) features taken into consideration (Fig. 3c). The HBD sphere (white) has a 0.50 radius, along with x = 49.89, y = 20.00, z = 45.40, θ = 38.557 and φ = 108.405. The HBA sphere (yellow) has a 0.50 radius, along with x = 44.42, y = 14.29, z = 46.48, θ = 59.686 and φ = 78.543. The hydrophobic sphere (green) has a 1.00 radius, along with x = 51.17, y = 14.39, z = 40.97, θ = 0.0 and φ = 0.0. The aromatic sphere1 (yellow white sphere) has a 1.10 radius, along with x = 48.35, y = 17.79, z = 44.79, θ = 0.0 and φ = 0.0. The other aromatic sphere2 (green) has a 1.10 radius, along with x = 44.04, y = 16.11, z = 43.27, θ = 0.0 and φ = 0.0. All the pharmacophore features play a vital role in screening for the best lead compounds.

Fig. 3

The pharmacophore features of TNF-α, TNFR1 and the TNF-α–TNFR1 complex. a The TNF-α pharmacophore contained six pharmacophoric features with default ligand (307), three aromatic spheres (pink), two hydrophobic spheres (green) and HBA (orange). b The TNFR1 pharmacophore contained five pharmacophoric features with Physcion-8-glucoside, one hydrophobic (green), two aromatic (pink), one HBA (yellow), HBD (white) features. c The TNF-α–TNFR1 complex pharmacophore contained five pharmacophoric features with Erythrosine B, one hydrophobic (green), two aromatic (pink), one HBA (yellow), HBD (white) features

Pharmacophore based virtual screening

The pharmacophore design models were used for screening molecules against the Zinc database (https://zinc.docking.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. Table 1 shows TNFα, TNF1 and the TNFα-TNF1 complex’s best top five small molecules’ ZINC IDs, popular name, SMILES (Simplified molecular-input line-entry system), binding energies, protein–ligand interaction residues, angles, distance between hydrogen bonds and number of hydrogen bonds. The ZINC09609430 molecule HN group interacts with one hydrogen bond to the Gly121 amino acid CH group with − 9.2 kcal/mol⁻¹ binding energy. The ZINC49467549 molecule may interact with electrostatic or Van der Waal bonds with − 9.0 kcal/mol⁻¹ binding energy to the Ile58, Leu120, Gly121, and Tyr151 active site amino acids of TNFα. The ZINC13113075 molecule OC group interacts with one hydrogen bond to the Tyr151 amino acid CO functional group with − 8.8 kcal/mol⁻¹ binding energy. The ZINC39907639 molecule may interact electrostatic or Van der Waal bonds with − 8.5 kcal/mol⁻¹ binding energy to the Ile58, Leu120, Gly121, and Tyr151 active site amino acids of TNFα. The ZINC25251930 molecule NC, NH and OC groups interacted with three hydrogen bonds to the Ile58, Leu120, Gly121, and Tyr151 amino acids CO, CO, and CN functional groups with − 8.1 kcal/mol⁻¹ binding energy respectively. The 307 (default ligand) molecule may interact with electrostatic or Van der Waal bonds with − 6.8 kcal/mol⁻¹ binding energy to the Ile58, Leu120, Gly121, and Tyr151 active site amino acids of TNFα (Fig. 4). Ile58, Leu120, Gly121, and Tyr151 are key residues for interacting with small molecules to inhibit TNFα trimer formation.

Table 1 The list of TNFα, TNFR1 and TNFα-TNFR1 complex and their best lead molecules interactions, binding energies, smiles, bond angles, bond lengths and number of hydrogen bonds

Fig. 4

The top five ligand molecular poses of docking protein–ligand interaction analysis of TNFα along with docking functional residues. The names of these compounds are ZINC09609430, ZINC49467549, ZINC13113075, ZINC39907639 and ZINC25251930

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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ binding energy respectively (Fig. 5).

Fig. 5

The top five ligand molecular poses of docking protein–ligand interaction analysis of TNFR1 along with docking functional residues. The names of these compounds are ZINC02968981, ZINC09544246, ZINC58047088, ZINC72021182 and ZINC08704414

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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ binding energy respectively (Fig. 6). The functional key residues play a vital role in TNFR1 and TNFα-TNFR1 complex formation and activation of the TNFα signaling pathway in proinflammation.

Fig. 6

The top five ligand molecular poses of docking protein–ligand interaction analysis of the TNFα-TNFR1 complex along with docking functional residues. The names of these compounds are ZINC05462670, ZINC35681945, ZINC23553920, ZINC05328058, and ZINC17206695

Chemical analysis of drug-likeness

All fifteen inhibitors were performed Lipinski “Rule of 5” and “drug-likeness” by MolSoft (https://www.molsoft.com/) and FAF-Drugs4 (http://fafdrugs4.mti.univ-paris-diderot.fr/) tools. The compounds showed Log P ≤ 5, relative molecular mass ≤ 500, range of HBA (hydrogen bond acceptors) ≤ 10, and range of HBD (hydrogen bond donors) ≤ 5 considering the best ligand molecules were used as drug leads for biological activity (Table 2). Lipinski’s Rule of five could be a rule of thumb designed for evaluating the drug likeness, or deciding whether a substance through a particular pharmacologic or biological action that may possibly create it a credible verbally energetic compound in humans. The results showed that fifteen TNFα, TNFR1 and TNFα–TNFR1 complex inhibitors, i.e., ZINC09609430, ZINC49467549, ZINC13113075, ZINC39907639, ZINC25251930, ZINC02968981, ZINC09544246, ZINC58047088, ZINC72021182, ZINC08704414, ZINC05462670, ZINC35681945, ZINC23553920, ZINC05328058 and ZINC17206695 satisfied the Rule of five and drug-likeness.

Table 2 TNF-α, TNFR1, TNF-α -TNFR1 complex inhibitors and their molecular properties and drug-likeness predicted by Molsoft and FAF-Drugs4

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.

Fig. 7

The ADMET properties (2D structure of each ligand atoms, physicochemical filter positioning, compound complexity, oral property space, oral absorption estimation and Pfizer 3/75 rule positioning) of a the best TNF-α inhibitor (ZINC09609430), b the TNFR1 inhibitor (ZINC02968981), and c the TNF-α–TNFR1 complex inhibitor (ZINC05462670)

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 anti-inflammatory 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 membrane-bound 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.