Targeting SARS-CoV-2 nucleocapsid oligomerization: Insights from molecular docking and molecular dynamics simulations.

The outbreak of COVID-19 caused by SARS-CoV-2 virus continually led to infect a large population worldwide. Currently, there is no specific viral protein-targeted therapeutics. The Nucleocapsid (N) protein of the SARS-CoV-2 virus is necessary for viral RNA replication and transcription. The C-terminal domain of N protein (CTD) involves in the self-assembly of N protein into a filament that is packaged into new virions. In this study, the CTD (PDB ID: 6WJI) was targeted for the identification of possible inhibitors of oligomerization of N protein. Herein, multiple computational approaches were employed to explore the potential mechanisms of binding and inhibitor activity of five antiviral drugs toward CTD. The five anti-N drugs studied in this work are 4E1RCat, Silmitasertib, TMCB, Sapanisertib, and Rapamycin. Among the five drugs, 4E1RCat displayed highest binding affinity (-10.95 kcal/mol), followed by rapamycin (-8.91 kcal/mol), silmitasertib (-7.89 kcal/mol), TMCB (-7.05 kcal/mol), and sapanisertib (-6.14 kcal/mol). Subsequently, stability and dynamics of the protein-drug complex were examined with molecular dynamics (MD) simulations. Overall, drug binding increases the stability of the complex with maximum stability observed in the case of 4E1RCat. The CTD-drug complex systems behave differently in terms of the free energy landscape and showed differences in population distribution. Overall, the MD simulation parameters like RMSD, RMSF, Rg, hydrogen bonds analysis, PCA, FEL, and DCCM analysis indicated that 4E1RCat and TMCB complexes were more stable as compared to silmitasertib and sapanisertib and thus could act as effective drug compounds against CTD.Communicated by Ramaswamy H. Sarma.

Introduction

