Computational Prediction of the Molecular Mechanism of Statin Group of Drugs Against SARS-CoV-2 Pathogenesis

The COVID-19 pandemic had raised a severe health concern across the globe, as most countries are ghting with the second wave of SARS-CoV-2 infection. Although several vaccines had cleared the clinical trials and are being administered in different countries, the degree of protection varies due to the emergence of novel viral strains. Moreover, no antiviral drug against SARS-CoV-2 was reported to date. Recently published reports on statin therapy on COVID-19 patients indicated that statin therapy is associated with a better clinical outcome and a signicant reduction in mortality risk. using docking genome sequences in using the Swiss Model via The percentage was determined by the Ramachandran The MD simulation studies of free protein templates and protein-ligand complex were performed using NAMD (Nanoscale Molecular Dynamics) program; v 2.8.


Abstract Background
The COVID-19 pandemic had raised a severe health concern across the globe, as most countries are ghting with the second wave of SARS-CoV-2 infection. Although several vaccines had cleared the clinical trials and are being administered in different countries, the degree of protection varies due to the emergence of novel viral strains. Moreover, no antiviral drug against SARS-CoV-2 was reported to date. Recently published reports on statin therapy on COVID-19 patients indicated that statin therapy is associated with a better clinical outcome and a signi cant reduction in mortality risk.

Methods
Blind docking of the critical structural and functional proteins of SARS-CoV-2 like RNA-dependent RNA polymerase, Mprotease of 3-CL-Pro, Helicase, and the Spike proteins ( wild type and mutants) with the statin molecules were performed using the Schrodinger docking tool. Prevalent mutations in the spike protein were determined by analyzing the SARS-CoV-2 genome sequences deposited in GISAID, using the NextAlign aligner tool. Wild type and mutant spike proteins were modeled using the Swiss Model server and subjected to energy minimization via YASARA server. The percentage of the allowed region was determined by the Ramachandran plot. The MD simulation studies of free protein templates and protein-ligand complex were performed using NAMD (Nanoscale Molecular Dynamics) program; v 2.8.

Results
We observed that uvastatin and pitavastatin showed fair binding a nities to RNA polymerase, 3-CL-Pro, and spike-double mutant with similar docking scores. But of all the target proteins, uvastatin showed the strongest binding a nity to the helicase with the glide score, emodel score, and glide energy values of -11.333, -66.511 and -58.72 (kcal/mol), respectively.
Pitavastatin having similar chemical properties like uvastatin (logP and pKa values) exhibited strong biding a nities to RdRp, 3-CL-Pro, and S-double mutant. Additionally, molecular dynamics simulation con rmed the formation of a stable drugprotein complex.

