Next Article in Journal
Rapid Analysis of Compounds from Piperis Herba and Piperis Kadsurae Caulis and Their Differences Using High-Resolution Liquid–Mass Spectrometry and Molecular Network Binding Antioxidant Activity
Next Article in Special Issue
SAR, Molecular Docking and Molecular Dynamic Simulation of Natural Inhibitors against SARS-CoV-2 Mpro Spike Protein
Previous Article in Journal
Chemical Composition, and Antioxidant and Antimicrobial Activity of Oregano Essential Oil
Previous Article in Special Issue
Structure–Activity Studies on Bis-Sulfonamide SHIP1 Activators
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Favipiravir Analogues as Inhibitors of SARS-CoV-2 RNA-Dependent RNA Polymerase, Combined Quantum Chemical Modeling, Quantitative Structure–Property Relationship, and Molecular Docking Study

by
Magdalena Latosińska
and
Jolanta Natalia Latosińska
*
Faculty of Physics, Adam Mickiewicz University, Uniwersytetu Poznańskiego 2, 61-614 Poznań, Poland
*
Author to whom correspondence should be addressed.
Molecules 2024, 29(2), 441; https://doi.org/10.3390/molecules29020441
Submission received: 31 October 2023 / Revised: 9 January 2024 / Accepted: 12 January 2024 / Published: 16 January 2024

Abstract

:
Our study was motivated by the urgent need to develop or improve antivirals for effective therapy targeting RNA viruses. We hypothesized that analogues of favipiravir (FVP), an inhibitor of RNA-dependent RNA polymerase (RdRp), could provide more effective nucleic acid recognition and binding processes while reducing side effects such as cardiotoxicity, hepatotoxicity, teratogenicity, and embryotoxicity. We proposed a set of FVP analogues together with their forms of triphosphate as new SARS-CoV-2 RdRp inhibitors. The main aim of our study was to investigate changes in the mechanism and binding capacity resulting from these modifications. Using three different approaches, QTAIM, QSPR, and MD, the differences in the reactivity, toxicity, binding efficiency, and ability to be incorporated by RdRp were assessed. Two new quantum chemical reactivity descriptors, the relative electro-donating and electro-accepting power, were defined and successfully applied. Moreover, a new quantitative method for comparing binding modes was developed based on mathematical metrics and an atypical radar plot. These methods provide deep insight into the set of desirable properties responsible for inhibiting RdRp, allowing ligands to be conveniently screened. The proposed modification of the FVP structure seems to improve its binding ability and enhance the productive mode of binding. In particular, two of the FVP analogues (the trifluoro- and cyano-) bind very strongly to the RNA template, RNA primer, cofactors, and RdRp, and thus may constitute a very good alternative to FVP.

Graphical Abstract

1. Introduction

1.1. Motivation of the Research

Over the last two centuries, there have been several pandemics/epidemics caused by different human-infective RNA viruses: in 1889–1890 (Asiatic/Russian flu, A subtype H2), in 1918–1920 (Spanish flu, H1N1), in 1957–1958 (Asian flu, H2N2), in 1968–1969 (Hong Kong flu, H3N2), in 1977–1979 (Russian flu, H1N1), since 1981 (AIDS, HIV), in 2009–2010 (Swine flu, H1N1/09), in 2002–2003 (SARS, SARS-CoV), and in 2012–2013 (MERS, MERS-CoV). The emergence of Swine flu 13 highlighted the urgent need for effective antiviral therapy against RNA viruses due to their pandemic nature. The coronavirus disease-2019 (COVID-19), a life-threatening infectious disease caused by the novel severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) [1,2], which emerged in Hubei Province in Central-Eastern China in December 2019, caused panic and accelerated these efforts. Since the global spread of the COVID-19 pandemic, many different SARS-CoV-2 variants (Alpha, Beta, Gamma, Delta, Epsilon, Etha, Iota, Kappa, Lambda Omicron, Theta, and Zeta) and subvariants, in particular, so-called variants of interests (XBB.1.5, XBB.1.16, EG.5, and JN.1) and variants under monitoring (BA.5, CH.1.1, XBB, and BA.2.86) appeared, which raised concerns about the possibility of prolonging the current pandemic and the occurrence of subsequent ones. Therefore, SARS-CoV-2 variants are monitored on an ongoing basis by the World Health Organization Technical Advisory Group on SARS-CoV-2 Virus Evolution [3].
Three recent pandemics, influenza [4], HIV [5], and SARS-CoV-2 [6], not only have common zoonotic origins but also are caused by RNA viruses (negative-sense single-stranded RNA ssRNA (−) or positive-sense single-stranded RNA ssRNA (+)), belonging to the kingdom Orthornavirae (realm Riboviria). The HIV group antigens (Gag) and influenza virus matrix (M1/M2) proteins also have even evolved from a common ancestor protein [7]. Two pandemic viruses, influenza and coronavirus, infect the respiratory tract, share similar symptoms, and use surface proteins to infect the host cells [8]. However, influenza uses surface glycoproteins hemagglutinin (HA) and neuraminidase (NA) for infection [9], while SARS-CoV-2 (similarly to MERS-CoV, SARS-CoV-1) uses the spike (S) protein [10]. In the case of HIV-1, the mechanisms of entry are different, but the strategy is similar to the one used by SARS-CoV-2 [11]. All three viruses use the viral RNA polymerase to express their proteins, although only SARS-CoV-2 has a proofreading mechanism [12], which has been considered an argument against the possibility of rapid mutation. Furthermore, in the case of the Ebola virus (Zaire Ebolavirus, EBOV), which is an ssRNA (−) enveloped virus from Riboviria [13], the trimeric transmembrane glycoprotein (GP) that plays a key role in infection mediates EBOV attachment and entry into host cells. However, EBOV, which is transmitted through body fluids and causes flu-like symptoms, was the cause of the Ebola hemorrhagic fever epidemic in West Africa in 2013–2016. Fortunately, EBOV does not mutate rapidly, which, contrary to initial fears, reduces the risk of a pandemic. In any case, viruses belonging to the same group of ssRNA viruses with a similar structure are considered potential culprits for the next global pandemic. In light of this, effective therapy targeting RNA viruses is highly desirable.
Most anti-viral drugs used in previous pandemics either inhibit the virion’s M2 ion channel, such as amantadine (Gocovri/Symadine/Symmetrel, 1964) and rimantadine (Flumadine, 1964), or inhibit the viral neuraminidase (e.g., oseltamivir (Tamiflu, 1997), zanamivir (Relenza, 1993), peramivir (Rapivab, 2000), and laninamivir (CS-8958, 2009). They prevent the virus from entering (entry/uncoating phase) and exiting cells (release phase). These treatments have demonstrated some effectiveness; however, there is potential for greater effectiveness with more specific treatment methods. The search for more effective drugs with a broad spectrum of activity and an innovative mode of action was dictated by three main factors: poor performance against new variants, the emergence of resistance to antiviral drugs (e.g., amantadine, rimantadine [14,15,16], and oseltamivir [17,18,19]), and the need to treat previously unknown or untreatable diseases. The most current first-line anti-viral therapies either focus on suppressing the cytokine storm (e.g., tocilizumab, a humanized anti-human IL-6 receptor antibody) or act as broad-spectrum antivirals (e.g., favipiravir, ribavirin, and remdesivir) which inhibit the viral RNA polymerase. The priority is to develop or improve broad-spectrum antivirals for the treatment of epidemic diseases caused by emerging or re-emerging viruses.

1.2. Favipiravir–State of the Art

Favipiravir (6-fluoro-3-hydroxypyrazine-2-carboxamide, FPV, T-705;), Figure 1, is an active pharmaceutical ingredient discovered by a Japanese company Toyama Co. Ltd. (now Toyama Kagaku Kōgyō, Toyama, Japan), which was approved for treating influenza strains unresponsive to adamantanes or neuraminidase inhibitors [20].
Food and Drug Administration (FDA) phase II and III clinical trials demonstrated the safety of FVP in humans and better efficacy against the influenza virus than oseltamivir [21]. In influenza virus-infected patients treated with FVP, the high barrier to resistance was observed [22]. FVP has been shown to have broad antiviral activity at higher concentrations against many other RNA viruses [23,24,25,26], such as OTV-resistant influenza A, B, and C viruses; as well as flavi-, alpha-, filo-, bunya, and arena-noroviruses [27,28]; Ebola [29]; Lassa virus [30]; West Nile Fever [31]; Zika [32]; Rift Valley fever [33]; Yellow fever [34]; Crimean-Congo hemorrhagic fever (CCHF) [35]; Nipah tick-borne encephalitis [36]; rabies [37]; and Argentine hemorrhagic fever (Junín) [38], including those for which there is currently no antiviral treatment available. After the identification of the new virus, SARS-CoV-2, FVP was among the first medications screened for its effectiveness and safety. Due to its mechanism of action, preclinical results, and high safety in humans, FVP shows promise in treating human diseases caused by various RNA viruses, including SARS-CoV-2. According to the ClinicalTrials.gov database, there are 67 completed, ongoing, or planned clinical trials around the world (Europe, North America, and Australia) assessing the efficacy and safety of FPV in monotherapies and combination therapies in the treatment of COVID-19. (Another 12 are efficacy studies for Influenza (7), Ebola (3), Lassa (1), and CCHF (1).) Currently, FVP is registered as a drug for COVID-19 in China, India, and Russia and is approved in Turkey, Hungary, Serbia, the KSA, Thailand, Egypt, Bangladesh, Pakistan, Jordan, and Saudi Arabia [39,40,41]. Oral doses of FVP have been shown to be safe in outpatients with mild to moderate infections. Furthermore, the reduction in viral load and improvement in radiological and clinical outcomes are significant [24,25]. Compared to Molnupiravir (MOL), Nirmatrelvir (NMVr) [18,19], or Remdesivir (RDV) [42], FVP exhibits reduced potency against SARS-CoV-2 in vitro but similar effectiveness in animal models.
Certain studies have indicated a potential association between FVP and oxidative stress/cardiotoxicity, [43] hepatotoxicity [44], teratogenicity [45,46,47], and embryotoxicity [48]. However, other anti-viral drugs, including MOL, NMVr, or RDV, also exhibit hepatotoxicity [49], nephrotoxicity [50], and mutagenicity [51]. Inconsistencies in clinical trial results, despite the synergistic effect of combined therapy with Oseltamivir (OSE) [52,53], Ivermectin (IVM) [54], Ribavirin (RBV) [55], Remdesivir (RDV), and Tocilizumab (TOZ) [56] in the fight against H1N1, H3N2, H5N1, and SARS-CoV-2 [57], have caused the use of FVP as a regular drug to remain under investigation. Even RDV, although introduced as a drug for the treatment of COVID-19 after obtaining its first emergency use authorization in May 2020 in the United States [58] and then later in Japan, has faced much criticism. The FDA approved RVD for the treatment of hospitalized patients, and the EU Commission granted conditional authorization, but the European Medicines Agency (EMA) approved it only for patients 12+ suffering from pneumonia who require oxygen supply [59]. The only anti-viral drug approved by the FDA is Pfizer’s PaxlovidTM, which is the combination of Nirmatrelvir (NMVR)—SARS-CoV-2 3-CL Mpro inhibitor and Ritonavir (RTV)—protease inhibitor. Thus, its mechanism of action is completely different. However, its use carries many risk factors, including severe liver impairment or liver disease. Recent studies show that SARS-CoV-2 has acquired phenotypic resistance to RDV and PaxlovidTM [60] but not to FVP. Thus, it is highly desirable to seek more suitable FVP alternatives, keeping in mind the issues mentioned above.

1.3. Favipiravir—Unique Mechanism of Action

Due to its relatively simple structure, FVP seems to be easily modifiable while maintaining its unique mode of action. It is therefore an attractive basis for investigating new drugs. However, the exact molecular mechanism behind the broad-spectrum antiviral activity of FVP has not yet been fully elucidated. It is known that FVP, a pro-drug, must be first converted by hypoxanthinguanine phosphoribosyl transferase (HGPRT) to FVP-ribosyl 5′-monophosphate (FVP-RMP) and then metabolized to the FVP-ribofuranosyl-5′-triphosphate (FVP-RTP) by cellular kinases [20,21] to inhibit viral RNA-dependent RNA polymerase (RdRp), as shown in Figure 2.
RdRp, a key enzyme regulating both the replication and transcription of viral RNA, is a component of the so-called minimal core of viral RNA replication. FVP-RTP acts as a purine nucleoside analogue, which pairs with either cytosine or uracil. It competitively inhibits the viral RdRp polymerase substrate, leading to chain termination [61] and ultimately inhibits the transcription and replication of the viral genome RNA. Another possible mechanism of its action is lethal mutagenesis, resulting in an increased frequency of, primarily, guanosine to adenine (G→A) and, secondarily, C→U mutations, which produce non-infectious progeny during replication [62]. In this case, the therapeutic effect of FPV will result from the accumulation of mutations in the replicated RNA of nascent viruses, which leads to a cascade of mistakes, the so-called “error catastrophe”, and loss of the virus’s ability to reproduce. The malfunction of RdRP in the presence of FPV depends on the intracellular phosphorylation of the drug to its active form (FPV triphosphate), a false nucleoside that is built by the viral RdRP into the nascent viral RNA, resulting in a “defective”, mutated RNA. Since RdRp structure and function are conserved among RNA viruses and no RdRp homolog has been found in human cells, its selection as a target is highly advantageous for the development of broad-spectrum antivirals. Furthermore, it was observed that the conserved active site of RdRp does not mutate as easily as other targets, such as the Spike protein. However, FVP is considered a poor substrate for HGPRT, and the requirement for its conversion to the active metabolite RTP is, in fact, a limiting factor for its antiviral activity [63]. Viral resistance to FVP may result from depletion or a lack of HGPRT [44]. In HGPRT-deficient cells, FVP is completely devoid of antiviral activity. To overcome the inadequate activation of FVP, the utilization of ribonucleoside pro-drugs, as well as di- and triphosphate analogues, seems to be a viable solution, which increases the biodistribution and therapeutic efficacy of FVP.

1.4. FVP Analogues—A Research Hypothesis and Concept

The motivation for our study is the hypothesis that new FVP analogues may provide more efficient nucleic acid recognition and binding processes while reducing side effects. Based on our previous observations from the analysis of FVP binding modes in solid, pre-catalytic, and active forms [64], we proposed a set of the FVP analogues (pro-drugs) along with their triphosphate forms as novel SARS-CoV-2 RdRp inhibitors. The pyrazine heterocyclic ring and three functional groups (halogen/CH3/CF3, hydroxide, and amide) in the FVP analogues allow them to engage in weak non-covalent interactions (e.g., hydrogen-bonds, van der Waals, steric, and stacking). The presence of the three donor (amide and hydroxide) or five/eight acceptor atoms (i.e., two aromatic ring nitrogens, oxygens of amide and hydroxide, and halogen/CF3/CN) in one molecule facilitates the formation of hydrogen bonds in the parental FVP. According to the Etter rule [65], FVP analogues can theoretically realize up to twenty-four different types of hydrogen bonds, which makes them highly appealing for pharmaceutical purposes. Enzymatic conversion to the active ribofuranose-5′-triphosphate (RTP) form leads to an increase in the number of oxygen atoms while facilitating the establishment of hydrogen bonds. However, both of these factors can be modified by altering the substituents at the C(6) position of the pyridazine ring. The key question is whether and to what extent alteration of the molecular structure via substituents will help in the effective recognition of FPV analogues, their incorporation into the RNA strand, and their binding to the RdRp. The effective screening of the halogenated and non-halogenated (R = I, Br, Cl, CF3, H, CH3, CF3, or CN) analogues of FVP-RTP is the objective of the current research. Using the Quantum Theory of Atoms in Molecules (QTAIM), Quantitative Structure–Property Relationship (QSPR), and Molecular Docking (MD) approaches, the differences in their reactivity, toxicity, binding efficiency, and ability to be incorporated by SARS-CoV-2 RdRp were assessed. These methods, supplemented by new global indices describing relative reactivity and new quantitative methods for the estimation and visualization of the differences in the binding modes of individual analogues with RdRp, provide insight into the set of desirable characteristics responsible for the inhibition of the SARS-CoV-2 RdRp. Molecular docking studies demonstrate the high binding affinity of the FPV analogues to SARS-CoV-2, indicating their potency as antiviral drugs against COVID-19. The proposed modification of the FVP structure seems to improve its binding ability to SARS-CoV-2 RdRp and enhance productive binding modes. If sufficient efficacy in inhibiting viral replication in cell culture is established, they could be explored as potential drugs against COVID-19. Our method for quantifying differences in binding mode holds promise for guiding future research on new anti-SARS-CoV-2 agents.

2. Results and Discussion

2.1. Characteristics of the Candidate Ligands

2.1.1. Physicochemical Profile (ADMET) and Key Pharmacokinetic Parameters

FVP analogues (Figure 1, R = I, Br, Cl, CF3, H, CH3, CF3, or CN) were selected as candidate ligands for SARS-CoV-2 RdRp and potential pro-drugs. The physicochemical profile parameters that describe the pharmacokinetic behavior of the candidate ligands and known drugs (MOL, RVD, and the only registered drug, PaxlovidTM, components: NMVr and RTV) have been evaluated in Table 1 and Table 2.
The molecular weights of the ligands (MW) ranged from 139.11 to 265.01 g/mol and did not exceed 500 g/mol. According to Swiss ADME [66], the predicted consensus LogP (lipophilicity) for the ligands ranged from −0.83 to 0.52. The CN analogue, like FVP, showed low and negative lipophilicity, while the CF3 analogue showed high and positive lipophilicity (closer to optimal for drugs). The substitution of –F with –CF3 increases the hydrophobicity (lipophilicity) of the ligand. Thus, the alteration of lipophilicity has the potential to modify hydrophobic targets. Indeed, in medicinal chemistry, –CF3 is often used as a substituent due to its strong electron-withdrawing nature, poor polarizability, and broad hydrophobic domain [67] (e.g., Tecovirimat [68], Doravirine [69], and Tipranavir [70]). However, ADMETlab 2.0 [71], which uses a different model, predicts negative LogP for both above-mentioned ligands. The predicted water solubility index (Solubility, SILICOS-IT) suggests that all the candidate ligands, like FVP, should be highly water soluble. Moreover, none of them showed lead-likeness violations and no violation of drug-like rules (Lipinski, Egan, or Veber). The predicted synthetic accessibility score (SAS) for the ligands ranges from 1.73 to 2.43, indicating the ease of synthesis of these compounds.
All the ligands showed high gastrointestinal absorption (GI), a key parameter in assessing the in vivo performance of an orally administered drug formulation. Abbot’s Bioavailability score for all ligands placed them within the 55% probability class. The topological polar surface area (TPSA), which describes the passive molecular transport across the membranes, is 88.84, except for the CN analogue, which showed a significantly higher TPSA of 112.63. However, since TPSA does not exceed 140, it is still an optimal level. Overall, the physicochemical profile of the candidate is better than those of MOL, RDV, NMVr, and RTV. The candidate ligands show no Pan Assay of Interference Structures (PAINS) or Structural (BRENK) alerts. None of the ligands are expected to cross the blood–brain barrier (BBB), inhibit cytochrome P450 (CYP) isoforms or have a potential to be a substrate of multidrug resistance protein (permeability glycoprotein, P-gp). The risk scores for carcinogenicity or genotoxic mutagenicity are very low for all of them, while the risk scores for hepatotoxicity (H-HT) and drug-induced liver injury (DILI) in humans are high. However, the CF3 analogue showed two times lower hepatotoxicity than FVP. The probability of genetic toxicity (Salmonella typhimurium reverse mutation assay, i.e., AMES assay) and the probability of rat acute oral toxicity (LD50) are relatively low for the CN and CH3 derivatives. Overall, the AMDE toxic profile for candidate ligands is non-inferior to FVP and better than MOL and RVD in terms of TPSA, genotoxic mutagenicity, genetic toxicity, H-HT, and DILI. The NVMr shows smaller toxicity, but it is combined with highly toxic RTV.
However, there is no set of ideal pharmacokinetic parameters that a given drug candidate should exhibit, as it depends on the specific requirements of the target. Therefore, our objective is to identify drugs that possess an optimal pharmacological profile, achieving the desired biological activity with minimal side effects (the last rows of Table 1 and Table 2 show the values that can be considered optimal).

2.1.2. Reactivity Profile (Quantitative Structure–Property Relationships)

Each ligand (potential pro-drug) and its active form, ribofuranosyl-5′-triphosphate (RTP), were constructed and optimized at the B3LYP/6–311G(d,p) level of theory. The enol tautomer (with intramolecular OH⋯O hydrogen bond), although more stable than the keto tautomer in solids and solutions, is further excluded by ribofuranosyl substitution. The planar conformation of the pyrazine moiety in active forms is maintained and supported by the intramolecular NH⋯O hydrogen bond in all FVP analogues, regardless of the type of substituent. The optimized molecular geometries served as the initial configurations for further QSPR and MD studies.
Determining the impact of molecular structure on reactivity is essential when designing the ligands with desired properties. A powerful tool for describing the chemical reactivity of the ligands based on their structure is the frontier molecular orbital (FMO) theory, which conceptualizes chemical bonding and reactivity in terms of the interactions between frontier orbitals. The highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) are very useful for assessing the chemical reactivity of molecules. LUMO accepts electrons, and its energy corresponds to an electron affinity (EA), while HOMO donates electrons, and its energy is related to ionization potential (IP). Low IP and high EA correspond to high nucleophilic and high electrophilic properties, respectively. Therefore, the HOMO-LUMO gap is useful for predicting charge transport and ligand. Theoretical global indices of ligand reactivity, such as absolute electronegativity, χ; absolute hardness, η; electrophilicity index (reactivity), ω; softness, S; electro-donating power, ω; electro-accepting power ω+; net electrophilicity, Δω; and a maximum number of electrons transferred in a chemical reaction, ΔNmax, provide further details regarding the reactivity of the ligands. The HOMO, LUMO, and global reactivity indices for the ligands were evaluated at MP2/6-311G(d,p) and M062X/6-311G(d,p) levels of the theory (in the gas phase—single molecule and aqueous solution) in Table 3 and Table 4.
Comparing reactivity parameters makes sense in relation to groups of structurally similar compounds; therefore, analyzing MOL, RVD, NVMr, or RTV will not bring anything in this respect.
The ligands studied can be ordered according to the decreasing HOMO-LUMO gap as follows: CF3 > H > CN > F > CH3 > Cl > Br > I. The results show that regardless of the phase (gas or aqueous solution), the CF3 and I analogues have the highest and lowest stability, respectively. The CF3 and CN analogues, which have a higher HOMO-LUMO energy gap, are more stable and are therefore chemically harder than the other ligands. However, each of these ligands is slightly less stable in aqueous solutions than in the gas phase.
The ordering of the ligands according to decreasing absolute electronegativity, χ, a measure of the ligand’s ability to attract electrons to itself, is correlated with the electronegativity of the halogen: CH3 < H < I < Br < Cl < F < CN < CF3. The CF3 group has a significantly strong electronegativity, typically intermediate between that of F and Cl, while the CN bond is strongly polarized toward nitrogen and more electronegative than Cl. Surprisingly, among the FVP analogues, the CF3 derivative is more electronegative than the CN derivative. The ordering of the ligands according to decreasing global hardness, a measure of the ligands’ resistance to change its electronic configuration, is as follows: CF3 > H > CN > F > CH3 > Cl > Br > I. The very high value of the absolute hardness for the CF3 derivative indicates its high degree of stability and low reactivity. The χ value describes the tendency to donate/accept electrons, while η measures the ease with which this can occur, which for CF3, are high and low, respectively.
The most important descriptor measuring electrophilic power (capacity of an electrophile to accept the maximal number of electrons in a neighboring reservoir of electron pool) is the global electrophilicity index, ω. Its values for the ligands are in the range of 1.318–1.885 eV and 1.323–1.670 eV in the gas phase and aqueous solution, respectively, as can be seen in Table 1. The parameter ω, which actually measures the reactivity of the ligand, revealed the following trend: CN > CF3 > F > Cl > Br > I > H > CH3. Therefore, reducing the inductive electron-withdrawing effect (F > Cl > Br > I) and electron-donating effect via resonance (F > Cl > Br > I) leads to a decrease in electrophilic activation (ω = 1.676, 1.648, 1.606, 1.592 eV for F, Cl, Br, and I, respectively). Very high reactivity values describing the system’s tendency to acquire electrons from the environment are observed for the CN and CF3 derivatives, while very low values for the H and CH3 derivatives. Thus, the CN and CF3 analogues seem most promising because highly electrophilic reagents lead to low substrate selectivity, which means they can inhibit a wide range of RdRps, not just SARS-CoV-2. Both analogues also have the highest local electro-donating power, ω+, electro-accepting power, ω, and overall electrophilicity, Δω. A larger ω+ value for CN corresponds to its better ability to accept charge, whereas a smaller ω value for CN makes this ligand a better electron donor. However, the unusually low-lying LUMO level for the CN and CF3 analogues suggests their easy participation in molecular reactions with nucleophiles, and the low-lying HOMO level for the CN and CF3 analogues suggests their easy participation in molecular reactions with electrophiles. LUMO is even slightly lower, while HOMO is slightly higher in the aqueous solution than in the gas phase, while ω+ and ω are higher in the gas phase.
Two new reactivity descriptors, the so-called relative electron donation power, R+, and the relative electro-accepting power, R, have been defined by us.
The relative electro-donating power, R+, is the quotient of the electro-donating power of the tested ligand and the reference ligand:
R + = ω l i g a n d + ω r e f e r e n c e   l i g a n d +
Similarly, the relative electro-accepting power is the quotient of the electro-accepting power of the tested ligand and the reference ligand:
R = ω l i g a n d ω r e f e r e n c e   l i g a n d
Both parameters, R+ and R, describe the ability of the ligands to accept and donate charge, as can be seen in Table 3 and Table 4.
The orderings of the ligands in descending order of R+ and R- are as follows:
CN > CF3 > F > Cl > I > Br > H > CH3 and CN > CF3 > F > Cl > Br > I > H > CH3 (the gas phase)
and
CN > CF3 > F > Cl > Br > I > H > CH3 and CF3 > CN > F > Cl > Br > I > H > CH3 (aqueous solution).
Both parameters, R+ and R, allow for the classification of ligands based on the reference ligand, as can be seen in Figure 3.
The R values in the gas phase and the aqueous solution are similar, while the R+ values differ significantly. Based on these parameters, the two ligands CN and CF3 are better as both acceptors and donors than FVP. These features are maintained in the aqueous solution (pH = 7, i.e., close to body pH). Moreover, the conclusions remain the same regardless of the calculation method used. However, the MP2 function appears to be more sensitive to changes in reactivity, as shown in Figure 3.
The degree of interaction, DOI [65], calculated for the most promising ligands in Table 5 confirms the observations regarding the differences between the two analogues CN and CF3.
The DOI characterizes the strength of an atom’s attachment to its molecular neighborhood, i.e., the degree of electron density sharing between the atom and its surroundings. The highest DOI was obtained for the acceptors –N(4) (39.26% in FVP, 37.964 for CN, and 37.963 for CF3), followed by –C(3) (20.942% for FVP, 20.651% for CN, and 20.300% for CF3), and –C(2) (22.819% for FVP, 22.352% for CN, and 22.349% for CF3). The lowest DOI was obtained for –NH2 (1.053% for FVP, 1.18% for CN, and 1.024% for CF3), –OH (2.626% for FVP, 2.449% for CN, and 2.571% for CF3) acting as a donor, and =O (0.141 for FVP, 0.143 for CN, and 0.156 for CF3%). The very low DOI for the –NH2 group is maintained regardless of the ligand type, making it highly suitable for binding ligands to the RNA strand. Based on the DOI parameters, the CF3 ligand appears more attractive than CN. After ribofuranosyl substitution, any R substituent loses the ability to share electron density with its surroundings in favor of RTP, and therefore, RTP plays a dominant role in the formation of bindings with the protein.

