Molecular dynamics simulation and docking studies reveal NF-κB as a promising therapeutic drug target for COVID-19

Coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is rampant worldwide and is a deadly disease for humans. Our current work emphasizes on molecular dynamics simulation (MDS) targeting nuclear factor-kappa B (NF-κB), the well-known human transcription factor controlling innate and adaptive immunity, to understand its mechanism of action during COVID-19 in humans. NF-κB was interacted with the SARS-CoV-2 spike protein in an in silico MDS experiment, revealing the NF-κB site at which the SARS-CoV-2 spike protein interacts. We screened some known drugs via docking studies on NF-κB used as a receptor. The MDS software Schrodinger generated more than 2000 complexes from these compounds and using the SMILES format of these complexes, 243 structures were extracted and 411 conformers were generated. The drug used as a ligand that docked with NF-κB with the best docking score and binding anity was Sulindac sodium as its trade name. Furthermore, RMSF data of sulindac sodium and NF-κB displayed minimal uctuations in the protein structures, and the protein-ligand complex had reduced exibility. Sulindac sodium is hence suggested as a suitable drug candidate for repurposing in clinical trials for SARS-CoV-2 infections. This drug potently blocked the spike protein’s interaction with NF-κB by inducing a conformational change in the latter. Arguably, NF-κB inaction is desired to have normal immunity and can possibly be retained using proposed drug. This work provides a signicant lead for drug repurposing to combat SARS-CoV-2 and its various mutant forms and reveals new approach for controlling SARS-CoV-2-induced disease.


