Small Molecules to Destabilize the ACE2-RBD Complex: A Molecular Dynamics Study for Potential COVID-19 Therapeutics.

The ongoing COVID-19 pandemic has infected millions of people, claimed hundreds of thousands of lives, and made a worldwide health emergency. Understanding the SARS-CoV-2 mechanism of infection is crucial in the development of potential therapeutics and vaccines. The infection process is triggered by direct binding of the SARS-CoV-2 receptor-binding domain (RBD) to the host cell receptor, Angiotensin-converting enzyme 2 (ACE2). Many efforts have been made to design or repurpose therapeutics to deactivate RBD or ACE2 and prevent the initial binding. In addition to direct inhibition strategies, small chemical compounds might be able to interfere and destabilize the meta-stable, pre-fusion complex of ACE2-RBD. This approach can be employed to prevent the further progress of virus infection at its early stages. In this study, Molecular docking is employed to analyze the binding of two chemical compounds, SSAA09E2 and Nilotinib, with the druggable pocket of the ACE2-RBD complex. The structural changes as a result of the interference with the ACE2-RBD complex are analyzed by molecular dynamics simulations. Results show that both Nilotinib and SSAA09E2 can induce significant conformational changes in the ACE2-RBD complex, intervene with the hydrogen bonds, and influence the flexibility of proteins. Moreover, essential dynamics analysis suggests that the presence of small molecules can trigger large-scale conformational changes that may destabilize the ACE2-RBD complex.

INTRODUCTION

The SARS-CoV-2 is the newest member of coronaviruses, a family of single-stranded RNA based viruses able to infect mammals such as bats and humans. While the SARS-CoV-2 preserves almost 70% of the structure of the previous SARS-CoV virus, the differences are enough to affect the Pathogenicity. Several studies show that the SARS-CoV-2 has a higher binding affinity with its host cell receptor, Angiotensin-Converting Enzyme 2 (ACE2). The higher risk of virus initial recognition can justify the higher infectivity of the novel coronavirus. However, the overall path of infection is similar to the SARS-CoV virus. The infection process is initialized by virus recognition through binding of the receptor-binding domain (RBD) of spike glycoproteins with ACE2 receptors. The recognition step is followed by the cleavage and activation processes induced by protease enzymes such as TMPRSS2 which paves the way for virus fusion and entry into the cell membrane. Since the initial days of the emergence of the novel coronavirus, huge scientific efforts have been made to target various stages of the infection process from inhibiting the initial host-cell binding process to intervene in the virus reproduction pathway. Here we focus on the early stages of infection and virus recognition process by host cell receptors.

