Supervised molecular dynamics for exploring the druggability of the SARS-CoV-2 spike protein.

The recent outbreak of the respiratory syndrome-related coronavirus (SARS-CoV-2) is stimulating an unprecedented scientific campaign to alleviate the burden of the coronavirus disease (COVID-19). One line of research has focused on targeting SARS-CoV-2 proteins fundamental for its replication by repurposing drugs approved for other diseases. The first interaction between the virus and the host cell is mediated by the spike protein on the virus surface and the human angiotensin-converting enzyme (ACE2). Small molecules able to bind the receptor-binding domain (RBD) of the spike protein and disrupt the binding to ACE2 would offer an important tool for slowing, or even preventing, the infection. Here, we screened 2421 approved small molecules in silico and validated the docking outcomes through extensive molecular dynamics simulations. Out of six drugs characterized as putative RBD binders, the cephalosporin antibiotic cefsulodin was further assessed for its effect on the binding between the RBD and ACE2, suggesting that it is important to consider the dynamic formation of the heterodimer between RBD and ACE2 when judging any potential candidate.

Introduction

In early 2020, the severe acute respiratory syndrome-related coronavirus (SARS-CoV-2) caused the spread of the coronavirus disease (COVID-19) on an unprecedented global scale. Symptoms comprising smell and taste dysfunction [1], fever, cough, fatigue, headache and dyspnea [2] usually appear after an average incubation period of just over 5 days [3], lasting for 6 to 41 days [4]. The clinical picture is worsened by severe pneumonia with high levels of pro-inflammatory cytokines [2] (a condition also known as “cytokine storm”), which triggers extensive organ damage. SARS-CoV-2 is a positive-stranded RNA Betacoronavirus [5] in which about one-third of the genome codes for structural proteins such as envelope proteins, nucleocapsid proteins, membrane proteins and the spike (S) glycoprotein [6].

The S glycoprotein gives the characteristic “crown” shape to the coronaviruses phenotype, by forming prominent structures on the virion envelope, and drives virulence by mediating the first interactions with the host cell [7] and inducing immune responses [8]. The trimeric S glycoprotein is composed of two subunits, S1 and S2, respectively mediating the binding to the host cell and membrane fusion [9, 10]. Specifically, S1 comprises the receptor-binding domain (RBD, Fig. 1a) for the angiotensin-converting enzyme 2 (ACE2), and stabilizes the prefusion state of the S2 subunit, which in turn drives the fusion of the viral particle into the host cell membrane [11]. ACE2 is a type 1 membrane protein with an extracellular peptidase domain (PD) responsible for the maturation of the vasoactive hormone angiotensin [12]. Organs with high expression level of ACE2 are lungs, testis, heart, kidneys, and intestine [13, 14].

Fig. 1

Supramolecular organization of the spike protein:ACE2 complex and the RBD:ACE2 binding surface. a Schematic representation showing the dimeric ACE2 protruding from the cytoplasmic membrane and acting as receptor for the trimeric SARS-CoV2 spike protein in prefusion open conformation, which binds through one of the three available RBD in up position. b two side-views of the contact interface between the RBD (orange) and ACE2 (purple). The three RBD sites (S1, S2, and S3) primarily responsible for binding ACE2 α1 helix are indicated within parentheses

An intriguing way to tackle the SARS-CoV-2 infection is to disrupt the entrance of the virus into host cells by blocking the molecular machinery implied in this fundamental step [15-18]. The structural basis of the interaction between ACE2 and the S glycoprotein RBD (Fig. 1) could drive the discovery of small molecules and monoclonal antibodies [19-21] able to slow down, or even prevent, the development of COVID-19 [22-25]. An ideal drug candidate should selectively target the RBD without interacting with ACE2, to avoid possible side effects linked to angiotensin physiology [26, 27]. The S glycoprotein RBD binds to the α1 helix of the ACE2 PD mainly through polar interactions (Fig. 1b). Three different sites on RBD can be distinguished (Fig. 1b, Figure S3): site 1 (residues Q498, T500, N501, and Y505) and site 3 (N487 and F486) interacts with ACE2 α1 helix C (residues Q24 and T27) and N (Y41, Q42, K353, and R357) termini, respectively. The RBD site 2 (residues R403, L455, F456, Y453, and Q493) interacts with the central segment of the α1 helix (D30, K31, H34, and D38).

