Virtual screening and molecular dynamics simulations provide insight into repurposing drugs against SARS-CoV-2 variants Spike protein/ACE2 interface

After over two years of living with Covid-19 and hundreds of million cases worldwide there is still an unmet need to find proper treatments for the novel coronavirus, due also to the rapid mutation of its genome. In this context, a drug repositioning study has been performed, using in silico tools targeting Delta Spike protein/ACE2 interface. To this aim, it has been virtually screened a library composed by 4388 approved drugs through a deep learning-based QSAR model to identify protein–protein interactions modulators for molecular docking against Spike receptor binding domain (RBD). Binding energies of predicted complexes were calculated by Molecular Mechanics/Generalized Born Surface Area from docking and molecular dynamics simulations. Four out of the top twenty ranking compounds showed stable binding modes on Delta Spike RBD and were evaluated also for their effectiveness against Omicron. Among them an antihistaminic drug, fexofenadine, revealed very low binding energy, stable complex, and interesting interactions with Delta Spike RBD. Several antihistaminic drugs were found to exhibit direct antiviral activity against SARS-CoV-2 in vitro, and their mechanisms of action is still debated. This study not only highlights the potential of our computational methodology for a rapid screening of variant-specific drugs, but also represents a further tool for investigating properties and mechanisms of selected drugs.


