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 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 averts the 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 (M^pro) 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 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 antiviral effects [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) inhibitors [41,42,43,44], these studies were looking for better alternatives to remdesivir. Few studies also reported the screening of analogues of different scaffolds such as quinolines [40], cytidines [45] and andrographolides [46]. It’s also common for repurposing efforts to propose drugs that are yet to be tested in vitro [13,32,38,39,40,47,48,49,50,51,52,53]. As RdRp is a key target, this paper looks at a class of small molecules called bisphosphonates (BPs) for five reasons: 1. The use of selected BPs is associated with a 3-5 fold reduction in the incidence of SARS-CoV-2 testing, COVID-19 diagnosis, and related hospitalization during the pandemic [54,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 different (in silico, in vitro and ex vivo) studies [21]; 4. Our recent study which has shown that alendronate, a BP drug, has promising in silico results compared to remdesivir [32]; and 5. BPs benefit the age groups that also experience severe forms of COVID-19. 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]. Such a range of immunomodulatory effects and the drug’s binding with RdRp inspired us to investigate BPs as prospective COVID-19 and/or PASC drug candidates [54,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 seven 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). To start screening BPs 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_SA) energies of RdRp-ligand complex structure and free RdRp and ligand molecules [65]. The protein-ligand complex 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 at a constant temperature of 310 K using Nosé -Hoover chain thermostat, and constant pressure of 1.013 bar using Martyna-Tobias-Klein barostat [67,68]. 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 ų [69]. 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 eight selected complexes (including the reference drug remdesivir, and the negative control cinnamaldehyde). The root mean 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 [70,71]. In this calculation, the molecular mechanism Poisson-Boltzmann surface area (MMPBSA) 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 [72]. 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 [73]. Since the initial frames during the dynamic simulation are involved in the equilibrium, the frames after 50 ns were considered for data analysis. 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 [74].

2.7. ADMET (Absorption, Distribution, Metabolism, Excretion, Toxicity) Studies

To determine the pharmacokinetic and pharmacodynamic profiles, the smiles format of the hit compounds was used and their parameters were predicted using the pkCSM web-tool [75]. All the absorption, distribution, metabolism, excretion and toxicity parameters obtained were compared with the standard reference drug remdesivir.

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 the SP mode and we obtained 1,992 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 molecular weight (<500 Daltons) to filter out large molecules. The number of drugs that passed the filter was 1,398. To further reduce this number, a filter corresponding to the number of rotatable bonds was applied, as the drugs having fewer rotatable bonds are acceptable [63,76]. With this filter, 628 molecules were obtained, which were further subjected to molecular docking using the extra precision (XP) mode.

All the molecules were binding to the active site with varying glide scores; hence, a cut-off docking score of -9.0 kcal/mol was chosen for down-selection as it is 2-3 times that of comparables (remdesivir -3.27 kcal/mol; favipiravir -3.44 kcal/mol; molnupiravir -4.93 kcal/mol). Concurrently we also enabled protomerization so that the uncharged BPs can be ionized at pH 7.0 ± 2.0. This is because BPs contain two phosphonate groups attached to a carbon or a nitrogen atom that can undergo ionization depending on the environment’s pH. For instance, alendronate, a bone resorption inhibitor used for treating osteoporosis exhibits multi-level ionization for the dissociation of 4-hydroxyl groups present on the phosphorous atom, resulting in a total of 4 dissociation constants at varying pH [77]. Similarly, at a pH ranging from 5 to 9, these BPs may exist as a mixture of protonated or ionized forms. Hence, considering the ionization behaviour at this stage is crucial to evaluate their binding mechanism with the target residues. Therefore, to improve the accuracy and reliability of our docking studies further, and to allow for a comprehensive exploration of ligand’s behavior in different biological environments, protomers at pH 7.0 ± 2.0 were generated using Epik module for the top 14 compounds that met the -9.0 kcal/mol docking score cut-off. This resulted in a total of 48 protomers, using which the XP mode of molecular docking was carried out. The list of molecules obtained after the cut-off and their corresponding best protomers are 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, and their respective protomers, 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. We see that seven candidates which are either uncharged molecules or charged protomers (highlighed in bold font in Table 2 along with their ChEMBL ID’s) have lower binding energy than that of remdesivir (-40.32 kcal/mol). For these 7 molecules with lower energy compared to remdesivir, we further analyzed their interaction pattern as depicted in Figure 2 and Table 3.

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 its corresponding protein were low, resulting in a stable RMSD and RMSF plots of the protein (Figure 3 and Figure 4). 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) (Supplementary Materials S1). 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 the 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. In the case of its protomer, the ligand has stable deviations for the initial 25 ns (RMSD 4.0 Å – 6.0 Å), later increasing to 16 Å till 45 ns and then attaining equilibrium until the end with RMSD ranging between 11 Å and 13 Å. The residues that significantly participated in the interactions include Tyr455 (73% π-π stacking, 68% H-bond), Lys551 (69%, H-bond), Arg553 (92%, H-bond), Arg555 (35%, H-bond), and Lys621 (33%, π-cation). It is noted that the protomer’s RMSF is similar to that of the uncharged molecule (Figure 4). 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%). 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 was observed until the end of the simulation (last frame RMSD is 5.4 Å) (Figure 3).