Given the fast rate of diffusion of the pandemic, an outstanding number of drug repurposing campaigns have started [28]. The advantages of repurposing an “old” drug to treat new diseases reside in the shorter development timelines and costs, due to the low-risk safety profile of already approved drugs [29]. In silico repurposing approaches usually rely on the virtual screening (VS) of small drug databases, followed by a computational validation through molecular dynamics (MD) simulations to evaluate the stability of the predicted complexes in a flexible and fully hydrated simulated environment. Some examples of potential disruptors of the S glycoprotein interaction with ACE2, identified by simple molecular docking of approved drugs include hesperidin [30], paritaprevir [31], cladribine, clofarabine, and fludarabine [32]. Differently, docking combined with MD simulations highlighted denopamine, bometolol, naminterol, rotigaptide, benzquercin [33], simeprevir, lumacaftor [34], tegobuvir, bromocriptine, baicalin [35], KT185, KT203, GSK1838705A, and BMS195614 [36] as potential binders of the RBD.

Assuming some redundancy between the libraries of compounds considered, the lack of consistency from different groups could be considered surprising. However, causes for this discrepancy reside in the RBD coordinates used for VS, the different protocols for ligand preparation and docking, the number of MD replicas and the total amount of post-docking MD sampling, the different force fields employed and the free-energy method used for rescoring binding modes. Here, we performed the VS of 2421 FDA-approved drugs against MD-derived coordinates of RBD. Post-docking MD simulations (in triplicate and up to 500 ns) of the top-ranked 23 complexes highlighted six promising compounds, which were further evaluated according to the molecular mechanics energies combined with the generalized Born and surface area continuum solvation (MMGBSA) energy. Cefsulodin, which stood out alongside Nilotinib as a potential disruptor of the RBD–ACE2 interface, was then evaluated through supervised molecular dynamics simulations [37-39] (SuMD) aimed to reconstruct the binding event between the RBD and ACE2. Overall, cefsulodin hindered the formation of the heterodimer. However, cefsulodin’s affinity for the RBD was affected by metastable interactions with ACE2, indicating the importance of including the heterodimer during the screening of potential candidates.

Methods

The sequential steps of the computational protocol followed in the present work (Figure S1) are here described.

Classic molecular dynamics simulation of the spike protein RBD

The coordinates of a monomer of the spike protein (S protein) RBD were retrieved from PDB entry 6M17 [25], which reports the pre-fusion open conformation of the trimeric S protein (one RBD in up position) in complex with one monomer of ACE2. Only residues C336-F515 from one monomer of the S protein were considered. After removing the N-acetylglucosamine residues (which are not located in the RBD pocket [22]), hydrogen atoms were added using pdb2pqr [40] and the propka [41] software. The resulting system was prepared for simulations with the Amber14SB [42] force field by creating a 90 Å × 92 Å × 73 Å TIP3P [43] solvated box. Overall charge neutrality was maintained by adding two Cl− ions.

ACEMD [44] was used for both equilibration and 200 ns of MD productive simulation. Restraints were applied to protein alpha carbon atoms and gradually released throughout 2 ns of equilibration before starting the MD production. Productive trajectories (Table S4) were computed with an integration time step of 4 fs in the canonical ensemble (NVT). The target temperature was set at 300 K, using a thermostat damping of 0.1 ps−1; the M-SHAKE algorithm [45, 46] was employed to constrain the bond lengths involving hydrogen atoms. The cut-off distance for electrostatic interactions was set at 9 Å, with a switching function applied beyond 7.5 Å. Long range Coulomb interactions were handled using the particle mesh Ewald summation method (PME) [47] by setting the mesh spacing to 1.0 Å.