Introduction
In a normal state of physiological activity, our immune system maintains homeostasis in the expression levels of NF-kappa B (NF-κB), which remains in a steady state by staying in a complex of proteins (such as IκBs), not allowing NF-κB to dissociate and translocate to the nucleus from the cytosol to exert its activity. However, stimuli trigger the activation of NF-κB through phosphorylation of IκBs by IκB kinase (IKK), which leads to nuclear translocation of NF-κB, where it binds to its cognate DNA and activates transcription of a wide variety of genes involved in host immunity, in ammation, cell proliferation and apoptosis (Oeckinghaus and Ghosh 2009). Various NF-κB inducers are known to be highly variable and include aggregates of bacterial lipopolysaccharides, ionizing radiation, reactive oxygen species (ROS), cytokines such as tumour necrosis factor alpha (TNF-α) and interleukin 1-beta (IL-1β) and viral DNA and RNA . Upon activation, NF-κB promotes the gene expression of an broad range of cytokines (e.g., IL-1, IL-2, IL-6, IL-12, TNF-α, LT-α, LT-β, and GM-CSF), chemokines (e.g., IL-8, MIP-1, MCP1, RANTES, and eotaxin), adhesion molecules (e.g., ICAM, VCAM, and E-selectin), acute phase proteins (e.g., serum amyloid A: SAA), and inducible effector enzymes (e.g., inducible nitric oxide synthase: iNOS and cyclooxygenase-2: COX-2). NF-κB can be considered a "quick action" primary transcription factor that is able to regulate myriad cellular responses, including the host's early innate immune response to infection, and is also associated with chronic in ammatory states, viral infections, septic shock syndrome and multiorgan failure Li and Verma 2002). The constitutive activation of NF-κB pathways has been found to be stimulated during in ammatory diseases, such as multiple sclerosis and rheumatoid arthritis (Liu et al. 2017).
Since coronavirus disease 2019  caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and became a pandemic worldwide, research on immediately available drugs that can block the causative virus or boost the human immune system to counter this virus has become rampant. With a lack of adequate speci c treatment for COVID-19, there is urgency to develop or repurpose drugs to help end the epidemic soon. The World Health Organization launched trials to investigate four potential treatments: remdesivir; chloroquine/hydroxychloroquine; lopinavir and ritonavir; and lopinavir and ritonavir plus interferon-β.1 (Neufeldt et al., 2020). However, these kinds of trials need to include thousands of patients from several countries and hence will consume a large amount of time and money. Furthermore, countries have developed their own vaccines, but proving their worth still requires further monitoring. We therefore need to urgently develop more focused therapy before the virus begins killing humans in an exponential manner.
The drug discovery process has received a great impetus with the emergence of computer-aided drug discovery (CADD) and has been instrumental over the last decade in exploring protein inhibitors in protein-drug interactions and protein-protein interactions (Keretsu et al., 2020). Since the process involved in the development of a candidate drug into an approved drug is lengthy and expensive, a combination of computational methodologies such as virtual screening, docking, molecular dynamics simulation (MDS), and binding free energy evaluation contributes to identifying potential drug candidates from compound libraries (Elmezayen et al., 2020). To date, MDS and structure-based virtual screening studies have been carried out to understand the SARS-CoV-2 spike protein function (Arantes et al., 2020;Selvaraj et al., 2020), role of repurposed protease inhibitors in blocking SARS-CoV-2 (Cardoso and Mendanha, 2021), effects of some antiviral drugs on SARS- CoV-2 (Zhao et al., 2020), and effect of temperature on the structure of the spike protein (Lipsa Rath and Kumar, 2020) and to perform molecular docking and simulation studies on SARS-CoV-2 elucidating some potent drugs that can be effective in controlling COVID-19 (Lokhande et al., 2020), as well as many other aspects of SARS-CoV-2 biology.
However, the rate of occurrence of different pathological human phenotypes and their heterogeneous lethality rates arising from constant mutations of SARS-CoV-2 (Martinez et al., (2020)) indicate serious problems that need to be appropriately addressed. A recent study revealed the importance of the NF-κB pathway in developing therapy regimens for critical COVID-19 patients (Hariharan et al., 2021). Other reports have also substantiated the above fact by arguing that therapeutic bene ts of NF-κB inhibitors, including dexamethasone, a synthetic form of glucocorticoid, have increasingly been realized (Elkhodary, 2020). Abnormal activation of NF-κB resulting from SARS-CoV-2 infection might be associated with the pathogenic pro le of immune cells, cytokine storms and multiorgan defects (Elkhodary, 2020).
Motivated by the mentioned facts, we currently explored the structural details of NF-κB via an in silico approach wherein we pinpointed the protein pockets for the binding of SARS-CoV-2 spike proteins. While performing its function upon entering the human body, this viral protein may trigger the action of the NF-κB transcription factor directly or may activate IKK, the kinase phosphorylating IκBs that otherwise remain bound with NF-κB to keep it inactive and in the cytoplasm. IκBα is known to be phosphorylated by speci c kinases at two sites near the N-terminus (Ser-32 and Ser-36), while the phosphorylated protein is then ubiquitinated at Lys-21 and Lys-22, leading to proteosome-mediated degradation .
Removal of IκB unmasks the nuclear localization sequence (NLS) of NF-κB and allows its translocation to the nucleus. Thus, IκBα is a multifunctional inhibitor of NF-κB that blocks nuclear translocation, DNA binding, and phosphorylation by PKA. Inhibition of NF-kB inhibits the NF-κB-mediated deactivation of immune and in ammatory responses to stimuli such as cytokines or bacterial/viral infection products . We further report here the NF-κB amino acids that interact with the viral spike protein to show the precise protein pocket location on NF-κB (Fig. 4B). Concomitantly, our selected drug named sulindac sodium (molecular formula: C20H16FO3S) from a total of 411 conformers was generated using the ligand preparation wizard and docked at the active site (used site map for identi cation) using Glide. The drug currently chosen has a prior history of use as therapy in various human diseases (Brunell et al., 2011;Marchetti et al., 2009) and was currently subjected to MDS together with the SARS-CoV-2 spike protein and NF-κB separately ( Fig. 1.). Interestingly, our molecular dynamics simulation data revealed some interesting features. Sulindac sodium was able to block the interaction between the spike protein and NF-κB at its designated pocket, where the former normally interacts with NF-κB and leads the spike protein to interact at another (possibly non-speci c) site on NF-κB. This could possibly keep NF-κB inactive in the cytoplasm. An MDS study on the NF-κB pocket showed the presence of most of the residues (Fig. 4C) that interact with the spike protein; thus, drug binding to the NF-κB-spike protein complex interferes with the interaction of the spike protein complex. In the presence of the drug molecule, the spike protein binds to the NF-κB-IkB complex using a pocket other than that used when it binds in the absence of the drug molecule. Furthermore, our analysis revealed that in the presence of the drug molecule, the spike protein failed to interact with IkB.
We further investigated the changes in the 3-D structure of the NF-κB-IkB complex due to the binding of the drug molecule and the binding of the spike protein in the presence of the drug molecule compared to the unbound complex. It was found that due to binding of the drug molecule, the compactness of the complex increased. The compactness further increased due to the binding of the spike protein in the presence of the drug molecule. The increase in compactness can be understood as an increase in the interaction between the chains of the NF-κB-IkB complex, which further supports the nding that in the presence of the drug molecule, the inactivity of NF-κB is further enhanced. We also employed MM-GBSA (molecular mechanics energies combined with the Poisson-Boltzmann or generalized Born and surface area continuum solvation) to estimate ligand binding a nities for sodium sulindac and receptor NF-κB and the values found as -2.15 kcal/mol for MMGBSA dG bind and -2.71 kcal/mol for MMGBSA dG bind(NS) show strong interaction between the two.

