The SARS-CoV-2 helicase as a target for antiviral therapy: Identification of potential small molecule inhibitors by in silico modelling.

Although vaccines that provide protection against severe illness from coronavirus disease (COVID-19) have been made available, emerging variant strains of severe acute respiratory syndrome 2 coronavirus 2 (SARS-CoV-2) are of concern. A different research direction involves investigation of antiviral therapeutics. In addition to structural proteins, the SARS-CoV-2 non-structural proteins are of interest and this includes the helicase (nsp13). In this study, an initial screen of 300 ligands was performed to identify potential inhibitors of the SARS-CoV-2 nsp13 examining the nucleoside triphosphatase site (NTPase activity) as the target region. The antiviral activity of polyphenols has been previously reported in the literature and as a result, the phenolic compounds and fatty acids from the OliveNet™ library were utilised. Synthetic compounds with antimicrobial and anti-inflammatory properties were also selected. The structures of the SARS-CoV and MERS-CoV helicases, as well as the human RECQ-like DNA helicase, DHX9 helicase, PcrA helicase, hepatitis C NS3 helicase, and mouse Dna2 nuclease-helicase were used for comparison. As expected, sequence and structural homology between the various species was evident. A number of broad-spectrum and well-known inhibitors interacted with the NTPase active site highlighting the need to potentially identify more specific inhibitors for SARS-CoV-2. Acetylcysteine, clavulanic acid and homovanillic acid were identified as potential lead compounds for the SARS-CoV-2 helicase. Molecular dynamics simulations were performed with the leads bound to the SARS-CoV-2 helicase for 200 ns in triplicate, with favourable binding free energies to the NTPase site. Given their availability, further exploration of their potential inhibitory activity could be considered.

Introduction

Coronavirus disease 19 (COVID-19), which is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has continued to spread around the globe and a number of public health measures have been implemented to prevent transmission [1]. In addition to having an impact on the health of individuals, the socio-economic consequences of the COVID-19 pandemic have become apparent. Coronaviruses are enveloped viruses that consist of a positive-sense single-stranded RNA genome and to date, seven human coronaviruses have been recorded [2,3]. Like severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV), SARS-CoV-2 belongs to the beta-genera and all three viruses are highly pathogenic [2,4]. They predominantly infect cells of the respiratory tract and the coronavirus spike protein plays an important role in the viral entry stage [5]. In terms of SARS-CoV-2, studies have shown that the spike protein binds to the angiotensin-converting enzyme 2 (ACE2) receptor and the structural basis of this interaction has been elucidated [[5], [6], [7]]. The spike protein is also the main target of COVID-19 vaccines and neutralising monoclonal antibodies [8,9].

The first open reading frame (ORF1a/b) of the SARS-CoV-2 RNA genome encodes for 16 non-structural proteins, while the remaining ORFs encode for structural and accessory proteins [[10], [11], [12]]. The non-structural proteins (nsps) have been found to assemble into a replication and transcription (RTC) complex that is essential for viral survival [4,13]. Nsp16 and nsp14, for example, are S-adenosylmethionine (SAM)-dependent methyltransferases that are responsible for capping the viral mRNA and nsp10 is also an important cofactor [14]. The coronavirus nsp13 functions as a helicase, as it unwinds deoxyribonucleic acid (DNA) and ribonucleic acid (RNA) in a nucleoside triphosphate (NTP) dependent manner [15]. Interestingly, nsp13 is highly conserved across coronavirus species and the binding domains have been established [[15], [16], [17], [18]]. Genome sequencing has resulted in the identification of mutations in structural proteins, namely the spike protein, and non-structural proteins [19]. This includes nsp13 and two missense mutations, P504L and Y541C, were found to occur in the 2A domain [11]. In a study by Ugurel et al., the interactions between FDA-approved drugs and the wildtype and mutant SARS-CoV-2 helicases were examined using in silico methods. Although the mutations were predicted to result in structural changes, the compounds cangrelor, fludarabine, folic acid, and polydatin were found to interact with both the wildtype and mutant helicases [11]. Several studies have examined the association between the helicase and the RNA-dependent RNA polymerase (RdRp), as well as the auxiliary factors nsp7 and nsp8 (mini-RTC) [20,21]. Nsp13 is consequently an attractive target protein and prior to the crystal structure of the SARS-CoV-2 helicase being released, homology models were generated and were used to discover potential drug candidates [12,15].

Drug repurposing has formed a significant part of the response to the COVID-19 pandemic and there is a need to identify compounds that may have therapeutic potential or be used as prophylaxis [22,23]. This process provides an opportunity for existing therapeutics to be re-evaluated and for new applications to be identified [23]. Remdesivir was the first drug to be approved by the FDA as a treatment for COVID-19 and this is a potent RdRp inhibitor that can be used in patients requiring hospitalisation [24,25]. Furthermore, compounds that have anti-inflammatory properties and are able to alleviate symptoms have been of particular interest [26]. This includes phytochemicals and the protective properties of natural products against viruses have been investigated.