The 200 ns trajectory was clustered into 10 groups of frames using the VMD plugin Clustering (at < https://github.com/luisico/clustering), according to the RMSD of the residues delimiting the RBD pocket (residues Y453, R403, and K417) with respect to the starting conformation. A visual inspection of the trajectory revealed that the breathing of the RBD pocket was primarily due to side chain movements. R403 and K417, in particular, heavily contributed to shaping the pocket because the side chains are longer than those of other neighboring residues. Moreover, position 417 is occupied by valine in the SARS-CoV RBD, indicating a potential element of resistance or selectivity for drug candidates. The frame from the most populated cluster with the highest RMSD value (representative for a “relaxed” and more open conformation of the pocket) was extracted and used for the subsequent steps. This double criterion should allow selection of a structure quite close to the RBD average conformation, taking into account also the possible induced-fit produced by ligand binding, which might favor a more open pocket state.

The RBD from both PDB 6M17 and the MD frame were used as input for the FTMap [48, 49] and DeepSite [50] webservers to evaluate the RBD binding site before and after the MD simulation of RBD (Figure S2). As shown in Figure S2, the binding pocket is larger in the MD structure than in the RBD structure from PDB 6M17.

Virtual Screening of the approved compounds library

The RBD from the MD frame was used to dock an SDF library of FDA-approved drugs [51] (available at https://www.selleckchem.com/screening/fda-approved-drug-library.html).

RDKit (https://www.rdkit.org/) was used to remove organometallic compounds and entries with more than 12 rotatable bonds, strip counterions from salt molecules, and generate initial 3D conformations. The resulting database was further processed through Chimera [52] and Open Babel [53] to check the protonation state of the compounds, assign MMFF94 [54] partial charges and generate final pdbqt files for the 2421 structures.

Docking simulations were performed employing AutoDock Vina [55] within a cubic 25 Å × 25 Å × 25 Å volume centered at the oxygen atom of the Y495 side chain, therefore focusing on the main RBD pocket shown in Figure S2b. Ten conformations for each run were generated.

Post docking molecular dynamics simulations

The 23 compounds (Table S1) with the best docking scores were prepared for classic MD with the Amber14SB [42]/GAFF [56, 57] force field by creating a 77 Å × 95 Å × 71 Å TIP3P [43] solvated box (a padding of 15 Å along the x, y, and z dimensions was considered). Overall charge neutrality was maintained by adding Cl− or Na+ counterions. ACEMD [44] was used for both equilibration and three replicas of duration 100 ns each. Restraints were applied to protein alpha carbon atoms and gradually released throughout 2 ns of equilibration before starting the MD production replicas.

After 100 ns, the RMSD of the compounds with respect to the docking pose was calculated with VMD [58] after aligning the RBD alpha carbon atoms to the initial conformation. The MD replicas in which the final ligand RMSD was lower than 7 Å were extended to 500 ns, for a total of 13 replicas involving 9 compounds (Table S1). The RMSD of the drugs was then computed again, highlighting the compounds still bound to the RBD pocket (Table S1): cefsulodin, cromoglycate, nafamostat, nilotinib, penfluridol, and radotinib.

For these six compounds, the binding energies were evaluated by computing the molecular mechanics energies combined with the generalized Born and surface area continuum solvation were computed with the MMPBSA.py [59] script, from the AmberTools17 suite (The Amber Molecular Dynamics Package, at https://ambermd.org/) employing the default settings and considering the last 50 ns of simulations out of 500 ns (1 frame every 100 ps of the simulation).

Classic MD simulations of ACE2:RBD complex

The coordinates of the spike protein RBD (residues C336-F515 from one monomer of the S protein) and human ACE2 (residues I21-N599) were retrieved from one monomer present in the PDB entry 6M17 [25]. After removing the N-acetylglucosamine residues (which are not located at the interface between the proteins) and the zinc ion, hydrogen atoms were added using pdb2pqr [40] and the propka [41] software. The resulting system was prepared for simulations with the Amber14SB [42] force field by creating a 124 Å × 117 Å × 160 Å TIP3P [43] solvated box (a padding of 15 Å along the x, y, and z dimensions was considered). Overall charge neutrality was maintained by adding 21 Na+ ions.

ACEMD [44] was used for both equilibration and 200 ns of MD productive simulation. Restraints were applied to protein alpha carbon atoms and gradually released throughout 2 ns of equilibration before starting the MD production.

MMGBSA binding energies were computed as described in Sect. “Post Docking molecular dynamics simulations” considering 1000 frames (1 frame every 200 ps of the simulation).

Supervised MD simulations of the ACE2:RBD complex

The supervised molecular dynamics (SuMD) is an adaptive sampling method for speeding up the simulation of binding [37, 39] and unbinding processes [38, 60]. Sampling is gained without the introduction of any energetic bias, by applying a tabu–like algorithm to monitor the distance between the centers of mass (or the geometrical centers) of the ligand and the predicted binding site or the receptor. However, the supervision of a second metric of the system can be considered [60]. A series of short unbiased MD simulations are performed, and after each simulation, the distances (collected at regular time intervals) are fitted to a linear function. If the resulting slope is negative (for binding) or positive (for unbinding) the next simulation is started from the last set of atom velocities and coordinates, otherwise, the simulation is discarded and a new one is started from the last productive coordinates by randomly assigning the atomic velocities.

The RBD (residues C336-F515 from one monomer of the S protein) and ACE2 (residues I21-N599) from one monomer present in the PDB entry 6M17 [25] were placed about 25 Å away from each other and prepared for simulations as reported in Sect. “Classic molecular dynamics simulation of the spike protein RBD”. Restraints were applied to protein alpha carbon atoms and gradually released throughout 2 ns of equilibration before starting the MD production simulation. Four SuMD replicas were collected by monitoring the distance between the centroid of RBD residue Q493 and the centroid of ACE2 residues K31, E35, throughout successive 2 ns-long windows. In the PDB structure 6M17, the interactions between Q493 (RBD) and K31, E35 (ACE2) are in a central position at the interface between proteins. The supervision was iterated until the distance was lower than 7 Å, then 200 ns of classic (unsupervised) MD was performed.

MMGBSA binding energies were computed as described in Sect. Post Docking molecular dynamics simulations considering 1 frame every 100 ps of the simulation.

SuMD simulations of the RBD:cefsulodin:ACE2 ternary complex

The RBD:cefsulodin complex resulting from the post docking MD (500 ns, Sect. “Classic MD simulations of ACE2:RBD complex”) was placed about 25 Å away from ACE2 (residues I21-N599). The resulting system was prepared for simulations as reported in Sect. Classic MD simulations of ACE2:RBD complex. In analogy with SuMD simulations in the absence of cefsulodin bound to RBD (Sect. SuMD simulations of the RBD:cefsulodin:ACE2 ternary complex), four SuMD replicas were collected by monitoring the distance between the centroid of RBD residue Q493 and the centroid of ACE2 residues K31, E35, throughout successive 2 ns-long windows. The supervision was iterated until the distance was lower than 7 Å, then 200 ns of classic (unsupervised) MD was performed.

Molecular mechanics energies combined with the generalized Born and surface area continuum solvation were computed with the MMPBSA.py [59] script, from the AmberTools17 suite (The Amber Molecular Dynamics Package, at https://ambermd.org/) employing the default settings and considering 1 frame every 100 ps of simulation.

Results and discussion

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 [30, 61], targeting the S glycoprotein could represent a promising way to treat and, most importantly, prevent COVID-19. Indeed, the ACE2:RBD binding interface (Fig. 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, can enhance the docking predictivity [62] 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 [63].

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 with respect 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 with respect 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 with respect 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, Fig. 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 MMGBSA binding 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 (Fig. 2a). Direct interactions with cromoglycate and nafamostat, for example, folded RBD site 3 towards the putative binding pocket occupied by the ligand. 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 (Fig. 2b). Within RBD site 1, the ligands formed hydrophobic π–π 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).

Fig. 3

Influence of RBD-bound cefsulodin on the intermolecular recognition between RBD and ACE2 during binding. Top panels plots show the RMSD with respect to PDB 6M17 of the whole RBD (cyan line) and the RBD residues located at the binding interface with ACE2 (green line), respectively; the red and black lines indicate the MMGBSA binding energy of RBD:ACE2 complex and cefsulodin, respectively. Bottom panels: RBD:ACE2 intermolecular contacts potted on the surface of the proteins; the RBD pocket, which was occupied by cefsulodin at the beginning of the simulations, is indicated in green in (c) and (d). a classic MD simulation of the RBD:ACE2 complex from PDB 6M17; b RBD:ACE SuMD binding; the dimer reached the 6M17 bound conformation (RMSD ≈ 2 Å) in less than 20 ns. c)) RBD:cefsulodin:ACE ternary complex SuMD binding replica 1; despite reaching the ACE2 surface, RBD was not stabilized due to the co-presence of cefsulodin on protein interfaces; d RBD:cefsulodin:ACE2 ternary complex SuMD binding replica 2; cefsulodin was displaced by ACE2 after 140 ns of simulation, However, the RBD did not reach the experimental conformation within 230 ns (RMSD ≈ 10–15 Å).). Vertical dashed lines indicate the start of 200 ns of classic MD after SuMD; the three RBD sites responsible for binding ACE2 are indicated with superscripts on the RBD residue names (S1, S2, and S3)

Fig. 2

Binding conformation of cefsulodin, cromoglycate, nafamostad, nilotinib, penfluridol, and radotinib after 500 ns of post-docking MD simulations. a comparison between the last MD frames highlighting the different RBD orientation; cromoglycate and nafamostat strongly affected the RBD site 3. b schematic representation of the interaction features characterizing the RBD pocket as determined by the binding modes of the six proposed ligands; the three sites responsible for binding ACE2 are indicated with dashed lines (site 1 in red, site 2 in blue, and site 3 in green respectively). c-h binding modes of the ligands after 500 ns of post-docking classic MD. Hydrogen bonds are shown as dashed red lines, while hydrophobic contacts are depicted as transparent cyan surfaces. The 2D structures of the compounds are reported for clarity. c) cefsulodin; d cromoglycate: e nafamostat; f nilotinib; g penfluridol; h radotinib. The three RBD sites responsible for binding ACE2 are indicated with superscripts on the RBD residue names (S1, S2, and S3)

Cefsulodin, cromoglycate, nafamostat, nilotinib, penfluridol and radotinib proposed binding modes

Cefsulodin (Fig. 2c) is a third-generation β-lactam cephalosporin antibiotic with a spectrum of activity restricted to Pseudomonas aeruginosa and Staphylococcus aureus [64]. 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 (MMGBSA binding 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 (Fig. 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 β-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 (Fig. 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, Fig. 2d, Fig. 2a, Video S2). The drug formed two persistent hydrogen bonds between the carboxylates and N501, Y453 side chains (Fig. 2d). Hydrophobic interactions involved Y505, Y495, and Y453 (Fig. 2d).

Nafamostat (Fig. 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 [65, 66]. 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 (Fig. 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 (Fig. 2e, Video S3).

Nilotinib (Fig. 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 [67] via ABL-mediated cytoskeletal rearrangement [68, 69]. During the 500 ns MD simulation (Video S4), nilotinib, which remained stably bound to the RBD pocket (MMGBSA binding 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 (Fig. 2f, Video S4). Hydrogen bonds were established between the drug and the N501 and R403 side chains (Fig. 2f, Video S4).

Penfluridol (Fig. 2g) is a long-acting antipsychotic agent [70]. 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 (Fig. 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 [71], which may indicate a common molecular mechanism.

Radotinib (Fig. 2h) is a second-generation BCR-ABL tyrosine kinase inhibitor with therapeutic indication for patients resistant or intolerant to imatinib [72]. Despite the high structural similarity with nilotinib (from which it differs only by the pyridine terminal ring, Fig. 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, Fig. 2h), radotinib displayed lower stability (MMGBSA binding 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 [63, 73] on the basis of the ranking provided by docking scoring functions.

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 (Figs. 1, 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, Fig. 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 (Fig. 2c), and the low toxicity [74].

The first batch of SuMD simulations sought 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, Fig. 3a, b). A further SuMD replica reproduced the complex with good agreement to the cryo-EM structure (Figure S6a), while two other 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 (Fig. 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) [34, 75, 76].

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, Fig. 3c, d, Figure S7). During two simulations, the RMSD of the RBD with respect to the experimental conformation increased, indicating non-productivity of the simulated binding event due to stochastic reasons (e.g. the supervision on the distance between the RBD residue Q493 and the ACE2 residues K31 and E35 failed to move the proteins closer), 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 Å (Fig. 3c, d), consistent with a possible productive contact between the proteins. In one of the two productive SuMD simulations (Video S8, Fig. 3c), despite the formation of a ternary complex involving the RBD, cefsulodin and ACE2, the MMGBSA binding energy between the RBD and ACE2 oscillated around zero, indicating the formation of a sub-optimal binding interface between the two proteins (Fig. 3c). Following a slight conformational rearrangement, cefsulodin remained bound to the interface between proteins for the whole simulation (MMGBSA binding energy in the − 20– − 40 kcal/mol range, Fig. 3c, cf. MMGBSA binding energy in the range of − 32– − 50 kcal/mol in Figure S3). The RBD formed limited contacts with ACE2 (residues D38, Y41, and K353, Fig. 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, Fig. 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 MMGBSA binding energy ≈ − 20 kcal/mol, against ≈ − 45 kcal/mol as computed for the experimental complex, Fig. 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 of 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.

Conclusion

We performed the VS and extensive MD and SuMD simulations of FDA-approved drugs. According to our computational protocol, six compounds (cefsulodin, cromoglycate, nafamostat, nilotinib, penfluridol, and radotinib) are predicted to bind the SARS-CoV-2 S glycoprotein RBD. Cefsulodin (a cephalosporin antibiotic) and nilotinib (a kinases inhibitor) were the most efficient binders in silico.

We further evaluated cefsulodin as a potential disruptor of the RBD:ACE2 complex by simulating the binding process in the absence and presence of the drug bound to the RBD. While cefsulodin hindered the formation of the complex in the simulated conditions, it also appears that ACE2 altered the affinity of cefsulodin for the RBD. According to SuMD, the presence of cefsulodin in the RBD pocket during the approach to ACE2 could modify the overall binding path between the two proteins, hampering the formation of the experimental arrangement observed in the cryo-EM structure. This is promising for the success of drug repurposing strategies targeting the early step of SARS-CoV-2 infection. However, our simulations suggested also the alternative scenario in which the dynamic protein–protein interface formed during the binding between RBD and ACE2 creates a supramolecular environment that modifies the affinity of the putative disruptor for the RBD pocket, with the final effect of partially or completely displacing it. This double outline raises the question whether disruptors should be in silico designed considering only complementarity and affinity for the RBD pocket or designed to stabilize unproductive metastable ternary complexes (e.g. involving the RBD, ACE2 and the ligand). Future work will aim to compare several other potential RBD binders (i.e. nilotinib) to deliver insights about what set of interactions with the RBD are more prone to stabilize the ligand during the ACE2 recognition.

If repurposing strategies fail to deliver an actual disrupter of the protein, it would be necessary to further investigate the binding process between the spike protein and ACE2, characterizing intermediate complexes that could be drugged to prevent the formation of the one which is productive for the cell invasion. These results extend the computational toolkit for the rational discovery of PPIs disruptors to adaptive or enhanced sampling methods, like SuMD, able to simulate the nonequilibrium formation of homo- and hetero dimers.

Below is the link to the electronic supplementary material.

Supplementary file1 (pdf 4867 kb)

Supplementary file2 (mpg 5714 kb)

Supplementary file3 (mpg 5348 kb)

Supplementary file4 (mpg 5484 kb)

Supplementary file5 (mpg 5086 kb)

Supplementary file6 (mpg 5560 kb)

Supplementary file7 (mpg 4978 kb)

Supplementary file8 (mp4 26467 kb)

Supplementary file9 (mpg 6500 kb)

Supplementary file10 (mpg 5112 kb)