In silico drug repurposing using molecular docking and dynamics to target the protein interaction between the SARS-CoV-2 S-glycoprotein and the ACE2 receptor

Introduction

The novel corona virus, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has had an unprecedented impact worldwide. It is the viral entity responsible for the Coronavirus disease 2019 (COVID-19) pandemic, which has reaped significant detrimental effects upon societies and economies worldwide. Up until present day, efforts remain ongoing towards effective infection management and control.¹

The SARS-CoV-2 virus is an enveloped single stranded RNA virus.² The entire viral genome has been sequenced and found to comprise of 29,881 bp encoding 9860 amino acids (GenBank no. MN908947).² The viral RNA codes for both structural and non-structural proteins, most notable of which are the structural S-glycoproteins or ‘spike proteins’- key to viral-host attachment and infection.² The S-glycoproteins are responsible for binding to the angiotensin converting enzyme 2 (ACE2) host cell receptor and initiating viral entry into the host cell.³ Both SARS-CoV and SARS-CoV-2 bind to the host ACE2 receptor to trigger viral entry and infection, however the binding affinity of the SARS-CoV-2 S-glycoprotein to ACE2 is shown to be over 20 times greater than that of SARS-CoV S-glycoprotein and ACE2 binding interaction.^2,4

The viral S-glycoprotein is a trimeric protein with a characteristic ‘stalk and halo’ like appearance.² Its peptide chain, composed of a total of 1237 amino acids, includes the S1 (aa residues 14–685) and S2 subunits (aa residues 686–1273), which have been characterized and found to be critical for viral attachment and membrane fusion.^5,6 In these regions, both the critical receptor binding domain (aa residues 319–541), as well as the fusion peptide domain (aa residues 788–806), were identified and characterized.^4,7

While vaccination efforts have greatly stemmed the spread of the infection, the main drawback of long-lasting vaccine efficacy appears to be the growing incidence of variants.⁸ Managing the disease with currently available and approved therapeutics, remains an essential approach for treatment. However, there appears to be is no consensus on effective management strategies yet.^9,10 Drug repurposing is an approach to identify novel indications and uses for approved medications. In the era of COVID-19 drug repurposing has been shown to be an essential path towards identifying potentially effective therapies for disease management.^11,12 Strategies towards this include in silico drug discovery tools and virtual screening, which have proven to be most valuable in identifying potential drugs that can be repurposed towards a new target.^11,13 Utilizing computational modelling platforms to visualize and identify ligand-target interactions is a powerful and efficient approach towards successful drug development.¹⁴

This study aims to apply a comprehensive virtual screening approach using multiple platforms to identify drugs and natural products that may potentially be repurposed to inhibit the binding interaction between the SARS-CoV-2 S-glycoprotein and the host ACE2 receptor. Furthermore, the interaction profile of all hit compounds is verified by molecular dynamic simulation studies to provide a detailed insight towards the stability and nature of the binding interaction at the viral-host interface.

Methods

Results

Docking studies

Preliminary docking studies to identify ligands that interact favorably with the receptor binding domain of the spike protein, involved screening a total of over 7000 approved medications and natural compounds (see Table 1). Initial docking studies were carried out on the S-RBD and ACE2 interface using Schrödinger Glide and Autodock Vina to identify hits that would prevent/disrupt the crucial protein-protein interaction. The results from both programs were analyzed and compared. The compounds that showed the most favorable binding profile and the best binding scores were shortlisted in Table 1.^29,30

From the initial hits 13 were identified that showed consistent good binding scores and binding poses. The compounds shortlisted bound to the protein-protein interface, and established interactions with key residues namely, His34, Asp38, and Lys353 from the ACE-2 binding region and Lys417, Gly496 and Tyr505 from the S-RBD. Hexoprenaline and tricocin, being highly flexible molecules and n-acetylglucosamine being a relatively small molecule were able to burrow deeply within the protein interface, while maintaining interactions with the above-mentioned key residues (Figure 1). The hit compounds all expressed high binding affinity scores using both platforms; ranging from -7.5 — -12.4 (Glide scoring), and -4.8 — -8.9 (Autodock vina scoring).