2.2. Characteristic of RNA-Directed RNA Polymerase

The genomic arrangement of SARS-CoV-2 is primarily composed of 4 structural proteins, nucleocapsid protein (N), spike protein(S), envelope protein (E), and membrane protein (M); 16 non-structural (NSPs); and 9 accessories (ORFs). Therefore, the SARS-CoV-2 virus consists of many proteins that can mutate rapidly. Two-thirds of the viral genome is occupied by the replicase gene referred to as two Open Reading Frames (ORFs), ORF 1a and ORF1ab, which encode the non-structural proteins (NSPs), the so-called pp1a and pp1ab polyproteins, respectively. The non-structural protein pp1ab includes RNA-directed RNA polymerase (RdRp), so-called nsp12 (chain 4393–5324). RdRp is a core component of viral replication and transcription [26] and exhibits significant catalytic activity, but only with the help of other cofactors: nsp7 and nsp8 [26,27]. Thus, nsp12-nsp7-nsp8 is defined as the minimal core component for viral RNA replication. It has been observed that the conserved active site of RdRp does not mutate as easily as other targets, such as the S protein.
The palm subdomain of SARS-CoV-2 RdRp (residues 585–625 and 680–807) forms the catalytic core of the polymerase, which contains the four highly conserved motifs (A–D). The principal target for SARS-CoV-2 is the active site of the RdRp polymerase, which is formed by two catalytic motifs: A, composed of the residues from 611 to 626, and C, containing residues from 753 to 767 [72]. In general, different RdRp polymerase inhibitors bind more or less strongly with the following residues: ASP760, ASP761, GLY616, TRP617, ASP618, TYR619, PRO620, LYS621, CYS622, LEU758, SER759, ALA762, ALA797, LYS798, CYS799, TRP800, HIS810, GLU811, PHE812, CYS813, SER814, and GLN815 [73,74]. Three hydrophilic and polar residues, ASP618, ASP760, and ASP761, play a key role in SARS-CoV-2 RdRp inhibition [74]. ASP618 is the most conserved residue in viral RdRp and, together with two strictly conserved residues, ASP760 and ASP761, is responsible for the formation of the RdRp catalytic center. ASP623 is involved in a hydrogen bond with the 2′-OH group of the nucleoside triphosphate and therefore appears to be important in sugar selection. The neutral SER759 residue is involved in the positioning of the priming nucleotide [75,76], while the hydrophilic and polar LYS798 stabilizes the RdRp core [77]. The multi-subunit RdRp binds nucleotide triphosphate substrates that enter the main enzyme channel via a hydrophilic cluster formed by the polar residues LYS545, ARG553, and ARG555. When the substrate enters the active site of the enzyme, the complex is formed.

2.3. Binding Modes of the Native and Candidate Ligands to RdRp in Different States

The procedure used to dock FVP-RTP analogues (Figure 1, R = I, Br, Cl, CF3, H, CH3, CF3, or CN) to the SARS-CoV-2 RdRp was nearly identical to that previously described [64,78,79]. The binding site (cavity) was identified, and the search space was defined as a subset region of about 9.0–15.0 Å.
Knowing that FVP-RTP can effectively mimic either guanosine or adenosine and bind to cytosine or uracil [66,80], both potential possible methods to bind to the SARS-CoV RdRp were explored. This may help resolve the source of bias in the spectrum of mutations induced by FVP, which is likely to be competition with adenosine and guanosine during nucleotide incorporation. Furthermore, the differences between desirable/productive and undesirable/unproductive binding modes were analyzed. To evaluate the quality of the docking process, we performed a redocking task. In each case, the actual ligand was removed from the parental structure and redocked in its own binding site. The redocking protocol was considered successful when the root-mean-square deviation (RMSD) of the pose relative to its conformation in the parental structure did not exceed 3 Å. The binding mode of the native ligand to the RdRp was described in detail. Then, in the same way, the new ligand was docked into the rigid protein structure, and its binding mode was characterized.

2.3.1. Pre-Catalytic State—Productive Mode I (Binding to Cytosine and Stacking to Adenosine)

Binding Mode of the Native FVP-RTP Ligand to RdRp

The structure of the replicating polymerase complex of SARS-CoV-2 with FVP-RTP in the pre-catalytic state (7CTT) [81] was retrieved from the PDB database. The pocket containing FVP-RTP has a surface area of 967.55 Å2 and a volume of 998.91 Å3, which is thus a surface/volume ratio of 0.97. Its hydrophobicity is 0.58. In this pre-catalytic state, one FVP-RTP molecule is incorporated into the RNA primer strand and forms a base-stacking interaction with adenosine in this strand. The conformation of the amide group in FVP-RTP is stabilized by the intramolecular N–H⋯O hydrogen bond of 2.612 Å (closing 6-member ring), which causes it to resemble guanine and facilitate π⋯π stacking. FVP-RTP is also involved in strong hydrogen bonds (two N–H⋯O of 2.4 Å and 3.2 Å and one N–H⋯N of 2.74 Å) with the cytosine from the RNA template strand. Furthermore, FVP-RTP binds to at least nine RdRp residues and one cofactor (magnesium Mg2+ ion) via various non-covalent interactions, as shown in Table 6 and Table 7 and Figure 4.
Non-covalent interactions include one N–H⋯O of 2.99 Å and two O–H⋯O hydrogen bonds of 2.44 Å and 3.38 Å binding FVP-RTP to serine SER682, and plenty of O–H⋯O and N–H⋯O hydrogen bonds between ribosyl and phosphate of FVP-RTP, and the neighboring residues, as shown in Table 1. FVP-RTP is also coordinated with ARG555 via π-cation interaction, with LYS545, which accepts hydrogen bonds, and with ARG555, LYS798, LYS621, and ARG553 via salt bridges. The interaction of FVP-RTP with Mg2+ is weak because it is long-range (of 6.94 Å). Fluorine participates in two F⋯O contacts of 4.957 Å and 4.691 Å and supports a C–H⋯O hydrogen bond of 2.694 Å to maintain the specific conformation of RTP. Moreover, the F⋯N contact of 4.151 Å supports the π⋯π stacking. This all adheres to the typical pattern of fluoride’s role in drug structures, specifically its impact on conformation. The binding mode of the candidate ligands should be consistent with that described above.

Binding Mode of the Candidate Ligands to RdRp

  • Docking Results
Triphosphorylated forms of the candidate ligands (with the F replaced by I, Br, Cl, H, CH3, CF3, or CN at the C(6) position) were prepared and docked to the binding site in RdRp that had been previously prepared by correcting protonation and atomic hybridization. The docking results are summarized in Table 8, and the best poses that led to the stabilization of the complex with the highest binding/docking score are shown in Figure 5.
The ordering of the ligands by descending docking score is as follows:
CN > CF3 > CH3 > Cl > Br > I > H > F.
As shown in Table 6, the total protein–ligand binding energy increased by 20–30% compared to the actual ligand, FVR-RTP. The docking score and total protein–ligand binding energy are highest for the CN analogue, followed by CF3 and Cl. According to decreasing total binding energy, the ligands can be ordered as follows:
CN > CF3 > Cl > I > CH3 > Br > H > F.
The sum of the binding energies of the ligand to the cofactor, RNA Template, and RNA primer is the highest for the CF3 analogue, followed by CN and Cl. Protein–ligand hydrogen bonds are the strongest for the CF3 analogue, followed by CN and Br. The orderings according to different energy characteristics, as presented in Table 9, show similar trends, i.e., the substitution of F to CF3 and CN leads to the most significant changes.
The ordering of the ligands according to the decreasing binding affinity is as follows:
CF3 > CN > F > H > Cl ≅ Br > I > CH3
The binding affinity is strongly but non-linearly correlated with the relative reactivity power R+ and R, Figure 6.
The highest values of R+ and R correspond to the strongest binding affinity. Thus, the CF3 and CN analogues indeed seem the most promising.
  • In-Depth Analysis of the Binding Mode
The in-depth inspection of the binding mode reveal that all FVP-RTP analogues bind to the same set of the RdRp residues (76 in total), RNA primer, RNA template, and one cofactor (Mg2+ ion), although via various non-covalent interactions. In this sense, the binding mode remains consistent among all analogues, while there are differences in the strength and nature of the interactions.
The binding mode of each ligand can be treated as a specific kind of “binding fingerprint”. Global differences in the binding modes of specific ligands relative to the actual ligand FVP-RTP (the entire complex, only protein residues, only cofactors, RNA template, and RNA primer) can be compared using different mathematical metrics (Euclidean or Manhattan) or simple summations, as shown in Table 10 and Table 11. Manhattan distance, which measures distance by pairwise aggregating the absolute difference between each variable, seems the most intuitive for relating the binding mode differences. Euclidean distance, which uses the squared difference in each variable, seems less convenient as it over-optimizes the result. Simple summing (additive method) shows only the balance of contributions.
Although each ligand interacts with as many as 76 residues of RdRp, most of these interactions are of minor importance, as shown in Table 8 and Table 9 and Figure S1. The CF3 and CN derivatives show the most significant differences in their binding modes compared to FVP, as Table 6, Table 8 and Table 9 show. Furthermore, these differences arise from different factors. While CN substitution modifies bindings to the RdRp residues, primarily those near the active site, CF3 significantly modifies interactions with the co-factors, RNA template, and RNA primer. Moreover, the docking results suggest that among the candidate ligands, the CF3 analogue should bind most strongly to both the RdRp, as well as the cofactor, RNA template, and RNA primer, and thus may be a very good alternative to FVP. Although the CN derivative has the highest total binding affinity, as shown in Table 6, it binds relatively weakly to the RNA template due to the unfavorable conformation of its amide group (tilted relative to the plane of the pyridazine ring).
Detailed insight into the binding mode using Ligplot+ [82,83] reveals a set of hydrogen bonds binding the ligands with ARG555, SER682, LYS98, LYS621, and LYS545 and hydrophobic interactions binding the ligands with ASP760, THR687, ASP623, VAL557, ASP618, adenosine, and uracil, as shown in Figure S2, which further stabilizes the protein–ligand complex. It should be noted that Ligplot+ suggests a similar set of hydrophobic interactions for all considered ligands binding to the RdRp. However, ASP623 forms NH⋯O or OH⋯O hydrogen bonds only with the CF3 and CN analogues and is therefore not recognized as hydrophobic, as shown in Figure S2 (middle and right).
The binding energies of the FVP-RTP analogues to the RdRp residues, RNA primer, RNA template, and cofactor (7CTT [81] target) are summarized in Table 12, and a radar plot comparing these data is shown in Figure 7.
As shown in Table 12 and Figure 7, the replacement of the F at C(6) with I, Br, Cl, CF3, H, CH3, CF3, and CN leads primarily to the changes in the binding strength of the ligand to the RNA primer, LYS545, ARG553, ARG555, ASP618, LYS621, CYS622, ASP623, SER682, ASN691, SER759, ASP760, ASP761, and LYS798. Moreover, the interactions of all FVP-RTP analogues with the conserved residue ASP618 and critical residues SER759, ASP760, and ASP761 are primarily electrostatic and repulsive; the bindings to the polar and hydrophilic LYS798, LYS621, SER682, ASP623, and ASN691 are strong and attractive, while their bindings to the hydrophilic cluster formed by LYS545, ARG553, and ARG555 are strong and mainly ionic (salt bridges). The bindings to the Mg2+ and Zn2+ cofactors are relatively weak, and the latter is negligible. However, the allocation of the binding energy between the individual residues is not uniform, Figure 7. Moreover, some residues, such as LYS798, LYS621, and ASP618, are highly sensitive to the type of ligand.
Regardless of whether only residues near the active site or all of them are taken into account, the most significant changes occur when F is replaced by CF3 or CN, as shown in Table 6, Table 8, Table 9 and Table 10. The most promising CF3 analogue interacts with the conserved residue ASP618 (11.267 kcal/mol) and three critical residues neutral SER759 (−0.752 kcal/mol), polar ASP761 (4.363 kcal/mol), and ASP760 (5.699 kcal/mol); these interactions are primarily electrostatic and repulsive in nature. Its bindings with LYS798 (−21.657 kcal/mol), LYS621 (−31.682 kcal/mol), SER682 (−10.147 kcal/mol), ASP623 (−5.553 kcal/mol), and ASN691 (−2.873 kcal/mol) are significantly stronger than the other ligands. The binding of the CF3 analogue to the hydrophilic cluster residues LYS545 (−9.728 kcal/mol), ARG553 (−16.868 kcal/mol), and ARG555 (−32.469 kcal/mol) is very strong, but its binding to the Mg2+ cofactor is relatively weak: only of −6.464 kcal/mol. Moreover, its binding to the RNA template and primer is the strongest among the candidate ligands (−48.652 and −16.309 kcal/mol for the RNA primer and template, respectively). The CN analogue binds to the RNA template and RNA primer less strongly than the CF3 analogue, but its binding to LYS621 and LYS798 is much stronger. The cosine distance in Table 10 shows the relatively low similarity of the CF3 and CN analogues to the FVP-RTP in terms of binding mode when considering the active site residues, RNA primer, RNA template, and cofactor and even slightly smaller when considering only active site residues. This effect is noticeable in all analogues, but it is noteworthy only in the CN and CF3 analogues. Importantly, the discrepancies in molecular fingerprints (Tanimoto distances, Figure 1) do not correspond to the discrepancies in the binding modes. Thus, the alteration in the binding mode resulting from a mere substitution of a group is complex and multifactorial.
  • Binding Mode Visualization