Computational methods have also enabled scientists to study various aspects of infectious diseases and virtual screening tools have made it possible for lead compounds to be discovered from large ligand databases [27]. As discussed by Wu et al., the main strategies that are being employed to combat the pandemic include testing broad-spectrum antivirals, screening existing databases for small molecules that may be effective against the virus and resulting disease, and developing novel drugs from scratch [28]. Here we investigated the SARS-CoV-2 helicase focusing on the NTPase active site. For comparison, the active site domain from various species including the microbial SARS-CoV, MERS-CoV, PcrA and hepatitis C NS3 helicases, and the mammalian RECQ-like DNA, DHX9 (human) helicases, and Dna2 nuclease-helicase (mouse). A library of 300 compounds consisting of a curated database of natural and synthetic compounds was utilised for initial screening. Following further modelling studies of interesting compounds, potential lead compounds including acetylcysteine, clavulanic and homovanillic acids, were identified for the SARS-CoV-2 helicase.

Materials and methods

Protein structures and ligands

The crystal structure of the SARS-CoV-2 helicase (PDB ID: 6ZSL) was obtained from the RCSB Protein Data Bank [29,30]. A single chain was isolated and the zinc ions were retained [31]. A second cryo-electron microscopy (cryo-EM) structure of the SARS-CoV-2 helicase in complex with the non-structural proteins nsp7, nsp8 and nsp12 (PDB ID: 6XEZ) was utilised and the two nsp13 chains were isolated (chains E and F) [21]. The zinc ions were retained, while the adenosine-5′-diphosphate (ADP) ligand was removed. Moreover, the SARS-CoV (PDB ID: 6JYT) and MERS-CoV (PDB ID: 5WWP) helicase crystal structures were examined [15,32]. A single chain was isolated from both structures and the zinc ions were retained. Likewise, the structures of the human RECQ-like DNA helicase (PDB ID: 2V1X), human RNA helicase DHX9 (PDB ID: 3LLM), mouse Dna2 nuclease-helicase (PDB ID: 5EAW), PcrA DNA helicase from Geobacillus stearothermophilus (PDB ID: 3PJR), and the hepatitis C virus NS3 helicase (PDB ID: 4OJQ) were obtained [[33], [34], [35], [36], [37]]. The monomer of each protein was generated by isolating a single chain and the zinc ions were kept for the human RECQ-like DNA helicase. ADP and adenosine-5′-triphosphate (ATP) ligands were removed.

The chemical structures of 300 ligands and ADP, which was used as the control compound, were obtained from the National Centre for Biotechnology Information (NCBI) PubChem Database [38]. The compounds were predominantly sourced from the OliveNet™ database and this included 215 phenolic compounds, as well 9 fatty acids [39]. OliveNet™ is a comprehensive library of 676 compounds from Olea Europaea and due to the bioactivities of phenolic compounds, it has been suggested that this subclass may provide a starting point for the development of novel treatments for viral infections [40,41]. The remaining ligands consisted of known drugs with antiviral, anti-inflammatory, antiparasitic, antioxidant, and antibiotic properties.

Docking to the nucleoside triphosphatase (NTPase) region

The Schrödinger Suite was used for molecular docking and the structures were imported into Maestro [42]. The Protein Preparation Wizard and LigPrep tool were used to prepare the proteins and ligands, respectively [[42], [43], [44], [45]]. The receptor grid was generated using the Receptor Grid Generation tool and the optimized potential for liquid simulations (OPLS3e) force field was used [46,47]. The receptor grid was 20 × 20 × 20 Å in size. For the SARS-CoV-2 (PDB ID: 6ZSL), SARS-CoV (PDB ID: 6JYT) and MERS-CoV (PDB ID: 5WWP) helicases, key residues of the NTPase region were used to define the binding site. They were K288, S289, D374, E375, Q404 and R567. The receptor grid was centred around ADP for the 6XEZ, 2V1X, 3LLM, and 5EAW nsp13 structures. The grid was generated around ATP for the PcrA DNA helicase (PDB ID: 3PJR) and residues T206, G207, S208, G209, K210, S211, and T212 for the hepatitis C virus NS3 helicase (PDB ID: 4OJQ) [37]. For the Hoechst 33342 ligand, the receptor grid dimensions were increased to 25 × 25 × 25 Å. The Glide Ligand Docking protocol was used for the initial screen of ADP and the 300 compounds. The Glide standard precision (SP) mode was selected for this stage and the results were refined further using the Quantum-Mechanics Polarized Ligand Docking (QPLD) protocol [48,49]. A selection of 30 compounds with a range of binding affinities was subsequently docked to the target binding site using the QPLD protocol and the extra precision (XP) mode was selected for improved docking accuracy [50,51]. The QM level was set to accurate and the GlideScore option was chosen for final selection. The protein-ligand interactions were observed using the Ligand Interaction Diagram tool.