The crystal structure of RBD-ACE2 complex with the resolution of 2.45 Å (PDB:6M0J) was determined by Lan et al.(1). Moreover, Yan et al.(2) presented the cryo-electron microscopy structures of full-length ACE2 receptors with and without the RBD (PDB:6M17). Various strategies have been evaluated to block RBD-ACE2 interactions. Monolocal antibodies show promising results in deactivating SARS-CoV-2 spike proteins (3–5). Moreover, ACE2-based peptides can be designed to compete with real cellular receptors and reduce the risk of infection (6–8). In the absence of approved drugs and vaccines, screening techniques have been widely employed in order to repurpose approved drugs and find chemical compounds that are able to inhibit various proteins that are essential for the infection process. Adedeji et al (9) carried out the screening of a large library of pharmacologically active small molecules and suggested three compounds with different mechanisms of action to inhibit the infection process of SARS-CoV. Interestingly, they reported that N-[[4-(4-methylpiperazin-1-yl)phenyl]methyl]-1,2-oxazole-5-carboxamide or SSAA09E2 has a novel inhibition mechanism that allows it to interfere with ACE2-RBD interactions and block the early stages of infections. However, details of this novel mechanism and its effectiveness in the case of the SARS-CoV-2 virus are yet to be understood completely. Wei et al. (10) employed virtual screening to evaluate the binding energy of approximately 15000 drugs with the RBD. Based on their analysis, among all the drugs listed in the DrugBank database (11), Digitoxin and Nilotinib have the best docking scores with the RBD. In addition to virtual screening techniques, all-atomistic (12, 13) and coarse-grained (14–16) molecular dynamics methods are effective computational tools for multiscale modeling of virus infection process and computer aided drug design (17). Deganutti et al. (18) employed structure-based screening and supervised molecular dynamics to find potential drugs that can bind to the druggable pockets of RBD and inhibit the binding. They reported Cefsulodin and Nilotinib as the potential disruptors of the RBD–ACE2 interface. Nilotinib is a member of Abelson (Abl) tyrosine-protein kinase (Abl-TK) inhibitor drugs that can potentially affect the infection process through inhibiting the virus fusion prior to hemifusion (19). In addition, these drugs can reduce the risk of infection through Abl-mediated cytoskeletal rearrangement and interfering with actin dynamics (20). The in-vitro observations of positive impacts of Abl-TK inhibitors have recently been strengthened by in-vivo studies and medical surveys. For instance, Foa et al.(21) and Breccia et al. (22) reported a lower percentage of COVID-19 infection among the Philadelphia-positive acute lymphoblastic leukemia and chronic myeloid leukemia patients under treatment by Abl-TK inhibitors.

In this study, we focused on the interference mechanism of small molecules that can potentially destabilize the ACE2-RBD complex. The disruption of protein-protein interactions can be a strategy for potential drugs to interfere with the early stages of infection. Based on in-vitro and virtual screening studies, we chose SSAA09E2 and Nilotinib to study their potential intervention effects. Instead of direct inhibition of ACE2 or RBD, this study aims to understand the structural changes induced by binding of SSAA09E2 and Nilotinib into the druggable pocket at the central part of the ACE2-RBD interface. Such intervention process can shed light on the mechanism of action of SSAA09E2 and Abl-TK drugs. Molecular docking is employed to find the most favorable binding poses. Molecular dynamics is then utilized to analyze the structural changes of the ACE2-RBD complex. We showed that both of these chemical compounds are able to make stable bonds with the ACE2-RBD complex and reduce the infection risk by interfering with the ACE2-RBD interactions and triggering large scale conformational changes.

METHODS

The crystal structure of SARS-CoV-2 receptor-binding domain complexed with ACE2 enzyme was obtained from protein data bank (PDB:6M0J) (1). The chemical structures of Nilotinib and SSAA09E2 are retrieved from the PubChem database (23). Avogadro software is used to prepare the ligand structures and adding hydrogen atoms (24). The CHARMM-GUI ligand reader and modeler module (25) is utilized to prepare parameter and topology files of ligand molecules based on the CGenFF (26) method. Molecular docking simulations are performed by Autodocktools 4 (27). A grid with the size of 60 Å × 60 Å × 60 Å centered on the druggable pocket at the interface of the ACE2-RBD complex is used for docking simulations. The docking process is carried out with the Lamarckian Genetic algorithm (28, 29). The best docking pose is used for further molecular dynamics simulations. The system topology is made by CHARMM-GUI solution builder module by using the CHARMM36(m) all-atom forcefield (30, 31) and TIP3P water model (32). An octahedral computational box with 12 Å distance between the edges is employed for the protein solvation. The system is neutralized by adding Potassium(K) and chloride(Cl) ions with the final concentration of 0.15 M. Systems are initially minimized by the steepest descent algorithm with 5000 steps. The minimization is followed by 125 ps of equilibration simulation in the canonical ensemble (NVT). Linear Constraint Solver (LINCS) algorithm is employed to constrain the hydrogen bonds(33). The Verlet cutoff scheme with the cut-off radius of 1.2 nm and the Particle Mesh Ewald method(PME) are used to compute long-range electrostatics. All the systems are simulated at the temperature of 310 K that is maintained by Nosé–Hoover temperature coupling method with the time constant for coupling (tau-t) of 1 ps (34). The 200 ns of production simulations are performed at the NPT (constant number of particles, pressure, and temperature) ensemble at the pressure of 1atm which is controlled by the Parinello-Rahman method with the time constant for pressure coupling (tau-p) of 5 ps. Gromacs 5.1 is employed for all the molecular dynamics simulations and trajectories are visualized by PyMol and VMD packages. g_RMS, g_RMSF, g_gyrate, g_hbond, g_covar, and g_anaeig packages of Gromacs software are employed to analyze the root mean square deviations (RMSD), root mean square fluctuations (RMSF), radius of gyration (R), number of hydrogen bonds, calculation of eigenvalues and eigenvectors, and 2D projection of eigenvectors respectively.