Figure 1. Surface representation of the 13 hits bound to the ACE-2 S-glycoprotein interface.

ACE-2 residues with carbon in green and S-glycoprotein residues with carbon in brown. The hit compounds are represented in yellow with hexoprenaline, n-acetyglucosamine, and tricocin highlighted in magenta.

Lymecycline, as an example, established key interactions via its side chain groups (Figure 2). The carboxylate formed a crossbridge at the interface; with strong charge assisted H bond between the Asp38 and Lys353 sidechains from ACE-2 and the Gly496 sidechain from the spike protein. The interaction was further strengthened by the ammonium ion of the terminal amino acid forming a charge assisted H bond with both the His34 backbone (ACE2) and Gly496 sidechain (S-RBD). A water mediated H bond for this ammonium ion with the Tyr453 sidechain of the S-RBD was also observed. Potential H bonds were also noted between the sidechain amide oxygen and Lys417, as well as the ring’s tertiary amine with Tyr505 from the spike protein. An increase of selectivity for lymecycline is expected owing to the presence of two opposite electrostatic interactions with both chains. The first a H bond between the sidechain amide and His34 sidechain (ACE-2) and the second is a H-π interaction between the central aliphatic ammonium and Tyr453 aromatic ring (S-RBD).

Figure 2. Binding interactions of lymecycline with the ACE-2 S-glycoprotein interface.

ACE-2 residues with carbon in green and S-glycoprotein residues with carbon in brown. Potential electrostatic interactions are represented as red dotted lines and distances are in Angstrom.

MD simulation results

The top hit molecules were selected based on a detailed visual analysis of the interaction profile of each compound with the target proteins. A hit was identified as any compound which had a favorable binding score and profile. Compounds were shortlisted based on reproducibility of the results using two different docking platforms (Glide and AutoDock Vina). 13 hit ligands were selected for molecular dynamic simulation studies (Table 2). Investigation of the optimal 2D and 3D docked positions reveled that each of the hit ligands form interactions with key residues of the binding pocket. Root mean square deviation (RMSD) plot analysis was used to measure the average displacement of atoms with respect to a reference frame. RMDS analysis revealed a stable binding profile for the top 7 ligands as shown in Figure 3 (RMSD fluctuation range of 2-4 Å). Of the 7 hit ligands, 3 exhibited a highly stable and robust binding profile within the interface binding pocket, these include: lymecycline, pentagalloylglucose and polydatin. The interaction profile for each individual ligand is depicted in Figures 4-6 (for MD results for the remaining 4 ligands see extended data). As the S-RBD-ACE2 interface structure was used for simulations, all interactions of the ligand with key binding residues of both chain A of the ACE2 binding domain and chain E of the S-RBD were considered as significant.

Table 2. The individual docking scores for the hit molecules using both AutoDock Vina and Schrödinger’s Glide at the interface of the S-glycoprotein and ACE2 receptor PDB 6M0J.

Ligand	Binding affinity (kcal/mol) Software: Schrödinger Glide	Binding affinity (kcal/mol) Software: AutoDock Vina
Colistin	-12.4	-4.8
Argipressin (vasopressin)	-9.7	-6.8
Lymecycline	-9.1	-8.3
Setmelanotide	-9.0	-6.1
Polydatin	-7.6	-8.5
Plazomicin	-10.6	-7.3
Hexoprenaline	-10.1	-6.5
Tricrocin	-9.1	-8.4
Ginsenosides	-11.3	-8.0
Ademetionine	-7.6	-7.4
N-Acetyl-D-Glucosamine	-7.5	-6.0
Pentagalloylglucose	-12.7	-8.9
Forsythiaside	-11.7	-8.7

Figure 3. Root mean square deviation (RMSD) graphs for the top hit compounds A: Lymecycline, B: Hexoprenaline, C: Pentagalloylglucose, D: Polydatin, E: Tricrocin, F: Setmelanotide, G: Forsythiaside.