Molecular dynamics simulations

Classical molecular dynamics (MD) simulations were performed with a selection of small molecules bound to the NTPase site of SARS-CoV-2 helicase using GROMACS 2018.4 software [52,53]. The structure of the SARS-CoV-2 helicase was obtained from the RCSB Protein Data Bank (ID: 6ZSL) [29,30]. Crystallographic water and solvent molecules were removed, and zinc ions were retained. Missing loops were modelled using Modeller 9.22 [54,55], generating five models. The structure with the lowest zDOPE score (−1.16) was selected [56]. Docked ligands served as starting structures with topologies generated using SwissParam [57]. Simulations were performed with the CHARMM36 force field [58] as previously described [59]. Protein systems were solved in a dodecahedral box containing 0.15 M NaCl with a minimum distance of 2.0 nm between protein atoms and the box edge. Simulations were performed for 200 ns in triplicate with a time-step of 2 fs, with random generation of velocities according to a Maxwell distribution.

Trajectories were visualised and analysed using Visual Molecular Dynamics 1.9.3 [60]. Root mean square deviation (RMSD), radius of gyration (Rg), and root mean square fluctuation (RMSF) calculations were performed using analysis tools within GROMACS. Binding free energy was calculated with molecular mechanics Poisson–Boltzmann surface area (MM-PBSA), using the g_mmpbsa tool [61] as previously described [59]. MM-PBSA calculations were performed on a 10 ns segment of the trajectory for each system preceding ligand detachment in any replicate: 135–145 ns for ADP-bound, 35–45 ns for acetylcysteine-bound, 20–30 ns for clavulanic acid-bound, and 60–70 ns for homovanillic acid-bound SARS-CoV-2 helicase (Fig. S2). Simulations and MM-PBSA calculations were carried out on the Spartan HPC GPU cluster [62].

PrankWeb, PyMOL alignment and blind docking

The PrankWeb server was used to predict potential binding sites in each structure using conservation analysis (Table S6) [63,64]. These regions were subsequently aligned in PyMOL and compared to the SARS-CoV-2 NTPase region in the 6ZSL structure [31]. The cryo-EM structure of the SARS-CoV-2 (PDB ID: 6XEZ chains E and F) was also aligned to the crystal structure (PDB ID: 6ZSL) and the RMSD values of the amino acids were recorded. The proteins were prepared as macromolecules and ligands using PyRx [65]. For blind docking, the receptor grid was generated around the entire protein. The exhaustiveness was set to 2048. The jobs were run using AutoDock Vina on the cloud computing server Galileo (Hypernet labs) [66,67].

Results

SARS-CoV-2 nsp13

A library of 300 compounds consisting of natural and synthetic ligands (Table S1) was screened against the NTPase site of the SARS-CoV-2 helicase (PDB ID: 6ZSL). Based on the PrankWeb binding site analysis, the active NTPase region was comprised of residues P284, G285, G287, K288, S289, H290, A312, A313, A316, K320, D374, E375, S377, G400, Q404, R442, R443, Q537, G538, E540 and R567 (Fig. 1 ). The GlideScores that were generated from molecular docking ranged from −2.3 to −9.3 kcal/mol. ADP was used as the positive control in this study and the GlideScore was found to be −6.5 kcal/mol. Based on the binding affinities, commercial availability, and safety profile, 30 compounds were selected for further analysis. This included the OliveNet™ phenolic compounds hellicoside, rutin, homovanillic acid, cyanidin-3-O-glucoside and oleocanthal. A number of antibiotics (amikacin, meropenem, cefuroxime, ciprofloxacin, levofloxacin, ceftazidime, moxifloxacin, cefotaxime, clavulanic acid, amoxicillin and ceftriaxone), protease inhibitors (nelfinavir, indinavir, saquinavir, darunavir, ritonavir, and lopinavir), DNA-binding ligands (Hoechst 33342 and 33258), chemotherapy drugs (doxorubicin), anti-inflammatory and antioxidant compounds (indomethacin, (−)-epigallocatechin gallate and acetylcysteine), sirtuin activators (SRT2104) and RNA-dependent RNA polymerase inhibitors (remdesivir) were also chosen for comparison.

Fig. 1

Structures of the SARS-CoV-2 helicase. The domains of the SARS-CoV-2 nsp13 are labelled (PDB ID: 6ZSL) and the complex that the helicase forms with nsp12, nsp7 and nsp8 is also shown (PDB ID: 6XEZ). The NTPase pocket that was identified through the PrankWeb server can be seen (coloured tan). The protein-ligand interactions of ADP with the NTPase region of the SARS-CoV-2 nsp13 structures (PDB ID: 6ZSL and PDB ID: 6XEZ chain E) are also provided.