The 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 (Supplementary Materials S1). In the case of CHEMBL387132, a 2-amino isoquinoline compound, the RMSD plot is comparatively lower than 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, Supplementary Materials S1).

The 2-amino thiazole containing BP (CHEMBL4569308) showed a stable RMSD plot for initial 45 ns (between 1.0 Å and 2.0 Å), had a slight increment over the next few time frames and then attained equilibrium from 60 ns until the end with deviations between 2.5 Å and 4.0 Å. The stability is observed due to the significant involvement of thiazole ring in the interactions with Asp623 (99%, H-bond; 71% water-mediated H-bond) and Arg624 (89%, π-cation; 52% and 40% H-bond). Other residues such as Lys545 (68%, H-bond), Arg553 (76%, 82%, H-bond) and Arg555 (68%, 41%, H-bond) showed moderate interactions with the oxygen atoms of the compound. The imidazole compound (CHEMBL98211) attained stability from 20 ns until the end with RMSD between 5.5 Å and 7.0 Å. However, there were slight deviations observed between 45 ns-55 ns and 85 ns-90 ns with RMSD ranging between 4.0 Å and 7.0 Å. These deviations could be due to the intermittent water-mediated interactions with Asp623 (72%). The significantly contributing residues in the interactions include Lys545 (96%, H-bond), Arg555 (77%, 67%, H-bond) and Arg624 (73%, 49%, H-bond) (Supplementary Materials S1). Other residues such as Thr556 and Lys621 contributed with 52% of water-mediated and 48% of H-bond interactions, respectively. The imidazo[1,2-α]pyridine-4-ium compound (CHEMBL608526) exhibited the most stable RMSD plot among all the studied complexes. For the initial 30 ns, its RMSD was between 3.0 Å and 4.0 Å, later there was a gradual decrease until 60 ns, and then equilibrium was attained until the end with RMSD ranging between 1.0 Å and 2.0 Å. There were several residues that significantly contributed to the interactions during the simulation which include Asp450 (91%, H-bond), Lys545 (99%, 68%, H-bond), Arg553 (100%, 99%, H-bond), Arg555 (100%, 93%, H-bond), Asp623 (102%, water-mediated H-bond) and Arg624 (101%, 84%, H-bond interactions). The other residues such as Lys551, Thr556, Ala558 and Ser682 contributed moderately with 63%, 43%, 33% and 60% interactions, respectively. The negative control cinnamaldehyde showed highly unstable RMSD plot with high deviations (up to 90.0 Å), and interactions with <10% contributions (Supplementary Materials S1).

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 5 out of our top 7 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 good probability that these compounds may emerge as potential RdRp inhibitors if evaluated in vitro and/or ex vivo.

3.3. Estimation of Entropic Contribution by gmx_MMPBSA

The MMPBSA analysis helps to understand entropic contribution between protein and ligand during the molecular dynamic simulations [78]. 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 [70]. 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)

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

(3)

$$\Delta E_{MM}=\left(\Delta E_{vdW}+\Delta E_{EL}+\Delta E_{PB}+\Delta E_{NP}\right)$$

(4)

In the above equation (1), the G_complex 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 adding the interaction entropy (I.E. = - TΔS) to the change in total energy (ΔE_MM) of the system [72]. Δ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_NP). Applying equation (3), the free energy of the system (ΔG_binding) was calculated and tabulated in Table 4 and represented graphically in Figure 5.

The results reveal that the protomer of CHEMBL608526 has the least value (-75.62 kcal/mol), while CHEMBL196676 (-33.41 kcal/mol), CHEMBL164344 (-26.77 kcal/mol), CHEMBL4291724 (-26.38 kcal/mol) and CHEMBL4569308 (-22.97 kcal/mol) have values lower than remdesivir (-19.91 kcal/mol). The more negative values signify that the free energy of the complex is lower than the sum of the individual free energies of the protein and the ligand, as seen from equation (2). Therefore, these five candidates, highlighted in bold in Table 4, are worth evaluating in vitro and/or ex vivo. The two other candidates, CHEMBL387132 (-1.16 kcal/mol) and CHEMBL98211 (0.42 kcal/mol), have values close to zero and to the negative control cinnamaldehyde (-0.99 kcal/mol), but they may also be worth evaluating for the following reasons: 1. They perform as well as the 5 highlighted compounds in molecular dynamics studies, showing several interactions with the protein’s active site that are absent in the case of negative control (c.f. Figure 3, Supplementary Materials S1); 2. ΔG_binding values can vary substantially, e.g. CHEMBL164344’s uncharged state value (-26.77 kcal/mol) is twice that of its protomer (-12.65 kcal/mol), therefore it may be premature to dismiss these two compounds; and 3. CHEMBL98211 resembles the approved drug zoledronate which has been observed to benefit COVID-19 human patients as described in section 3.4 below. Although charged molecules are associated with electrostatic repulsions, higher degrees of freedom, and generally lower binding energies leading to more spontaneous binding, we have to exercise caution. This is because the relative contribution of the first term (enthalpy) versus the second term (entropy) in equation (3) is not straightforward in the case of proteins due to conformational effects in the latter and specific properties of a given system [79]. More studies are required with a range of molecules wherein the in silico predictions can be experimentally validated to determine whether free energy calculations for charged or uncharged molecules are better predictive and more useful for drug selection.