The green graph shows fluctuations in the protein backbone from the initial reference point while the red shows the ligand fluctuations. The RMSD profile of the ligand is with respect to its initial fit to the protein binding pocket indicates that all ligands did not fluctuate beyond a 2-4 Å range.

Figure 4. Interaction diagram of Lymecycline with S-RBD-ACE2 interface binding pocket.

(A) Interaction of lymecycline with residues in each trajectory frame. The depth of color indicating the higher the interaction with contact residues; (B) The protein-ligand contacts showing the binding interactions fraction; (C) Lymecycline interactions with the protein residues during MD simulation. Interactions shown are occurring more than 30% during the simulation time. A: chain A of the ACE2 binding domain, E: chain E of the S-glycoprotein receptor binding domain.

Figure 5. Interaction diagram of Pentagalloylglucose with S-RBD-ACE2 interface binding pocket.

(A) Interaction of Pentagalloylglucose with residues in each trajectory frame. The depth of color indicating the higher the interaction with contact residues; (B) The protein-ligand contacts showing the bonding interactions fraction; (C) Pentagalloylglucose interaction with the protein residues during MD simulation. Interactions shown are occurring more than 30% during the simulation time.

Figure 6. Interaction diagram of polydatin with S-RBD-ACE2 interface binding pocket.

(A) Interaction of polydatin with residues in each trajectory frame. The depth of color indicating the higher the interaction with contact residues; (B) The protein-ligand contacts showing the bonding interactions fraction; (C) Polydatin interaction with the protein residues during MD simulation. Interactions shown are occurring more than 30% during the simulation time.

Lymecycline shows an exceptional binding profile with key binding residues of the binding pocket as shown in Figure 4. Key binding residues of the S-RBD Lys417 and Gly496 form a H bonding interaction with the lymecycline, furthermore it exhibited water bridge interactions with Tyr 505 and Ser494. While at the interface lymecycline was also shown to form interactions with key binding residues of the ACE2 binding domain including a strong H bond with His34 and mixed interactions (ionic and hydrogen bonds) with residues Asp38 and Lys353. All significant interactions are represented in Figures 4 A and B, which highlight residues with the strongest ligand interactions that are stable over the entire simulation period. Figure 4 C depicts all significant interactions displayed by the ligand and interacting residues occurring for over 30% of the simulation period. Of significance Lys417 and Gly496 are bound to lymecycline over 60% of the simulation period.

MD simulation results for pentagalloylglucose in Figure 5 show significant interactions with key residues of the binding pocket. Key binding residues of the S-RBD: Lys417 and Gly496, form H bonds, and Lys417 was also shown to form a water bridge interaction with the ligand. Additionally, a stable hydrophobic interaction with residue Tyr505 was observed. Strong interactions with key residues of the ACE2 binding domain were also observed in the simulation period, notably H bonds and water bridges with residues Asp38 and Arg393. All significant interactions are represented in Figures 5A‐C. Of significance Lys417 interacts with pentagalloylglucose approximately 85% of the simulation period.

Polydatin exhibits several stable interactions throughout the 100 nanosecond simulation period as depicted in Figure 6. Polydatin binds with key residues of the S-RBD Lys417 and Tyr505. A mixed binding profile is observed including H bond interactions, water bridges and for Lys417 hydrophobic interactions. Polydatin was also shown to form interactions with key residues of the ACE2 binding domain including a strong hydrogen bond with His37 and mixed H bond and water bridge interaction with residue Asp30. All significant interactions are represented in Figure 6.

Discussion

Over the last two decades, different docking tools and programs have been developed that use different algorithms in which the conformation of the ligand is extensively evaluated in a binding pocket until an energy minimum is reached. Most programs treat the ligand as a flexible component and the receptor as rigid while others treat both interacting components (ligand and receptor) as flexible.¹⁶ Such programs not only differ in the type of docking, but also in their ligand placement strategies.¹⁹ In this study, we selected Schrödinger’s Maestro and Autodock vina in order to assess the docking accuracy and mode of binding. Autodock vina was used to perform rigid docking while maestro was used for induced fit (flexible) docking. Autodock uses the genetic algorithm while maestro uses systematic search techniques for ligand placement.