Conclusion
Thus our study shows that of all the statins, uvastatin can bind to multiple target proteins of SARS-CoV-2, including the spike-mutant proteins. This property might contribute to the potent antiviral e cacy of this drug.
The ORFs 1a and 1b code for the proteases 3Cl-Pro (also known as Mpro or the Main protease), and PL-Pro respectively, that further cleave the polypeptide into 16 nonstructural proteins. These structural proteins play pivotal roles in the entry and assembly of the virus particles in the host cells [6,7]. The replication/transcription machinery of the virus is mediated by two enzymes, namely the RNA-dependent R.N.A. polymerase (RdRp) or nsp12 and helicase or nsp13. While RdRP mediates the viral replication, the helicase catalyzes the unwinding of double-stranded RNA formed during replication and allows the next round of viral replication [8,9]. Most of these structural and nonstructural proteins of SARS-COV-2 share highly conserved functional domains with its predecessor strain SARS-CoV-2(5) and hence may be targeted with the existing or novel antiviral agents (6).
Genomic analysis of the available SARS-CoV-2 sequences across the globe revealed that the viral genome had acquired a speci c mutation in S protein, which facilitated its entry and infectivity in host cells. Those variants were found to possess the D614G mutation, i.e. the replacement of the aspartate (D) with glycine (G) at the 614th amino acid of S protein, and became predominant across the globe outcompeting the wild type strain (7). Also, it was further revealed that the D614G mutation enhanced the infectivity of the virus in the host cells (8,9). With the progression of pandemic, SARS-CoV-2 acquired additional mutations in S protein on the background of D614G mutation in different geographical regions and classi ed as B.1.1.7or N501Y.V1 variant (the UK variant), B.1.351or N501Y.V2 variant (the South African variant), and P.1 or N501Y.V3 variant (the Brazilian variant). All these three variants possess the mutation in RBD of S protein, characterized by the replacement of asparagine (N) with tyrosine (Y) at position 501 of RBD, on the background of D614G mutation (10,11). The N501Y mutation has been reported to enhance the binding of viral S protein to the host ACE2 receptors and reduce the e cacy of the neutralizing antibodies targeting the RBD (12,13). This N501Y mutation hence is known as mutation of concern (MOC), and the new SARS-CoV-2 strains harboring this mutation and also additional mutations are termed as the variants of concern (VOC) (14). Thus these accumulating variations had enhanced the immune escape potential of new SARS-CoV-2 variants. Since the prevailing vaccines were designed against the wild-type SARS-CoV-2 discovered in 2019, concerns have been raised about whether these vaccines will be effective against the new VOCs. Hence the search for novel drug candidates which can target both the wild type as well as mutant proteins warrants investigation on utmost priority.
Presently there are no well-de ned effective therapies against COVID-19. The only medicine prescribed by the doctors is dexamethasone, which has been shown to reduce 28-day mortality in COVID-19 (15). Although Remdesivir, a novel nucleotide analogue, initially showed some promising results in lowering the oxygen requirement of hospitalized COVID-19 patients (16), later WHO has recommended against the use of remdesivir (17). Currently, a plethora of randomized trials are going on to investigate the possible therapeutic remedies against COVID-19, but no conclusive outcome has been found to date.
The bene cial effects of the statin group of drugs are well documented. The statins have shown a promising therapeutic role against various autoimmune in ammatory induced disorders such as multiple sclerosis, systemic lupus erythematosus, rheumatoid arthritis, etc. (18). They also possess anti-hyperlipidemic and cardioprotective properties (19,20). Moreover, statins have also been effective against several viral infections such as Avian in uenza, Zika virus, hepatitis C virus, H1N1 pandemic, and also the Ebola outbreak in West Africa (21)(22)(23)(24)(25). In an interesting case study published in August 2020, a positive association between statin usage and reduced mortality of COVID-19 patients was rst reported (26). It was observed that the antecedent use of statin in hospitalized COVID-19 patients is associated with the lowering of mortality. At a similar time frame, another retrospective case-control study with the COVID-19 patients hospitalized at Hubei province, China, reported that the overall mortality risk had signi cantly reduced in the patients with in-hospital use of statins, mainly atorvastatin (27). In a retrospective single-center study, Daniels et al. published that statin use before hospitalization was associated with a substantially lower risk of COVID-severity and was also associated with a faster recovery (28). Another multi-centric retrospective cohort study on COVID-19 positive old population living in Belgium revealed a statistically signi cant association between statin usage and reduced COVID-19-severity (29). In a very recent preprint, Gerold et al.
demonstrated the direct effect of statins on coronavirus infection in ex vivo conditions using the human lung cells (30). They observed that among all the statins, uvastatin signi cantly reduced the entry of the low pathogenic coronavirus 229E into the human respiratory cells by modulating the host gene expressions (30). Interestingly, they also observed that uvastatin treatment signi cantly decreased SARS-CoV-2 genome copy numbers. But none of the above-mentioned studies had demonstrated the direct effect of statins on SARS-CoV-2 target proteins. Hence the mechanism by which statins may inhibit viral pathogenesis remains inconclusive. In the present study by in silico molecular docking and molecular dynamics simulation analysis of the interactions of the statin group of drugs with the essential structural and functional proteins of SARS-CoV-2, we have aimed to predict a possible mechanism by which the statins may inhibit SARS-CoV-2 infection.