A. Protein and Ligand Preparation
Coordinates of the X-ray diffraction-based crystal structure of the I-κ-B-α/NF-κB complex with solvents (PDB ID: 1NFI) at a resolution of 2.70 Å were downloaded from the Protein Data Bank (PDB). The PDB IDs of the structures used were recorded. Solvent molecules and chains A, B, and F were removed during protein preparation (Fig. 4A). using the Protein Preparation wizard from Schrodinger Maestro (Schrödinger, LLC, New York, NY). The remaining structures were rst preprocessed using the Protein Preparation wizard, which included proper assignment of bond order; adjustment of ionization states; orientation of disoriented groups; creation of disul de bonds; removal of unwanted water molecules, metal and cofactors; capping of the termini; assignment of partial charges; and addition of missing atoms and side chains. In the case of the SARS-CoV structure, a loop (residues 376-381) missing in the PDB structure was modelled using Schrödinger's Prime module. Hydrogen atoms were incorporated, and a standard protonation state at pH 7 was used. Bond orders were assigned using the chemical components dictionary (CCD) database. The heterostate was generated using Epik with a pH of 7 (+/-2.0), and Prime was used to ll the missing chains and loops. Furthermore, all water molecules were removed in updated structures. While re ning the structure, PROPKA with pH 7.0 was used along with the sample water orientation. While minimizing the structure, a root mean square deviation (RMSD) of 0.30 Å with the OPLS_2005 force eld was used. The ligand was selected based on published information regarding its use as a therapeutic molecule for various human diseases and its e ciency in yielding positive results with AutoDock-and Swiss Dock-based docking with NF-κB (Kumar et al., 2021, unpublished data). The ligand was prepared using the Ligprep wizard of Schrodinger with the OPLS_2005 force eld; the Epik job was submitted with the metal-binding site and included the original state. Thirty-two stereoisomers were assigned per ligand with speci ed chirality. Glide was used to lter the search to locate the ligand in the active-site region of the receptor. The shape and properties of the receptor were represented on a grid to provide a more accurate scoring of the ligand poses. The docked complexes were superimposed on the original crystal structure to calculate the RMSD.

B. Virtual Screening and Molecular Docking
The Glide grid generator was used for generating the grid for blind docking. A Schrodinger virtual screening work ow (VSW) was used to score the virtual screening with default parameters employing the Glide program of Schrödinger. Use of HTVS mode allowed the elimination of most of the stereoisomers, and only 10% of the total ligands could be retained that passed the screening. Ligands following QikProp and Lipinski's rule were ltered for docking. Docking in the VSW was performed with Epik state penalties and further passed through HTVS, SP and XP docking modes. The ligands that showed the best a nity towards the main protease quali ed for our study. Finally, the interactions of the selected ligand and protein docked complexes were analysed by a poseviewer.

