Interaction of surface glycoprotein of SARS-CoV-2 variants of concern with potential drug candidates: A molecular docking study

Background: COVID-19 has become a global threat. Since its first outbreak from Wuhan, China in December 2019, the SARS-CoV-2 virus has gone through structural changes arising due to mutations in its surface glycoprotein. These mutations have led to the emergence of different genetic variants threatening public health due to increased transmission and virulence. As new drug development is a long process, repurposing existing antiviral drugs with potential activity against SARS-CoV-2 might be a possible solution to mitigate the current situation. Methods: This study focused on utilizing molecular docking to determine the effect of potential drugs on several variants of concern (VOCs). The effect of various drugs such as baricitinib, favipiravir, lopinavir, remdesivir and dexamethasone, which might have the potential to treat SARS-CoV-2 infections as evident from previous studies, was investigated for different VOCs. Results: Remdesivir showed promising results for B.1.351 variant (binding energy: -7.3 kcal/mol) with residues Gln319 and Val503 facilitating strong binding. Favipiravir showed favorable results against B.1.1.7 (binding energy: -5.6 kcal/mol), B.1.351 (binding energy: -5.1 kcal/mol) and B.1.617.2 (binding energy: -5 kcal/mol). Molecular dynamics simulation for favipiravir/B.1.1.7 was conducted and showed significant results in agreement with our findings. Conclusions: From structural modeling and molecular docking experiments, it is evident that mutations outside the receptor binding domain of surface glycoprotein do not have a sharp impact on drug binding affinity. Thus, the potential use of these drugs should be explored further for their antiviral effect against SARS-CoV-2 VOCs.


Introduction
The spread of human coronavirus SARS-CoV-2 has been increasing since it was first detected in the Chinese city of Wuhan in December 2019. 1 Several efforts have been taken to prevent its spread after the World Health Organization (WHO) declared it a public health emergency on January 31, 2020. However, its continual spread across the world compelled the WHO to declare it a pandemic. 2 Different genetic variants of this novel coronavirus have appeared and been transmitted across the world amidst the pandemic. 3 WHO classified some of these genetic variants into three categories: variant of interest (VOI), variant of concern (VOC) and variants of high consequence (VOHC). 4 Five different genetic variants have been placed in the VOC category, due to their increased transmission rates, more severe disease, or significant reduction in antibodies generated due to previous infection or vaccination namely by B.1.1.7 (Alpha), B.1.351 (Beta), B.1.617.2 (Delta), AY.1 (Delta plus) and P.1 (Gamma). Certain anti-retroviral drugs have been used owing to their promising results for the emergency treatment of COVID-19 patients. 5 Remdesivir was the first drug to be approved by the United States Food and Drug Administration (USFDA) for the treatment of COVID-19 patients. 6 However, emerging mutations in drug targets (such as the receptor binding domain [RBD] of the surface glycoprotein) are likely to affect the binding affinity through altered drug-receptor interaction. 7 Considering the emergence of several VOCs in different parts of the world, it is important to ascertain the effect of their signature genomic variants such as single nucleotide polymorphisms (SNPs) and insertions/deletions (InDels). To the best of our knowledge, there is a paucity of such information, especially regarding the newly emerged Delta variant. 8 Docking studies are very helpful and serve as the first starting point for such investigations. 9 Therefore, in this study, surface glycoprotein sequences of different VOCs were modeled in silico and their interactions with the drugs baricitinib, dexamethasone, favipiravir, lopinavir and remdesivir were studied using molecular docking. These drugs have shown promising results in various clinical studies and thus have been considered to determine their binding affinity on SARS-CoV-2 VOCs. [10][11][12] The main aim of this study was to utilize an in silico docking approach to estimate relative changes in the binding affinity of these potential drugs against the characteristic mutational profile of the spike protein sequence in different VOCs, to predict their potential therapeutic efficacy against VOCs.

Spike mutations
The surface glycoprotein (Spike) allows the virus to bind to hACE2 receptors and thereby promotes the virus's entry into the host cell. 13,14 It is divided into two subunits, S1 and S2. The S1 subunit consists of the RBD which directly binds to hACE2. It is also the target of neutralizing antibodies. Thus, it is the region with most mutations with clinical significance in terms of viral transmissibility and virulence. 15,16 Major mutations reported in different VOCs are shown in Table 1.

