Transcription factor NF-κB as target for SARS-CoV-2 drug discovery efforts using inflammation-based QSAR screening model.

NF-κB is a central regulator of immunity and inflammation. It is suggested that the inflammatory response mediated by SARS-CoV-2 is predominated by NF-κB activation. Thus, NF-κB inhibition is considered a potential therapeutic strategy for COVID-19. The aim of this study was to identify potential anti-inflammation lead molecules that target NF-κB using a quantitative structure-activity relationships (QSAR) model of currently used and investigated anti-inflammatory drugs as the basis for screening. We applied an integrated approach by starting with the inflammation-based QSAR model to screen three libraries containing more than 220,000 drug-like molecules for the purpose of finding potential drugs that target the NF-κB/IκBα p50/p65 (RelA) complex. We also used QSAR models to rule out molecules that were predicted to be toxic. Among screening libraries, 382 molecules were selected as potentially nontoxic and were analyzed further by short and long molecular dynamics (MD) simulations and free energy calculations. We have discovered five hit ligands with highly predicted anti-inflammation activity and nearly no predicted toxicities which had strongly favorable protein-ligand interactions and conformational stability at the binding pocket compared to a known NF-κB inhibitor (procyanidin B2). We propose these hit molecules as potential NF-κB inhibitors which can be further investigated in pre-clinical studies against SARS-CoV-2 and may be used as a scaffold for chemical optimization and drug development efforts.

Introduction

Transcription factor NF-κB or nuclear factor kappa B is a central regulator of immunity and inflammation. Dysregulation of this remarkable family of proteins has been associated with the development of many diseases including those linked to autoimmunity, inflammation and viral infections [[1], [2], [3]]. Several human viruses such as HIV-1 and hepatitis B and C viruses have evolved strategies to target and modulate NF-κB signaling pathways to evade the immune response and to promote their survival [4,5]. Recently, there has been growing evidence that suggests that severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) mediates its hyperinflammatory systemic response by inducing specific activation of NF-κB in infected lung epithelial cells [6]. It is suggested that this predominant NF-κB-based inflammation leads way to an imbalance in the immune response and hence inability to mount an effective anti-viral control response, especially in patients with comorbidities and weaker immunity systems [6]. According to the latest WHO coronavirus data, there were over 153 million people confirmed as COVID-19 cases since the start of the global pandemic, leading to over 3 million deaths worldwide [7].

The role of NF-κB in SARS-CoV-2 infection has been highlighted by several researchers. A large-scale study which aimed at identifying potential druggable proteins for COVID-19 stated that there are two SARS-CoV-2 proteins –NSP13 and ORF9c– which target the NF-κB pathway [8]. Another study that analyzed the human coronavirus-host interactome network revealed that the NF-κB signaling pathway is considered a significant pathway for SARS-CoV-2 [9]. One study reported a potential agent that was shown to have anti-viral activity against the novel SARS-CoV-2 through inhibition of NF-κB in cell lines [10]. The link of NF-κB to the cytokine release syndrome and its potential for use as a therapeutic strategy for COVID-19 has been reviewed by Hirano and Murakami [11].

The major role of NF-κB activation in severe acute respiratory syndrome (SARS) was reported earlier for SARS-CoV, where NF-κB inhibition led to decreased inflammation and increased survival in mice infected with SARS-CoV [12]. Additionally, it was previously shown that the spike (S) protein of the SARS-CoV induced activation of the innate immunity and release of cytokines such as IL-8 via activation and nuclear translocation of NF-κB p65 in human peripheral monocyte macrophages [13]. Use of a NF-κB inhibitor suppressed the release of IL-8 in cells infected with human coronavirus 229E [13]. The up-regulation of IL-6 and TNF-α post-S protein treatment in murine macrophages was also dependent on NF-κB activation, specifically via the degradation of IκBα [14].