Using the QPLD protocol of the Schrödinger Suite for greater docking accuracy, ADP and a selection of 30 ligands were docked to the NTPase site (Table S2). Amikacin had the strongest GlideScore (−12.4 kcal/mol), while ritonavir had a weaker binding affinity (−3.2 kcal/mol). The GlideScore for ADP was −9.5 kcal/mol (Fig. 1). Blind docking was also performed on the SARS-CoV-2 helicase to investigate whether the ligands would preferentially bind to the NTPase region or another site in the protein, which may be a potential allosteric binding site. The compounds that were predicted to have conformations positioned in the NTPase region were amikacin, meropenem, cefuroxime, indomethacin, acetylcysteine, ciprofloxacin, levofloxacin, ceftazidime, (−)-epigallocatechin gallate, cefotaxime, clavulanic acid, doxorubicin, amoxicillin, Hoechst 33258, Hoechst 33342, cyanidin-3-O-glucoside, homovanillic acid and oleocanthal. The ligands acetylcysteine (6 out of 9 poses), cefotaxime (3 out of 9 poses), clavulanic acid (7 out of 9 poses), doxorubicin (2 out of 9 poses), amoxicillin (3 out of 9 poses), Hoechst 33342 (3 out of 9 poses), Hoechst 33258 (2 out of 9 poses) and homovanillic acid (7 out of 9 poses) were selected for the subsequent stage (Table S3). The number of conformations to be generated for these compounds from blind docking was increased and the results can be seen in Fig. 2 .

Fig. 2

Molecular docking and blind docking results for the SARS-CoV-2 crystal structure (PDB ID: 6ZSL). The compounds were docked to the NTPase site of the SARS-CoV-2 helicase using the QPLD protocol and the GlideScores (kcal/mol) are provided. The protein-ligand interactions are depicted for acetylcysteine, amoxicillin, cefotaxime, clavulanic acid, doxorubicin, Hoechst 33258, Hoechst 33342, and homovanillic acid. The polar residues are coloured dark blue, the positively charged residues are coloured purple, the negatively charged residues are coloured maroon, the hydrophobic residues are coloured green, and the glycine residues are coloured dark yellow. The blind docking results are also provided for these compounds and the number of poses that were found to be in the NTPase region can be seen. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

In terms of protein-ligand interactions, ADP formed hydrogen bonds with residues E540, Q404, R567, and K288 (Fig. 1). Salt bridges were also present with residues R443, R567, and K288. Acetylcysteine had a stronger GlideScore (−9.8 kcal/mol) than ADP, and this was followed by homovanillic acid, cefotaxime, clavulanic acid, doxorubicin, amoxicillin, Hoechst 33258, and Hoechst 33342. The interactions that each of the ligands made with the protein residues are described in Fig. 2. It was evident that six ligands formed intermolecular bonds with K288 including amoxicillin, acetylcysteine, clavulanic acid, cefotaxime, Hoechst 33258, and homovanillic acid. With the exception of amoxicillin and homovanillic acid, the ligands formed hydrogen bonds with G538. In addition to K288, several compounds formed intermolecular bonds with S289 and R567. This included amoxicillin, acetylcysteine, and clavulanic acid. Based on the docking results for the SARS-CoV-2 helicase NTPase region, ADP and the ligands homovanillic acid, clavulanic acid, and acetylcysteine were selected for further analysis.

The PrankWeb analysis revealed that pocket 2 and pocket 1 were the NTPase sites on chains E and F of the 6XEZ SARS-CoV-2 cryo-EM structure, respectively (Table S6). Molecular docking was conducted on the NTPase site, and ADP had a GlideScore of −7.6 kcal/mol for chain E and −6.6 kcal/mol for chain F (Fig. 1 and Fig. S1). On chain E, ADP formed interatomic contacts with D374 (H-bond), T286 (H-bond), E540 (H-bonds), K320 (salt bridges), and K569 (H-bond). On chain F, ADP formed intermolecular bonds with G282 (H-bond), Q404 (H-bond), K320 (H-bonds and salt bridge), R567 (H-bond), and E540 (H-bonds). ADP was the strongest binding ligand for chain E and this was followed by homovanillic acid, clavulanic acid, and acetylcysteine (Fig. S1). The GlideScores for chain E ranged from −3.8 to −7.6 kcal/mol. For chain F, ADP was the strongest binding ligand (−6.6 kcal/mol) and this was followed by acetylcysteine, homovanillic acid, and clavulanic acid. The GlideScores for chain F ranged from −3.7 to −6.6 kcal/mol (Fig. S1). The blind docking results revealed that acetylcysteine had 11 out of 20 poses, clavulanic acid had 10 out of 20 poses, homovanillic acid had 2 out of 18 poses, and ADP had 16 out of 20 poses in the NTPase region of chain E. For chain F, ADP had 15 out of 19 poses, and clavulanic acid had 4 out of 20 poses within the NTPase site (Fig. S1). The cryo-EM structure of the SARS-CoV-2 helicase (PDB ID: 6XEZ) was aligned to the crystal structure (PDB ID: 6ZSL) and when examining several of the residues that were within 5 Å of the docked ligands, some of the amino acids had greater RMSD values in the 6XEZ cryo-EM structure (Table S4).