The sign[λ2(r)]ρ(r) surface mapped to the reduced density gradient, RDG(r), isosurface in the red–green–blue scheme for the RdRp-ligand complex visualizes the binding mode and reveals the nature of the non-covalent interactions between the RdRp and the ligands, as can be seen in Figure 8.
The overlap of the isosurfaces of RDG with sign(λ2BCP mapped over the surface shows the differences in binding mode for the candidate and actual ligand. Weak van der Waals-type interactions clearly dominate, which is indicated by the green color of the surface for each ligand. The large and nearly flat green area above the six-membered ring of the FVP ring proves π⋯π stacking between the ligands and adenosine. The N–H⋯O hydrogen bond linking the ligands to the RNA strand is clearly visible and depicted by a small light cyan disc-shaped area near =O and -NH2. The CF3 analogue binds more strongly to RdRp than the CN analogue, as confirmed by the significantly larger blue-green surfaces depicting attractive interactions, as shown in Figure 8.

2.3.2. Pre-Catalytic State—Productive Mode II (Binding to Uracil and Stacking to Guanosine)

Binding of the Native RVD Ligand to RdRp

The structure of the replicating polymerase complex of SARS-CoV-2 in the pre-catalytic state bound to RVD (7UO4) [80] was retrieved from the PDB database. The pocket containing RDV has a surface area of 1364.91 Å2, a volume of 1364.99 Å3, and a surface/volume ratio of 1.0. Its hydrophobicity factor is slightly higher and equal to 0.61. The backbone of 7UO4 [80], which holds the protein together and gives it its tertiary structure, differs from 7CTT [81] but only by 0.475% (backbone residues) and 1.817% (all residues). In this form, one RVD molecule incorporates into the RNA primer strand and forms a base-stacking interaction with the guanosine. It also participates in strong hydrogen bonds with the uracil moiety from the RNA template strand. Thus, this structure actually represents a slightly different binding mode than FVP-RTP in 7CTT [81].
Considering the possibility of alternative binding of FVP-RTP to the RNA strand, we used 7UO4 [80] to simulate this particular variant. However, the actual ligand, RVD, binds to a similar set of the RdRp residues (ASP623, ASN691, ARG555, LYS798, LYS621, ARG553, LYS551 ASP618, ASP760, and SER682) and one magnesium ion (cofactor), as shown in Table 13 and Table 14 and Figure 9.
The binding of RVD to the hydrophilic cluster residues LYS545, ARG553, ARG555, and LYS551 and the Mg2+ cofactor is very strong and involves π-cation (ARG555), salt bridge (LYS551), or electrostatic/charge–charge (ARG553 and Mg2+) interactions. The interaction between RVD and LYS798 is mixed in nature and can be described as a combination of salt bridges and electrostatic (charge–charge, attractive) interactions. RVD interactions with ASP618, ASP760, and ASP761 are primarily electrostatic but repulsive. Ligand bonds with guanosine from the RNA primer and SER682 are of the π⋯π type (hydrophobic). Conventional strong and attractive hydrogen bonds bind RDV to ASP623, ASN691, and uracil. In-depth analysis reveals that the actual ligand, RVD, is linked to the uracil in the RNA template via two hydrogen bonds: N–H⋯O of 2.878 Å and N–H⋯N of 3.10 Å, involving the amine group and nitrogen N(4), respectively. In addition, three N–H⋯O hydrogen bonds link –OH moieties from RDV to ASN691 (of 2.80 Å), LYS798 (of 2.89 Å), CYS622 (of 3.44 Å), and two O–H⋯O hydrogen bonds link RDV with ASP623 (of 3.21 Å). The atoms of the RDV heterocyclic ring and nitrogen atom from the –NH2 group participate in the π⋯π stacking with guanosine (of 3.99–5.0 Å). The RVD ligand binds to the NH3+ group of LYS798 using three oxygen atoms of its PO3 moiety.

Binding Mode of the Candidate Ligands to RdRp

  • Docking Results
The triphosphorylated forms of the candidate ligands in which the F at the C(6) position was replaced by I, Br, Cl, H, CH3, CF3, or CN were docked to the active site instead of the parental RDV. The docking results are summarized in Table 15, and the best poses are shown in Figure 10.
As shown in Table 15, the total protein–ligand binding energy increased by approximately 30% compared to the energy of the actual ligand, RDV. The highest docking score and total protein–ligand binding energy are observed for FVP-RTP and the CF3 and CN analogues are shown in Table 16. The ligands’ ordering according to the decreasing docking score is as follows:
F > CF3 > CN > Cl > Br > I > CH3 > H
According to the decreasing total energy of binding, the ligands can be ordered as follows:
CF3 > F > CN > CH3 > Cl > I > H > Br
The sum of the energies of binding ligand to cofactor, RNA Template, and RNA primer is the highest for the CH3 analogue, followed by F and CN. The protein–ligand hydrogen bonds are the strongest for FVP-RTP, followed by the H and Cl analogues.
  • In-Depth Analysis of the Binding Mode
Overall, the original FVP-RTP ligand, as well as its analogues, bind to the same set of RdRp residues, RNA primer, RNA template, and one cofactor (Mg2+ ion) but via different non-covalent interactions. In-depth analysis shows that the bindings of RVD and FVP-RTP to RdRp differ by as much as 20 residues (the Tanimoto distance 0.73). Differences in binding modes of specific ligands compared to the actual ligand RVD (the entire complex, only protein residues, only cofactors, RNA template, and RNA primer) calculated using different metrics (Euclidean or Manhattan) and additively are listed in Table 17 and Table 18.
In terms of overall interaction pattern, the CF3 analogue is closest to RDV, while the H and I analogues are the most different. However, the differences are mainly due to binding to the co-factors, RNA primer, and RNA template. The three analogues, F, CN, and CF3, appear to be the most preferred because they exhibit the strongest total protein–ligand binding, as well as a strong binding ability to both the RNA primer and RN template. The differences between the interactions within the entire complex and active site residues, as shown in Table 17 and Table 18, are subtle, so only residues near the ligand play a key role. Although each ligand interacts with as many as 83 residues of RdRp, most of them are of minor importance, Figure S3.
Docked ligands interact mainly with the active site residues via hydrogen bonding and hydrophobic interactions. Detailed insight into the binding mode using Ligplot+ [82,83] confirms the presence of the hydrogen bonds and reveals the hydrophobic interactions binding the ligands and ASP760, THR687, ASP623, VAL557, ASP618, adenosine, and uracil, as shown in Figure S4. Note that Ligplot+ suggests a similar set of hydrophobic interactions for the ligands binding to the same protein. However, only FVP-RTP forms an NH⋯O hydrogen bond with ASP623 and ASP760. In the case of the CF3 analogue, the hydrophobic contacts replace the hydrogen bonds. The CF3 analogue has a high docking score but relatively low binding affinity, so it may seem most convenient when it is desired to exclude binding to uracil.
The binding energies of FVR-RTP analogues with the RdRp residues, RNA primer, RNA template, and cofactor are summarized in Table 19, and a radar plot comparing these data is shown in Figure 11.
As follows from Table 19 and Figure 11, the replacement of the F at C(6) with I, Br, Cl, CF3, H, CH3, CF3, and CN leads primarily to the changes in the binding strength of the ligand to the RNA primer, Mg2+, RNA template, and residues: LYS545, LYS551, ARG553, ARG555, ASP618, LYS621, CYS622, ASP623, SER682, THR687, ASP760, ASP761, and LYS798. The nature of these interactions is the same as previously described for 7CTT [81], but their strengths are different, e.g., the bond with Mg2+ is very strong. Therefore, the distribution of the binding energy between individual residues is relatively uniform, Figure 11. Moreover, some residues, such as ARG553, ARG555, LYS621, and THR687, are highly sensitive to the type of the ligand. The high similarity of the binding modes in candidate ligands to RDV, greatest for the CF3 analogue, is confirmed by the high cosine distance values.
  • Binding Mode Visualization
The sign[λ2(r)]ρ(r) surface mapped to the RDG(r) isosurface in the red–green–blue scheme visualizes the binding mode for the RdRp-ligand complex and reveals the nature of the non-covalent interactions between RdRp and candidate ligands, shown in Figure 12.
Weak van der Waals-type interactions clearly dominate, as indicated by the green color of the surface. The N–H⋯O hydrogen bond linking the ligands to the RNA strand is clearly visible and depicted by a small light cyan disc-shaped region near =O and -NH2. The large and nearly flat green surfaces above the six-membered ring of the FVP ring prove π⋯π stacking between the ligands and guanosine. The differences between the modes of binding of the CF3 and CN analogues are small, as shown in Figure 12.
However, the differences between the binding modes in both pre-catalytic states (binding to cytosine and stacking to adenosine vs. binding to uracil and stacking to guanosine) for the same ligands are significant as the cosine distance ranges from 0.550 to 0.610. The most significant difference concerns the binding of ligand to Mg2+, which is six times stronger. Three residues, ARG555, ARG553, and ASP618, are important, but their ligand binding is only twice as strong. Ligand binding to the RNA template is only slightly stronger. The docking results suggest that among the analogues studied, the CF3 derivative should bind most strongly to RdRp, the cofactor Mg2+, the RNA template, and the RNA primer and may therefore be a very good alternative to FVP-RTP and RVD. The key factor for improving FVP effectiveness seems to be increasing the binding strength of the ligand to the RNA template.

2.3.3. Pre-Catalytic State—Non-Productive Mode

The structure of the complex of the SARS-CoV-2 replicating polymerase bound to FVP-RTP in the pre-catalytic and non-productive state (7AAP) [84] was retrieved from the PDB database. The pocket containing FVP-RTP has a surface area of 1120.59 Å2, a volume of 1205.75 Å3, and therefore a surface/volume ratio of 0.93. Its hydrophobicity is 0.61. The 7AAP [84] backbone differs from 7UO4 [80] by 0.737 (backbone residues) and 1.638 (all residues), which is almost twice as much as 7CTT [81]. The FVP-RTP molecule in 7AAP [84] weakly interacts with ASP618 (5.238 kcal/mol), SER759 (−0.79 kcal/mol), ASP760 (1.779 kcal/mol), and ASP761 (6.171 kcal/mol) but binds strongly to LYS545 (−9.426 kcal/mol), ARG553 (−4.621 kcal/mol), and ARG555 (−8.070 kcal/mol). In the non-productive mode, the phosphoribosyl conformation is different than in the productive mode and ensures a very strong binding of FVP-RTP to two Mg2+ ions (−28.57 and −44.10 kcal/mol).
FVP-RTP interactions with the hydrophilic cluster residues LYS545, ARG553, ARG555, and LYS551, as well as the Mg2+ cofactor, are very strong and include π-cation (ARG555), salt bridge (LYS551) or electrostatic/charge–charge interactions (ARG553, Mg2+). The bonding between FVP-RTP and LYS798 exhibits a mixed nature and can be described as a combination of salt bridges and electrostatic (charge–charge, attractive) interactions. The FVP-RTP interactions with ASP618, ASP760, and ASP761 are mainly electrostatic but repulsive. The interactions of the ligands with guanosine from RNA primer and SER682 are of the π⋯π type (hydrophobic). Conventional strong and attractive hydrogen bonds bind RDV to ASP623, ASN691, and uracil. In-depth analysis reveals that the actual ligand, FVP-RTP, is linked to uracil in the RNA template via two hydrogen bonds, N–H⋯O of 2.878 Å and N–H⋯N of 3.10 Å, involving the amine group and nitrogen N(4), respectively, Figure 13.

Docking Results

The docking results of FVP-RTP and its seven analogues are summarized in Table 20. Negative values of the docking scores suggest generally favorable docking of the ligands at the binding site.
As shown in Table 20, the docking score values are the highest for the CF3 analogue, followed by CN and I.
The ordering of the ligands in descending order of the docking scores is as follows:
CF3 > CN > I > Br > F > H > Cl > CH3
According to the decreasing total energy of binding, the ligands can be ordered as follows:
CF3 > Br > Cl > I > CH3 > CN > H > F
The total protein–ligand binding energy for the CF3 analogue was increased by 20% compared to that of the actual ligand, FVR-RTP. The sum of the binding energies of the ligand to the cofactor, RNA Template, and RNA primer is highest for the CH3 analogue, followed by CF3 and Br (this order is mainly due to the steric effect). Hydrogen bonds between protein and ligand are strongest for the CF3 analogue, followed by CH3 and H. The CF3 analogue’s highest binding affinity suggests it is the most promising ligand. On the contrary, the lowest binding affinity of the CN analogue suggests that it may offer a non-productive mode exclusion.

In-Depth Analysis of the Binding Mode

Differences in the binding modes of specific ligands compared to the actual ligand, RVD, calculated using different metrics (Euclidean or Manhattan) and additively are summarized in Table 21 and Table 22. Calculations were performed for three variants: the entire complex; protein residues; cofactors, RNA template and RNA primer.
The binding of the cofactors, RNA primer, and RNA template is the key factor that determines changes in the entire complex. The distinctions among the interactions throughout the entire complex and the active site residues are subtle, validating that only specific residues in close proximity to the ligand are crucial. In terms of the total interaction pattern, the CH3 analogue is the closest to the actual ligand FVR-RTP, while the CF3 and I analogues are the most different, as can be seen in Table 21 and Table 22. However, the differences in the strength of the residues’ interactions with individual ligands in 7AAP [84] are significantly greater than those observed in 7CTT [81] (especially for CF3 and I analogues). Overall, of the analogues tested, the CF3 derivative should bind most strongly to both RdRp and the cofactor, RNA template, and RNA primer, while the CN analogue provided the weakest binding.
The docked ligands interact mainly with the active site residues via hydrogen bonding and hydrophobic interactions. A thorough examination of binding modes using Ligplot+ uncovers the hydrophobic interactions that bind the ligands with ASP760, THR687, SER682, adenosine, and uracil, as shown in Figure S5.
The CF3 analogue additionally interacts with Asp618, while CN interacts with Asp555 and Asp623. Moreover, the CF3 derivative forms NH⋯O, NH⋯N, and NH⋯F hydrogen bonds with LYS545, while the CN analogue forms only one NH⋯N hydrogen bond. These bonds are absent in the structure of the 7AAP [84] complex. Thus, the binding mode in non-productive 7AAP is heavily influenced by the substituent.
The binding energies of the FVR-RTP analogues to the RdRp residues, RNA primer, RNA template, and cofactor are summarized in Table 23, and the best poses are shown in Figure 14.
A radar plot illustrating the discrepancies in binding modes is shown in Figure 15.
The largest variations in binding strength caused by the substituent effect are observed in the bonds connecting the ligand and the cofactor. Most residues, except LYS545 and ARG555, are insensitive to the type of ligand. The strong resemblance of binding modes in the complexes of RdRp with candidate ligands and the complex with FVP-RTP is reaffirmed by the elevated cosine distance values. The highest cosine distance value for the CF3 analogue revealed its high similarity to FVP-RTP.

Binding Mode Visualization

The sign[λ2(r)]ρ(r) surface mapped to the RDG(r) isosurface in the red–green–blue scheme for RdRp-ligand complex visualizes the binding mode and reveals the nature of the non-covalent interactions between the RdRp and the candidate ligands, as shown in Figure 16.
There are extensive green areas indicating the presence of weak van der Waals interactions, especially in the case of the CN derivative. The N–H⋯O hydrogen bond linking the ligands to cytosine is clearly visible and depicted as a small light cyan disc-shaped area near =O and -NH2. Both its shape and color reveal that the CF3 analogue binds more strongly to cytosine than the CN analogue. The large and nearly flat green area above the six-membered ring of the pyridazine ring proves π⋯π stacking between the ligands and guanosine. The CN analogue shows larger surfaces representing the interactions of phosphate groups with RdRp. However, the mostly green color indicates that these interactions are weaker than those observed for CF3.
The binding modes in 7AAP [84] and 7OU4 [80] differ in the binding strength within the complex components, i.e., the RNA primer, LYS621, LYS798, and Mg2+. The most significant differences between the non-productive (7AAP [84]) and productive modes (7CTT [81] and 7OU4 [80]) concern the binding strength of LYS798 and Mg2+ (three magnesium ions involved) with the ligand, which are six times weaker and very strong, respectively. The inorganic Mg2+ cofactor, which binds strongly to each ligand in the non-productive mode, changes the structure and chemical potential of the active site. The LYS798 residue, which stabilizes the core of the RdRp, is weakly bound by the ligand in this mode. The CF3 analogue shows the strongest binding to LYS798 among the FVP-RTP analogues and forms a strong binding to the cofactor. Moreover, the binding of the ligand to the RNA primer in the productive state is weaker, but to the RNA template, it is stronger. Importantly, the binding affinity in the non-productive mode is not correlated with the R+ and R-values.
The differences between the binding modes in productive and non-productive states can be compared using the radar plots in Figure 7, Figure 12 and Figure 15. In the productive modes, the bindings between the ligands and LYS621, LYS798, ARG555, and ARG553 are stronger than in the non-productive mode. Furthermore, the allocation of interaction energy between the active-site residues is more evenly distributed in the non-productive mode (i.e., the interactions are less directional).

2.3.4. Active State

The development of nucleotide-based medicines for combating COVID-19 relies on the understanding of the structure in the active state. The structure of the replicating polymerase complex of SARS-CoV-2 in the active state bound to FVP-RMP (7DFG) [85] was retrieved from the PDB database. The pocket containing FVP-RMP incorporated into the RNA strand has a volume of 2765 Å3. The two most promising ligands (the CF3 or CN analogues of FVP-RMP) were prepared, incorporated into the RNA strand, and docked.

Docking Results

The docking results are summarized in Table 24, and a radial plot comparing these data is shown in Figure 17.

In-Depth Analysis of the Binding Mode

The docking score is higher for the CN analogue than for the CF3 analogue, but the total energy of protein–ligand binding and hydrogen bond strength are higher for CF3 than for CN, as shown in Table 25.
Furthermore, the CF3 analogue has a much higher (−75.135 kcal/mol), while the CN analogue has a lower (−69.515 kcal/mol) binding affinity than the actual ligand (−72.457 kJ/mol). The radial plot presented in Figure 17 demonstrates uniformity and symmetry similar to that observed in the pre-catalytic state (target from 7CTT [81]) that was previously described. Substitution leads to significant differences in binding strengths of ARG836, ARG858, ARG513, ARG555, and LYS545, with ARG836 being most strongly affected. However, determining which ligand has a binding mode closest to FVP-RMP from examining just the plot shapes poses a challenge. In this task, the cosine distance is helpful. The cosine distance between the energetic profiles of binding of both ligands is high (0.996), which confirms their high similarity and binding efficiency. Moreover, the cosine distance between the energetic profiles of binding for FVP-RMP and the CN analogue incorporated in the RNA strand is slightly higher than that of FVP-RMP and the CF3 analogue (0.993 vs. 0.991, respectively). Thus, in terms of the overall interaction pattern, the CN analogue is closer to FVP-RMP than the CF3 analogue. The distinctions between the two ligands primarily stem from their bindings to the co-factors and RNA template. There is relatively little variation in the strength of the bindings within the active site residues. The best poses with the greatest binding/docking scores, illustrated in Figure 18, appear to be insignificantly distinct.
The disparities in the binding modes of particular ligands compared to the actual ligand, ascertained using various metrics (Euclidean and Manhattan) and additively, are summarized in Table 26 and Table 27. Three variants, the entire complex, protein residues, and cofactors and the RNA template, were taken into account.
Detailed insight into the structures of protein–ligand complexes using Ligplot+ reveals the hydrophobic interactions between the candidate ligands incorporated into the RNA strand and RdRp residues, Figure S6. As can be seen from Figure S6, all ligands interact with the residues ASP623, CYS813, ALA840, SER861, and MET855. Additionally, FVP-RTP interacts hydrophobically with ASP760 and ASP865; the CF3 analogue interacts with ASP691, LYS849, and GLU857; and CN interacts with LEU862, LYS849, and GLU857. Furthermore, both CN and CF3 analogues form OH⋯O hydrogen bonds with ASP865, while ASP760 forms metallic bonds with Mg2+ instead of the ligand. These bonds are not present in the current active state complex. Thus, the nature of the binding mode in the active state is strongly modulated by the substituent. Additionally, the hydrogen bond NH⋯F of 3.299Å between Lys545 and the CF3 ligand results in the RNA strand with the CF3 analogue being held more securely in the hydrophilic pocket. The CN analogue does not form an extra bond. Nevertheless, thanks to the -CN group, it facilitates water binding.

Binding Mode Visualization

The mapping of sign[λ2(r)]ρ(r) onto the RDG(r) isosurface demonstrates the significant similarity between the binding modes of the CF3 and CN analogues, as shown in Figure 19.
Green areas in Figure 19 indicate weak non-covalent interactions. The green, flat, and almost parallel surfaces in the single-stranded RNA template reveal numerous instances of hydrophobic π⋯π stacking between aromatic rings.

2.3.5. Allosteric Effect

A possible alternative mechanism of FVP action, which may explain the scattering of the results of clinical trials [86,87] or the synergistic effect observed in combined treatment against SARS-CoV-2 [88], can be associated with the allosteric effect [64]. The cryo-electron microscopy structure of SARS-CoV-2 virus full-length nsp12 in complex with cofactors nsp7 and nsp8 (6M71) [74] was retrieved from the PDB database. The allosteric sites were detected with Allosite-Pro [89] and PASSer [90] using a scheme [64,91,92].
Near the active site (2943 Å3 in volume), there are two allosteric pockets on either side of the RNA strand. One of them is composed of VAL557, ARG555, LYS545, THR680, SER682, THR556, ASP623, THR687, SER759, ASN691, ALA688, ASP760, and ARG553.

Docking Results

Ligands in which the F at the C(6) position was replaced by I, Br, Cl, H, CH3, CF3, or CN were docked to this site. The docking results are summarized in Table 28, and the best poses are shown in Figure 20.
The ligands’ ordering according to the decreasing docking score is as follows:
F > I > CN > Br > CF3 > CH3 > Cl > H
According to the decreasing total energy of binding, the ligands can be ordered as follows:
F > I > CN > Br > CF3 > H > CH3 > Cl
The ligands’ ordering according to the decreasing binding affinity is as follows:
CF3 > Br > I > F > Cl > H > CN > CH3
The CF3 analogue binds directly to the RdRp, much weaker than the F or CN ones. This direct binding of the CF3 analogue is also much weaker than when it is incorporated into the RNA strand.

In-Depth Analysis of the Binding Mode

Comparison of the cosine distance between the binding energy profiles of the potential ligands indicates that the CF3 analogue deviates most from others in terms of overall interaction pattern.
The conformation of FVP-RTP docked directly into the allosteric pocket is very similar to its conformation in the RNA template, Figure 18.
Docking of the FVP-RTP analogues directly into the allosteric pocket leads to their very strong bindings to hydrophilic cluster residues LYS545, ARG555, and ARG553, but much weaker to conserved residues ASP623 and ASP618 and critical catalytic residues SER682, SER759, ASP760 and ASP761, as shown in Table 29 and Figure 21.
Changing the ligand has no effect on binding to LYS798 and ASP761, while bindings to LYS545, SER682, and ASN691 are strongly modulated by the ligand, as shown in Figure 21. Surprisingly, despite the high binding affinity for the CF3 derivative, the cosine distance to the F analogue is the largest, which is mainly due to the weak binding to SER682.
Detailed insight into the structures of protein–ligand complexes using Ligplot+ reveals the hydrophobic interactions between the candidate ligands incorporated into the RNA strand and RdRp residues, Figure S7. As can be seen from Figure S7, all ligands interact hydrophobically with the residues ARG555, VAL557, and ALA688. Additionally, FVP-RTP interacts hydrophobically with CYS622, ASP623, ASP691, and ASP760, while the CN analogue interacts with SER682, CYS622, and THR556. Furthermore, FVP-RTP forms OH⋯O hydrogen bonds with SER759 and NH⋯O with ASN691, ASP623, TYR556, ARG553, and LYS545. The CF3 analogue forms OH⋯O hydrogen bonds with SER759 and THR676, and NH⋯O bonds with Arg553 and LYS545, while the CN analogue forms OH⋯O bonds with SER759 and NH⋯O bonds with ASP760, ASP623, THR687, ARG555, and LYS545. Thus, their binding modes differ significantly.

Binding Mode Visualization

The mapping of sign[λ2(r)]ρ(r) onto the RDG(r) isosurface in Figure 22 demonstrates above-mentioned differences between the binding modes of the CF3 and CN analogues.
Green areas in Figure 22 indicate non-covalent interactions. The green and flat regions reveal hydrophobic π⋯π stacking, while cyan discs indicate hydrogen bonds.

2.4. Comparison of the Productive, Non-Productive, Active State and Allosteric Site Binding Modes

Upon comparison of all protein–ligand complexes described above, it is apparent that the main channel proteins (LYS545, ARG553, and ARG555) are relatively loosely bound to the ligand in the non-productive and active states. Nevertheless, in the pre-catalytic states and in instances of direct binding, the binding is significantly stronger. The ligand binding to LYS798 firmly stabilizes the RdRp core in pre-catalytic states, but the binding is feeble in non-productive and direct binding states.
Radar plots, as shown in Figure 23, aid in quick screening, and highly bonded components are easily identified with pop-up peaks.
The CN analogue exhibits lower selectivity, as it promotes comparable binding modes in pre-catalytic state I and II. Conversely, the CF3 analogue primarily reinforces pre-catalytic state I. The latter, combined with a better ADME profile, increases CF3 attractiveness in terms of selective anti-virial activity. Strong direct binding of FVP-RTP to the active site residues suggests a possible alternative mechanism of FVP action, which may explain the scattering of the results of clinical trials [25,63] or the synergistic effect observed in combined treatment against SARS-CoV-2 [64]. The replacement of F with CF3 may help to overcome this limitation.

3. Methods

3.1. ADME and Drug-Likeness Prediction

In silico ADME drug-likeness evaluation was performed using the SwissADME tool developed by the Swiss Institute of Bioinformatics [66]. Drug-likeness was tested according to Lipinski, Veber, and Egan’s rules of 5 (RO5). The Abbot Bioavailability scores were computed to predict the probability of a compound having at least 10% oral bioavailability. Lipophilicity was predicted with iLOGP, XLOGP3, WLOGP, MLOGP, and SILICOS-IT models, from which a consensus log Po/w was determined [66]. The solubility (log S) of the ligands was predicted using SILICOS-IT [66]. The mutagenicity and carcinogenicity were predicted using models implemented in ADMET2.0 [71]. Synthetic accessibility (SA) was predicted based on the assumption that the frequency of molecular fragments in “really” obtainable molecules correlates with the ease of synthesis [66].

3.2. Density Functional Theory

Quantum chemical calculations were carried out within the Density Functional Theory (DFT) approach rooted in the Kohn–Sham [93] theorem, generalized by Levy [94]. Becke’s hybrid exchange-correlation functional, B3LYP [95,96]; Møller-Plesset, MP2 [97]; Minnesota hybrid meta M062X a high-nonlocality exchange-correlation functional with double the amount of nonlocal exchange (2X) [98], and the all-electron split-valence basis set 6-311G(d,p) were used to optimize the geometry of the ligands and their triphosphate forms. Geometry optimization was performed in the internal coordinates system. Following geometry optimization, a Hessian matrix (a matrix of second derivatives of the energy with respect to atomic coordinates) was calculated for all optimized ligand structures to verify the absence of imaginary frequencies (i.e., true stationary points). The polarizable continuum model (PCM) [99] was used to model solvation effects, required for the estimation of the reactivity indices in solution. The calculations were performed using Gaussian 16 rev. C01 [100].
The theoretical reactivity indices defined by Par and Pearson [101], Gazquez [102], and Chattaraj [103] are as follows: the absolute electronegativity [χ = −(ELUMO + EHOMO)/2; eV]; absolute hardness [η = (ELUMO − EHOMO)/2; eV)]; electrophilicity index (reactivity) [ω = χ2/2η; eV]; softness [S = 1/η; 1/eV]; electro-donating power [ω; eV]; electro-accepting power [ω+; eV]; net electrophilicity [Δω = ω+ + ω; eV]; and maximum number of electrons transferred in a chemical reaction [ΔNmax]. These indices were evaluated at MP2/6-311G(d,p) and M062X/6-311G(d,p), and describe the effect of donating (HOMO) and accepting (LUMO) electrons in the frozen core approximation (the Koopmans theorem [104]).