3.4. Results of in Silico Predicted ADMET Profiles

All 7 compounds have satisfactory pharmacokinetic and toxicity (ADMET) profiles compared to remdesivir. The intestinal absorption rate ranges from 34% to 72% for all the molecules except for CHEMBL164344 (16.77%). All compounds exhibited similar skin permeability; none acted as a P-glycoprotein (P-gp) substrate. A drug that acts as a substrate to P-gp implies low bioavailability, as evident from the studies on remdesivir [80]. The volume of distribution in steady state condition (Vd_ss) for all 7 compounds is less except for CHEMBL196676, as the latter 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, higher the fraction unbound (Fu) value, the higher is its distribution [81]. All 7 compounds had better Fu values than remdesivir (0.005 Fu), with CHEMBL608526 showing higher distribution pattern as its value was 0.811 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 slightly higher permeability profiles than remdesivir. None of the molecules are likely to be metabolized by cytochrome P450, and all showed a good 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 above results, summarized in Table 5, are very promising and indicate that these 7 compounds should be evaluated further, while two of them have a very close similarity to approved drugs. CHEMBL602586 is very similar to minodronic acid, while CHEMBL98211 is close to zoledronic acid (Supplementary Materials S1). Users of BPs, especially alendronate/alendronic acid and zoledronic acid, had “lower odds ratios (OR) of testing for SARS-CoV-2 infection (OR = 0.22; 95%CI:0.21–0.23; p<0.001) COVID-19 diagnosis (OR = 0.23; 95%CI:0.22–0.24; p<0.001), and COVID-19-related hospitalization (OR = 0.26; 95%CI:0.24–0.29; p<0.001)” in a recent study on “whether prior use of BPs is associated with reduced incidence and/or severity of COVID-19” [54,55]. These authors decided to investigate BPs (such as alendronate/alendronic acid and zoledronate) because pre-COVID-19 observational studies had reported “decreased in-hospital mortality for patients in the ICU” [82] and “reduced incidence of pneumoniae and pneumonia-related mortality in patients treated with amino-BPs versus controls” [83]. The in silico basis for these observations in people using alendronate/alendronic acid is provided by our previous work [32]; and by this follow up paper for zoledronate and a new drug minodronate. The latter drug and related ChEMBL ID’s (Table 4) are worthy of further investigation pertaining to COVID-19 and/or PASC. An in vitro/ex vivo analysis similar to McAuley et al. [21], or an in situ analysis similar to Thompson et al. [54], could be faster and useful ahead of in vivo studies and clinical trials.

4. Conclusions

In the current study, in silico investigation of BP-containing molecules from the CHEMBL database was carried out. Initial analysis identified 48 molecules with docking scores 2-3 times superior to remdesivir. Subsequent MM-GBSA analysis led to the selection of 7 molecules with notably better free energy scores than remdesivir. These 7 candidates underwent further assessment of their dynamics profiles and entropy calculations, resulting in the identification of 5 promising candidates for in vitro evaluation, comprising of 3 uncharged and 2 charged molecules. It is important to consider the charged nature of a molecule to assess protein binding and interactions in biological systems; while its uncharged state could be relevant in the content of stability and metabolism. Notably, our investigation, exemplified by CHEMBL164344, revealed that despite charged and uncharged molecules having similar MM-GBSA scores, dynamic behavior, and total internal energies, their entropy contributions can differ significantly (for example, CHEMBL164344’s values were -12.65 kcal/mol for charged and -26.77 kcal/mol for uncharged). Whether it is important to only consider the charged state for down-selection remains to be seen from future experimental and clinical trial observations. The agreement between entropy and dynamics results, combined with in silico predicted ADMET analysis, lead us to predict that seven compounds (CHEMBL196676, CHEMBL164344, CHEMBL4291724, CHEMBL608526, CHEMBL4569308, CHEMBL98211, CHEMBL387132) are suitable candidates for synthesis, in vitro and in vivo evaluation against RdRp. Of these seven, CHEMBL608526 very closely resembles the approved drug minodronate, while CHEMBL98211 resembles the approved drug zoledronate.