NF-κB activation can enhance the expression of hundreds of target genes in response to a diverse set of stimuli that include cellular stress such as acidic pH, physical stress including UV-light and ionizing radiation, proinflammatory cytokines namely IL-1β and TNFα, bacterial components such as lipopolysaccharides (LPS), viral molecules, parasites, and receptor ligands such as CD40 ligand [15,16]. These stimuli act via different pathways which have been described as the classical or canonical pathway, the alternative or non-canonical pathway and other atypical activation pathways [[17], [18], [19], [20]]. Activation of the classic pathway leads to the partial proteasomal degradation of the inhibitor of NF-κB, IκBα, allowing for the nuclear translocation of the p50-p65 subunits of NF-κB and expression of target genes [3,18,21,22]. This process leads to the upregulation of hundreds of important target genes that play a fundamental role in mounting an effective immune response in particular. This response includes upregulating adhesion molecules which are necessary for the recruitment of leukocytes; inducing the production of cytokines such as IL-1 [23] that are critical for the differentiation of cells such as T cell lymphocytes and macrophages into M1 or M2 subtypes; release of antimicrobial products such as neutrophilic extracellular traps; and hence effective initiation and orchestration of the innate and adaptive immune responses [15,16].

This diverse biological function can be explained by (i) the different protein structures of NF-κB which come from multiple families and their structural combinations; (ii) the different cell types they are expressed in; (iii) the diverse stimuli that induce their activation, and (iv) the complex crosstalk with a wide variety of other transcription factors, signaling proteins and signaling molecules such as reactive oxygen species and micro RNAs [15,19,20]. This complex array of interconnections creates a network that can enhance or inhibit the transcription of NF-κB genes and/or directly control the expression of NF-κB target genes. NF-κB is further tightly regulated by various mechanisms such as epigenetic modifiers and positive and negative feedback loops [15,20,24]. IκBα is tightly bound to NF-κB inhibiting its nuclear translocation and further transcriptional activity, but leaves room for constitutional activity of NF-κB [15]. NF-κB also inhibits its overactivation via certain molecules such as the A20 deubiquitinase [16], and via elimination of damaged macrophages via inducing intrinsic mitochondrial autophagy [23]. Any dysregulation of the homeostasis of this remarkable family of proteins has been associated with the development of inflammatory and autoimmune diseases, viral infections and even the progression to cancer [3,16,19,20,[25], [26], [27], [28], [29], [30]].

Thus, downregulation of NF-κB-mediated cellular responses is critically sought as a therapeutic strategy for inflammatory diseases in particular and may be used as a potential therapeutic strategy for COVID-19. It was recognized that the potential use of NF-κB inhibitors is in treatment of corticosteroid-resistant asthma and chronic obstructive pulmonary disease [31], as exacerbation of asthma attacks following bacterial or viral infections may be related to induced activation of NF-κB via the stimulation of TLRs [32]. Drugs such as corticosteroids and non-steroidal anti-inflammatory drugs (NSAIDs) have been found to be inhibitors of the NF-κB pathways [32,33]. Other antioxidants such Vitamins C and E, B-carotene and omega-3 fatty acids were found to indirectly inhibit the NF-κB pathways [15]. NF-κB inhibition was found beneficial in early sepsis as an example of an inflammatory-thrombotic disease state [15]. This stresses the fact that chronic inflammation is closely associated to increased risk of thrombosis and coagulopathy disorders [15].

The search for potent and safe NF-κB inhibitors is of paramount importance as it may represent a novel therapeutic strategy for many diseases including inflammatory diseases and viral infections. While the full structural details of NF-κB have not been characterized yet, computational studies have thus far used mathematical modeling and simulation tool to investigate NF-κB and its signaling pathways [34]. Predictions made by various in silico studies which were successfully complemented by wet lab studies have advanced our understanding of NF-κB signaling pathways [35]. In our previous work, we developed fragment-based energy-optimized pharmacophore models for the NF-κB/IκBα which can be used for library screening purposes [36]. The use of quantitative structure-activity relationships (QSAR) models as a basis of screening to identify lead ligands have recently been successful in our previous studies [37], and shows potential as a powerful virtual screening methodology that should be investigated for promising molecular targets. Thus, the aim of this study was to apply an integrated computational approach by starting with a QSAR model of anti-inflammation drugs to screen a comprehensive library of drug-like molecules with the goal of finding lead ligands that target the NF-κB/IκBα p50/p65 (RelA) heterodimer complex. These screening efforts were followed by molecular dynamics (MD) simulations and free energy calculations.

Methods

Library and protein preparation