Dataset collection and mutation analysis
A total of 24 full-length sequences of SARS-CoV-2 genomes categorized into five VOCs from different geographical regions were selected and retrieved from the Global Initiative on Sharing All Influenza Data (GISAID) database. 18,19 The first sequence of SARS-CoV-2 originating from Wuhan was retrieved from the National Center for Biotechnology Information (NCBI) nucleotide database as a reference (NCBI reference sequence: NC_045512.2). Mutation analysis was carried out by multiple sequence alignment of the retrieved sequences using the ClustalW algorithm in MEGA-X software v 10.2.3 20 and mutation positions were determined. This analysis was performed to check the frequency of mutations across sequences from different VOCs. After this, a particular mutation was inserted in the sequence of a VOC, followed by its modeling. As these sequences were derived from COVID-19 infected patients, they represent the actual frequency of genetic mutations acquired by SARS-CoV-2. It reveals that once a specific mutation in any variant has evolved, it remains conserved in the descendant population, which may again acquire new characteristic mutations.
3D structure prediction, model quality assessment and validation The novel SARS-CoV-2 surface glycoprotein nucleotide wild type (WT) gene sequence NC_045512.2 was retrieved from the NCBI 21 nucleotide database. Reported mutations were induced in the retrieved sequence. A homology model was built for the surface glycoprotein of SARS-CoV-2 VOCs using SWISS-MODEL software. 22  The modeled surface glycoprotein of all SARS-CoV-2 VOCs was served as binding target and five approved drugs as ligand. All the compounds were first optimized in their active forms in physiological conditions.

Protein and ligand preparation
The structure of investigated drugs, namely baricitinib (PubChem CID 44205240), favipiravir (PubChem CID 492405), lopinavir (PubChem CID 92727), remdesivir (PubChem CID 121304016) and dexamethasone (PubChem CID 5743) were retrieved from PubChem database. 34 AutoDockTools 1.5.6, a free graphical user interface of MGL software package was used for all the required file conversions needed for the docking study. 35,36 The rotatable bonds present on the ligands were treated as non-rotatable for performing the docking. All the water molecules and hetero atoms present on the receptor surface were removed, followed by the addition of Kollman charges and polar hydrogen atoms using AutoDockTools 1.5.6. The Gasteiger charge calculation method and partial charges were also applied to the ligand molecules. 37 Grid box preparation and docking Molecular Docking was performed with modeled surface glycoproteins of different VOCs as receptors and selected drugs as ligands. Grid box parameters were selected using AutoDockTools 1.5.6 ( Table 2). The Lamarckian Genetic Algorithm was used for performing docking to explore the conformational space required for the ligand with a population size of 150 individuals. The total number of current grid points per map was 64,000. Other parameters were set at default.

Molecular dynamics (MD) simulation
To check the validity of molecular docking results for favipiravir against the B.1.1.7 (Alpha) variant, a molecular dynamics simulation was conducted. The simulation was conducted in GROMACS 2018, 38,39 with CHARMM36 as all-atom force field. 40 For the ligand-receptor complex, all receptors missing hydrogen atoms were added using Table 1. Spike mutations reported in variants of concern VOCs. Mutations in bold represent its presence in different variants.
Variant of concern (VOCs) H655Y T1027I V1176F *AY.1 is commonly referred as Delta plus, 17 although this is not as per WHO classification which considers it as one of the types within Delta lineage.
Chimera. 41,42 The protein-ligand complex was placed in an isotonic box with a dodecahedron cell. The box contained a neutralizing number of sodium (Na + ) and chloride (Cl -) ions based on the total charge of the protein. Topology parameters for the ligand were built using the CHARMM General Force Field (CGenFF) tool 43 to generate the CHARMM36 parameters. The solvation step was followed by energy minimization, equilibration number of particles volume temperature (NVT) ensemble and number of particles pressure temperature (NPT) ensemble; then, MD simulation with 2 femtoseconds (fs) integration steps for 20 ns were conducted. The output trajectory was then subjected to Periodic Boundary Conditions (PBC) correction and the system was fitted to its start position based on the backbone of the receptor. Further analysis was performed to plot the root-mean-square deviation (RMSD) of the ligand and Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) energy computation.