3.3. Quantum Theory of Atoms in Molecules

Theoretical analysis of the topology of intermolecular interactions was performed within Bader’s Quantum Theory of Atoms in Molecules (QTAIM) [105] approach using Gaussian 16 rev. C01 [100]. The DFT wavefunction was calculated at the M062X/6-311+G(d,p) level of the theory. The so-called topological descriptors calculated at the Bond Critical Points (BCPs) include the electron density at BCP, ρ(rBCP); three eigenvalues of the principal components of Hessian matrix composed of second partial derivatives of the electron density describing curvature of the electron density according to the principal axes λ1, λ2, and λ3; and Laplacian (the sum of Hessian eigenvalues), Δρ. Some of these descriptors (λ2 and Δρ) provide information about the nature and strength of the interactions. In many molecular systems, the reduced density gradient (RDG), which is highly sensitive to the electron density distribution and its variations, helps to detect weak or specific non-covalent interactions that were uncertain or not revealed by classical QTAIM [105,106].
The reduced density gradient is defined as follows:
R D G ( r ) = | ρ ( r ) | 2 ( 3 π ) 1 / 3 ρ ( r ) 4 / 3
where ∇ρ(r) is the electron density gradient, and ρ(r) is electron density. The RDG(r) vs. sign(λ2)ρ(r) plot reveals characteristic spikes in the low-gradient and low-density regions as long as the non-covalent contacts are present in the structure. The nature of these contacts can be determined using the sign of λ2, which allows their classification as attractive (stabilizing; λ2 < 0) or repulsive (destabilizing; λ2 > 0). A spike in the low-gradient, low-density region at λ2 < 0 suggests a stabilizing interaction such as a hydrogen bonding, a smaller spike accompanied with only slightly negative λ2 indicates a weakly stabilizing interaction, and a spike associated with λ2 > 0 indicates the absence of non-covalent interactions. Thus, the reduced density gradient (RDG) isosurface analysis allows for the analysis and visualization of non-covalent repulsive and attractive interactions. The results of the RDG analysis were visualized using Visual Molecular Dynamics (VMDs) [107]. The isosurfaces have been colored using a blue–green–red scale, where blue indicates strong attraction, red indicates strong repulsive interactions, and green indicates weak van der Waals interactions.

3.4. Molecular Docking