Three different small molecule libraries from Specs (specs.net) were used in our high throughput virtual screening research study. The Natural Products (834 compounds), World Diversity Set (10,000 compounds) and Specs SC library (210,550 compounds) were downloaded from https://www.specs.net. The crystal structure of the p50-p65 NF-κB and IκBα heterodimer complex was taken from our previous study, where the N-terminal signal receiving domain of IκBα was structurally modeled [38]. The protein complex and small molecule libraries were prepared according to the protein preparation module of the Schrodinger's Maestro suite, where hydrogen bonds were added, side chains and loops fixed and disulfide bonds generated [39,40]. To mimic physiological conditions, PROPKA was used to generate protonation states of amino acids at pH 7.4 [41,42]. The protein complex structure was optimized using the OPLS3 forcefield [43].

Binary QSAR model analysis

MetaCore/MetaDrug from Clarivate Analytics® is a comprehensive platform used to derive several biochemical, physical and pharmacological properties about chemical compounds. The QSAR models have been developed using a set of compounds that exhibit experimental activity and function (positives) and compounds which do not exhibit any activity or function (negatives) in approximately equal numbers. For each QSAR model developed, a set of compounds (the training set) is used to build the model, whereas other compounds are used to test the model for its validity (the test/validation set). The training set includes compounds that are known to exhibit activity against a particular disease and/or condition including drugs in the market. Only QSAR models that possessed accuracy, estimated from its correlation coefficient (R2) and root mean squared error (RMSE), and which were found to have the highest statistical powers (i.e. specificity, sensitivity, accuracy and the Matthews Correlation Coefficient) were selected for use by the platform. In this study, the inflammation-based QSAR model was used to screen ~210,000 small molecules from the SPECS library in addition to ~834 from the Natural Products and 10,000 from the WorldDiversity libraries. The inflammation-QSAR mathematical model can predict and calculate potential anti-inflammatory activity of any compound based on its chemical structural descriptors with a sensitivity of 0.86, specificity of 0.84 and accuracy of 0.85. This QSAR model was constructed with the use of 598 training set and 93 test set compounds. A calculated predicted value greater than 0.5 indicates a potentially active target ligand which can be evaluated further. We also used QSAR models covering over 26 toxicities to predict the toxicity of our ligands. Similarly, a calculated predicted value greater than 0.5 indicates a potentially toxic compound which can be eliminated. More details about this platform can be found in our previous paper [36].

Molecular docking

The Glide Standard Precision (Glide/SP) is a grid-based docking methodology was used in this study [[44], [45], [46]]. Glide/SP settings were used as default. For this purpose, we generated a grid box for our target protein and allowed the rotatable amino acids to rotate their side chains in order to provide additional flexibility. The grid box was centered around the interaction site between NF-κB and IκBα (see Figure S16). In grid-box generation, inner and outer box sizes were used as 10 and 30 Å, respectively. Glide/SP was used to dock the top-382 filtered ligands into target protein.

Molecular dynamics (MD) simulations

MD simulations were prepared using the TIP3P solvent model. 0.15 M NaCl solution was used to mimic the physiological conditions. OPLS2005 force field and RESPA integrator were used in the simulations [47,48]. We conducted MD simulations at 310 K with Nose-Hoover temperature coupling [49] and constant pressure of 1.01 bar via Martyna−Tobias−Klein pressure coupling [50] in the Desmond program [47,51]. Simulation box shape was selected as orthorhombic and box size calculation methods was used as “buffer”. The simulation box size is calculated by using the given buffer distance between the solute structures and the simulation box boundary which was 10 Å from each dimension. For the 50 ns MD simulations, our systems consisted of around 91,000 atoms and around 26,660 water molecules. Other settings were used as default. We performed 1 ns (short) MD simulations for 382 molecules (in total 382 ns), followed by 50 ns (long) MD simulations for 11 selected lead molecules. We also ran both 1 ns and 50 ns MD simulations for one of the positive control molecules, procyanidin B2.

The molecular mechanics-generalized Born surface area (MM/GBSA) continuum solvation calculations