C. Molecular Dynamics Simulation
To study the dynamic behaviour of the protein complex under simulated physiological conditions, MDSs of the protein-ligand (PL) complex were performed using Desmond, which is available with Schrodinger Maestro (v12.5). The PL complex (9873 atoms) was solvated in a 10 × 10 × 10 Å orthorhombic periodic box built with TIP3P water molecules (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983). The whole system was neutralized by adding an appropriate number of 6 Na + counterions. This solvated system was energy minimized and position restrained with OPLS_2005 as the force eld (Jorgensen and Tirado-Rives 1988). Furthermore, 100 ns of MDS was carried out at 1 atm pressure and 300 K temperature, implementing the NPT ensemble with a recording interval of 100 ps, resulting in a total of 1000 read frames. Finally, various parameters of the MDS, such as ligand-binding site analysis, RMSD, root mean square uctuation (RMSF), PL contacts, secondary structure element (SSE) analysis, etc., were also analysed to check the stability, compactness, structural uctuations and PL interactions in the solvated system.

Structural Data:
The PDB was used to extract all 3-D structural information. The X-ray diffraction-based crystal structure of the I-κ-B-α/NF-κB complex with solvents (PDB ID: 1NFI) at a resolution of 2.70 Å was used for this study. Solvents and chains A, B, and F were removed during protein preparation (Fig. 4A). The structure of the SARS-CoV-2 spike glycoprotein (closed state) solved using electron microscopy (PDB ID: 6VXX) with a resolution of 2.80 Å was obtained (Fig. 4B.). Furthermore, all water molecules were removed from both structural data sets. The educational version of PyMOL (The PyMOL Molecular Graphics System, v1.2r3pre, Schrödinger, LLC) was used to generate the images derived to understand and analyse the structure and interchain interaction information.

Pocket Analysis:
NF-κB has multiple pockets, as calculated using the CASTp 3.0 server (Liang et al., 2018). Out of these pockets, the three largest pockets (based on the accessible area) were considered for analysis. The residues present in these pockets were considered interacting residues (Fig. 4). The area of the pockets and the interacting residues (exposed) changed signi cantly due to interaction with the spike protein and the drug molecule or the combination of the two.

Molecular Docking:
The Schrodinger VSW was used for virtual screening with default parameters using the Glide program. HTVS mode eliminated most of the stereoisomers, and only a few of the isomers passed after screening. QikProp and Lipinski's rule were further used to select the ligands without the functional group that were screened beforehand. Ligands were already prepared using the Ligprep wizard of Schrodinger with the OPLS_2005 force eld. Less than or equal to 32 stereoisomers were assigned per ligand with speci ed chirality. Docking in the VSW was performed with Epik state penalties and further passed through HTVS, SP and XP docking modes. The ligands that showed the best a nity towards the main protease quali ed. Finally, we found 1 unique ligand that showed the best docking score and binding a nity (Fig.  3). The complex of the drug molecule docked in the pocket of NF-κB was designated S0 (Fig. 4A) and was taken as the reference for further interaction comparative analysis.

Protein-Protein Docking:
Protein-protein docking studies were carried out using ClusPro 2.0 (Vajda et al., 2017) to understand the interaction of the SARS-CoV-2 spike glycoprotein with NF-κB in the presence and absence of the drug molecule. The best pose of the spike protein binding with the NF-κB-IkB complex was selected based on the cluster size, which is how the models are ranked in ClusPro. S1 (Fig. 4B) represents the complex of NF-κB with the SARS-CoV-2 spike glycoprotein, and S2 (Fig. 4C) represents the complex of S0 with the SARS-CoV-2 spike glycoprotein.

Interaction Studies
LigPlot+ v.2.2 (Laskowski and Swindells, 2011) was used to nd the ligand-protein interaction of complex S0 (Fig. 2). PDBsum was used to calculate the interaction between NF-κB and IkB (Fig. 5B), the interaction between the spike protein and NF-κB (Fig. 5B) and the interaction between the spike protein and NF-κB docked with drug molecules. PyMOL was used to analyse the structural differences imposed on NF-κB when it interacted with the drug molecule and the spike protein in the presence of the drug molecule. NS here refers to binding or interaction energy without any involvement of conformational changes occurring in the receptor or ligand. The potential energy of the complex will be reported in kcal/mol.