RESULTS AND DISCUSSION

Our docking simulations show that both Nilotinib and SSAA09E2 have high binding affinities with the druggable pocket (as shown in figure 1.a) of the ACE2-RBD complex . Figure 1.b and figure 1.c show the two-dimensional diagram of ligand-protein interactions (diagram is made by LigPlot+ software (35)). It can be seen that at this pocket, Nilotinib has non-bonded interactions with residues Lys26, Thr27, Leu29, Asn33, His34, Glu37, ASP38, Lys353, Arg393 of the ACE2 protein which is shown by chain A in the diagram. Two potential hydrogen bonds between the Glu37 and Gln96 are identified between the Nilotinib and ACE2. The RBD (shown by chain E) non-bonded interactions with Nilotinib are through Arg 403, Ser494, Tyr495, Gly496, and Tyr505. The best docking score for the pose shown in figure 1 is −8.34 kcal/mol. Figure 1.c shows the interaction diagrams and best docking pose of SSAA09E2. Residues Asn33, His34, Glu35, Glu37, Asp38 of ACE2 have non-bonded interactions with SSAA09E2. Moreover, the Lys353 is able to make a hydrogen bond with the ligand. In contrast to Nilotinib, SSAA09E2 can form a hydrogen bond with RBD by interacting with Gly496. Arg403, Tyr449, Tyr453, Ser494, Tyr495, Gly496, Tyr505 from the RBD have non-bonded interactions with SSAA09E2. Three dimensional poses of drugs in the ACE2-RBD complex are shown in Figure 1.d and 1.e. The best docking pose has a score of −8.04 kcal/mol which is close to the case of Nilotinib with the score −8.34 kcal/mol. Our docking results show that both Nilotinib and SSAA09E2 can potentially interfere with some of the important ACE2-RBD interactions such as Tyr449–Asp38, Tyr453–His34, and Gln493–Glu35 (1, 2, 36, 37).

Figure 1:

(a) Electrostatic surface of ACE2-RBD complex and the pocket at the interface of ACE2 and RBD, Ligand-protein interactions for (b) Nilotinib, (c) SSAA09E2, and best docking poses of (d) Nilotinib, and (e) SSAA09E2.

To analyze the stability and conformational changes of protein structures, 200 ns of MD simulation is performed for the control group without drugs, and ACE2-RBD complexes bound with Nilotinib and SSAA09E2, respectively. The RMSD is a measure of protein drift from a reference structure. To exclude less meaningful fluctuations of often flexible side chain atoms, we only focus on the RMSD of backbone atoms. Figure 2.a shows the RMSD of the ACE2 proteins in various cases. In all simulations, the ACE2 protein reaches an equilibrium after 25ns of simulation. However, the SSAA09E2 has a jump at 50ns followed by a new equilibration state at a higher RMSD. The average RMSD is 2.79 Å, 2.73 Å, and 3.50 Å for the control, Nilotinib, and SSAA09E2 cases respectively. Thus, adding SSAA09E2 results in the greatest structural changes in comparison with the initial reference structure. The RMSF shows the time-averaged fluctuations of individual residues. As shown in figures 2.b and 2.c, both drugs can induce some changes in the RMSF of key residues. For instance, the RMSF of residues number 20–80 (figures 2.d), which are at the interface of ACE2 interaction with RBD, is increased for both drugs. The higher estimated flexibility of interface residues can lead to higher entropy penalty that negatively influences the binding energy, especially at higher temperatures (38). The average RMSF over all residues is 1.43 Å, 1.52 Å, and 1.77 Å for the control, Nilotinib, and SSAA09E2 cases respectively. Thus, adding drug molecules leads to an increase in the flexibility of the ACE2 protein and its random fluctuations.

Figure 2:

(a) RMSD of the ACE2 protein, (b) RMSF of the ACE2 protein, (c) RMSF of residues number 20–80 of the ACE2 protein, (d) snapshot of the ACE2 protein (blue), RBD (red) and interface residues number 20–80 (yellow)

Figure 3.a shows the RMSD of the RBD. For the control case, the RMSD converges at ~100ns after a jump at ~80ns. However, a mild positive slope can be seen in the Nilotinib case which seems to be flattened after ~150ns. The final RMSD for the RBD with Nilotinib is close to the control case. The RBD has the most unstable behavior in the SSAA09E2 case in which the RMSD is fluctuating around ~4 Å. The average RMSD values over 200 ns of the simulation are 2.40 Å, 2.60 Å, and 3.92 Å for the control, Nilotinib, and SSAA09E2 cases respectively. However, by analyzing the RMSF plot of RBD, it can be seen that the higher instabilities in the SSAA09E2 case are mainly stemmed from the large fluctuations of the terminal residues. Figure 3.b shows the RMSF plots for the RBD protein. The presence of the SSAA09E2 molecule has significant effects on at least four regions of the RBD protein. The most notable increase in the flexibility is at the terminal residues His519-Gly526 in which the RMSF is much larger in the SSAA09E2 case. While the whole spike protein is not simulated in our simulations, very large fluctuations and flexibility of the terminal residues in the case of SSAA09E2 may be effective in triggering conformational changes in the rest of the spike protein structure. On the other hand, the RMSF of the Nilotinib case is close to the control system in this region. Moreover, adding SSAA09E2 leads to a substantial drop in the flexibility of the loop region between Ser477-Gly485. This can be effective in the ACE2-RBD binding process. For instance, the lower flexibility of this loop can affect the formation of hydrogen bonds between the residues Asn487 of the RBD and residue Gln24 of the ACE2 protein (1) (Figure 3.c). On the other hand, residues Asn437-Leu441 have higher flexibility in the SSAA09E2 case. In comparison with the control case, less significant increases in the RMSF can be observed for both Nilotinib and SSAA09E2 at residues Ser359–369 and Thr333-Glu340. The average of RMSF over all residues is 1.70 Å for the SSAA09E2 case which is significantly larger than 1.28 Å and 1.36 Å for the control and Nilotinib cases respectively, mainly due to huge fluctuations at the terminal residues.

Figure 3:

(a) RMSD of the RBD, (b) RMSF of the RBD, (c) Residues with more significant changes in the RMSF are shown with yellow sticks, Gly485 and Gln24 are shown by gray sticks (d) RMSD of ligand molecules