The MM/GBSA was used for calculating the binding free energies of studied compounds [52]. Hou et al. in their succeeding studies represented that re-scoring by MM/GBSA is an effective procedure to improve the predictions of docking methods [53]. As our aim was to propose novel compounds as NF-κB inhibitors in a relatively efficient and inexpensive technique, MM/GBSA approach was preferred [54]. The MM/GBSA tool was conducted in this study using Schrodinger's Prime module [48]. 100 trajectory frames were used to calculate the MM/GBSA from short (1 ns) MD simulations and 2000 trajectory frames were considered from the 50 ns (long) MD simulations. More details about the simulation protocol can be found in our previous studies [37,[55], [56], [57]].

The study's full methodology is summarized in Scheme 1 .

Scheme 1

The study's computational approach for screening libraries for the purpose of identifying drug-like molecules with potential anti-inflammatory activity.

Results and discussion

In this study, we used a QSAR model of inflammation to screen nearly 223,000 drug-like molecules from three different small molecule libraries. This mathematical model can predict the anti-inflammation therapeutic activity of any molecule based on their chemical and structural characteristics. The lead ligands from the QSAR model screening of all three libraries were selected based on having predicted anti-inflammation activity with a value greater than 0.80. This high threshold was selected to increase the likelihood of including ligands with highly predicted anti-inflammation activity based on their structure. This process was followed by filtering out predicted toxic compounds and including compounds that are nontoxic or have low toxicities. Compounds with predicted toxicity with a value less than 0.60 in a maximum of 4 toxicities were included for further analysis. 382 molecules were selected as potentially nontoxic molecules and were analyzed further. These ligands were docked into the prepared inflammation target protein (NF-κB/IkBα) and advanced to MD simulations. In total, we ran 1 ns MD simulations and performed MM/GBSA free energy calculations for 382 molecules (Fig. 1 , left). The top-11 ligands (2 standard deviation away from the mean of their MM/GBSA scores) were selected to undergo 50 ns MD simulations (Fig. 1, right). This approach allowed us to expand our pool of molecules to be evaluated in the MD simulations so that their stability, behaviors and nonbonded interactions with the NF-κB/IκBα complex structure can be analyzed over time, opposed to using pure rigid docking protocols. The MM/GBSA tool is also a powerful tool to predict the free binding energy of biological systems allowing for an analysis of protein-protein and protein-ligand interactions [54].

Fig. 1

(left) The average MM/GBSA scores for the 382 hit compounds identified from the virtual screening analysis. The mean value is −40.52 kcal/mol. The value of −58.29 kcal/mol (i.e., 2 SD (standard deviation) from the mean) was used cutoff value for the filtration of compounds. These compounds were used in long MD simulations. Plots were produced using GraphPad Prism v6.01, www.graphpad.com. (right) The MM/GBSA scores for the top-11 lead compounds. Average MM/GBSA scores of short and long MD simulations of compounds were compared.

All of the 11 selected ligands were predicted to have a high therapeutic activity against inflammation (Figure S1). Based on the analysis of other QSAR models that are available for other diseases in the MetaCore/MetaDrug platform, the majority of these 11 ligands were also predicted to have potential therapeutic activity in arthritis, asthma, cancer, depression, HIV, hyperlipidemia, migraines and obesity, disease states that are strongly linked to inflammation. This highlights that using one model to describe inflammation is not enough to encompass the other diseases that are considered inflammatory in their pathogenesis. However, the predicted anti-inflammatory activity of the procyanidin B2 by the model was noted to be low. There are several reasons that could explain this observation such as the nonspecific nature of the polyphenolic structure. Hence, procyanidin B2's nonspecific chemical structure may not fit the chemical descriptors used by the construction of QSAR model. It is also plausible that the difference seen could be due to procyanidin B2 acting at a different step in the NF-κB pathway, at a different binding site and/or having a different or additional mechanism of action in its inhibition of NF-κB [58]. Nevertheless, other QSAR models predicted its therapeutic activity for arthritis and asthma which are conditions in which NF-κB plays an important role in the pathophysiology of these diseases [3]. In comparison with our 11 hit ligands which had essentially no predicted toxicities, the control had significant predicted toxicities namely cardiotoxicity, liver necrosis, liver cholestasis, nephron toxicity as well as liver weight gain (Figure S2).

