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

Materials and Software

The molecular modelling software Maestro by Schrödinger¹⁵ and AutoDock Vina,¹⁶ were both used for virtual screening and molecular dynamics simulation studies. The desktop workstation was equipped with Intel® Core™ i7-10700F Processor, Linux Ubuntu 22.10 operating system and a RTX 5000 graphics card.

Crystal structures

The protein crystal structures were retrieved from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) (RRID:SCR_012820). For virtual screening and molecular dynamic simulations, the structure of the SARS-CoV-2 S-glycoprotein receptor binding domain (S-RBD) bound to the ACE2 receptor interface was used (PDB code: 6M0J).

Virtual databases

A ligand library was compiled comprised of Food and Drug Administration (FDA) approved and worldwide approved drugs as well as approved nutraceuticals and natural compounds with verified in vivo efficacy. A total of 1100 FDA approved drugs were retrieved in SDF format from DrugBank database (RRID:SCR_002700),¹⁷ along with a library of 3440 worldwide approved medications and 74 approved nutraceuticals. Furthermore, 1537 natural compounds with verified physiological effects in vivo, were retrieved from Zinc database (RRID:SCR_006082).¹⁸ Table 1 summarizes the ligand categories of the compiled virtual library.

Maestro

Protein preparation

All crystal structures were prepared using the Schrödinger Maestro’s (RRID:SCR_016748) protein preparation wizard tool. Structure preparation and minimization was done at a pH of 7.4 with corrected ionization states, polar hydrogens were added, and non-essential water molecules were removed. The entire structure was minimized and optimized with the OPLS3 force field and the default value for the RMSD of 0.30 Å was used for non-hydrogen atoms.

Ligand library preparation

Downloaded SDF structures were prepared for docking studies using the maestro Ligprep tool (RRID:SCR_016748). Structures were converted into 3D maestro format, and ionization states and chirality were optimized at physiological pH (7.4) using OPLS3 force field. The final 3D conformations were utilized for virtual screening.

Binding pocket determination and docking studies

In Schrödinger’s Maestro (RRID:SCR_016748) the binding pocket was identified using Schrödinger’s Sitemap, a single binding pocket between the interface of the RBD and ACE2 binding region (PBD 6M0J) was used. The selected binding pockets were verified by ensuring the presence of essential binding residues identified in previous studies.⁶ The site score value was within the range of 1-1.1 for each of the binding pockets indicating high accessibility and druggability of the selected binding pocket. The selected binding pockets were used to generate a docking grid using Maestro’s Glide module for docking studies. The receptor grids were generated using the prepared proteins, with the docking grids centered on the identified receptor binding pocket for each protein. A receptor grid was generated using a 1.00 van der Waals (vdw) radius scaling factor and 0.25 partial charge cut-off. The binding sites were enclosed in a grid box of 20 ų without constraints using default parameters. Docking was repeated and verified using three screening settings. All compounds were screened under a high-throughput docking setting and the top 200 compounds with the highest binding scores where then selected for standard precision docking, of these verified hits the top 80 compounds where further verified using extra precision (XP) docking settings. The ligands were docked using the extra precision mode (XP) without using any constraints and a 0.80 van der Waals (vdw) radius scaling factor and 0.15 partial charge cut-off. Induced fit docking was carried out with flexibility of the residues of the pocket in close proximity to the ligand.

GlideScore implemented in Glide (RRID:SCR_016748), was used to estimate binding affinity and rank ligands. The XP Pose Rank was used to select the best-docked pose for each ligand. The final list of thrice verified compounds was then analysed in detail based on binding scores and a detailed study of all binding interactions.

Molecular dynamics simulation studies

Molecular dynamic (MD) simulation studies were carried out using the Desmond Module on Schrödinger’s Maestro platform (RRID:SCR_016748). The protein preparation wizard was used to minimize the hit protein-ligand complex and the simulation environment was built using the system builder application of Desmond. A water based solvent system: TIP3P was employed to generate the simulation environment contained within an orthorhombic simulation box with 10 Å buffer parameter from the protein surface. The system was neutralized and isotonic conditions attained via the addition of counter ions and 0.15 M NaCl. All MD simulations were conducted at a temperature of 300 K and a pressure of 1.013 bar. A simulation period of 100 nanoseconds was run for each of the hit ligand-protein complexes. Analysis calculations were subsequently run and results presented using the simulation interaction diagram tool of Desmond.

Binding free energy calculations MM/GBSA