Molecular docking is a widely used technique in drug development for the screening of novel therapeutic agents. It requires knowledge of the reliable 3D structures of the ligand and target. The structures of the candidate ligands were optimized at the M062X/6-311+G(d,p) level, while the target RdRp structures were downloaded from the Protein databank PDB database (http://www.rcsb.org/pdb, accessed on 30 September 2023). The structures of the replicating polymerase complex of SARS-CoV-2 with FVP-RTP (7CTT) [81], RVD (7UO4) [80], FVP-RMP(7AAP) [84], and FVP-RTP (7DFG) [85] were retrieved from the PDB database. Prior to docking the candidate ligands, the native ligand that co-crystallized with the RNA was removed from the structure. The docking was performed using candidate ligands incorporated into the RNA strand. The cryo-electron microscopy structure of SARS-CoV-2 full-length nsp12 in complex with cofactors nsp7 and nsp8 (6M71) [74] was used exclusively for the studies of allosteric effect. The allosteric sites were detected with software dedicated for this purpose, i.e., Allosite-Pro [89] and PASSer [90], using schemes [64,91,92]. Allosite-Pro 2016 predicts allosteric sites using Support Vector Machines (SVMs), while PASSer 2021 uses eXtreme gradient boosting (XGBoost) and graph convolutional neural networks (GCNNs). However, both methods are based on a set of topological and physicochemical parameters. The results obtained from both methods were very similar. The active site differs only slightly in shape and size but not in position.
Conversion of the files with receptor and ligand structures to the .pdbqt format was carried out using MGLTools ver. 1.5.7. Molecular docking was performed using the automated docking tools AutoDock ver. 4.2.6 [108] and AutoDock Vina ver. 1.2.3 [109]. Both are designed to predict the binding of small molecules, such as substrates or drug candidates, to a receptor. AutoDock calculations were performed in three steps: (1) the preparation of coordinate files using AutoDockTools, (2) the pre-calculation of atomic affinities using AutoGrid, and (3) the docking of ligands using AutoDock. AutoDock AutoGrid was used to pre-calculate the grids for each type of atom in the ligand, as well as grids of electrostatic and desolvation potentials. AutoDock Vina automatically performs this task. The pre-calculated grid was treated as a template. Two techniques were used for docking: Template docking based on the native ligand and docking within a defined searched space around the active site. The search space was a grid box with step sizes covering 9–15 Å centered on the active site. The Autodock was applied for the docking of the ligand to a set of grids describing the target protein. Some specific sidechains (in the nearest neighborhood of the cavity) were treated as flexible. For docking, we applied the most computationally efficient method: a Lamarckian genetic algorithm (LGA) with only a slightly modified set of parameters.
To evaluate the quality of the docking process, we performed a redocking task. The actual ligand was removed from the parental structure and redocked in its own binding site. The redocking protocol was considered successful when the root-mean-square deviation (RMSD) of pose relative to its conformation in the parental structure did not exceed 3 Å. Then, in the same way, the new ligand was docked into the rigid protein structure. The number of rotatable bonds in the ligands was calculated, and the ligands were treated as flexible bodies. Due to the high flexibility of the ligands, as many as 50 conformational poses were searched for each ligand. AutoDock was performed many times to obtain several docked poses. After docking, the best poses leading to the stabilization of the protein–ligand complex with the highest docking score were selected and further investigated. The best pose was selected on the basis of the scoring function estimating the protein–ligand binding energy. In the final post-docking step, the residues (side chains) closest to the active ligand were minimized with respect to the pose found, the ligand was energy-minimized using standard potentials, and the positions of the protons in the entire complex were optimized. The docking results were further verified (redocking task) using a Genetic Evolutionary Method for molecular DOCKing (GEMDOCK) [110], which, similar to AutoDock, uses the genetic algorithm but a different evolution operator (limited Gaussian and Cauchy mutations and introduced rotamer-based mutations). The docking score is defined as a combination of the pairwise inter- and intramolecular energy with penalty term, the energy of binding of pharmacophore summarizing (electrostatic and hydrophilic terms) for all hot-spot atoms and penalty function which limits non-active conformations. Further post-docking exploration of protein–ligand non-covalent and hydrophobic interactions was performed using PoseEdit [111], Discovery Studio Visualizer [112], LigPlot+ [82,83], QTAIM/RDG, and GEMDOCK. Several techniques were employed because detecting weak, non-covalent interactions can be challenging. The GEMDOCK docking score was calculated as a sum of the intramolecular and intermolecular interactions corrected by a penalty value keeping the ligand inside the box (ligand–protein term), the pharmacophore-based interaction between the ligand and protein (binding site term) and penalty value for the electrostatic and hydrophilic preferences of the ligand (penalty term). The energy of the intramolecular interactions was described using the sum of the piecewise linear potential (PLP) term, which describes steric (van der Waals), hydrogen bonding interactions and the Coulomb electrostatic term. The energy of intramolecular interactions was estimated using the Gehlhaar model [113] with original parameterization, which employs piecewise linear potentials and a steric term. The binding affinity was estimated using a multiple linear regression model. The final 3D visualization of the binding modes was made using VMD [107].

3.5. Comparison of the Differences in the Binding Modes—Mathematical Metrics

A similarity of the molecules was calculated using the Tanimoto distance, defined as the ratio of the intersection of A and B to the union of A and B:
d ( A , B ) = A B A + B A B
where A and B are sets characterizing molecular structures.
A similarity of the binding modes was calculated using a few metrics defined below:
  • Euclidian distance d p , q = i ( p i q i ) 2 ;
  • Manhattan distance d p , q = i p i q i ;
  • the additive formula d p , q = i ( p i q i ) ;
  • cosine distance d p , q = i p i q i i p i i q i ;
  • Tanimoto distance d ( p , q ) = P Q P + Q P Q ;
where pi and qi are the binding interactions in each structure, P = {pi}, and Q = {qi}.
The plots were generated using OriginLab 2024 (OriginLab corporation, Northampton, MA, USA).

4. Conclusions

One of the primary objectives of small molecule drug discovery is to find innovative chemical compounds that possess desired properties. For this purpose, new global quantum chemical indices describing relative reactivity and new quantitative methods for the estimation and visualization of the differences in the binding modes of individual analogues with RdRp have been developed and applied, allowing the candidate ligands to be conveniently screened. The use of new relative reactivity indices can aid in the preliminary screening of candidate ligands. Euclidean, Manhattan, and cosine distances are helpful in quantifying the differences in binding modes between native and candidate ligands. Isosurfaces and radial plots can be used to visualize these differences, which is helpful in recognizing their location. This combined approach is particularly useful when comparing modifications or screening candidate drugs.
Quantitative analysis of the binding modes provided important clues to design high-affinity RdRp ligands. The proposed modification of the FVP structure seems to improve its binding ability to the SARS-CoV-2 RdRp and enhance productive binding mode.
In particular, two of the candidate ligands (the trifluoro- and cyano- derivatives) bind very strongly to both the RNA template and RdRp, so they may constitute a very good alternative to FVP. The CN analogue, like FVP, shows low and negative lipophilicity, while the CF3 analogue shows high and positive lipophilicity, which is more optimal for drugs. Replacing –F with –CF3, which has a strong electron-withdrawing nature, poor polarizability, and a broad hydrophobic domain, results in an increase in the hydrophobicity (lipophilicity) of the ligand and thus modifies the hydrophobic targets. Both CN and CF3 analogues seem promising because, as highly electrophilic reagents, they should have low substrate selectivity and inhibit a wide range of RdRps, not only SARS-CoV-2. Moreover, their relative electro-donating and electro-accepting powers are high; thus, both should act more efficiently than FVP-RTP. The CF3 analogue appears to be a better choice because it does not affect the preferred arrangement of the amide and triphosphate groups required for the effective binding to RNA primer and Mg2+, respectively. However, the CN analogue has a disadvantage because it is conformationally more labile. The CF3 analogue incorporated into the RNA strand binds to the RNA template and RNA primer more strongly than FVP, even in monophosphate form. Moreover, it shows the strongest binding to LYS798 among the FVP-RMP analogues and forms a strong binding to the cofactor. On the other hand, the CF3 analogue binds directly to RdRp much more weakly than the F or CN derivatives. This direct binding of the CF3 analogue is also much weaker than when it is incorporated into the RNA strand.
Upon comparison of all protein–ligand complexes described above, it is apparent that the main channel proteins (LYS545, ARG553, and ARG555) are relatively loosely bound to the ligand in the non-productive and active states. Nevertheless, in the pre-catalytic states and in instances of direct binding, the binding is significantly stronger. The ligand binding to LYS 798 firmly stabilizes the RdRp core in pre-catalytic states, but the binding is feeble in non-productive and direct binding states. The CN analogue exhibits lower selectivity, as it promotes comparable binding modes in pre-catalytic state I and II. Conversely, the CF3 analogue primarily reinforces pre-catalytic state I. The latter combined with a better ADME profile, increases CF3 attractiveness in terms of selective anti-virial activity. Replacing F with CF3 may help to overcome the limitations caused by the alternative mechanism of binding responsible for the scattering of the results of clinical trials or the synergistic effect observed in combined treatment against SARS-CoV-2. Increasing the number of acceptors seems to be a good step toward greater efficiency of FVP analogues. The docking results suggest that among the analogues studied, the CF3 derivative should bind most strongly to the cofactor Mg2+, the RNA template, and the RNA primer and therefore may be a very good alternative to FVP-RTP and RVD. Additionally, the incorporation of the CF3 group is expected to considerably reduce the likelihood of the formation of an alternate metabolite, in which the OH group is attached to the C(5) position due to the steric hindrance.
Our method for quantifying differences in binding mode holds promise for guiding future research on new anti-SARS-CoV-2 agents. If sufficient effectiveness in inhibiting viral replication in cell culture is established, they could be explored as potential drugs against COVID-19.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/molecules29020441/s1. Figure S1. The comparison of the total binding energy of FVR-RTP analogues to RdRp, RNA primer, RNA template, and cofactor (the RdRp target from 7CTT [81]); Figure S2. A LigPlot+ schematic 2D representation of the protein–ligand interactions (the RdRp target from 7CTT [81]). The hydrogen bonds are shown as green dashed lines, and residues involved in hydrophobic contacts are shown as red arcs and labeled; Figure S3. The comparison of the total binding affinity of FVR-RTP analogues with RdRp, RNA primer, RNA template, and cofactor (the RdRp target from 7UO4); Figure S4. A LigPlot+ schematic 2D representation of the protein–ligand interactions (the RdRp target from 7UO4 [80]). The hydrogen bonds are shown as green dashed lines, and residues involved in hydrophobic contacts are shown as red arcs and labeled; Figure S5. A LigPlot+ schematic 2D representation of the protein–ligand interactions (target from 7AAP [84]). The hydrogen bonds are shown as green dashed lines, and residues involved in hydrophobic contacts are shown as red arcs and labeled.; Figure S7. A LigPlot+ schematic 2D representation of the protein–ligand interactions (the RdRp target from 6M71 [74]). The hydrogen bonds are shown as green dashed lines, and residues involved in hydrophobic contacts are shown as red arcs and labeled.

Author Contributions

Conceptualization, J.N.L.; methodology, J.N.L.; investigation, quantum chemical calculations and modeling, J.N.L. and M.L.; visualizations, M.L.; writing—original draft preparation, J.N.L.; writing—review and editing, M.L. and J.N.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in article and Supplementary Materials.

Acknowledgments

Allotment of computer time from Wrocław Supercomputing and Networking Center (WCSS) is gratefully acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

AbbreviationMeaning
ADMEAbsorption, Distribution, Metabolism, and Excretion
BBBblood–brain barrier
CYPcytochrome P450
DFTDensity Functional Theory
DILIdrug-induced liver injury
DOIdegree of interaction
FVPfavipiravir
GIgastrointestinal absorption
HGPRThypoxanthinguanine phosphoribosyl transferase
H-HThepatotoxicity
MDMolecular Docking
MOLMolnupiravir
NMVrNirmatrelvir
PAINSPan Assay of Interference Structures
PCMpolarizable continuum model
P-gppermeability glycoprotein
QSPRQuantum Quantitative Structure–Property Relationship
QTAIMQuantum Theory of Atoms in Molecules
RDGReduced Density Gradient
RdRpRNA-dependent RNA polymerase
RDPribonucleoside-5′-diphosphate
RMPribonucleoside-5′-monophosphate
RTPribonucleoside-5′-triphosphate
RTVRitonavir
RNAribonucleic acid
RVDRemdesivir
RBVRibavirin
SASsynthetic accessibility score
ssRNAsingle-stranded RNA viruses
TPSAtopological polar surface area

References

  1. Perlman, S. Another Decade, Another Coronavirus. N. Engl. J. Med. 2020, 382, 760–762. [Google Scholar] [CrossRef]
  2. Zhang, T.; Wu, Q.; Zhang, Z. Probable Pangolin Origin of SARS-CoV-2 Associated with the COVID-19 Outbreak. Curr. Biol. 2020, 30, 1578. [Google Scholar] [CrossRef]
  3. WHO. Tracking SARS-CoV-2 Variants. Available online: https://www.who.int/activities/tracking-SARS-CoV-2-variants (accessed on 27 September 2023).
  4. Taubenberger, J.K.; Morens, D.M. 1918 Influenza: The Mother of All Pandemics. Emerg. Infect. Dis. 2006, 12, 15–22. [Google Scholar] [CrossRef] [PubMed]
  5. Hartman, A.L.; Zadeh, V.R.; Afowowe, T.O.; Abe, H.; Urata, S.; Yasuda, J. Potential and action mechanism of favipiravir as an antiviral against Junin virus. PLoS Pathog. 2022, 18, e1010689. [Google Scholar] [CrossRef]
  6. Konda, M.; Dodda, B.; Konala, V.M.; Naramala, S.; Adapa, S. Potential Zoonotic Origins of SARS-CoV-2 and Insights for Preventing Future Pandemics Through One Health Approach. Cureus 2020, 12, e8932. [Google Scholar] [CrossRef] [PubMed]
  7. Harris, A.; Sha, B.; Luo, M. Structural similarities between influenza virus matrix protein M1 and human immunodeficiency virus matrix and capsid proteins: An evolutionary link between negative-stranded RNA viruses and retroviruses. J. Gen. Virol. 1999, 80, 863–869. [Google Scholar] [CrossRef]
  8. Bar-On, Y.M.; Flamholz, A.; Phillips, R.; Milo, R. SARS-CoV-2 (COVID-19) by the numbers. eLife 2020, 9, e57309. [Google Scholar] [CrossRef] [PubMed]
  9. Kosik, I.; Yewdell, J.W. Influenza Hemagglutinin and Neuraminidase: Yin-Yang Proteins Coevolving to Thwart Immunity. Viruses 2019, 11, 346. [Google Scholar] [CrossRef]
  10. Letko, M.; Marzi, A.; Munster, V. Functional assessment of cell entry and receptor usage for SARS-CoV-2 and other lineage B betacoronaviruses. Nat. Microbiol. 2020, 5, 562–569. [Google Scholar] [CrossRef]
  11. Fantini, J.; Chahinian, H.; Yahi, N. Convergent Evolution Dynamics of SARS-CoV-2 and HIV Surface Envelope Glycoproteins Driven by Host Cell Surface Receptors and Lipid Rafts: Lessons for the Future. Int. J. Mol. Sci. 2023, 24, 1923. [Google Scholar] [CrossRef]
  12. Robson, F.; Khan, K.S.; Le, T.K.; Paris, C.; Demirbag, S.; Barfuss, P.; Rocchi, P.; Ng, W.-L. Coronavirus RNA Proofreading: Molecular Basis and Therapeutic Targeting. Mol. Cell 2020, 79, 710–727. [Google Scholar] [CrossRef]
  13. Kuhn, J.H.; Becker, S.; Ebihara, H.; Geisbert, T.W.; Johnson, K.M.; Kawaoka, Y.; Lipkin, W.I.; Negredo, A.I.; Netesov, S.V.; Nichol, S.T.; et al. Proposal for a revised taxonomy of the family Filoviridae: Classification, names of taxa and viruses, and virus abbreviations. Arch. Virol. 2010, 155, 2083–2103. [Google Scholar] [CrossRef] [PubMed]
  14. Cheung, C.L.; Rayner, J.M.; Smith, G.J.D.; Wang, P.; Naipospos, T.S.P.; Zhang, J.; Yuen, K.Y.; Webster, R.G.; Peiris, J.S.M.; Guan, Y.; et al. Distribution of Amantadine-Resistant H5N1 Avian Influenza Variants in Asia. J. Infect. Dis. 2006, 193, 1626–1629. [Google Scholar] [CrossRef]
  15. Deyde, V.M.; Xu, X.; Bright, R.A.; Shaw, M.; Smith, C.B.; Zhang, Y.; Shu, Y.; Gubareva, L.V.; Cox, N.J.; Klimov, A.I. Surveillance of Resistance to Adamantanes among Influenza A(H3N2) and A(H1N1) Viruses Isolated Worldwide. J. Infect. Dis. 2007, 196, 249–257. [Google Scholar] [CrossRef] [PubMed]
  16. Nelson, M.I.; Simonsen, L.; Viboud, C.; Miller, M.A.; Holmes, E.C. The origin and global emergence of adamantane resistant A/H3N2 influenza viruses. Virology 2009, 388, 270–278. [Google Scholar] [CrossRef] [PubMed]
  17. Hata, M.; Tanaka, S.; Yasui, Y.; Fujiwara, N.; Kobayashi, S.; Minagawa, H. Rapid Genotypic Assay for Detection of Oseltamivir-Resistant Influenza A (H1N1) Viruses. J. Clin. Microbiol. 2010, 48, 1983–1984. [Google Scholar] [CrossRef] [PubMed]
  18. Earhart, K.C.; Elsayed, N.M.; Saad, M.D.; Gubareva, L.V.; Nayel, A.; Deyde, V.M.; Abdelsattar, A.; Abdelghani, A.S.; Boynton, B.R.; Mansour, M.M.; et al. Oseltamivir resistance mutation N294S in human influenza A(H5N1) virus in Egypt. J. Infect. Public Health 2009, 2, 74–80. [Google Scholar] [CrossRef]
  19. Meijer, A.; Lackenby, A.; Hungnes, O.; Lina, B.; van der Werf, S.; Schweiger, B.; Opp, M.; Paget, J.; van de Kassteele, J.; Hay, A.; et al. Oseltamivir-Resistant Influenza Virus A (H1N1), Europe, 2007–2008 Season. Emerg. Infect. Dis. 2009, 15, 552–560. [Google Scholar] [CrossRef]
  20. Report on the Deliberation Results March 4, 2014 Evaluation and Licensing Division, Pharmaceutical and Food Safety Bureau Ministry of Health, Labour and Welfare, Pharmaceuticals and Medical Devices Agency. Available online: https://www.pmda.go.jp/files/000210319.pdf (accessed on 27 September 2023).
  21. Hayden, F.G.; Lenk, R.P.; Stonis, L.; Oldham-Creamer, C.; Kang, L.L.; Epstein, C. Favipiravir Treatment of Uncomplicated Influenza in Adults: Results of Two Phase 3, Randomized, Double-Blind, Placebo-Controlled Trials. J. Infect. Dis. 2022, 226, 1790–1799. [Google Scholar] [CrossRef]
  22. Takashita, E.; Ejima, M.; Ogawa, R.; Fujisaki, S.; Neumann, G.; Furuta, Y.; Kawaoka, Y.; Tashiro, M.; Odagiri, T. Antiviral susceptibility of influenza viruses isolated from patients pre- and post-administration of favipiravir. Antivir. Res. 2016, 132, 170–177. [Google Scholar] [CrossRef]
  23. Jordan, P.C.; Stevens, S.K.; Deval, J. Nucleosides for the treatment of respiratory RNA virus infections. Antivir. Chem. Chemother. 2018, 26, 2040206618764483. [Google Scholar] [CrossRef] [PubMed]
  24. Watanabe, M. The COVID-19 Pandemic in Japan. Surg. Today 2020, 50, 787–793. [Google Scholar] [CrossRef] [PubMed]
  25. Kumar, R.; Gupta, N.; Kodan, P.; Mittal, A.; Soneja, M.; Wig, N. Battling COVID-19: Using old weapons for a new enemy. Trop. Dis. Travel Med. Vaccines 2020, 6, 6. [Google Scholar] [CrossRef] [PubMed]
  26. Şimşek Yavuz, S.; ÜNal, S. Antiviral treatment of COVID-19. Turk. J. Med. Sci. 2020, 50, 611–619. [Google Scholar] [CrossRef]
  27. Delang, L.; Abdelnabi, R.; Neyts, J. Favipiravir as a potential countermeasure against neglected and emerging RNA viruses. Antivir. Res. 2018, 153, 85–94. [Google Scholar] [CrossRef]
  28. Furuta, Y.; Komeno, T.; Nakamura, T. Favipiravir (T-705), a broad spectrum inhibitor of viral RNA polymerase. Proc. Jpn. Acad. Ser. B 2017, 93, 449–463. [Google Scholar] [CrossRef]
  29. Bai, C.-Q.; Mu, J.-S.; Kargbo, D.; Song, Y.-B.; Niu, W.-K.; Nie, W.-M.; Kanu, A.; Liu, W.-W.; Wang, Y.-P.; Dafae, F.; et al. Clinical and Virological Characteristics of Ebola Virus Disease Patients Treated With Favipiravir (T-705)—Sierra Leone, 2014. Clin. Infect. Dis. 2016, 63, 1288–1294. [Google Scholar] [CrossRef]
  30. Rosenke, K.; Feldmann, H.; Westover, J.B.; Hanley, P.W.; Martellaro, C.; Feldmann, F.; Saturday, G.; Lovaglio, J.; Scott, D.P.; Furuta, Y.; et al. Use of Favipiravir to Treat Lassa Virus Infection in Macaques. Emerg. Infect. Dis. 2018, 24, 1696–1699. [Google Scholar] [CrossRef]
  31. Escribano-Romero, E.; Jiménez de Oya, N.; Domingo, E.; Saiz, J.C. Extinction of West Nile Virus by Favipiravir through Lethal Mutagenesis. Antimicrob. Agents Chemother. 2017, 61, 10–1128. [Google Scholar] [CrossRef]
  32. Marlin, R.; Desjardins, D.; Contreras, V.; Lingas, G.; Solas, C.; Roques, P.; Naninck, T.; Pascal, Q.; Behillil, S.; Maisonnasse, P.; et al. Antiviral efficacy of favipiravir against Zika and SARS-CoV-2 viruses in non-human primates. Nat. Commun. 2022, 13, 5108. [Google Scholar] [CrossRef]
  33. Williams, M.; Caroline, A.L.; Powell, D.S.; Bethel, L.M.; Oury, T.D.; Reed, D.S.; Hartman, A.L. Broad Spectrum Antiviral Activity of Favipiravir (T-705): Protection from Highly Lethal Inhalational Rift Valley Fever. PLoS Negl. Trop. Dis. 2014, 8, e2790. [Google Scholar] [CrossRef]
  34. Julander, J.G.; Shafer, K.; Smee, D.F.; Morrey, J.D.; Furuta, Y. Activity of T-705 in a Hamster Model of Yellow Fever Virus Infection in Comparison with That of a Chemically Related Compound, T-1106. Antimicrob. Agents Chemother. 2009, 53, 202–209. [Google Scholar] [CrossRef]
  35. Singh, S.K.; Oestereich, L.; Rieger, T.; Neumann, M.; Bernreuther, C.; Lehmann, M.; Krasemann, S.; Wurr, S.; Emmerich, P.; de Lamballerie, X.; et al. Evaluation of Antiviral Efficacy of Ribavirin, Arbidol, and T-705 (Favipiravir) in a Mouse Model for Crimean-Congo Hemorrhagic Fever. PLoS Negl. Trop. Dis. 2014, 8, e2804. [Google Scholar] [CrossRef]
  36. Bologheanu, R.; Schubert, L.; Thurnher, M.; Schiefer, J.; Santonja, I.; Holzmann, H.; Oesterreicher, Z.; Tobudic, S.; Winkler, S.; Faybik, P.; et al. Unexpected complete recovery of a patient with severe tick-borne encephalitis treated with favipiravir. Antivir. Res. 2020, 184, 104952. [Google Scholar] [CrossRef]
  37. Banyard, A.C.; Mansfield, K.L.; Wu, G.; Selden, D.; Thorne, L.; Birch, C.; Koraka, P.; Osterhaus, A.D.M.E.; Fooks, A.R. Re-evaluating the effect of Favipiravir treatment on rabies virus infection. Vaccine 2019, 37, 4686–4693. [Google Scholar] [CrossRef] [PubMed]
  38. Veliziotis, I.; Roman, A.; Martiny, D.; Schuldt, G.; Claus, M.; Dauby, N.; Van den Wijngaert, S.; Martin, C.; Nasreddine, R.; Perandones, C.; et al. Clinical Management of Argentine Hemorrhagic Fever using Ribavirin and Favipiravir, Belgium, 2020. Emerg. Infect. Dis. 2020, 26, 1562–1566. [Google Scholar] [CrossRef] [PubMed]
  39. Agrawal, U.; Raju, R.; Udwadia, Z.F. Favipiravir: A new and emerging antiviral option in COVID-19. Med. J. Armed Forces India 2020, 76, 370–376. [Google Scholar] [CrossRef] [PubMed]
  40. Joshi, S.; Parkar, J.; Ansari, A.; Vora, A.; Talwar, D.; Tiwaskar, M.; Patil, S.; Barkate, H. Role of favipiravir in the treatment of COVID-19. Int. J. Infect. Dis. 2021, 102, 501–508. [Google Scholar] [CrossRef]
  41. Mallhi, T.H.; Fathi, M.; Vakili, K.; Sayehmiri, F.; Mohamadkhani, A.; Hajiesmaeili, M.; Rezaei-Tavirani, M.; Eilami, O. The prognostic value of comorbidity for the severity of COVID-19: A systematic review and meta-analysis study. PLoS ONE 2021, 16, e0246190. [Google Scholar] [CrossRef]
  42. Wang, Y.; Li, P.; Rajpoot, S.; Saqib, U.; Yu, P.; Li, Y.; Li, Y.; Ma, Z.; Baig, M.S.; Pan, Q. Comparative assessment of favipiravir and remdesivir against human coronavirus NL63 in molecular docking and cell culture models. Sci. Rep. 2021, 11, 23465. [Google Scholar] [CrossRef]
  43. Gunaydin-Akyildiz, A.; Aksoy, N.; Boran, T.; Ilhan, E.N.; Ozhan, G. Favipiravir induces oxidative stress and genotoxicity in cardiac and skin cells. Toxicol. Lett. 2022, 371, 9–16. [Google Scholar] [CrossRef]
  44. Kumar, P.; Kulkarni, A.; Sharma, M.; Rao, P.N.; Reddy, D.N. Favipiravir-induced Liver Injury in Patients with Coronavirus Disease 2019. J. Clin. Transl. Hepatol. 2021, 9, 276. [Google Scholar] [CrossRef] [PubMed]
  45. Madelain, V.; Mentré, F.; Baize, S.; Anglaret, X.; Laouénan, C.; Oestereich, L.; Nguyen, T.H.T.; Malvy, D.; Piorkowski, G.; Graw, F.; et al. Modeling Favipiravir Antiviral Efficacy Against Emerging Viruses: From Animal Studies to Clinical Trials. CPT Pharmacomet. Syst. Pharmacol. 2020, 9, 258–271. [Google Scholar] [CrossRef]
  46. Tırmıkçıoğlu, Z. Favipiravir exposure and pregnancy outcome of COVID-19 patients. Eur. J. Obstet. Gynecol. Reprod. Biol. 2022, 268, 110–115. [Google Scholar] [CrossRef] [PubMed]
  47. Ertem, O.; Guner, O.; Incir, C.; Kalkan, S.; Gelal, A. The outcomes of favipiravir exposure in pregnancy: A case series. Arch. Gynecol. Obstet. 2022, 307, 1385–1395. [Google Scholar] [CrossRef] [PubMed]
  48. Nagata, T.; Lefor, A.K.; Hasegawa, M.; Ishii, M. Favipiravir: A New Medication for the Ebola Virus Disease Pandemic. Disaster Med. Public Health Prep. 2014, 9, 79–81. [Google Scholar] [CrossRef] [PubMed]
  49. Casalini, G.; Giacomelli, A.; Antinori, S. Liver tests abnormalities with licensed antiviral drugs for COVID-19: A narrative review. Expert Opin. Drug Saf. 2023, 21, 1483–1494. [Google Scholar] [CrossRef] [PubMed]
  50. Dufour, I.; Devresse, A.; Scohy, A.; Briquet, C.; Georgery, H.; Delaey, P.; Greef, J.D.; Goffin, E.; Labriola, L. Safety and efficiency of molnupiravir for COVID-19 patients with advanced chronic kidney disease. Kidney Res. Clin. Pract. 2023, 42, 275–278. [Google Scholar] [CrossRef]
  51. Zhou, S.; Hill, C.S.; Sarkar, S.; Tse, L.V.; Woodburn, B.M.D.; Schinazi, R.F.; Sheahan, T.P.; Baric, R.S.; Heise, M.T.; Swanstrom, R. β-d-N4-hydroxycytidine Inhibits SARS-CoV-2 Through Lethal Mutagenesis But Is Also Mutagenic To Mammalian Cells. J. Infect. Dis. 2021, 224, 415–419. [Google Scholar] [CrossRef]
  52. Smee, D.F.; Hurst, B.L.; Wong, M.-H.; Bailey, K.W.; Tarbet, E.B.; Morrey, J.D.; Furuta, Y. Effects of the Combination of Favipiravir (T-705) and Oseltamivir on Influenza A Virus Infections in Mice. Antimicrob. Agents Chemother. 2010, 54, 126–133. [Google Scholar] [CrossRef]
  53. Smee, D.F.; Tarbet, E.B.; Furuta, Y.; Morrey, J.D.; Barnard, D.L. Synergistic combinations of favipiravir and oseltamivir against wild-type pandemic and oseltamivir-resistant influenza A virus infections in mice. Future Virol. 2013, 8, 1085–1094. [Google Scholar] [CrossRef]
  54. Sarojvisut, P.; Apisarnthanarak, A.; Jantarathaneewat, K.; Sathitakorn, O.; Pienthong, T.; Mingmalairak, C.; Warren, D.K.; Weber, D.J. An Open Label Randomized Controlled Trial of Ivermectin Plus Favipiravir-Based Standard of Care versus Favipiravir-Based Standard of Care for Treatment of Moderate COVID-19 in Thailand. Infect. Chemother. 2023, 55, 50. [Google Scholar] [CrossRef]
  55. Oestereich, L.; Rieger, T.; Lüdtke, A.; Ruibal, P.; Wurr, S.; Pallasch, E.; Bockholt, S.; Krasemann, S.; Muñoz-Fontela, C.; Günther, S. Efficacy of Favipiravir Alone and in Combination With Ribavirin in a Lethal, Immunocompetent Mouse Model of Lassa Fever. J. Infect. Dis. 2016, 213, 934–938. [Google Scholar] [CrossRef] [PubMed]
  56. Zhao, H.; Zhu, Q.; Zhang, C.; Li, J.; Wei, M.; Qin, Y.; Chen, G.; Wang, K.; Yu, J.; Wu, Z.; et al. Tocilizumab combined with favipiravir in the treatment of COVID-19: A multicenter trial in a small sample size. Biomed. Pharmacother. 2021, 133, 110825. [Google Scholar] [CrossRef] [PubMed]
  57. Abdel-Halim, H.; Hajar, M.; Hasouneh, L.; Abdelmalek, S.M.A. Identification of Drug Combination Therapies for SARS-CoV-2: A Molecular Dynamics Simulations Approach. Drug Des. Dev. Ther. 2022, 16, 2995–3013. [Google Scholar] [CrossRef]
  58. Rubin, D.; Chan-Tack, K.; Farley, J.; Sherwat, A. FDA Approval of Remdesivir—A Step in the Right Direction. N. Engl. J. Med. 2020, 383, 2598–2600. [Google Scholar] [CrossRef] [PubMed]
  59. Gérard, A.O.; Laurain, A.; Fresse, A.; Parassol, N.; Muzzone, M.; Rocher, F.; Esnault, V.L.M.; Drici, M.D. Remdesivir and Acute Renal Failure: A Potential Safety Signal From Disproportionality Analysis of the WHO Safety Database. Clin. Pharmacol. Ther. 2021, 109, 1021–1024. [Google Scholar] [CrossRef]
  60. Stevens, L.J.; Pruijssers, A.J.; Lee, H.W.; Gordon, C.J.; Tchesnokov, E.P.; Gribble, J.; George, A.S.; Hughes, T.M.; Lu, X.; Li, J.; et al. Mutations in the SARS-CoV-2 RNA-dependent RNA polymerase confer resistance to remdesivir by distinct mechanisms. Sci. Transl. Med. 2022, 14, eabo0718. [Google Scholar] [CrossRef]
  61. Sangawa, H.; Komeno, T.; Nishikawa, H.; Yoshida, A.; Takahashi, K.; Nomura, N.; Furuta, Y. Mechanism of Action of T-705 Ribosyl Triphosphate against Influenza Virus RNA Polymerase. Antimicrob. Agents Chemother. 2013, 57, 5202–5208. [Google Scholar] [CrossRef]
  62. Baranovich, T.; Wong, S.-S.; Armstrong, J.; Marjuki, H.; Webby, R.J.; Webster, R.G.; Govorkova, E.A. T-705 (Favipiravir) Induces Lethal Mutagenesis in Influenza A H1N1 Viruses In Vitro. J. Virol. 2013, 87, 3741–3751. [Google Scholar] [CrossRef]
  63. Naesens, L.; Guddat, L.W.; Keough, D.T.; van Kuilenburg, A.B.P.; Meijer, J.; Vande Voorde, J.; Balzarini, J. Role of Human Hypoxanthine Guanine Phosphoribosyltransferase in Activation of the Antiviral Agent T-705 (Favipiravir). Mol. Pharmacol. 2013, 84, 615–629. [Google Scholar] [CrossRef] [PubMed]
  64. Latosińska, J.N.; Latosińska, M.; Seliger, J.; Žagar, V.; Apih, T.; Grieb, P. Elucidating the Role of Noncovalent Interactions in Favipiravir, a Drug Active against Various Human RNA Viruses; a 1H-14N NQDR/Periodic DFT/QTAIM/RDS/3D Hirshfeld Surfaces Combined Study. Molecules 2023, 28, 3308. [Google Scholar] [CrossRef] [PubMed]
  65. Etter, M.C. Encoding and decoding hydrogen-bond patterns of organic compounds. Acc. Chem. Res. 2002, 23, 120–126. [Google Scholar] [CrossRef]
  66. Daina, A.; Michielin, O.; Zoete, V. SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 2017, 7, 42717. [Google Scholar] [CrossRef]
  67. Jia, H.; Häring, A.P.; Berger, F.; Zhang, L.; Ritter, T. Trifluoromethyl Thianthrenium Triflate: A Readily Available Trifluoromethylating Reagent with Formal CF3+, CF3•, and CF3− Reactivity. J. Am. Chem. Soc. 2021, 143, 7623–7628. [Google Scholar] [CrossRef]
  68. Mucker, E.M.; Goff, A.J.; Shamblin, J.D.; Grosenbach, D.W.; Damon, I.K.; Mehal, J.M.; Holman, R.C.; Carroll, D.; Gallardo, N.; Olson, V.A.; et al. Efficacy of Tecovirimat (ST-246) in Nonhuman Primates Infected with Variola Virus (Smallpox). Antimicrob. Agents Chemother. 2013, 57, 6246–6253. [Google Scholar] [CrossRef]
  69. Wang, Z.; Yu, Z.; Kang, D.; Zhang, J.; Tian, Y.; Daelemans, D.; De Clercq, E.; Pannecouque, C.; Zhan, P.; Liu, X. Design, synthesis and biological evaluation of novel acetamide-substituted doravirine and its prodrugs as potent HIV-1 NNRTIs. Bioorg. Med. Chem. 2019, 27, 447–456. [Google Scholar] [CrossRef] [PubMed]
  70. Turner, S.R.; Strohbach, J.W.; Tommasi, R.A.; Aristoff, P.A.; Johnson, P.D.; Skulnick, H.I.; Dolak, L.A.; Seest, E.P.; Tomich, P.K.; Bohanon, M.J.; et al. Tipranavir (PNU-140690): A Potent, Orally Bioavailable Nonpeptidic HIV Protease Inhibitor of the 5,6-Dihydro-4-hydroxy-2-pyrone Sulfonamide Class. J. Med. Chem. 1998, 41, 3467–3476. [Google Scholar] [CrossRef]
  71. Xiong, G.; Wu, Z.; Yi, J.; Fu, L.; Yang, Z.; Hsieh, C.; Yin, M.; Zeng, X.; Wu, C.; Lu, A.; et al. ADMETlab 2.0: An integrated online platform for accurate and comprehensive predictions of ADMET properties. Nucleic Acids Res. 2021, 49, W5–W14. [Google Scholar] [CrossRef]
  72. Gao, Z.-G.; Jacobson, K.A. Emerging adenosine receptor agonists. Expert Opin. Emerg. Drugs 2007, 12, 479–492. [Google Scholar] [CrossRef]
  73. Xu, X. Molecular model of SARS coronavirus polymerase: Implications for biochemical functions and drug design. Nucleic Acids Res. 2003, 31, 7117–7130. [Google Scholar] [CrossRef] [PubMed]
  74. Gao, Y.; Yan, L.; Huang, Y.; Liu, F.; Zhao, Y.; Cao, L.; Wang, T.; Sun, Q.; Ming, Z.; Zhang, L.; et al. Structure of the RNA-dependent RNA polymerase from COVID-19 virus. Science 2020, 368, 779–782. [Google Scholar] [CrossRef] [PubMed]
  75. Tao, Y.; Farsetta, D.L.; Nibert, M.L.; Harrison, S.C. RNA Synthesis in a Cage—Structural Studies of Reovirus Polymerase λ3. Cell 2002, 111, 733–745. [Google Scholar] [CrossRef] [PubMed]
  76. Butcher, S.J.; Grimes, J.M.; Makeyev, E.V.; Bamford, D.H.; Stuart, D.I. A mechanism for initiating RNA-dependent RNA polymerization. Nature 2001, 410, 235–240. [Google Scholar] [CrossRef]
  77. Jacobo-Molina, A.; Ding, J.; Nanni, R.G.; Clark, A.D.; Lu, X.; Tantillo, C.; Williams, R.L.; Kamer, G.; Ferris, A.L.; Clark, P. Crystal structure of human immunodeficiency virus type 1 reverse transcriptase complexed with double-stranded DNA at 3.0 A resolution shows bent DNA. Proc. Natl. Acad. Sci. USA 1993, 90, 6320–6324. [Google Scholar] [CrossRef]
  78. Latosińska, J.N.; Latosińska, M.; Orzeszko, A.; Maurin, J.K. Synthesis and Crystal Structure of Adamantylated 4,5,6,7-Tetrahalogeno-1H-benzimidazoles Novel Multi-Target Ligands (Potential CK2, M2 and SARS-CoV-2 Inhibitors); X-ray/DFT/QTAIM/Hirshfeld Surfaces/Molecular Docking Study. Molecules 2022, 28, 147. [Google Scholar] [CrossRef]
  79. Latosińska, J.N.; Latosińska, M.; Maurin, J.K.; Orzeszko, A.; Kazimierczuk, Z. Quantum-Chemical Insight into Structure-Reactivity Relationship in 4,5,6,7-Tetrahalogeno-1H-benzimidazoles: A Combined X-ray, DSC, DFT/QTAIM, Hirshfeld Surface-Based, and Molecular Docking Approach. J. Phys. Chem. A 2014, 118, 2089–2106. [Google Scholar] [CrossRef]
  80. Malone, B.F.; Perry, J.K.; Olinares, P.D.B.; Lee, H.W.; Chen, J.; Appleby, T.C.; Feng, J.Y.; Bilello, J.P.; Ng, H.; Sotiris, J.; et al. Structural basis for substrate selection by the SARS-CoV-2 replicase. Nature 2023, 614, 781–787. [Google Scholar] [CrossRef]
  81. Peng, Q.; Peng, R.; Yuan, B.; Wang, M.; Zhao, J.; Fu, L.; Qi, J.; Shi, Y. Structural Basis of SARS-CoV-2 Polymerase Inhibition by Favipiravir. Innovation 2021, 2, 100080. [Google Scholar] [CrossRef]
  82. Wallace, A.C.; Laskowski, R.A.; Thornton, J.M. LIGPLOT: A program to generate schematic diagrams of protein-ligand interactions. Protein Eng. Des. Sel. 1995, 8, 127–134. [Google Scholar] [CrossRef]
  83. Laskowski, R.A.; Swindells, M.B. LigPlot+: Multiple Ligand–Protein Interaction Diagrams for Drug Discovery. J. Chem. Inf. Model. 2011, 51, 2778–2786. [Google Scholar] [CrossRef] [PubMed]
  84. Naydenova, K.; Muir, K.W.; Wu, L.-F.; Zhang, Z.; Coscia, F.; Peet, M.J.; Castro-Hartmann, P.; Qian, P.; Sader, K.; Dent, K.; et al. Structure of the SARS-CoV-2 RNA-dependent RNA polymerase in the presence of favipiravir-RTP. Proc. Natl. Acad. Sci. USA 2021, 118, e2021946118. [Google Scholar] [CrossRef] [PubMed]
  85. Li, Z.; Zhou, Z.; Yu, X. Structure of COVID-19 RNA-Dependent RNA Polymerase Bound to Favipiravir; RCSB Protein Data Bank: Piscataway, NJ, USA, 2021. [Google Scholar] [CrossRef]
  86. Golan, Y.; Campos, J.A.S.; Woolson, R.; Cilla, D.; Hanabergh, R.; Gonzales-Rojas, Y.; Lopez, R.; Finberg, R.; Balboni, A. Favipiravir in Patients With Early Mild-to-moderate Coronavirus Disease 2019 (COVID-19): A Randomized Controlled Trial. Clin. Infect. Dis. 2023, 76, e10–e17. [Google Scholar] [CrossRef] [PubMed]
  87. Sirijatuphat, R.; Manosuthi, W.; Niyomnaitham, S.; Owen, A.; Copeland, K.K.; Charoenpong, L.; Rattanasompattikul, M.; Mahasirimongkol, S.; Wichukchinda, N.; Chokephaibulkit, K. Early treatment of Favipiravir in COVID-19 patients without pneumonia: A multicentre, open-labelled, randomized control study. Emerg. Microbes Infect. 2022, 11, 2197–2206. [Google Scholar] [CrossRef]
  88. Eloy, P.; Le Grand, R.; Malvy, D.; Guedj, J. Combined treatment of molnupiravir and favipiravir against SARS-CoV-2 infection: One + zero equals two? eBioMedicine 2021, 74, 103663. [Google Scholar] [CrossRef]
  89. Huang, W.; Lu, S.; Huang, Z.; Liu, X.; Mou, L.; Luo, Y.; Zhao, Y.; Liu, Y.; Chen, Z.; Hou, T.; et al. Allosite: A method for predicting allosteric sites. Bioinformatics 2013, 29, 2357–2359. [Google Scholar] [CrossRef]
  90. Tian, H.; Jiang, X.; Tao, P. PASSer: Prediction of allosteric sites server. Mach. Learn. Sci. Technol 2021, 2, 035015. [Google Scholar] [CrossRef] [PubMed]
  91. Faisal, S.; Badshah, S.L.; Kubra, B.; Sharaf, M.; Emwas, A.-H.; Jaremko, M.; Abdalla, M. Identification and Inhibition of the Druggable Allosteric Site of SARS-CoV-2 NSP10/NSP16 Methyltransferase through Computational Approaches. Molecules 2022, 27, 5241. [Google Scholar] [CrossRef]
  92. Kijewska, M.; Sharfalddin, A.A.; Jaremko, Ł.; Cal, M.; Setner, B.; Siczek, M.; Stefanowicz, P.; Hussien, M.A.; Emwas, A.-H.; Jaremko, M. Lossen Rearrangement of p-Toluenesulfonates of N-Oxyimides in Basic Condition, Theoretical Study, and Molecular Docking. Front. Chem. 2021, 9, 662533. [Google Scholar] [CrossRef]
  93. Kohn, W.; Sham, L.J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. [Google Scholar] [CrossRef]
  94. Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. Proc. Natl. Acad. Sci. USA 1979, 76, 6062–6065. [Google Scholar] [CrossRef]
  95. Becke, A.D. A new mixing of Hartree–Fock and local density-functional theories. J. Chem. Phys. 1993, 98, 1372. [Google Scholar] [CrossRef]
  96. Becke, A.D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. [Google Scholar] [CrossRef]
  97. Møller, C.; Plesset, M.S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618–622. [Google Scholar] [CrossRef]
  98. Zhao, Y.; Truhlar, D.G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2007, 120, 215–241. [Google Scholar] [CrossRef]
  99. Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilizaion of AB initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117–129. [Google Scholar] [CrossRef]
  100. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.; Nakatsuji, H.; et al. Gaussian 16 Rev. C.01; Gaussian, Inc.: Wallingford, CT, USA, 2016. [Google Scholar]
  101. Parr, R.G.; Szentpály, L.v.; Liu, S. Electrophilicity Index. J. Am. Chem. Soc. 1999, 121, 1922–1924. [Google Scholar] [CrossRef]
  102. Gázquez, J.L.; Cedillo, A.; Vela, A. Electrodonating and Electroaccepting Powers. J. Phys. Chem. A 2007, 111, 1966–1970. [Google Scholar] [CrossRef]
  103. Chattaraj, P.K.; Chakraborty, A.; Giri, S. Net Electrophilicity. J. Phys. Chem. A 2009, 113, 10068–10074. [Google Scholar] [CrossRef]
  104. Koopmans, T. Über die Zuordnung von Wellenfunktionen und Eigenwerten zu den Einzelnen Elektronen Eines Atoms. Physica 1934, 1, 104–113. [Google Scholar] [CrossRef]
  105. Bader, R. Atoms in Molecules: A Quantum Theory (International Series of Monographs on Chemistry); Oxford University Press: Cary, NC, USA, 1994. [Google Scholar]
  106. Johnson, E.R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A.J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498–6506. [Google Scholar] [CrossRef] [PubMed]
  107. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 14, 33–38. [Google Scholar] [CrossRef] [PubMed]
  108. Trott, O.; Olson, A.J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2009, 31, 455–461. [Google Scholar] [CrossRef] [PubMed]
  109. Eberhardt, J.; Santos-Martins, D.; Tillack, A.F.; Forli, S. AutoDock Vina 1.2.0: New Docking Methods, Expanded Force Field, and Python Bindings. J. Chem. Inf. Model. 2021, 61, 3891–3898. [Google Scholar] [CrossRef] [PubMed]
  110. Yang, J.M.; Chen, C.C. GEMDOCK: A Generic Evolutionary Method for Molecular Docking. Proteins 2004, 55, 288–304. [Google Scholar] [CrossRef]
  111. Diedrich, K.; Krause, B.; Berg, O.; Rarey, M. PoseEdit: Enhanced ligand binding mode communication by interactive 2D diagrams. J. Comput. Aided Mol. Des. 2023, 37, 481–503. [Google Scholar] [CrossRef]
  112. Dassault Systèmes, Discovery Studio Visualizer, version v21.1.0.20298; Dassault Systèmes; BIOVIA: San Diego, CA, USA, 2021.
  113. Gehlhaar, D.K.; Bouzida, D.; Rejto, P.A. Fully automated and rapid flexible docking of inhibitors covalently bound to serine proteases. Lect. Notes Comput. Sci. 1198, 1447, 449–461. [Google Scholar] [CrossRef]
Figure 1. Structural formula of Favipiravir (6-fluoro-3-hydroxypyrazine-2-carboxamide, FVP, R = F) and its analogues R = H, Cl, Br, I, CH3, CF3 and CN (pro-drug (left), ribofuranosyl-5′-triphosphate active forms (right)) and the Tanimoto coefficient (molecular similarity distance).
Figure 1. Structural formula of Favipiravir (6-fluoro-3-hydroxypyrazine-2-carboxamide, FVP, R = F) and its analogues R = H, Cl, Br, I, CH3, CF3 and CN (pro-drug (left), ribofuranosyl-5′-triphosphate active forms (right)) and the Tanimoto coefficient (molecular similarity distance).
Molecules 29 00441 g001
Figure 2. Conversion of favipiravir and its analogues to its various metabolites (a series of ribosylation and phosphorylation steps to form the active triphosphate FVP-RTP).
Figure 2. Conversion of favipiravir and its analogues to its various metabolites (a series of ribosylation and phosphorylation steps to form the active triphosphate FVP-RTP).
Molecules 29 00441 g002
Figure 3. The ligand’s ability to accept (R+) and donate (R) charge in relation to FVP: (a) MP2/6-311G(d,p) and (b) M062X/6-311G(d,p) (the ordering of the substituents according to group electronegativity).
Figure 3. The ligand’s ability to accept (R+) and donate (R) charge in relation to FVP: (a) MP2/6-311G(d,p) and (b) M062X/6-311G(d,p) (the ordering of the substituents according to group electronegativity).
Molecules 29 00441 g003
Figure 4. The binding mode of FVP-RTP to RdRp in the 7CTT [81] complex.
Figure 4. The binding mode of FVP-RTP to RdRp in the 7CTT [81] complex.
Molecules 29 00441 g004
Figure 5. The docked poses of the ligands in the RdRp binding site. The protein backbone is represented as a cartoon, the binding cavity residues are shown as thin sticks, and the docked ligands are shown as color sticks (FVP in cyan, CN in yellow, and CF3 in red).
Figure 5. The docked poses of the ligands in the RdRp binding site. The protein backbone is represented as a cartoon, the binding cavity residues are shown as thin sticks, and the docked ligands are shown as color sticks (FVP in cyan, CN in yellow, and CF3 in red).
Molecules 29 00441 g005
Figure 6. The binding affinity versus ligand’s ability to accept (R+) and donate (R) charge.
Figure 6. The binding affinity versus ligand’s ability to accept (R+) and donate (R) charge.
Molecules 29 00441 g006
Figure 7. The comparison of the binding energies of FVR-RTP analogues with the RdRp split to active site residues, RNA primer, RNA template, and cofactor. (FVP*—native ligand).
Figure 7. The comparison of the binding energies of FVR-RTP analogues with the RdRp split to active site residues, RNA primer, RNA template, and cofactor. (FVP*—native ligand).
Molecules 29 00441 g007
Figure 8. The overlap of the isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) FVP (red) and the CF3 derivative (green) and (b) FVP (red) and the CN derivative (green). (Pre-Catalytic State—Productive Mode I).
Figure 8. The overlap of the isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) FVP (red) and the CF3 derivative (green) and (b) FVP (red) and the CN derivative (green). (Pre-Catalytic State—Productive Mode I).
Molecules 29 00441 g008
Figure 9. The binding mode of RDV (7UO4 [80]).
Figure 9. The binding mode of RDV (7UO4 [80]).
Molecules 29 00441 g009
Figure 10. The docked poses of the ligands in the binding site of RdRp (target from 7UO4 [80]). The protein backbone is represented as a cartoon, the binding cavity residues are shown as thin sticks, and the docked ligands are shown as sticks.
Figure 10. The docked poses of the ligands in the binding site of RdRp (target from 7UO4 [80]). The protein backbone is represented as a cartoon, the binding cavity residues are shown as thin sticks, and the docked ligands are shown as sticks.
Molecules 29 00441 g010
Figure 11. The comparison of the binding energies of FVR-RTP analogues with RdRp (selected residues), RNA primer, RNA template, and cofactor (target from 7UO4 [80]). (RDV*—native ligand).
Figure 11. The comparison of the binding energies of FVR-RTP analogues with RdRp (selected residues), RNA primer, RNA template, and cofactor (target from 7UO4 [80]). (RDV*—native ligand).
Molecules 29 00441 g011
Figure 12. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 analogue and (b) the CN analogue.(Pre-Catalytic State—Productive Mode II).
Figure 12. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 analogue and (b) the CN analogue.(Pre-Catalytic State—Productive Mode II).
Molecules 29 00441 g012
Figure 13. The binding mode of FVP-RTP in 7AAP [84].
Figure 13. The binding mode of FVP-RTP in 7AAP [84].
Molecules 29 00441 g013
Figure 14. The docked poses of the ligands in the binding site of RdRp (target from 7AAP [84]). The protein backbone is represented as a cartoon; the binding cavity residues are shown as thin sticks, and the docked ligands are shown as sticks (FVP-RTP in cyan, the CF3 analogue in pink, and the CN analogue in white).
Figure 14. The docked poses of the ligands in the binding site of RdRp (target from 7AAP [84]). The protein backbone is represented as a cartoon; the binding cavity residues are shown as thin sticks, and the docked ligands are shown as sticks (FVP-RTP in cyan, the CF3 analogue in pink, and the CN analogue in white).
Molecules 29 00441 g014
Figure 15. The comparison of the binding energies of FVR-RTP analogues to RdRp (selected residues), RNA primer, RNA template, and cofactor (target from 7AAP [84]) (FVP-RTP*—native ligand).
Figure 15. The comparison of the binding energies of FVR-RTP analogues to RdRp (selected residues), RNA primer, RNA template, and cofactor (target from 7AAP [84]) (FVP-RTP*—native ligand).
Molecules 29 00441 g015
Figure 16. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 derivative and (b) the CN derivative (Pre-Catalytic State—Non-Productive Mode).
Figure 16. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 derivative and (b) the CN derivative (Pre-Catalytic State—Non-Productive Mode).
Molecules 29 00441 g016
Figure 17. The comparison of the binding energies of the FVR-RMP analogues to RdRp (selected residues) and cofactor (target from 7DFG [85]). (FVP-RMP*—native ligand).
Figure 17. The comparison of the binding energies of the FVR-RMP analogues to RdRp (selected residues) and cofactor (target from 7DFG [85]). (FVP-RMP*—native ligand).
Molecules 29 00441 g017
Figure 18. The docked poses of the ligands in the binding site of RdRp (Active state). The protein backbone is represented as a cartoon; the binding cavity residues are shown as thin sticks and docked ligands are shown as sticks (FVP in pink, the CF3 analogue in yellow, and the CN analogue in green).
Figure 18. The docked poses of the ligands in the binding site of RdRp (Active state). The protein backbone is represented as a cartoon; the binding cavity residues are shown as thin sticks and docked ligands are shown as sticks (FVP in pink, the CF3 analogue in yellow, and the CN analogue in green).
Molecules 29 00441 g018
Figure 19. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 derivative and (b) the CN derivative. (Active state).
Figure 19. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 derivative and (b) the CN derivative. (Active state).
Molecules 29 00441 g019
Figure 20. The docked poses of the ligands in the binding site of RdRp (the RdRp target from 6M71 [74]). The protein backbone is represented as a cartoon; the binding cavity residues are shown as thin sticks and docked ligands are shown as sticks (FVP in pink, the CF3 analogue in yellow, and the CN analogue in green).
Figure 20. The docked poses of the ligands in the binding site of RdRp (the RdRp target from 6M71 [74]). The protein backbone is represented as a cartoon; the binding cavity residues are shown as thin sticks and docked ligands are shown as sticks (FVP in pink, the CF3 analogue in yellow, and the CN analogue in green).
Molecules 29 00441 g020
Figure 21. The comparison of the binding energies of the FVR-RTP analogues to RdRp (the RdRp target from 6M71 [74]).
Figure 21. The comparison of the binding energies of the FVR-RTP analogues to RdRp (the RdRp target from 6M71 [74]).
Molecules 29 00441 g021
Figure 22. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 derivative and (b) the CN derivative. (Allosteric effect).
Figure 22. The isosurfaces of RDG (isovalue 0.5a.u.) with sign(λ2BCP mapped over the surface for (a) the CF3 derivative and (b) the CN derivative. (Allosteric effect).
Molecules 29 00441 g022
Figure 23. The comparison of the binding energies of the FVR-RTP analogues to RdRp of all residues, cofactors Mg2+ and Zn2+, RNA Template, and RNA Primer; (a) pre-catalytic productive mode I; (b) pre-catalytic productive mode II; (c) non-productive mode. Pop-up peaks represent components with high binding energy.
Figure 23. The comparison of the binding energies of the FVR-RTP analogues to RdRp of all residues, cofactors Mg2+ and Zn2+, RNA Template, and RNA Primer; (a) pre-catalytic productive mode I; (b) pre-catalytic productive mode II; (c) non-productive mode. Pop-up peaks represent components with high binding energy.
Molecules 29 00441 g023
Table 1. The physicochemical profile parameters that describe the pharmacokinetic behavior of the candidate ligands.
Table 1. The physicochemical profile parameters that describe the pharmacokinetic behavior of the candidate ligands.
RMWLogP *LogP **Solubility (SILICOS-IT)Abbot’s BioavailabilityTPSA [A2]SASGILead-likeness (Lipinski, Egan, Veber)
H139.11−0.65−1.539Soluble0.5588.841.73HighYes
F (FVP)157.1−0.27−0.934Soluble0.5588.842.08HighYes
Cl173.56−0.05−0.472Soluble0.5588.842.07HighYes
Br218.010.04−0.293Soluble0.5588.842.12HighYes
I265.010.06−0.313Soluble0.5588.842.43HighYes
CF3207.110.52−0.142Soluble0.5588.842.03HighYes
CH3153.14−0.31−0.982Soluble0.5588.841.88HighYes
CN164.12−0.83−0.878Soluble0.55112.632.10HighYes
MOL329.31−0.89−0.402Soluble0.55143.144.49LowNo; 1 violation (MW > 300)
RVD602.230.521.961Moderately soluble0.17203.556.43LowNo; 2 violations (MW > 350, rotors > 7, TPSA > 140)
NMVr499.241.64 2.013Moderately soluble0.55131.44.82HighNo; 1 violation (rotors > 10, MW > 350)
RTV720.945.04 4.5Insoluble0.17202.266.45LowNo, 3 violations
(MW > 350, XLOGP3 > 3.5, rotors > 7)
Optimal values<5001 < logP < 4Soluble1.0<140<10HighYes
* SwissADME, ** ADMET2.0.
Table 2. The parameters describing the toxicity of the candidate ligands.
Table 2. The parameters describing the toxicity of the candidate ligands.
RBBB, PAINSBRENKP-gp CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4 CarcinogenicityGenotoxic
Carcinogenicity/Mutagenicity
AMESLD50 RatH-HTDILI
HNoNoNoNo0.13400.0220.2690.1770.956
F (FVP)NoNoNoNo0.2300.0310.5220.8330.935
ClNoNoNoNo0.50100.0270.6360.210.92
BrNoNoNoNo0.70800.0550.9620.0820.835
INoNoNoNo0.16900.0320.4870.1010.603
CF3NoNoNoNo0.1400.3260.9670.4910.948
CH3NoNoNoNo0.15800.0180.3640.3310.893
CNNoNoNoNo0.16800.0280.4340.8930.941
MOLNoNoNoNo0.3034 rules violation0.5120.0210.4860.976
RVDNo 1 alert (phosphorus)YesNo CYP3A4 Yes 0.3145 rules violation0.8220.7190.2570.946
NMVrNoNoYes No CYP3A4 Yes0.11900.0030.9620.2930.556
RTVNoNoYesNo CYP3A4 Yes0.0211 rule violation0.0230.2450.9650.978
Optimal valuesNoNoYesNo000000
Table 3. The global reactivity descriptors calculated at the MP2/6-311G(d,p) level of the theory.
Table 3. The global reactivity descriptors calculated at the MP2/6-311G(d,p) level of the theory.
R *HOMO [eV]LUMO [eV]Gap
[eV]
Χ
[eV]
η
[eV]
ω
[eV]
S
[1/eV]
ω+
[eV]
ω
[eV]
Δω
[eV]
ΔNmax
[−]
R+
[−]
R
[−]
single molecule (gas)
H−9.7031.67211.3754.0165.6881.4180.1764.1360.1214.2570.7060.5670.914
F−9.8641.23711.1014.3135.5501.6760.1804.5260.2134.7400.7771.0001.000
Cl−9.7631.24511.0074.2595.5041.6480.1824.4650.2064.6720.7740.9680.987
Br−9.6411.27010.9114.1865.4551.6060.1834.3800.1954.5750.7670.9140.968
I−9.4651.21810.6834.1245.3411.5920.1874.3210.1984.5190.7720.9270.955
CF3−10.3661.13811.5054.6145.7521.8510.1744.8770.2635.1390.8021.2321.077
CH3−9.3661.72311.0903.8225.5451.3170.1803.9210.0994.0200.6890.4660.866
CN−10.2011.00811.2094.5975.6041.8850.1784.8840.2875.1710.8201.3481.079
aqueous solution, pH = 7
H−9.5821.67111.2543.9555.6271.3900.1784.0710.1164.1870.7030.6100.929
F−9.6861.30310.9894.1925.4941.5990.1824.3820.1904.5720.7631.0001.000
Cl−9.6261.34510.9714.1405.4861.5620.1824.3180.1784.4960.7550.9370.986
Br−9.5441.35810.9024.0935.4511.5370.1834.2650.1724.4360.7510.9040.973
I−9.3151.39010.7053.9625.3521.4670.1874.1170.1554.2710.7400.8140.940
CF3−10.1041.32611.4304.3895.7151.6850.1754.5940.2054.8000.7681.0811.049
CH3−9.2911.67410.9653.8095.4821.3230.1823.9120.1044.0160.6950.5470.893
CN−9.9371.28111.1694.3285.6091.6700.1784.5350.2074.7420.7721.0891.035
* R—substituent.
Table 4. The global reactivity descriptors calculated at the M062X/6-311G(d,p) level of the theory.
Table 4. The global reactivity descriptors calculated at the M062X/6-311G(d,p) level of the theory.
R *HOMO [eV]LUMO [eV]Gap
[eV]
Χ
[eV]
η
[eV]
ω
[eV]
S
[1/eV]
ω+
[eV]
ω
[eV]
Δω
[eV]
ΔNmax
[−]
R+
[−]
R
[−]
single molecule (gas)
H−8.750−1.0807.6694.9153.8353.1500.2616.0871.1727.2581.2820.8050.926
F−8.799−1.4327.3675.1153.6843.5520.2716.5701.4558.0251.3891.0001.000
Cl−8.722−1.4777.2465.1003.6233.5890.2766.5921.4928.0841.4081.0261.003
Br−8.574−1.4637.1125.0183.5563.5410.2816.4951.4777.9721.4111.0150.989
I−8.422−1.5016.9204.9613.4603.5570.2896.4701.5097.9791.4341.0370.985
CF3−9.258−1.5827.6765.4203.8383.8270.2617.0171.5978.6141.4121.0981.068
CH3−8.434−0.9977.4374.7153.7182.9900.2695.8121.0976.9091.2680.7540.885
CN−9.188−1.7117.4785.4493.7393.9710.2677.1631.7148.8771.4571.1781.090
aqueous solution, pH = 7
H−8.675−1.0997.5764.8873.7883.1520.2646.0691.1827.2521.2900.8300.940
F−8.668−1.4007.2685.0343.6343.4860.2756.4581.4247.8811.3851.0001.000
Cl−8.634−1.4137.2205.0233.6103.4950.2776.4581.4357.8931.3911.0081.000
Br−8.533−1.4107.1224.9723.5613.4700.2816.4011.4307.8311.3961.0040.991
I−8.343−1.3806.9634.8613.4823.3940.2876.2601.3987.6581.3960.9820.969
CF3−9.098−1.4437.6555.2713.8283.6290.2616.7431.4728.2151.3771.0341.044
CH3−8.402−1.0697.3334.7353.6663.0580.2735.8841.1497.0331.2920.8070.911
CN−8.977−1.4927.4855.2353.7433.6610.2676.7461.5118.2571.3991.0611.045
* R—substituent.
Table 5. The degree of interaction, DOI, for FVP and its two analogues, CN and CF3.
Table 5. The degree of interaction, DOI, for FVP and its two analogues, CN and CF3.
FVPCF3CN
1N(4)3.531N(4)3.545N(4)3.537
2C(2)3.964C(2)4.279C(2)4.251
3C(6)5.151C(6)5.715C(6)5.897
4N(1)3.695N(1)3.783N(1)3.711
5C(2)5.124C(2)5.266C(2)5.048
6C(7)5.294C(7)5.453C(7)5.390
7N(7)3.037N(7)3.034N(7)3.049
8=O2.293=O2.291=O2.284
9C(3)5.154C(3)5.378C(3)5.277
10OH2.492OH2.503OH2.502
11F1.495C5.128C5.164
11a--F1.555N2.684
11b--F1.511--
11c--F1.510--
12H0.795H0.810H0.796
13H0.858H0.854H0.852
14H(OH)0.896H(OH)0.895H(OH)0.899
15H1.013H1.007H0.998
Table 6. The list of hydrogen bonds binding FVP-RTP to RdRp in the 7CTT [81] complex.
Table 6. The list of hydrogen bonds binding FVP-RTP to RdRp in the 7CTT [81] complex.
No.Residue ChainResidue NameH⋯A
Distance [Å]
D⋯A
Distance [Å]
<DHATypeDonorAcceptorMoiety
1621ALYS3.304.05133.53D4327 [Nam]8910 [O3]Phosphate
2622ACYS3.033.63120.57D4336 [Nam]8913 [O2]Phosphate
3623AASP2.603.46147.96A8902 [O3]4349 [O2]Ribosyl
4623AASP2.703.59153.00A8903 [O3]4349 [O2]Ribosyl
5623AASP2.823.74156.06D4342 [Nam]8913 [O2]Phosphate
6682ASER2.092.99150.84A8898 [Nam]4792 [O2]Pyridazine
7691AASN2.723.19109.79D4855 [Nam]8902 [O3]Ribosyl
8760AASP3.003.77140.59D5420 [O3]8901 [O3]Ribosyl
Table 7. The list of salt bridges (opposite charges), π-cation interactions (a positive charge and an aromatic ring), and metallic interactions in the 7CTT [81] complex.
Table 7. The list of salt bridges (opposite charges), π-cation interactions (a positive charge and an aromatic ring), and metallic interactions in the 7CTT [81] complex.
No.Residue ChainResidue NameDistance [Å]Ligand
Moiety
Interaction Type
1555AARG4.12 Aromaticπ-cation
2553AARG5.42Phosphateionic/salt bridge
3553AARG4.78Phosphateionic/salt bridge
4555AARG5.05Phosphateionic/salt bridge
5621ALYS4.36Phosphateionic/salt bridge
6798ALYS5.32Phosphateionic/salt bridge
7798ALYS5.44Phosphateionic/salt bridge
8798ALYS3.73Phosphateionic/salt bridge
9 Mg2+6.94Phosphatemetal
Table 8. The docking results for FVR-RTP analogues (the RdRp target from 7CTT [81]).
Table 8. The docking results for FVR-RTP analogues (the RdRp target from 7CTT [81]).
ParameterFVP *FClBrIHCH3CF3CN
binding score **, kcal/mol-−6.92−6.95−6.94−6.90−6.54−6.96−10.51−11.58
docking score, kcal/mol-−125.466−126.032−125.752−125.343−118.565−126.244−190.616−209.746
protein–ligand, kcal/mol−112.71−132.906−134.130−134.054−134.117−130.160−134.035−135.452−167.177
protein–ligand hydrogen bonds, kcal/mol−13.545−11.008−11.438−11.455−11.444−11.005−11.453−13.023−12.502
cofactor-RNA template/primer–ligand, kcal/mol−54.877−56.192−54.984−54.854−54.579−54.348−59.124−60.590−55.773
binding affinity, kJ/mol−24.66−29.05−25.31−25.31−25.30−25.60−25.26−76.80−58.57
* The actual ligand without redocking; ** AutoDock result.
Table 9. The ordering of the ligands according to the different parameters describing the binding strength.
Table 9. The ordering of the ligands according to the different parameters describing the binding strength.
ParameterOrdering
docking scoreCN > CF3 > CH3 > Cl > Br > I > H > F
total energy of bindingCN > CF3 > Cl > I > CH3 > Br > H > F
hydrogen bonds linking the ligand to the proteinCF3 > CN > Br > CH3 > I > Cl > F > H
binding of the ligand to the RNA templateCF3 > H > F > CN > Br > Cl > I > CH3
binding of the ligand to the RNA primerCF3 > F > CH3 > Cl > H > Br > I > CN
binding of the ligand to both the RNA template and RNA primerCF3 > F > CH3 > Cl > Br > H > I > CN
binding of the ligand to magnesium ionCN > CF3 > H ≅ F > Br > I > CH3 > Cl
binding affinityCF3 > CN > F > H > Cl ≅ Br > I > CH3
Table 10. Comparison of the similarity of binding patterns of the individual ligands with respect to the FPV-RTP ligand (Pre-Catalytic State—Productive Mode I, all residues).
Table 10. Comparison of the similarity of binding patterns of the individual ligands with respect to the FPV-RTP ligand (Pre-Catalytic State—Productive Mode I, all residues).
RThe Entire ComplexProtein ResiduesMg2+, Zn2+, RNA Template/Primer
EuclidianManhattanAdditiveEuclidianManhattanAdditiveEuclidianManhattanAdditive
H1.5670.6470.2401.5780.6270.2421.4070.9340.203
Cl1.5810.6360.2851.5510.9340.2031.5510.5880.289
Br1.5760.6350.2831.5490.5890.2891.9191.2940.199
I1.5660.6310.2781.5490.5900.2881.7911.2280.131
CF31.8770.7560.3781.7180.6910.3103.4261.6961.370
CH31.5900.6400.2881.5500.5880.2882.0911.3840.284
CN2.8451.0080.6612.9221.0080.7311.2861.005−0.352
FVP-RTP *1.5850.6360.2991.5290.5870.2792.2341.3360.585
* Redocked vs. actual.
Table 11. Comparison of the similarity of binding patterns of the individual ligands with respect to the FPV-RTF ligand (Pre-Catalytic State—Productive Mode I, selected residues).
Table 11. Comparison of the similarity of binding patterns of the individual ligands with respect to the FPV-RTF ligand (Pre-Catalytic State—Productive Mode I, selected residues).
RThe Entire Complex (Selected Residues, Mg2+, Zn2+, RNA Template, and RNA Primer)Selected Residues
(555, 798, 621, 553, 682, 545, 691, 622, 623, 760, 761, and 618)
EuclidianManhattanAdditiveEuclidianManhattanAdditive
H1.4920.4800.2131.4980.4490.214
Cl1.5230.4730.2731.4880.9340.203
Br1.5170.4720.2691.4850.4150.274
I1.5070.4690.2651.4860.4160.274
CF31.8030.5770.3771.6320.5000.308
CH31.5330.4770.2761.4860.4140.275
CN2.7290.7370.6342.8010.7190.703
FVP-RTP *1.5120.4700.2741.4480.4090.252
* Redocked vs. actual.
Table 12. The binding mode of the FVR-RTP analogues to RdRp (the RdRp target from 7CTT [81]); kcal/mol units. (The active site residues, RNA primer, RNA template, and cofactor are taken into account).
Table 12. The binding mode of the FVR-RTP analogues to RdRp (the RdRp target from 7CTT [81]); kcal/mol units. (The active site residues, RNA primer, RNA template, and cofactor are taken into account).
ResidueFVP-RTP *FClBrIHCH3CF3CN
RNA primer−41.019−45.808−44.858−44.732−44.396−43.848−45.166−48.652−42.647
RNA template−13.292−11.988−11.262−11.263−11.257−12.039−11.254−16.309−11.66
LYS545−5.379−5.136−4.984−4.986−4.991−4.686−4.981−9.728−5.070
ARG553−15.047−18.025−18.809−18.756−18.773−18.002−18.773−16.868−17.259
ARG555−31.286−30.126−30.431−30.376−30.360−27.792−30.436−32.469−31.030
ASP61814.76717.18915.05515.15815.18617.22515.11611.2674.697
LYS621−20.662−30.887−30.866−30.872−30.869−30.827−30.880−31.682−40.246
CYS622−4.298−6.501−6.469−6.472−6.462−6.520−6.463−6.449−7.593
ASP623−4.074−4.861−5.687−5.657−5.637−4.869−5.661−5.553−7.966
SER682−6.868−10.886−11.051−11.051−11.059−10.859−11.041−10.147−9.894
ASN691−4.639−5.911−6.279−6.276−6.266−5.905−6.276−2.873−6.288
SER759−0.740−1.690−1.857−1.861−1.852−1.686−1.861−0.752−1.061
ASP7602.6170.6551.8061.7761.6990.5701.7865.6992.603
ASP7614.6594.2424.2034.2064.2044.2404.2064.3633.125
LYS798−23.722−21.897−20.299−20.352−20.363−21.910−20.329−21.657−29.618
Mg2+−7.022−6.450−6.306−6.314−6.313−6.450−6.312−6.464−5.305
cosine distance (total) ** 1.00.9880.9850.9850.9850.9850.9850.9850.956
cosine distance (residues) **1.00.9800.9770.9770.9770.9750.9770.9740.939
* The actual ligand without redocking ** Cosine distance between the binding mode of FVP-RTP and its analogues.
Table 13. The list of most important hydrogen bonds binding RVD to RdRp in 7UO4 [80] complex.
Table 13. The list of most important hydrogen bonds binding RVD to RdRp in 7UO4 [80] complex.
No.Reside Chain Residue NameH⋯A
Distance [Å]
D⋯A
Distance [Å]
<DHADonor/AcceptorDonorAcceptorMoiety
1621ALYS3.174.09155.61D4992 [Nam]12374 [O3]Phosphate
2622ACYS3.124.04156.57D5001 [Nam]12373 [O3]Phosphate
3623AASP2.643.21117.47D5007 [Nam]12378 [O3]Phosphate
4623AASP2.153.11168.63A12,378 [O3]5014 [O]Phosphate
5687ATHR2.743.70168.20A12,379 [O3]5489 [O3]Phosphate
6691AASN2.002.80136.69D5520 [Nam]12379 [O3]Ribosyl
7759ASER2.913.59128.06D6082 [O3]12381 [N1]Ribosyl
Table 14. The list of salt bridges (opposite charges) and metallic interactions in 7UO4 [80].
Table 14. The list of salt bridges (opposite charges) and metallic interactions in 7UO4 [80].
No.Residue ChainResidue NameDistance [Å] Ligand MoietyInteraction Type
1551ALYS5.08Phosphatesalt bridge
2798ALYS3.58Phosphatesalt bridge
3 Mg2+2.34Phosphatemetal
4 Mg2+2.82Phosphatemetal
Table 15. The docking results for FVR-RTP analogues (the RdRp target from 7UO4 [80]).
Table 15. The docking results for FVR-RTP analogues (the RdRp target from 7UO4 [80]).
RDV *FClBrIHCH3CF3CN
Tanimoto binding fingerprint distance1.00.7350.7350.7350.7350.7350.7530.7530.735
binding score **, kcal/mol-−8.201−7.673−7.663−7.651−7.448−7.292−7.887−7.786
docking score, kcal/mol-−251.895−235.698−235.415−235.02−228.801−223.990−242.268−239.205
protein–ligand, total, kcal/mol−114.144−156.685−155.140−153.131−155.108−154.503−155.543−158.402−155.860
protein–ligand hydrogen bonds, kcal/mol−6.768−26.465−24.459−22.111−24.378−25.941−18.922−22.176−23.164
cofactor-RNA template/primer–ligand, total, kcal/mol−57.037 (−7.328)−55.375 (−23.783)−50.282 (−15.983)−53.411 (−14.170)−53.411 (−14.055)−45.374 (−15.276)−59.404 (−15.487)−53.785 (−10.571)−54.901 (−15.961)
hydrogen bonds, total, kcal/mol−13.439−26.484−23.114−22.107−22.250−22.756−18.649−18.649−22.189
binding affinity, kJ/mol-−112.947−84.893−82.002−83.860−84.413−82.009−70.516−83.455
* The actual ligand without redocking, 7 = ** AutoDock result.
Table 16. The ligands’ ordering according to the decreasing binding strength.
Table 16. The ligands’ ordering according to the decreasing binding strength.
ParameterOrdering
docking scoreF > CF3 > CN > Cl > Br > I > CH3 > H
total energy of bindingCF3 > F > CN > CH3 > Cl > I > H > Br
hydrogen bonds linking the ligand to the proteinF > H > Cl > I > CN > CF3 > Br > CH3
binding of the ligand to the RNA templateF > CF3 > CN > Br > I > H > Cl > CH3
binding of the ligand to the RNA primerCF3 > F > CN > Br > CH3 > Cl > I > H
binding of the ligand to both the RNA template and RNA primerF > CF3 > CN > Br > CH3 > Cl > I > H
binding of the ligand to magnesium ionF > H > Cl > Br > CN > I > CF3 > CH3
binding affinityF > Cl > H > I > CN > CH3 > Br > CF3
Table 17. The differences in binding mode of the particular FVP-RTP analogues relative to actual RVD ligand (Pre-Catalytic State—Productive Mode II, all residues).
Table 17. The differences in binding mode of the particular FVP-RTP analogues relative to actual RVD ligand (Pre-Catalytic State—Productive Mode II, all residues).
RComplexProteinMg2+, Zn2+, RNA Template/Primer
EuclidianManhattanAdditiveEuclidianManhattanAdditiveEuclidianManhattanAdditive
H3.0741.2540.1762.4191.0470.3598.0984.489−2.678
Cl2.8181.2080.2392.4491.0540.3706.1803.618−1.809
Br2.6751.1540.1472.4021.0230.3455.3623.192−2.934
I2.8331.2010.2342.4501.0470.3696.2913.615−1.871
CF32.4391.1520.1722.0710.9920.2205.6433.644−0.568
CH32.5531.0740.0601.5670.7930.4018.3585.454−5.257
CN2.6651.1690.3022.4971.0630.3814.5462.814−0.928
F2.7411.2600.4122.5091.0870.3925.1483.9560.733
Table 18. The comparison of the similarity of binding patterns of the individual ligands in relation to the actual RVD ligand (Pre-Catalytic State—Productive Mode II, selected residues).
Table 18. The comparison of the similarity of binding patterns of the individual ligands in relation to the actual RVD ligand (Pre-Catalytic State—Productive Mode II, selected residues).
RThe Entire Complex (Selected Residues, Mg2+, Zn2+, RNA Template, RNA Primer)Active Site Residues
(621, 553, 555, 798, 622, 551, 682, 691, 620,
687, 759, 545, 623, 760, 761, 618)
EuclidianManhattanAdditiveEuclidianManhattanAdditive
H3.0350.9820.1312.3660.7570.311
Cl2.7740.9300.2032.3960.7570.332
Br2.6320.8800.1072.3510.7310.302
I2.7920.9290.1912.3990.7560.323
CF32.1050.7700.1142.0050.6850.157
CH32.5120.7960.0021.4960.4980.340
CN2.6190.8890.2682.4440.7650.345
F2.7411.2600.4122.5091.0870.392
Table 19. The binding mode of the FVP-RTP analogues with RdRp (the RdRp target from 7UO4 [80]); kcal/mol units.
Table 19. The binding mode of the FVP-RTP analogues with RdRp (the RdRp target from 7UO4 [80]); kcal/mol units.
ResidueRDV *FClBrIHCH3CF3CN
RNA primer −46.462−46.702−41.536−42.079−41.228−37.102−41.599−48.096−45.449
RNA template−18.993−26.878−19.374−20.311−19.535−19.460−19.361−23.264−19.867
LYS545−4.088−3.222−3.223−3.367−3.268−3.089−3.158−3.124−3.540
LYS551−7.972−7.551−7.516−7.679−7.632−7.848−9.866−9.473−7.655
ARG553−8.612−22.002−22.002−21.063−21.955−21.964−18.994−9.599−22.066
ARG555−15.802−28.985−28.669−28.606−28.995−18.856−22.904−23.205−22.065
ASP6186.9616.8716.7276.9766.8646.8755.8946.4676.764
LYS621−13.844−20.934−21.102−20.496−21.156−28.162−26.845−19.236−20.700
CYS622−8.369−7.207−7.561−7.608−7.477−7.797−8.526−10.510−7.238
ASP623−0.638−2.238−2.176−1.930−1.997−1.692−3.646−1.565−2.476
SER682−10.710−14.323−13.310−12.666−13.161−13.968−13.635−13.314−12.987
THR687−4.667−13.176−13.316−12.602−12.985−12.841−11.555−10.987−13.311
ASP7602.6114.8294.5714.5054.5614.7841.2781.9874.550
ASP7614.1583.8253.8043.8453.8133.8223.3374.1563.804
LYS798−23.603−23.419−23.449−23.872−23.432−23.473−18.078−22.486−23.559
Mg2+−44.712−45.197−45.136−45.079−44.971−45.138−27.898−44.230−45.109
cosine distance (total) ** 1.00.9680.9640.9680.9630.9560.9360.9900.976
cosine distance (residues) ** 1.00.9400.9400.9470.9410.9340.9250.9760.946
cosine distance to 7CTT0.6280.5620.5510.5500.5460.5590.6100.5940.602
* Actual ligand without redocking ** Cosine distance between the binding mode of RDV and the FVP-RTP analogues.
Table 20. The docking results for the FVR-RTP analogues (the RdRp target from 7AAP [84]).
Table 20. The docking results for the FVR-RTP analogues (the RdRp target from 7AAP [84]).
ParameterFVP−RTP *FClBrIHCH3CF3CN
binding score **, kcal/mol−8.85−10.40−9.85−10.49−10.80−10.05−9.84−11.74−10.81
docking score, kcal/mol −161.033−189.604−179.185−190.741−196.385−182.759−179.006−213.435−196.565
protein–ligand, total, kcal/mol−44.038−42.534−46.801−50.122−46.242−44.699−46.246−58.546−44.98
protein–ligand, hydrogen bonds, kcal/mol−5.089−4.152−7.057−7.162−4.529−7.481−7.829−8.902−5.461
cofactor-RNA template/primer–ligand, total, kcal/mol−66.056−39.391−57.662−47.412−21.649−37.123−58.932−57.830−30.120
binding affinity, kJ/mol−25.460−79.409−79.795−85.965−83.816−82.637−76.086−100.253−72.851
* The actual ligand without redocking; ** AutoDock result.
Table 21. The differences in binding mode of the particular FVP-RTP analogues relative to actual ligand FVP-RTP (Pre-Catalytic State—Non-Productive Mode, all residues).
Table 21. The differences in binding mode of the particular FVP-RTP analogues relative to actual ligand FVP-RTP (Pre-Catalytic State—Non-Productive Mode, all residues).
RComplexProtein ResiduesMg2+, Zn2+, RNA Template/Primer
EuclidianManhattanAdditiveEuclidianManhattanAdditiveEuclidianManhattanAdditive
H2.2730.5890.2230.5970.242−0.0257.7624.6263.045
Cl1.8870.4900.2170.5740.2070.0226.3703.7782.449
Br2.0120.6090.5290.6310.2790.2796.7744.4553.460
I3.4460.8680.4780.6450.309−0.00411.9557.3415.985
CF34.0761.0650.6711.0420.5260.16013.9417.3836.557
CH31.5420.4350.1920.5910.2080.0275.0673.0812.091
CN3.0070.8010.3301.1160.337−0.0779.9276.2064.956
F2.6430.6770.3070.5550.248−0.0279.1315.6484.112
Table 22. Comparison of the similarity of binding patterns of individual ligands in relation to the current FVP-RTP ligand (Pre-Catalytic State—Non-Productive Mode, selected residues).
Table 22. Comparison of the similarity of binding patterns of individual ligands in relation to the current FVP-RTP ligand (Pre-Catalytic State—Non-Productive Mode, selected residues).
RThe Entire Complex (Selected Residues, Mg2+, Zn2+, RNA Template, RNA Primer)Active Site Residues
(545, 682, 555, 553, 814, 551, 691, 798,
687, 622, 760, 623, 618, 761)
EuclidianManhattanAdditiveEuclidianManhattanAdditive
H2.2710.5450.2210.5880.193−0.027
Cl1.8850.4520.2220.5680.1650.028
Br2.0080.5410.4610.6160.2040.204
I3.4390.7570.4710.6050.186−0.012
CF34.0620.8880.6470.9800.3270.134
CH31.5410.4090.1920.5870.1800.026
CN3.0040.7190.3281.1050.246−0.079
F2.6390.6010.3020.5350.163−0.032
Table 23. The binding affinity of the FVP-RTP analogues to RdRp (the RdRp target from 7AAP); kcal/mol units.
Table 23. The binding affinity of the FVP-RTP analogues to RdRp (the RdRp target from 7AAP); kcal/mol units.
ParameterFVP-RTP *FClBrIHCH3CF3CN
RNA primer −36.368−37.596−36.257−39.591−42.655−37.085−35.759−44.212−39.227
RNA template−26.551−21.173−22.008−23.071−21.805−21.016−23.694−23.658−22.178
LYS545−9.973−9.094−11.129−12.049−8.987−12.001−11.715−13.715−10.082
LYS551−4.067−3.951−3.622−3.900−3.746−3.623−3.796−3.826−4.172
ARG553−4.621−4.323−3.967−5.312−4.452−4.059−4.527−5.341−4.098
ARG555−8.070−6.249−7.762−6.433−6.640−5.964−6.792−9.788−8.180
ASP6185.2385.3315.2815.6373.6165.5755.5174.3445.941
LYS621−2.203−2.357−2.096−2.537−2.465−2.199−2.298−2.864−2.408
CYS622−0.565−0.819−0.771−0.734−0.822−0.818−0.581−0.960−0.995
ASP6231.8110.812−2.265−1.5511.5441.209−1.191−0.3391.148
SER682−8.537−8.580−7.740−7.216−9.875−8.078−8.172−6.880−8.805
THR687−2.717−1.879−1.163−0.515−2.18300−1.189−0.872
ASP7601.7890.3330.8780.6530.741−0.84300.301−0.847
ASP7616.1716.3776.7196.2676.3486.7526.4506.1986.789
LYS798−3.326−4.516−3.313−3.689−4.602−3.495−3.498−5.348−4.271
Mg2+−28.566−47.517−35.764−40.192−56.691−36.563−35.608−64.101−44.510
Mg2+−0.308−0.325−0.313−0.318−0.329−0.317−0.314−0.345−0.327
Mg2+−44.103−58.033−58.648−56.893−56.244−62.176−55.132−49.372−64.308
cosine distance (total) **0.9730.9830.9850.9630.9780.9870.9450.9770.973
cosine distance (residues) **0.9450.9390.9310.9430.9360.9360.9230.9460.945
* The actual ligand without redocking, ** Cosine distance between the binding mode of RDV and the FVP-RTP analogues.
Table 24. The docking results for FVP-RMP analogues (the RdRp target from 7DFG [85]).
Table 24. The docking results for FVP-RMP analogues (the RdRp target from 7DFG [85]).
ParameterFVP-RMP *CF3CN
binding score **, kcal/mol−29.20−35.51−35.90
docking score, kcal/mol−529.700−648.051−654.057
protein–ligand, total, kcal/mol−192.623−219.376−212.312
steric, total, kcal/mol−155.733−155.733−153.221
van der Waals, total, kcal/mol−55.040−61.936−61.538
protein–ligand hydrogen bonds, total, kcal/mol−21.747−26.440−23.089
cofactors, total, kcal/mol−1.076−3.126−1.662
binding affinity, kJ/mol−72.457−75.135−69.515
* The actual ligand, ** AutoDock.
Table 25. The binding mode of the FVP-RMP analogues to RdRp (target from 7DFG [85]); kcal/mol units.
Table 25. The binding mode of the FVP-RMP analogues to RdRp (target from 7DFG [85]); kcal/mol units.
ResidueFVP-RMP *CF3CN
LYS500−6.213−5.970−6.190
ARG513−9.294−4.533−4.611
LYS545−6.555−9.620−8.610
ARG555−5.812−9.247−7.042
SER682−13.242−13.258−13.578
SER759−6.137−5.961−6.014
CYS813−16.024−15.159−15.478
SER814−9.724−11.143−10.050
ARG836−32.449−38.565−36.865
LYS849−23.107−23.623−23.229
ARG858−24.517−28.346−29.649
SER861−8.122−8.890−10.14
Mg2+−23.746−23.381−23.695
Mg2+−12.684−13.075−12.743
cosine distance1.00.9910.993
* The actual ligand without redocking.
Table 26. The differences in binding mode of the particular FVP-RTP analogues relative to actual ligand FVP-RTP (Active stat, all residues).
Table 26. The differences in binding mode of the particular FVP-RTP analogues relative to actual ligand FVP-RTP (Active stat, all residues).
RComplexProtein ResiduesMg2+, Zn2+, RNA Template
EuclidianManhattanAdditiveEuclidianManhattanAdditiveEuclidianManhattanAdditive
CF30.9190.3810.1400.9330.3800.1660.8600.466−0.266
CN0.6980.3110.1180.7460.3250.1240.0520.1380.044
Table 27. Comparison of the similarity of binding patterns of individual ligands in relation to the current FVP-RTP ligand (Active state, selected residues).
Table 27. Comparison of the similarity of binding patterns of individual ligands in relation to the current FVP-RTP ligand (Active state, selected residues).
RThe Entire Complex (Selected Residues, Mg2+, Zn2+, RNA Template)Active Site Residues
(836, 858, 849, 813, 682, 814, 513, 861, 545, 500, 759, 555, 862)
EuclidianManhattanAdditiveEuclidianManhattanAdditive
CF30.6520.1840.0740.5030.1260.058
CN0.5060.1460.0790.4380.1130.052
Table 28. The binding mode of the FVP-RTP analogues to RdRp (the RdRp target from 6M71 [74]).
Table 28. The binding mode of the FVP-RTP analogues to RdRp (the RdRp target from 6M71 [74]).
FClBrIHCH3CF3CN
binding score *, kcal/mol−8.42−6.00−6.80−8.05−5.60−6.15−6.52−7.97
docking score, kcal/mol−152.723−108.945−123.513−146.582−101.65−111.641−118.237−144.458
protein–ligand, total, kcal/mol−157.797−121.337−143.673−156.914−127.477−124.751−129.888−149.934
protein–ligand hydrogen bonds, total−21.207−39.023−18.336−21.327−19.156−14.831−19.482−20.409
binding affinity, kJ/mol−34.145−34.102−34.983−34.245−32.689−30.110−40.309−31.148
* AutoDock.
Table 29. The binding mode of the FVP-RTP analogues to RdRp (the RdRp target from 6M71 [74]); kcal/mol units.
Table 29. The binding mode of the FVP-RTP analogues to RdRp (the RdRp target from 6M71 [74]); kcal/mol units.
ResidueFClBrIHCH3CF3CN
LYS545−33.460−26.092−26.185−33.609−27.282−24.831−32.410−26.185
ARG553−14.026−13.621−9.488−13.234−9.836−13.100−15.893−9.488
ARG555−20.555−21.398−23.292−19.977−19.827−22.241−20.148−23.292
ASP6181.0691.0881.0871.0611.1141.1251.110-
ASP623−10.604−6.813−9.559−10.063−4.021−4.222−9.262−9.559
SER682−17.852−6.008−19.222−17.973−15.860−9.322−0.381−19.222
THR687−12.302−12.053−12.409−12.342−7.133−10.746−10.406−12.409
ASN691−8.817−1.631−9.319−9.208−6.982−3.638−5.182−9.319
SER759−6.907−8.066−6.738−6.907−6.410−1.590−6.602−6.738
ASP760−3.664-−3.616−3.616−3.339−3.885−0.814−3.616
ASP7611.4841.5491.5491.4851.5721.6131.5671.550
LYS798−0.960−0.965−0.959−0.953−0.969−0.967−0.998−0.951
cosine distance **1.00.9530.9851.0000.9870.9640.9300.985
** Cosine distance calculated relative to the F analogue.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Latosińska, M.; Latosińska, J.N. Favipiravir Analogues as Inhibitors of SARS-CoV-2 RNA-Dependent RNA Polymerase, Combined Quantum Chemical Modeling, Quantitative Structure–Property Relationship, and Molecular Docking Study. Molecules 2024, 29, 441. https://doi.org/10.3390/molecules29020441

AMA Style

Latosińska M, Latosińska JN. Favipiravir Analogues as Inhibitors of SARS-CoV-2 RNA-Dependent RNA Polymerase, Combined Quantum Chemical Modeling, Quantitative Structure–Property Relationship, and Molecular Docking Study. Molecules. 2024; 29(2):441. https://doi.org/10.3390/molecules29020441

Chicago/Turabian Style

Latosińska, Magdalena, and Jolanta Natalia Latosińska. 2024. "Favipiravir Analogues as Inhibitors of SARS-CoV-2 RNA-Dependent RNA Polymerase, Combined Quantum Chemical Modeling, Quantitative Structure–Property Relationship, and Molecular Docking Study" Molecules 29, no. 2: 441. https://doi.org/10.3390/molecules29020441

APA Style

Latosińska, M., & Latosińska, J. N. (2024). Favipiravir Analogues as Inhibitors of SARS-CoV-2 RNA-Dependent RNA Polymerase, Combined Quantum Chemical Modeling, Quantitative Structure–Property Relationship, and Molecular Docking Study. Molecules, 29(2), 441. https://doi.org/10.3390/molecules29020441

Article Metrics

Back to TopTop