MD Simulation Studies
MDSs of the PL complex were carried out using Desmond 6.1 (Maestro v12.5), which provides information about the receptor-ligand complex over time. Thus, we performed MDSs for 100 ns on the complex after docking. After the simulation, we analysed the trajectory les for RMSD, RMSF, PL interactions, etc.

Root Mean Square Deviation (RMSD)
The RMSD is an indicator that describes the average displacement of an atom in a speci c molecular conformation concerning a reference conformation. In trajectory analysis, the complex RMSD in the initial phase at 1.63 Å went to 4.11 Å while stabilizing the structure for 100 ns of simulation. Initially, at up to 50 ns, the RMSD did not uctuate much, but after that, the stability uctuated for 10 ns and again started stabilizing. The c-α RMSD uctuated initially up to 1.47, after which a small uctuation was noticed (Fig.   1). Local changes along the protein chain were calculated throughout the simulation deliberations, and the RMSF was studied to characterize such changes. The RMSD of each residue of a conformation in a frame relative to the average conformation was studied to determine the exibility of a region of the protein. The peak in the RMSF plot indicated the protein's maximum uctuation during the simulation work, and the lower RMSF values represent fewer conformational changes. Overall, the RMSF plot analysis displayed minimal uctuations in the protein structures (Fig. 2). Furthermore, minimal uctuations were also observed in the PL complex, and the RMSF plot showed uctuations in some regions of the protein. The results of the residue interaction analysis after docking are highlighted with speci c colours based on the features (Fig. 3).
We assumed that the spike protein of SARS-CoV-2 directly interacts with IkB and dislocates this molecule from the IkB and NF-κB complex. However, this hypothesis was not found to be true based on the docking outcomes of the spike protein and NFkB (Fig. 4A). We further identi ed the ligand that could strengthen the interaction between IkB and NF-κB so that the spike protein of SARS-CoV-2 could not activate NF-κB for nuclear translocation, remaining in the inactive form in the cytoplasm, as shown in Fig. 4C. We studied the docking interactions of the identi ed drug and found that it docked in the largest pocket of the NF-κB-IkB complex that contains most of the residues interacting with the spike protein. It is possible that drug binding to the complex interferes with the above interaction of the spike protein. In the presence of the drug molecule, the spike protein bound to the NF-κB-IkB complex using a pocket other than that used when it binds in the absence of the drug molecule (Fig. 4C). Additionally, the drug molecule did not allow the spike protein to interact with IkB.
The changes in the 3-D structure of the NF-κB-IkB complex due to the binding of the drug molecule and spike protein in the presence of the drug molecule were compared to that of the unbound complex (Fig.  5A). The binding of the drug molecule to the NF-κB-IkB complex led to an increase in the compactness of the complex. The area and volume data regarding the pockets presented in Fig. 5A corroborated this observation. The binding of the spike protein in the presence of the drug molecule further increased the compactness of this complex. The increase in compactness can be understood as the increase in the interaction between the chains of the NF-κB-IkB complex, which further supports the nding that in the presence of the drug molecule, the inactivity of NF-κB is enhanced.
MM-GBSA studies are performed for docked complexes to accurately measure relative binding energy/a nity. MMGBSA dG bind values have a reasonable agreement with ranking based experimental scores, i.e., binding values (Uttarkar and Niranjan, 2021 Fig. 7. At present, our understanding of the involvement of the NF-κB signalling pathway in COVID-19 is limited, but NF-κB inhibitors have been increasingly utilized to gain therapeutic bene ts in many human diseases (Park and Hong, 2016;Yu et al., 2020). Abnormal activation of NF-κB in the manifestation of SARS-CoV-2 infection is considered associated with the pathogenic pro le of immune cells, cytokine storms and multiorgan defects. Thus, pharmacological inactivation of the NF-κB, such as proposed currently in the use of sulindac sodium, NF-κB-kinase interaction either directly or by inhibiting the dephosphorylation of its inhibitors, such as IKKα, IKKβ and other associated molecules, strongly represents a potential therapeutic target to treat the symptoms of COVID-19 in clinical trials.
Declarations Figure 1 Root mean square deviation (RMSD) of the protein and ligand after the initial RMSD values were stabilized. This plot shows RMSD values for the protein on the left Y-axis, whereas for the ligand, these values are indicated on the right Y-axis. The RMSD graph for c-αis shown in blue colour, and for ligand t on the protein, it is in red colour.   Interactions that occur more than 7.0% of the simulation time in the selected trajectory from 0.00 through 100.00 ns. The highlighted structure represents the ligand (drug) molecule. The interacting residues are shown in circular shapes with the representation of chains as C (representing residues of the NFkB molecule) and E (representing residues of IkBα). Circles with H2O represent the water-mediated H-bond interactions. Fig.   3B. A timeline representation of the interactions and contacts (H-bond, hydrophobic, ionic, and water bridge) summarized in the above gure. The top panel shows the total number of speci c contacts that the protein makes with the ligand over the course of the trajectory. The bottom panel shows speci c residues interacting with the ligand in each trajectory frame. Some residues make more than one speci c contact with the ligand, which is represented by a darker shade of orange according to the scale shown to the right of the plot. The interaction made by the drug molecule with the residues of the active site cavity of NF-kB, IkB and P50 peptide chains is dynamic in nature and formed, broke and reformed during the simulation duration (Fig. 3C). A brief duration of 20 ns was observed when no interaction occurred between the drug molecule and the NF-kB active site residues. Although there was an interaction between the ligand (drug) molecule and residues of NF-kB (chain C) and IkB (chain E) in the rst 50 ns, the interaction between the drug and IkB (chain E) disappeared after 50 ns. This loss of interaction strengthened the ligand-NFkB interaction. The binding pocket is represented by the labelled residues (G209, D210, E211, D257, R253, R259, S269, E285, etc., of the NFkB molecule and I102, N182, Q183, H184, etc., of the IkBα molecule). Fig. 4B.
Structure of the complex formed after protein-protein docking of the SARS-CoV-2 spike glycoprotein (closed state) (PDB ID: 6VXX, shown in orange colour at a resolution of 2.80 Å) with the NF-κB molecule. The highlighted region represents the binding pocket of the spike protein in the NF-κB molecule. The binding pocket is represented by the labelled residues (G209, D210, E211, D257, R253, R259, S269, E285, etc., of the NFkB molecule and N182, Q183, H184, etc., of the IkBα molecule). This information indicates that the drug molecule binds in the same pocket that the spike protein used to interact with NF-κB, possibly to activate its translocation to the nucleus. Fig. 4C. Structure of the complex formed after protein-protein docking of the SARS-CoV-2 spike glycoprotein (shown in orange colour) with the NF-κB molecule (already docked with the drug molecule). The highlighted region represents the binding pocket of the spike protein in the NF-κB molecule in the presence of the drug molecule. The drug molecule occupies the same binding pocket on NF-κB that is also used by the spike protein during docking interaction, thus inhibiting the spike protein from interacting with NF-κB at its designated site. In the interim, the spike protein interacts with NF-κB at a different site and hence may not be able to activate it.
Currently, these assumptions need validating experimental proof, which is lacking because of the strict restraint on experimentation with live samples from active human COVID-19 cases.

Figure 5
A 3-D structural alignment of the NF-κB -drug complexes (represented in a red cartoon model) and spikeprotein-NF-κB-drug complexes (represented in a green cartoon model), with the NF-κB molecule represented in an orange cartoon model. Fig. 5A. These alignments show the conformational changes imposed in the structure of the NF-κB molecule upon binding of the drug and spike protein. Fig. 5B.
Sequence-based comparison of the interacting residues of the NF-κB molecule with the drug sulindac sodium (yellow coloured bar), IkB molecule (green coloured bar) and spike protein (red coloured bar).

Figure 6
Comparison of the residues of the spike protein, MAP kinase and pLC gamma interacting with the residues of the NF-κB molecule. The underlined residues are in common (black: common in the s protein and pLC gamma, red: common in the spike protein and MAP kinase, and blue: common in all three).