Stability of ligands bound to SARS-CoV-2 nsp13 NTPase site

The stability of potential small molecule SARS-CoV-2 helicase inhibitors in a dynamic aqueous system was examined using classical MD simulations. Overall, five separate systems were studied and simulated for 200 ns in triplicate: the apo SARS-COV-2 helicase, and the protein bound with ADP, acetylcysteine, clavulanic acid, and homovanillic acid. Out of the ligands studied, the natural substrate ADP demonstrated the greatest stability in binding to the NTPase site of SARS-CoV-2 helicase (Fig. S2). Two out of three replicates remained bound to the protein throughout the entire trajectory, with unbinding occurring in one replicate after approximately 150 ns. Similarly, while acetylcysteine also remained bound to the protein for two out of three replicates, unbinding occurred in the remaining replicate much earlier at approximately 50 ns. For both clavulanic and homovanillic acids, one out of three replicates were bound to the protein for 200 ns. Unbinding occurred for clavulanic acid at 30 and 120 ns, and for homovanillic acid at 70 and 80 ns into the trajectory.

Frames were extracted from a single trajectory at 100, 150, and 200 ns to examine protein-ligand hydrogen bonds. When examining the intermolecular bonds that were predicted to occur, several key active site residues formed hydrogen bonds with the ligands (Table 1 ). They included K288, S289, and R567. Residues K320, R443, and E540 were also predicted to form hydrogen bonds (Table 1).

Table 1

The SARS-CoV-2 helicase residues that were predicted to form hydrogen bonds with the ligands based on frames that were extracted from a single trajectory at 100, 150, and 200 ns are provided.

Ligand	Trajectory frame	Hydrogen bonds
ADP	100 ns	G285, S289, K320, R443, E540
	150 ns	K320, R442, R443, E540, R567
	200 ns	G285, S289, R442, R443, K320, E540, R567
Clavulanic acid	100 ns	G285, T286, G287, K288, H290
	150 ns	G285, T286, K288, S289, H290
	200 ns	G285, T286, G287, K288, S289, H290, R443
Homovanillic acid	100 ns	K320
	150 ns	K320
	200 ns	K320
Acetylcysteine	100 ns	Q404, R443, G538, R567
	150 ns	Q404, R443, G538, R567
	200 ns	None

The overall protein structure was not greatly affected by small molecule binding to the NTPase site of helicase. Average root mean square deviation (RMSD) of the protein backbone was slightly higher for apo (0.40 nm) and homovanillic-acid bound (0.43 nm) helicase (Fig. 3 A). ADP-, acetylcysteine-, and homovanillic acid-bound helicase had similar average RMSD values of 0.34, 0.32, and 0.34 nm, respectively. Radius of gyration (Rg) was largely similar, with an average Rg of 2.8 nm for all systems throughout the trajectory (Fig. 3B). Protein flexibility was examined by calculating the root mean square fluctuation (RMSF) of both the protein backbone and sidechains (Fig. 3C and D, respectively). The zinc-binding domain (ZBD) of the helicase was shown to be the most flexible, with large fluctuations in both the protein backbone and sidechains. The 1B domain was also shown to be flexible. Peaks in RMSF were observed in regions around residues 340 and 480 of the RecA1 and RecA2 domains, corresponding to flexible loops on the outer surface on the protein.

Fig. 3

Classical MD simulations of SARS-CoV-2 nsp13 in its apo (grey), and in complex with ADP (purple), acetylcysteine (blue), clavulanic acid (green), and homovanillic acid (yellow) were performed for 200 ns in triplicate. A) Average root mean square deviation (RMSD) and B) average radius of gyration (Rg) of the protein backbone. Average root mean square fluctuation (RMSF) of C) protein backbone and D) sidechains. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Binding free energy of small molecules to SARS-CoV-2 nsp13 NTPase site

To quantitate binding free energy of the small molecules to the NTPase active site of SARS-CoV-2 helicase, trajectories were further analysed using molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) calculations. ADP was the strongest binding ligand with an average ΔGbinding of −109.0 kcal/mol (Table 2 ). Acetylcysteine, clavulanic acid and homovanillic acid were found to have similar average ΔGbinding, with values of −27.9, −27.2, and −27.2 kcal/mol respectively (Table 2). Electrostatic forces were the greatest contributor to favourable ligand binding, along with smaller contributions from van der Waal and non-polar interactions.

Table 2

Average binding energy of small molecules bound to the NTPase site of SARS-COV-2 helicase. Energies are shown as average ± standard error in kcal/mol.