Results
Exploring the druggable binding sites on Delta and Omicron Spike at ACE2-RBD interface. Our in silico druggability analysis using SiteMap identified in both Delta and Omicron S protein a pocket at the ACE2-RBD interface embedding residues important for the stabilization of the complex (hotspot). Table 1 reports the list of the hotspot residues which were identified through computational alanine scanning mutagenesis, i.e., substituting each residue at the interface with alanine and calculating the corresponding binding free energy change relative to the wild type (∆∆G bind ). Mutation of hotspot residues to alanine changes the binding free energy of at least 1.0 kcal/mol. In agreement with what has been previously observed for several other PPIs 21,22 , identified ACE2-RBD site had lower score than typical binding pockets (0.6 and 0.5 in Delta and Omicron, respectively). The site was lined by Arg403, Glu406, Tyr495, Gly496 (Ser in Omicron), Phe497, Asn501 (Tyr in Omicron) and Tyr505 (His in Omicron) ( Fig. 1) where Asn501 (Tyr in Omicron) and Tyr505 were main hotspots www.nature.com/scientificreports/ (Table 1). Edge of the larger Delta pocket also included Tyr453, Gln493 and Ser494. The predicted site overlapped with that of Spike RBD identified using FTMap and DeepSite and used as target for docking 23 . Analysis of the non-bonded interactions showed that in both Spike proteins residues at 496, 501 and 505 positions interact with Lys353 and Tyr41 of ACE2. In addition, Delta Gln493 and Omicron Tyr453 of Spike contact ACE2 Glu35 and His34, respectively. Given the challenges of identifying druggable sites on PPIs, pockets detected by SiteMap which rely on static protein conformations were assessed by molecular dynamics simulations (MD). The binding site volume was calculated along the MD trajectory and plotted over the simulation time using the script trajectory_binding_site_ volumes.py available in the Schrödinger suite of software. As evident in the Fig. 2, no significant differences existed in the binding-pocket volumes of the two simulated systems. The pockets were stable over the simulation time and an average binding-pocket volume of 53.0 Å 3 and 56.5 Å 3 was calculated for Delta and Omicron variants, respectively. A major fluctuation was observed in Omicron suggesting that its pocket structure may acquire a more expanded conformation as compared to Delta.    20 to screen the library of 4388 approved drugs, we were able to identify 505 potential PPIs modulators. Structure-based virtual screening, with different levels of increasing docking precision, was used to rank the 505 compounds based on their binding affinity. To obtain further accuracy for our protocol 24 , poses of Delta RBD-ligand complexes obtained in Glide XP docking were reranked with MM-GBSA which demonstrated to better distinguish between real and decoy poses of a ligand and assess the energetically preferred pose 25 . The ΔG bind values lower than − 47 kcal/mol were considered to retrieve the final set of molecules, leading to the identification of 20 top-ranked hits. Table 2 reports the calculated ΔG bind values of the selected compounds and the various energy components. Inspection of the free energy components in this Table revealed that for all the compounds the van der Waals and the lipophilic energies (ΔG vdW and ΔG Lipo ) contribute most to the ligands binding energy. As shown in Fig. 3, an interaction fingerprint (IFP) matrix was generated by Maestro to analyze the binding site residues occurring in any type of contact with the selected 20 compounds. Considering that a cutoff of 1 was used for the number of residue interactions, we found that these 20 ligands made interactions mostly with Arg403, Tyr453, Gly496 and Tyr 505, the last being the main hotspot residue (Table 1).   www.nature.com/scientificreports/ Crystal structure of Delta RBD-ACE2 complex was stable in the MD simulation. The stability of the RBD-ACE2 complex was analyzed through the root mean square deviation (RMSD) of the simulated protein structure's atomic positions from its native coordinates. The RMSD profile (showed in Fig. 4A) of the Cα atoms remained stable around average values of 2.4 ± 1.7 Å along the entire trajectory indicating that the complex was stable throughout the simulation. Root mean square fluctuation (RMSF) helped to understand the flexibility of each amino acid residue during the simulation time. As already described 26 , the binding interface of Delta RBD-ACE2 complex is made up of two patches of residues. Patch 1 consists of Ser19, Gln24, Phe28, Asp30, Lys31, His34, Leu79, Met82, Tyr83 on ACE2 and of Ala475, Asn487, Gln493, Tyr453, Phe486, Tyr489, Gln498 on RBD. Patch 2 comprises Glu37, Asp38, Tyr41, Gln42, Lys353 on ACE2 and Tyr449, Gly496, Thr500, Tyr505 on RBD.
The RMSF analysis showed that slight fluctuations occurred in patches 1 and 2 with RMSF values less than 1.5 Å and 1.3 Å for patch 1 and 2, respectively ( Fig. 4B, C). The RMSF of unbound RBD will provide a baseline for comparing the fluctuations with different ligand bound complexes.
Analysis of structural stability and binding free energy for Delta complexes. To better investigate the correct binding mode and estimate the stability of the predicted complexes, we performed unconstrained 100 ns MD simulations of the most promising 20 candidates within Delta RBD, followed by free energy of binding calculations using MM-GBSA. First, we monitored the conformational stability of the protein through MD simulations calculating Ca-RMSDs (  www.nature.com/scientificreports/ Compounds #4, #15 and #17 are very stable within the predicted binding site for the total simulation time. Compound #2 remained bound to the binding pocket, but while its hydroxy(diphenyl)methyl moiety ( Fig. 5) forms stable hydrophobic and hydrogen bonds interactions with Tyr505 and Gln493, respectively, all over the simulation time, the hydrogen bonding interactions of the methylpropionyl group of #2 with Arg452 or Ser494 were broken between 20 and 65 ns (giving rise to the increase of the RMSD lig ) and reformed between 65 and 90 ns (Fig. S1). Superimposition of RMSD and RMSF plots for the stable complexes is shown in Fig. 6. The Cα-RMSD of unbound RBD is also shown for comparison indicating that the protein's structure did not change during the simulation time (Fig. 6A). The variation of RMSD lig highlighted the considerable stability of the ligand position inside the binding pocket, and the major flexibility of fexofenadine (ZINC000003824921) (Fig. 6B).
In addition, the Cα-RMSF was also calculated to monitor structural fluctuations and get further insight about the binding mode analysis of selected ligands. It is important to point out that also the RMSF plots of the complexes were highly similar to that of RBD (Fig. 6C). This was especially remarkable for those residues critically involved in molecular recognition, where least fluctuations were observed (green bars in Fig. 6C). Analysis of trajectories indicated that the selected drugs remained close to their initial positions obtained from docking analysis, confirming the formation of stable complexes.
Simulation interactions diagrams during the entire MD run provided insights into the interaction pattern of the four ligands with Delta RBD. Figure 7 shows those interactions occurring during the MD simulation. Ligand #2, #4 and #15 exhibited hydrogen bonding and water bridges with Spike RBD through Gln493 and Ser494 and hydrogen bonding interactions with Tyr453. Hydrophobic interactions of #2, #4, #17 ligands with Tyr505 were also observed. Important hydrophobic interactions were observed for #17 with Phe456, Tyr489 and Tyr505. As it can be seen from Table 3 #17 exerted the most favored MM-GBSA energy (− 74.8 kcal/mol), with almost 18 kcal/mol, 24 kcal/mol and 27 kcal/mol difference comparing to #2 (− 56.4 ± 12.0), #4 (− 50.8 ± 16.7) and #15 (− 46.9 ± 4.6), respectively. The favored binding of hydroxy itraconazole can mainly be attributed to the lipophilic term of the binding energy (ΔG Lipo ) and to the non-polar (ΔG SolvSA ) contribution to the solvation free energy (Table 3).   Table 2) in comparison with the benchmark ligands suggest that our screening approach can achieve efficient inhibitor identification. Of note, more favorable binding energy was shown by fexofenadine when compared to the whole set of benchmark ligands. Molecular dynamics simulations of the best two benchmark ligands (Pixatimod and AB-00011778, Table 4) bound at the RBD binding site were also performed. Plots showing the stability of the ligands within the predicted site and comparison with the selected four drugs (fexofenadine, 3'-hydroxy repaglinide, RPR121056-d3, and hydroxy itraconazole) are presented in Supplementary Information (Figs. S2-S4). In all these systems the distance between the center of mass of the ligands and the center of mass of the Delta Spike RBD site remains constant over the simulation time ( Supplementary Information, Figs S2 and S4).   www.nature.com/scientificreports/ Å and 58.6 ± 9.7 Å, respectively). Compound #15 showed an average value of 9.1 ± 1.6 Å in the first 50 ns, 6.3 ± 0.7 Å between 50 and 80 ns and 10.8 ± 1.7 Å from 80 ns to the end of simulation (Fig. 8). Compound #15 was stable at the Omicron RBD-ACE2 interface, by means of hydrogen bonds with Arg403 and Gly502 and hydrophobic interactions with Tyr501 (Fig. 9). These data were corroborated by binding free energies values of the four selected compounds to the Omicron RBD calculated by the MM-GBSA method using the snaphots extracted from the MD trajectories. Table 5 shows that the first ranked compound was compound #15 and that hydrophobic interactions provide a favorable factor for association with Omicron RBD.

Discussion
The management of SARS-CoV-2 post-pandemic needs general health strategies also accounting for the possibility that this infection will become endemic 32,33 . In this context the vast research effort done so far, that must be celebrated for the great human achievements, should now give way to a new research approach, by means of fast pipelines for drug discovery and repositioning and social policies. The continuous variation of the Spike protein in fact causes concern of risk in terms of transmissibility and/or pathogenicity, also in population with high prevalence of vaccinated individuals. In this regards an unmet need is represented by the availability of fast  www.nature.com/scientificreports/ methodological approach able to promptly deepen the features of the putative novel variants and hopefully to predict the efficacy of molecules/drugs to be repurposed. The most successful drugs currently approved for clinical use, have been predicted and deepened in their mechanisms of action also through computational methodologies, such as in silico screenings 34,35 , immunoinformatic and structure-based drug design [36][37][38] . Several strategies have been approached to then identify and confirm computationally predicted molecules or to repurpose FDA-approved drugs to counteract the SARS-CoV-2 pandemic. Since the emergence of this infectious disease, a variety of in vitro methodology has been performed, mainly based on infected Vero 39 , HEK 37 , Calu-3, Huh7 cells 40 and hPSC lung organoids 41 , among the wide amount of other cell types. Very few in vitro tested molecules progressed to in vivo models and clinical studies and only a very small amount of them were successful at different phases of clinical trials.
Computational methods allowed not only to deepen the interacting properties of SARS-CoV-2 Spike with the ACE2, but also to evaluate a variety of natural and synthetic compounds as potential small-molecule specific inhibitors of Spike protein-ACE2 binding. Starting from these results and aware to transitioning to an era of endemicity with SARS-CoV-2 research, in this work we purpose a computational strategy for predicting if the therapeutic agents identified for past variants may continue to be useful against new ones. To the best of our knowledge few studies reported the use of molecular docking for repurposing already approved drugs as inhibitors of the interaction between Spike-RBD and the ACE2 receptor. For example, in 2020 a high-throughput virtual screening campaign of the University of Tennessee identified nitrofurantoin and isoniazid as potential SARS-CoV-2 treatment agents 44 . Targeting Spike-RBD and combining molecular docking with dynamics simulations few good anti-COVID-19 candidates (gonadorelin 17 , fondaparinux 17 , atorvastatin 17 and the antiviral  www.nature.com/scientificreports/ atovaquone 45 and praziquantel 45 ) were identified. In silico screening against ACE2 identified cefpiramide as a potential inhibitor of Spike protein-ACE2 binding 46 . The advantage of our approach is the use of a subset of the library of approved drugs which includes compounds with high probability to modulate the PPI interface and inhibit the infection in the host cells. From the site identification analysis, one pocket was identified at the Delta-and Omicron-ACE2 interface embedding hotspot residues. This pocket is shallow, and, from the MD simulations, we found that most predicted ligands were unable to form stable interactions with cavity-lining residues, this result confirming the importance of corroborating docking data by dynamic analysis. Here we used different approaches to select the compounds including rescoring, MM-GBSA, and a machine learning model to try to distinguish molecules able to interfere with the Spike RBD/ACE2 interaction. The computational pipeline gave us four top-hits, which exhibited good binding affinity towards Delta RBD. The compound with the most interesting features is fexofenadine, an antihistaminic drug already in use for the treatment of seasonal allergies (hay fever), skin itching and chronic idiopathic urticaria in adults and children 47 . Remarkably, antihistaminic drugs were found to exhibit direct antiviral activity against SARS-CoV-2 in vitro 48 . Fexofenadine, together with other histamine agonists, was also proposed as a potential COVID-19 drug by molecular docking studies against main protease in 2020 by Singh and co-workers 49 . Moreover, it has been demonstrated that patients hospitalized with COVID-19 have had benefits after the administration of a single or combined antihistaminic drugs 50,51 , correlating with increased survival and halted COVID-19 symptom progression, especially in the ones with critical conditions 52 . It has been suggested that histamine agonists could be able to down-regulate the excessive cytokine release, although the potential mechanism(s) of action about their impact on COVID-19 symptoms is(are) still unclear. Similarly, it is debated their usage in persons experiencing long-Covid symptoms (PASC, Post-Acute sequelae of SARS-CoV-2), where anecdotal descriptions have been reported 53 . Our data on the one hand corroborate these clinical reports, suggesting and giving new insights into the putative mechanism(s) of action of antihistaminic drugs on SARS-CoV-2, and, on the other hand, highlight the potential of this methodology for a rapid screening of variant-specific drugs.
Other compounds showed a good inhibitory ability against Spike-ACE2 complex both for Delta and Omicron variants: repaglinide, a hypoglycemic drug, approved for the treatment (in combination with rosiglitazone or pioglitazone) of type 2 diabetes 54 ; RPR121056-d3, a metabolite of Irinotecan (CPT-11), already approved as antineoplastic agent 55 ; itraconazole, already approved for the treatment of certain fungal infections (candidiasis and histoplasmosis) 56 .
Fexofenadine has been purposed as a potential target drug for COVID-19 treatment 49,57 . This substance has antihistamine properties and is indicated for the symptoms of allergies, through blocking of H1 receptors. In SARS-CoV-2 we agree with previous works that hypothesized that fexofenadine could impact indirectly on the cytokine storm by inhibiting histamine and, consequently, interleukin-6 (IL-6) production, directly on the viral replication, interacting on the protease enzyme MPro 49,58 . A similar direct mechanism of action has been proposed for repaglinide, an antidiabetic drug, that in analogy to other antiviral drugs (e.g., nelfinavir) has been candidate as a potential inhibitor of the MPro. Further experiments would be needed to determine the mechanism(s) of action of these two drugs. While RPR121056-d3 has not yet been reported as a potential treatment of COVID-19, itraconazole showed an in vitro activity against SARS-CoV-2, possibly by inhibiting oxysterol-binding protein that interferes with intracellular lipid transfer 59 and by favoring the cholesterol accumulation in the endosomal membrane preventing the viral transfer into the target cells 60 .
Interestingly all these drugs, although with different mechanism of actions and with different composition, are already used to modulate immunopathologic processes, including fexofenadine, which is recommended for immune-mediated respiratory diseases, more specifically in upper airway diseases and allergic rhinitis 61 . Our results from prediction and in silico studies have thus succeeded in selecting drugs that, with distinct suggested or proven mechanisms of action, appear to be able to interfere with the virus' ability to infect and damage host cells. To date one of the unmet needs in counteracting the spread of SARS-CoV-2 and the generation of new variants is represented by the control of novel or recurrent infections, both for adults and children and for fragile and elderly patients [62][63][64][65] . In this regard, among the selected drugs, fexofenadine in a ready-to-use spray formulation, could be a valid option, as suggested also for other drugs 66 . These drugs, identified via computational approaches can be assessed experimentally and provide a basis to develop novel treatment options against COVID-19 variants.

Materials and methods
Proteins preparation. Crystallographic 26 and cryo-EM 26 structures of SARS-CoV-2 Delta and Omicron RBD, bound to ACE2, were downloaded from the protein data bank (PDB) (PDB codes: 7WBQ and 7WBL, respectively) and used as targets for virtual screening. RBDs were extracted from the PDB structures of ACE2-RBD complexes and prepared using the "Protein Preparation Wizard" tool of the Schrödinger suite (Schrödinger Release 2021-4: Protein Preparation Wizard; Epik, Schrödinger, LLC, New York, NY, 2021). The protocol included, following water molecules and cofactors removal, correcting mislabeled elements, adding hydrogen atoms, assigning bond orders, hydrogen bond optimization, and restrained energy minimization using OPLS4 force field 67 . The prepared proteins were then considered for grid generation using the "Receptor Grid Generation" panel of the Glide module of the Schrödinger suite (Schrödinger Release 2021-4: Glide, Schrödinger, LLC, New York, NY, 2021). The center of each grid (size 15 × 15 × 10 Å) was arranged at the centroid of the predicted druggable site using SiteMap 20,68 . Identification of druggable pockets. SiteMap was used to explore possible druggable pockets on Delta and Omicron RBD surfaces 69,70 . The program uses OPLS4 force field to estimate the interaction energies of probes placed at all points along a three-dimensional grid that encompasses the entire RBD with a minimum www.nature.com/scientificreports/ of 15 site points per site. SiteMap ranks identified sites using two druggability assessment scores: SiteScore and Dscore, which characterize the binding site in terms of size, exposure to solvent, hydrophobicity and hydrophilicity, degree of hydrogen bond donation and acceptance. The highest ranking regions were further selected for embedding hotspot residues at RBD-ACE2 surface 20  were placed in a cubic water box at a buffer distance of 10 Å and solvated with SPC water models. A 0.15 M NaCl salt concentration was added and additional Na + /Cl − ions were added to neutralize the systems. The particlemesh Ewald method was used to calculate the long-range electrostatic interactions. A cut-off radius of 9.0 Å was applied for short-range van der Waals and Coulomb interactions. Each solvated system was minimized and equilibrated using the default protocol of Desmond in Maestro which includes 2 NVT and 2 NPT restrained short simulations. All equilibrated systems were then subjected to a MD run with periodic boundary conditions in the NPT ensemble using OPLS4 force field 78 for 100 ns. The temperature of 300 K and the pressure of 1 atm of the systems were maintained by the Nosè-Hoover chain thermostat and Martyna-Tobiase-Klein barostat methods, respectively. The binding energy between the Spike RBDs and the docked ligands was calculated over the 100 ns period with thermal_mmgbsa.py python script provided by Schrödinger which takes a Desmond trajectory file, splits it into individual snapshots, runs the Prime-MMGBSA calculations on each frame, and yields the average calculated binding energy. Best compounds for Delta RBD were also analyzed by molecular dynamics simulations on Omicron RBD. 100 ns MD simulations of the best benchmark ligands bound at the RBD site and 200 ns MD simulations of Delta RBD/ACE2 complex were also performed using the same conditions. Details of the starting systems for MD simulations are reported in Supplementary Materials (Table TS1).

Conclusions
The possibility to identify and reposition drugs via computational analyses can address the current challenges related to the succession of waves of the infection of SARS-CoV-2 variants worldwide. This strategy allowed the prediction of selective compounds consistent with already suggested and prescribed drugs and is confirmed as a useful tool in support to clinical therapeutic approaches.

Data availability
Docking results and simulation trajectories datasets are freely accessible at zenodo.org as https:// doi. org/ 10. 5281/ zenodo. 68685 48. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.