In Silico Exploration of Bisphosphonate Scaffolds as Potential Inhibitors of SARS-CoV-2 RdRp for COVID-19 and PASC

1. Introduction

Coronaviruses are responsible for causing seasonal respiratory tract infections (RTIs) in people and are associated with common cold symptoms [1]. The highly pathogenic human coronaviruses (HCoVs) such as severe acute respiratory syndrome (SARS) associated coronavirus (SARS-CoV), Middle East respiratory syndrome–related coronavirus (MERS-CoV) and the novel SARS-CoV-2 cause infection to epithelial cells of the bronchi and pneumocytes, which could lead to life-threatening lung injuries [2]. Among these three, the SARS-CoV-2 virus which emerged in December 2019 exhibits faster human-to-human transmission, and resulted in over 780 million confirmed cases and 7 million reported deaths due to the novel coronavirus disease (COVID-19) [3]. Long COVID, also known as postacute sequelae of SARS-CoV-2 infection or COVID-19 (PASC) [4,5], may affect typically 3.1% of the population, especially those aged 35 to 69 years, females, people living in more deprived areas, those working in social care, those aged 16 years or over who were not working and not looking for work, and those with another activity-limiting health condition or disability [6].

SARS-CoV-2 employs a multi-subunit machinery for replication and transcription. Non-structural proteins (Nsp's) produced as cleavage products due to the open reading frame 1a and 1b (ORF1a and ORF1b) facilitate viral replication and transcription [7]. One of these, known as Nsp12 or RNA-dependent RNA polymerase (RdRp), catalyzes the synthesis of viral RNA and plays a central role in the replication and transcription cycle of SARS-CoV-2 with Nsp7 and Nsp8 as co-factors [8,9]. Therefore, Nsp12 is considered a primary target for antiviral agents, with the potential for treating COVID-19 [10], and possibly other coronaviral diseases because it is a highly conserved motif. For example, sequence alignment results from the literature shows 96% common identity between SARS-CoV and SARS-CoV-2 [11].

In this work, as we shall be using RdRp as a key target let us first describe its components: it has the Nsp12 catalytic subunit, two accessory subunits (Nsp8 and Nsp7), and more than two turns of RNA template-product duplex [12]. The RdRp domain is analogous to a cupped right hand, consisting of the finger ‘F’ (amino acid residues 398–581, 628–687), palm (amino acid residues 582–627, 688–815) and thumb (amino acid residues 816–919) subdomains found in all single-subunit polymerases [13]. An in-depth structural analysis depicts the Nsp12 subunit binding to the first turn of RNA between its F and thumb sub-domains. The core protein consists of a single chain of 942 amino acids. The active site comprises of five conserved Nsp12 elements that are found in the palm motif. The amino acids Asp760 and Asp761 (Figure 1) are necessary for synthesis, which binds to the 3′ end of the RNA. The RNA template is positioned by the supplementary Nsp12 finger motif. The second turn is positioned by two copies of Nsp8 that bind to the cleft on the opposite sides. As RNA exits, large helical extensions of Nsp8 protrude and create positively charged sliding poles that are necessary for coronaviruses to replicate their lengthy genomes [14]. Structural stability requires two Zinc (Zn) ions interacting with the residues present in the N-terminal domain (His295, Cys301, Cys306, Cys310) and finger domain (Cys487, His642, Cys645, Cys646). The presence of Zn in this site indicates its crucial role in stabilizing the overall 3D structure of the protein [15]. The binding of drugs to the amino acid residues in motif F of RdRp results in the aversion of entry of the substrate and divalent cations into the central active site cavity, thereby inhibiting the catalytic activity of the enzyme and preventing the RNA replication [16].

So far, the standard drugs of choice for treatment have been the emergency use authorized (EUA) drugs, nirmatrelvir and remdesivir. Nirmatrelvir is an orally available main protease inhibitor [17], while remdesivir is an RdRp inhibitor that is administered parenterally. The EUA status of remdesivir was revoked on April 2022, and approval as a supplemental new drug application (sNDA) has been given in December 2022 [16]. However, a recent study published in The Lancet shows that while remdesivir could reduce the risk it had an insignificant effect in ventilated COVID-19 patients [17]. Nirmatrelvir, a main protease inhibitor, exhibited promising antiviral effects but was susceptible to rapid degradation. To mitigate this issue, ritonavir, a protease inhibitor, was incorporated in a combined formulation which received FDA approval as the first oral antiviral pill PAXLOVID [18]. Other drugs that have shown promise include molnupiravir [19], favipiravir [20] and fluvoxamine [21,22]. However, molnupiravir and favipiravir are no longer recommended as the former drug has poor clinical outcomes [23,24,25] and the latter showed ineffective in viral clearance [26,27,28]. Although fluvoxamine has been shown to have both immunomodulatory effects [29] as well as anti-viral effect [21,30,31], it is clear that this drug can only be used in combination and not in its own [21,29].

Therefore, the quest to identify potential molecules to treat COVID-19 and PASC is still ongoing [32]. In the literature, when databases of molecules including synthetic and natural origin were screened against RdRp [33,34,35,36,37,38,39,40] together with other non-structural proteins (Nsp's) [41,42,43,44] inhibitors, such studies suggest the use of natural substances as alternative to remdesivir. Few studies also reported the screening of analogues of different scaffolds such as quinolines [40], cytidines [45] and andrographolides [46]. Notably, drug repurposing efforts across the world have proposed several drugs that are yet to be tested in vitro [13,32,38,39,40,47,48,49,50,51,52,53,54]. As RdRp is a key target, this paper looks at a class of small molecules called bisphosphonates (BPs) for four reasons: 1. The use of BPs is associated with a significant three- to five-fold reduction in the incidence of SARS-CoV-2 testing, COVID-19 diagnosis, and COVID-19-related hospitalization during the pandemic [55]; 2. BP scaffolds exhibit competitive binding to the nidovirus RdRp-associated nucleotidyl (NiRAN) transferase domain [56]; 3. Our prior experience with an immunomodulatory drug, fluvoxamine in in silico and ex vivo studies [21]; and 4. Our recent study which has shown that alendronate, a BP drug, shows more promising in silico results compared to remdesivir [32]. Therefore, we shall further explore BPs in detail as it is the focus of this follow-up study.

BPs are a class of small-molecule drugs that have two phosphonate groups. They are categorized into nitrogen-containing (amino-BPs) and nitrogen-free BPs (non-amino-BPs). They are mainly used to treat osteoporosis, Paget's disease of the bone, and to lower high calcium levels in people with cancer [57]. In addition, amino-BPs control the activation, expansion, and function of a significant portion of human γδT cells (i.e., it reduces the amount of circulation γδT cells), as well as neutrophils, monocytes and macrophages. They can also modify the dendritic cell's ability to present antigens to the immune system. Results from animal experiments indicate that both amino-BPs and non-amino-BPs have strong adjuvant-like effects of increasing antibody and T-cell responses to viral antigens. [58]. Acknowledging the range of immunomodulatory effects and their binding with RdRP inspired us to investigate these as prospective COVID-19 drug candidates [55,56,58].

2. Materials and Methods

In this study, a series of compounds having BPs were downloaded from the ChEMBL database [59] and molecular docking study was carried out to get all the “hit” molecules. These hits were further down-selected using molecular mechanics with generalized Born and surface area (MM-GBSA) studies and the top four ligands were analyzed for their dynamic behavior via molecular dynamics studies followed by the entropic calculations using the software “gromacs” (v2022.4).

2.1. Protein Preparation

While the literature shows a high sequence identity (96%) between the RdRp proteins of SARS-CoV and SARS-CoV-2, we couldn’t find any report with MERS-CoV, so the latter’s RdRp domain (residues 4845–5412 from accession code A0A023YA54) was taken and multiple sequence alignment was carried out to reveal the conserved motifs (75%) for these three coronaviruses (Figure 1). Further, to provide a valuable starting point for BPs screening against SARS-CoV-2, RdRp protein (PDB: 6M71) [8] was selected and downloaded from the protein data bank (https://www.rcsb.org/). The protein was prepared using the protein preparation wizard of the glide module of the Schrödinger software (Schrödinger LLC., NY, v2022) [60]. Missing hydrogens and residues were added using the software’s prime module and pre-processed using the Epik module at pH 7.0±2.0 [61]. The protein was optimized by removing water molecules beyond 3 Å and subjected to minimization using the optimized potentials for liquid simulations 4 (OPLS4) force field [61].

2.2. Ligand preparation

ChEMBL database has been used to download molecules containing BPs based on the Tanimoto similarity check (95%) [44]. These molecules were imported and prepared using Ligprep module of maestro (Schrödinger LLC., NY, v2022)). The ionization state was set to neutral and chirality was determined from the 3D structure. The force field OPLS4 was employed to prepare the ligands.

2.3. Molecular docking

For carrying out molecular docking studies, a grid box of 10 ų was generated using the maestro receptor grid generation wizard (Schrödinger) by specifying the binding (active) site residue Arg555 as no co-crystal ligand was available [49]. The ligands were docked using the ligand docking wizard of the glide module of Schrödinger with standard-precision (SP) mode initially. The molecules that were able to bind were subjected to filtration by applying the criterion of molecular weight (<500 Da) and number of rotatable bonds (<10) [63]. Molecules obtained from this filtration process were then subjected to docking with extra-precision (XP) mode [64].

2.4. MM-GBSA calculations

Molecular mechanics with generalised Born and surface area solvation (MM-GBSA) studies help to calculate the ligand's binding-free energy (∆G_bind) value. The calculations were mainly based on the summation of differences in the minimization (∆E_MM), solvation (∆G_Solv), and surface area (∆G_S.A.) energies of RdRp-ligand complex structure and free RdRp and ligand molecules [65]. The protein-ligand complex (PLC) from the docking studies was used to calculate the binding free energy. The analysis was carried out in the prime module of Schrödinger software. OPLS4 force field with dielectric surface generalized Born (VSGB) continuum solvation model was used. Based on the following formula, the binding free energy of the ligand was calculated:

$${∆G}_{bind}={∆E}_{MM}+{∆G}_{solv}+{∆G}_{SA}$$

(1)

2.5. Molecular dynamics studies

The molecular dynamics (MD) simulation was accomplished to determine the ligand molecules' stability, confirmation and intermolecular interaction with the target RdRp protein (PDB ID: 6M71) [66]. The time-dependent modification of the complexes was estimated over 100 ns using the Desmond module. The MD simulation was executed under thermodynamic conditions (applied volume, density, pressure and temperature). The complete system was annealed and equilibrated using ensembles.

The down-selected ligand complexes were imported and prepared using the protein preparation wizard of the Desmond module (Schrödinger LLC., NY, v2020). The complex was solvated using the transferable intermolecular potential with 3 points (TIP3P) model and the grid boundary dimensions was set to 10 ų [67]. The complex model was electrically neutralized with Na⁺/Cl⁻ ions, and built using a system builder wizard. The complete solvated model was minimized and molecular dynamics was carried out for 100 ns on all five selected complexes (including the reference drug remdesivir). The root means square deviation (RMSD), root mean square fluctuation (RMSF) and interaction plots were used to interpret the stability of appropriate complexes.

2.6. Entropy calculation for molecular dynamics trajectories

The binding free energy of the protein-ligand complex was determined using the gmx_mmPBSA tool [68,69]. In this calculation, the Poisson-Boltzmann surface area (PBSA) method with a dielectric model (ipb = 2) and a non-polar solvation model (inp = 1) were employed. The ionic strength of the surrounding medium was maintained at 0.15 M, and the temperature was set at 310 K. To calculate the entropy (-TΔS), the interaction entropy (I.E.) method was used [70]. This allowed us to evaluate the change in binding free energies with entropy contributions for the protein, ligand, and their complexes. The trajectories from the protein-ligand MD simulation in explicit water from the Desmond module were used to generate the gromacs trajectory file required for calculations using visual molecular dynamics (VMD) software.

Additionally, topology files for protein and ligand were obtained separately by converting the *.cms files to *.gro and *.top files using the intermol software [71]. Since the initial frames during the dynamic simulation are involved in the equilibrium, the frames after 50 ns were considered for analysis data. To speed up the procedure and to get better averaging, we ran five independent MMPBSA calculations (10 ns each) for every complex from 50-100 ns and reported the average and standard deviation of these calculations [72].

3. Results and Discussion

All the drug molecules (synthetic and approved) obtained from the ChEMBL database were docked into the binding site of the RdRP protein (PDB ID: 6M71) using the ligand docking wizard of the glide module of Schrödinger software. The test molecules were docked using SP mode and we obtained 1992 molecules, signifying that all test molecules have occupied the active site pocket. Larger molecules >500 Daltons find it harder to absorb as smaller molecules are more absorbable. Hence, we used a molecular weight (<500 Daltons) filter to filter out large molecules. The number of drugs that passed the filter was 1398, leaving behind 594 molecules. To further reduce the number of molecules, a filter corresponding to the number of rotatable bonds (<10) was applied, as the drugs having fewer rotatable bonds are acceptable (<10) [73,63]. With this filter, six hundred twenty-eight (628) molecules were obtained, which were further subjected to molecular docking using extra precision (XP) mode. All the molecules were binding to the active site with varying glide scores; hence, a cut-off score of -9.0 kcal/mol was considered to down-select the compounds. The list of molecules obtained after the cut-off is shown in Table 1.

3.1. MM-GBSA studies

These top 14 molecules obtained using a cut-off value of <-9.0 kcal/mol from molecular docking studies were subjected to MM-GBSA studies to evaluate the binding free energy of the ligand, and the molecules with the least score when compared to remdesivir were selected. The binding energy of remdesivir was found to be -40.32 kcal/mol, and we see that four molecules have lower binding energy. The 14 molecules (Table 2) were divided into three categories, i.e., low, moderate, and high binding energy, with 6, 4 and 4 each. As the last category of molecules showed the least energy compared to remdesivir, we further analyzed their interaction pattern (depicted in Figure 2 and Table 3) and extrapolated our study using these molecules for molecular dynamics studies.

3.2. Molecular dynamics results

The top ligands with the least binding energy from MM-GBSA studies were analyzed by molecular dynamics studies (Figure 3). The clinically used drug remdesivir was also subjected to dynamics simulations for comparative analysis, and showed a stable RMSD plot for the entire duration (5.0 Å-7.0 Å from the initial 5 ns until the end). Additionally, the fluctuations in the active site region of the protein were low, resulting in a stable RMSD and RMSF plots of the protein (Figure 3 and Figure 4). During this, significant interactions were observed by the residues of the palm domain that are necessary to bind with the RNA. The residual interactions include Asp618 (63%, water-mediated), Asp623 (70%, water-mediated), and Asp760 (90%, H-bond) (Figure 5). In the case of compound CHEMBL164344, the RMSD plot has deviations till 10.0Å for the initial 18 ns, during which strong interactions were observed with the residues of palm domain, Thr556, Asp623, and Arg624; later, these interactions gradually decreased, and the interactions with two residues of F domain, Asp452 (81%, H-bond) and Tyr455 (61%, π-π stacking) increased over the time. This has caused a decrease in the RMSD from 12.0 Å at 20^th ns to 7.5 Å at 100^th ns. Change of phenoxymethyl (as in CHEMBL164344) with biphenyl group (as in CHEMBL196676) results in significantly more interactions with the residues Asp452 (81%, H-bond) of F domain and Asp623 (85%, H-bond) of palm domain, as well as inconsistent water-mediated interactions with Asp760 (27%) (Figure 5). Initially, the interactions with Arg553, Thr556, and Asp623, as observed in the docked complex, remained for 18 ns, during which no significant deviations were observed (RMSD between 1.0 Å – 3.0 Å); later, the interactions with Asp452, Lys621 along with Asp623 had changed the ligand's conformation, resulting in the increase in the RMSD ranging between 8.0 Å and 7.0 Å which retained till 70 ns, and then a gradual decrease until the end of the simulation (last frame RMSD is 5.4 Å) (Figure 3). The N-methyl quinazoline-4-amine compound (CHEMBL4291724) shows greater stability in terms of its RMSD plot as equilibrium is attained after 20 ns until the end of the simulation with RMSD ranging between 10.0 Å and 12.0 Å (Figure 3). Most of the interactions are between the phosphonic acid groups and the residues of the palm domain, which include Asp623, Thr680, and Asp760, with a contribution of 50%, 39%, and 86%, respectively (Figure 5). In the case of CHEMBL387132, the RMSD plot is comparatively lower than all other complexes, as the overall deviations for 100 ns were between 3.0 Å and 6.4 Å. Additionally, strong interactions were observed with Asp623 for the first 30 ns and then with Asp760 of the palm domain, contributing to an overall interaction of 30% and 64%, respectively. A residue from the F domain, Lys545, has a total of 33% H-bond interactions with the ligand. (Figure 3 and Figure 5). The protein RMSD plot shows a similar pattern of deviations with all the hit ligands, indicating their stability with the ligands during the MD simulation. Similarly, the RMSF plot shows that interacting residues in the active site have acceptable fluctuations (<2.0 Å) (Figure 4). From this, we infer that the selected top ligands behave in a similar fashion to that of remdesivir, with the active site residues forming strong bonds, and most of the interactions by our “hit compounds” attributable to bisphosphonic acid groups. Therefore, there is a high probability that these compounds may emerge as potential RdRp inhibitors if evaluated in vitro.

3.3. Estimation of entropic contribution by gmx_MMPBSA

The molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) analysis helps to understand entropic contribution between protein and ligand during the molecular dynamic simulations [74]. It refers to the degree of randomness in a system and can guide to understand the entropic contributions of the ligand in the active site of the protein. Since MMPBSA analysis module is not available in Schrödinger software, the gmx_MMPBSA tool was used to determine the entropic contributions of the ligands, protein, and the protein-ligand complexes [68]. The change in free energy (ΔG) of the complex is then calculated using the following equation:

$$\Delta G_{binding}=G_{complex}-\left(G_{protein}+G_{ligand}\right)$$

(2)

In the above equation (1), the G_binding represents the free energy of the protein-ligand complex in water, while G_protein and G_ligand represent the free energies of the protein and ligand respectively in water. The free energy of the total system (ΔG_binding) can be obtained by subtracting the interaction entropy (I.E.) from the change in total energy (ΔE_MM) of the system. ΔE_MM is the summation of various change in energies such as van der Waals (ΔE_vdw), electrostatic columbic (ΔE_EL), electrostatic potential (ΔE_PB) and non-polar (ΔE_NPOLAR).

$$\Delta G_{binding}=\left(\Delta E_{MM}-T\Delta S\right)$$

(3)

Whereas,

$$\Delta E_{MM}=\left(\Delta E_{vdw}+\Delta E_{EL}+\Delta E_{PB}+\Delta E_{NPOLAR}\right)$$

(4)

Applying equation (3), the free energy of the systems (ΔG_binding) was calculated and tabulated in Table 4. From the results, it is clear that the compound CHEMBL196676 had entropy contributions less than remdesivir as their interaction entropies (I.E.) were 9.476 kcal/mol and 10.856 kcal/mol respectively, whereas, their ΔG_binding were -33.404 kcal/mol and -19.914 kcal/mol respectively. The higher negative values indicate the free energy of the complex is lower than the sum of individual free energies of the protein and ligand. Similarly, the I.E. values in the case of CHEMBL164344 (13.764 kcal/mol) and CHEMBL4291724 (15.428 kcal/mol) were slightly greater than remdesivir, however, due to lower total energy (ΔE_MM) the binding free energies were better than remdesivir (Figure 6). The compound CHEMBL387132 shows higher I.E. and ΔE_MM values, resulting in poor binding free energy. These results are in agreement with our molecular dynamics studies wherein the 4 hit compounds showed significant interactions with the active site residues except for CHEMBL387132.

3.4. Results of in silico predicted ADMET profiles

This study estimated the pharmacokinetic and toxicological predictions of the top-hit BP molecules. The predicted results showed satisfactory results for three compounds (CHEMBL4291724, CHEMBL387132, CHEMBL196676). The compound CHEMBL164344 showed a low Caco-2 permeability score, and the intestinal absorption is <30%, indicating the compound may have a low absorption profile. The other three compounds exhibited similar skin permeability; none acted as a P-glycoprotein (P-gp) substrate. The drug that acts as a substrate to P-gp showed low bioavailability, as evident from the studies of remdesivir [54]. The volume of distribution in steady state condition (Vd_ss) for all the compounds is less except for CHEMBL196676, which possessed better tissue distribution value than remdesivir. The fraction of the drug that remains unbound with blood plasma protein is essential to pass through cell membranes; therefore, the higher the fraction unbound (Fu) value, the higher its distribution [75]. All four compounds have better Fu values than the reference drug (0.005 Fu). Due to the presence of phosphonate groups, the compounds are likely to have poor permeability to the blood-brain barrier (BBB) or the central nervous system (CNS); however, CHEMBL387132 and CHEMBL196676 exhibited better permeability profiles than remdesivir. None of the molecules are likely to be metabolized by cytochrome P450, and all showed a better clearance rate. The maximum tolerated dose is higher than remdesivir and does not inhibit the human ether-a-go-go gene (hERG). The rat oral toxicity levels are higher than remdesivir, indicating high concentrations are required to cause toxicity. The other toxicity parameters are better compared to remdesivir (Table 5).

4. Conclusions

In the current study, an in silico investigation of BP-containing molecules obtained from the CHEMBL database was carried out. The preliminary studies produced 14 compounds with the best docking score compared to remdesivir, but from MM-GBSA analysis, four molecules were observed to be better than remdesivir with superior molecular dynamics. The agreement between entropy and dynamics results, combined with in silico predicted ADMET analysis, lead us to predict that three of these compounds (CHEMBL196676, CHEMBL164344 and CHEMBL4291724) are suitable candidates for synthesis and in vitro evaluation against RdRp. These candidates also have better water solubility, Caco₂ permeability, intestinal absorption, distribution, metabolism and excretion profiles. They have poor bioavailability and toxicity profile but are better than the reference drug, remdesivir, which is currently under clinical use.