|
article-title
|
Dual inhibition of SARS-CoV-2 spike and main protease through a repurposed drug, rutin
|
|
alt-title
|
A. Kumari et al.
|
|
alt-title
|
Journal of Biomolecular Structure and Dynamics
|
|
abstract
|
Abstract
The global health emergency caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has led to alarming numbers of fatalities across the world. So far the researchers worldwide have not been able to discover a breakthrough in the form of a potent drug or an effective vaccine. Therefore, it is imperative to discover drugs to curb the ongoing menace. In silico approaches using FDA approved drugs can expedite the drug discovery process by providing leads that can be pursued. In this report, two drug targets, namely the spike protein and main protease, belonging to structural and non-structural class of proteins respectively, were utilized to carry out drug repurposing based screening. The exposed nature of the spike protein on the viral surface along with its instrumental role in host infection and the involvement of main protease in processing of polyproteins along with no human homologue make these proteins attractive drug targets. Interestingly, the screening identified a common high efficiency binding molecule named rutin. Further, molecular dynamics simulations in explicit solvent affirmed the stable and sturdy binding of rutin with these proteins. The decreased Rg value (4 nm for spike-rutin and 2.23 nm for main protease-rutin) and stagnant SASA analysis (485 nm/S2/N in spike-rutin and 152 nm/S2/N in main protease-rutin) for protein surface and its orientation in the exposed and buried regions suggests a strong binding interaction of the drug. Further, cluster analysis and secondary structure analysis of complex trajectories validated the conformational changes due to binding of rutin.
|
|
title
|
Abstract
|
|
p
|
The global health emergency caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has led to alarming numbers of fatalities across the world. So far the researchers worldwide have not been able to discover a breakthrough in the form of a potent drug or an effective vaccine. Therefore, it is imperative to discover drugs to curb the ongoing menace. In silico approaches using FDA approved drugs can expedite the drug discovery process by providing leads that can be pursued. In this report, two drug targets, namely the spike protein and main protease, belonging to structural and non-structural class of proteins respectively, were utilized to carry out drug repurposing based screening. The exposed nature of the spike protein on the viral surface along with its instrumental role in host infection and the involvement of main protease in processing of polyproteins along with no human homologue make these proteins attractive drug targets. Interestingly, the screening identified a common high efficiency binding molecule named rutin. Further, molecular dynamics simulations in explicit solvent affirmed the stable and sturdy binding of rutin with these proteins. The decreased Rg value (4 nm for spike-rutin and 2.23 nm for main protease-rutin) and stagnant SASA analysis (485 nm/S2/N in spike-rutin and 152 nm/S2/N in main protease-rutin) for protein surface and its orientation in the exposed and buried regions suggests a strong binding interaction of the drug. Further, cluster analysis and secondary structure analysis of complex trajectories validated the conformational changes due to binding of rutin.
|
|
abstract
|
Graphical Abstract
Communicated by Ramaswamy H. Sarma
|
|
title
|
Graphical Abstract
|
|
p
|
Communicated by Ramaswamy H. Sarma
|
|
body
|
1. Introduction
Pandemics involving pathogenic human coronaviruses have wreaked havoc on the human populace. During the severe acute respiratory syndrome coronavirus (SARS-CoV) epidemic in 2002–2003, a high fatality rate of 10% was observed in the approximately 8,000 individuals infected (Marra et al., 2003; Peiris et al., 2003; Rota et al., 2003; SARS Working Group, 2003). In 2012, the middle east respiratory syndrome coronavirus (MERS-CoV) had a much higher fatality rate of 36% and infected more than 1,700 individuals (de Groot et al., 2013; Zaki et al., 2012). However, the recent SARS-CoV outbreak of 2019, also referred to as SARS-CoV-2 or COVID-19 (Chen et al., 2020), began in Wuhan in China and spread rapidly due to close human-to-human proximity (Li et al., 2020). The SARS-CoV-2 has caused significant devastation globally due to the sheer number of people infected. According to statistics published by the World Health Organization (WHO), there were 4,334,451 known cases of Covid-19 worldwide, out of which 19,737,794 cases alone are reported from America, while Europe attributes 9,664,042 cases, causing death casualties of over 1,157,509 globally as of Oct 27, 2020.
Similar to SARS-CoV and MERS-CoV, SARS-CoV-2 (exhibiting 80% and 50% homology, respectively) (Kim et al., 2020; Zhu et al., 2020) belongs to the genera beta coronavirus (family Coronaviridae in the order Nidovirales) (Enjuanes et al., 2006; Perlman & Netland, 2009), which is known to infect mammals (Li, 2016) and has manifested an illustrious capability of cross-species transmission, including humans (Menachery et al., 2017). The SARS-CoV-2 is an enveloped virus characterized by a positive-sense, single-stranded RNA genome of almost 30 kb, surrounded by a helical capsid comprised of nucleo capsid protein (N). The viral envelope is associated with three primary structural proteins viz. membrane proteins (M) and envelope proteins (E), which perform virus assembly; and spike proteins (S), which facilitate the viral attachment and thus, the virus entrance into the host cells. The large protrusions formed by spikes on the surface of the virus give it the appearance of having crowns, which, in Latin, translates to corona (Li, 2016). These shorter, sgRNA-encoded structural proteins and several accessory proteins are known to be conserved (Kim et al., 2020). The large ectodomain of the spike consists of receptor binding S1 and membrane fusion S2 subunits. For many CoVs, the S1 and S2 domains remain non-covalently linked. In β coronaviruses, the cleavage between the S1 and S2 regions is not obligatory. However, the host proteases have been observed to cause cleavage within the fusion domain (S2), which leads to irreversible conformational changes, activating the protein for membrane fusion (Zhou et al., 2020). The S2 subunit contains two regions with 4, 3 hydrophobic heptad repeats (HR), HR1 and HR2, which are conserved in sequence and position. The HR2 is located adjacent to the trans-membrane domain, while the HR1 occurs ∼170 residues upstream of HR2. The HR region has been observed to be a common motif in several viral fusion proteins. Structurally, the HR1 domain exists as a homo-trimeric coil, packing the HR2 domain (in an anti-parallel manner) in its hydrophobic grooves, thereby bringing the N-terminal fusion peptide closer to the trans-membrane anchor, further facilitating fusion due to proximity (Bosch et al., 2003). Various host receptors viz. angiotensin-converting enzyme 2 (ACE2), amino peptidase N (APN), dipeptidyl peptidase 4 (DPP4), carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1), and sugar recognize S1 domain of the spike protein to render virus entry into cells (Li, 2016). Several reports have been published recently discovering inhibitors against the target via a computational methodologies, which involve the usage of FDA approved compounds and natural compounds (Wei et al., 2020), spice molecules (Rout et al., 2020), small-molecule compounds of ZINC Drug Database (Kadioglu, 2020; Wu et al., 2020) along with traditional Chinese medicine and natural products and derivatives (Wu et al., 2020) and medicinal compounds (Salman et al., 2020).
Apart from structural proteins which play a crucial role in the lifecycle of the virus, remarkable number of functional proteins categorized as non-structural proteins (Nsp’s), involved in viral replication, transcription, translation and protein modifications, and host infection are equally important (Wu et al., 2020). Among many such proteins, a well characterized distinguished drug target for SARS-Cov-2 is the 33.8 kDa, 3 C-like main protease, also called as the main protease or Nsp5 (3CLpro or Mpro or Nsp5) (Wu et al., 2020; Zhang et al., 2020). Together with the papain like proteases (Plpro), the enzyme is responsible for the processing of polyproteins produced by the viral RNA (Zhang et al., 2020). The Mpro along with Plpro cleaves the 790 kDa polyprotein1ab to generate 15 Nsp’s. The Mpro processes the polyprotein1ab at 11 sites to produce mature Nsp4-Nsp16 (Wu et al., 2020). Recently, the crystal structure of Mpro has also been reported at a resolution of 1.75 Å. The Mpro structure consists of three domains, domains I and II comprising of six-stranded anti-parallel β barrels spanning from residues 10 to 99 and 100 to 182, respectively forming the substrate-binding site between them and the domain III which is responsible for modulating the Mpro dimerization is a cluster of five helices ranging from residues 198 to 303 (Zhang et al., 2020). Due to the availability of the crystal structure and the absence of any human homologue, the Mpro is a lucrative drug target to work upon for the discovery of novel antiviral agents (Jin et al., 2020; Wu et al., 2020; Zhang et al., 2020). Recently, numerous studies have emerged against this target, discovering novel inhibitors by utilizing computational approaches which include the use of natural molecules (Aanouz et al., 2020; Das et al., 2020), commercially available drugs and zinc library (Das et al., 2020; Elmezayen, 2020; Ton et al., 2020), antiviral compounds (Khan, 2020; Kumar et al., 2020; Muralidharan, 2020), peptide molecules (Pant, 2020), generative chemistry approaches for drug design (Alex, 2020), combinatorial strategy of using anti-virals, natural products, anti-fungals, anti-protozoals and anti-nematodes (Das et al., 2020), and spice molecules (Rout et al., 2020).
Many pharmacological agents such as antiviral agent’s viz. remdesvir, ribavirin, favipiravir, chloroquine, hydroxychloroquine, oseltamivir and umifenovir; immunomodulatory agents including toclizumab, inteferons; adjunctive agents like azithromycin, corticosteroids have been repurposed for treating COVID-19. However, so far no therapy has been (Lam et al., 2020) recognized as an effective treatment against the deadly disease. It has, therefore, become the need of the hour to adopt computational approaches for drug discovery and repurpose other existing marketed drugs for its treatment. Drug repurposing via high-throughput screening of various databases is a cost-effective and time-saving approach. The present report pertains to the above-described drug discovery effort using a repurposing screening strategy to target the SARS-CoV-2 spike protein and main protease. Together, these proteins play a key role in binding to host cell receptors and causing fusion to render viral entry, determining the host range and tissue tropism, eliciting an immune response by the host, viral replication, transcription and protein processing (Li, 2016).
2. Material and methods
2.1. Preparation of protein and ligand structures
The spike protein and main protease regulatory enzyme of SARS-CoV-2 have been revealed to play a vital role in COVID-19 infection in the host (Ke et al., 2020; Xue et al., 2008). In order to delineate the drug treatment, the deduced three dimensional structures of target SAR-CoV-2 spike protein in open state (PDB ID: 6VYB) and main protease (PDB ID: 6LU7) were retrieved from RCSB Protein Data Bank (Burley et al., 2019). The crystal structures of target proteins were imported to Schrödinger Maestro (Moore, 2015) for molecular docking analysis. Pre-processing and preparation of target proteins were performed using Schrödinger Protein Preparation Wizard (Schrödinger, 2011) by adding the missing hydrogen atoms, correcting the bond orders, capping the protein termini, creating disulfide bonds and removing the water molecules. Also, the H-bonds were optimized and further the target proteins were equilibrated using OPLS 2005 force field (Robertson et al., 2015).
The library of ligands was downloaded from DrugBank (Wishart et al., 2018) and the coordinate files were prepared and optimized using LigPrep module (Release, 2017). The energy minimization of ligand structures was carried using OPLS 2005 force field and ionization states were generated at pH 7.0 ± 2.0 using Epik. This generated an output file with various stereoisomers and tautomeric conformers of each ligand producing chemical and structural diversity.
2.2. Determination of binding pocket and generation of receptor grid
The active sites enclosed in the binding pocket of spike protein and main protease were predicted using Computed Atlas of Surface Topography of proteins (CASTp) webserver (Tian et al., 2018) which calculates geometric and topological measurements of protein structures. It measures the volume and surface area of the binding pocket and also predicts the atoms involved in formation of this binding cavity. Although the active sites of spike protein and the main protease are known but still research is undergoing in which it is predicted that new amino acids might be involved in the inhibition of virus activity. In order to undertake all the possibilities, binding pocket analysis was performed through CastP. Further for the ease of molecular docking, a receptor grid was generated around the active sites of target proteins using Receptor Grid generation module. The receptor grid file was created with a grid box around the centroid of selected active residues.
2.3. Protein ligand interactions and molecular docking
Protein ligand docking was conducted using GLIDE module of Schrödinger (Release, 2018) which performs all the docking calculations using OPLS 2005 force field. A total of 2652 FDA approved drugs were screened against spike protein and main protease of SAR-CoV-2 through molecular docking approaches. The virtual screening of compounds against the target protein was conducted using two docking methodologies, initially through a rapid screening of large set of compounds using high throughput virtual screening (HTVS) method followed by further screening a subset of top scoring HTVS compounds with a more precise and accurate extra precise (XP) docking method. Protein ligand complexes were ranked using GlideScore function to predict the binding efficacy of ligands with the protein. The docked complex with the lowest docking score was selected for further molecular dynamic simulation analysis to elucidate the inhibition mechanism against SARS-CoV-2 infection.
2.4. Molecular dynamics trajectory analysis of protein and their docked complex
The SARS-CoV-2 spike protein and main protease were used separately for molecular dynamics (MD) simulations alone and in complex with the repurposed drug, rutin. All the four MD simulations were executed through GROMACS v5.0 under the force field GROMOS96 54a7 having water model SPC216 along with the time step of 1 fs for 100 ns (Abraham et al., 2015; Darden et al., 1993). Varied sizes of the simulation box were created for each MD simulation event, which were further loaded with about respective amount of water molecules using the SPC model. The total charge on spike protein, main protease and the two in complex with rutin were neutralized by adding –6, –4, –6, and –4 charges, respectively, and were incorporated into the simulation system by compensating the water molecules in the arbitrary locations inside the simulation box. The NPT ensembles, along with periodic boundary conditions, were utilized to carry out MD simulations. A cut-off of about 12 Å was used in order to manage the Vander Waals forces. The Particle Mesh Ewald model manifesting a cut-off of 14 Å was further utilized to calculate the electrostatic interactions (Darden et al., 1993). The spike protein, main protease and the two in complex with rutin were solvated through a slab of about 10 Å in every direction. The neighbor list was updated to a frequency of 10 ps.
The MD simulations were achieved for each system employing the four major steps. The first step deals with the energy minimization of the entire system utilizing the integrator of steepest descent in continuation with second integrator of conjugate gradients algorithms. The second step involves the minimization and molecular dynamics of NVT and NPT ensembles for 500 ps and 1000 ps, respectively allowing the solvents and ions to evolve by keeping the same starting configuration for the structures. In the third step the systems were heated using a lower temperature coupling (τ = 0.1 ps) along with pressure coupling (τ = 0.5 ps) to attain equilibrium at 300 K and 1 atm of temperature and pressure. In the equilibration phase, the thermostat and barostat were evaluated through the Berendsen algorithm (Berendsen et al., 1984). The hydrogen-containing bond lengths were constrained with the help of the LINCS algorithm (Hess et al., 1997). Finally, the last step also called as the production step was carried out, where the MD simulation for 100 ns at 300 K temperature with 2 fs of time step were performed for both systems, and the last final structures were achieved. The Maxwell Boltzmann distribution was utilized in order to reassign the velocities at every step. The Nose Hoover thermostat and Parrinello Rahman barostat were the respective thermostat and barostat for the final MD or production run (Berendsen et al., 1984).
Various analyses were performed with the help of inbuilt analysis commands of GROMACS. The root mean square deviation (RMSD) is a magnitude of the dimensional disparity among the two stagnant structures, and RMSD calculation is achieved depending upon the native structure and each consecutive trajectory frames in the simulation. In addition, root mean square fluctuation (RMSF) profile measures the affability of every protein residue depending on the fluctuation about an average location within all MD simulations (Knapp et al., 2011). Therefore, RMSD and RMSF of each simulation system were determined to examine the stability and residual fluctuations. Further, the radius of gyration (Rg) analysis was performed to evaluate the compactness of both the simulation systems separately. Also, the hydrogen bond analysis was performed to check the neighboring interactions with both simulation systems separately, including the hydrophobic interactions with the help of the LigPlot tool for both spike protein and main protease complexes with rutin before and after simulation.
Additionally, the solvent accessibility surface area (SASA) was also computed to examine the solvent attributable areas of all simulation system. The cluster analysis having a cut-off value of 0.25 and 0.2 nm for spike protein and main protease respectively, depending upon the RMSD profile were utilized to demonstrate the conformations found utmost intermittently throughout the trajectory. Here, all the structures having RMSD values of below 0.25 nm for all components within a cluster are incorporated to the initial cluster. It is rare that a molecule having a higher value for RMSD than 0.25 nm from other cluster supposedly would be treated as a structure. The secondary structure analysis was also performed using the DSSP program (Martin et al., 2005). The visualization of protein nature during the entire simulation was accomplished by using Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) and UCSF Chimera (Pettersen et al., 2004).
3. Results and discussion
SARS-CoV-2 is an enveloped non-segmented positively stranded RNA infectious virus that mainly affects enteric and respiratory systems (Graham et al., 2013). Its infection and amelioration in the host are majorly governed by several structural proteins (Joshi et al., 2020; Prajapat et al., 2020). Amidst numerous structural proteins, spike protein and main protease of SARS-CoV-2 are stated to show essential impacts on viral replication through proteolytic machinery and involvement in transcription, translation, and amplification of viral proteins (Paules et al., 2020). Spike protein plays a critical role in binding to host cell receptor and is thought to represent a key determinant of the host range restriction (de Wilde, 2017). Main protease on the other hand exclusively cleaves polypeptide sequences after the glutamine residue and to the best of our knowledge is known as an ideal drug target due to the absence of human host cell protease showing this substrate specificity (Ullrich & Nitsche, 2020). The 3 D structures of spike protein and main protease were preprocessed, optimized and minimized to obtain the refined structures for drug repurposing.
3.1. Binding pocket analysis
The binding pocket volume and surface area were determined through CASTp webserver. Binding pocket was forecasted at the surface as well as in the interior of proteins. Binding pocket volume of spike protein and main protease was 22,908 and 319 (SA), respectively which signifies an optimum space for ligand binding. The binding residues as predicted in the binding cavity for spike protein were Y38-S45, V47-L54, P85-D88,T108,K195-Y200,K202-Y204,P225-D228,I233-R237,Q271-T274,C291-A292,E298,C301-K304,Y313-N317,R319-V327,D364-S366,Y369-N370,Y380-T393,P412,D427-T430,F515-H519,A522,P527-T531,N540-T549,Q564-D568,A570-T573,V576-R577,P579,I587,P589-F592,G594,M740-Y741,G744,F855-N856,Q954,Q957-A958,N960-D979,R983,V991,D994-R995,T998-R1000,Q1002-S1003,T1006-Y1007,Q1010, and R1014 around which the receptor grid was formed. The receptor grid was also created around the catalytic region (T24-L27, H41, C44-S46, M49, P52, Y54, F140-C145, H163-P168, H172, D187-T190, and Q192) of main protease.
3.2. Probable drugs against spike protein and main protease of SARS-CoV-2
The molecular docking studies were performed to determine the ligand interactions and their binding affinity with SARS-CoV-2 spike protein and main protease. The prepared ligands structures were docked against the binding sites of the target proteins. The HTVS docking approach filtered out a large number of compounds on the basis of their binding interactions within the binding sites of the target proteins. The compounds displaying a HTVS glide score of more than −6.0 kJ/mol were employed for XP docking. The docking simulations evaluated a high Glide score for the compounds with best binding affinity and their interactions within the binding pocket of the binding site residues. The docking analysis revealed a common compound, rutin (DrugBank ID: DB01698) which binds strongly to both the major targets of SARS-CoV-2. Our study here presents rutin exhibiting high binding efficiency against spike protein and main protease with an XP Glide score of −8.367 kcal/mol and −11.553 kcal/mol respectively. Rutin (3,3′,4′,5,7-pentahydroxyflavone-3-rhamnoglucoside) belongs to the flavonol class of compound, generously present in plants including buckwheat, apple, passion flower and tea. It is well known to exhibit many bioactivities, such as anti-viral, anti-bacterial, ant-larvicidal, cytoprotective, anticarcinogenic, antioxidant, vasoprotective, cardioprotective and neuroprotective activities (Saluja & Ganeshpurkar, 2016). The list of best ten candidates derived from the docking calculations with their corresponding docking score against spike protein and main protease are shown in supporting information (supporting material Tables S2 and S3). To determine the interaction pattern between rutin and target proteins, their docked complexes were visualized using Ligplot (Figure 1). The interaction pattern between rutin and spike protein revealed the formation of five hydrogen bonds and four hydrophobic interactions at the binding site, thereby contributing to the stability of the complex. The residues involved in the formation of hydrogen bonds between rutin and spike protein were F970, N969, H49, Q52, and T274. While the hydrophobic contact between the two involved residues T51, S50, S967 and S968. The interactions between rutin and main protease was stabilized by residues E166, T190, Y54, D187, T26, N142 forming hydrogen bond and various hydrophobic contacts with residues P168, L167, H41, G143, T25, L27, M49, M165, H164, R188 and Q189 (Tables 1 and 2). A similar study exploiting a blind molecular docking approach utilizing the Swiss Dock server also highlights main protease inhibition by rutin where by common amino acids viz. E166, T190, N142, H41 are found to be interacting with the inhibitor. In our study however, the inhibitor displayed slightly better binding efficiency than their report (Das et al., 2020). From the contact analyses of both proteins, it can be attributed that rutin has high affinity and wide molecular contacts for both spike protein and main protease. Also, it was further scrutinized by molecular dynamics simulations to attain insights towards the inhibitory aspects and efficacy to combat spike protein and main protease drug targets of SARS-CoV-2.
Figure 1. Protein-ligand interactions diagram of spike protein and main protease in complex with rutin. 2D Ligplot representation showing hydrogen bond and hydrophobic interactions, and 3D structure representation showing rutin in the binding pocket of spike protein (A and C) and main protease (B and D).
Table 1. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between spike protein with Rutin.
Spike protein + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
His 49.A NE2–Rutin.het O15 3.28 HNE2-H--O15 Ser50, Thr51, Ser967, Ser968
GLN 52.A OE1–Rutin.het O14 2.89 HOE1-H--O14
THR 274.A OG1–Rutin.het O14 3.30 HOG1-H--O14
ASN 969.A ND2–Rutin.het O10 3.17 HND2-H--O10
PHE 970.A N–Rutin.het O9 3.10 HN-H--O9
Table 2. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between main protease with Rutin.
Main protease + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
THR 26.A O–Rutin.het O15 3.05 HO-H--O15 Thr25, Leu27, His41, Met49, Gly143, His164,
THR 26.A O–Rutin.het O16 3.02 HO-H--O16 Met165, Leu167, Pro168, Arg188,
TYR 54.A OH–Rutin.het O14 2.91 HOH-H--O14 Gln189
ASN 142.A OD1–Rutin.het O5 3.01 HOD1-H--O5
ASN 142.A OD1–Rutin.het O7 2.91 HOD1-H--O7
GLU 166.A O–Rutin.het O8 3.02 HO-H--O8
ASP 187.A O–Rutin.het O14 3.31 HO-H--O14
THR 190.A O–Rutin.het O10 2.87 HO-H--O10
3.3. Molecular dynamics simulation of spike protein and main protease in absence and presence of rutin
The binding affinity of rutin against both spike and main protease were investigated through MD simulations using GROMACS. Separate MD simulations for spike protein and main protease with and without rutin for 100 ns were assessed through trajectory analysis. For the course of 100 ns MD simulation, the stable trajectory was observed and the representative structures were obtained. The deviation of the backbone atoms for simulated structures relative to the starting structures used as a reference was evaluated through RMSD for the backbone atoms (Figures 2 and 3) and C-alpha atoms (supporting material Figure S1). A steep magnitude RMSD variation during the entire simulation can be an implication of a malleable and free instinctive protein or the alteration of the force field. Depending upon the outcomes of the RMSD evaluation, Figure 2(A) represents that the RMSD fluctuation stabilize at about 50 ns MD simulations for both spike protein and spike protein in presence of rutin, and the simulation time was acceptable. In the time from 50–100 ns, the RMSD for spike protein and spike-rutin complex have approximate values about 1.50–1.65 nm and 1.20–1.26 nm, respectively. Similarly, the main protease in absence and presence of rutin for time from 20–70 ns and 20–100 ns showed RMSD values of about 0.25–0.30 nm and 0.20–0.25 nm (Figure 3(A)). Interestingly, RMSD for both the complex systems were more stable than the native state of target receptor spike protein and main protease (Figures 2(A) and 3(A)). The RMSD of C-alpha atoms were also attained and was found to be similar to the RMSD of backbone atoms (supporting material Figure S1). The average RMSD for backbone and C-alpha atoms of both the protein and their complex with rutin are shown in supporting material Table S1.
Figure 2. Representation of the MD analysis plots for spike protein (purple color) and spike-rutin complex (red color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg), and (D) solvent accessible surface area (SASA) analysis.
Figure 3. Representation of the MD analysis plots for main protease (blue color) and main protease-rutin complex (green color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg) and (D) solvent accessible surface area (SASA) analysis. The integral extent of the residual dynamic criterion is achieved by assessing the variations arising from shifts of each of the protein residues that majorly feature the most flexible chain frames. Hence, we validated the residual fluctuations by calculating the mean fluctuation for stable trajectories of each simulation. The RMSF evaluation of all the protein residues were achieved in order to check the residues that may have tend towards an enhancement in the RMSD results. Compelling fluctuations existed in terminal residues and in various loops along β sheets and coils of the spike protein (residues 852-1000) to about 0.47 nm and spike-rutin (residues 62–90, 400–520 and 675–730) to about 0.65 nm (Figure 2(B)). Similarly, in case of main protease and itself in complex with rutin showed fluctuations in the terminal residues and in some regions of β sheet with some coils and loops in main protease (residues 136-148) to about 0.37 nm and main protease-rutin (residues 39-52, 150-156) to about 0.37 nm, respectively (Figure 3(B)). Also, we noticed that residues 270 to 350 in spike protein and 160 to 190 in main protease showed comparatively reduced discrepancy in RMSF values, which were the corresponding binding segment of both the complexes after MD simulation. Interestingly, the binding region for spike-rutin and main protease-rutin before MD simulation were found to be the same after simulation, suggesting as a strong binding interaction of the rutin towards spike protein and main protease and indicating stability of both the complexes.
The radius of gyration analysis was evaluated to determine the change in compactness of protein systems used throughout the MD simulations. The Rg plots for spike, main protease and its complex with rutin show slight fluctuations at the initial frame and attain compactness after 30 ns with Rg score of 3.85 nm and 4 nm (spike and spike-rutin) and after 20 ns with Rg score of 2.27 nm and 2.23 nm (main protease and main protease-rutin), respectively (Figures 2(C) and 3(C)). When compared with spike protein, Rg value for spike-rutin is stabilized and remains constant, suggesting strong binding interaction of the inhibitor, and the same is observed in case of main protease and rutin. Similar observations were determined through SASA analysis representing the solvent defined protein surface and its orientation through folding, making the alterations in the exposed and buried regions of the surface area of proteins. The SASA values for all the simulation systems were about 460 nm/S2/N (spike), 485 nm/S2/N (spike-rutin), 154 nm/S2/N (main-protease), and 152 nm/S2/N (main protease-rutin) (Figures 2(D) and 3(D)). Here, spike-rutin and main protease-rutin solvation profile shows a convincing SASA value suggesting a stable structure and strong binding interaction with the rutin.
Further, the cluster analysis with a RMSD based cut-off value of 0.25 nm demonstrated the development of 13,28,3 and 2 distinctive clusters for spike, spike-rutin, main protease, and main protease-rutin simulation systems respectively (Figure 4(A–D)). The most dominant cluster attained after 100 ns of MD simulation for all four systems are shown in Figure 4(E–H). Also, the secondary structure analyses of the stable trajectory for both the simulation systems were performed using the DSSP tool of GROMACS. Both, spike and main protease were formed mainly of conserved β-sheet along with some connecting loops and α-helix with small coil regions as secondary structure elements infused with various small segments of bend, turn, and β-bridge (Figure 5). Both cluster analysis and secondary structure analysis reveals the conformational changes before and after simulations for spike and main protease with and without rutin structures. A noticeable observation through these analyses supports the rationale that spike protein and main protease binds with rutin firmly and attain similar interactions after MD simulations also.
Figure 4. Representation of the MD analysis plots for the RMSD matrix cluster formation showing number of clusters in (A) spike protein, (B) spike-rutin complex, (C) main protease, and (D) main protease-rutin complex. The most dominant cluster conformation after MD cluster analysis of (E) spike protein, (F) spike-rutin complex, (G) main protease and (H) main protease-rutin complex.
Figure 5. The secondary structure analysis of (A) spike protein, (B) spike-rutin complex, (C) main protease and (D) main protease-rutin complex. Further, hydrogen bond analysis revealed a higher number of hydrogen bonds with the surrounding water molecules in case of spike protein (2300 hydrogen bonds) than in complex with rutin (approximately, 2200 hydrogen bonds). Similarly, in case of main protease, 700 hydrogen bonds were formed while 560 hydrogen bonds in case of main protease complex with rutin (Figure 6(A,B)). In addition, the hydrogen bond landscape assessed against the inhibitor revealed the dynamic equilibration of the complex trajectories with a high number of hydrogen bonds, as shown in Figure 6(C,D). The average number of hydrogen bonds within spike-rutin and main protease-rutin are 3.72 and 2.29, respectively. The consistent numbers of hydrogen bonds were observed which contributed significantly to the proximal binding of the potent drug rutin with the spike and main protease receptor. The superimposed structures of spike protein and main protease both in complex with rutin before and after simulations is shown in Figure 7(A,B). Further, these results were strengthened by the vital contribution of the complex binding energies throughout the simulation run. These calculations with consistent high binding energies and large hydrogen bonds involvement demonstrated the stable binding of rutin with both spike and main protease proteins.
Figure 6. Representation of the MD analysis plots for the hydrogen bond with the surrounding solvent against (A) spike protein (purple color) and spike-rutin complex (red color), and (B) main protease (blue) and main protease-rutin complex (green). The MD analysis plots representing hydrogen bonds with the potent inhibitor i.e. rutin against (C) main protease and (D) main protease-rutin complex.
Figure 7. Superimposed structures of (A) spike protein before simulation (green) and after simulation (blue), and (B) main protease before simulation (yellow) and after simulation (pink).
4. Conclusion
SARS-CoV-2 which causes the deadly and highly infectious COVID-19 remains an unmet medical issue that requires urgent attention, owing to the number of deaths it is causing globally. We report here a computational approach for a repurposing screening program against two crucial SARS-CoV-2 drug targets viz. the spike protein and the main protease. We discovered that the drug rutin, showed a high binding potency against both the drug targets and may therefore be capable of eliciting reduction in viral load. The current study delineates that the highly effective repurposed drug rutin has the potential to obtain strong inhibition against SARS-CoV-2 spike protein and main protease. Rutin showed high binding efficiency against spike protein and main protease having an XP Glide score of –8.367 kcal/mol and –11.553 kcal/mol.
Further, the molecular interactions of rutin with both the drug targets were analysed thoroughly with the help of molecular dynamics simulation studies. The simulations for spike protein and main protease provided us with the stable trajectory analysis for both the targets. We analyzed various parameters which include RMSD, RMSF, radius of gyration, SASA, binding energies, hydrogen bond interactions, secondary structure contents, and cluster formation throughout the simulation, which proposed insights about the stable binding and compactness due to strong interaction of the rutin with both the drug targets. Through this in silico approach, the repurposed drug, rutin depicted reliable inhibition of two vital proteins and therefore may lead to elicitation of anti-SARS-CoV-2 activity.
Supplementary Material
Supporting_Information_Revised.pdf
Click here for additional data file.
|
|
sec
|
1. Introduction
Pandemics involving pathogenic human coronaviruses have wreaked havoc on the human populace. During the severe acute respiratory syndrome coronavirus (SARS-CoV) epidemic in 2002–2003, a high fatality rate of 10% was observed in the approximately 8,000 individuals infected (Marra et al., 2003; Peiris et al., 2003; Rota et al., 2003; SARS Working Group, 2003). In 2012, the middle east respiratory syndrome coronavirus (MERS-CoV) had a much higher fatality rate of 36% and infected more than 1,700 individuals (de Groot et al., 2013; Zaki et al., 2012). However, the recent SARS-CoV outbreak of 2019, also referred to as SARS-CoV-2 or COVID-19 (Chen et al., 2020), began in Wuhan in China and spread rapidly due to close human-to-human proximity (Li et al., 2020). The SARS-CoV-2 has caused significant devastation globally due to the sheer number of people infected. According to statistics published by the World Health Organization (WHO), there were 4,334,451 known cases of Covid-19 worldwide, out of which 19,737,794 cases alone are reported from America, while Europe attributes 9,664,042 cases, causing death casualties of over 1,157,509 globally as of Oct 27, 2020.
Similar to SARS-CoV and MERS-CoV, SARS-CoV-2 (exhibiting 80% and 50% homology, respectively) (Kim et al., 2020; Zhu et al., 2020) belongs to the genera beta coronavirus (family Coronaviridae in the order Nidovirales) (Enjuanes et al., 2006; Perlman & Netland, 2009), which is known to infect mammals (Li, 2016) and has manifested an illustrious capability of cross-species transmission, including humans (Menachery et al., 2017). The SARS-CoV-2 is an enveloped virus characterized by a positive-sense, single-stranded RNA genome of almost 30 kb, surrounded by a helical capsid comprised of nucleo capsid protein (N). The viral envelope is associated with three primary structural proteins viz. membrane proteins (M) and envelope proteins (E), which perform virus assembly; and spike proteins (S), which facilitate the viral attachment and thus, the virus entrance into the host cells. The large protrusions formed by spikes on the surface of the virus give it the appearance of having crowns, which, in Latin, translates to corona (Li, 2016). These shorter, sgRNA-encoded structural proteins and several accessory proteins are known to be conserved (Kim et al., 2020). The large ectodomain of the spike consists of receptor binding S1 and membrane fusion S2 subunits. For many CoVs, the S1 and S2 domains remain non-covalently linked. In β coronaviruses, the cleavage between the S1 and S2 regions is not obligatory. However, the host proteases have been observed to cause cleavage within the fusion domain (S2), which leads to irreversible conformational changes, activating the protein for membrane fusion (Zhou et al., 2020). The S2 subunit contains two regions with 4, 3 hydrophobic heptad repeats (HR), HR1 and HR2, which are conserved in sequence and position. The HR2 is located adjacent to the trans-membrane domain, while the HR1 occurs ∼170 residues upstream of HR2. The HR region has been observed to be a common motif in several viral fusion proteins. Structurally, the HR1 domain exists as a homo-trimeric coil, packing the HR2 domain (in an anti-parallel manner) in its hydrophobic grooves, thereby bringing the N-terminal fusion peptide closer to the trans-membrane anchor, further facilitating fusion due to proximity (Bosch et al., 2003). Various host receptors viz. angiotensin-converting enzyme 2 (ACE2), amino peptidase N (APN), dipeptidyl peptidase 4 (DPP4), carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1), and sugar recognize S1 domain of the spike protein to render virus entry into cells (Li, 2016). Several reports have been published recently discovering inhibitors against the target via a computational methodologies, which involve the usage of FDA approved compounds and natural compounds (Wei et al., 2020), spice molecules (Rout et al., 2020), small-molecule compounds of ZINC Drug Database (Kadioglu, 2020; Wu et al., 2020) along with traditional Chinese medicine and natural products and derivatives (Wu et al., 2020) and medicinal compounds (Salman et al., 2020).
Apart from structural proteins which play a crucial role in the lifecycle of the virus, remarkable number of functional proteins categorized as non-structural proteins (Nsp’s), involved in viral replication, transcription, translation and protein modifications, and host infection are equally important (Wu et al., 2020). Among many such proteins, a well characterized distinguished drug target for SARS-Cov-2 is the 33.8 kDa, 3 C-like main protease, also called as the main protease or Nsp5 (3CLpro or Mpro or Nsp5) (Wu et al., 2020; Zhang et al., 2020). Together with the papain like proteases (Plpro), the enzyme is responsible for the processing of polyproteins produced by the viral RNA (Zhang et al., 2020). The Mpro along with Plpro cleaves the 790 kDa polyprotein1ab to generate 15 Nsp’s. The Mpro processes the polyprotein1ab at 11 sites to produce mature Nsp4-Nsp16 (Wu et al., 2020). Recently, the crystal structure of Mpro has also been reported at a resolution of 1.75 Å. The Mpro structure consists of three domains, domains I and II comprising of six-stranded anti-parallel β barrels spanning from residues 10 to 99 and 100 to 182, respectively forming the substrate-binding site between them and the domain III which is responsible for modulating the Mpro dimerization is a cluster of five helices ranging from residues 198 to 303 (Zhang et al., 2020). Due to the availability of the crystal structure and the absence of any human homologue, the Mpro is a lucrative drug target to work upon for the discovery of novel antiviral agents (Jin et al., 2020; Wu et al., 2020; Zhang et al., 2020). Recently, numerous studies have emerged against this target, discovering novel inhibitors by utilizing computational approaches which include the use of natural molecules (Aanouz et al., 2020; Das et al., 2020), commercially available drugs and zinc library (Das et al., 2020; Elmezayen, 2020; Ton et al., 2020), antiviral compounds (Khan, 2020; Kumar et al., 2020; Muralidharan, 2020), peptide molecules (Pant, 2020), generative chemistry approaches for drug design (Alex, 2020), combinatorial strategy of using anti-virals, natural products, anti-fungals, anti-protozoals and anti-nematodes (Das et al., 2020), and spice molecules (Rout et al., 2020).
Many pharmacological agents such as antiviral agent’s viz. remdesvir, ribavirin, favipiravir, chloroquine, hydroxychloroquine, oseltamivir and umifenovir; immunomodulatory agents including toclizumab, inteferons; adjunctive agents like azithromycin, corticosteroids have been repurposed for treating COVID-19. However, so far no therapy has been (Lam et al., 2020) recognized as an effective treatment against the deadly disease. It has, therefore, become the need of the hour to adopt computational approaches for drug discovery and repurpose other existing marketed drugs for its treatment. Drug repurposing via high-throughput screening of various databases is a cost-effective and time-saving approach. The present report pertains to the above-described drug discovery effort using a repurposing screening strategy to target the SARS-CoV-2 spike protein and main protease. Together, these proteins play a key role in binding to host cell receptors and causing fusion to render viral entry, determining the host range and tissue tropism, eliciting an immune response by the host, viral replication, transcription and protein processing (Li, 2016).
|
|
label
|
1.
|
|
title
|
Introduction
|
|
p
|
Pandemics involving pathogenic human coronaviruses have wreaked havoc on the human populace. During the severe acute respiratory syndrome coronavirus (SARS-CoV) epidemic in 2002–2003, a high fatality rate of 10% was observed in the approximately 8,000 individuals infected (Marra et al., 2003; Peiris et al., 2003; Rota et al., 2003; SARS Working Group, 2003). In 2012, the middle east respiratory syndrome coronavirus (MERS-CoV) had a much higher fatality rate of 36% and infected more than 1,700 individuals (de Groot et al., 2013; Zaki et al., 2012). However, the recent SARS-CoV outbreak of 2019, also referred to as SARS-CoV-2 or COVID-19 (Chen et al., 2020), began in Wuhan in China and spread rapidly due to close human-to-human proximity (Li et al., 2020). The SARS-CoV-2 has caused significant devastation globally due to the sheer number of people infected. According to statistics published by the World Health Organization (WHO), there were 4,334,451 known cases of Covid-19 worldwide, out of which 19,737,794 cases alone are reported from America, while Europe attributes 9,664,042 cases, causing death casualties of over 1,157,509 globally as of Oct 27, 2020.
|
|
p
|
Similar to SARS-CoV and MERS-CoV, SARS-CoV-2 (exhibiting 80% and 50% homology, respectively) (Kim et al., 2020; Zhu et al., 2020) belongs to the genera beta coronavirus (family Coronaviridae in the order Nidovirales) (Enjuanes et al., 2006; Perlman & Netland, 2009), which is known to infect mammals (Li, 2016) and has manifested an illustrious capability of cross-species transmission, including humans (Menachery et al., 2017). The SARS-CoV-2 is an enveloped virus characterized by a positive-sense, single-stranded RNA genome of almost 30 kb, surrounded by a helical capsid comprised of nucleo capsid protein (N). The viral envelope is associated with three primary structural proteins viz. membrane proteins (M) and envelope proteins (E), which perform virus assembly; and spike proteins (S), which facilitate the viral attachment and thus, the virus entrance into the host cells. The large protrusions formed by spikes on the surface of the virus give it the appearance of having crowns, which, in Latin, translates to corona (Li, 2016). These shorter, sgRNA-encoded structural proteins and several accessory proteins are known to be conserved (Kim et al., 2020). The large ectodomain of the spike consists of receptor binding S1 and membrane fusion S2 subunits. For many CoVs, the S1 and S2 domains remain non-covalently linked. In β coronaviruses, the cleavage between the S1 and S2 regions is not obligatory. However, the host proteases have been observed to cause cleavage within the fusion domain (S2), which leads to irreversible conformational changes, activating the protein for membrane fusion (Zhou et al., 2020). The S2 subunit contains two regions with 4, 3 hydrophobic heptad repeats (HR), HR1 and HR2, which are conserved in sequence and position. The HR2 is located adjacent to the trans-membrane domain, while the HR1 occurs ∼170 residues upstream of HR2. The HR region has been observed to be a common motif in several viral fusion proteins. Structurally, the HR1 domain exists as a homo-trimeric coil, packing the HR2 domain (in an anti-parallel manner) in its hydrophobic grooves, thereby bringing the N-terminal fusion peptide closer to the trans-membrane anchor, further facilitating fusion due to proximity (Bosch et al., 2003). Various host receptors viz. angiotensin-converting enzyme 2 (ACE2), amino peptidase N (APN), dipeptidyl peptidase 4 (DPP4), carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1), and sugar recognize S1 domain of the spike protein to render virus entry into cells (Li, 2016). Several reports have been published recently discovering inhibitors against the target via a computational methodologies, which involve the usage of FDA approved compounds and natural compounds (Wei et al., 2020), spice molecules (Rout et al., 2020), small-molecule compounds of ZINC Drug Database (Kadioglu, 2020; Wu et al., 2020) along with traditional Chinese medicine and natural products and derivatives (Wu et al., 2020) and medicinal compounds (Salman et al., 2020).
|
|
p
|
Apart from structural proteins which play a crucial role in the lifecycle of the virus, remarkable number of functional proteins categorized as non-structural proteins (Nsp’s), involved in viral replication, transcription, translation and protein modifications, and host infection are equally important (Wu et al., 2020). Among many such proteins, a well characterized distinguished drug target for SARS-Cov-2 is the 33.8 kDa, 3 C-like main protease, also called as the main protease or Nsp5 (3CLpro or Mpro or Nsp5) (Wu et al., 2020; Zhang et al., 2020). Together with the papain like proteases (Plpro), the enzyme is responsible for the processing of polyproteins produced by the viral RNA (Zhang et al., 2020). The Mpro along with Plpro cleaves the 790 kDa polyprotein1ab to generate 15 Nsp’s. The Mpro processes the polyprotein1ab at 11 sites to produce mature Nsp4-Nsp16 (Wu et al., 2020). Recently, the crystal structure of Mpro has also been reported at a resolution of 1.75 Å. The Mpro structure consists of three domains, domains I and II comprising of six-stranded anti-parallel β barrels spanning from residues 10 to 99 and 100 to 182, respectively forming the substrate-binding site between them and the domain III which is responsible for modulating the Mpro dimerization is a cluster of five helices ranging from residues 198 to 303 (Zhang et al., 2020). Due to the availability of the crystal structure and the absence of any human homologue, the Mpro is a lucrative drug target to work upon for the discovery of novel antiviral agents (Jin et al., 2020; Wu et al., 2020; Zhang et al., 2020). Recently, numerous studies have emerged against this target, discovering novel inhibitors by utilizing computational approaches which include the use of natural molecules (Aanouz et al., 2020; Das et al., 2020), commercially available drugs and zinc library (Das et al., 2020; Elmezayen, 2020; Ton et al., 2020), antiviral compounds (Khan, 2020; Kumar et al., 2020; Muralidharan, 2020), peptide molecules (Pant, 2020), generative chemistry approaches for drug design (Alex, 2020), combinatorial strategy of using anti-virals, natural products, anti-fungals, anti-protozoals and anti-nematodes (Das et al., 2020), and spice molecules (Rout et al., 2020).
|
|
p
|
Many pharmacological agents such as antiviral agent’s viz. remdesvir, ribavirin, favipiravir, chloroquine, hydroxychloroquine, oseltamivir and umifenovir; immunomodulatory agents including toclizumab, inteferons; adjunctive agents like azithromycin, corticosteroids have been repurposed for treating COVID-19. However, so far no therapy has been (Lam et al., 2020) recognized as an effective treatment against the deadly disease. It has, therefore, become the need of the hour to adopt computational approaches for drug discovery and repurpose other existing marketed drugs for its treatment. Drug repurposing via high-throughput screening of various databases is a cost-effective and time-saving approach. The present report pertains to the above-described drug discovery effort using a repurposing screening strategy to target the SARS-CoV-2 spike protein and main protease. Together, these proteins play a key role in binding to host cell receptors and causing fusion to render viral entry, determining the host range and tissue tropism, eliciting an immune response by the host, viral replication, transcription and protein processing (Li, 2016).
|
|
sec
|
2. Material and methods
2.1. Preparation of protein and ligand structures
The spike protein and main protease regulatory enzyme of SARS-CoV-2 have been revealed to play a vital role in COVID-19 infection in the host (Ke et al., 2020; Xue et al., 2008). In order to delineate the drug treatment, the deduced three dimensional structures of target SAR-CoV-2 spike protein in open state (PDB ID: 6VYB) and main protease (PDB ID: 6LU7) were retrieved from RCSB Protein Data Bank (Burley et al., 2019). The crystal structures of target proteins were imported to Schrödinger Maestro (Moore, 2015) for molecular docking analysis. Pre-processing and preparation of target proteins were performed using Schrödinger Protein Preparation Wizard (Schrödinger, 2011) by adding the missing hydrogen atoms, correcting the bond orders, capping the protein termini, creating disulfide bonds and removing the water molecules. Also, the H-bonds were optimized and further the target proteins were equilibrated using OPLS 2005 force field (Robertson et al., 2015).
The library of ligands was downloaded from DrugBank (Wishart et al., 2018) and the coordinate files were prepared and optimized using LigPrep module (Release, 2017). The energy minimization of ligand structures was carried using OPLS 2005 force field and ionization states were generated at pH 7.0 ± 2.0 using Epik. This generated an output file with various stereoisomers and tautomeric conformers of each ligand producing chemical and structural diversity.
2.2. Determination of binding pocket and generation of receptor grid
The active sites enclosed in the binding pocket of spike protein and main protease were predicted using Computed Atlas of Surface Topography of proteins (CASTp) webserver (Tian et al., 2018) which calculates geometric and topological measurements of protein structures. It measures the volume and surface area of the binding pocket and also predicts the atoms involved in formation of this binding cavity. Although the active sites of spike protein and the main protease are known but still research is undergoing in which it is predicted that new amino acids might be involved in the inhibition of virus activity. In order to undertake all the possibilities, binding pocket analysis was performed through CastP. Further for the ease of molecular docking, a receptor grid was generated around the active sites of target proteins using Receptor Grid generation module. The receptor grid file was created with a grid box around the centroid of selected active residues.
2.3. Protein ligand interactions and molecular docking
Protein ligand docking was conducted using GLIDE module of Schrödinger (Release, 2018) which performs all the docking calculations using OPLS 2005 force field. A total of 2652 FDA approved drugs were screened against spike protein and main protease of SAR-CoV-2 through molecular docking approaches. The virtual screening of compounds against the target protein was conducted using two docking methodologies, initially through a rapid screening of large set of compounds using high throughput virtual screening (HTVS) method followed by further screening a subset of top scoring HTVS compounds with a more precise and accurate extra precise (XP) docking method. Protein ligand complexes were ranked using GlideScore function to predict the binding efficacy of ligands with the protein. The docked complex with the lowest docking score was selected for further molecular dynamic simulation analysis to elucidate the inhibition mechanism against SARS-CoV-2 infection.
2.4. Molecular dynamics trajectory analysis of protein and their docked complex
The SARS-CoV-2 spike protein and main protease were used separately for molecular dynamics (MD) simulations alone and in complex with the repurposed drug, rutin. All the four MD simulations were executed through GROMACS v5.0 under the force field GROMOS96 54a7 having water model SPC216 along with the time step of 1 fs for 100 ns (Abraham et al., 2015; Darden et al., 1993). Varied sizes of the simulation box were created for each MD simulation event, which were further loaded with about respective amount of water molecules using the SPC model. The total charge on spike protein, main protease and the two in complex with rutin were neutralized by adding –6, –4, –6, and –4 charges, respectively, and were incorporated into the simulation system by compensating the water molecules in the arbitrary locations inside the simulation box. The NPT ensembles, along with periodic boundary conditions, were utilized to carry out MD simulations. A cut-off of about 12 Å was used in order to manage the Vander Waals forces. The Particle Mesh Ewald model manifesting a cut-off of 14 Å was further utilized to calculate the electrostatic interactions (Darden et al., 1993). The spike protein, main protease and the two in complex with rutin were solvated through a slab of about 10 Å in every direction. The neighbor list was updated to a frequency of 10 ps.
The MD simulations were achieved for each system employing the four major steps. The first step deals with the energy minimization of the entire system utilizing the integrator of steepest descent in continuation with second integrator of conjugate gradients algorithms. The second step involves the minimization and molecular dynamics of NVT and NPT ensembles for 500 ps and 1000 ps, respectively allowing the solvents and ions to evolve by keeping the same starting configuration for the structures. In the third step the systems were heated using a lower temperature coupling (τ = 0.1 ps) along with pressure coupling (τ = 0.5 ps) to attain equilibrium at 300 K and 1 atm of temperature and pressure. In the equilibration phase, the thermostat and barostat were evaluated through the Berendsen algorithm (Berendsen et al., 1984). The hydrogen-containing bond lengths were constrained with the help of the LINCS algorithm (Hess et al., 1997). Finally, the last step also called as the production step was carried out, where the MD simulation for 100 ns at 300 K temperature with 2 fs of time step were performed for both systems, and the last final structures were achieved. The Maxwell Boltzmann distribution was utilized in order to reassign the velocities at every step. The Nose Hoover thermostat and Parrinello Rahman barostat were the respective thermostat and barostat for the final MD or production run (Berendsen et al., 1984).
Various analyses were performed with the help of inbuilt analysis commands of GROMACS. The root mean square deviation (RMSD) is a magnitude of the dimensional disparity among the two stagnant structures, and RMSD calculation is achieved depending upon the native structure and each consecutive trajectory frames in the simulation. In addition, root mean square fluctuation (RMSF) profile measures the affability of every protein residue depending on the fluctuation about an average location within all MD simulations (Knapp et al., 2011). Therefore, RMSD and RMSF of each simulation system were determined to examine the stability and residual fluctuations. Further, the radius of gyration (Rg) analysis was performed to evaluate the compactness of both the simulation systems separately. Also, the hydrogen bond analysis was performed to check the neighboring interactions with both simulation systems separately, including the hydrophobic interactions with the help of the LigPlot tool for both spike protein and main protease complexes with rutin before and after simulation.
Additionally, the solvent accessibility surface area (SASA) was also computed to examine the solvent attributable areas of all simulation system. The cluster analysis having a cut-off value of 0.25 and 0.2 nm for spike protein and main protease respectively, depending upon the RMSD profile were utilized to demonstrate the conformations found utmost intermittently throughout the trajectory. Here, all the structures having RMSD values of below 0.25 nm for all components within a cluster are incorporated to the initial cluster. It is rare that a molecule having a higher value for RMSD than 0.25 nm from other cluster supposedly would be treated as a structure. The secondary structure analysis was also performed using the DSSP program (Martin et al., 2005). The visualization of protein nature during the entire simulation was accomplished by using Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) and UCSF Chimera (Pettersen et al., 2004).
|
|
label
|
2.
|
|
title
|
Material and methods
|
|
sec
|
2.1. Preparation of protein and ligand structures
The spike protein and main protease regulatory enzyme of SARS-CoV-2 have been revealed to play a vital role in COVID-19 infection in the host (Ke et al., 2020; Xue et al., 2008). In order to delineate the drug treatment, the deduced three dimensional structures of target SAR-CoV-2 spike protein in open state (PDB ID: 6VYB) and main protease (PDB ID: 6LU7) were retrieved from RCSB Protein Data Bank (Burley et al., 2019). The crystal structures of target proteins were imported to Schrödinger Maestro (Moore, 2015) for molecular docking analysis. Pre-processing and preparation of target proteins were performed using Schrödinger Protein Preparation Wizard (Schrödinger, 2011) by adding the missing hydrogen atoms, correcting the bond orders, capping the protein termini, creating disulfide bonds and removing the water molecules. Also, the H-bonds were optimized and further the target proteins were equilibrated using OPLS 2005 force field (Robertson et al., 2015).
The library of ligands was downloaded from DrugBank (Wishart et al., 2018) and the coordinate files were prepared and optimized using LigPrep module (Release, 2017). The energy minimization of ligand structures was carried using OPLS 2005 force field and ionization states were generated at pH 7.0 ± 2.0 using Epik. This generated an output file with various stereoisomers and tautomeric conformers of each ligand producing chemical and structural diversity.
|
|
label
|
2.1.
|
|
title
|
Preparation of protein and ligand structures
|
|
p
|
The spike protein and main protease regulatory enzyme of SARS-CoV-2 have been revealed to play a vital role in COVID-19 infection in the host (Ke et al., 2020; Xue et al., 2008). In order to delineate the drug treatment, the deduced three dimensional structures of target SAR-CoV-2 spike protein in open state (PDB ID: 6VYB) and main protease (PDB ID: 6LU7) were retrieved from RCSB Protein Data Bank (Burley et al., 2019). The crystal structures of target proteins were imported to Schrödinger Maestro (Moore, 2015) for molecular docking analysis. Pre-processing and preparation of target proteins were performed using Schrödinger Protein Preparation Wizard (Schrödinger, 2011) by adding the missing hydrogen atoms, correcting the bond orders, capping the protein termini, creating disulfide bonds and removing the water molecules. Also, the H-bonds were optimized and further the target proteins were equilibrated using OPLS 2005 force field (Robertson et al., 2015).
|
|
p
|
The library of ligands was downloaded from DrugBank (Wishart et al., 2018) and the coordinate files were prepared and optimized using LigPrep module (Release, 2017). The energy minimization of ligand structures was carried using OPLS 2005 force field and ionization states were generated at pH 7.0 ± 2.0 using Epik. This generated an output file with various stereoisomers and tautomeric conformers of each ligand producing chemical and structural diversity.
|
|
sec
|
2.2. Determination of binding pocket and generation of receptor grid
The active sites enclosed in the binding pocket of spike protein and main protease were predicted using Computed Atlas of Surface Topography of proteins (CASTp) webserver (Tian et al., 2018) which calculates geometric and topological measurements of protein structures. It measures the volume and surface area of the binding pocket and also predicts the atoms involved in formation of this binding cavity. Although the active sites of spike protein and the main protease are known but still research is undergoing in which it is predicted that new amino acids might be involved in the inhibition of virus activity. In order to undertake all the possibilities, binding pocket analysis was performed through CastP. Further for the ease of molecular docking, a receptor grid was generated around the active sites of target proteins using Receptor Grid generation module. The receptor grid file was created with a grid box around the centroid of selected active residues.
|
|
label
|
2.2.
|
|
title
|
Determination of binding pocket and generation of receptor grid
|
|
p
|
The active sites enclosed in the binding pocket of spike protein and main protease were predicted using Computed Atlas of Surface Topography of proteins (CASTp) webserver (Tian et al., 2018) which calculates geometric and topological measurements of protein structures. It measures the volume and surface area of the binding pocket and also predicts the atoms involved in formation of this binding cavity. Although the active sites of spike protein and the main protease are known but still research is undergoing in which it is predicted that new amino acids might be involved in the inhibition of virus activity. In order to undertake all the possibilities, binding pocket analysis was performed through CastP. Further for the ease of molecular docking, a receptor grid was generated around the active sites of target proteins using Receptor Grid generation module. The receptor grid file was created with a grid box around the centroid of selected active residues.
|
|
sec
|
2.3. Protein ligand interactions and molecular docking
Protein ligand docking was conducted using GLIDE module of Schrödinger (Release, 2018) which performs all the docking calculations using OPLS 2005 force field. A total of 2652 FDA approved drugs were screened against spike protein and main protease of SAR-CoV-2 through molecular docking approaches. The virtual screening of compounds against the target protein was conducted using two docking methodologies, initially through a rapid screening of large set of compounds using high throughput virtual screening (HTVS) method followed by further screening a subset of top scoring HTVS compounds with a more precise and accurate extra precise (XP) docking method. Protein ligand complexes were ranked using GlideScore function to predict the binding efficacy of ligands with the protein. The docked complex with the lowest docking score was selected for further molecular dynamic simulation analysis to elucidate the inhibition mechanism against SARS-CoV-2 infection.
|
|
label
|
2.3.
|
|
title
|
Protein ligand interactions and molecular docking
|
|
p
|
Protein ligand docking was conducted using GLIDE module of Schrödinger (Release, 2018) which performs all the docking calculations using OPLS 2005 force field. A total of 2652 FDA approved drugs were screened against spike protein and main protease of SAR-CoV-2 through molecular docking approaches. The virtual screening of compounds against the target protein was conducted using two docking methodologies, initially through a rapid screening of large set of compounds using high throughput virtual screening (HTVS) method followed by further screening a subset of top scoring HTVS compounds with a more precise and accurate extra precise (XP) docking method. Protein ligand complexes were ranked using GlideScore function to predict the binding efficacy of ligands with the protein. The docked complex with the lowest docking score was selected for further molecular dynamic simulation analysis to elucidate the inhibition mechanism against SARS-CoV-2 infection.
|
|
sec
|
2.4. Molecular dynamics trajectory analysis of protein and their docked complex
The SARS-CoV-2 spike protein and main protease were used separately for molecular dynamics (MD) simulations alone and in complex with the repurposed drug, rutin. All the four MD simulations were executed through GROMACS v5.0 under the force field GROMOS96 54a7 having water model SPC216 along with the time step of 1 fs for 100 ns (Abraham et al., 2015; Darden et al., 1993). Varied sizes of the simulation box were created for each MD simulation event, which were further loaded with about respective amount of water molecules using the SPC model. The total charge on spike protein, main protease and the two in complex with rutin were neutralized by adding –6, –4, –6, and –4 charges, respectively, and were incorporated into the simulation system by compensating the water molecules in the arbitrary locations inside the simulation box. The NPT ensembles, along with periodic boundary conditions, were utilized to carry out MD simulations. A cut-off of about 12 Å was used in order to manage the Vander Waals forces. The Particle Mesh Ewald model manifesting a cut-off of 14 Å was further utilized to calculate the electrostatic interactions (Darden et al., 1993). The spike protein, main protease and the two in complex with rutin were solvated through a slab of about 10 Å in every direction. The neighbor list was updated to a frequency of 10 ps.
The MD simulations were achieved for each system employing the four major steps. The first step deals with the energy minimization of the entire system utilizing the integrator of steepest descent in continuation with second integrator of conjugate gradients algorithms. The second step involves the minimization and molecular dynamics of NVT and NPT ensembles for 500 ps and 1000 ps, respectively allowing the solvents and ions to evolve by keeping the same starting configuration for the structures. In the third step the systems were heated using a lower temperature coupling (τ = 0.1 ps) along with pressure coupling (τ = 0.5 ps) to attain equilibrium at 300 K and 1 atm of temperature and pressure. In the equilibration phase, the thermostat and barostat were evaluated through the Berendsen algorithm (Berendsen et al., 1984). The hydrogen-containing bond lengths were constrained with the help of the LINCS algorithm (Hess et al., 1997). Finally, the last step also called as the production step was carried out, where the MD simulation for 100 ns at 300 K temperature with 2 fs of time step were performed for both systems, and the last final structures were achieved. The Maxwell Boltzmann distribution was utilized in order to reassign the velocities at every step. The Nose Hoover thermostat and Parrinello Rahman barostat were the respective thermostat and barostat for the final MD or production run (Berendsen et al., 1984).
Various analyses were performed with the help of inbuilt analysis commands of GROMACS. The root mean square deviation (RMSD) is a magnitude of the dimensional disparity among the two stagnant structures, and RMSD calculation is achieved depending upon the native structure and each consecutive trajectory frames in the simulation. In addition, root mean square fluctuation (RMSF) profile measures the affability of every protein residue depending on the fluctuation about an average location within all MD simulations (Knapp et al., 2011). Therefore, RMSD and RMSF of each simulation system were determined to examine the stability and residual fluctuations. Further, the radius of gyration (Rg) analysis was performed to evaluate the compactness of both the simulation systems separately. Also, the hydrogen bond analysis was performed to check the neighboring interactions with both simulation systems separately, including the hydrophobic interactions with the help of the LigPlot tool for both spike protein and main protease complexes with rutin before and after simulation.
Additionally, the solvent accessibility surface area (SASA) was also computed to examine the solvent attributable areas of all simulation system. The cluster analysis having a cut-off value of 0.25 and 0.2 nm for spike protein and main protease respectively, depending upon the RMSD profile were utilized to demonstrate the conformations found utmost intermittently throughout the trajectory. Here, all the structures having RMSD values of below 0.25 nm for all components within a cluster are incorporated to the initial cluster. It is rare that a molecule having a higher value for RMSD than 0.25 nm from other cluster supposedly would be treated as a structure. The secondary structure analysis was also performed using the DSSP program (Martin et al., 2005). The visualization of protein nature during the entire simulation was accomplished by using Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) and UCSF Chimera (Pettersen et al., 2004).
|
|
label
|
2.4.
|
|
title
|
Molecular dynamics trajectory analysis of protein and their docked complex
|
|
p
|
The SARS-CoV-2 spike protein and main protease were used separately for molecular dynamics (MD) simulations alone and in complex with the repurposed drug, rutin. All the four MD simulations were executed through GROMACS v5.0 under the force field GROMOS96 54a7 having water model SPC216 along with the time step of 1 fs for 100 ns (Abraham et al., 2015; Darden et al., 1993). Varied sizes of the simulation box were created for each MD simulation event, which were further loaded with about respective amount of water molecules using the SPC model. The total charge on spike protein, main protease and the two in complex with rutin were neutralized by adding –6, –4, –6, and –4 charges, respectively, and were incorporated into the simulation system by compensating the water molecules in the arbitrary locations inside the simulation box. The NPT ensembles, along with periodic boundary conditions, were utilized to carry out MD simulations. A cut-off of about 12 Å was used in order to manage the Vander Waals forces. The Particle Mesh Ewald model manifesting a cut-off of 14 Å was further utilized to calculate the electrostatic interactions (Darden et al., 1993). The spike protein, main protease and the two in complex with rutin were solvated through a slab of about 10 Å in every direction. The neighbor list was updated to a frequency of 10 ps.
|
|
p
|
The MD simulations were achieved for each system employing the four major steps. The first step deals with the energy minimization of the entire system utilizing the integrator of steepest descent in continuation with second integrator of conjugate gradients algorithms. The second step involves the minimization and molecular dynamics of NVT and NPT ensembles for 500 ps and 1000 ps, respectively allowing the solvents and ions to evolve by keeping the same starting configuration for the structures. In the third step the systems were heated using a lower temperature coupling (τ = 0.1 ps) along with pressure coupling (τ = 0.5 ps) to attain equilibrium at 300 K and 1 atm of temperature and pressure. In the equilibration phase, the thermostat and barostat were evaluated through the Berendsen algorithm (Berendsen et al., 1984). The hydrogen-containing bond lengths were constrained with the help of the LINCS algorithm (Hess et al., 1997). Finally, the last step also called as the production step was carried out, where the MD simulation for 100 ns at 300 K temperature with 2 fs of time step were performed for both systems, and the last final structures were achieved. The Maxwell Boltzmann distribution was utilized in order to reassign the velocities at every step. The Nose Hoover thermostat and Parrinello Rahman barostat were the respective thermostat and barostat for the final MD or production run (Berendsen et al., 1984).
|
|
p
|
Various analyses were performed with the help of inbuilt analysis commands of GROMACS. The root mean square deviation (RMSD) is a magnitude of the dimensional disparity among the two stagnant structures, and RMSD calculation is achieved depending upon the native structure and each consecutive trajectory frames in the simulation. In addition, root mean square fluctuation (RMSF) profile measures the affability of every protein residue depending on the fluctuation about an average location within all MD simulations (Knapp et al., 2011). Therefore, RMSD and RMSF of each simulation system were determined to examine the stability and residual fluctuations. Further, the radius of gyration (Rg) analysis was performed to evaluate the compactness of both the simulation systems separately. Also, the hydrogen bond analysis was performed to check the neighboring interactions with both simulation systems separately, including the hydrophobic interactions with the help of the LigPlot tool for both spike protein and main protease complexes with rutin before and after simulation.
|
|
p
|
Additionally, the solvent accessibility surface area (SASA) was also computed to examine the solvent attributable areas of all simulation system. The cluster analysis having a cut-off value of 0.25 and 0.2 nm for spike protein and main protease respectively, depending upon the RMSD profile were utilized to demonstrate the conformations found utmost intermittently throughout the trajectory. Here, all the structures having RMSD values of below 0.25 nm for all components within a cluster are incorporated to the initial cluster. It is rare that a molecule having a higher value for RMSD than 0.25 nm from other cluster supposedly would be treated as a structure. The secondary structure analysis was also performed using the DSSP program (Martin et al., 2005). The visualization of protein nature during the entire simulation was accomplished by using Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) and UCSF Chimera (Pettersen et al., 2004).
|
|
sec
|
3. Results and discussion
SARS-CoV-2 is an enveloped non-segmented positively stranded RNA infectious virus that mainly affects enteric and respiratory systems (Graham et al., 2013). Its infection and amelioration in the host are majorly governed by several structural proteins (Joshi et al., 2020; Prajapat et al., 2020). Amidst numerous structural proteins, spike protein and main protease of SARS-CoV-2 are stated to show essential impacts on viral replication through proteolytic machinery and involvement in transcription, translation, and amplification of viral proteins (Paules et al., 2020). Spike protein plays a critical role in binding to host cell receptor and is thought to represent a key determinant of the host range restriction (de Wilde, 2017). Main protease on the other hand exclusively cleaves polypeptide sequences after the glutamine residue and to the best of our knowledge is known as an ideal drug target due to the absence of human host cell protease showing this substrate specificity (Ullrich & Nitsche, 2020). The 3 D structures of spike protein and main protease were preprocessed, optimized and minimized to obtain the refined structures for drug repurposing.
3.1. Binding pocket analysis
The binding pocket volume and surface area were determined through CASTp webserver. Binding pocket was forecasted at the surface as well as in the interior of proteins. Binding pocket volume of spike protein and main protease was 22,908 and 319 (SA), respectively which signifies an optimum space for ligand binding. The binding residues as predicted in the binding cavity for spike protein were Y38-S45, V47-L54, P85-D88,T108,K195-Y200,K202-Y204,P225-D228,I233-R237,Q271-T274,C291-A292,E298,C301-K304,Y313-N317,R319-V327,D364-S366,Y369-N370,Y380-T393,P412,D427-T430,F515-H519,A522,P527-T531,N540-T549,Q564-D568,A570-T573,V576-R577,P579,I587,P589-F592,G594,M740-Y741,G744,F855-N856,Q954,Q957-A958,N960-D979,R983,V991,D994-R995,T998-R1000,Q1002-S1003,T1006-Y1007,Q1010, and R1014 around which the receptor grid was formed. The receptor grid was also created around the catalytic region (T24-L27, H41, C44-S46, M49, P52, Y54, F140-C145, H163-P168, H172, D187-T190, and Q192) of main protease.
3.2. Probable drugs against spike protein and main protease of SARS-CoV-2
The molecular docking studies were performed to determine the ligand interactions and their binding affinity with SARS-CoV-2 spike protein and main protease. The prepared ligands structures were docked against the binding sites of the target proteins. The HTVS docking approach filtered out a large number of compounds on the basis of their binding interactions within the binding sites of the target proteins. The compounds displaying a HTVS glide score of more than −6.0 kJ/mol were employed for XP docking. The docking simulations evaluated a high Glide score for the compounds with best binding affinity and their interactions within the binding pocket of the binding site residues. The docking analysis revealed a common compound, rutin (DrugBank ID: DB01698) which binds strongly to both the major targets of SARS-CoV-2. Our study here presents rutin exhibiting high binding efficiency against spike protein and main protease with an XP Glide score of −8.367 kcal/mol and −11.553 kcal/mol respectively. Rutin (3,3′,4′,5,7-pentahydroxyflavone-3-rhamnoglucoside) belongs to the flavonol class of compound, generously present in plants including buckwheat, apple, passion flower and tea. It is well known to exhibit many bioactivities, such as anti-viral, anti-bacterial, ant-larvicidal, cytoprotective, anticarcinogenic, antioxidant, vasoprotective, cardioprotective and neuroprotective activities (Saluja & Ganeshpurkar, 2016). The list of best ten candidates derived from the docking calculations with their corresponding docking score against spike protein and main protease are shown in supporting information (supporting material Tables S2 and S3). To determine the interaction pattern between rutin and target proteins, their docked complexes were visualized using Ligplot (Figure 1). The interaction pattern between rutin and spike protein revealed the formation of five hydrogen bonds and four hydrophobic interactions at the binding site, thereby contributing to the stability of the complex. The residues involved in the formation of hydrogen bonds between rutin and spike protein were F970, N969, H49, Q52, and T274. While the hydrophobic contact between the two involved residues T51, S50, S967 and S968. The interactions between rutin and main protease was stabilized by residues E166, T190, Y54, D187, T26, N142 forming hydrogen bond and various hydrophobic contacts with residues P168, L167, H41, G143, T25, L27, M49, M165, H164, R188 and Q189 (Tables 1 and 2). A similar study exploiting a blind molecular docking approach utilizing the Swiss Dock server also highlights main protease inhibition by rutin where by common amino acids viz. E166, T190, N142, H41 are found to be interacting with the inhibitor. In our study however, the inhibitor displayed slightly better binding efficiency than their report (Das et al., 2020). From the contact analyses of both proteins, it can be attributed that rutin has high affinity and wide molecular contacts for both spike protein and main protease. Also, it was further scrutinized by molecular dynamics simulations to attain insights towards the inhibitory aspects and efficacy to combat spike protein and main protease drug targets of SARS-CoV-2.
Figure 1. Protein-ligand interactions diagram of spike protein and main protease in complex with rutin. 2D Ligplot representation showing hydrogen bond and hydrophobic interactions, and 3D structure representation showing rutin in the binding pocket of spike protein (A and C) and main protease (B and D).
Table 1. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between spike protein with Rutin.
Spike protein + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
His 49.A NE2–Rutin.het O15 3.28 HNE2-H--O15 Ser50, Thr51, Ser967, Ser968
GLN 52.A OE1–Rutin.het O14 2.89 HOE1-H--O14
THR 274.A OG1–Rutin.het O14 3.30 HOG1-H--O14
ASN 969.A ND2–Rutin.het O10 3.17 HND2-H--O10
PHE 970.A N–Rutin.het O9 3.10 HN-H--O9
Table 2. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between main protease with Rutin.
Main protease + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
THR 26.A O–Rutin.het O15 3.05 HO-H--O15 Thr25, Leu27, His41, Met49, Gly143, His164,
THR 26.A O–Rutin.het O16 3.02 HO-H--O16 Met165, Leu167, Pro168, Arg188,
TYR 54.A OH–Rutin.het O14 2.91 HOH-H--O14 Gln189
ASN 142.A OD1–Rutin.het O5 3.01 HOD1-H--O5
ASN 142.A OD1–Rutin.het O7 2.91 HOD1-H--O7
GLU 166.A O–Rutin.het O8 3.02 HO-H--O8
ASP 187.A O–Rutin.het O14 3.31 HO-H--O14
THR 190.A O–Rutin.het O10 2.87 HO-H--O10
3.3. Molecular dynamics simulation of spike protein and main protease in absence and presence of rutin
The binding affinity of rutin against both spike and main protease were investigated through MD simulations using GROMACS. Separate MD simulations for spike protein and main protease with and without rutin for 100 ns were assessed through trajectory analysis. For the course of 100 ns MD simulation, the stable trajectory was observed and the representative structures were obtained. The deviation of the backbone atoms for simulated structures relative to the starting structures used as a reference was evaluated through RMSD for the backbone atoms (Figures 2 and 3) and C-alpha atoms (supporting material Figure S1). A steep magnitude RMSD variation during the entire simulation can be an implication of a malleable and free instinctive protein or the alteration of the force field. Depending upon the outcomes of the RMSD evaluation, Figure 2(A) represents that the RMSD fluctuation stabilize at about 50 ns MD simulations for both spike protein and spike protein in presence of rutin, and the simulation time was acceptable. In the time from 50–100 ns, the RMSD for spike protein and spike-rutin complex have approximate values about 1.50–1.65 nm and 1.20–1.26 nm, respectively. Similarly, the main protease in absence and presence of rutin for time from 20–70 ns and 20–100 ns showed RMSD values of about 0.25–0.30 nm and 0.20–0.25 nm (Figure 3(A)). Interestingly, RMSD for both the complex systems were more stable than the native state of target receptor spike protein and main protease (Figures 2(A) and 3(A)). The RMSD of C-alpha atoms were also attained and was found to be similar to the RMSD of backbone atoms (supporting material Figure S1). The average RMSD for backbone and C-alpha atoms of both the protein and their complex with rutin are shown in supporting material Table S1.
Figure 2. Representation of the MD analysis plots for spike protein (purple color) and spike-rutin complex (red color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg), and (D) solvent accessible surface area (SASA) analysis.
Figure 3. Representation of the MD analysis plots for main protease (blue color) and main protease-rutin complex (green color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg) and (D) solvent accessible surface area (SASA) analysis. The integral extent of the residual dynamic criterion is achieved by assessing the variations arising from shifts of each of the protein residues that majorly feature the most flexible chain frames. Hence, we validated the residual fluctuations by calculating the mean fluctuation for stable trajectories of each simulation. The RMSF evaluation of all the protein residues were achieved in order to check the residues that may have tend towards an enhancement in the RMSD results. Compelling fluctuations existed in terminal residues and in various loops along β sheets and coils of the spike protein (residues 852-1000) to about 0.47 nm and spike-rutin (residues 62–90, 400–520 and 675–730) to about 0.65 nm (Figure 2(B)). Similarly, in case of main protease and itself in complex with rutin showed fluctuations in the terminal residues and in some regions of β sheet with some coils and loops in main protease (residues 136-148) to about 0.37 nm and main protease-rutin (residues 39-52, 150-156) to about 0.37 nm, respectively (Figure 3(B)). Also, we noticed that residues 270 to 350 in spike protein and 160 to 190 in main protease showed comparatively reduced discrepancy in RMSF values, which were the corresponding binding segment of both the complexes after MD simulation. Interestingly, the binding region for spike-rutin and main protease-rutin before MD simulation were found to be the same after simulation, suggesting as a strong binding interaction of the rutin towards spike protein and main protease and indicating stability of both the complexes.
The radius of gyration analysis was evaluated to determine the change in compactness of protein systems used throughout the MD simulations. The Rg plots for spike, main protease and its complex with rutin show slight fluctuations at the initial frame and attain compactness after 30 ns with Rg score of 3.85 nm and 4 nm (spike and spike-rutin) and after 20 ns with Rg score of 2.27 nm and 2.23 nm (main protease and main protease-rutin), respectively (Figures 2(C) and 3(C)). When compared with spike protein, Rg value for spike-rutin is stabilized and remains constant, suggesting strong binding interaction of the inhibitor, and the same is observed in case of main protease and rutin. Similar observations were determined through SASA analysis representing the solvent defined protein surface and its orientation through folding, making the alterations in the exposed and buried regions of the surface area of proteins. The SASA values for all the simulation systems were about 460 nm/S2/N (spike), 485 nm/S2/N (spike-rutin), 154 nm/S2/N (main-protease), and 152 nm/S2/N (main protease-rutin) (Figures 2(D) and 3(D)). Here, spike-rutin and main protease-rutin solvation profile shows a convincing SASA value suggesting a stable structure and strong binding interaction with the rutin.
Further, the cluster analysis with a RMSD based cut-off value of 0.25 nm demonstrated the development of 13,28,3 and 2 distinctive clusters for spike, spike-rutin, main protease, and main protease-rutin simulation systems respectively (Figure 4(A–D)). The most dominant cluster attained after 100 ns of MD simulation for all four systems are shown in Figure 4(E–H). Also, the secondary structure analyses of the stable trajectory for both the simulation systems were performed using the DSSP tool of GROMACS. Both, spike and main protease were formed mainly of conserved β-sheet along with some connecting loops and α-helix with small coil regions as secondary structure elements infused with various small segments of bend, turn, and β-bridge (Figure 5). Both cluster analysis and secondary structure analysis reveals the conformational changes before and after simulations for spike and main protease with and without rutin structures. A noticeable observation through these analyses supports the rationale that spike protein and main protease binds with rutin firmly and attain similar interactions after MD simulations also.
Figure 4. Representation of the MD analysis plots for the RMSD matrix cluster formation showing number of clusters in (A) spike protein, (B) spike-rutin complex, (C) main protease, and (D) main protease-rutin complex. The most dominant cluster conformation after MD cluster analysis of (E) spike protein, (F) spike-rutin complex, (G) main protease and (H) main protease-rutin complex.
Figure 5. The secondary structure analysis of (A) spike protein, (B) spike-rutin complex, (C) main protease and (D) main protease-rutin complex. Further, hydrogen bond analysis revealed a higher number of hydrogen bonds with the surrounding water molecules in case of spike protein (2300 hydrogen bonds) than in complex with rutin (approximately, 2200 hydrogen bonds). Similarly, in case of main protease, 700 hydrogen bonds were formed while 560 hydrogen bonds in case of main protease complex with rutin (Figure 6(A,B)). In addition, the hydrogen bond landscape assessed against the inhibitor revealed the dynamic equilibration of the complex trajectories with a high number of hydrogen bonds, as shown in Figure 6(C,D). The average number of hydrogen bonds within spike-rutin and main protease-rutin are 3.72 and 2.29, respectively. The consistent numbers of hydrogen bonds were observed which contributed significantly to the proximal binding of the potent drug rutin with the spike and main protease receptor. The superimposed structures of spike protein and main protease both in complex with rutin before and after simulations is shown in Figure 7(A,B). Further, these results were strengthened by the vital contribution of the complex binding energies throughout the simulation run. These calculations with consistent high binding energies and large hydrogen bonds involvement demonstrated the stable binding of rutin with both spike and main protease proteins.
Figure 6. Representation of the MD analysis plots for the hydrogen bond with the surrounding solvent against (A) spike protein (purple color) and spike-rutin complex (red color), and (B) main protease (blue) and main protease-rutin complex (green). The MD analysis plots representing hydrogen bonds with the potent inhibitor i.e. rutin against (C) main protease and (D) main protease-rutin complex.
Figure 7. Superimposed structures of (A) spike protein before simulation (green) and after simulation (blue), and (B) main protease before simulation (yellow) and after simulation (pink).
|
|
label
|
3.
|
|
title
|
Results and discussion
|
|
p
|
SARS-CoV-2 is an enveloped non-segmented positively stranded RNA infectious virus that mainly affects enteric and respiratory systems (Graham et al., 2013). Its infection and amelioration in the host are majorly governed by several structural proteins (Joshi et al., 2020; Prajapat et al., 2020). Amidst numerous structural proteins, spike protein and main protease of SARS-CoV-2 are stated to show essential impacts on viral replication through proteolytic machinery and involvement in transcription, translation, and amplification of viral proteins (Paules et al., 2020). Spike protein plays a critical role in binding to host cell receptor and is thought to represent a key determinant of the host range restriction (de Wilde, 2017). Main protease on the other hand exclusively cleaves polypeptide sequences after the glutamine residue and to the best of our knowledge is known as an ideal drug target due to the absence of human host cell protease showing this substrate specificity (Ullrich & Nitsche, 2020). The 3 D structures of spike protein and main protease were preprocessed, optimized and minimized to obtain the refined structures for drug repurposing.
|
|
sec
|
3.1. Binding pocket analysis
The binding pocket volume and surface area were determined through CASTp webserver. Binding pocket was forecasted at the surface as well as in the interior of proteins. Binding pocket volume of spike protein and main protease was 22,908 and 319 (SA), respectively which signifies an optimum space for ligand binding. The binding residues as predicted in the binding cavity for spike protein were Y38-S45, V47-L54, P85-D88,T108,K195-Y200,K202-Y204,P225-D228,I233-R237,Q271-T274,C291-A292,E298,C301-K304,Y313-N317,R319-V327,D364-S366,Y369-N370,Y380-T393,P412,D427-T430,F515-H519,A522,P527-T531,N540-T549,Q564-D568,A570-T573,V576-R577,P579,I587,P589-F592,G594,M740-Y741,G744,F855-N856,Q954,Q957-A958,N960-D979,R983,V991,D994-R995,T998-R1000,Q1002-S1003,T1006-Y1007,Q1010, and R1014 around which the receptor grid was formed. The receptor grid was also created around the catalytic region (T24-L27, H41, C44-S46, M49, P52, Y54, F140-C145, H163-P168, H172, D187-T190, and Q192) of main protease.
|
|
label
|
3.1.
|
|
title
|
Binding pocket analysis
|
|
p
|
The binding pocket volume and surface area were determined through CASTp webserver. Binding pocket was forecasted at the surface as well as in the interior of proteins. Binding pocket volume of spike protein and main protease was 22,908 and 319 (SA), respectively which signifies an optimum space for ligand binding. The binding residues as predicted in the binding cavity for spike protein were Y38-S45, V47-L54, P85-D88,T108,K195-Y200,K202-Y204,P225-D228,I233-R237,Q271-T274,C291-A292,E298,C301-K304,Y313-N317,R319-V327,D364-S366,Y369-N370,Y380-T393,P412,D427-T430,F515-H519,A522,P527-T531,N540-T549,Q564-D568,A570-T573,V576-R577,P579,I587,P589-F592,G594,M740-Y741,G744,F855-N856,Q954,Q957-A958,N960-D979,R983,V991,D994-R995,T998-R1000,Q1002-S1003,T1006-Y1007,Q1010, and R1014 around which the receptor grid was formed. The receptor grid was also created around the catalytic region (T24-L27, H41, C44-S46, M49, P52, Y54, F140-C145, H163-P168, H172, D187-T190, and Q192) of main protease.
|
|
sec
|
3.2. Probable drugs against spike protein and main protease of SARS-CoV-2
The molecular docking studies were performed to determine the ligand interactions and their binding affinity with SARS-CoV-2 spike protein and main protease. The prepared ligands structures were docked against the binding sites of the target proteins. The HTVS docking approach filtered out a large number of compounds on the basis of their binding interactions within the binding sites of the target proteins. The compounds displaying a HTVS glide score of more than −6.0 kJ/mol were employed for XP docking. The docking simulations evaluated a high Glide score for the compounds with best binding affinity and their interactions within the binding pocket of the binding site residues. The docking analysis revealed a common compound, rutin (DrugBank ID: DB01698) which binds strongly to both the major targets of SARS-CoV-2. Our study here presents rutin exhibiting high binding efficiency against spike protein and main protease with an XP Glide score of −8.367 kcal/mol and −11.553 kcal/mol respectively. Rutin (3,3′,4′,5,7-pentahydroxyflavone-3-rhamnoglucoside) belongs to the flavonol class of compound, generously present in plants including buckwheat, apple, passion flower and tea. It is well known to exhibit many bioactivities, such as anti-viral, anti-bacterial, ant-larvicidal, cytoprotective, anticarcinogenic, antioxidant, vasoprotective, cardioprotective and neuroprotective activities (Saluja & Ganeshpurkar, 2016). The list of best ten candidates derived from the docking calculations with their corresponding docking score against spike protein and main protease are shown in supporting information (supporting material Tables S2 and S3). To determine the interaction pattern between rutin and target proteins, their docked complexes were visualized using Ligplot (Figure 1). The interaction pattern between rutin and spike protein revealed the formation of five hydrogen bonds and four hydrophobic interactions at the binding site, thereby contributing to the stability of the complex. The residues involved in the formation of hydrogen bonds between rutin and spike protein were F970, N969, H49, Q52, and T274. While the hydrophobic contact between the two involved residues T51, S50, S967 and S968. The interactions between rutin and main protease was stabilized by residues E166, T190, Y54, D187, T26, N142 forming hydrogen bond and various hydrophobic contacts with residues P168, L167, H41, G143, T25, L27, M49, M165, H164, R188 and Q189 (Tables 1 and 2). A similar study exploiting a blind molecular docking approach utilizing the Swiss Dock server also highlights main protease inhibition by rutin where by common amino acids viz. E166, T190, N142, H41 are found to be interacting with the inhibitor. In our study however, the inhibitor displayed slightly better binding efficiency than their report (Das et al., 2020). From the contact analyses of both proteins, it can be attributed that rutin has high affinity and wide molecular contacts for both spike protein and main protease. Also, it was further scrutinized by molecular dynamics simulations to attain insights towards the inhibitory aspects and efficacy to combat spike protein and main protease drug targets of SARS-CoV-2.
Figure 1. Protein-ligand interactions diagram of spike protein and main protease in complex with rutin. 2D Ligplot representation showing hydrogen bond and hydrophobic interactions, and 3D structure representation showing rutin in the binding pocket of spike protein (A and C) and main protease (B and D).
Table 1. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between spike protein with Rutin.
Spike protein + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
His 49.A NE2–Rutin.het O15 3.28 HNE2-H--O15 Ser50, Thr51, Ser967, Ser968
GLN 52.A OE1–Rutin.het O14 2.89 HOE1-H--O14
THR 274.A OG1–Rutin.het O14 3.30 HOG1-H--O14
ASN 969.A ND2–Rutin.het O10 3.17 HND2-H--O10
PHE 970.A N–Rutin.het O9 3.10 HN-H--O9
Table 2. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between main protease with Rutin.
Main protease + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
THR 26.A O–Rutin.het O15 3.05 HO-H--O15 Thr25, Leu27, His41, Met49, Gly143, His164,
THR 26.A O–Rutin.het O16 3.02 HO-H--O16 Met165, Leu167, Pro168, Arg188,
TYR 54.A OH–Rutin.het O14 2.91 HOH-H--O14 Gln189
ASN 142.A OD1–Rutin.het O5 3.01 HOD1-H--O5
ASN 142.A OD1–Rutin.het O7 2.91 HOD1-H--O7
GLU 166.A O–Rutin.het O8 3.02 HO-H--O8
ASP 187.A O–Rutin.het O14 3.31 HO-H--O14
THR 190.A O–Rutin.het O10 2.87 HO-H--O10
|
|
label
|
3.2.
|
|
title
|
Probable drugs against spike protein and main protease of SARS-CoV-2
|
|
p
|
The molecular docking studies were performed to determine the ligand interactions and their binding affinity with SARS-CoV-2 spike protein and main protease. The prepared ligands structures were docked against the binding sites of the target proteins. The HTVS docking approach filtered out a large number of compounds on the basis of their binding interactions within the binding sites of the target proteins. The compounds displaying a HTVS glide score of more than −6.0 kJ/mol were employed for XP docking. The docking simulations evaluated a high Glide score for the compounds with best binding affinity and their interactions within the binding pocket of the binding site residues. The docking analysis revealed a common compound, rutin (DrugBank ID: DB01698) which binds strongly to both the major targets of SARS-CoV-2. Our study here presents rutin exhibiting high binding efficiency against spike protein and main protease with an XP Glide score of −8.367 kcal/mol and −11.553 kcal/mol respectively. Rutin (3,3′,4′,5,7-pentahydroxyflavone-3-rhamnoglucoside) belongs to the flavonol class of compound, generously present in plants including buckwheat, apple, passion flower and tea. It is well known to exhibit many bioactivities, such as anti-viral, anti-bacterial, ant-larvicidal, cytoprotective, anticarcinogenic, antioxidant, vasoprotective, cardioprotective and neuroprotective activities (Saluja & Ganeshpurkar, 2016). The list of best ten candidates derived from the docking calculations with their corresponding docking score against spike protein and main protease are shown in supporting information (supporting material Tables S2 and S3). To determine the interaction pattern between rutin and target proteins, their docked complexes were visualized using Ligplot (Figure 1). The interaction pattern between rutin and spike protein revealed the formation of five hydrogen bonds and four hydrophobic interactions at the binding site, thereby contributing to the stability of the complex. The residues involved in the formation of hydrogen bonds between rutin and spike protein were F970, N969, H49, Q52, and T274. While the hydrophobic contact between the two involved residues T51, S50, S967 and S968. The interactions between rutin and main protease was stabilized by residues E166, T190, Y54, D187, T26, N142 forming hydrogen bond and various hydrophobic contacts with residues P168, L167, H41, G143, T25, L27, M49, M165, H164, R188 and Q189 (Tables 1 and 2). A similar study exploiting a blind molecular docking approach utilizing the Swiss Dock server also highlights main protease inhibition by rutin where by common amino acids viz. E166, T190, N142, H41 are found to be interacting with the inhibitor. In our study however, the inhibitor displayed slightly better binding efficiency than their report (Das et al., 2020). From the contact analyses of both proteins, it can be attributed that rutin has high affinity and wide molecular contacts for both spike protein and main protease. Also, it was further scrutinized by molecular dynamics simulations to attain insights towards the inhibitory aspects and efficacy to combat spike protein and main protease drug targets of SARS-CoV-2.
|
|
figure
|
Figure 1. Protein-ligand interactions diagram of spike protein and main protease in complex with rutin. 2D Ligplot representation showing hydrogen bond and hydrophobic interactions, and 3D structure representation showing rutin in the binding pocket of spike protein (A and C) and main protease (B and D).
|
|
label
|
Figure 1.
|
|
caption
|
Protein-ligand interactions diagram of spike protein and main protease in complex with rutin. 2D Ligplot representation showing hydrogen bond and hydrophobic interactions, and 3D structure representation showing rutin in the binding pocket of spike protein (A and C) and main protease (B and D).
|
|
p
|
Protein-ligand interactions diagram of spike protein and main protease in complex with rutin. 2D Ligplot representation showing hydrogen bond and hydrophobic interactions, and 3D structure representation showing rutin in the binding pocket of spike protein (A and C) and main protease (B and D).
|
|
table-wrap
|
Table 1. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between spike protein with Rutin.
Spike protein + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
His 49.A NE2–Rutin.het O15 3.28 HNE2-H--O15 Ser50, Thr51, Ser967, Ser968
GLN 52.A OE1–Rutin.het O14 2.89 HOE1-H--O14
THR 274.A OG1–Rutin.het O14 3.30 HOG1-H--O14
ASN 969.A ND2–Rutin.het O10 3.17 HND2-H--O10
PHE 970.A N–Rutin.het O9 3.10 HN-H--O9
|
|
label
|
Table 1.
|
|
caption
|
Intermolecular hydrogen bonds and hydrophobic residues showing close contact between spike protein with Rutin.
|
|
p
|
Intermolecular hydrogen bonds and hydrophobic residues showing close contact between spike protein with Rutin.
|
|
table
|
Spike protein + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
His 49.A NE2–Rutin.het O15 3.28 HNE2-H--O15 Ser50, Thr51, Ser967, Ser968
GLN 52.A OE1–Rutin.het O14 2.89 HOE1-H--O14
THR 274.A OG1–Rutin.het O14 3.30 HOG1-H--O14
ASN 969.A ND2–Rutin.het O10 3.17 HND2-H--O10
PHE 970.A N–Rutin.het O9 3.10 HN-H--O9
|
|
tr
|
Spike protein + rutin
|
|
th
|
Spike protein + rutin
|
|
tr
|
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
|
|
th
|
Interacting residue
|
|
th
|
H-bond distance (Å)
|
|
th
|
H-bond (D-H--A)
|
|
th
|
Hydrophobic residues
|
|
tr
|
His 49.A NE2–Rutin.het O15 3.28 HNE2-H--O15 Ser50, Thr51, Ser967, Ser968
|
|
td
|
His 49.A NE2–Rutin.het O15
|
|
td
|
3.28
|
|
td
|
HNE2-H--O15
|
|
td
|
Ser50, Thr51, Ser967, Ser968
|
|
tr
|
GLN 52.A OE1–Rutin.het O14 2.89 HOE1-H--O14
|
|
td
|
GLN 52.A OE1–Rutin.het O14
|
|
td
|
2.89
|
|
td
|
HOE1-H--O14
|
|
td
|
|
|
tr
|
THR 274.A OG1–Rutin.het O14 3.30 HOG1-H--O14
|
|
td
|
THR 274.A OG1–Rutin.het O14
|
|
td
|
3.30
|
|
td
|
HOG1-H--O14
|
|
td
|
|
|
tr
|
ASN 969.A ND2–Rutin.het O10 3.17 HND2-H--O10
|
|
td
|
ASN 969.A ND2–Rutin.het O10
|
|
td
|
3.17
|
|
td
|
HND2-H--O10
|
|
td
|
|
|
tr
|
PHE 970.A N–Rutin.het O9 3.10 HN-H--O9
|
|
td
|
PHE 970.A N–Rutin.het O9
|
|
td
|
3.10
|
|
td
|
HN-H--O9
|
|
td
|
|
|
table-wrap
|
Table 2. Intermolecular hydrogen bonds and hydrophobic residues showing close contact between main protease with Rutin.
Main protease + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
THR 26.A O–Rutin.het O15 3.05 HO-H--O15 Thr25, Leu27, His41, Met49, Gly143, His164,
THR 26.A O–Rutin.het O16 3.02 HO-H--O16 Met165, Leu167, Pro168, Arg188,
TYR 54.A OH–Rutin.het O14 2.91 HOH-H--O14 Gln189
ASN 142.A OD1–Rutin.het O5 3.01 HOD1-H--O5
ASN 142.A OD1–Rutin.het O7 2.91 HOD1-H--O7
GLU 166.A O–Rutin.het O8 3.02 HO-H--O8
ASP 187.A O–Rutin.het O14 3.31 HO-H--O14
THR 190.A O–Rutin.het O10 2.87 HO-H--O10
|
|
label
|
Table 2.
|
|
caption
|
Intermolecular hydrogen bonds and hydrophobic residues showing close contact between main protease with Rutin.
|
|
p
|
Intermolecular hydrogen bonds and hydrophobic residues showing close contact between main protease with Rutin.
|
|
table
|
Main protease + rutin
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
THR 26.A O–Rutin.het O15 3.05 HO-H--O15 Thr25, Leu27, His41, Met49, Gly143, His164,
THR 26.A O–Rutin.het O16 3.02 HO-H--O16 Met165, Leu167, Pro168, Arg188,
TYR 54.A OH–Rutin.het O14 2.91 HOH-H--O14 Gln189
ASN 142.A OD1–Rutin.het O5 3.01 HOD1-H--O5
ASN 142.A OD1–Rutin.het O7 2.91 HOD1-H--O7
GLU 166.A O–Rutin.het O8 3.02 HO-H--O8
ASP 187.A O–Rutin.het O14 3.31 HO-H--O14
THR 190.A O–Rutin.het O10 2.87 HO-H--O10
|
|
tr
|
Main protease + rutin
|
|
th
|
Main protease + rutin
|
|
tr
|
Interacting residue H-bond distance (Å) H-bond (D-H--A) Hydrophobic residues
|
|
th
|
Interacting residue
|
|
th
|
H-bond distance (Å)
|
|
th
|
H-bond (D-H--A)
|
|
th
|
Hydrophobic residues
|
|
tr
|
THR 26.A O–Rutin.het O15 3.05 HO-H--O15 Thr25, Leu27, His41, Met49, Gly143, His164,
|
|
td
|
THR 26.A O–Rutin.het O15
|
|
td
|
3.05
|
|
td
|
HO-H--O15
|
|
td
|
Thr25, Leu27, His41, Met49, Gly143, His164,
|
|
tr
|
THR 26.A O–Rutin.het O16 3.02 HO-H--O16 Met165, Leu167, Pro168, Arg188,
|
|
td
|
THR 26.A O–Rutin.het O16
|
|
td
|
3.02
|
|
td
|
HO-H--O16
|
|
td
|
Met165, Leu167, Pro168, Arg188,
|
|
tr
|
TYR 54.A OH–Rutin.het O14 2.91 HOH-H--O14 Gln189
|
|
td
|
TYR 54.A OH–Rutin.het O14
|
|
td
|
2.91
|
|
td
|
HOH-H--O14
|
|
td
|
Gln189
|
|
tr
|
ASN 142.A OD1–Rutin.het O5 3.01 HOD1-H--O5
|
|
td
|
ASN 142.A OD1–Rutin.het O5
|
|
td
|
3.01
|
|
td
|
HOD1-H--O5
|
|
td
|
|
|
tr
|
ASN 142.A OD1–Rutin.het O7 2.91 HOD1-H--O7
|
|
td
|
ASN 142.A OD1–Rutin.het O7
|
|
td
|
2.91
|
|
td
|
HOD1-H--O7
|
|
td
|
|
|
tr
|
GLU 166.A O–Rutin.het O8 3.02 HO-H--O8
|
|
td
|
GLU 166.A O–Rutin.het O8
|
|
td
|
3.02
|
|
td
|
HO-H--O8
|
|
td
|
|
|
tr
|
ASP 187.A O–Rutin.het O14 3.31 HO-H--O14
|
|
td
|
ASP 187.A O–Rutin.het O14
|
|
td
|
3.31
|
|
td
|
HO-H--O14
|
|
td
|
|
|
tr
|
THR 190.A O–Rutin.het O10 2.87 HO-H--O10
|
|
td
|
THR 190.A O–Rutin.het O10
|
|
td
|
2.87
|
|
td
|
HO-H--O10
|
|
td
|
|
|
sec
|
3.3. Molecular dynamics simulation of spike protein and main protease in absence and presence of rutin
The binding affinity of rutin against both spike and main protease were investigated through MD simulations using GROMACS. Separate MD simulations for spike protein and main protease with and without rutin for 100 ns were assessed through trajectory analysis. For the course of 100 ns MD simulation, the stable trajectory was observed and the representative structures were obtained. The deviation of the backbone atoms for simulated structures relative to the starting structures used as a reference was evaluated through RMSD for the backbone atoms (Figures 2 and 3) and C-alpha atoms (supporting material Figure S1). A steep magnitude RMSD variation during the entire simulation can be an implication of a malleable and free instinctive protein or the alteration of the force field. Depending upon the outcomes of the RMSD evaluation, Figure 2(A) represents that the RMSD fluctuation stabilize at about 50 ns MD simulations for both spike protein and spike protein in presence of rutin, and the simulation time was acceptable. In the time from 50–100 ns, the RMSD for spike protein and spike-rutin complex have approximate values about 1.50–1.65 nm and 1.20–1.26 nm, respectively. Similarly, the main protease in absence and presence of rutin for time from 20–70 ns and 20–100 ns showed RMSD values of about 0.25–0.30 nm and 0.20–0.25 nm (Figure 3(A)). Interestingly, RMSD for both the complex systems were more stable than the native state of target receptor spike protein and main protease (Figures 2(A) and 3(A)). The RMSD of C-alpha atoms were also attained and was found to be similar to the RMSD of backbone atoms (supporting material Figure S1). The average RMSD for backbone and C-alpha atoms of both the protein and their complex with rutin are shown in supporting material Table S1.
Figure 2. Representation of the MD analysis plots for spike protein (purple color) and spike-rutin complex (red color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg), and (D) solvent accessible surface area (SASA) analysis.
Figure 3. Representation of the MD analysis plots for main protease (blue color) and main protease-rutin complex (green color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg) and (D) solvent accessible surface area (SASA) analysis. The integral extent of the residual dynamic criterion is achieved by assessing the variations arising from shifts of each of the protein residues that majorly feature the most flexible chain frames. Hence, we validated the residual fluctuations by calculating the mean fluctuation for stable trajectories of each simulation. The RMSF evaluation of all the protein residues were achieved in order to check the residues that may have tend towards an enhancement in the RMSD results. Compelling fluctuations existed in terminal residues and in various loops along β sheets and coils of the spike protein (residues 852-1000) to about 0.47 nm and spike-rutin (residues 62–90, 400–520 and 675–730) to about 0.65 nm (Figure 2(B)). Similarly, in case of main protease and itself in complex with rutin showed fluctuations in the terminal residues and in some regions of β sheet with some coils and loops in main protease (residues 136-148) to about 0.37 nm and main protease-rutin (residues 39-52, 150-156) to about 0.37 nm, respectively (Figure 3(B)). Also, we noticed that residues 270 to 350 in spike protein and 160 to 190 in main protease showed comparatively reduced discrepancy in RMSF values, which were the corresponding binding segment of both the complexes after MD simulation. Interestingly, the binding region for spike-rutin and main protease-rutin before MD simulation were found to be the same after simulation, suggesting as a strong binding interaction of the rutin towards spike protein and main protease and indicating stability of both the complexes.
The radius of gyration analysis was evaluated to determine the change in compactness of protein systems used throughout the MD simulations. The Rg plots for spike, main protease and its complex with rutin show slight fluctuations at the initial frame and attain compactness after 30 ns with Rg score of 3.85 nm and 4 nm (spike and spike-rutin) and after 20 ns with Rg score of 2.27 nm and 2.23 nm (main protease and main protease-rutin), respectively (Figures 2(C) and 3(C)). When compared with spike protein, Rg value for spike-rutin is stabilized and remains constant, suggesting strong binding interaction of the inhibitor, and the same is observed in case of main protease and rutin. Similar observations were determined through SASA analysis representing the solvent defined protein surface and its orientation through folding, making the alterations in the exposed and buried regions of the surface area of proteins. The SASA values for all the simulation systems were about 460 nm/S2/N (spike), 485 nm/S2/N (spike-rutin), 154 nm/S2/N (main-protease), and 152 nm/S2/N (main protease-rutin) (Figures 2(D) and 3(D)). Here, spike-rutin and main protease-rutin solvation profile shows a convincing SASA value suggesting a stable structure and strong binding interaction with the rutin.
Further, the cluster analysis with a RMSD based cut-off value of 0.25 nm demonstrated the development of 13,28,3 and 2 distinctive clusters for spike, spike-rutin, main protease, and main protease-rutin simulation systems respectively (Figure 4(A–D)). The most dominant cluster attained after 100 ns of MD simulation for all four systems are shown in Figure 4(E–H). Also, the secondary structure analyses of the stable trajectory for both the simulation systems were performed using the DSSP tool of GROMACS. Both, spike and main protease were formed mainly of conserved β-sheet along with some connecting loops and α-helix with small coil regions as secondary structure elements infused with various small segments of bend, turn, and β-bridge (Figure 5). Both cluster analysis and secondary structure analysis reveals the conformational changes before and after simulations for spike and main protease with and without rutin structures. A noticeable observation through these analyses supports the rationale that spike protein and main protease binds with rutin firmly and attain similar interactions after MD simulations also.
Figure 4. Representation of the MD analysis plots for the RMSD matrix cluster formation showing number of clusters in (A) spike protein, (B) spike-rutin complex, (C) main protease, and (D) main protease-rutin complex. The most dominant cluster conformation after MD cluster analysis of (E) spike protein, (F) spike-rutin complex, (G) main protease and (H) main protease-rutin complex.
Figure 5. The secondary structure analysis of (A) spike protein, (B) spike-rutin complex, (C) main protease and (D) main protease-rutin complex. Further, hydrogen bond analysis revealed a higher number of hydrogen bonds with the surrounding water molecules in case of spike protein (2300 hydrogen bonds) than in complex with rutin (approximately, 2200 hydrogen bonds). Similarly, in case of main protease, 700 hydrogen bonds were formed while 560 hydrogen bonds in case of main protease complex with rutin (Figure 6(A,B)). In addition, the hydrogen bond landscape assessed against the inhibitor revealed the dynamic equilibration of the complex trajectories with a high number of hydrogen bonds, as shown in Figure 6(C,D). The average number of hydrogen bonds within spike-rutin and main protease-rutin are 3.72 and 2.29, respectively. The consistent numbers of hydrogen bonds were observed which contributed significantly to the proximal binding of the potent drug rutin with the spike and main protease receptor. The superimposed structures of spike protein and main protease both in complex with rutin before and after simulations is shown in Figure 7(A,B). Further, these results were strengthened by the vital contribution of the complex binding energies throughout the simulation run. These calculations with consistent high binding energies and large hydrogen bonds involvement demonstrated the stable binding of rutin with both spike and main protease proteins.
Figure 6. Representation of the MD analysis plots for the hydrogen bond with the surrounding solvent against (A) spike protein (purple color) and spike-rutin complex (red color), and (B) main protease (blue) and main protease-rutin complex (green). The MD analysis plots representing hydrogen bonds with the potent inhibitor i.e. rutin against (C) main protease and (D) main protease-rutin complex.
Figure 7. Superimposed structures of (A) spike protein before simulation (green) and after simulation (blue), and (B) main protease before simulation (yellow) and after simulation (pink).
|
|
label
|
3.3.
|
|
title
|
Molecular dynamics simulation of spike protein and main protease in absence and presence of rutin
|
|
p
|
The binding affinity of rutin against both spike and main protease were investigated through MD simulations using GROMACS. Separate MD simulations for spike protein and main protease with and without rutin for 100 ns were assessed through trajectory analysis. For the course of 100 ns MD simulation, the stable trajectory was observed and the representative structures were obtained. The deviation of the backbone atoms for simulated structures relative to the starting structures used as a reference was evaluated through RMSD for the backbone atoms (Figures 2 and 3) and C-alpha atoms (supporting material Figure S1). A steep magnitude RMSD variation during the entire simulation can be an implication of a malleable and free instinctive protein or the alteration of the force field. Depending upon the outcomes of the RMSD evaluation, Figure 2(A) represents that the RMSD fluctuation stabilize at about 50 ns MD simulations for both spike protein and spike protein in presence of rutin, and the simulation time was acceptable. In the time from 50–100 ns, the RMSD for spike protein and spike-rutin complex have approximate values about 1.50–1.65 nm and 1.20–1.26 nm, respectively. Similarly, the main protease in absence and presence of rutin for time from 20–70 ns and 20–100 ns showed RMSD values of about 0.25–0.30 nm and 0.20–0.25 nm (Figure 3(A)). Interestingly, RMSD for both the complex systems were more stable than the native state of target receptor spike protein and main protease (Figures 2(A) and 3(A)). The RMSD of C-alpha atoms were also attained and was found to be similar to the RMSD of backbone atoms (supporting material Figure S1). The average RMSD for backbone and C-alpha atoms of both the protein and their complex with rutin are shown in supporting material Table S1.
|
|
figure
|
Figure 2. Representation of the MD analysis plots for spike protein (purple color) and spike-rutin complex (red color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg), and (D) solvent accessible surface area (SASA) analysis.
|
|
label
|
Figure 2.
|
|
caption
|
Representation of the MD analysis plots for spike protein (purple color) and spike-rutin complex (red color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg), and (D) solvent accessible surface area (SASA) analysis.
|
|
p
|
Representation of the MD analysis plots for spike protein (purple color) and spike-rutin complex (red color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg), and (D) solvent accessible surface area (SASA) analysis.
|
|
figure
|
Figure 3. Representation of the MD analysis plots for main protease (blue color) and main protease-rutin complex (green color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg) and (D) solvent accessible surface area (SASA) analysis.
|
|
label
|
Figure 3.
|
|
caption
|
Representation of the MD analysis plots for main protease (blue color) and main protease-rutin complex (green color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg) and (D) solvent accessible surface area (SASA) analysis.
|
|
p
|
Representation of the MD analysis plots for main protease (blue color) and main protease-rutin complex (green color) (A) root mean square deviation (RMSD) for backbone atoms, (B) root mean square fluctuation (RMSF), (C) radius of gyration (Rg) and (D) solvent accessible surface area (SASA) analysis.
|
|
p
|
The integral extent of the residual dynamic criterion is achieved by assessing the variations arising from shifts of each of the protein residues that majorly feature the most flexible chain frames. Hence, we validated the residual fluctuations by calculating the mean fluctuation for stable trajectories of each simulation. The RMSF evaluation of all the protein residues were achieved in order to check the residues that may have tend towards an enhancement in the RMSD results. Compelling fluctuations existed in terminal residues and in various loops along β sheets and coils of the spike protein (residues 852-1000) to about 0.47 nm and spike-rutin (residues 62–90, 400–520 and 675–730) to about 0.65 nm (Figure 2(B)). Similarly, in case of main protease and itself in complex with rutin showed fluctuations in the terminal residues and in some regions of β sheet with some coils and loops in main protease (residues 136-148) to about 0.37 nm and main protease-rutin (residues 39-52, 150-156) to about 0.37 nm, respectively (Figure 3(B)). Also, we noticed that residues 270 to 350 in spike protein and 160 to 190 in main protease showed comparatively reduced discrepancy in RMSF values, which were the corresponding binding segment of both the complexes after MD simulation. Interestingly, the binding region for spike-rutin and main protease-rutin before MD simulation were found to be the same after simulation, suggesting as a strong binding interaction of the rutin towards spike protein and main protease and indicating stability of both the complexes.
|
|
p
|
The radius of gyration analysis was evaluated to determine the change in compactness of protein systems used throughout the MD simulations. The Rg plots for spike, main protease and its complex with rutin show slight fluctuations at the initial frame and attain compactness after 30 ns with Rg score of 3.85 nm and 4 nm (spike and spike-rutin) and after 20 ns with Rg score of 2.27 nm and 2.23 nm (main protease and main protease-rutin), respectively (Figures 2(C) and 3(C)). When compared with spike protein, Rg value for spike-rutin is stabilized and remains constant, suggesting strong binding interaction of the inhibitor, and the same is observed in case of main protease and rutin. Similar observations were determined through SASA analysis representing the solvent defined protein surface and its orientation through folding, making the alterations in the exposed and buried regions of the surface area of proteins. The SASA values for all the simulation systems were about 460 nm/S2/N (spike), 485 nm/S2/N (spike-rutin), 154 nm/S2/N (main-protease), and 152 nm/S2/N (main protease-rutin) (Figures 2(D) and 3(D)). Here, spike-rutin and main protease-rutin solvation profile shows a convincing SASA value suggesting a stable structure and strong binding interaction with the rutin.
|
|
p
|
Further, the cluster analysis with a RMSD based cut-off value of 0.25 nm demonstrated the development of 13,28,3 and 2 distinctive clusters for spike, spike-rutin, main protease, and main protease-rutin simulation systems respectively (Figure 4(A–D)). The most dominant cluster attained after 100 ns of MD simulation for all four systems are shown in Figure 4(E–H). Also, the secondary structure analyses of the stable trajectory for both the simulation systems were performed using the DSSP tool of GROMACS. Both, spike and main protease were formed mainly of conserved β-sheet along with some connecting loops and α-helix with small coil regions as secondary structure elements infused with various small segments of bend, turn, and β-bridge (Figure 5). Both cluster analysis and secondary structure analysis reveals the conformational changes before and after simulations for spike and main protease with and without rutin structures. A noticeable observation through these analyses supports the rationale that spike protein and main protease binds with rutin firmly and attain similar interactions after MD simulations also.
|
|
figure
|
Figure 4. Representation of the MD analysis plots for the RMSD matrix cluster formation showing number of clusters in (A) spike protein, (B) spike-rutin complex, (C) main protease, and (D) main protease-rutin complex. The most dominant cluster conformation after MD cluster analysis of (E) spike protein, (F) spike-rutin complex, (G) main protease and (H) main protease-rutin complex.
|
|
label
|
Figure 4.
|
|
caption
|
Representation of the MD analysis plots for the RMSD matrix cluster formation showing number of clusters in (A) spike protein, (B) spike-rutin complex, (C) main protease, and (D) main protease-rutin complex. The most dominant cluster conformation after MD cluster analysis of (E) spike protein, (F) spike-rutin complex, (G) main protease and (H) main protease-rutin complex.
|
|
p
|
Representation of the MD analysis plots for the RMSD matrix cluster formation showing number of clusters in (A) spike protein, (B) spike-rutin complex, (C) main protease, and (D) main protease-rutin complex. The most dominant cluster conformation after MD cluster analysis of (E) spike protein, (F) spike-rutin complex, (G) main protease and (H) main protease-rutin complex.
|
|
figure
|
Figure 5. The secondary structure analysis of (A) spike protein, (B) spike-rutin complex, (C) main protease and (D) main protease-rutin complex.
|
|
label
|
Figure 5.
|
|
caption
|
The secondary structure analysis of (A) spike protein, (B) spike-rutin complex, (C) main protease and (D) main protease-rutin complex.
|
|
p
|
The secondary structure analysis of (A) spike protein, (B) spike-rutin complex, (C) main protease and (D) main protease-rutin complex.
|
|
p
|
Further, hydrogen bond analysis revealed a higher number of hydrogen bonds with the surrounding water molecules in case of spike protein (2300 hydrogen bonds) than in complex with rutin (approximately, 2200 hydrogen bonds). Similarly, in case of main protease, 700 hydrogen bonds were formed while 560 hydrogen bonds in case of main protease complex with rutin (Figure 6(A,B)). In addition, the hydrogen bond landscape assessed against the inhibitor revealed the dynamic equilibration of the complex trajectories with a high number of hydrogen bonds, as shown in Figure 6(C,D). The average number of hydrogen bonds within spike-rutin and main protease-rutin are 3.72 and 2.29, respectively. The consistent numbers of hydrogen bonds were observed which contributed significantly to the proximal binding of the potent drug rutin with the spike and main protease receptor. The superimposed structures of spike protein and main protease both in complex with rutin before and after simulations is shown in Figure 7(A,B). Further, these results were strengthened by the vital contribution of the complex binding energies throughout the simulation run. These calculations with consistent high binding energies and large hydrogen bonds involvement demonstrated the stable binding of rutin with both spike and main protease proteins.
|
|
figure
|
Figure 6. Representation of the MD analysis plots for the hydrogen bond with the surrounding solvent against (A) spike protein (purple color) and spike-rutin complex (red color), and (B) main protease (blue) and main protease-rutin complex (green). The MD analysis plots representing hydrogen bonds with the potent inhibitor i.e. rutin against (C) main protease and (D) main protease-rutin complex.
|
|
label
|
Figure 6.
|
|
caption
|
Representation of the MD analysis plots for the hydrogen bond with the surrounding solvent against (A) spike protein (purple color) and spike-rutin complex (red color), and (B) main protease (blue) and main protease-rutin complex (green). The MD analysis plots representing hydrogen bonds with the potent inhibitor i.e. rutin against (C) main protease and (D) main protease-rutin complex.
|
|
p
|
Representation of the MD analysis plots for the hydrogen bond with the surrounding solvent against (A) spike protein (purple color) and spike-rutin complex (red color), and (B) main protease (blue) and main protease-rutin complex (green). The MD analysis plots representing hydrogen bonds with the potent inhibitor i.e. rutin against (C) main protease and (D) main protease-rutin complex.
|
|
figure
|
Figure 7. Superimposed structures of (A) spike protein before simulation (green) and after simulation (blue), and (B) main protease before simulation (yellow) and after simulation (pink).
|
|
label
|
Figure 7.
|
|
caption
|
Superimposed structures of (A) spike protein before simulation (green) and after simulation (blue), and (B) main protease before simulation (yellow) and after simulation (pink).
|
|
p
|
Superimposed structures of (A) spike protein before simulation (green) and after simulation (blue), and (B) main protease before simulation (yellow) and after simulation (pink).
|
|
sec
|
4. Conclusion
SARS-CoV-2 which causes the deadly and highly infectious COVID-19 remains an unmet medical issue that requires urgent attention, owing to the number of deaths it is causing globally. We report here a computational approach for a repurposing screening program against two crucial SARS-CoV-2 drug targets viz. the spike protein and the main protease. We discovered that the drug rutin, showed a high binding potency against both the drug targets and may therefore be capable of eliciting reduction in viral load. The current study delineates that the highly effective repurposed drug rutin has the potential to obtain strong inhibition against SARS-CoV-2 spike protein and main protease. Rutin showed high binding efficiency against spike protein and main protease having an XP Glide score of –8.367 kcal/mol and –11.553 kcal/mol.
Further, the molecular interactions of rutin with both the drug targets were analysed thoroughly with the help of molecular dynamics simulation studies. The simulations for spike protein and main protease provided us with the stable trajectory analysis for both the targets. We analyzed various parameters which include RMSD, RMSF, radius of gyration, SASA, binding energies, hydrogen bond interactions, secondary structure contents, and cluster formation throughout the simulation, which proposed insights about the stable binding and compactness due to strong interaction of the rutin with both the drug targets. Through this in silico approach, the repurposed drug, rutin depicted reliable inhibition of two vital proteins and therefore may lead to elicitation of anti-SARS-CoV-2 activity.
|
|
label
|
4.
|
|
title
|
Conclusion
|
|
p
|
SARS-CoV-2 which causes the deadly and highly infectious COVID-19 remains an unmet medical issue that requires urgent attention, owing to the number of deaths it is causing globally. We report here a computational approach for a repurposing screening program against two crucial SARS-CoV-2 drug targets viz. the spike protein and the main protease. We discovered that the drug rutin, showed a high binding potency against both the drug targets and may therefore be capable of eliciting reduction in viral load. The current study delineates that the highly effective repurposed drug rutin has the potential to obtain strong inhibition against SARS-CoV-2 spike protein and main protease. Rutin showed high binding efficiency against spike protein and main protease having an XP Glide score of –8.367 kcal/mol and –11.553 kcal/mol.
|
|
p
|
Further, the molecular interactions of rutin with both the drug targets were analysed thoroughly with the help of molecular dynamics simulation studies. The simulations for spike protein and main protease provided us with the stable trajectory analysis for both the targets. We analyzed various parameters which include RMSD, RMSF, radius of gyration, SASA, binding energies, hydrogen bond interactions, secondary structure contents, and cluster formation throughout the simulation, which proposed insights about the stable binding and compactness due to strong interaction of the rutin with both the drug targets. Through this in silico approach, the repurposed drug, rutin depicted reliable inhibition of two vital proteins and therefore may lead to elicitation of anti-SARS-CoV-2 activity.
|
|
sec
|
Supplementary Material
Supporting_Information_Revised.pdf
Click here for additional data file.
|
|
title
|
Supplementary Material
|
|
caption
|
Supporting_Information_Revised.pdf
|
|
title
|
Supporting_Information_Revised.pdf
|
|
caption
|
Click here for additional data file.
|
|
p
|
Click here for additional data file.
|
|
back
|
Acknowledgements
AK is grateful to Jawaharlal Nehru University and Teri School of Advanced Studies for the management and other facilities. AK is contended to Indian Council of Medical Research (ICMR), India for the SRF position.
Disclosure statement
The authors declare no conflicting interests.
|
|
ack
|
Acknowledgements
AK is grateful to Jawaharlal Nehru University and Teri School of Advanced Studies for the management and other facilities. AK is contended to Indian Council of Medical Research (ICMR), India for the SRF position.
|
|
title
|
Acknowledgements
|
|
p
|
AK is grateful to Jawaharlal Nehru University and Teri School of Advanced Studies for the management and other facilities. AK is contended to Indian Council of Medical Research (ICMR), India for the SRF position.
|
|
sec
|
Disclosure statement
The authors declare no conflicting interests.
|
|
title
|
Disclosure statement
|
|
p
|
The authors declare no conflicting interests.
|