Results and discussion
Multiple sequence alignment and structure analysis The mutation analysis was carried out for the surface glycoprotein of SARS-CoV-2 VOC, and random clinical samples from different geographical regions were retrieved from the GISAID database 18  . However, no novel mutation was reported. The reported mutations which had a frequency of occurrence of 50% and above were considered for surface glycoprotein modeling for VOCs (Table 1). A high-quality model was constructed for the surface glycoprotein of SARS-CoV-2 based on the matching templates (7KRS, 7N1U and 7N1Q) for different SARS-CoV-2 VOCs. The 7KRS PDB accession number is a viral protein complex characterized by the mutation D614G solved by electron microscopy with a resolution of 3.20 Å. 45   residues Asn315, Arg317 and His617, one carbon-hydrogen bond with Gln319, and one alkyl bond with Ala290, Val503 and Ala621. Two amino acid residues, i.e., Gln319 and Val503 are present in the RBD of SARS-CoV-2 which facilitate its binding to the hACE2. 50 Similarly, baricitinib showed a strong binding affinity of 105.97% (binding energy: The emergence of new SARS-CoV-2 variants leads to the need for new treatment drugs development, which is a long process. However, with an increased patient mortality ratio, it is of utmost importance to repurpose existing drugs used to treat other viral diseases. 51 Failure to target the gene encoding the surface glycoprotein has been observed as SARS-CoV-2 variants are detected. 7 According to the latest guidelines of the Indian Council of Medical Research, approved test kits must employ multiplex RT-PCR assays, as the tests assessing only the surface glycoprotein may fail and produce false negative results. The viral entry into the host cell is facilitated by its successful binding to the angiotensin converting enzyme (ACE2) receptor. 52 Overexpression of ACE2 may lead to disease severity as observed in mice. 53 Lung damage can be reversed by blocking the renin-angiotensin pathway. 54 A recent study had shown that the surface glycoprotein of SARS-CoV-2 binds to ACE2 with a 10-to 20-fold higher affinity than other SARS-CoV surface glycoproteins, 55 which might be the reason for the high infectivity of SARS-CoV-2. Thus, viral entry into the host cell is a vital step which must be exploited for an efficient therapeutic development. There is a rapid ongoing search for therapeutic agents against SARS-CoV-2. Various computational studies have been conducted to discover potential drugs against SARS-CoV-2. 51,56,57 Recent studies have been based on the drugs targeting either surface glycoprotein or main protease of SARS-CoV-2. These approaches have led to the discovery of molecules with high binding affinities to these proteins. 58 The molecular docking analysis of surface glycoproteins with selected drugs for different VOCs along (Supplementary Figure 2, Extended data 44 ) revealed promising results for B.1.351 (Beta variant). Three drugs, namely baricitinib, lopinavir and remdesivir, showed maximum binding affinities against the Beta variant as compared to the WT and other VOCs. Other variants also expressed significant binding energies. As per a recent study, 10 the combination of baricitinib and remdesivir was more effective than Remdesivir alone and thus helped to lower the recovery time and accelerate the clinical status of patients suffering from COVID-19, especially those requiring high-inflow oxygen ventilation. Remdesivir efficiently inhibits the replication of SARS-CoV-2 by causing delayed chain termination when getting incorporated into the viral RNA. 59 However, it also showed considerable binding affinity when docked with the surface glycoprotein. 60 Here, our molecular docking study revealed its potential as an effective drug against SARS-CoV-2 VOCs. The high energy score resulting from these docking analyses suggests that these drugs may be recommended for administration to patients with B.1.351 infection. Molecular docking revealed that two RBD residues, namely Gln319 and Val503 facilitated a strong binding. Upon comparison with the WT, favipiravir showed a significant binding affinity with B. 1.1.7, B.1.617.2 and B.1.351. Favipiravir is a purine analog which inhibits the elongation phase of RNA synthesis. Favipiravir was proven to be effective in viral clearance and fast clinical improvement. 61 It has shown positive results in COVID-19 patients by improving patient's health. 62 In concordance with our findings, favipiravir was successfully docked with the surface glycoprotein of B.1.1.7 (Alpha variant). Dexamethasone has been widely used as a therapeutic intervention to treat COVID-19 patients. The docking score of dexamethasone with the surface glycoprotein of B.1.617.2 (Delta variant) was the highest (binding energy: -7.7 Kcal/mol) compared to other variants of concern. In contrast, the docking score for AY.1 (Delta plus variant) showed the lowest affinity with dexamethasone (binding energy: -4.4 Kcal/ mol). This observation shows that dexamethasone binding to the surface glycoprotein of the SARS-CoV-2 Delta variant (which has spread as one of the most dominant lineage worldwide) may represent an additional contribution to its efficacy in treating COVID-19. 63 Lopinavir is a drug approved by the FDA and serves as a protease inhibitor commonly used in the treatment of the Human Immunodeficiency Virus (HIV); it may be considered useful in the treatment against SARS-CoV-2 infection. 64 Our findings reveal that it may be a suitable choice for treatment as it shows significant binding with different VOCs, especially with B.1.351 (Beta variant; binding energy: -10.1 Kcal/mol). Remdesivir has been found to be a more potent drug than lopinavir, both in vitro and in MERS-CoV infected mice. 65 In concordance with our findings, remdesivir has shown more significant binding with B.1.351 than lopinavir as compared to WT. There is evidence of lopinavir being selective against other coronaviruses. 12 Despite the high binding energy with surface glycoproteins, our results encourage further in vitro and in vivo investigations. In comparison to the WT, the binding residues for different VOCs vary and some of them lie outside the receptor binding domain of the surface glycoprotein which does not have a direct role in drug affinity. However, they may impact the interaction of drug with surface glycoprotein through weak molecular interactions.

MD simulation and RMSD
The studied favipiravir-Alpha variant (B.1.1.7) ligand heavy atoms complex showed a RMSD of 1.028 Å. The Lennard-Jones potential and the binding potential of the complex was calculated to be -129.178 kJ/mol and -137.227 kJ/mol respectively. During the initial simulation run, up to 4 ns, the ligand-receptor showed a high RMSD value, which may be due to the stearic changes the protein underwent. The ligand managed to follow the change, keeping a stable number of H-bonds. The ligand initially stabilized at the binding pocket with some additional hydrophobic interactions. At t =14 ns, the pocket was obliterated, however, the ligand was kept near to the pocket until the pocket was opened. Subsequently, the ligand was restored to its original position ( Figure 2). MD simulation for favipiravir/Alpha variant was performed as it showed a significant binding affinity compared to other drugs under study. A similar procedure can be incorporated to conduct the simulation studies on other variants with drugs showing significant binding.

Conclusions
Drug repurposing may help to discover and identify the potential therapeutic effect of existing drugs against the genomic targets of SARS-CoV-2 virus. This study shows that the mutations (except Gln319 and Val503) outside the RBD of the surface glycoprotein of several VOCs do not largely affect the binding affinity of these drugs. No drastic structural change has been observed in variants irrespective of binding with the residues occurring outside the RBD of the surface glycoprotein. However, favipiravir showed the highest binding affinity against the Alpha variant, whereas dexamethasone showed approximately a 50% reduction in its binding affinity with the Delta plus variant when compared to the Delta variant, revealing that dexamethasone bound to surface glycoprotein of the Delta variant more strongly than the Delta plus variant. These residual fluctuations may play a role in antibody evasion and their molecular roles should be explored further. 66 However, the candidate drugs besides favipiravir and dexamethasone showed no significant alteration in the surface glycoprotein structure when compared to WT, implying that the current regimen of approved drugs can be continued in patients infected with these SARS-CoV-2 strains.
Further, molecular docking approaches offer great promise for predicting, shortlisting and quickly evaluating the anti-SARS-CoV-2 potential of candidate and existing drugs which can help timely effective interventions. The workflow depicting the study has been summarized in Figure 3.

Data availability
Underlying data Zenodo: Interaction of Surface Glycoprotein of SARS-CoV-2 variants of concern with Potential Drug Candidates:

Akinwunmi Oluwaseun Adeoye
Department of Biochemistry, Federal University Oye-Ekiti, Oye-Ekiti, Nigeria This study describes the interaction of surface glycoprotein of SARS-CoV-2 variants of concern with potential drug candidates: a molecular docking study. The authors concluded that the mutations which occurred outside the receptor-binding domains of the surface glycoprotein of several variants of concern do not affect the binding affinity of these drugs. The study is well-designed and well written. The materials and methods section is clear and detailed. Information on multiple sequence alignment and structural analysis, as well as molecular dynamics simulation, is sufficient. The results presented are good as well as the discussion section. It is an interesting study that I believe falls within the scope of this journal. The study contains stimulating information for the readers of F1000Research with an interest in repurposing existing antiviral drugs with potential activity against SARS-CoV-2 variants of concern. This paper is suitable for publication in this journal.
I have reviewed this article with keen interest. In general, the study design seems perfect. However, I observed in the structured abstract that the objective of the study was stated in the methods. From the study, the authors reported that the binding affinity of each drug for different variants was computed by assuming 100% binding affinity with SARS-CoV-2 wild type. Also, they considered the lowest binding energy and RMSD conformation as the most suitable docking pose. I believe that the drug molecule with the highest binding affinity to different variants of SARS-CoV-2 surface glycoproteins may not necessarily have the lowest binding energy and vice versa. The authors should have reported the binding affinities of the selected drugs with the SARS-CoV-2 surface glycoproteins separately in Kcal/mol from the docking study. Probably, the affinity of the drug molecules to the surface glycoprotein of SARS-CoV-2 variants of concern could be presented or represented in fold with respect to the SARS-CoV-2 wild type. The computation of the binding affinity of each drug concerning the assumption of 100% binding affinity with SARS-CoV-2 WT led to over 100% affinity of Baricitinib with beta variant, Favipiravir with alpha, beta, and delta variants as well as Remdesivir with the beta variant. What is the justification for these increases in affinity over 100% when compared with the reference SARS-CoV-2 wild type? Do any of the drugs alter the structures of SARS-CoV-2 surface glycoproteins?

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and is the work technically sound? Yes

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate?