Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a rapidly progressing pandemic disease, which developed in China at the end of 2019 and then spread rapidly around the world [1, 2]. This virus is highly contagious, and according to WHO statistics in June 2022 it infected 532 M people and more than 6.3 M deaths in total [3]. Accordingly, COVID-19 had put all researchers in an alert situation, to gain a deep understanding of viral pathophysiology and find effective therapies.

SARS-CoV2 may cause moderate to severe diseases ranging from colds to organ disturbances and acute respiratory distress syndrome (ARDS) [4, 5]. Treatment ranges from symptomatic therapy to oxygen therapy to mechanical ventilation [6]. Moreover, Food and Drug Administration (FDA) granted emergency use authorization (EUA) of Paxlovid™ (nirmatrelvir/ritonavir tablets) against SAR-COV-2 [7]. In addition, many control strategies have emerged since 2020, including vaccination, and some drugs such as remdesivir, dexamethasone, plasma therapy and monoclonal antibodies have shown effective results in clinical trials [4, 8]. Although the field of research is actively continuing to understand the molecular basis of this virus to a clear and effective targeted treatment protocol [8].

Coronavirus contains single-positive RNA with genome sizes 26–32 kilobases and belongs to the subgenus Sarbecovirus of the family Coronaviridae member of the genus Betacoronavirus [1, 4, 5, 9]. The SARS-CoV-2 genome contains open reading frames (ORFs) that encode two types of proteins: structural and non-structural proteins (nsps) that play a vital role in viral genome replication and gene transcription [4, 10, 11]. Spike glycoprotein (S), envelope protein (E), membrane protein (M) and nucleocapside protein (N) are the main structural proteins of SARS-CoV-2. While RNA-polymerase-dependent RdRp (RdRp), also known as nsp12, is a structural protein that represents an attractive drug selective target in SARS-COV-2, as it provides drug selectivity as it is absent in its human counterparts [12,13,14,15]. The size of RdRp is 240–250 KD, and its architecture consists of a catalytic core containing domains of palm, thumb and fingers with various residues [15, 16]. RdRp plays a central role in the replication and transcription process of SARS-CoV-2 with the help of cofactors nsp7 and nsp8 [15, 16]. The C-terminal domain of RdRp is similar to that of other viruses, but the selectivity of targeting SARS-CoV-2 resides in its specific N-terminal nidovirus domain used for replication [2]. For this, remdesivir and favipiravir are nucleoside analogues that have been used for COVID-19 treatment, since they can inhibit RdRp, more specifically the specific domain of N-terminal nidovirus, and block the activity of nucleotide transferase [2, 10, 17,18,19].

The life cycle of SARS-CoV-2 infection is briefly mediated through its spike proteins to allow its binding to the host cell and then releases its RNA which is then processed by the papain-like protease and the main protease (Mpro) to produce 16 nsps. Released nsps help in viral replication and translation, where released RNA and proteins will infect another human through exocytosis after being collected as prognostic virions [9].

Viral proteases, in particular the main protease (Mpro), are a cysteine protease, which is an important enzyme in the production of structural and non-structural proteins and, therefore, RNA translation by special processing of polyproteins. Blocking of Mpro leads to effective inhibition of viral growth and replication, since viral protease is responsible for infection and the underlying disease associated with SARS-CoV-2 [5, 12, 20, 21]. Mpro SARS-CoV-2 is found in the regions 3264–3569 (containing 306 residues) of polyprotein pp1a and pp1ab [22, 23]. The importance of Mpro is due to the cleavage of polypeptides together with papaine-like protease (PLpro) to form 11 pp1a and 5 pp1ab, which are used in viral replication [5, 6, 9, 22, 24]. Lopinavir and ritonavir are useful inhibitors of Mpro in the treatment of COVID-19 by binding to Thr24, Thr26 and Asn119 Mpro residues [6, 9].

Since drug discovery is a long process that can take up to 10–15 years to approve a new drug, and costs can reach a billion dollars, CDD helps greatly in advancing drug discovery and development as well as cost and time savings [25, 26]. There are several approaches to CDD; there is a drug repurposing approach that is an elegant method to help find new uses for already approved drugs, especially in unexpected emerging diseases such as the COVID-19 outbreaks, which have placed an additional burden on medical researchers to find effective treatment in a short time to save lives [16, 27, 28].

Drug repurposing is considered to be a fast and cost-effective process because long-term trials are not required to approve drug safety [12, 24, 25]. Three repurposing approaches are used, which can be a computational approach, an experimental approach or mixed approaches [25]. Great advances in the development of computational drugs (CDDs) have led to more advanced algorithms that play a huge role in advances in drug discovery processes [6, 29, 30].

Despite great effort in the research literature published in COVID-19 using in silico approaches, no study has combined Mpro and RdRp as dual targets using approved drugs from the ChEMBL database. Consequently, this study aims to find dual target inhibitors against these two targets as promising new candidates for COVID-19 treatment using molecular coupling, MM-GBSA calculations and molecular dynamics.

Materials and methods

All in silico studies were performed on the Maestro v12.8 of Schrödinger Suite and academic Desmond v6.5 by D.E. Shaw Research for molecular dynamics simulations.

Protein preparation

The RNA-dependent RNA polymerase (RdRp, PDB ID: 7BV2) and SARS-CoV-2 Mpro (PDB ID: 5R81) protein structures were downloaded from the RCSB protein databank (PDB). The crystal structure of proteins can have problems such as incorrect bonding orders, missing side chains, missing loops and atoms, so they need to be prepared before docking. The crystal structures of the two proteins were prepared using the Protein Preparation Wizard; in this process, the purpose was to check the protein structures and verify the assignment of bonds/binding orders, the addition of hydrogens, the detection of disulfide bonds, the completion of missing loops or side chains and the rectification of any mislabeled components. In addition, crystal protein structures were subjected to a process of minimization in the Imperf utility using the OPLS3e force field [31]. During this process, heavy atoms in the structures were restricted to relieve torsional tension with a harmonic potential of 25 kcal/mol while hydrogen remained unrestrained.

Ligand preparation

A ChEMBL library containing more than nine thousand approved drugs was used for repurposing in this study. First, ligand preparation was performed using the Ligprep tool of the Maestro software to obtain three-dimensional structures with correct chirality, ionization states, stereochemistry, tautomers and ring conformations for each input structure. Second, the exact degree of ionization was evaluated at pH 7.0. Also, an OPLS3e force field was used to minimize ligand energy. Finally, the resulting ligand structures were subjected to a docking process.

Molecular docking and MM-GBSA calculations

To measure the binding affinity of receptor-ligand complexes against the two targets (RdRp and Mpro), a molecular docking process was carried out. The Grid-based docking was carried out with Schrödinger Glide, the grid box was produced for each receptor and specific amino acids which interact with the ligand were determined. Later, the docking process was carried out with two approaches: high throughput virtual screening (HTVS) with extra precision (XP) to identify the strongest molecular interactions.

MM-GBSA bonding-free energies were calculated with the Prime of Schrodinger [32]. The binding-free energy of protein and ligand can be calculated using the following equation:

$$\mathrm{\Delta G}(\mathrm{bind})={\mathrm{G}}_{\mathrm{complex}}-{\mathrm{G}}_{\mathrm{protein}}-{\mathrm{G}}_{\mathrm{ligand}}$$

where Gcomplex is the energy of the protein–ligand complex, Gprotein is the energy of the protein and Gligand energy of the ligand.

Molecular dynamics simulation

MD studies were performed using academic Desmond v6.5. MD simulations were to assess binding stability and the interaction profile of the three best scoring drugs in complex with the protein Mpro. Desmond’s System Builder tool was used to build the system for the MD process. The complex was solvated in 9658 TIP3P molecules in an orthorhombic box (10 Å × 10 Å × 10 Å). Salt was added in a specific concentration of Na+ and Cl− charge, as 73.420 mM (total charge + 39) for Na+, while 50.829 mM (total charge −27) for Cl− for the system. Energy minimization was performed utilizing the OPLS3e force field. NPT ensemble was used to equilibrate the minimized system, and the MD simulation was conducted at an atmospheric pressure of 1.013 bar at a temperature of 300 K for 100 nN. Various parameters were used to analyse the MD results, such as root mean square deviation (RMSD), squire mean square fluctuation (RMSF) and 2D ligand protein interaction [33].

Results

The study workflow is summarized in Fig. 1.

Fig. 1
figure 1

The overall works of the study

Full size image

Molecular docking and MM-GBSA calculations

The docking process was first validated by redocking the cocrystal inhibitors remdesivir (RTP) and Z1367324110 in RdRp (PDB ID: 7BV2) and Mpro (PDB ID: 5R81), respectively. The redocking of these native ligands has shown that they bind at the same location as the co-crystal ligands in the original structures and that the root agent square deviation (RMSD) values were 1.9862 and 0.9658 Å respectively (Figs. 2 and 3).

Fig. 2
figure 2

Redocking of co-crystallized inhibitors. Z1367324110 over SARs COV-2 MPro (5R81) (A) Three-dimensional structure (3D) of SARS-CoV-2 main protease complexed with co-crystallized ligand Z1367324110. (B) Superimposition of the co-crystal ligand and re-docked native ligand, the co-crystal ligand shown in white and the redocked ligand in pink (RMSD = 0.9658 Å)

Full size image
Fig. 3
figure 3

Redocking of co-crystallized inhibitors. Remdesivir over RNA-dependent RNA polymerase (RdRp) (7BV2) (A) Three-dimensional structure (3D) of RNA-dependent RNA polymerase complexed with co-crystallized ligand Remdesivir. (B) Superimposition of the co-crystal ligand and re-docked native ligand, the co-crystal ligand shown in white and the redocked ligand in pink (RMSD = 1.9862 Å)

Full size image

The ChEMBL library containing 9923 approved drugs was prepared with Maestro’s LigPrep module, resulting in 23,334 conformations that were then used in the docking process to measure the binding affinity of these drugs against RdRp and Mpro of SARS-CoV-2. HTVS docking mode from Glide was used on RdRp, and the top 24 drugs from HTVS with docking scores ≤ −6.00 kcal/mol were further docked with RdRP with extra precision (XP) docking filter. Drugs with docking results ≤ −7.00 kcal/mol were obtained. The five uppermost compounds were then docked to Mpro, and the docking scores obtained were analysed and compared with their cocrystal ligands remdesivir and Z1367324110. The short-listed drugs, natamycin, nadide, cangrelor, denufosol and fomidacillin, have shown good affinity for both targets with docking results between −7.46 to −10.54 and −6.65 to −9.002 kcal/mol in RdRp and Mpro, respectively, while the reference ligands remdesivir and Z1367324110 had lower scores of −5.978 and −5.969 respectively (Table 1).

Table 1 XP docking scores and MM-GBSA free binding energies for the top five approved drugs with dual activity against RdRp and Mpro
Full size table

The receptor-ligand interactions of the five best hits with the binding pockets of RdRp and Mpro were further investigated in detail. At the active site of RdRp, the reference remdesivir, nadide, denufosol and fomidacillin showed common hydrogen binding interactions with residues Tyr619 and Arg553, while the interaction with Arg555 was observed with the reference ligand together with nadide, cangrelor, denufosol and fomidacillin (Figs. 4 and 5). Furthermore, a hydrogen bond with Asp623 with remdesivir, natamycin and fomidacillin was observed. The interaction with the key residues Arg551 and Asp760 was produced by cangrelor or denufosol, respectively. Additional polar, charge, salt bridge, π cation and hydrophobic interactions were observed as shown in Table 2.

Fig. 4
figure 4

Two-dimensional (2D) interactions of reference compounds in targeted SARs-COV-2 proteins. (A) SARS-CoV-2 main protease with co-crystallized ligand Z1367324110. (B) RNA-dependent RNA polymerase (RdRp) with co-crystallized ligand Remdesivir (PDB id: 7BV2). The hydrogen-bond interactions with residues are represented by a purple dashed arrow directed towards the electron donor. The pi-pi interactions are represented by a green line

Full size image
Fig. 5
figure 5

Two-dimensional (2D) interactions of short-listed drugs in the active sites of SARS-CoV-2 RNA-dependent RNA polymerase (RdRp) (PDB id: 7BV2). (A) Natamycin, (B) Nadide, (C) Cangrelor, (D) Denufosol, (E) Fomidacillin. The hydrogen-bond interactions with residues are represented by a purple dashed arrow directed towards the electron donor. The pi-pi interactions are represented by a green line. The pi-cation interactions are represented by a red line. The salt bridge interactions are represented by a red/purple line

Full size image
Table 2 Intermolecular interaction of tested compounds with RdRp and Mpro enzymes
Full size table

In the Mpro binding site, Glu166 was the main residue in the hydrogen binding of reference ligand Z1367324110 as well as natamycin and denufosol. The latter formed three further hydrogen bonds with Gly143, Arg188 and Gln189, while nadides and fomidacillin have common hydrogen bonds with Thr26, Ser46 and Gly143 (Figs. 4 and 6). Interestingly, cangrelor and fomidacillin showed hydrogen bonds and π-π stacking interactions with the key residue His41. Moreover, the tested drugs build common hydrophobic (Met49, Cys145, Met165, Leu167, Pro168), polar (Thr25, His41, Ser46, Gln189, Thr190) and charged interactions (Glu166, Asp187, Arg188) in the binding pocket (Table 2 and Fig. 6).

Fig. 6
figure 6

Two-dimensional (2D) interactions of short-listed drugs in the active sites of SARS-CoV-2 main protease (Mpro) (PDB id: 5r81). (A) Natamycin, (B) Nadide, (C) Cangrelor, (D) Denufosol, (E) Fomidacillin. The hydrogen-bond interactions with residues are represented by a purple dashed arrow directed towards the electron donor. The pi-pi interactions are represented by a green line. The pi-cation interactions are represented by a red line. The salt bridge interactions are represented by a red/purple line

Full size image

MM-GBSA for the top five drugs has been calculated to augment docking results. The five medicines showed dual target activity and MM-GBSA binding free energy ranged from −29.78 to −73.01 kcal/mol (Table 1).

Molecular Dynamics simulations (MD)

During the docking procedure, the receptor is considered a rigid structure, while ligand and receptor binding is dynamic and receptor flexibility is critical for reliable prediction of drug binding and associated thermodynamics and kinetics [34]. Therefore, drugs with better docking scores, nadide, cangrelor and denufosol-Mpro complexes, were subjected to molecular dynamics (MD) for 100 ns, to investigate their stability, dynamics and conformational changes [5]. In this regard, we analysed root mean square deviation (RMSD), root mean square variation (RMSF) and a histogram that describes the interaction of ligand and amino acids and other properties of the ligand.

The stability of the three complexes was evaluated via RMSD. The RMSD gives insight into the structural conformation throughout the simulation. It measures the average change between Cα backbones from their initial conformation to their final position during the simulation [5]. The RMSD values for the three complexes are depicted in the plot (RMSD (nm) vs. Time (ns)) in Fig. 7 which revealed that Mpro has attained an equilibrium state around 2 Å after 5 ns of the simulation and persisted up to 40 ns. After that, minor fluctuations for the next 20 ns were observed; later, equilibrium was re-established around 2.4 Ả to the end of the simulation. This indicates the high stability of the Mpro backbone during the simulation. Whereas ligands RMSD have revealed that nadide was stabilized around 1.6 Ả in the first 58 ns of the simulation, small drift was then noticed around 60 ns, then equilibrium was restored around 2.4 Ả to the end of the simulation (Fig. 7A).

Fig. 7
figure 7

The protein–ligand root mean square deviation (RMSD) plot of short-listed compounds complexed SARS-CoV-2 main protease (Mpro) (PDBid: 5R81). (A) Nadide. (B) Denufosol. (C) cangrelor

Full size image