Figure 3.d illustrates the RMSD of the ligand molecules at the binding pocket. SSAA09E2 is stable during the 200 ns simulation with the average RMSD of 1.53 Å. On the other hand, Nilotinib shows more unstable behavior by fluctuating between the RMSD values of ~2–4 Å with the average RMSD of 2.84 Å. Nilotinib is a larger molecule in comparison with SSAA09E2. While at the best docking pose almost half of the Nilotinib molecule is in the druggable pocket, another half of the molecule can not make stable bonds with ACE2 or RBD and is freer to fluctuate (see the supporting movie S1). This leads to higher total deviations from the reference structure. These higher fluctuations of the Nilotinib may help this molecule to interfere with the hydrogen bonds between ACE2 and RBD. Despite higher RMSD, visual analysis of trajectory shows that the Nilotinib molecule is fixed at the pocket during the 200 ns of simulation and does not escape from its initial location. Figure 4 shows the radius of gyration (R) of the ACE2 and RBD proteins in all simulations. R is a measure for the compactness of proteins and shows the stability of protein folding. It can be seen that in all cases the proteins are stably folded and the differences in the compactness are negligible. Thus, drugs do not influence the consistent shape and size of ACE2-RBD complex during the 200 ns of simulation.

Figure 4:

R of the (a) ACE2, and (b) RBD proteins.

Hydrogen bonds have a significant role in the strength of protein-protein interactions (39). A recent study shows that the average number of hydrogen bonds formed between the RBD and cell receptors can be almost two times more for the novel SARS-CoV-2 in comparison with the SARS-CoV (40). This can be a reason for the greater binding energy of ACE2-RBD and more infectivity of the SARS-CoV-2 virus. We studied the hydrogen bonds between the ACE2 and RBD proteins in various cases. To make plots more clear, in addition to the number of hydrogen bonds at each frame, running averages of every 50 frames are plotted. As illustrated in figure 5.a, the SSAA09E2 and Nilotinib have different effects on the number of hydrogen bonds between the ACE2 and RBD proteins. In the case of SSAA09E2, more hydrogen bonds form between the two proteins. On the other hand, Nilotinib weakens the structural integrity by reducing the number of hydrogen bonds. The average number of hydrogen bonds during the 200 ns of simulation is 7, 5.4, and 8.5 for the control, Nilotinib, and SSAA09E2 cases respectively. Thus, on average, the Nilotionib case has ~1–2 fewer hydrogen bonds in comparison with the control system. Moreover, we analyzed the number of hydrogen bonds between ligand molecules and ACE2 or RBD. figures 5.b to 5.f show the number of hydrogen bonds formed between SSAA09E2 and Nilotinib with ACE2 and RBD proteins. Both ligands form more hydrogen bonds with the ACE2 protein rather than the RBD. On average, the Nilotinib has ~1.35 hydrogen bonds during the 200 ns of simulation with the ACE2 protein which is slightly greater than ~1.10 bonds made between the SSAA09E2 and ACE2. On the other hand, the average number of hydrogen bonds of SSAA09E2 with the RBD is ~0.30 which is larger than ~0.09 for the Nilotinib case. Analysis of hydrogen bonds shows that Nilotinib and SSAA09E2 may have two different mechanisms of action to interfere with the ACE2-RBD complex. Fluctuations of the Nilotinib molecule intervene with the hydrogen bonds formed between ACE2 and RBD leading to weaker complex integrity. However, SSAA09E2 does not decrease the ACE2-RBD hydrogen bonds and may even increase them.

Figure 5:

Number of Hydrogen bonds for the (a) control ,(b) Nilotinib , and (c) SSAA09E2 and rolling average of hydrogen bonds for (d) Control, (e) Nilotinib, and (f) SSAA09E2.

We also employed the Principal Component Analysis (PCA) or essential dynamics to identify dominant motion modes of the ACE2 and RBD in various cases. The covariance matrix of fluctuations is calculated for backbone atoms and diagonalized to find the eigenvalues and eigenvectors. The dominant motions of each system are then identified by projecting the trajectory along the calculated eigenvectors. Usually, a few numbers of low-frequency eigenvectors are responsible for the majority of overall fluctuations. As shown in figure 6, our results indicate that the first three eigenvectors account for 53%, 71%, 69% of overall fluctuations in the control, Nilotinib, and SSAA09E2 case respectively.

Figure 6:

Eigenvectors and corresponding eigenvalues for various systems.