Ligand	ΔEelec	ΔEvdW	ΔGpolar	ΔGnonpolar	ΔGbinding
ADP	−394.9 ± 48.2	−13.5 ± 2.6	302.8 ± 43.2	−3.3 ± 0.3	−109.0 ± 9.6
Acetylcysteine	−111.3 ± 5.7	−11.1 ± 0.8	96.8 ± 5.9	−2.3 ± 0.1	−27.9 ± 2.0
Clavulanic acid	−85.4 ± 15.1	−13.7 ± 2.7	74.2 ± 20.3	−2.3 ± 0.2	−27.2 ± 4.0
Homovanillic acid	−92.2 ± 3.4	−7.1 ± 2.3	74.0 ± 2.2	−2.0 ± 0.1	−27.2 ± 1.3

Binding energy was decomposed into per-residue basis. While residues across the entire protein contributed both favourably and unfavourably to ligand binding (Fig. S3), the residues that were most prominent in terms of energy contribution were located within the NTPase active site highlighted in Fig. 4 . Conserved residues involved in NTPase hydrolysis were among these, with K288 and R567 contributing very favourably to the binding of all ligands studied. Energy contributions from K288 were −25.9 kcal/mol for ADP, −8.9 kcal/mol for acetylcysteine, −7.7 for clavulanic acid, and −7.3 kcal/mol for homovanillic acid. Similarly, respective energy contributions for R567 were −25.7, −6.7, −5.9, and −6.4 kcal/mol.

Fig. 4

Energy contribution of residues to binding of ligands to the NTPase site of SARS-CoV-2 nsp13. Residues are coloured according to their energy contribution in kcal/mol, with favourable contributions in red and unfavourable energy contributions in blue. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

While D374 and E375 are also conserved NTPase residues of helicase, these residues were found to contribute unfavourably to ligand binding: +23.3 and + 23.5 kcal/mol for ADP, +7.3 and + 10.7 kcal/mol for acetylcysteine, +6.5 and + 7.2 kcal/mol for clavulanic acid, and +6.1 and + 7.4 kcal/mol for homovanillic acid. E540 was another residue within the NTPase site contributing unfavourably to ligand binding with energy contributions of +30.3 kcal/mol for ADP, +7.5 kcal/mol for acetylcysteine, +9.2 kcal/mol for clavulanic acid, and +7.6 kcal/mol for homovanillic acid. Favourable energy contributions were also apparent from K320 and R443 located in the NTPase site: 22.0 and −28.1 kcal/mol for ADP, −7.3 and −8.9 kcal/mol for acetylcysteine, −6.0 and −6.5 kcal/mol for clavulanic acid, and −6.8 and −4.7 kcal/mol for homovanillic acid.

SARS-CoV and MERS-CoV nsp13

Doxorubicin, cefotaxime, homovanillic acid, Hoechst 33258, clavulanic acid, amoxicillin, acetylcysteine, and Hoechst 33342 were docked to the SARS-CoV (PDB ID: 6JYT) and MERS-CoV (PDB ID: 5WWP) helicase crystal structures. Based on the PrankWeb results, pocket 3 for the SARS-CoV helicase and pocket 1 for the MERS-CoV helicase contained the residues from the conserved NTPase region (Fig. S8 and Table S6). The GlideScores for the SARS-CoV-2 nsp13 ranged from −3.3 to −5.6 kcal/mol, while the GlideScores for the MERS-CoV nsp13 ranged from −3.7 to −10.4 kcal/mol.

For the SARS-CoV structure, the strongest binding ligand was doxorubicin and this was followed by cefotaxime, homovanillic acid, Hoechst 33258, clavulanic acid, amoxicillin, acetylcysteine, and Hoechst 33342 (Fig. 5 ). In terms of the interactions with the protein residues, several ligands formed intermolecular bonds with K288. This included amoxicillin (H-bond), acetylcysteine (H-bond and salt bridge), clavulanic acid (H-bond and salt bridge), cefotaxime (H-bond and salt bridge), Hoechst 33258 (π−π cation) and homovanillic acid (π−π cation and salt bridge). Amoxicillin and Hoechst 33342 formed hydrogen bonds with E375, Hoechst 33258 formed hydrogen bonds with S289 and D374, and homovanillic acid formed a salt bridge with R567. The residues R443, K320, S310, Q537, D534, E201 and A316 were also involved in interatomic atomics with the ligands (Fig. 5).

Fig. 5