The ligand RMSD for denufosol-bound Mpro was fluctuating throughout the simulation, and stability was observed (2.2 Å) around the duration of 80–95 ns (Fig. 7B), while in Mpro bound cangrelor complex, the ligand RMSD was stabilized around 1.9 Å after 20 ns of the simulation, and minor fluctuation was displayed between 40 and 60 ns (Fig. 7C). Overall, a comparison of the three complex ligand RMSD revealed that nadide is superior to the other drugs followed by cangrelor which makes them interesting candidates for further investigations.

The root mean square fluctuations (RMSF) plot (Fig. 8) was used to analyse the residual fluctuation during the simulation. RMSF lower than 2.4 Å was observed with most of the residues except spike value displayed by residues 209–303, 166–170 and 191–192 which correspond to N- and C-terminal or loop regions that adopt a variety of conformations. The overall low fluctuation indicates a strong attachment of ligands with the protein.

Fig. 8
figure 8

The root mean square fluctuation (RMSF) of SARS-CoV-2 main protease (Mpro) (PDBid: 5R81). (A) Nadide. (B) Denufosol. (C) Cangrelor

Full size image

Protein–ligand interaction analysis revealed that the complex of Mpro and nadide was maintained by direct and indirect hydrogen bonds throughout the simulation. Thr24, His 41, Gly 143 and His 164 interacted with nadide via hydrogen bonds with 35%, 30%, 25% and 35% occupancy, respectively, while hydrogen bonds with Thr26 and Asn 142 persisted for 60% of the simulation time. Glu 166 displayed ultimate contact with nadide via hydrogen bond (100%) (Fig. 9A). Denufosol was mainly stabilized via direct and indirect hydrogen bonds with Thr 24 (80%), Cys 44 (45%) and Lys 61 (50%), while the ultimate contact was observed with Ser46 (Fig. 9B). In Mpro-cangrelor complex, hydrogen and bridged hydrogen bonds were observed with residues His 41(70%) and Gln 189 (25%), while Glu166 displayed ionic and bridged hydrogen bonds (25%) (Fig. 9C).

Fig. 9
figure 9

Protein–ligand contact histogram of SARS-CoV-2 main protease (Mpro) (PDBid: 5R81). (A) Nadide. (B) Denufosol. (C) Cangrelor

Full size image

As shown in Fig. 10, the following ligand properties, the radius of gyration (ROG), intramolecular H-bonds, molecular surface area (MolSA), solvent accessible surface area (SASA) and polar surface area (PSA), were further analysed.

Fig. 10
figure 10

Variation of ligands properties generated during molecular dynamic simulation (100 ns) of SARS-CoV-2 main protease (Mpro) (PDBid: 5R81). (A) Nadide. (B) Denufosol. (C) Cangrelor

Full size image

The radius of gyration (ROG) is used to calculate the spread of mass of ligands around their geometrical central axis [35], whereas firmer structures are associated with low ROG values, while high value indicates a more labile character of the complex [36]. Nadide, denufosol and cangrelor-Mpro complexes displayed a rGyr varied between 4.3–7, 4.5–7, 4.8–6.9 and 4.8–6.9Ả, respectively.

The formation of hydrogen bonds (Hb) in the active site of the protein plays a significant role in molecular recognition of the ligand [5, 29]. MD results showed that intramolecular Hbs were detected in the three complexes (Fig. 10).

With respect to the molecular surface (MolSA), which corresponds to the Van der Waals surface, nadide and cangrelor showed more constant fluctuations in MolsA values compared to denufosol. While solvent accessible surface area (SASA) shows the solvent-like behaviour of the protein–ligand complex, a lower value indicates the stability of the system [36]. Results obtained showed that the lowest SASA value was observed with nadide. Polar surface area (PSA) measures the solvent-accessible surface in a molecule that is only contributed by oxygen and nitrogen atoms. As shown in Fig. 10, the cangrelor-Mpro complex showed the least fluctuation.

Discussion