Analysis of the MD simulations and Gibbs free energy calculations allowed us to consider five of our 11 ligands as most promising as potential NF-κB inhibitors (Fig. 2 , Table 1 ). These ligands were: AF-399/32354064, AG-690/12890456, AK-968/11841158, AK-968/41926571 and AP-064/41252894 from the Specs small molecule library. The control inhibitor showed an average MM/GBSA score that is significantly weaker (ΔG of −35.14 kcal/mol) than the identified top-5 lead ligands throughout the 50 ns MD simulations indicating that it is less energetically favorable and structurally stable compared to the lead ligands. Additionally, all of the 11 hit ligands had better MM/GBSA scores than the control. This was also true throughout the 1 ns MD simulations compared to the lead molecules. These results highlight the higher level of conformational stability attained by the NF-κB complex with the discovered lead ligands. The MM/GBSA free energy analysis for all 11 hit ligands along with NF-κB inhibitor, Procyanidin B2, at the binding pocket of NF-κB/IκBα throughout the 50 ns MD simulations can be found in Figure S3.

Fig. 2

MM/GBSA binding free energy analysis for the hits along with NF-κB inhibitor (Procyanidin B2) at the active site of the NF-κB/IκBα during the long MD simulations.

Table 1

Structures and the average free energy scores of the top selected 11 molecules. These 11 molecules were selected based on a QSAR-based screening methodology followed by 1 ns MD simulations for all filtered molecules. The average MM/GBSA scores are taken from the 50 ns MD simulations.

The compound that has the lowest MM/GBSA score at the binding pocket of the NF-κB/IκBα complex was AG-690/12890456 as confirmed by the least root mean squared deviation (RMSD) fluctuations seen throughout the 50 ns MD simulations (Figure S4). The RMSD is primarily used to investigate the structural stability of the protein-ligand complex throughout the MD simulations. This molecule also had the lowest average Gibbs free energy score (average ΔG of −60.18 kcal/mol) throughout the 50 ns MD simulation, suggesting a more stable ligand conformation and a predicted higher affinity at the binding site. Additionally, as is seen in Figure S13 which shows the protein-ligand interactions, this molecule had made significant interactions with amino acid residues Asn276, Leu277 and Gln278 (observed in over 40% of the simulation time) and with Met279 (observed in over 80% of the simulation time). Other interactions with the protein residues are also shown in the diagrams, but only when the interaction occurs greater than 30% of the simulation time. Additionally, Met279 was noted to make several contacts with the protein throughout the 50 ns MD simulations, leading to an interaction fraction greater than 1 when compared to other residues (Fig. 3 ). When analyzing the protein-ligand interactions for AG-690/12890456, it is noted that these contacts were made mostly via hydrogen bonding, although hydrophobic interactions and water bridges were also observed (Fig. 3). As described in our previous study, these residues are part of the PEST sequence that spans residues 276–287 in IκBα and which is considered important for identifying high affinity ligands [36]. Further details on the interactions between the protein and ligands at each of the 1 ns and 50 ns poses of the MD simulations are provided in the Supplementary Materials for the lead five ligands (see Figures S14 and S15).

Fig. 3

The protein-ligand interaction diagrams of the identified 5 compounds and the reference compound in the active site of the NF-κB/IkBα complex throughout the 50 ns MD simulations. The interaction fraction indicates the percentage of time the contact is made during the simulation. The profiles are displayed for ligands AF-399/32354064, AG-690/12890456, AK-968/11841158, AK-968/41926571 and AP-064/41252894 and Procyanidin B2, the control molecule.

The RMSD value of Cα atoms away from the initial positions is used to describe the flexibility of all possible protein conformers throughout the MD simulations. Our analysis showed that all the complex systems have less than average RMSD of 5 Å. Ligand AK-968/11841158-bound target protein had the highest fluctuations in RMSD values with an average of 6.69 Å although it reaches a plateau towards the last 30 ns. (Figure S5).

The “fit on protein” and “fit on ligand” modes represent the RMSD of the nonhydrogen atoms of ligands referenced to protein backbone and ligand itself, respectively. In the fit on protein mode, ligand AG-690/12890456 had the least observed mean values of RMSD (Figure S4). Although AF-399/32354064 had a higher fit on ligand RMSD values (Figure S6), it was relatively stable throughout the simulation unlike ligands which had more fluctuations seen over the first 20 ns. Although ligand AF-399/32354064 had an abrupt change in its rotational motion around 10 ns, it was one of the most stable ligands with the least fluctuations observed throughout the 50 ns MD simulations. Indeed, in the fit on ligand mode shown in Figure S6, all of the investigated compounds as well as the control molecule had average RMSD values less than 2 Å. This may indicate a limited rotational motion of the molecules at the investigated binding site.