Interactions of the ligands with the SARS-CoV and MERS-CoV helicase crystal structures. The intermolecular bonds that formed between the compounds and the protein residues in the NTPase site are depicted. The GlideScores (kcal/mol) of the ligands acetylcysteine, amoxicillin, cefotaxime, clavulanic acid, doxorubicin, Hoechst 33258, Hoechst 33342, and homovanillic acid are provided. The polar residues are coloured dark blue, the positively charged residues are coloured purple, the negatively charged residues are coloured maroon, the hydrophobic residues are coloured green, and the glycine residues are coloured dark yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Homovanillic acid had the strongest GlideScore for the MERS-CoV helicase, and this was followed by clavulanic acid, doxorubicin, acetylcysteine, amoxicillin, cefotaxime, Hoechst 33258, and Hoechst 33342 (Fig. 5). A number of ligands formed intermolecular bonds with the residues E375 (amoxicillin, acetylcysteine, clavulanic acid, doxorubicin, and Hoechst 33258: H-bonds), K288 (clavulanic acid, cefotaxime, homovanillic acid, amoxicillin, and cefotaxime: H-bonds and π−π cations), R567 (amoxicillin, cefotaxime, acetylcysteine, clavulanic acid, and homovanillic acid: H-bonds and salt bridges), Q404 (amoxicillin and acetylcysteine: H-bonds) and S289 (doxorubicin: H-bonds). The blind docking results showed that Hoechst 33258 and Hoechst 33342 had poses within the NTPase region of the SARS-CoV helicase, whereas cefotaxime, doxorubicin, Hoechst 33258, and Hoechst 33342 had poses within this site for the MERS-CoV helicase (Fig. S4).

Helicase from other species

The PrankWeb analysis for each crystal structure revealed that pocket 1 of the human RECQ-like DNA helicase, human RNA helicase, mouse Dna2-nuclease helicase, and PcrA DNA helicase were conserved (Fig. S8 and Table S6). These pockets were the binding sites of ADP in the human RECQ-like DNA helicase, human RNA helicase, and mouse Dna2-nuclease helicase, and ATP in the PcrA DNA helicase. Pocket 3 for the hepatitis C virus NS3 helicase consisted of residues T206, E291, H293, A323, T324, V329, V331, V456 and T459, and this was adjacent to the ATP phosphate binding site. The binding site residues for each helicase structure were aligned to the NTPase region in the SARS-CoV-2 nsp13 crystal structure (PDB ID: 6ZSL). The alignment results can be seen in Fig. S8 and the RMSD values can be found in the Supplementary Information (Table S5).

The ligands were docked to each crystal structure using the Schrödinger Suite and blind docking was also performed. For the human RECQ-like DNA helicase and DHX9 helicase, clavulanic acid was predicted to be the strongest binding ligand and this was followed by acetylcysteine and homovanillic acid (Fig. S5). Furthermore, the crystal structures of the mouse Dna2-nuclease helicase, the PcrA DNA helicase and hepatitis C virus NS3 helicase were examined. The binding affinities and intermolecular bonds that the compounds formed with the protein residues can be seen in Figs. S6 and S7.

Discussion

Helicases can be divided into six superfamilies based on several characteristics including the presence of conserved motifs in their sequences, the type of nucleic acid they can bind and unwind, and the polarity of unwinding [[68], [69], [70]]. The SARS-CoV-2 nsp13 has been found to unwind RNA in the 5′ to 3′ direction and this protein belongs to the superfamily 1 (SF1) helicases [71]. Shu et al. demonstrated that the binding and hydrolysis of NTPs plays a crucial role in this process, as it provides a source of energy [71]. Although the SARS-CoV-2 helicase was able to hydrolyse four kinds of NTPs, there was a preference for ATP and guanosine-5′-triphosphate (GTP) [71]. The NTPase and unwinding activities of the protein were also affected by the presence of divalent metallic ions and Mg2+ was most efficient at supporting these activities [71]. This is also the case for the SARS-CoV and MERS-CoV helicases [32,70,72,73]. The RNA 5′-triphosphatase activity of coronavirus helicases has also been explored [73].

The SARS-CoV-2, SARS-CoV and MERS-CoV helicases have a triangular pyramidal shape, and five domains have been identified [21,32]. There are two Rec-A like domains (1A and 2A domains), a 1B domain, a stalk domain, and an N-terminal zinc binding domain (ZBD) [21]. There are six key residues that are involved in NTP hydrolysis and they are K288, S289, D374, E375, Q404, and R567 [15,74]. The NTPase region can be targeted by small molecules and some of the compounds that have been investigated for the SARS-CoV-2 nsp13 include bismuth salts (bismuth potassium citrate and ranitidine bismuth citrate), simeprevir, scutellarin, ertapenem, cangrelor, cepharanthine and ergoloid [11,12,16,17,74]. Antiviral compounds such as bananin and its derivatives have also been identified as potential inhibitors, as they have been previously found to be effective against the SARS-CoV helicase [75,76]. As aforementioned, in silico tools allowed for the structure of the SARS-CoV-2 nsp13 to be predicted before the crystal structure was made available and this information has formed the foundation for the drug discovery process [12].