This study was a comprehensive in silico investigation of a key target of the SARS-CoV-2 virus; the surface S-glycoprotein or spike protein. The interaction between the S-glycoprotein and target ACE2 host receptor is essential to the virulent nature of the virus. Agents that can disrupt this binding interaction may potentially inhibit viral entry into the host cell and infection. Structures selected for this study were all solved by x-ray crystallography or EM with a minimum accepted resolution of 2.45 Angstrom. For each of the target structures, the selected binding pockets from Schrödinger’s Sitemap were verified by ensuring the presence of essential binding residues identified in previous studies (Table 2). A detailed study of the viral host interface depicted in the 6M0J crystal structure revealed the key residues involved in strong charge assisted H bonds, ionic bonds, and strong H bonds, judged by the bond length, charge, and orientation (Table 3). Binding to these residues is necessary to disrupt the natural interaction between the two proteins in favor of the hit ligand.

In our study, a comprehensive screening library of over 7000 compounds was compiled, comprised of both FDA and worldwide approved drugs and nutraceuticals in addition to all natural products with established in vivo activity (Zinc)¹⁸ (Drugbank).¹⁷ In order to ensure robust and reproducible results virtual screening was preformed using two different platforms: Glide¹⁵ and Autodock Vina.¹⁶ The thorough protocol employed within this study ensures that the hit ligands identified had verified binding profiles. Further molecular dynamic simulation studies were vital towards verifying and visualizing the nature of the binding interaction of the hit ligands within the binding interface. Hits, that were selected for MD simulation studies, were determined to be the most likely to disrupt the strongest interactions between the viral and host proteins (Table 2).

MD simulation studies identified the following hit ligands: lymecycline, hexoprenaline, pentagalloylglucose, polydatin, tricrocin, setmelanotide and forsythiaside. All hits expressed a stable interaction profile as indicated by RMSD below 4 Å for both proteins and ligand position. It is important to highlight that all compounds were shortlisted from initial screening results based not only on their binding profile but also their suitability to be translated towards a clinical setting. Several drugs used in the management of COVID-19 have had detrimental effects owing to the adverse drug reaction profile of the employed therapeutic. The pathogenesis of COVID-19 culminates in several immune and cardiovascular manifestations ranging from hypercoagulability to kidney failure, as of such agents that are not appropriately selected may cause more harm than benefit within the overall scope of disease management.²⁰ Four hits were selected including lymecycline, pentagalloylglucose, polydatin, and hexoprenaline, which expressed very good RMSD profiles, indicating that the ligand remained bound in a stable manner within the binding pocket throughout the entirety of the simulation period. Other shortlisted ligands such as setmelanotide, and forsythiaside appear to have fluctuated with respect to their initial position relative to the protein backbone at some point within the simulation period although, the range of the fluctuation did not exceed beyond the range of 2-4 Å. The three final hits; lymecycline, pentagalloylglocose and polydatin, maintain strong bonds with residues within the binding pocket and in some cases vital binding residues throughout the entire simulation period.

Lymecycline, a broad-spectrum second-generation tetracycline antibiotic commonly used in the management of acne, gynecological and respiratory tract infections, was shown to exhibit stable binding with key binding residues of both the spike and ACE2 RBD. Lymecycline maintained its interaction with its side chain throughout the simulation period. The stable charge assisted H bond with Lys 353 from ACE-2 and Gly 496 from spike was conserved during the dynamic simulation which suggests that these interactions are more energetically favored over the initial H bond between the two mentioned residues. Other crucial interactions are shown with residues; Glu 37 and Asp 38 from ACE-2 and Tyr 505 from the S-RBD. These residues were bound via a H bond in the original protein-protein interaction. Lymecycline has been reported in previous studies to bind to additional viral targets in the SARS-CoV-2 virus, namely the main protease (M^pro).²¹ Like many members of its class it expresses anti-inflammatory properties,²² this makes it a particularly attractive subject that can be further investigated for repurposing studies in the treatment of SARS-CoV-2.