We further analyzed the RMSD evolution of the backbone atoms covering the PEST sequence as this set of residues have been suggested to be important for targeting in NF-κB inhibition [17,59,60] (Figure S7). The selected hit molecules as well as the control molecule were stable throughout the MD simulations with respect to the PEST sequence (residues 276–287). Ligands AK-968/41926571 and AP-064/41252894 showed some fluctuations in the last 20 ns of the simulations, but this increase in the RMSD from a baseline of nearly 3.5 to 4.5 Å may be insignificant. Again, ligand AK-968/11841158 showed prominent fluctuations when compared to the other molecules, and continued oscillating throughout the 50 ns MD simulation, however its RMSD on average (3.27 Å) was similar with both ligands AK-968/41926571 (3.15 Å) and AP-064/41252894 (3.56 Å).

Analysis of the MD simulations in different parameters is represented in Figures S8-S11. Figure S8 shows the RMSD evolution for the side chains over a 50 ns MD simulations with the NF-κB/IκBα (p50/p65) complex. Highest side chain RMSD values was observed for ligand AK-968/11841158, which may indicate the level of perturbation of these side chain residues as they construct their interactions at the binding pocket. Figures S9 and S10 shows root mean square fluctuations (RMSF) evolution over time for the Cα atoms and side chains over a 50 ns MD simulation with the NF-κB/IκBα (p50/p65) complex, respectively. Figure S11 represents solvent accessible surface area (SASA) in Å2 of the 5 hit ligands over the 50 ns MD simulations with the NF-κB/IκBα (p50/p65) complex. The control procyanidin B2 showed the highest SASA value (328.76 Å2) compared to other ligands.

Based on our thorough analysis of the simulations and post-MD Gibbs free energy calculations, it can be suggested that ligand AG-690/12890456 is the most energetically favorable ligand at the binding site of NF-κB/IκBα. This ligand along with AF-399/32354064, AK-968/41926571 and AP-064/41252894 may show effective inhibition of NF-κB in in vitro studies and may be investigated further. ADME properties of these compounds were also calculated using MetaCore/MetaDrug and it is found that selected ligands have proper drug-like profiles (Figure S12). Use of these ligands as lead molecules and expanding their activity and purpose by chemical optimization and structural methods is also predicted to lead to the development of more potent inhibitors.

Conclusions

The COVID-19 pandemic highlights the critical need for rapidly and efficiently identifying important molecular targets in diseases and developing safe and effective drugs that could be used as novel therapies. In this regard, computer-aided drug design and discovery plays a critical role in efficiently searching massive databases of drug-like molecules and analyzing the chemical and physical characteristics of these potential compounds. The diversity and complexity of NF-κB pathways and interconnections necessitate a powerful understanding of its signaling pathways and connection to disease states in order to deduce the optimal patients that could benefit from targeting NF-κB, as well as to derive the optimal timing, dosing and mode of administration of these potential drugs [15,20].

Several research studies have presented the potential use of NF-κB inhibitors in isolation or in combination with other drugs for its maximum intended anti-inflammation benefits [20]. Thus, the search for potent and nontoxic NF-κB inhibitors is considered a highly desired strategy for treatment of inflammatory and inflammation-associated diseases and represents a novel way that has a strong potential for advancing to clinical use in the future. In our study, we have utilized an integrated computational approach to identify molecules predicted to have anti-inflammatory activity that interacts favorably with the NF-κB complex. We identified five hit molecules that could potentially act as NF-κB inhibitors. However, it is important to highlight that our study is limited by its in silico nature. Further studies to test these molecules in vitro and subsequently in vivo are necessary to evaluate their effectiveness in reducing inflammation marked by increased cellular overstimulation of NF-κB and its pathway mediators. NF-κB inhibitors discovered via in silico approaches can be investigated in pre-clinical studies and may be further considered as a scaffold for further structural optimization and drug development efforts.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.