In the current study, the PrankWeb ligand binding site prediction server recognised the NTPase site as a ligand binding pocket and this was the region of interest in the SARS-CoV-2, SARS-CoV and MERS-CoV helicases [63]. The 6ZSL structure from the PDB was used for the initial screen of the 300 ligands against the SARS-CoV-2 nsp13. By analysis of docking scores, bioactivity profiles, and commercial availability 30 compounds including antibiotics, antiviral agents, immunomodulators, polyphenols with antioxidant and anti-inflammatory properties, and antineoplastic drugs, were selected.

Inflammation and tissue damage have been observed in the lungs of patients with severe COVID-19, and there is also the risk of complications and secondary infections [77]. Although antibiotic resistance is a major concern, research is being conducted into the possible therapeutic effects of various antimicrobials and whether they are able to interfere with the virus lifecycle [77,78]. In saying this, the mechanisms of action of synthetic and natural compounds that may have antiviral properties and may assist with regulating the immune system are still being explored [79,80].

ADP was the ligand present in the cryo-EM structure of nsp13 (PDB ID: 6XEZ) and was consequently used as the positive control. ADP had a strong GlideScore for both structures of the SARS-CoV-2 helicase and interacted with several protein residues. When taking into consideration the protein-ligand interactions, binding affinities and number of poses that were predicted to be in the NTPase region from blind docking, the ligands acetylcysteine, clavulanic acid and homovanillic acid were identified as potential lead compounds. The data obtained from the MD simulations also demonstrated that ADP was binding strongly to the NTPase site of the SARS-CoV-2 helicase and that the lead compounds were binding favourably to this region. Amoxicillin with clavulanic acid is a widely used antimicrobial formulation used to treat a number of bacterial infections, including pneumonia [81]. Acetylcysteine acts as a potent antioxidant and mucolytic agent, and has been proposed as a potential treatment of SARS-CoV-2 [82,83]. Homovanillic acid is a phenolic compound that belongs in the hydroxyphenylacetic acid subclass, as characterised in the OliveNet™ database [39].

Additionally, DNA and RNA helicases from other species were examined for comparison. When developing drugs, specificity is a challenge and it is important to differentiate between viral and mammalian cellular helicases [84]. The lead compounds were docked to the NTP binding site in each of the structures and the results were compared to that of SARS-CoV-2. This included the human RECQ1 helicase and the conserved domain 1 of the DEIH-motif-containing helicase DHX9 [34,85]. Both of these proteins belong to the superfamily 2 helicases and the unwinding of nucleic acids is driven by the hydrolysis of NTPs [34,85]. RecQ helicases are critical for DNA replication, recombination and repair, while DHX9 participates in numerous cellular processes involving RNA [86]. Interestingly, the involvement of DHX9 in viral biology and its potential as a therapeutic target are being explored [87,88]. It is well known that the Hoechst ligands and doxorubicin are able to interact with DNA and can inhibit the catalytic activity of DNA helicase [89,90]. These compounds were also docked to the SARS-CoV-2 helicase for comparison. The results revealed that these ligands had weaker GlideScores and fewer blind docking poses positioned within the NTPase region.

The substrate complex of the PcrA DNA helicase, the NS3 helicase from the hepatitis C virus, and the mouse Dna2 nuclease/helicase were also analysed. Although the sequences differ between the species, there may be structural similarities and this could also be seen in the alignment results. Kwong et al. also highlighted how there are conserved structural elements in helicase structures and they showed an overlay of the HCV and PcrA helicases [84]. In terms of the ligand binding site analysis server results, the NTPase site was identified as a ligand binding pocket for all structures except PDB: 4OJQ. For the 4OJQ protein, the third ranked pocket was the closest to the ATP phosphate binding site. The first ranked pocket was the region that the co-crystallised inhibitor, (5-bromo-1H-indol-3-yl)acetic acid, was bound to and this site may be of interest in future studies.

In general, given the need of functional helicases for viral replication investigation of specific inhibitors is important. Despite difficulties associated with the overlapping activities between viral and cellular enzymes, there is an interest in developing targeted small molecules. In this context, compounds that target the conserved NTPase region of helicase among coronaviruses may be developed as antiviral drugs that act in a specific manner. Regarding the SARS-CoV-2 helicase, the small molecules acetylcysteine, homovanillic acid, and clavulanic acid were identified, in our studies, as potential lead compounds for further examination to investigate specificity and in antiviral assays. Given their availability, known bioactivity, and safety profiles these compounds are promising candidates for future studies.

Author contributions statement

TCK and AH conceptualized the aims and methodology and were involved in supervision. EP performed data analysis, data curation, and was involved in production of the first draft of the manuscript. JL was involved in data analysis and curation and was involved in production of the first draft of the manuscript. All authors contributed to editing and reviewing the manuscript.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Epigenomic Medicine Program (TCK) is supported financially by McCord Research (Iowa, USA), which has a financial interest in dietary compounds described in this work. However, there is no conflict of interest with respect to the inhibition of the SARS-CoV-2 helicase. The remaining co-authors also have no conflicts of interest.