3.1 Virtual screening and post-docking molecular dynamics simulations suggest putative SARS-CoV-2 RBD binders amongst FDA-approved drugs
Virtual screening (VS) of approved compounds represents the first step of drug repurposing strategies in the pursuit of candidates able to target SARS-CoV-2. Amongst the viral proteins that could be therapeutically exploted[24, 72], targeting the S glycoprotein could represent a promising way to treat and, most importantly, prevent COVID-19. Indeed, the ACE2:RBD binding interface (Figure 1), partially overlaps a pocket (Figure S2) that can be targeted with small molecules to disrupt the early stage of SARS-CoV-2 infection.
The use of MD structures, instead of cryo-EM or crystal structures, usually enhances the docking predictivity[47] thanks to the relaxation of the system from the solid phase to the fully hydrated one and the inclusion of full molecular flexibility. Extracting the RBD coordinates from the cryo-EM structure of the RBD-ACE2 complex (PDB 6M17) and processing the RBD through MD simulations allowed the RBD pocket to expand (Figure S1a,b) and to expose residues located within site 2 (Figure S1c,d). On this optimized structure, we docked 2421 approved drugs. The 23 top-ranked compounds from the VS (Table S1) were subjected to three MD replicas to validate the stability of the predicted complexes. Performing MD simulations from docking poses represents a valuable approach for evaluating the stability of docking-predicted complexes in a fully hydrated and flexible environment, therefore overcoming well-known limits of molecular docking[58].
Notably, none of these 23 drugs remained bound to RBD in all the three MD replicas (Table S1). The best-ranked compound according to Autodock, dihydroergotamine, completely dissociated during one replica, and moved away in the other two replicas, transiently interacting with other regions of the RBD (final RMSDs to the docked pose were 25.6 Å and 14.5 Å respectively). The second-ranked compound, nilotinib, displayed good stability with slight conformational rearrangement throughout the three MD simulations (final RMSDs to the docked pose 3.5 Å, 3.8 Å and 5.4 Å respectively). Cefsulodin and cromoglycate remained stably bound during two MD replicas out of three, while fluralaner, nafamostat, naringin, penfluridol, radotinib, and regorafenib were stable during one replica (Table S1).
The 13 simulations characterized by good stability were extended to 500 ns for further evaluating the stability of the complexes. Fluralaner, naringin, and regorafenib dissociated from the RBD (final RMSD to the docked pose 14.4 Å, 11.8 Å and 12.4 Å, respectively). Six compounds, on the other hand, remained bound: cefsulodin, cromoglycate, nafamostat, nilotinib, penfluridol, and radotinib (Table S1, Figure S3a). MMGBSA binding energies (Table S1, Figure 3Sb) indicated nilotinib and cefsulodin as the most stable ligands (-53.2 ± 4.1 kcal/mol and -41.3 ± 6.7 kcal/mol, respectively), while the other compounds displayed a similar interaction energy (e.g. values close to -30 kcal/mol). A comparison of the last frames from the six 500 ns MD simulations reveals similarities and differences in the conformations of the bound ligands and the RBD (Figure 2a). Direct interactions with cromoglycate and nafamostat, for example, folded RBD site 3 towards the putative binding pocket. As a general view, the six compounds occupied a pocket region partially overlapping RBD site 1 and site 2, and formed four main groups of interactions with the protein (Figure 2b). Within RBD site 1, the ligands formed hydrophobic p-p interactions with the side chain of Y505 and hydrogen bonds with the hydrophilic region in the proximity of backbone and the side chain of N501. Site 2, instead, stabilized the ligands through a mix of hydrophobic interactions deep in the pocket and hydrogen bonds with residues on the surface of the binding site (which is comprised of R403, K417, Y453, and G496).
3.2 Cefsulodin, cromoglycate, nafamostat, nilotinib, penfluridol and radotinib proposed binding modes
Cefsulodin (Figure 2c) is a third-generation b-lactam cephalosporin antibiotic with a spectrum of activity restricted to Pseudomonas aeruginosa and Staphylococcus aureus[48]. Interestingly, during the 500 ns MD simulation, cefsulodin rapidly dissociated from the starting docked pose before rebinding the RBD pocket after 100 ns with a different stable conformation (interaction energy -41.3 ± 6.7 kcal/mol, Table S1), Figure S3a, Video S1). In this binding mode, the drug formed hydrogen bonds with G502, G496, R403, and K417, and hydrophobic contacts with Y505, K417, I418, Y453, and Y495 side chains (Figure 2c, Video S1).
Interestingly, the other cephalosporins that were docked all obtained modest scores. The analog closer to cefsulodin in term of docking results was cefonicid (docking score = -7.0 kcal/mol, Figure S4). A direct comparison between cefsulodin, cefonicid and cefazolin (a further b-lactam cephalosporin, modest docking score of -6.2 kcal/mol, Figure S4) shows how different structures affected docking results: cefsulodin and cefonicid, indeed, bear a terminal phenyl ring which appeared to drive a common docking pose, while cefazolin’s tetrazole moiety produced a divergent flipped predicted conformation (Figure S4).
Cromoglycate (Figure 2d) is an anti-inflammatory drug indicated for the prevention of acute episodes in patients with mild to moderate asthma. During the 500 ns post-docking MD simulation cromoglycate experienced remarkable fluctuations within the binding site (Figure S3, Video S2) before forming interactions with residues that are part of RBD site 3 (i.e. G482, Figure 2d, Figure 2a, Video S2). The drug formed two persistent hydrogen bonds between the carboxylates and N501, Y453 side chains (Figure 2d). Hydrophobic interactions involved Y505, Y495, and Y453 (Figure 2d).
Nafamostat (Figure 2e) is a serine protease inhibitor currently in clinical trials for SARS-CoV2 infection in the light of its proposed ability to inhibit the transmembrane protease serine 2 (TMPRSS2)-dependent host cell entry[29, 30]. During our simulation (Video S3) nafamostat remained bound to the RBD pocket thanks to hydrogen bonds between the guanidinium group and N501, G502 and G496 (Figure 2e), as well as hydrophobic contacts with Y505 (MMGBSA binding energy = -29.8 ± 3.9). Transient interactions were formed between the molecule and the RBD site 3, in particular the Y489 side chain (Figure 2e, Video S3).
Nilotinib (Figure 2f) is a tyrosine kinases inhibitor with antineoplastic activity. In has been proposed to act in two different stages of the coronavirus infection cycle: in the early phases of infection, it has been proposed to inhibit the virion fusion with the endosome1, while at a later stage it has been proposed to inhibit the viral replication[59] via ABL-mediated cytoskeletal rearrangement[11, 20]. During the 500 ns MD simulation (Video S4), nilotinib, which remained stably bound to the RBD pocket (MMGBSA energy = -53.2 ± 4.1 kcal/mol), positioned the trifluoromethyl group deep in the binding pocket and formed hydrophobic interactions with I418, Y495, F497, Y505, and the alkyl chain of R403 (Figure 2f, Video S4). Hydrogen bonds were established between the drug and the N501 and R403 side chains (Figure 2f, Video S4).
Penfluridol (Figure 2g) is a long-acting antipsychotic agent[10]. During the 500 ns MD simulation (Video S5) the drug was anchored to the RBD through interactions between the piperidin-4-ol group and residues D406, Y453 (Figure 2g, Video S5), while the flexible bis(4-fluorophenyl)butyl and 4-chloro-3-(trifluoromethyl)phenyl moieties explored several conformations. Notably, penfluridol shares the phenylpiperidin-4-ol scaffold with the SARS-CoV inhibitor K22[42], which may indicate a common molecular mechanism.
Radotinib (Figure 2h) is a second-generation BCR-ABL tyrosine kinase inhibitor with therapeutic indication for patients resistant or intolerant to imatinib[4]. Despite the high structural similarity with nilotinib (from which it differs only by the pyridine terminal ring, Figure 2f,h) and the good interaction network within the RBD pocket (hydrogen bonds with R403 and Y453 side chains, hydrophobic interactions with Y505, I418, and Y495, Figure 2h), radotinib displayed lower stability (MMGBSA energy = -28.6 ± 3.8) and higher conformational fluctuation (Video S6) than nilotinib. The reason for this could reside in the different conformations predicted by docking (Figure S5). Indeed, the nilotinib top-ranked pose oriented the trifluoromethyl group toward the inner side of the pocket, while radotinib’s best pose was predicted to be almost planar on the pocket surface. This binding mode of radotinib did not evolve towards more stable states during the MD simulation (e.g. it did not reorient the trifluoromethyl group), resulting in sub-optimal interactions with the RBD. The different outcomes from almost identical ligands highlights the importance of the docking pose choice. The common practice to pick the top-ranked pose can lead to overlooking other suitable conformations that are penalized by lower docking scores computed without considering the explicit solvent contribution to the binding. However, the MD post-processing of docking poses can be time-consuming, and therefore applicable only to a small set of compounds[56, 58] on the basis of the ranking provided by docking scoring functions.
3.3 Simulating the RBD:ACE2 binding mechanism to retrieve insights into the effect of potential disruptors
During post-docking MD simulations, putative binders of the RBD interacted with protein residues that are responsible for direct contacts with ACE2 (Figure 1, Figure 2). Cefsulodin and nilotinib, for example, made stable contacts with residues located on both RBD site 1 (Y505, N501, G502, and F497) and site 2 (L455, Y495 and G496) that contributes to the stabilization of the complex with ACE2 (Table S2). This indicates the potential ability of these drugs to disrupt key interactions of the SARS-CoV-2 spike protein. To further test this hypothesis, we simulated the binding between the RBD and ACE2 in the absence of any small molecule or in presence of cefsulodin bound to RBD (conformation sampled after 500 ns of post-docking MD, Figure 2c). We focused on cefsulodin rather than nilotinib in the light of the spontaneous re-binding event observed during the simulation (Video S1, Figure S3), the numerous hydrogen bonds formed with the RBD (Figure 2c), and the low toxicity[71].
The first batch of SuMD simulations aimed to reproduce a putative binding path of the RBD:ACE2 hetero dimer (PDB 6M17). In one out of four replicas, the RBD was able to bind in close agreement with the experimental conformation reported in the PDB structure 6M17 (Video S7, Figure 3a,b). A further SuMD replica reproduced the complex with good agreement to the cryo-EM structure (Figure S6a), while other two simulations were not productive after 230 ns (Figure S6b,c). SuMD reproduced the stability and intermolecular contacts observed between the RBD and ACE2 during MD of the PDB structure 6M17 (Figure 3a,b), highlighting RBD site 1 residues T500, Y505, N501, and Q498 as the most involved in interactions alongside Q493, L455, K417 (site 2) and F486 (site 3)[5, 60, 63].
These data represented a robust reference for evaluating the potential effect of a ligand on the formation of the RBD:ACE complex. We then performed four SuMD replicas to study the binding between the RBD in complex with cefsulodin and ACE2 (Video S8, Video S9, Figure 3c,d, Figure S7). During two simulations, the RMSD of the RBD to the experimental conformation increased, indicating non-productivity of the simulated binding event due to stochastic reasons, rather than the presence of cefsulodin (Figure S7). The other two replicas, instead, were characterized by RMSD values decreasing to a plateau of about 10-15 Å (Figure 3c,d), consistent with a possible productive contact between the proteins. In one of the two productive SuMD simulations (Video S8, Figure 3c), despite the formation of a ternary complex involving the RBD, cefsulodin and ACE2, the interaction energy between the RBD and ACE2 oscillated around zero, indicating the formation of a sub-optimal binding interface between the two proteins (Figure 3c). Following a slight conformational rearrangement, cefsulodin remained bound to the interface between proteins for the whole simulation (GBSA interaction energy in the -20 to -40 kcal/mol range, Figure 3c, cf. GBSA interaction energy in the range of -32 - -50 kcal/mol in Figure S3). The RBD formed limited contacts with ACE2 (residues D38, Y41, and K353, Figure 3c) through Y449 and Q498 (site 1), and Q494 (site 2) side chains. A divergent scenario was sampled during the other productive SuMD replica (Video S9, Figure 3d). During this simulation, cefsulodin was displaced by ACE2 after 140 ns, allowing the formation of a metastable binary complex between RBD and ACE2 (RBD –ACE2 binding energy » -20 kcal/mol, against » -45 kcal/mol as computed for the experimental complex, Figure 3a,b); RBD residues F486, N487, S477 (site 3), Y489 Q493 (site 2), and T449 (site 1) were highly involved. Despite the displacement of cefsulodin, the RBD did not reach the experimental conformation over the simulated 230 ns (RMSD » 10-15 Å). However, over a longer timescale, it is plausible that the RMD eventually would rearrange to the fully bound conformation.
Taken together, these results suggest that a compound selected through VS and post-docking MD simulations may be less effective when the dynamic protein-protein formation is taken into account. SuMD (and other adaptive or enhanced MD sampling methods) may represent a further step in silico to the identification of protein-protein interaction (PPI) disruptors that may ultimately make the design process more effective.