Materials And Methods
Selection of target proteins of SARS-CoV-2 and its sequence homology with other coronaviruses As discussed earlier, the spike (S)-protein, RNA dependent RNA polymerase (RdRp), 3-Chymotrypsin-like protease (3CL-Pro) or the main protease and helicase are the critical proteins that regulate the various stages of SARS-CoV-2 infection, including the entry and replication of the viral genome in the host cells [5-9, 26, 27]. Hence for this study, the crystalline structures of RdRp (PDB ID: 7BV2), 3CL-Pro (PDB ID: 6LU7), and helicase (PDB ID: 6ZSL) were obtained from the RCSC Protein Data Bank (https://www.rcsb.org). The S-protein and its two mutants, such as the single mutant having D614G, and the double mutant having D614G/ N501Y, were modeled and subsequently processed via energy minimization (Sect. 2.2).
The FASTA sequences of the target proteins related to SARS-CoV and coronaviruses from other species were derived from UniProt database (https://www.uniprot.org/) for the sequence homology analysis. Multiple sequence alignment studies were performed with Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) to get the sequence homology data and generate the phylogenetic tree.
Determination of the mutational landscape of SARS-CoV-2 genome SARS-CoV-2 genome sequences deposited to GISAID until April 2021 from the human host were obtained and processed further. After excluding the sequences having an inadequate length (< 25000 nt) and duplicated entries, a total of 1299807 sequences were considered for downstream analysis. Sequences were mapped against the SARS-CoV-2 reference sequence (Wuhan/WH01/2019) using NextAlign aligner tool (https://github.com/nextstrain/nextclade). Thereafter, nonsynonymous single nucleotide substitutions in each of the aligned sequences were detected using the Biopython tool. Additionally, to predict the time-dependent emergence of a speci c mutation, phylodynamic analysis of global subsamples of SARS-CoV-2 sequences (N = 3461) was performed using nextstrain/ncov pipeline (https://github.com/nextstrain/ncov). Entropy for each position on SARS-CoV-2 Spike protein, helicase, RdRp and main protease were calculated using a method encoded in the nextstrain/ncov pipeline (https://github.com/nextstrain/ncov).
Homology modeling, energy minimization, and validation of wild type and spike mutants The sequence of SARS-CoV-2 spike protein was downloaded from the NCBI protein database (Accession No: YP_009724390). We prepared three sequence les by changing the amino acid in the speci ed region, i.e. wild type (no alteration), single mutant (D614G), and double mutant (N501Y and D614G). The template sequence was identi ed via alignment of the available PDB sequences using the BLASTp program (NCBI ) for all three types of protein. Due to the unavailability of the crystal structure of the region of interest, we have taken the EM structure as a template (PDB ID: 7BBH) from the PDB database on the basis of query coverage, highest GMQE score, and percent (%) identity. The structures were then modeled using Swiss Model Server (31) and subjected to further analysis in SWISS-MODEL for the calculation of QMEAN Z-score, which includes cumulative Z-score of Cβ, atomic con guration, solvation, and torsion values. The predictions were nally validated by PROCHECK. The summation of the number of residues in favored regions and in additionally allowed regions was considered for percent (%) quality assessment using the Ramachandran plot in PROCHECK. For each of the models, energy minimization is done using YASARA protein minimization server (32). After the minimization of the predicted model, all the above parameters were checked. The percentage of the allowed region was given under the Ramachandran plot in Table 1.  Table S1. To make the statins suitable for docking, we performed the following steps such as salt removal, the addition of H-atoms, and deprotonation using the LigPrep v4.7 module of the Schrodinger suite. Epik v4.7 module was used for charge neutralization of the drug candidates for attaining the biological relevant pH (pH 7.0 ± 2). The high-energy tautomeric states were excluded to retain up to four stereoisomers, and generation of only one stereoisomer per ligand was allowed. The best drug candidate(s) for different target proteins was predicted by the virtual screening work ow module of Maestro (33).

Preparation of the target proteins
The target proteins obtained from PDB database or modelled in other platforms required further processing for being used for the docking analysis. The non-standard residues (residues other than amino acids) were rst removed from the target proteins using DS client, and the resulted protein structures were processed using the 'protein preparation wizard' of Maestro, Schrodinger (34). The target proteins were processed by the addition of H-atoms, assigning the bond order using CCD.
(Cambridge Crystallographic Datacenter) database, and introducing the missing disul de linkages. The structures were further optimized for H-bond assignment at pH 7.0 using PROPKA function. It was followed by a restrained energy minimization module by applying OPLS3e force eld to converge heavy atoms to an RMSD value of 0.30 Å (34).

Grid generation
The grid generation is a very crucial step for the analysis of drug-target protein interaction. In this process, a 3-D boundary for the ligand-binding site was generated in the target protein using Glide, version 8.2 of Maestro, Schrodinger (35). The top-ranked potential binding sites were identi ed by keeping at least 15 points per site for the generation of the standard grid.
The 'receptor grid generation' module was used to generate grids of 20 Å size across the binding site.

Molecular Docking
The docking of the statin molecules with the target proteins was performed in three steps, namely, high throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP). The processed ligands (statins) were rst screened for their binding a nities to the target proteins using the High-Throughput-Virtual Screening ( (41), and generalized Born Implicit solvent model was applied (42). The system was heated from 300 K to 310 K for 100 ps, and then equilibrated at 310 K for 100 ps, followed by the production stage for 2 ns, and nally cooling to 300 K for 100 ps. The equation of motion was discretized at a 2 fs interval. Here, Langevin dynamics model (43,44) was used for the proposed molecular dynamics simulation described above. Frames were collected at 0.5 ps interval. The trajectory les were analyzed with VMD.
(38) and VEGA ZZ (45) for generating the data of RMSD, RMSF, SASA, and radius of gyration plot. The radius of gyration is a signi cant factor for the structural integrity of the crystal structure of protein and protein-ligand complexes. Here, the student t-test was analyzed for statistical evaluation (times = 2fs; for 4600 frames), and the results were statistically signi cant (p < 0.0001). Solvent accessible surface area (SASA) is a measurement of the surface area of the ligand that is accessible to solvent molecules and also referred to as PSA (Polar solvent accessible area). SASA versus MD simulation time plots were used throughout the study to investigate changes in the position and conformation of the ligands within the protein complex corresponding through 4600 frames (for 2ns). Moreover, the gyration radius (rg) and solvent accessible surface area (SASA) for the free protein and statin-protein complexes were calculated. The radius of gyration describes the overall spread of the molecule and is de ned as the root mean square distance of the collection of atoms from their common centre of gravity.
SASA represents the solvent accessable surface area of the protein and determines wheather the amino acid residues are buried or exposed.

Evolutionary relationship of therapeutically targeted SARS-CoV-2 proteins with other viral strains
The evolutionary conservancy is a critical feature of a gene or protein which guides it to act as a potential drug target. Hence the determination of interspecies evolutionary distance is necessary to predict the effectiveness of the candidate protein to become a potential drug target. In this study, we investigated the interspecies divergence of four selected target proteins of SARS-CoV-2, namely the spike protein, RdRp, helicase, and main protease among the Coronavirinae family members (Fig.   S1). Protein sequence (FASTA) alignment of SARS-CoV-2 proteins with other coronavirus species from bat, civet, pangolin, including SARS-CoV were performed using P-BLAST. The analysis showed high sequence similarity for RdRp (> 95%), helicase (> 99%), and main protease (> 95%) proteins (Fig. S1). SARS-CoV-2 spike (S) protein showed maximum sequence homology with pangolin coronavirus (92.43%) and over 76% homology with coronaviruses from other species.
It was observed that the spike D614G mutation became predominant (~ 70%) since March 2020, just within two months of its origin, both P681H and N501Y started gaining frequency from December 2020 (~ 17-18%) and achieved 50% frequency in two months down the line. Both these mutations occurred in the background of the D614G mutation. Previous studies have already shown the functional impact of spike D614G and N501Y mutations that provided a selective advantage to these mutations in viral entry, either through the interaction of ACE2 or by S1-S2 cleavage (9,(46)(47)(48). Moreover, the N501Y mutation in S protein is present in more than 90% of sequences of B.1.1.7 variant (the UK variant), variant (the South African variant), and P.1 variant (the Brazilian variant) and exhibited stronger interaction with the host ACE-2 receptor (49). Hence it is also known as the MOC or mutation of concern.
Additionally, to remove biases due to oversampling of SARS-CoV-2 sequences from a speci c geographic region, a phylodynamic time tree analysis from a global subsample of 3461 SARS-CoV-2 sequences was performed ( Homology modeling, re nement, and structure validation of wild type and mutant S proteins The template sequence for modeling the wild type and spike mutants was identi ed through searching on PDB using the BLASTp program. The BLAST search against the PDB database provided homology sequences against several structures. The PDB 7BBH showed maximum query coverage (96%), highest GMQE score (0.71), and present identity (89%), along with the resolution (2.90A) of the structure. Hence 7BBH was selected as the template structure for the modeling of wild-type and mutant S proteins. Also, the GMQE score is generated by the target-template alignment of the modeled protein in SWISS-MODEL, and the value is expressed as a number between zero and one. Higher GMQE scores indicate increased structural reliability. The GMQE score for the modeled protein was found to be 0.98, which indicates good model accuracy. We prepared three sequence les by changing the amino acid in the speci ed region, i.e. wild type (no alteration), single mutant (Aspartate at 614th position altered to Glycine-S D614G , and double mutant (Aspartate at 614th position altered to Glycine-S D614G and Asparagine at 501st position changed to Tyrosine-S D614G/N501Y ). After the modeling and re nement, each structure was validated through the generation of the Ramachandran plot using the PROCHECK tool (Table 1). Ramachandran plot con rmed that most of the amino acids of the modeled structures lie in the allowed region (> 95%), which indicates the stability of the structures (Table 1) residues: 398-581, 628-687), a 'palm' domain (amino acid residues: 582-627, 688-815), and a 'thumb' domain (amino acid residues: 816-919) and also an additional nidovirus-unique N-terminal extension (amino acid residues: 1-397) [32]. As per PDBsum record, 7BV2 has seven beta-sheets, forty-six α -helices, fty-eight beta-turns, and two gamma-turns. 7BV2 comprises six active sites such as AC1, AC2, AC3, AC4, AC5, and AC6, respectively, of which AC3 and AC6 constitute the remdesivir binding site.
Molecular docking of 7BV2 with nine statin molecules revealed that only uvastatin, pitavastatin, pravastatin, rosuvastatin, and simvastatin quali ed with a docking score while the other statin molecules did not qualify the screening (Table 2). Among those statins, uvastatin and pitavastatin were the best candidate molecules with higher binding a nities (Fig. 3,  Fig. S3 and Fig. S4, Table 2). In contrast, the other statin molecules did not exhibit strong binding ( Table 2). The 3D structures of RdRp-ligand complexes and their binding sites were shown in Fig. 3 and Fig. S3-S4 (A&B). Both the molecules bind at the remdesivir binding site of the enzyme. The glide score, emodel score, and energy values for uvastatin were − 7.441, -60.23, and − 57.076 kcal/mol (Table 2), and the corresponding values for pitavastatin were − 7.5, -62.341, and − 58.654 (kcal/mol), respectively ( Table 2). In our study, we observed that the glide score, emodel score, and energy values for remdesivir were − 7.312, -56.672, and − 71.38 kcal/mol, respectively. The binding pockets for uvastatin and pitavastatin were further analyzed using the Discovery Studio client. The uvastatin binding pocket was found to be hydrophilic in nature (Fig. S3C) (Fig. S3D). The presence of both the H-bond donor residues (ASN 496 and ARG 569) and H-bond acceptor residues (GLY 683) has made this site favorable for the formation of H-bond with the ligand ( uvastatin). These residues, present in 32nd and 47th β-turn formed conventional H-bonds with the carboxyl and hydroxyl groups of uvastatin (Fig. 3B). Moreover, the alkyl group of uvastatin had hydrophobic interactions with the VAL 557 residue (Fig. 3B). All these interactions are favorable for forming a stable uvastatin-RdRP complex.  (Fig. S4C). The presence of basic amino acids imparts a slight basic character to the binding pocket of pitavastatin to RdRp (Fig. S4D).  Table 2). Thus the data obtained from the docking analysis indicated the strong a nity of uvastatin and pitavastatin towards the enzyme.
Characterization of the binding cavity of uvastatin revealed that it is partly hydrophilic (Fig. S5C, Blue color). The binding cavity is comprised of the residues like THR 25, THR 26, TYR 54, ASP 187, and GLN 189, which impart a polar (hydrophilic) character to the binding pocket (Fig. S5C). As per ionizibility, the binding site mainly was neutral due to the presence of the residues like THR 26, HIS 41, TYR 54, and CYS 145 (fade white color, Fig. S5D Moreover, the alkyl group of uvastatin had hydrophobic interactions with CYS 145 and HIS 163, present on the 12th β strand of ß sheet and 17th β turn, respectively (Fig. 4B). All these interactions contributed to a stable uvastatin-3CL-Pro complex.
Due to the similar LogP values (Table S1) the binding pocket for pitavastatin is similar to that of uvastatin. The binding site of pitavastatin was slightly basic due to the occurrence of the basic amino acid residues like HIS 41, HIS 163, HIS 164, HIS 141, and ARG 188 (Fig. S 6D). HIS 41 acts as H-donor while GLU 166 and ARG 188 act as H-acceptor and participate in the formation of conventional H-bond with hydroxyl groups of pitavastatin (Fig. S6E). The electron cloud of the uorophenyl group of pitavastatin contributed to π-π stacked interaction with LEU 141, present in 16th β turn of the protein (Fig. 4D Molecular docking of 6ZSL with the statins revealed that uvastatin exhibited the strongest binding a nity with the glide score, emodel, and glide energy values of -11.333, -66.511 and − 58.72 (kcal/mol), respectively ( Table 2). The 3D structure of the uvastatin-helicase complex is shown in Fig. 5A and 5B, while the major 2D interactions are shown in Fig. 5C. additional H-bond with the -OH group attached to the aliphatic chain. Moreover, ten hydrophobic interactions, such as alkylalkyl, pi-alkyl, pi-sigma, amide-pi stacked, were observed between the π electron-rich aromatic rings of uvastatin and the surrounding amino acid residues like GLY 538, ALA 312, ALA 313, and ALA 316. All these interactions had synergistic effects to make helicase the preferable target protein for uvastatin. Interestingly, the previous study had revealed that the ATP binding site of helicase is comprised of amino acid residues like K288, E375, Q404, R443, and R567 (55). In our study, we have observed that uvastatin binds to helicase in a similar region, and the residues like LYS 288 and ARG 443 play important roles in this interaction (Fig. 5B & Fig. S7). Hence it may be concluded that uvastatin may interfere with the ATP binding site of helicase and inhibit the activity of the enzyme.
Wild type (Spike wt ) and mutant spike proteins (Spike D614G and spike D614G/N501Y ) Spike protein is one of the essential proteins required for the entry of SARS-CoV2 into the host cell. The viral entry is mediated by the interaction of the RBD (receptor binding domain) of S1 subunit of spike with the ACE2 receptor of the host cell membrane is mainly responsible for the entry of SARS-CoV2 into the host cell. The structural insights of the different conformational states of the S-protein and S1(RBD)-ACE2 complex have been reported (56). Further S1-subunit is composed of four domains, namely the N-terminal domain (NTD), RBD, and two C-terminal domains (CTDs), CTD-1, and CTD-2 [45]. In our modeled S-proteins (wild type S protein and the mutants), the amino acids 14 to 685 constitute the S1 region, while the residues 686 to 1273 comprise the S2 site, which corroborates with the published report (57).
Molecular docking of the wild-type S-protein with the statins revealed that uvastatin exhibited the highest binding a nity to the S-protein while the other statins did not show strong binding a nity (Table 2). Fluvastatin binds with Spike wt protein at NTD. (N-Terminal Domain) region (the 3D structure is shown in Fig. 6A and 2D interaction diagram is shown in Fig. 6B), with the glide score and glide energy were − 5.854 and − 38.632 Kcal/mole, respectively ( Table 2). The binding cavity of uvastatin on S-protein was composed of hydrophilic residues like SER 98, THR 259, HIS 245, and LYS 97; and hydrophobic amino acid residues like TYR 144, TRP 258, ALA 260, and ALA 262 (Fig. S8C). The presence of HIS 245, TYR 144, and LYS 97 shifted the pH of the binding pocket to a slightly basic range that is suitable for the interaction of ionized state of uvastatin (Fig. S8D). The result of the molecular docking study of uvastatin to single mutant spike D614G was shown in Fig. 6C-6D & Fig. S9. The docking analysis revealed that the binding site of uvastatin on spike D614G was altered than that on the wild-type protein.
Additionally, the polar amino acids like ASN 955, THR 732, HIS 1058; the basic amino acids like LYS 854 and LYS 835; and the acidic amino acid ASP 830 were also present, contributing to the hydrophilicity of the binding pocket. The pH was favorable for docking of uvastatin molecules at a negatively charged state (Fig. S9D) Thus, the docking site at Spike D614G is more favorable for a higher number of noncovalent interactions compared to Spike wt. Therefore, the binding a nity was signi cantly higher, as re ected with the calculated glide score (-7.378 kcal/mol) and glide energy (-47.041 kcal/mol) ( Table 2).
The molecular docking study of uvastatin to double mutant spike spike D614G/N501Y protein was shown in Fig. 6E-6F & Fig.   S10. Analysis of the binding site revealed that uvastatin binds at the S2 region of spike D614G/N501Y protein. The binding pocket was found to be rich in hydrophobic amino acid residues like ILE 850, VAL 860, PRO 862, PRO 863, LEU 959, ALA 956, VAL 952, LEU 828, and ALA 829, which imparts a hydrophobic character to the binding site (Fig. S10C). But polar amino acids like ASN 955, THR 732, HIE 1058, GLN 853, THR 859; basic acid amino acids like LYS 854, LYS 835; and acidic amino acid ASP 830 contributed to the hydrophilicity of the binding pocket. The pH was suitable for docking of uvastatin in a negatively charged state (Fig. S10D). The O-atom of the carboxylic group (-COOH) of uvastatin formed a salt bridge with LYS854. The presence of both anionic and cationic amino acid residues made the pocket rich with H-bond acceptor and Hbond donor points (Fig. S10E). The carbonyl O-atom formed H-bond with LYS 854. Additionally, the aliphatic hydroxyl groups formed H-bond with VAL 860 and ASP 830. Moreover, the F-tom of uvastatin developed an H-bond with GLN 853. It also interacted with ALA 956. Additionally, the pose of uvastatin was suitable to form hydrophobic interactions like alkyl-alkyl, pi-alkyl between the electron cloud of aromatic rings and ILE 850, ALA 829, and LEU 828. All these interactions had synergistic effects of increasing the a nity of uvastatin to spike D614G/N501Y ( glide score − 7.8 kcal/mol, and glide energy − 45.217 kcal/mole) in comparison to Spike wt and Spike D614G (Table 2).

Molecular dynamics simulation of Fluvastain-target protein complexes
To validate the conformational changes of the target proteins in the free state or in complex with the statins, we performed two nanoseconds of MD simulation using NAMD (Nanoscale Molecular Dynamics) program; v 2.8. The parameters such as RMSD, RMSF, SASA, and radius of gyrations were analyzed throughout the simulation trajectory to determine the stability of the drug-protein complex. The docking analysis demonstrated that among all the nine statin molecules, uvastatin showed a fair, binding a nity to the target proteins, and the highest a nity was observed for the helicase protein. Hence we have performed the MD simulation for uvastatin-target protein complexes.

Helicase-Fluvastatin complex
The molecular dynamic simulation analyses of the uvastatin-helicase complex were shown in Fig. 7A-7D and summarized in Table 3. The RMSD plot (Fig. 7A) . 7C & 7D). Moreover, upon complexation with uvastatin, the average RMSF value helicase was reduced from 1.096 ± 0.0113 Å to 1.031 ± 0.0128 Å respectively, which indicates the formation of a thermodynamically favourable complex.
Moreover, The SASA value was reduced from 43091 ± 2.887 Å 2 for the free protein to 38259 ± 2.827 Å 2 in the uvastatinhelicase complex without any major change in the gyration radius (r g ) (Fig. S11A, S11B). These results indicate a strong interaction of uvastatin with the helicase without any unfolding of the protein structure.

RdRp-Fluvastatin complex
The MD simulation analysis of Fluvastatin-RdRp interaction revealed that the RMSD value of the ligand-protein complex had decreased signi cantly compared to the free protein ( Fig. 7E-7F, Table 3). The free RdRP and uvastatin-RdRp complex had RMSD values of 3.713 ± 0.006144 and 2.584 ± 0.00357, respectively (Fig. 7E). This signi cant reduction in RMSD is in favor of better stability of the uvastatin-RdRp complex. The RMSD difference was 0.9302 Å (< 2 Å), suggesting the stability of uvastatin at the active site of RdRp. Furthermore, the RMSF value was decreased from 0.9124 ± 0.01 Å (free RdRp) to 0.841 ± 0.01 Å due to complexation with uvastatin. The major regions of uctuations were highlighted with the color ballot ( Fig. 7G-7H). Also, the SASA value of the uvastatin-RDRP complex (33766 ± 1.987 Å 2 ) was signi cantly reduced in comparison to free RDRP (38565 ± 2.269 Å 2 ) without any alteration in the radius of gyrations (r g ) values (indicating NO unfolding of the complex) (Fig. S11C-S11D). Thus, the docking of uvastatin is strong enough so that water can not replace it.

3Cl-Pro-Fluvastatin complex
The MD analysis of uvastatin-3CL-Pro complex showed the reduction in the average RMSD values of uvastatin-3CL-Pro complex (2.576 ± 0.004 Å) in comparison to the free 3CL-Pro (2.888 ± 0.005Å) (Fig. 7I-7J, Table 3). The conformations of uvastatin at the beginning and end of MD simulation were superimposed in Fig. 7J, and the difference in RMSD value was 1.78 Å (< 2 Å), suggesting the stability of the uvastatin-protein complex. The average RMSF value of the amino acid residues was reduced from 0.845 ± 0.017 Å (free 3 CL-Pro) to 0.7794 ± 0.016 Å due to complexation with uvastatin. The major regions of uctuation were highlighted with a color ballot (Fig. 7K-7L). Moreover, there were no signi cant changes in the SASA values and r g values ( Fig. S11E-S11F, Table 3). These results indicate the formation of a stable complex between uvastatin and 3 CL-Pro.
S WT -Fluvastatin complex, S D614G -Fluvastatin complex and S D614G/N501Y -Fluvastatin complex As per the RMSD plot, complexation of Spike wt with uvastatin increased the average RMSD value from 5.95 ± 1.38 Å (free protein) to 6.25 ± 1.77 Å (complex) (Fig. 8A, Table 3 As per the molecular docking study, the uvastatin binds at the NTD. region (amino acid residue 240 to 250). The RMSF at this region (labeled as G1 in Fig. 8C) was reduced due to Fluvastatin binding, indicating the formation of a stable complex. The SASA plot and the gyration of radius (r g ) did not show any signi cant change, indicating the formation of a stable complex of uvastatin with spike wt protein ( Fig. S12A-S12B).
The comparative RMSD plot of free Spike D614G and uvastatin-Spike D614G complex is shown in Fig. 8E. Up to 1100 ns, the free protein and uvastatin-protein complex had very close RMSD values. After that, the RMSD increased gradually and reached a maximum value of 8.72 Å. Thus, there was 1.45 Å (< 2 Å) enhancement of RMSD at the end of the simulation period. If we consider the poses of uvastatin at the beginning and end of the simulation period (Fig. 8F), the RMSD value was 1.92 Å (< 2 Å). Again, there was no signi cant difference in the average RMSD value of free protein (5.97 ± 1.13 Å) and uvastatin-protein bound complex (6.22 ± 1.68 Å). Thus, it is evident that uvastatin-Spike D614G is stable, which is further supported by the fact that the average RMSF value after complexation was nearly the same (free protein: 0.98 ± 0.46 Å, uvastatin bound complex: 0.98 ± 0.44 Å). But major reductions in RMSF values of a few amino acid residues (C1, and C2 regions) were observed (Fig. 8G-8H). C2 comprises 800-850 amino acid residues that constitute the uvastatin binding (S2).
No signi cant changes in The SASA value and gyration of radius (r g ) were observed (Fig. S12C-S12D), which suggested forming the stable complex of uvastatin with the Spike D614G .
The molecular dynamics studies of uvastatin-Spike D614G/N501Y complex showed no signi cant deviation in the RMSD plot for the free protein or the ligand-protein complex (Fig. 8I-8J). The free protein had an RMSD value of 6.284 ± 0.027 Å, whereas the uvastatin bound complex had an RMSD value of 6.369 ± 0.025 Å. The curves nearly overlap with each other (Fig. 8I). The RMSD value of poses of uvastatin at the beginning and end of the simulation period (Fig. 8J) was 1.63 Å (< 2 Å). The average RMSF value of uvastatin-Spike D614G/N501Y complex (2.575 ± 0.022 Å) was less than that of free protein (2.788 ± 0.023 Å), suggesting better stability of the complex. The four regions of Spike D614G/N501Y, having reduced RMSF, were C1: 350-550, C2: 700-780, C3: 820-900, and C4: 980 to 1100, respectively, as highlighted in Fig. 8K-8L. The region C3 constitutes the uvastatin binding pocket. Hence, the binding of uvastatin at the docking site of Spike D614G/N501Y is stable.
This 2% increase in the r g also suggested the opening up of loops in the C1 region. This C1 region belongs to the RBD section of the spike protein and is involved in the interaction with the ACE2 receptor. Thus the conformational change in this region of the S-double mutant due to Fluvastatin binding may also interfere with spike-ACE-2 interactions in variant SARS-CoV-2 strains.

Discussions
Since the emergence of the COVID-19 pandemic, there have been many efforts to repurpose the existing drugs, already known to be bene cial to the patients to treat this disease. Several case-control studies showed a positive association between statin usage and reduced mortality of COVID-19 patients hospitalized with this disease (8, 10, 58). Also, it was hypothesized that since statins can up-regulate ACE2 expression (59), they may prevent coronavirus infection. The direct effect of statins on SARS-CoV-2 was rst demonstrated by Gerold et al. (60). They observed that selected statins, especially uvastatin, signi cantly reduced the entry of SARS-CoV-2 into the human respiratory cells and genome copy numbers of SARS-CoV-2 in the infected cells. But the mechanism by which statins can inhibit the SARS-Cov-2 progression is not well understood.
As most of the countries are facing the second or third wave of the pandemic, various vaccines have been commercialized for the people to confer protection against SARS-Cov-2. But, the emergence of VOCs acquiring unique mutations in essential viral genes had raised concern about the effectiveness of the vaccines against these novel variants. In a very recent study, it was reported that the neutralizing antibodies from the individuals, those who received one or two doses of either BNT162b2 or mRNA-1273 vaccines, showed limited e cacy against the pseudoviruses representing the globally predominant VOCs (61). This observation was supported by other reports revealing the limited effectiveness of the mRNA vaccines or convalescent plasma against the circulating variants(62, 63). Although there are parallel efforts of improving the vaccine e ciency against these variants, there is also the necessity of antiviral medications that can target these variants, harboring the concerned mutations.
In our study, nine statin molecules were screened against four selected target proteins of SARS-CoV2 and the mutated S protein with D614G and N501Y (double mutant) by in silico molecular docking and molecular dynamics study. As per the molecular docking studies, the details of the best-selected candidate drug molecules are tabulated in Table 2. The uvastatin exhibited a good binding a nity to all the selected target proteins. It binds with the active site of RdRp (remdesivir binding site), 3 CL-Pro (inhibition site), and helicase (ATP binding site). Along with uvastatin, pitavastatin had shown an equivalent binding a nity for RdRp, and 3-CL-Pro ( Table 2). But of all the target proteins, the interaction of uvastatin with helicase was the best, as a docking score of -11.3 kcal/mole was observed. Analysis of the binding site revealed that uvastatin binds to Chemical homology, thermodynamic parameters, polarity, and favorable interactions may lead to multiple target sites for a single drug molecule(66, 67). Analysis of the binding sites of uvastatin in all the target proteins (Table S2), revealed that amino acids like threonine, serine, asparagine, arginine, lysine, and aspartic acid residues were common. Pitavastatin having similar chemical properties like uvastatin (logP and pKa values) exhibited strong biding a nities to RdRp, 3-CL-Pro, and Sdouble mutant.

Conclusion
Although the major limitation of this work is that it is based on a computational prediction without any laboratory validation, but such studies are signi cant for the research groups performing wet-lab experiments intending to identify novel antiviral compounds. The most encouraging result obtained from our study is that uvastatin emerged as the best drug candidate from the blind docking studies with all the nine statin molecules (Fig. 9), which corroborates with the recently published functional study (30). Thus our research will help to predict the molecular mechanism by which this drug inhibits SARS-CoV-2 pathogenesis.