Figure 6 also illustrates that the first three eigenvalues are larger for the cases with Nilotinib and SSAA09E2 which means that systems with ligand molecules have greater movements along the first three eigenvectors. Figure 7 shows the conformational space sampled by the ACE2-RBD complex in various cases. Each point shows a frame of 200 ns trajectory projected onto the two-dimensional space of the first two eigenvectors. While there is an overlapping area for all cases, figure 7 clearly shows the effect of ligand molecules on intensifying dynamics of the ACE2-RBD complex which is consistent with previous results of the RMSD and RMSF. The largest covered area is for the Nilotinib case which shows more powerful random fluctuations. Moreover, the area covered by the ACE2-RBD complex bound with SSAA09E2 illustrates the higher flexibility of this case in comparison with the control case. For a better analysis of the dominant dynamics, porcupine plots (41) are employed to show characteristics of dynamic fluctuations for the first three eigenvectors. In figure 8, vectors indicate the motion of each backbone atom along the eigenvector direction. The length of each arrow represents the magnitude of the corresponding eigenvalue. As expected, the largest motions happen at the loop regions. For the first eigenvector, the overall motion of the ACE2-RBD complex bound with the SSAA09E2 molecule is similar to the control case. However, unlike the control case, for the SSAA09E2 case, residues S443 to N448 of the RBD (shown by the black circle in figure 8) are getting away from the ACE2 protein. The difference between the dominant motions is more significant in the Nilotinib case. In contrast to the control case, the ACE2 and RBD proteins bound with Nilotinib moves away from each other by a counter-clockwise rotation in RBD and a clockwise rotation in ACE2. In addition to different motions, the magnitude of eigenvalues is greater in the systems with ligand molecules. Results of essential dynamics analysis clarify the difference in the mechanisms of action of SSAA09E2 and Nilotinib. While SSAA09E2 is changing the motion dynamics especially at the terminal residues and loops close to the ACE2 interface, Nilotinib changes the large-scale motions of both proteins.

Figure 7:

Dynamic conformations projected onto the two first principal vectors for (a) control, (b) Nilotinib, (c) and SSAA09E2 cases. All cases are compared in part (d).

Figure 8:

Porcupine plots of characteristic dynamic fluctuations for the three lowest-frequency principal components

CONCLUSION

Molecular docking and molecular dynamics methods are employed to study the intervention and destabilization mechanisms of Nilotinib and SSAA09E2. Our structural analysis disovers on two different mechanisms of action. The stable binding of SSAA09E2 at the pocket between ACE2 and RBD leads to higher flexibility in the residues Ser438 to Ser443 and terminal residues of the RBD, less flexibility at the residues Ser477 to Gly484, and higher RMSD of ACE2 and RBD from the reference structure. Moreover, while the PCA analysis shows that the dominant motions of SSAA09E2 are similar to the control case, some large scale motions such as the outward rotation of residues S443 to N449 can affect the ACE2-RBD binding. On the other hand, the Nilotinib is less stable at the binding pocket and half of its structure which is posed parallel to the ACE2 interface with RBD can freely fluctuate out of the druggable pocket of the ACE2-RBD complex. However, Nilotinib is able to weaken the interactions between the ACE2 and RBD proteins by decreasing the total number of hydrogen bonds. More interestingly, our essential dynamics analysis illustrates that greater conformational spaces are sampled by the ACE2-RBD complexes bound with Nilotinib and SSAA0E2. In other words, both drugs are destabilizing the ACE2-RBD complex by triggering more powerful random fluctuations and larger-scale motions. In the case of Nilotinib, filtering the trajectory by the first three eigenvectors shows that the most dominant motions of ACE2 and RBD proteins are completely different from the control system. In contrast to the control system, Nilotinib induces a divergent large scale motion that rotates ACE2 and RBD in opposite directions. Our structural analysis broadens the understanding of small molecules destabilization and intervention mechanisms that would be an effective strategy to prevent infection at its early stages.