Binding free energy; molecular mechanics generalized Born surface area (MM/GBSA) was calculated using the PRIME module in Maestro (RRID:SCR_016748). MM/GBSA quantifies the difference in binding free energy calculated between the stabilized complex and individual ligand and receptor conformation.¹⁹ The final complex configuration from compiled trajectories of a 100 ns MD simulation was used to calculate MMGBSA free energy, calculations were preformed using the OPLS_2005 force field and VSGB 2.0 solvation model.²⁰

Autodock vina

Docking calculations were carried out using the AutoDock vina software version 1.1.2 (RRID:SCR_011958).¹⁶ All hydrogens were added to the ligand PDB file and Gasteiger charges were computed and all the torsion angles of the ligand were defined using the autodock-tools program. A grid box generated with the following dimensions: 36×24×-4 Å, with a grid spacing of 1 Å was used. The Lamarckian genetic algorithm was used as a search method with a total of 30 runs (maximum of 20 000 000 energy evaluations; 27 000 generations; initial populations of 150 conformers). The binding affinity calculation in AutoDock vina together with analysis of binding interactions were used to select hits for molecular dynamic simulation studies.

In-silico pharmacokinetic calculations

The Absorption, Distribution, Metabolism, and Excretion–Toxicity (ADMET) properties of the top seven hit drugs and nutraceuticals were derived and calculated from the PubChem database²¹ and using SwissADME²² and ADMETlab2.0²³ respectively.

Data availability

The crystal structures analysed during the current study are available in the PDB database, [PDB ID: 6M0J]. The molecular structures of the compound’s datasets used were obtained from DrugBank database (https://go.drugbank.com/drugs) and Zinc database (https://zinc.docking.org/substances/subsets/fda/).

Docking studies

Upon designing the study and to this date there are no known solved crystal structures for the SARS-Cov-2 RBD and ACE2 complex bound to a ligand within the interacting binding pocket, as of such for the purpose of validation two different docking platforms were employed in the initial docking investigations to screen and shortlist suitable candidates targeting our region of interest in the viral-host interface. The Site Score tool on Maestro was used to identify the region with the highest druggability. A score of 1.1 for the region of interest identified revealed high druggability and accessibility of this binding pocket and most importantly was selected to encompasses the viral-host interactions of interest. The key interactions between the viral host and receptor have been determine and verified in previous studies (Table 5). Validation drawn from prior experimental evidence revealed a consistent pattern associating good docking scores with experimental compounds that have been shown to interrupt the SARS-Cov-2 RBD and ACE2 interaction, complementing the findings of our study and supporting this approach of identifying agents capable of disrupting key viral-host interactions (see extended data).

Preliminary docking studies to identify ligands that interact favourably 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 disrupt crucial protein-protein interactions. The results from both programs were analysed and compared. The compounds that showed the most favourable binding profiles and the best binding scores were shortlisted in Table 1.^34,35

A total of 13 were identified that consistently showed 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 to -12.4 (Glide scoring), and -4.8 to -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 group formed a crossbridge at the interface; with strong charge assisted H bonds shown 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 being a H bond between the sidechain amide and His34 sidechain (ACE-2), and the second 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 favourable binding score and profile. Compounds were shortlisted based on reproducibility of the results using two different docking platforms (Glide and AutoDock Vina). A total of 13 hit ligands were selected for molecular dynamic simulation studies (Table 2). Investigation of the optimal 2D and 3D docked positions revealed 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.

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.

Binding free energy calculations MM/GBSA

Results for MM/GBSA binging free energy calculations are reported as a change in thermodynamic binding free energy (ΔGBind) of the ligand-protein complex in its most stable complexed configuration relative to the free energy of individual system binding partners and parameters in noncomplexed forms. The change in free energy was reported for the following compounds: lymecyclin, polydatin, setmelanotide, hexoprenaline, tricrocin, forsythiaside and pentagalloylglucose in Table 3. Results of MM/GBSA calculations complement the findings of the docking experiments and MD simulation. The ΔGBind is a composite of multiple components including the change of energy estimated from the formation of covalent bonds (ΔGCovalent), hydrogen bonds (ΔGHbond), lipophillic interactions (ΔGLipo), solvation energy (ΔGSolvation), and Van del Wal forces (ΔGvdW), each shedding light on distinct interaction dynamics between the ligands and the S-RBD-ACE2 binding pocket.

In-silico pharmacokinetic calculations

The Absorption, Distribution, Metabolism, and Excretion–Toxicity (ADMET) profiles for the top seven hit drugs and nutraceuticals are listed in Table 4. The majority of the ligands are FDA approved drugs or nutraceuticals with known and reported physiochemical properties which were derived from the PubChem database, for unknown properties SwissADME and ADMETlab 2.0 were used to predict the parameters of interest. All ligands were found to have favorable ADMET profiles with no significant pharmacokinetic or toxicity parameters of concern.