Since the declaration of the COVID-19 pandemic, several therapeutic options have been investigated to contain SARs-COV-2-related morbidity and mortality. However, limited specific treatment was considered for use [37]. There are multiple proteins involved in the virus life cycle considered to be potential targets to inhibit SARs-COV-2. Amongst them is SARS- CoV-2 main-protease (Mpro) [38,39,40,41,42,43,44,45,46,47,48]. Mpro is a primary enzyme responsible for the cleavage and activation of viral polyproteins (pp) into functional units that are pivotal in maintaining the viral life cycle [49]. Moreover, the active site of this protein is highly conserved, and there is no similarity between cysteine proteases in humans and that of SARS-CoV-2, thereby it is less possible to have candidates against SARS-CoV-2 proteases with cross-reactivity to those of humans [50]. All these features make Mpro a promising antiviral target.

RNA-dependent RNA polymerase (RdRp) is another relevant target that plays a critical role in the viral RNA transcription and synthesis step, hence, proliferation and growth of the virus [51]. The history of antiviral inhibitors that act against HIV and HCV suggests that multi-targeting approach is advantageous in terms of effectiveness and prevention of the emergence of resistance [52, 53]. In this regard, this study aims to repurpose compounds approved by international agencies as promising candidates to inhibit SARs-COV-2 by interfering with Mpro mainly and RdRp.

This study revealed that five approved drugs, natamycin, nadide, cangrelor, denufosol and fomidacillin, showed good binding affinity against RdRp and Mpro. In RdRp, the results revealed that the five short-listed compounds showed significant affinity at the binding site higher than that observed with the reference inhibitor. In addition, nadide, cangrelor, denufosol and fomidacillin are expected to inhibit RdRp in a similar way to remdesivir, by interacting with hydrophilic residues Lys545, Arg553 and Arg 555 in the NTP entry channel [37].

Top scoring drugs have also shown better binding affinity at the Mpro binding site compared to the non-covalent inhibitor Z1367324110. Non-covalent inhibition depends on hydrogen binding, and hydrophobic and ionic interactions to create a reversible interaction with the active site that may or may not involve the conserved catalytic residue Cys145 [54]. However, it is necessary to establish an interaction with one or more of the main interactive residues: Glu166, His41, Gly143, Ser144 and Cys145 [39]. As observed in Glide docking, the reference, natamycin and denufosol interacted with Glu 166 through hydrogen binding, while the final contact between nadide and Glu 166 was observed during molecular dynamic simulation. Thus, targeting this residue will block the catalytic activity of the enzyme [55].

Representative hits including nadide, cangrelor and denufosol showed good binding affinity in both target enzymes and formed interactions with their key residues. In particular, a nadide that is superior to other hits in the context of binding stability may serve as a drug candidate for COVID-19. Nadide is an adenine nicotinamide dinucleotide that exists widely in nature and participates in numerous enzymatic reactions, amongst which it acts as an electron carrier through alternating oxidation (NAD +) and reduction (NADH) [56]. Interestingly, it is suggested that nadide can inhibit SARs-COV-2 infection by targeting both the virus (Mpro and RdRp) and the host, since SARs-COV-2 infection causes a significant reduction in the host nicotinamide adenine dinucleotide (NAD), which is required to respond to innate immunity. Increased NAD via nadide may support innate immunity against SARs-COV-2 infection [56, 57].

Cangrelor is another promising potential drug candidate for COVID-19. Cangrelor is an analogue of ATP that directly blocks the P2Y12 receptor preventing platelet activation and aggregation [58]. P2Y12 inhibitors were used during the crisis of COVID-19 to prevent thrombotic events and to support the host’s defence against infection [59].

Conclusion

In this study, we applied different in silico approaches to drug repurposing for COVID-19. We propose three drugs, namely nadide, cangrelor and denufosol, which can target SARS-CoV-2 Mpro and RdRp. We selected these drugs based on their docking score, MMGBSA and MD simulation studies. In brief, these drugs demonstrated a very good binding affinity to the two targets and formed a stable interaction with Mpro, as evidenced by 100-ns MD. Further in vitro and in vivo studies will be required to confirm these results. These drugs have the potential to target SARS-CoV-2 Mpro and RdRp.