The current pandemic of COVID-19 (Coronavirus Disease-2019) caused by a novel coronavirus strain, 2019-nCoV/SARS-CoV-2 has led to over 27 million confirmed cases and more than 8 lakhs deaths in over 200 countries since its emergence in late 2019, in Wuhan, China, as of early-September, 2020 (https://covid19.who.int/). As the number of COVID-19 cases are continually increasing, there is a global need for antiviral drug development at an accelerated pace which brought drug repurposing as a promising strategy (Gil et al., 2020; Prasad et al., 2020). SARS-CoV-2 has a single-stranded, plus-sense, 30 kb RNA genome, which contains five major open reading frames encoding nonstructural replicase polyproteins (NSPs 1-16) and structural proteins, that is, spike (S), envelope (E), membrane (M), and nucleocapsid (N) (Zhu et al., 2020). SARS-CoV-2 nucleocapsid (N) protein is an essential conserved structural protein that binds to viral genomic RNA (gRNA) and plays an important role in viral transcription and translation (Baric et al., 1988). The N protein is also involved in binding other viral proteins including the E and M proteins to help the virus envelope formation and particle assembly (He, Leeson, et al., 2004; Kuo & Masters, 2002; Tseng et al., 2014). The N protein has an N-terminal RNA-binding domain (NTD), a C-terminal dimerization domain (CTD), and an intrinsically disordered central Ser/Arg (SR)-rich linker (McBride et al., 2014) (Figure 1).

Figure 1.

Graphical representation of the Nucleoprotein embedded N- & C- terminal domains, amino acid residues responsible for the dimerization of C-terminal domain and the In Silico work flow followed in the current study.

Both RNA binding and oligomerization of N protein are important for viral assembly (He et al., 2004; Nelson et al., 2000). The oligomerization of N protein is necessary for viral ribonucleoprotein (RNP) assembly that is crucial to protect the viral genome and sustain viral replication (Luo et al., 2006; Narayanan et al., 2003). The N-CTD has been shown to directly contribute to N protein dimerization and oligomerization (Ye et al., 2020; Yu et al., 2006). Furthermore, the N protein regulates the host innate immune response by inhibiting interferon β (IFN-β) production (Lu et al., 2011), and binds to host translation elongation factor 1α (EF1α) to inhibit host cell proliferation, mainly through the N-CTD (Zhou et al., 2008). N protein is also involved in different cellular activities, including the formation of stress granule and host translation shutoff (Raaben et al., 2007), inhibition of nonsense-mediated mRNA decay (Wada et al., 2018) and cell cycle regulation (Surjit et al., 2005). N protein is highly conserved and stable, with 90% homology and fewer mutations over time (Grifoni et al., 2020; Holmes & Enjuanes, 2003). It is highly immunogenic and is abundantly expressed during infection (Cong et al., 2020).

The first three-dimensional structure of the SARS-CoV-2 N-CTD was determined at 2.05 Å resolution and at pH 7.5 (PDB ID 6WJI; Center for Structural Genomics of Infectious Diseases (CSGID), unpublished) (Figure 1). Three-dimensional structures of N-CTD reveal a compact, tightly intertwined dimer with a central four-stranded β-sheet comprising the bulk of the dimer interface. This interface is composed of two β-strands and a short α-helix from each monomer that extend toward the opposite monomer and pack against its hydrophobic core. Many recent studies also demonstrated that N proteins of coronaviruses can be useful antiviral drug targets for its importance in the RNA genomic packing, viral transcription and replication (Chenavas et al., 2013; Lin et al., 2014; Lo et al., 2013; Monod et al., 2015; Sarma et al., 2020). These antiviral drugs either target the conserved RNA-binding sites to specifically block the formation of ribonucleoprotein during genome replication (Lejal et al., 2013; Tarus et al., 2015), or disrupt the oligomerization of N protein (Gerritz et al., 2011; Hung et al., 2012).

Very recently, Gordon et al (Gordon et al., 2020) presented a SARS-CoV-2-Human interactome to reveal potential drug targets against COVID-19. The authors identified 332 interactions between viral and host proteins, largely targeting the innate immune signaling pathway. The viral N protein targets stress granule protein G3BP1/2, the mTOR translational repressors LARP1, and the protein kinases CK2. The suppression of stress granules and host translation machinery is beneficial for viral replication, as stress granules inhibit the replication of viruses (Ivanov et al., 2019; Nakagawa et al., 2018; Raaben et al., 2007). Using this information, they identified a series of anti-N drugs with a high potential to fight COVID-19. These anti-N drugs inhibit the SARS-COV-2 proliferation in-vitro and the mechanism of action of these drugs are different like casein kinase-2 (CK-2) inhibitor, mTOR inhibitor, translation inhibitor, SG inhibitor, multi-targeted protein kinase inhibitor, etc. Moreover, these anti-N drugs could modulate the nucleocapsid assembly by interfering with CTD dimerization or oligomerization, and thus might be considered as valid therapeutic strategies against SARS‐CoV‐2.

Drug-repurposing is one of the well-versed techniques to understand the importance of active compounds to block the viral proteins during the process of infection. To explore the molecular mechanism of antiviral activity of anti-N drugs, here, a number of computational methods like molecular docking and all-atom molecular dynamics (MD) simulation has been performed on the SARS-CoV-2 CTD structure (PDB ID: 6WJI) (Figure 1). The crystal structure was prepared and the selected antiviral drugs such as 4E1RCat, Silmitasertib, TMCB, Sapanisertib, and Rapamycin were subjected to docking to acquire the better binding orientation. The overall results from the MD simulations revealed that the compounds 4E1RCat and TMCB have high stable confirmations and thus indicating their inhibitory activity toward the active cleft of the SARS-CoV-2 CTD, further hindering the virion assembly and interaction with the viral genome, which in turn halts the translation and replication processes. Thus, this mechanism-relevant study explored the binding conformations, spotted the important amino acid residues during the binding process, and also elucidated the detailed molecular level interactions between the protein-ligand complexes. The results anticipated insightful information into the interactive mechanism of the antiviral drugs toward CTD of N protein, which may pave a way for future exploration of efficient drug targets in COVID-19.

Methods

Molecular docking

Molecular docking approaches have broadly been used to know the binding modes of ligand and receptor and is widely used in drug discovery. We have used the crystal structure of the C-terminal domain of SARS-CoV-2 nucleocapsid protein (CTD) protein (PDB ID: 6WJI). The chemical structures of the drugs, 4E1RCat (ZINC000009311298), Silmitasertib (ZINC000058638454), TMCB (ZINC000058638742), Sapanisertib (ZINC000073069271), and Rapamycin (ZINC000169289388) were taken from the ZINC Database (Sterling & Irwin, 2015) as MOL2 files and converted into three dimensional PDB files and energy minimized using MM2 force field. Energy minimized 3 D structures were then converted into docking input file format (i.e. PDBQT). We have used Autodock (version 4.2), which is an extensively used automated tool for protein-ligand docking (G. M. Morris et al., 2009). The 3 D-grid box was set to be with dimensions 112 × 100 x 92 Å in XYZ axes, with a grid point spacing of 0.375 Å. The binding poses were generated by the Lamarckian Genetic Algorithm (Garrett M Morris et al., 1998). Ten binding poses were generated for each ligand. After docking, the 2 D receptor-ligand complexes were constructed using Schrödinger visualizer and visualized using PyMOL (DeLano, 2002).

Molecular dynamics (MD) simulations

The general methodology of MD simulations was carried out on unbound CTD and CTD-drug complexes using GROMACS-Version 5.18.3. The topology of CTD was generated by using GROMOS 9643a1 force field (Van Der Spoel et al., 2005). Due to the lack of suitable force field parameters for drug-like molecules in the GROMACS software package, PRODRG server (SchuÈttelkopf & Van Aalten, 2004) was used for generation of molecular topologies and coordinate files. The CTD systems were solvated using a simple point charge model (SPC/E) in a cubic box. The initial system of the protein and drug was set to a simulation grid size of 46.24, 42.89 and 33.14 on X-, Y- and Z-axis, respectively. The solvation system resulted 28,356 water molecules for protein and 29,486 water molecules for protein-drug complexes. To neutralize the system, 0.15 M counter ions (Na+ and Cl-) were added. The energy minimization of all the neutralized systems was performed using the steepest descent and conjugate gradients (50000 steps for each). The equilibration phase was carried out separately in NVT (constant volume) and NPT (constant pressure), each for 500 ps. The temperature of all systems was maintained at 300 K using Berendsen weak coupling method (Berendsen et al., 1987) and the pressure was maintained at 1 bar using Parrinello-Rahman barostat (Parrinello & Rahman, 1980). The SHAKE algorithm was used to confine the H atoms at their equilibrium distances, and periodic boundary conditions. Moreover, the Particle Mesh Ewald (PME) method used to define long-range electrostatic forces (Darden et al., 1993). The cutoffs for van der Waals and columbic interactions were set at 1.0 nm. LINC algorithm was used to constrain the bonds and angles. Using the NPT ensemble, production runs were performed for the period of 100 ps, with time integration. The energy, velocity, and trajectory were updated at the time interval of 10 ps. MD analysis of the root mean square deviations (RMSD), radius of gyration (Rg), solvent accessible surface area (SASA) and root mean square fluctuation (RMSF) was done using modules such as g_rmsd, g_gyrate, g_sasa, and g_rmsf, respectively.

Principal component analysis (PCA) and free-energy landscape

The principal component analysis (PCA) is one of the well-known statistical techniques that extracts the rigorous motion information in simulations which is essentially correlated for biological function. PCA works on the variance/covariance matrix which will be constructed based on the trajectories after removal of the rotational and translational movements. The eigenvalues that are identified by diagonalizing the matrix with their projection along with the first two principal components (PC1 and PC2) were calculated using essential dynamics (ED) method by g_covar, g_anaeig modules.

PC1 and PC2 were selected as reaction coordinates for the calculation of free energy, Gα using where k is the Boltzmann constant, T is the simulation temperature. P(qα) corresponds to the probability density function and was constructed using a joint probability distribution of reaction coordinates (PC1 and PC2). Pmax(q) signifies the probability of the most probable state. The protein movements and the dynamic motions of N-CTD and ligands bound to N-CTD throughout the trajectories in the subspace were further identified by Cartesian trajectory coordinates projecting most important eigenvectors from the complete analysis.

Results and discussion

Binding site exploration through molecular docking analysis

To predict the binding site for performing the molecular docking, a ligand-independent binding site search was done on CTD using CASTp Server (Tian et al., 2018). The server identifies a putative binding pocket with a solvent accessible surface area of 907 Å2 and volume of 1166 Å3. This binding pocket of CTD contains residues 259-264 (310 helix), 270, 274 (helix H1), 281-287 (turn and coil region between H1 and H5), H2 (291, 292, 295-296), H3 (301, 304-306), H4 (310-311, 314-318), β2 (329-334), and C-terminal coil region (335-349) (Figure S1, ). We have also used the other two web-servers, PrankWeb (Jendele et al., 2019) and COACH (Yang et al., 2013) for predicting the binding sites. These two servers also predict similar residues in the binding pocket, with minor differences.

Interestingly, we found that most of the binding pocket residues were predicted to be potential dimeric interfacial residues through the InterProSurf analysis (Negi et al., 2007) and the consensus Protein-Protein Interaction Site Predictor (cons-PPISP) server (Chen & Zhou, 2005). The consensus results from these two web-servers revealed that the residues involved in the dimerization of CTD are located in the H4 (Ser310-Ser318), β1 strand (Arg319-Val324), β2 strand (Gly328-Thr334) and in the coil regions adjacent to β2 strand (Gly335-Leu339).

As per the information about the binding site and the dimeric interfaces of the CTD, it is rational to look for the interaction of these target residues with the ligand. For this, we performed the focused molecular docking on the predicted binding site involving H4-β2 and nearby regions in the three-dimensional structure of CTD. For the molecular docking study, we utilized AutoDock, which is the most cited software for docking. Through the molecular docking, we evaluated the binding of five drugs to the N protein of SARS-CoV-2 as shown by Gordon et al. (Gordon et al., 2020). Most of these drugs are currently in clinical trials for the COVID-19, though their antiviral mechanisms are still elusive. Given the importance of N protein oligomerization, we here investigated whether the inhibition of protein dimerization/oligomerization may be a part of the mechanism of action. The five drugs studied in this work are 4E1RCat (ZINC000009311298), Silmitasertib (ZINC000058638454), TMCB (ZINC000058638742), Sapanisertib (ZINC000073069271), and Rapamycin (ZINC000169289388).

All these five drugs were docked against CTD and were ranked based on their binding affinity and the maximum number of interactions with the target residues (Table 1). The proposed binding mode for these drugs is presented in Figure 2. Among the five drugs, 4E1RCat displayed highest binding affinity (-10.95 kcal/mol), followed by rapamycin (-8.91 kcal/mol), silmitasertib (-7.89 kcal/mol), TMCB (-7.05 kcal/mol), and sapanisertib (-6.14 kcal/mol). The binding energy suggests that these drugs bind with the CTD with a strong affinity with KD (KD = expΔG/RT) values in the nanomolar to low micromolar range.

Table 1.

Molecular docking results of Nucleoprotein with compounds.

Compound	Binding affinity (kcal/mol)	Hydrogen bond interactions	Weak interactions
ZINC000009311298	−10.95	Arg277, Lys261 and Tyr333	Thr282, Gln283, Gly284, Asn286, Gly287, Gly275, Phe274, Ala273, Thr271, Val20, Ala264, Thr263, Arg262, Lys26, Gln260, Arg259, Pro258, Thr332, Leu331, Trp330, Ile320, Phe315, Phe315, Phe314, Ser312, Ala311 and Ser310
ZINC000058638454	−7.89	Arg277, Asn285 and Phe286	Phe274, Gly275, Gln283, Gly284, Gly287, Ile357, His356, Leu353,Arg259, Ile292, Lue291, Trp301, Glun303, Ile304, Ala305, Gln306, Phe307, Ala308, Pro309, Ala336, Ile337, Lys338 and Leu339
ZINC000058638742	−7.05	Arg277	Phe274, Gly275, Gln283, Gly284, Asn285, Phe286, Gly287, Val270, Lue291, Ile292, Gly295, Thr296, Trp301, Arg259, Gln260, Lys261. Arg262, Thr263, Ala265 and Thr265
ZINC000073069271	−6.14	Lue331 and Tyr333	Pro258, Arg259, Gln260, Gln281, Thr282, Gln283, Phe314, Arg319, Gly321, Glu323, Thr329, Trp330, Thr332, Thr334 and Gly335
ZINC000169289388	−8.91	Arg277, Glu283, Asn285 and Tyr333	Arg259, Val270, Phe274, Gly278, Pro279, Glu280, Gln281, Thr282, Gly284, Phe286, Gly287, Ile357, Asp358, Lys361, Trp330, Leu331, Thr332, Thr334, Gly335, Ala336, Ser318, Arg319 and Glu323

Figure 2.

Molecular docking interactions and orientations of five anti-N drugs with SARS-CoV-2 CTD (PDB ID: 6WJI). Two-dimensional (2 D) diagrams of ligand-protein interactions for (A) 4E1RCat, (B) Silmitasertib, (C) TMCB, (D) Sapanisertib, and (E) Rapamycin. The hydrogen-bond interactions with residues are represented by a purple arrow directed toward the electron donor. The amino acids are colored according to their properties.

As shown in Figure 2A, 4E1RCat, a translation inhibitor by inhibiting the eIF4E-eIF4G interactions, displayed the highest binding affinity and forms four H-bonds with Lys261, Arg277, and Tyr333 apart from thirteen hydrophobic and thirteen polar interactions (Table 1).

Rapamycin is an mTOR inhibitor and has been shown to reduce MERS infection by ∼60% in vitro (Kindrachuk et al., 2015). It interacted with four target amino acid residues with H-bond i.e. Arg277, Glu283, Asn285, and Tyr333 (Figure 2(B)). Other amino acids involved in interaction comprising nine hydrophobic interactions and ten polar interactions.

Silmitasertib inhibits CK2 and promotes the formation of SGs (Reineke et al., 2017) forms four H-bonds with residues Arg277, Asn285 and Phe286 along with ten hydrophobic and five polar interactions (Figure 2(C)). Similar to silmitasertib, TMCB targets CK2 protein and interacts with CTD through two H-bonds with Arg277 along with six hydrophobic interactions and ten polar interactions (Figure 2(D)). Sapanisertib, like rapamycin targets LARP1 and mTOR inhibitor, interacted with target proteins by forming three H-bonds with Lue331 and Tyr333 along with five hydrophobic and ten polar interactions (Figure 2(E)).

Thus, the binding mode of these five drugs indicates that hydrogen bonding, hydrophobic and polar interactions were the driving forces for binding. The results obtained after docking calculations for all drugs indicate that the most common residues formed H-bonds with CTD are Arg277, Asn285, and Tyr333 (Figure 3).

Figure 3.

The potential binding poses for the drugs. The yellow dotted line represents intermolecular hydrogen bond interactions.

Evaluating the stability and conformational dynamics of CTD-ligand complexes through MD simulation

To further characterize the binding induced structural and conformational changes, all-atom MD simulations for 100 ns were performed to obtain the conformational sampling of the unbound CTD and CTD-drugs complexes.

The selected docking conformations of CTD in complex with drugs were sampled by 100 ns MD simulation, and the dynamic stability of the complex was elucidated by calculating the Cα-RMSD values of the protein as the function of simulation time (Figure 4(A)). The Cα-RMSD values for both free and drug-bound protein changed during the initial 20 ns simulation, and became stable after 30 ns. For CTD-Rapamycin and CTD-TMCB, there were some fluctuations in the beginning and then these complexes gradually reached equilibrium after 30 ns simulation. After 30 ns of simulation, the fluctuation curve of CTD-bound 4E1Rcat, TMCB, and Sapanisertib was lower than the others, suggesting the overall stability of the complex. For CTD-Silmitasertib, fluctuations in RMSD was observed and reached to the maximum value of 0.8 nm at 70 ns of simulation after which it levels back similar to unbounded CTD. Overall, the drug binding induces stability to the CTD, and the decreasing order of the stability is arranged as: TMCB > Sapanisertib > 4E1Rcat > Rapamycin > Silmitasertib.

Figure 4.

Molecular dynamics simulations of drug-protein complexes during 100 ns at 300 K. (A) RMSD of the Cα backbone (B) Radius of gyration, Rg (C) SASA, and (D) RMSF of residues. In all panels the color code is- CTD (black), 4E1RCat (red), Silmitasertib (green), TMCB (blue), Sapanisertib (yellow), and Rapamycin (brown).

The analysis of the one-dimensional probability distribution function (PDF) of Cα-RMSD for CTD-drug complexes indicated that except silmitasertib, all other drugs impart the stability to CTD upon binding (Figure S2A, supplementary material).

For prediction of the compactness of the protein-ligand complexes, the Rg and SASA was calculated for CTD systems (Figure 4(B, C)). The Rg is an effective parameter to evaluate the structural integrity and compactness of the studied systems. It is defined as the mass-weighted root mean square distance of a collection of atoms from their common center of mass. SASA determines the accessibility of protein to the solvent and thus indicates the unfolding of the protein during MD simulation. The time evolution plot of Rg showed that all systems were compact. In the case of CTD-drug complexes, Rapamycin seems to be making the protein more compact with lesser Rg and SASA values, as compared to other drugs. With lesser Rg and SASA values, Rapamycin and 4E1Rcat are observed to be stably bound to CTD in comparison to other drugs. Silmitasertib and TMCB bound CTD complex has increased Rg and SASA values, indicating the potential unfolding effects of silmitasertib and TMCB on CTD. Overall, the drug binding induces a decrease in compactness and increases the flexibility in the order as: Silmitasertib > TMCB > Sapanisertib > 4E1Rcat > Rapamycin.

The PDF of Rg and SASA indicates that the silmitasertib and TMCB have a potentially unfolding effect on CTD and making it unstable. (Figure S2B and S2C, supplementary material).

The RMSF displays the flexibility/mobility of each amino acid residue in the CTD and CTD-drug complexes (Figure 4(D)). For the unbound CTD, high fluctuations were significantly found in the β-sheet residues, 320-330 which is also a predicted binding pocket.

Also, these residues are the most conserved residues of N protein according to multiple sequence alignment results obtained from ESPript 3.0 (Robert & Gouet, 2014) (Figure S3, supplementary material). According to a recent study, the amino acids, 331-333 were the most conserved yet dynamic site of N and involved in the multimerization of N protein (Gupta et al., 2020). The multimeric structure of N revealed that the conserved amino acids Val270, Phe274, Arg277, Asn285, Gly287, Phe286, and Asp288 are clustered and surface exposed, and thus could contribute to protein-protein interactions. Recently, hydrogen-deuterium exchange mass spectroscopy of CTD (PDB ID: 6WZO) results indicated that the residues 323-329 (located within a loop) are relatively exposed, while the residues 330-336 (within a β-strand) is strongly protected from H-D exchange (Troyer & Brady, 2020). Our RMSF results of CTD are consistent with these findings.

Compared to the unbound CTD structure, the CTD-drug complexes showed decreased RMSF values for the C-terminal region, especially for β1 (318-325), β2 (328-334), and α-helix, H7 (345-355) residues, indicating significant stabilization of these regions upon drug binding. Thus, the predicted binding site (329-339) has become one of the most stable regions in the CTD-drug complex system. However, silmitasertib binding to CTD increases the fluctuations of extreme N-and C-termini of the protein, but not to a greater extent and the average RMSF is similar to that of unbounded CTD. Thus, according to RMSF results, the overall fluctuations of the CTD decreased upon drug binding and the maximum decrease were observed in the case of TMCB bound CTD complex.

Hydrogen bonding is one of the most important interactions for stabilizing the protein-ligand complex. Figure 5 shows the number of hydrogen bonds formed in CTD-drug complexes during the last 60 ns trajectory. Sapanisertib, 4E1Rcat, and Silmitasertib bound CTD showed a greater number of average hydrogen bonds, which proves them as good inhibitor molecules. Whereas, hydrogen bonds formation in CTD-TMCB was unstable and absent most of the time during the simulation. For the last 20 ns of simulation, the number of intermolecular H-bonds is greater in the case of CTD-Rapamycin and CTD-Silmitasertib system, indicating the stable binding.

Figure 5.

Number of hydrogen bond interactions formed during MD simulation in the case of CTD-4E1RCat (red), CTD-Silmitasertib (green), CTD-TMCB (blue), CTD-Sapanisertib (yellow), and CTD-Rapamycin (brown).

Secondary structure dynamics of CTD-drug complexes

To further understand the drug binding induced changes in secondary structure, the time evolution of the secondary structure profiles calculated through DSSP are shown in Figure 6.

Figure 6.

Secondary structure changes during MD simulation of CTD-drug complexes. Time evolution of the secondary structure profiles (A) CTD, (B) 4E1RCat, (C) Silmitasertib, (D) TMCB, (E) Sapanisertib, and (F) Rapamycin.

As can be seen, CTD adopts a compact fold rich in helical structures (five α-helices and two 310-helix) along with two β-strands. The detailed changes of the secondary structure of unbounded CTD (Figure 6(A)) showed that H4 was the most unstable secondary structure segment, which began to unfold after ∼20 ns of simulation. The helices, H3 and H4 are present in the dimer interface of a subunit, and also involved in drug binding. A little distortion of β-sheets observed during the initial 10 ns and then 30-60 ns simulation. Notable changes in the secondary structure profile were observed in CTD-drug complexes. In the case of 4E1Rcat bound CTD complex, the interactions between the main-chain of strands β1 and β2 were decreased at the beginning of the simulation and were lost significantly after ∼20 ns (Figure 6(B)). Moreover, the drug binding induces the stability of the H4 structure which remained stable in the presence of 4E1RCat during the simulation. Also, the turn residues between β2 to α5 were lost predominantly to coil.

Silmitasertib binding induces greater changes in the helicity of the protein (Figure 6(C)). The drug-binding overall decreases the helicity of the protein. The helix, H1 lost after 50 ns of simulation, H2 also lost after 40 ns of simulation and then reappear after 45 ns. H3 also partially lost during the last ∼30 ns of simulation while H4 lost after ∼25 ns. However, the helix H6 and the two β-strands remained stable throughout the simulation.

Most significant changes were observed in the case of TMCB bound CTD complex structure (Figure 6(D)), where the residues 25-30 form helix for half of the simulation period and then transformed to coil. However, H4 gets disrupted after ∼15 ns of simulation, and the strands lost after ∼75 ns. For Sapanisertib bound complex (Figure 6(E)), the N-terminal helices, H1 and H2 remained unstable with slight distortion during the first 40 ns and then reappear for the rest of the simulation. Interestingly, the helix H4 which remained lost in CTD and other CTD-drug complexes, appears stable throughout the simulation. The other helices, H3, and H5 also remained stable during the simulation. The β-sheets also remained stable with transient loss during the simulation. Overall, the significant increase in helicity and structural stability was observed in the CTD-Sapanisertib complex. Similar to Sapanisertib, rapamycin binding also induces the stabilization of helices, H3, and H4 along with sheets (Figure 6(F)). The N-terminal helices, H1 and H2 disappear and reappear at the end of the simulation. The formation of other minor secondary structures like turn and bends are more prominent in this complex.

Thus, DSSP information showed that the secondary structure profile of the residues involved in drug binding has some differences among the complexes. For example, binding of the most of the drug (except silmitasertib) to CTD increases the helicity of the residues 54-62 (H4: 310-318 in original PDB). However, the drug binding residues of β1 (319-325) and β2 (328-334) behave differently in the presence of the drugs.

The drug silmitasertib, and rapamycin binding increases the sheet content while the drug 4E1RCat, TMCB, and Sapanisertib binding leads to the loos of sheets and thus decrease the rigidity of the protein complex. Overall, the secondary structure content increases in the case of rapamycin and sapanisertib binding while, 4E1RCat, silmitasertib, and TMCB binding induces a decrease in secondary structure content of CTD (Figure S4, supplementary material).

Evaluation of structural motions and conformational sampling in CTD-drug complexes

Principal component analysis (PCA) was carried out to investigate the important motions during drug binding. The covariance matrix of atomic fluctuations was diagonalized for predicting the eigenvalues. As the first few eigenvectors play a central role in the motions of protein, we have shown here the eigenvalue in decreasing order versus the corresponding eigenvector for CTD systems (Figure 7(A)). As can be seen, the first five principal components account for more than 70% of motions observed for the last 40 ns trajectories of the six CTD systems. The correlated motions showed that unbound CTD showed lesser correlated motions as compared to the CTD-drug complexes. The results showed that the binding of the drug reduced the motion during binding as in the order: TMCB> 4E1RCat > Rapamycin > Sapanisertib > Silmitasertib. Thus, it can be concluded that TMCB and 4E1RCat binding significantly reduced the motion while silmitasertib binding moderately increases the motion.

Figure 7.

Principal component analysis of CTD-drug complexes during 100 ns MD simulation. (A) Plot of eigenvalues vs eigenvector index. First 20 eigenvectors were considered. (B) Projection of the motion of the protein and drug-protein complexes in phase space along the PC1 and PC2. (C) Average Eigen RMSF values for CTD systems were predicted for PC1 and PC2. The color code for the figure is: CTD (black), 4E1RCat (red), Silmitasertib (green), TMCB (blue), Sapanisertib (yellow), and Rapamycin (brown).

Further, the conformational sampling of CTD systems in the essential subspace is shown in Figure 7B which shows the global motions along with the PC1 and PC2 projected by the Cα atom. The figure clearly indicates that unbounded CTD showed fewer stable clusters as compared to CTD-drug complexes, especially in PC1. In the case of the CTD-drug complexes, CTD-TMCB, and CTD-4E1RCat complex showed a more stable cluster than the CTD-Silmitasertib complex. This has also been confirmed with the investigations of the residue displacements along with PC1 and PC2 for the CTD systems (Figure 7(C, D)). Average fluctuation for CTD-drug complexes were lesser along both the eigenvectors compared to the unbounded CTD. Both CTD-4E1RCat and CTD-TMCB bound complex induced fewer motions during binding while CTD-Silmitasertib and CTD-Sapanisertib complex showed a moderate increase in fluctuations and thus could influence the protein motions during binding (Figure 7(C, D)).

Overall, the PCA results indicated that among the five drugs, TMCB and 4E1RCat induces lesser motion and thus more stable complex. The reduced conformational dynamics of the CTD-drug complex may be due to a highly stable and optimized interaction network, which effectively limits the backbone dynamics of both the protein and drugs.

Free energy landscape (FEL)

To demonstrate the effects of drug binding on conformational redistributions, free energy landscapes (FEL) for the six systems were determined as a function of the top two principal components, PC1 and PC2, and are shown in Figure 8. FEL can be used to effectively describe conformational changes induced by ligand binding (Kumar et al., 2018; Prakash et al., 2019). The FEL plots revealed ΔG value 0 to 16.7, 16.1, 13.8, 13.6, 14.7 and 16.3 kJ/mol (Figure 8) for unbound CTD, CTD-4E1RCat, CTD-Silmitasertib, CTD-TMCB, CTD-Sapanisertib, and CTD-Rapamycin complex, respectively.

Figure 8.

The free-energy landscape (FEL) of the simulated CTD and CTD-drug complexes based on the principal component analysis. (A) CTD, (B) CTD-4E1RCat, (C) CTD-Silmitasertib, (D) CTD-TMCB, (E) CTD-Sapanisertib, and (F) CTD-Rapamycin. The color bar represents the free energy value according to kcal mol−1. Dark blue spots indicate the energy minima and energetically favored protein conformations, and yellow spots indicate the unfavorable high-energy conformations.

In reference to Figure 8, we could visualize that unbound CTD explores a single large conformational space along with a single global energy minima (Figure 8(A)). The FEL of CTD -4E1RCat complex was similar to that of unbound CTD, suggesting minimal changes in conformational rearrangements (Figure 8(B)). However, for CTD bound to silmitasertib or TMCB, we observed that the protein explores a huge conformational space and displays several small local minima along with a reversed population shift relative to the unbound CTD (Figure 8(C)). This indicates that these drugs induce diverse ensembles of flexible conformations during the 100 ns simulation. Similarly, CTD-Sapanisertib and CTD-Rapamycin bound complex explores different intermediate conformations linked by low lying energy barriers (Figure 8(D)). Thus, compared to unbound CTD free energy surfaces, there is a considerable population shift in all the CTD-drug complex systems, suggesting global conformational switching to a number of metastable states. These results lead us to the conclusion that silmitasertib and TMCB strongly modulates CTD protein conformations, and thus could regulate the function of the protein.

Density distribution analysis

The atomic density distribution was further investigated to observe the changes in the atomic orientation and distribution due to drug binding. We plotted the partial density distribution for CTD systems calculated using the densmap analysis available in the Gromacs package (Figure S5 a–f, supplementary material). The partial density area of CTD was observed to be stable with a value of 2.29 nm−3, whereas high density was noticed for drugs 4E1Rcat (2.61 nm−3), silmitasertib (2.88 nm−3), and TMCB (2.47 nm−3). However, the drugs sapanisertib (1.98 nm−3) and rapamycin (2.05 nm−3) showed a lesser density area as compared to CTD. The density maps thus indicated that drug binding induced conformation changes to CTD which is much more evident in the case of silmitasertib and sapanisertib. This result was also in accordance with the DSSP results, further indicating a major conformational loss in the presence of these drugs (Figure 6).

Dynamic cross-correlation matrices analysis

To monitor the correlated conformational motions of CTD-drug complexes, the dynamic cross-correlation matrices (DCCM) analysis was further analyzed (Figure 9). Here, highly positive regions (red) signify a strong correlation in the movement of residues in the same direction, while the negative regions (blue) indicate strong anti-correlated motion of the residues (residues move in the opposite direction). The color depth of diagonal could signify the movement degree of atoms. As shown in Figure 9, among the five CTD-drug complex systems, silmitasertib bound CTD has relatively stronger correlated motions followed by CTD-4E1RCat complex. Also, the more correlated and anti-correlated motions between residues are observed in CTD-Silmitasertib complex, followed by CTD-4E1RCat, CTD-Rapamycin, CTD-Sapanisertib, and CTD-TMCB complex. These noticeable differences suggested the presence of more and stronger cross-correlation motions in CTD-Silmitasertib and CTD-4E1RCat, indicating more strong interaction and better stability of these complexes. Also, the lack of correlation in the motion of the CTD-TMCB and CTD-Sapanisertib complex hinted toward the loss of contact among the residues, which lead to the destabilization of the complex.

Figure 9.

Dynamic cross-correlation map (DCCM) of the Cα atoms around their mean positions during MD simulation. (A) CTD, (B) CTD-4E1RCat, (C) CTD-Silmitasertib, (D) CTD-TMCB, (E) CTD-Sapanisertib, and (F) CTD-Rapamycin. The degrees of the correlation motions and anti-correlation motions are represented in blue and red, respectively.

Tertiary contact map

Figure 10 shows the pairwise contact maps for the CTD and CTD-drug complexes This enables us to obtain insights into the details of residual contacts and thus structural changes due to drug binding. The contact map analysis indicates that binding of drugs leads to the loss of tertiary contacts especially in the N-terminal region of the protein. The drug 4E1RCat (Figure 10(B)) induces maximum changes in residual contacts as the drug-binding decreases the contacts of N-terminal residues (residues 1-10) to H1 (24-30), H3 (45-50) and H4 (54-62). The N-terminal helix H1 also lost its contacts with H3 and H4. The contacts between C-terminal β2 (72-80) and H6 (104-109) also gets disrupted. However, these drug binding leads to the gain of some long-range contacts between N-terminal 310 helix residues to C-terminal β-sheets. Silmitasertib binding also weakens the short-range interactions of H1(24-30) to N-terminal residues (1-10), H3 (45-50), and H4 (54-62) along with minor loss of long-range contacts between N-and C-termini (Figure 10(C)). The drug TMCB binding also disrupts the contacts similar to silmitasertib, however it gains some non-native contacts between N-terminal 310-helix and C-terminal β2 (Figure 10(D)). The change in contact map in the case of sapanisertib binding is similar to that of 4E1RCat binding with the maximum disruption of H3 and H4 helices (Figure 10(E)). Moreover, rapamycin binding doesn't show significant changes in the contact map and behave very much similar to unbounded CTD (Figure 10(F)). Overall, the drug binding results in the loss of tertiary contacts of the protein and thus implies the structural loss of the protein which is more significant in 4E1RCat and sapanisertib.

Figure 10.

Tertiary contact map of CTD changes on drug binding. Contact map of (A) CTD, (B) CTD-4E1RCat, (C) CTD-Silmitasertib, (D) CTD-TMCB, (E) CTD-Sapanisertib, and (F) CTD-Rapamycin. The loss of residual contacts due to drug binding is represented in the red circle, while the gain of contacts is shown in the black circle.

Conclusions

The involvement of N protein in many crucial functions of the viral life cycle makes it a candidate drug-target in COVID-19. However, the molecular mechanism for drug binding and inhibition of antivirals against newly emerged novel SARS-CoV-2 N protein remain largely elusive. Understanding these aspects should facilitate the discovery of compounds that specifically inhibits CoV genome replication and assembly. Here, we analyzed the binding mechanism of anti-N drugs against SARS-CoV-2 CTD at the molecular level with the help of structure-based computational methods.

Our study indicated that these five drugs bind CTD with strong affinity values in the nanomolar to low micromolar range. It can be observed that the drug-binding residues also form the dimeric interface and become one of the most stable regions in the CTD-drug complex system. The changes in structural and conformational dynamics of CTD upon drug binding are different among the drugs. Overall, the MD simulation parameters like RMSD, RMSF, Rg, hydrogen bonds analysis, PCA, FEL, and DCCM analysis indicated that 4E1RCat and TMCB complexes were more stable as compared to silmitasertib and sapanisertib. In conclusion, the study elucidated the detailed interaction mechanism of drugs binding to CTD and the associated structural and dynamical changes upon drug-binding. These results will significantly enhance our understanding of the working mode of these drugs at the molecular and structural level, and will contribute to the future rational drug design for COVID-19.

Click here for additional data file.