Pentagalloylglucose is a polyphenolic compound that has been shown in in vitro and in vivo studies to have anti-viral effects against the hepatitis C virus.²³ A number of recent in vitro studies have shown a dose dependent effect in inhibiting viral spike association with the host ACE2 receptor.²⁴ Our study reveals the detailed nature of the binding interaction of pentagalloylglucose at the S-RBD-ACE2 binding interface. Key residues Lys417, Gly496 and Tyr505 of the viral spike protein and residues Asp38 and Arg393 of the ACE2 binding domain, form a stable binding interaction with pentagalloylglucose. Interference with viral-host binding at these key residues is very likely to be vital for the inhibition of viral engagement with the host ACE2 receptor and subsequent cell entry.

Polydatin, another polyphenolic compound, is a glycoside precursor of resveratrol. In vitro studies have shown its potential to inhibit viral spike protein binding with the ACE2 receptor.²⁵ The findings of this study corroborate our own where similar binding scores were reported. Furthermore, our studies identify strong interactions of the key binding residues Lys417 of the S-RBD and Asp30 and Glu37 of the ACE2 binding region with polydatin. It appears that interfering with the binding of the viral Lys417 residue is a key inhibitory pattern detected in a number of docking profiles supported by in vitro findings.^24,25

Hexoprenaline, a β2 adrenoceptor agonist, is mentioned as a drug of interest although it was not identified from the top three- hits that showed the most stable interactions during MD simulations. It is for the first time reported to exhibit a favorable binding profile with a number of key binding residues of the S-glycoprotein RBD as well as the ACE2 binding domain. Hexoprenaline was shown to bind key residues and to interrupt key interactions. One of its side-chain nitrogen atoms formed two charge assisted H bonds with the two initially bound residues Glu35 from ACE-2 and Gln493 from S-RBD. The other chain nitrogen formed a charge assisted H bond with Glu37 (see extended data for MD results of remaining shortlisted compounds). The stability of these critical bonds during dynamics was slightly lower than those noted for lymecycline during the simulation period. Yet, considering its role as a bronchodilator and its pharmacological profile it may be considered a good candidate for drug repurposing in the management of COVID-19 infection where there is a high incidence of respiratory distress. However, caution must be taken, considering its potential nonselective activity on β1 receptors which may result in unwanted cardiovascular effects.²⁶

In this study lymecycline, pentagalloylglucose, and polydantin were identified as potential inhibitors of the S-RBD-ACE2 binding interface. Hexoprenaline could be also considered as a promising hit, due to its favorable docking and dynamic profile and taking into account its relevance and suitability for clinical testing. Of the nutraceuticals, forsythiaside also appears promising in its ability to potentially disrupt key binding interactions at the viral-host interface and is a prime candidate for further in vitro and in vivo studies. From a clinical perspective, both lymecycline and hexoprenaline may be considered as possible candidates for preliminarily clinical studies to assess their therapeutic potential in the management of COVID-19 infections.

Conclusions

Virtual screening is a highly attractive approach to identify potential compounds with therapeutic efficacy against target proteins, while efficient, it is not without its limitations. There have been a large number of docking studies published in the literature that have identified agents that may potentially be repurposed to inhibit SARS-CoV-2 targets. Our study is unique in its comprehensive approach to identify agents that can bind to the viral S-glycoprotein-ACE2 binding interface using multiple platforms. Molecular dynamic simulation studies were essential to identify a consistent binding pattern that appears to be common in the most effective agents that have the potential to inhibit the S-glycoprotein-ACE2 interaction. The agents identified in this study were additionally shortlisted for their suitability to be translated to a clinical COVID-19 setting by understanding their toxicity profile and identifying agents with verified anti-inflammatory and anti-viral capacity. The hit compounds identified are prime agents for further in vitro and clinical investigations to verify their efficacy in the potential treatment of COVID-19.

Author contributions

D.H. designed the research project. D.H., M.G., A.A., and A.H. conducted the computational studies. D.H. and M.G. analyzed the results. S.M., H.A., M.A., helped draft the manuscript. All authors were involved in writing, editing and revision of the manuscript.