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

Favipiravir (6-fluoro-3-hydroxypyrazine-2-carboxamide, FPV), an active pharmaceutical component of the drug discovered and registered in March 2014 in Japan under the name Avigan, with an indication for pandemic influenza, has been studied. The study of this compound was prompted by the idea that effective processes of recognition and binding of FPV to the nucleic acid are affected predominantly by the propensity to form intra- and intermolecular interactions. Three nuclear quadrupole resonance experimental techniques, namely 1H-14N cross-relaxation, multiple frequency sweeps, and two-frequency irradiation, followed by solid-state computational modelling (density functional theory supplemented by the quantum theory of atoms in molecules, 3D Hirshfeld Surfaces, and reduced density gradient) approaches were applied. The complete NQR spectrum consisting of nine lines indicating the presence of three chemically inequivalent nitrogen sites in the FPV molecule was detected, and the assignment of lines to particular sites was performed. The description of the nearest vicinity of all three nitrogen atoms was used to characterize the nature of the intermolecular interactions from the perspective of the local single atoms and to draw some conclusions on the nature of the interactions required for effective recognition and binding. The propensity to form the electrostatic N−H···O, N−H···N, and C−H···O intermolecular hydrogen bonds competitive with two intramolecular hydrogen bonds, strong O−H···O and very weak N−H···N, closing the 5-member ring and stiffening the structure, as well as π···π and F···F dispersive interactions, were analysed in detail. The hypothesis regarding the similarity of the interaction pattern in the solid and the RNA template was verified. It was discovered that the -NH2 group in the crystal participates in intermolecular hydrogen bonds N–H···N and N–H···O, in the precatalytic state only in N–H···O, while in the active state in N–H···N and N–H···O hydrogen bonds, which is of importance to link FVP to the RNA template. Our study elucidates the binding modes of FVP (in crystal, precatalytic, and active forms) in detail and should guide the design of more potent analogues targeting SARS-CoV-2. Strong direct binding of FVP-RTP to both the active site and cofactor discovered by us suggests a possible alternative, allosteric mechanism of FVP action, which may explain the scattering of the results of clinical trials or the synergistic effect observed in combined treatment against SARS-CoV-2.


Introduction
The coronavirus disease-2019 (COVID- 19), an infectious disease associated with novel severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) [1,2], originated in December 2019 in the Hubei Province of China. Since the spread of the COVID-19 pandemic throughout the world, efforts to find effective treatments have been persistently undertaken. The recent endeavours in the field of drug design are aimed at discovering old or newly designed drugs which could be effective for the treatment of this contagious disease caused by different SARS-CoV-2 variants (Alpha, Beta, Gamma, Delta, Omicron) and subvariants, especially so-called variants of concern (VOC) lineages under WHO monitoring (BA.5, BA.2.75, BA.4.6, XBB, BA.2.3.20) [3].
Favipiravir (6-fluoro-3-hydroxypyrazine-2-carboxamide, FPV, T-705), Figure 1, is an active pharmaceutical ingredient of the drug discovered by the Japanese company Toyama Kagaku Kōgyō and registered in March 2014 in Japan under the name Avigan, with an indication for pandemic influenza [4][5][6][7]. It is well known for its in vitro activity towards OTV-resistant influenza A, B, and C viruses as well as flavi-, alpha-, filo-, bunya-, arena-, and noroviruses [8,9]. Preclinical trials performed before and after its first registration have shown strong FPV activity against several viruses containing negative-strand RNA as genetic material and causing serious human diseases [10], such as Ebola [11], Lassa virus [12], West Nile Fever [13], Zika [14], tick-borne encephalitis [15], and even rabies [16]. The antiviral activity of FPV consists of disrupting the normal function of the RNA-dependent RNA polymerase (RdRP), which produces a "mirror image" of a viral RNA, later used as a positive-sense RNA needed to synthesize multiple copies of nascent virus RNAs. RdRP malfunction 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. Thus, FVP is, in fact, a prodrug. It has only recently been discovered that the dominant mechanism of action of FVP is unlikely to be related to the delayed chain termination but to its action as a strong viral mutagen not inducing mutations beneficial for the virus [17], likely because these mutations cannot be effectively repaired by the viral proteins performing a "technical check" on the replication of virus' genetic material. Importantly, FVP is not used by human RNA polymerases, making it a low-toxicity drug. The therapeutic effect of FPV is a result of the accumulation of mutations in the replicated RNA of nascent viruses which lose their ability to proliferate. In Africa, FPV (along with the US-made drug Remdesivir, RSV, brand name Veklury) was clinically tested a few years ago as a potential drug for Ebola. Initially, results were promising, but the Ebola epidemic subsided and the trials were discontinued due to the lack of patients. However, in 2022, about 60 new cases of Ebola virus disease (EVD) were discovered in Mubende, Kassanda, and Kampala Districts in the Central Region, and Bunyangabu, Kagadi, and Kyegegwa Districts in the Western Region of Uganda. It is well known for its in vitro activity towards OTV-resistant influenza A, B, and C viruses as well as flavi-, alpha-, filo-, bunya-, arena-, and noroviruses [8,9]. Preclinical trials performed before and after its first registration have shown strong FPV activity against several viruses containing negative-strand RNA as genetic material and causing serious human diseases [10], such as Ebola [11], Lassa virus [12], West Nile Fever [13], Zika [14], tick-borne encephalitis [15], and even rabies [16]. The antiviral activity of FPV consists of disrupting the normal function of the RNA-dependent RNA polymerase (RdRP), which produces a "mirror image" of a viral RNA, later used as a positive-sense RNA needed to synthesize multiple copies of nascent virus RNAs. RdRP malfunction 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. Thus, FVP is, in fact, a prodrug. It has only recently been discovered that the dominant mechanism of action of FVP is unlikely to be related to the delayed chain termination but to its action as a strong viral mutagen not inducing mutations beneficial for the virus [17], likely because these mutations cannot be effectively repaired by the viral proteins performing a "technical check" on the replication of virus' genetic material. Importantly, FVP is not used by human RNA polymerases, making it a low-toxicity drug. The therapeutic effect of FPV is a result of the accumulation of mutations in the replicated RNA of nascent viruses which lose their ability to proliferate. In Africa, FPV (along with the US-made drug Remdesivir, RSV, brand name Veklury) was clinically tested a few years ago as a potential drug for Ebola. Initially, results were promising, but the Ebola epidemic subsided and the trials were discontinued due to the lack of patients. However, in 2022, about 60 new cases of Ebola virus disease (EVD) were discovered in Mubende, Kassanda, and Kampala Districts in the Central Region, and Bunyangabu, Kagadi, and Kyegegwa Districts in the Western Region of Uganda.
After the identification of the new virus, SARS-CoV-2, FVP was one of the first drugs whose effectiveness and safety were tested on the COVID-19 disease caused by this virus. This drug, along with the aforementioned RSV, raised the greatest hopes for the effective treatment of COVID-19. Compared to Molnupiravir or Nirmatrelvir, FVP is less potent against the SARS-CoV-2 virus in vitro [18,19] but similarly potent in animal models of possible competition between interactions raises the question of whether the experimentally observed local hydrogen-bonded environment in the solid state corresponds to the most likely interactions required for effective processes of FPV recognition and binding to the proteins or nucleic acids (DNA or RNA). We anticipated that our study would clarify these aspects and help in outlining the directions of synthesis of new active pharmaceutical ingredients that would be promising in the fields of the chemical industry and drug design of more effective analogues that are promising for COVID-19 and Ebola treatments.

1H-14 N NQDR Spectrum
The 14 N NQR frequencies in FVP lie in the low-frequency range below 4 MHz, and direct NQR measurement on a small sample (approx. 0.5g) is not possible due to a very low signal-to-noise ratio. Instead, the 14 N NQR frequencies were determined indirectly by measuring the proton NMR signal using nuclear quadrupole double resonance techniques. A multistep procedure was used.
In the first step, the 1 H-14 N cross-relaxation (CR) [39] spectrum was measured using fast field cycling (FFC) relaxometry. In a standard CR experiment, three successive time intervals with different values of a static magnetic field are used: (i) Polarization field B P for time t P . Typically, t P ≥ 3T 1H (B P ), where T 1H (B P ) is the proton spin-lattice relaxation time in the field B P so that an equilibrium proton magnetization is established. (ii) Relaxation (also called "mixing") field B R for time τ. Here the proton magnetization relaxes towards the equilibrium value with a time constant T 1H (B R ). (iii) Acquisition field B A , at which the proton NMR signal amplitude, which is proportional to the proton magnetization at the end of B R interval, is measured.
By repeating the experiment with a different relaxation field B R and relaxation time τ, the full dispersion of spin-lattice relaxation time T 1H (B R ) is determined. At values of magnetic field B R , where the proton NMR Larmor frequency ν H coincides with one of the 14 N NQR frequencies ν Q or ν H = γ H 2π B R = ν Q , a decrease in T 1H (B R ) is often observed. This phenomenon is well known as "quadrupole dips" [40] and is often used to indirectly determine the NQR frequencies.
In the case of FVP, the initial measurements have shown that T 1H (B P ) is longer than 100 s in the polarization field B P = 0.5 T, which is typically used in our laboratory. At such long values, the power consumption of the relaxometer is too high, and the sample cannot be properly polarised. Additionally, long T 1H makes measurement very time-consuming. Therefore, we modified the technique in such a way that the polarization period was removed from the experiment (B P = 0), and the proton magnetization was zero at the beginning of the B R interval. Then, the proton magnetization grows towards the equilibrium value as M = M 0 1 − e − τ T 1H (B R ) . In addition, for values of B R at which the 1 H proton NMR transition frequency ν H matches one of the 14 N NQR transition frequencies ν Q , the proton system is additionally polarised by the transfer of polarization from nitrogen NQR levels to proton NMR levels with a characteristic cross-relaxation time T CR . If T CR < T 1H (B R ), additional "quadrupole peaks" thus arise in the CR spectrum.
The 1 H-14 N cross-relaxation spectrum of FVP, as measured by the modified crossrelaxation technique, is presented in Figure 2. Here, a fixed value of τ = 3 s was used, and the measurement was repeated so that the "relaxation" of the magnetic field B R was scanned in steps of 0.585 mT, which corresponds to a step in the proton NMR frequency ∆ν H equal to 25 kHz. The proton NMR frequency in the acquisition field B A was 9.25 MHz. Three strong peaks centred at the frequencies of 3.870 MHz, 3.610 MHz, and 2.940 MHz are clearly observed. These peaks are observed in the frequency range in which the higher 14 N NQR frequencies ν+ and ν− of pyrazine nitrogen are usually observed. Since there are two non-equivalent pyrazine nitrogen atoms in the molecule, we expect to observe four peaks. It is possible that two peaks overlap within the experimental resolution.
Two somewhat weaker peaks are observed at the frequencies of 2.140 MHz and 1.610 MHz. In this frequency range, the amide nitrogen NQR frequencies ν+ and ν− are usually observed. Below 1 MHz, there is a well-observed peak at 830 kHz, but the signal-to-noise ratio at low is too low to observe two additional peaks, which are expected close to 500 kHz. Although the peaks are well located in the spectrum, the resolution is relatively low due to the Zeeman-powder broadening of the quadrupole peaks. Additionally, the CR peaks are Zeeman-shifted with respect to the pure NQR frequencies in the zero field [41].
In the second step we, therefore, used the techniques of multiple frequency sweeps and two-frequency irradiation [42,43] for a more precise determination of the frequencies of the NQR triplets ν+, ν-, and ν0 at the three nitrogen positions in an FVP molecule. In the search for the 14 N NQR frequencies at the pyrazine nitrogen positions, we applied multiple frequency sweeps of the rf magnetic field in the frequency range 2.8-4 MHz and performed the νH scan. We clearly observed two dips in the νL dependence of the magnitude of the proton NMR signal at the end of the magnetic field cycle at the frequencies of 620 kHz and 900 kHz. These are the lowest frequencies ν0 at the two pyrazine nitrogen positions. Then, we fixed the low magnetic field within the dip at νH = 620 kHz and selectively varied the upper and the lower-frequency limit. When the upper-frequency limit passes 3570 kHz from above, the magnitude of the proton NMR signal increases. Similarly, when the lower-frequency limit passes 2.950 MHz from below, the proton NMR signal increases. Therefore, the frequencies of 3.570 MHz, 2.950 MHz, and 620 kHz are within the experimental resolution of this technique, i.e., the 14 N NQR frequencies at the one pyrazine nitrogen position. We tried to refine these values with the application of the two-frequency irradiation technique; however, within the experimental resolution of this technique, ±5 kHz, we obtained the same results. Similarly, we determined the 14 N NQR frequencies at the second pyrazine position and the amide nitrogen position. The obtained results are presented in Table 1. Three strong peaks centred at the frequencies of 3.870 MHz, 3.610 MHz, and 2.940 MHz are clearly observed. These peaks are observed in the frequency range in which the higher 14 N NQR frequencies ν + and ν − of pyrazine nitrogen are usually observed. Since there are two non-equivalent pyrazine nitrogen atoms in the molecule, we expect to observe four peaks. It is possible that two peaks overlap within the experimental resolution.
Two somewhat weaker peaks are observed at the frequencies of 2.140 MHz and 1.610 MHz. In this frequency range, the amide nitrogen NQR frequencies ν + and ν − are usually observed. Below 1 MHz, there is a well-observed peak at 830 kHz, but the signal-tonoise ratio at low ν H is too low to observe two additional peaks, which are expected close to 500 kHz. Although the peaks are well located in the spectrum, the resolution is relatively low due to the Zeeman-powder broadening of the quadrupole peaks. Additionally, the CR peaks are Zeeman-shifted with respect to the pure NQR frequencies in the zero field [41].
In the second step we, therefore, used the techniques of multiple frequency sweeps and two-frequency irradiation [42,43] for a more precise determination of the frequencies of the NQR triplets ν + , ν − , and ν 0 at the three nitrogen positions in an FVP molecule. In the search for the 14 N NQR frequencies at the pyrazine nitrogen positions, we applied multiple frequency sweeps of the rf magnetic field in the frequency range 2.8-4 MHz and performed the ν H scan. We clearly observed two dips in the ν L dependence of the magnitude of the proton NMR signal at the end of the magnetic field cycle at the frequencies of 620 kHz and 900 kHz. These are the lowest frequencies ν 0 at the two pyrazine nitrogen positions. Then, we fixed the low magnetic field within the dip at ν H = 620 kHz and selectively varied the upper and the lower-frequency limit. When the upper-frequency limit passes 3570 kHz from above, the magnitude of the proton NMR signal increases. Similarly, when the lower-frequency limit passes 2.950 MHz from below, the proton NMR signal increases. Therefore, the frequencies of 3.570 MHz, 2.950 MHz, and 620 kHz are within the experimental resolution of this technique, i.e., the 14 N NQR frequencies at the one pyrazine nitrogen position. We tried to refine these values with the application of the two-frequency irradiation technique; however, within the experimental resolution of this technique, ±5 kHz, we obtained the same results. Similarly, we determined the 14 N NQR frequencies at the second pyrazine position and the amide nitrogen position. The obtained results are presented in Table 1. The experimental 1 H-14 N NQDR spectrum of FVP consists of nine signals, Figure 2, which must be grouped into three groups for each of the three 14 N sites, as shown in Table 1. The NQR parameters e 2 Qqh −1 and η, calculated from these data using Equations (2) and (3), are listed in Table 1. The number of resonance lines indicates that the molecules in the FPV crystal unit cell are equivalent and symmetry-related. Moreover, the NQR spectrum of FVP is normally broadened, with no evidence of structural disorder.
The electronic environment of all three nitrogen nuclei in the FVP molecule is different, so they are chemically non-equivalent. All nitrogen sites have distinctly different e 2 Qqh −1 but very close η; the latter is rather unusual. The low ν + and e 2 Qqh −1 as well as the characteristic resonant frequencies observed in the spectrum are typical of amide nitrogen (-NH 2 ). The ν + frequency of -NH 2 is substantially lower than that of aniline (3.244 MHz) and close to those in acetamide (2.105 MHz) or formamide (1.920 MHz) [45], which suggests a strong conjugation between the amide nitrogen and the adjacent carbon in FVP. Additionally, the value of η calculated from the 14 N NQR spectrum, closer to that of -N= than the -NH 2 type of nitrogen, suggests that the electron density at -NH 2 is significantly shifted due to the involvement of this moiety in the intermolecular interactions. Both NQR parameters in FVP are not only close to those observed for acetamide or formamide having hydrogen bonds N-H···O but also closer to that of cytosine [46] having two kinds of hydrogen bonds N-H···O and N-H···N [30], than to that of aminopyridines with only weaker N-H···N bonds. It points to the participation of -NH 2 in FVP in at least two kinds of hydrogen bonds, most likely N-H···O and N-H···N. Thus, a basic analysis of the NQR spectrum suggests that the set of resonant lines resulting in the lowest e 2 Qqh −1 can be assigned to the amide nitrogen -NH 2 , while the sets resulting in the two higher ones can be assigned to the pyrazine nitrogen atoms -N(1)= and -N(4)=. An unambiguous assignment of pyrazine nitrogen signals based solely on the NQR spectra is not possible.
The substituents effect, analysed directly using the NQR parameters for the structurally related compounds or modelled using the calculations for the single molecules, should be helpful in the assignment of the pyrazine signals. Due to the lack of experimental data (NQR spectra are available only for unsubstituted pyrazine [44]), the substituent effect could be only modelled. The results are summarized in Table 2 and shown in Figure 3.   The fluorine and hydroxyl groups substituted directly into the pyrazine ring at the C(6) and C(3) positions, respectively, are both σ electron-withdrawing and π-electron-donating. However, due to their symmetrical position in relation to the pyrazine nitrogen atoms, they have a competitive effect and act very similarly on the -N(1)= and -N(4)= nitrogen atoms. As a consequence, the shifts in the NQR frequencies and the values of parameters for both -N(1)= and -N(4)= differ only slightly, as shown in Table 2. The amide group substituted at the C(2) position of the ring also acts as an σ electron-withdrawing group, and its nitrogen is an excellent π-donor primarily to the carboxylic carbon. Therefore, it affects the electron density distribution in the ring but has a slightly stronger influence on -N(4)= in the meta position, than -N(1)= in the ortho position. Surprisingly, the amide group is quite an important factor in determining the position of the resonance lines in the spectrum, as shown in Figure 3. It means that the effect of the resultant substituent in the FVP molecule is complicated by a combination of nonadditive overlap and drift effects, resulting in a significant decrease in both NQR parameters (e 2 Qqh −1 and η) on both pyrazine nitrogen atoms. The single-molecule calculations, performed for the optimized and "solid" molecule (i.e., cut out of the crystal), reproduce this effect only roughly. The Pearson correlation coefficient is not very high, the scattering is quite large, and the NQR frequencies derived from these parameters are overestimated. This is a result of neglecting the further surroundings, in particular the nearest neighbouring molecules. The simulation of the 14 N NQR spectra of the other heterocyclic compounds showed that NQR parameters were very sensitive, even to the influence of the distant proximity of nitrogen atoms [29,38]. This raises a suspicion that despite reasonable values of NQR parameters, the assignment of the signals based on single-molecule calculations in FVP may be questionable.
The solid-state calculations that take into account all intermolecular interactions require a high-quality reliable crystal structure. The Cambridge structural database (CSD) contains four FVP crystalline structures for two different polymorphs (three polymorph I and one polymorph II; all enol-like forms). Polymorph I of FPV crystallizes in the orthorhombic space group Pna2 1 , with four equivalent molecules in the unit cell whose dimensions vary with temperature from a= 9.06680 Å, b = 14.85080 Å, and c = 4.57550 Å at 100 K (CCDC 2048895 [47]) or a = 9.14400 Å, b = 14.87000 Å, and c = 4.55730 Å at 100 K (CCDC 2106318 [48]), to a = 9.11060 Å, b = 14.76190 Å, and c = 4.69100 Å at 296 K (CCDC 969968 [49]). Polymorph II of FPV crystallizes in the tetragonal space group P42/n with four equivalent molecules in the unit cell of dimensions a = 9.28050 Å, b = 9.28050 Å, and c= 14.64910 Å at 298 K (CCDC 2047143 [50]). The 14 N NQDR spectrum indicates the presence of only one polymorphic form in the sample; however, it is not possible to infer from it which one exactly. Therefore, all of these structures were used as inputs for the DFT/GGA calculations. The results of the calculations assuming the structures CCDC 2048895, 2106318, and 969968 as inputs, are significantly better than those obtained assuming CCDC 2047143 (the Pearson correlation is higher and the scattering is smaller), as shown in Tables 3 and 4, and Figure 4. Thus, the simulated spectrum, much closer to the experimental one, clearly indicates the presence of polymorph I in the studied sample.   The differences between the experimental and calculated NQR parameters are rather small for the pyrazine nitrogen but quite large for the amide nitrogen, regardless of the level of theory and type of structure assumed. While e 2 Qq/h is quite well reproduced, the asymmetry parameter is significantly, almost twice, overestimated. This leads to the res- The differences between the experimental and calculated NQR parameters are rather small for the pyrazine nitrogen but quite large for the amide nitrogen, regardless of the level of theory and type of structure assumed. While e 2 Qq/h is quite well reproduced, the asymmetry parameter is significantly, almost twice, overestimated. This leads to the resonance frequencies ν − and ν 0 for -NH 2 calculated from the NQR parameters and Equations (4)-(6) being far from the experimental values. The ν + frequency for -N(4)= also deviates from the experimental value. There are two main reasons for this phenomenon: one is related to the quality of the crystallographic structure and associated with thermal effects, and the other is related to the nature of intermolecular interactions. X-ray crystallography at the usual resolution does not permit the exact assessment of the positions of light atoms, mainly because their contributions to diffraction are weak. In addition, the atoms oscillate around their equilibrium positions and the amplitude of the oscillations strongly depends on temperature (the lower the temperature, the smaller the oscillations).
At higher temperatures, the electron density is spread out over a larger volume, resulting in a more rapid decrease in atomic scattering power and, consequently, all of the atomic positions, bond lengths, and bond angles are not sufficiently well determined. Actually, the results closer to the experimental ones were obtained assuming the X-ray structure was determined at 100 K rather than by assuming the X-ray structure at 296 K, as seen in Table 3, although the NQR experiment was performed at 295 K. Moreover, the NQR parameters were reproduced significantly better under the assumption of the optimized structure, which minimizes the influence of the aforementioned error (Table 5). Optimization of the low-temperature structures (CCDC 2048895 or 2106318) with fixed unit cell parameters and at the GGA/RBPE level of theory yields more reasonable C-H, O-H, and N-H bond lengths, and amide nitrogen hybridization is a bit closer to the real value (sp 2 ). As a consequence, even the asymmetry parameter for -NH 2 is closer to the experimental value, although still slightly overestimated (by about 20%).
The cluster calculations that take into account all intermolecular interactions only in the immediate vicinity of the molecule, performed with noncovalent forces dedicated to hybrid functional M062X, gave better results. The spread for all nuclei, including -NH 2 , was smaller; however, the NQR parameters were overestimated by an average of 24% and, consequently, the resonance lines were upshifted in Table 6. Such calculations confirmed the validity of the assignment made based on the solid-state calculations and delivered a high-quality function for further QTAIM analysis. It is noteworthy that the assignment of signals to pyrazine nitrogen atoms obtained from the cluster and solid-state calculations are the opposite of those obtained from the single-molecule calculations (Tables 3-6). Moreover, the same conclusion was drawn regardless of the level of calculations as well as the assumption of each of the X-ray or optimized structures. The detailed analysis of this effect reveals a significant difference between both pyrazine nitrogen atoms. The e 2 Qq/h value of the -N(1)= nitrogen atom increases, while that of the -N(4)= nitrogen decreases, and the decrease in η of -N(4)= is more significant than that of -N(1)=. Moreover, a comparison of the calculated NQR frequencies shows that regardless of the assumed geometry (X-ray or optimized), the changes in the ν + of pyrazine nitrogen -N(4)= and ν 0 and ν − of amide nitrogen -NH 2 are the highest. Since ν 0 and ν − are highly sensitive to hydrogen bonding, the scattering of the two latter clearly indicates the involvement of the -NH 2 group in such bonds. The map of the static deformation density, Figure 5, reveals the redistribution of valence electron density near nitrogen sites caused by both chemical bonding and intermolecular interactions.
The significant charge rearrangement in the solid is concentrated on the atomic nuclei rather than in the internuclear region. The areas of charge accumulation in the bonding regions and the lone-pair positions of the atoms are visible. Upon the N-H···N and N-H···O hydrogen bond formation, the electron density is shifted from the lone-pair orbital of -N= (and =O) to the N-H bond of the acceptor. This effect is reflected in an increase in ν + , ν −, e 2 Qqh −1 , and η of -N(4)=, which is the acceptor, as well as a decrease in ν − , e 2 Qqh −1 , and an increase in ν + , ν 0 , and η of -NH 2 , which is the donor. Upon the hydrogen bonding formation, the population of the π orbital of -NH 2 , which is normal to the planar arrangement of N-C and N-H bonds, increases when compared to the populations of the σ orbitals (σ NC and σ NH ). Therefore, the Z-axis of the EFG tensor is oriented along the π orbital and the X and Y-axes of the EFG tensor are rotated around the Z-axis ( Figure 6). creases, while that of the -N(4)= nitrogen decreases, and the decrease in η of -N(4)= is more significant than that of -N(1)=. Moreover, a comparison of the calculated NQR frequencies shows that regardless of the assumed geometry (X-ray or optimized), the changes in the + of pyrazine nitrogen -N(4)= and 0 and of amide nitrogen -NH2 are the highest. Since 0 and are highly sensitive to hydrogen bonding, the scattering of the two latter clearly indicates the involvement of the -NH2 group in such bonds. The map of the static deformation density, Figure 5, reveals the redistribution of valence electron density near nitrogen sites caused by both chemical bonding and intermolecular interactions.  The significant charge rearrangement in the solid is concentrated on the atomic nuclei rather than in the internuclear region. The areas of charge accumulation in the bonding regions and the lone-pair positions of the atoms are visible. Upon the N-H···N and N-H···O hydrogen bond formation, the electron density is shifted from the lone-pair orbital of -N= (and =O) to the N-H bond of the acceptor. This effect is reflected in an increase in  + , −, e 2 Qqh −1 , and η of -N(4) = , which is the acceptor, as well as a decrease in -, e 2 Qqh −1 , and an increase in  + , 0, and η of -NH2, which is the donor. Upon the hydrogen bonding formation, the population of the π orbital of -NH2, which is normal to the planar arrangement of N-C and N-H bonds, increases when compared to the populations of the σ orbitals (σNC and σNH). Therefore, the Z-axis of the EFG tensor is oriented along the π orbital and the X and Y-axes of the EFG tensor are rotated around the Z-axis ( Figure 6).  An increase in σ NH , caused by charge transfer from the electron pair of -NH 2 , results in a decrease in e 2 Qq/h. Because σ NH > σ NC , the deviation from the axial symmetry along the Z-axis increases, and as a result, η also increases. The fact that e 2 Qq/h is considerably lower for the amide nitrogen in the solid FVP than in a single molecule of FVP supports the above reasoning and is in good agreement with the results obtained, for example, for imidazole. The NQR parameters of -N(4)= are related to the populations of the lone-pair orbital σ LP and the nitrogen orbitals involved in the σ NC and π bonds to carbon. Upon hydrogen bond formation, the X, Y and Z-axes of the EFG tensor at -N(4)= keep their orientation unchanged ( Figure 6). However, a decrease in σ LP of -N(4)= results in an increase in e 2 Qq/h and η. The frequencies ν 0 , ν + , and ν − are proportional to the populations of the nitrogen atom bonds and both NQR parameters; therefore, this effect results in frequency shifts. The solid-state shift of the NQR parameters for nitrogen -N(1)= is small: η is well reproduced, while e 2 Qq/h is only slightly overestimated. The latter suggests that -N(1)= participates in intermolecular interactions, but not in strong ones. Upon hydrogen bond formation, the X, Y and Z-axes of the EFG tensor at -N(1)= keep their orientation unchanged ( Figure 6), which supports this conclusion. However, a possible competition between acceptors (nitrogen, oxygen, and fluorine) raises the question of whether the experimentally observed local nitrogen environment corresponds to the most likely hydrogen bonds, and if the -NH and -OH groups can form intra-and intermolecular interactions.

Interactions Pattern
For a detailed evaluation of the pattern of intermolecular interactions in FVP crystals, we applied 3D HS, 2D molecular fingerprints, and an ESP approach, which helps to identify their type and the QTAIM with RDS for their strength estimation.

Hirshfeld Surfaces
A combined analysis of 3D Hirshfeld Surfaces (3D HS) with d norm , shape index, or curvedness mapped over these surfaces and 2D molecular fingerprints (2D FP) derived from 3D HS, which summarizes the distribution of interactions of a molecule with its environment, delivers an insight into the mix of the 15 types of the homo-and heteronuclear intermolecular interactions in the FVP crystal. As much as 40% of the FVP molecular surface is generated by H atoms (only four), while the other contributing elements C, F, N, and O generate much smaller but quite similar percentages of the molecular surface-from about 13 to 18 % (Table 7 and Figure 7).  Approximately 81% of all interactions of the complex surface consist of seven types O···H, F···H, N···H, H···H, C···H, C···N, and C···O. There are only small differences in the percentages of these interactions for different structures of polymorph I (X-ray, proton, and fully optimized), as shown in Figure 7.  Approximately 81% of all interactions of the complex surface consist of seven types O···H, F···H, N···H, H···H, C···H, C···N, and C···O. There are only small differences in the percentages of these interactions for different structures of polymorph I (X-ray, proton, and fully optimized), as shown in Figure 7.
The distribution of interactions is weakly affected by temperature. The whole distribution is, to a small extent, changed by optimization; for heavy atoms within the error limit, by 0.1 %, only for protons, slightly more, by 0.9 %. Thus, the only significant changes after optimization are the positions of the protons, resulting in shorter and less nonlinear and, thus, stronger hydrogen bonds.
The enrichment ratio, E XY , of the main intermolecular contacts, which reveals privileged (E XY > 1) and disfavoured (E XY < 1) contacts between every two atomic species, X and Y, are collected in Table 8.
For all examined FVP crystalline structures of polymorph I (X-ray, proton, and fully optimized), the percentages of the molecular surface and enrichment ratios are very close (the differences do not exceed 2%) (Figure 8). Four C and two N atoms in the six-member heterocyclic pyrazine ring are exposed to the molecular surface because each of them is three (C) and two (N) bonded (trigonal planar sp 2 hybridization), and their p orbitals form π bonds above and below the ring. The C···N contacts in FVP are significantly enriched (ECN  1.9), thus, the C species are the preferred interaction partners for N and more favoured compared to N···N species (ENN  1.5). However, both contacts have enrichment ratios much higher than unity due to the  Four C and two N atoms in the six-member heterocyclic pyrazine ring are exposed to the molecular surface because each of them is three (C) and two (N) bonded (trigonal planar sp 2 hybridization), and their p orbitals form π bonds above and below the ring. The C···N contacts in FVP are significantly enriched (E CN ≈ 1.9), thus, the C species are the preferred interaction partners for N and more favoured compared to N···N species (E NN ≈ 1.5). However, both contacts have enrichment ratios much higher than unity due to the high nitrogen but low hydrogen abundance in the heterocyclic ring and the high propensity of π stacking with polarised atoms as interacting partners. This suggests that in π· · · π stacking, the C(5)···N(NH 2 ), C(3)···N(1), C(6)···N(4) as well as N(1)···N(4) contacts play an important role in the crystalline structure. For the O···C contacts, the enrichment ratio is also high (E OC ≈ 1.5), which suggests that in π· · · π stacking, the O atoms interact with positively charged C atoms, which are themselves covalently bonded to O atoms. Thus, the possibility of the contacts C(6)···O(OH) and C(CONH 2 )···O(OH) is raised. Among C···C and H···C interactions, which compete in the crystal, only C· · · H contacts are favoured, but their enrichment ratio E CH < 1 as H atoms are involved in hydrogen bonds. Concurrently, the C···C contacts are significantly impoverished (E CC ≈ 0), so their participation in π···π stacking is electrostatically unfavourable. A very high E CH /E CC suggests the classification of the packing motifs as herringbone [51]. The propensity of F···F and F···H/H···F contacts is significantly enriched (E FF ≈ 2.0 and E FH ≈ 1.4), thus, the F and H species are the interaction partners for F. Among them, the F···F contacts, which form zig-zag chains, are preferred. However, the F···H/H···F contacts are also favoured in the crystal packing and raise the probability of F···HX hydrogen bonds (F···H-C(3) and F···H-O(OH)). The H···H contacts appear with an enrichment ratio significantly smaller than unity (E HH ≈ 0.8) due to their strong competition with both H···F contacts (E FH ≈ 1.4) and the hydrogen acceptors H···N (E NH ≈ 0.8) and H···O (E OH ≈ 1.3), additionally forced by the presence of low hydrogen content. The O···H/H···O contacts show enrichment values (E OH ≈ 1.3) because they are favoured in the crystal packing, and, notably, the O=C oxygen atoms form intramolecular hydrogen bonds of the type O-H···O (closing the five-member ring) and intermolecular bonds of the type N-H···O (-NH from -NH 2 ). For the same reason, the O···O contacts are significantly impoverished (E OO ≈ 0.1). The N···H/H···N contacts show enrichment values smaller than unity (E OH ≈ 0.9), and, thus, are disfavoured, and the existence of the N-H···N hydrogen bonds linking -N(4)= and -NH 2 is not raised. In general, the favoured interaction partners for N, whose contribution to the 3D HS molecular surface is about 16%, are C and N species, while H, followed by O and F are disfavoured.
A little more information about the types of interactions and their location is provided by 3D HS, which, with the normalized contact distance d norm , shape index and curvedness mapped over this surface, are shown in Figure 9. interaction partners for N, whose contribution to the 3D HS molecular surface is about 16%, are C and N species, while H, followed by O and F are disfavoured.
A little more information about the types of interactions and their location is provided by 3D HS, which, with the normalized contact distance dnorm, shape index and curvedness mapped over this surface, are shown in Figure 9.  Optimization of the geometry with fixed parameters of the unit cell results in the shortening of hydrogen bonds in relation to crystallographic bonds, which affects the 3D crystal pattern. This effect is poorly reflected by the shape of the 3D HS surface but well described by d norm , shape index, and curvedness. The N-H···O is described by values of d norm = −0.6313 a.u., shape index = −0.9428 a.u., and curvedness = −3.3128 a.u., while the N-H···N is characterized by values of d norm = −0.4436 a.u., shape index = −0.9957 a.u., and curvedness = −1.4333 a.u. Both are shorter (R O···N = 2.796 Å and R N···N = 2.915 Å) and more linear (<NHO = 170.87 • and <NHN = 142.20 • ), and, thus, stronger. The C-H···O hydrogen bond, which is described by values of d norm = −0.1367 a.u., shape index = 0.8807 a.u., and curvedness = −1.8613 a.u. is longer (R C···O = 3.542 Å) and weaker. The C(5) -H···F hydrogen bond of length 3.542 Å is characterized by values of d norm = 0.2989 a.u., shape index = 0.7926 a.u., and curvedness = −1.5158 a.u., i.e., weaker. The N···C contacts of 3.351-3.411 Å are described by values of d norm = 0.068 a.u., shape index = −0.629 a.u., and curvedness = −1.738 a.u.
In order to verify the presence of interactions, mainly directional interactions, the characteristic features of the 2D fingerprint (2D FP) de/di plots (Figure 10), obtained from the 3D HS, were further analysed. In the 2D FP plot, the O···H/H···O contacts describing intramolecular O-H···O and intermolecular N-H···O hydrogen bonds are represented by symmetric spikes (''wings''), which are sharp and small, at de + di ~ 2.0 Å (de + di ~ 1.75 Å for the optimized geometry). These interactions cover a small area of 17.8% (17% for the optimized geometry) of the In the 2D FP plot, the O···H/H···O contacts describing intramolecular O-H···O and intermolecular N-H···O hydrogen bonds are represented by symmetric spikes ("wings"), which are sharp and small, at d e + d i~2 .0 Å (d e + d i~1 .75 Å for the optimized geometry). These interactions cover a small area of 17.8% (17% for the optimized geometry) of the total 3D HS. The N···H/H···N contacts are represented by the most external symmetric spikes sharp, and d e + d i~2 .2 Å (d e + d i~2 .0 Å for the optimized geometry). These interactions cover a small area of 11.5% (11.7% for the optimized geometry) of the total 3D HS. The F···H/H···F contacts are represented in 2D FP by the widest spikes at d e + d i~3 .00 Å. They cover 15.6% (14.9% for the optimized geometry) of the total 3D HS. The C···O/O···C contacts bring a contribution of 7.1% (7.7% for the optimized geometry) and they are represented by sharp wings, d e + d i~3 .00 Å (3.10 Å for the optimized geometry), located in the middle area of the whole fingerprint. The O···N, O···F, and F···F contacts bring a small contribution of 4.8, 4.2, and 3.4% (5.0, 4.3, and 4.1% for the optimized geometry), respectively. The relatively high contribution is brought by weak H···H interactions, which cover 12.0 % (12.1% for the optimized geometry) of the total 3D HS. In the 2D FP, they are reflected by the cloud of scattered points. The contribution brought by C···N/N···C and N· · · N contacts, which cover 8.1 and 3.8 % (8.1 and 3.9% for the optimized geometry), respectively, mainly comes from C(3)· · · N(4) and N(1)···N(4) contacts, i.e., interlayer π···π stacking interactions.
The local 2D molecular fingerprints of -N(1)=, -N(4)=, and -NH 2 =, generated based on 3D HS for the particular nitrogen atoms (Figure 11), vary significantly from one nitrogen site to another, in both shape and percentage. While the upper narrow wing attributed to the NC covalent bond is common and featureless, the lower one matters because it shows a considerable variation in the number (three in -N(4)=, four in -N(1)=, and two in -NH2) and type of interactions.
The local 2D FPs reveal the participation of -NH2 in two covalent NH bonds being While the upper narrow wing attributed to the NC covalent bond is common and featureless, the lower one matters because it shows a considerable variation in the number (three in -N(4)=, four in -N(1)=, and two in -NH 2 ) and type of interactions.
The local 2D FPs reveal the participation of -NH 2 in two covalent NH bonds being donors for two strong hydrogen bonds NH···O and NH···N and the participation of -N(4)= in intermolecular N-H···N hydrogen bond and N···N, N···O contacts, Table 9. Moreover, it suggests the involvement of -N(1)= in the intramolecular hydrogen bond N-H···N (which closes the five-member ring but was not detected in the global 3D HSbased analyses) and in three N···N, N···O, and N···F contacts. The wing describing the intramolecular hydrogen bond N-H···N, very wide and scattered at d e + d i =2.4 Å, is highly different from the intermolecular NH···N, which is very narrow and at d e + d i = 2.15 Å. A similar analysis of the local 2D FP for the oxygen atoms reveals a sharp upper wing attributed to the OC covalent bond at d e + d i = 1.3 Å and a bottom wider wing attributed to hydrogen bonds (intermolecular O-H···N at d e + d i = 1.65 Å for =O and intramolecular O-H···O at d e + d i = 2.15 Å for -OH) and some heavy-atom contacts, Table 9.
The wing describing the intramolecular hydrogen bond O-H···O is typical of this kind of bond. It is worth noting that N-H···N has a very similar shape of a wide and scattered cloud of points. The local 2D FP of the fluorine atom reveals a sharp upper wing attributed to the FC covalent bond and a wide bottom one indicating the involvement of fluorine in the F···H-C hydrogen bond and three kinds of contacts F···N, F···O, and F···F (Table 9).
Thus, the analysis of local, atomic, 2D molecular fingerprints helps to discover/confirm the presence of the very weak hydrogen bonds N···H-N (intramolecular) and F···H-C (intermolecular) which meet the Gilli criterion [52] and are nicely complementary to the total 3D HS and 2D FP. Moreover, they explain the source of the differences between particular nitrogen nuclei detected by NQR.

Characterization of the Strength of the Interactions The Total Interaction Energy Partitioning
The analysis of the total lattice energy and its terms (electrostatic, polarization, dispersion, and repulsion) derived from DFT calculations for a cluster of 10 Å in diameter has revealed that irrespective of temperature, the major contributions come from the electrostatic or dispersion terms and could not be cancelled by the repulsion. The total lattice energy (−105.6 kJ/mol) is dominated by the dispersion (−77.17 kJ/mol) and electrostatic (−70.82 kJ/mol) terms, which are stronger than the polarization (−13.76 kJ/mol) or repulsion terms (56.79 kJ/mol). Upon a temperature increase, all negative terms increase dispersion (−72.82 kJ/mol), electrostatic (−63.21 kJ/mol), and polarization (−12.14 kJ/mol), while repulsion decreases (48.32 kJ/mol), which results in slightly lower total energy (−100.5 kJ/mol). The optimization with the fixed unit cell results in an increase in each term of energy (electrostatic (−78.32 kJ/mol), dispersion (−77.95 kJ/mol), polarization (−15.98 kJ/mol), and repulsion (65.26 kJ/mol), which in turn causes a decrease in the total energy (−95.2 kJ/mol). The total interaction energy partitioning indicates that the crystalline packing is determined by mainly electrostatic and repulsive strong N-H· · · O hydrogen bonds followed by N−H···N and weaker C-H· · · O, and dispersive strong π· · · π stacking and relatively weak F· · · F contacts (Table 10 and Figure 12). (However, a partial contribution to the repulsion term in N-H· · · N may, in fact, come from N· · · C contacts). H···O at de + di = 2.15 Å for -OH) and some heavy-atom contacts, Table 9. The wing describing the intramolecular hydrogen bond O-H···O is typical of this kind of bond. It is worth noting that N-H···N has a very similar shape of a wide and scattered cloud of points. The local 2D FP of the fluorine atom reveals a sharp upper wing attributed to the FC covalent bond and a wide bottom one indicating the involvement of fluorine in the F···H-C hydrogen bond and three kinds of contacts F···N, F···O, and F···F (Table 9).
Thus, the analysis of local, atomic, 2D molecular fingerprints helps to discover/confirm the presence of the very weak hydrogen bonds N···H-N (intramolecular) and F···H-C (intermolecular) which meet the Gilli criterion [52] and are nicely complementary to the total 3D HS and 2D FP. Moreover, they explain the source of the differences between particular nitrogen nuclei detected by NQR.

Characterization of the Strength of the Interactions
The Total Interaction Energy Partitioning The analysis of the total lattice energy and its terms (electrostatic, polarization, dispersion, and repulsion) derived from DFT calculations for a cluster of 10 Å in diameter has revealed that irrespective of temperature, the major contributions come from the electrostatic or dispersion terms and could not be cancelled by the repulsion. The total lattice energy (−105.6 kJ/mol) is dominated by the dispersion (−77.17 kJ/mol) and electrostatic (−70.82 kJ/mol) terms, which are stronger than the polarization (−13.76 kJ/mol) or repulsion terms (56.79 kJ/mol). Upon a temperature increase, all negative terms increase dispersion (−72.82 kJ/mol), electrostatic (−63.21 kJ/mol), and polarization (−12.14 kJ/mol), while repulsion decreases (48.32 kJ/mol), which results in slightly lower total energy (−100.5 kJ/mol). The optimization with the fixed unit cell results in an increase in each term of energy (electrostatic (−78.32 kJ/mol), dispersion (−77.95 kJ/mol), polarization (−15.98 kJ/mol), and repulsion (65.26 kJ/mol), which in turn causes a decrease in the total energy (−95.2 kJ/mol). The total interaction energy partitioning indicates that the crystalline packing is determined by mainly electrostatic and repulsive strong N-HO hydrogen bonds followed by N−H···N and weaker C-HO, and dispersive strong ππ stacking and relatively weak FF contacts (Table 10 and Figure 12). (However, a partial contribution to the repulsion term in N-HN may, in fact, come from NC contacts).
Molecules 2023, 28, x FOR PEER REVIEW 24 of 39 Figure 12. Energy framework of diagrams for electrostatic, dispersion, and total energy visualised in the FVP structure by red, green, and blue "sticks", respectively (view along the a-axis-top, along the b-axis-middle, and along the c-axis-bottom).  The ordering of the most important interactions according to the decreasing energy is as follows: N-H· · · O ≈ π· · · π stacking (parallelly oriented molecules) > >N-H· · · N ≈ π· · · π stacking (flipped molecules) > C-H· · · O >> F· · · F. The N-H· · · O hydrogen bonds and π· · · π stacking (slightly different between parallel or flipped molecules) forcing planar arrangement and F· · · F contacts forming zig-zag chains, are the main driving forces acting cooperatively (Figure 12).

Electrostatic Potential
The electrostatic nature of the hydrogen bonds, that is, the match of positive and negative regions of the molecular surfaces (electrostatic complementarity), is also well pronounced in the molecular electrostatic potential (ESP). The positive (blue) regions of ESP are clearly separated from the negative (red) regions in Figure 13. The N-HO hydrogen bonds and ππ stacking (slightly different between parallel or flipped molecules) forcing planar arrangement and FF contacts forming zig-zag chains, are the main driving forces acting cooperatively (Figure 12).

Electrostatic Potential
The electrostatic nature of the hydrogen bonds, that is, the match of positive and negative regions of the molecular surfaces (electrostatic complementarity), is also well pronounced in the molecular electrostatic potential (ESP). The positive (blue) regions of ESP are clearly separated from the negative (red) regions in Figure 13. The negative ESP on both sides of the molecule occupies O and double O infinity sign and hood-shaped regions. The highly positive region is localised on almost the whole moiety, including -NH2 and -N(1)=, while the negative one is only on nitrogen -N(4)= and both oxygens. The arrangement of the positive and negative ESP sites helps to understand hydrogen bond formation. The negatively charged regions located at -N(4)=, -OH, and =O are linked by bonds to the positively charged regions located at -NH2, -N(1)=, and -C(5)H (they perfectly match each other). As a result, N-HO and C-HO hydrogen bonds and N(4)O contacts are formed. Moreover, the N-HN hydrogen bond, although dominated by electrostatic forces, is not pronounced in the ESP due to the lack of electronegativity difference between the acceptor and donor. These observations are extremely important from the docking point of view as -N(4)= plays a key role in the prodrug conversion to the active substance (phosphoribosylation).

Quantum Theory of Atoms in Molecules
The quantum theory of atoms in molecules (QTAIM) [31] approach was used for the verification and strength estimation of the inter-and intramolecular interactions revealed by 3D HS and 2D FP. The topological descriptors of electron density (ρ(r) and ρ(r), the HBCP(r) and its components GBCP(r) and VBCP(r)), as well as the strength of all intra-and The negative ESP on both sides of the molecule occupies O and double O infinity sign and hood-shaped regions. The highly positive region is localised on almost the whole moiety, including -NH 2 and -N(1)=, while the negative one is only on nitrogen -N(4)= and both oxygens. The arrangement of the positive and negative ESP sites helps to understand hydrogen bond formation. The negatively charged regions located at -N(4)=, -OH, and =O are linked by bonds to the positively charged regions located at -NH 2 , -N(1)=, and -C(5)H (they perfectly match each other). As a result, N-H· · · O and C-H· · · O hydrogen bonds and N(4)· · · O contacts are formed. Moreover, the N-H· · · N hydrogen bond, although dominated by electrostatic forces, is not pronounced in the ESP due to the lack of electronegativity difference between the acceptor and donor. These observations are extremely important from the docking point of view as -N(4)= plays a key role in the prodrug conversion to the active substance (phosphoribosylation).

Quantum Theory of Atoms in Molecules
The quantum theory of atoms in molecules (QTAIM) [31] approach was used for the verification and strength estimation of the inter-and intramolecular interactions revealed by 3D HS and 2D FP. The topological descriptors of electron density (ρ(r) and ∆ρ(r), the H BCP (r) and its components G BCP (r) and V BCP (r)), as well as the strength of all intra-and intermolecular bonds and contacts, are collected in Table 11. In general, QTAIM is, in some cases, not sufficient to detect weak noncovalent interactions [32], thus the isosurface of the reduced density gradient RDG(r) was examined, and the corresponding plot was generated for the FVP structure. The molecular graphs cut off the crystalline structure with BCP, RCP, and reduced density gradient isosurfaces (RDG) with sign(λ 2 )ρ BCP mapped over the surface, as shown in Figure 14. Table 11. Topological parameters describing the interactions in FVP (electron density at the bond critical point, BCP, and (ρ BCP (r)), its Laplacian (∆ρ BCP (r)), the potential electron energy density (V BCP (r)), the kinetic electron energy density (G BCP (r)), the total electron energy density (H BCP (r)), the bonding energy according to Espinosa (E E ), Matta (E M ), Emamian (E EM ), Afonin (E A ), Gilli (E G ) and Nikolaienko (E N ) calculated at the M062X/6-311 + G** level.

Interaction
Nitrogen Site  Results of QTAIM confirmed the presence of the hydrogen bonds and contacts detected by the total 2D FP. The energies estimated using QTAIM and Equations (7)-(12), Table 11, are in good agreement with those estimated using the energy framework. In contrast to the local 2D FP, QTAIM did not detect intramolecular N-H···N hydrogen bonds. However, the disc-shaped green surface of the reduced density gradient located between the -N(1)= and H (from -NH 2 ) indicated the presence of a weak intramolecular N-H···N hydrogen bond, although no critical points were detected ( Figure 14). In the plot of RDG(r) versus sign(λ 2 )ρ(r), Figure 15, the noncovalent interactions-hydrogen bonds (strong and attractive), van der Waals contacts (weak), and steric interactions (ring closures, repulsive)-are visible as spikes.
Their positions were identified based on the analysis of the closely related pyrazine structures. The intramolecular N-H···N hydrogen bond spike is shown in light blue. The analysis of the bonds in which each of the nitrogen atoms participates within QTAIM and RDG reveals significant disproportions between them. The -NH 2 site participates in two intermolecular hydrogen bonds of different strengths N-H···N and N-H···O; -N(4)= participates in N-H···N intermolecular hydrogen bond, but -N(1)= only in a very weak putative intramolecular hydrogen bond N-H···N. This is reflected in the asymmetry parameters, the highest for asymmetrically bonded -NH 2 , and the lowest for -N(4)=, participating only in one bond coaxially with the Z-axis of the EFG tensor. The position of the spike in the plot of RDG(r) versus sign(λ 2 )ρ(r), Figure 15, is very well correlated (85%) with the e 2 qQ/h.   Results of QTAIM confirmed the presence of the hydrogen bonds and contacts detected by the total 2D FP. The energies estimated using QTAIM and Equations (7)-(12), Table 11, are in good agreement with those estimated using the energy framework. In contrast to the local 2D FP, QTAIM did not detect intramolecular N-H···N hydrogen bonds. However, the disc-shaped green surface of the reduced density gradient located between the -N(1)= and H (from -NH2) indicated the presence of a weak intramolecular N-H···N hydrogen bond, although no critical points were detected ( Figure 14). In the plot of RDG(r) versus sign(λ2)ρ(r), Figure 15, the noncovalent interactions-hydrogen bonds (strong and attractive), van der Waals contacts (weak), and steric interactions (ring closures, repulsive)-are visible as spikes. Their positions were identified based on the analysis of the closely related pyrazine structures. The intramolecular N-H···N hydrogen bond spike is shown in light blue. The analysis of the bonds in which each of the nitrogen atoms participates within QTAIM and RDG reveals significant disproportions between them. The -NH2 site participates in two Figure 15. A 2D plot of RDS(r) vs. sign(λ 2 )ρ(r) revealing spikes describing the weak non-covalent interactions of hydrogen bonds, van der Waals contacts, and steric interactions in the FVP crystal.

Binding Mode of Biologically Active Form
For FVP to be biologically active, it must undergo intracellular phosphorylation to its active form (FPV triphosphate). A false nucleoside, which is built by the viral RdRP into the nascent viral RNA, results in a "defective", mutated RNA [17,18]. Three structures of COVID-19 RNA-dependent RNA polymerase bound to FVP (7DFG [53]), Nsp7-Nsp8-Nsp12 SARS-CoV-2 RNA-dependent RNA polymerase in complex with template primer dsRNA, and FVP-RTP (7AAP [54]), and replicating polymerase complex of SARS-CoV-2 in the precatalytic state bound to FVP (7CTT [55]), were retrieved from the PDB database.
In 7DFG, the FVP moiety has adopted the same amide group conformation as in the solid state, while in 7AAP and 7CTT, the amide group of FVP-RTP is rotated by 180 • around the -C(2)-CONH 2 bond, although it is energetically unfavoured (by 5.8 kJ/mol). Consequently, the two types of structures differ significantly in the binding mode within the FVP moiety. In the 7DFG structure, FVP mimics guanosine, and through its amide group, forms a noncanonical base pair (the so-called wobble base pair) with uracil in the template RNA strand bonded via two intermolecular N-H· · · O hydrogen bonds of 2.452 Å and 2.268 Å. Two O-H· · · O hydrogen bonds of 2.909 Å (10.46 kJ/mol) and 2.939 Å (10.02 kJ/mol) link =O from FVP with SER682. The intramolecular hydrogen bond N-H· · · N of 2.711 Å closes a 6-membered ring, while the π· · · π stacking between nucleoside pairs (uracil· · · FVP and adenosine· · · uracil), where an amide group links N(1) with adenosine and C(5) with uracil, additionally stabilize the structure. In the 7AAP structure, FVP mimics the guanosine base, and through its amide group, forms a noncanonical base pair with cytosine in the template RNA strand. Three intermolecular hydrogen bonds, N-H· · · O of 2.569 Å (=O from FVP), N-H· · · O of 3.170 Å, and N-H· · · N of 2.740 Å (-NH 2 from FVP) bind the cytosine, and two other bonds, strong N-H· · · O (-NH from FVP) of 2.989 Å and weak C-H· · · O (=O from FVP) of 3.539 Å, bind SER682. This system is additionally stabilised by an intramolecular hydrogen bond, N-H· · · O, of 2.549 Å (linking =O with -NH 2 ) closing a 6-member ring. Moreover, the van der Waals-type side chain contact C· · · O of 3.24 Å and the hydrogen bond N-H· · · O of 3.28 Å (involving =O with C(6) and N(6) from adenosine) supports the π· · · π stacking interaction realised by C(5) of FVP and uracil (average stack distance 4.06 Å) and the N· · · N contact (linking -N(1)= of FVP and N(3) of uracil) of 3.22 Å.
In the 7CTT structure, two N-H· · · O hydrogen bonds of 2.4 Å and 3.2 Å and N-H· · · N of 2.74 Å, involving amide group through =O and -NH 2 , respectively, link the FVP moiety with cytosine in the template RNA. Two O-H· · · O hydrogen bonds of 2.44 Å and 3.38 Å link =O from FVP and one N-H· · · O of 2.99 Å links -NH 2 from FVP with SER682. This system is additionally stabilised by an intramolecular N-H· · · O hydrogen bond of 2.612 Å, closing a 6-member ring. The nitrogen atom from the -NH 2 group forms the π· · · π stacking with uracil (average stack distance 3.97 Å), while C(5) forms π· · · π with adenosine (average stack distance 4.32 Å). In addition to base pairing with cytosine or uracil, FVP is coordinated by Lys545, which is positioned to accept hydrogen bonds from the nitrogen -N1= in the pyrazine ring or is less likely to donate to the fluorine atom (3.4 Å and 3.7 Å or 3.6 Å, and 4.6 Å for cytosine or uracil base pairing, respectively). While in the solid F···H, contacts are preferred, and weak C-H···F bonds are present in FVP-RTP bound with the RNA strand, which fluorine cannot form because the nearest protons are too far. Fluorine also does not appear to be directly involved in the base pairing but still plays a role in the binding as it makes contact with neighbouring N, C, or O atoms in the stack. In 7DEF, it participates in the F· · · O contact of 3.321 Å and supports the C-H· · · O hydrogen bond of 3.065 Å, both linking the FVP moiety with the ribosyl ring and forcing specific conformation of RMP. Additionally, F· · · C (C from uracil) of 3.544 Å supports the π· · · π stacking. In 7AAP, fluorine participates only in the F· · · O contact of 3.674 Å, but supports the C-H· · · O hydrogen bond of 3.082 Å forcing the specific RTP conformation. In 7CCT, it participates in two F· · · O contacts of 4.957 Å and 4.691 Å, supports the C-H· · · O hydrogen bond of 2.694 Å to keep the specific RTP conformation, and, additionally, supports the π· · · π stacking by F· · · N (N from adenosine) of 4.151 Å. All of this fits the typical pattern of fluorine's role within drug structures, i.e., the impact on conformation.
The Considering the cofactor binding strength, the ligands can be ordered as FVP-RTP (active state) > FVP-RTP (precatalytic state)>FVP-RMP. The total binding affinity (between target protein and ligand) is also highest in 7AAP (−494.30 kJ/mol), followed by 7CTT (−471.58 kJ/mol), and 7DEF (−247.02 kJ/mol). However, the conformational change in RTP (precatalytic vs. active) may be dictated not only by its binding to the cofactor but may also result from an allosteric effect.
The allosteric sites in 7AAP were detected with Allosite-Pro [57] and PASSer [58] using a scheme [59,60]. Near the active site, there are two allosteric pockets on either side of the RNA strand. One of the pockets is composed of VAL557, ARG555, LYS545, THR680, SER682, THR556, ASP623, C10, THR687, SER759, ASN691, ALA688, ASP760, ARG553, A11, and U20, and, thus, includes critical catalytic and hydrophilic cluster residues. The procedure used to dock the FVP-RTP ligand in the protein pocket was identical to that previously described [61,62]. The search space was defined as a subset region of 9.0-15.0 Å around the allosteric site. After docking, the best pose that leads to the stabilization of the complex with the highest binding/docking score was selected.
Docking of the FVP-RTP directly into this pocket leads to its very strong bindings to hydrophilic cluster residues ARG555 (−86.32 kJ/mol), ARG553 (−70.46 kJ/mol), and LYS545 (−41.71 kJ/mol), conserved residues ASP623 (−69.16 kJ/mol) and ASP618 (20.84 kJ/mol), critical catalytic residues ASP761 (21.38 kJ/mol; repulsion), and two Mg cofactors (−216.38 and −215.77 kJ/mol). The total binding affinity is −402.96 kJ/mol. The conformation of FVP-RTP is very similar to that of the RNA template (root-mean-square deviation of atomic positions, RMSD, is only 4.98%). The contribution of the cofactor-ligand electrostatic interaction to the total protein-ligand binding is almost twice as high as the sum of the steric and hydrogen bond contributions. Such a strong direct binding of FVP-RTP to the active site and cofactor 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 binding modes calculated for the nearest vicinity of the single FVP molecule in crystal and FVP-RMP or FVP-RTP molecules in proteins are shown in Figure 16. For clarity, only single-molecule structures are shown.
Molecules 2023, 28, x FOR PEER REVIEW 30 of 39 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 binding modes calculated for the nearest vicinity of the single FVP molecule in crystal and FVP-RMP or FVP-RTP molecules in proteins are shown in Figure 16. For clarity, only single-molecule structures are shown. There is one more nonobvious difference between the solid and RMP/RTP, precatalytic/active states. The intramolecular hydrogen bond O-HO observed in the solid FVP does not occur in any substituted form. However, the weak intramolecular hydrogen bond N-HN is still present in the 7DEF and is replaced by a much stronger N-HO of 2.55 Å and 2.61 Å in active 7AAP and precatalytic 7CTT, respectively. Thanks to the N-HO intramolecular bond, FVP-RTP successfully mimics guanosine. However, this bond strongly affects the intermolecular binding capacity of -NH2 and =O. In any case, the intramolecular bonds and, consequently, the planar conformation of the FVP moiety, is not only solid-specific but is also retained in solution during the interaction of the molecule with RNA protein. Furthermore, unlike the solid state with the enol-like conformation of the amide, each active/bonded form of FVP is deprived of a proton at -OH, which implies There is one more nonobvious difference between the solid and RMP/RTP, precatalytic/active states. The intramolecular hydrogen bond O-H· · · O observed in the solid FVP does not occur in any substituted form. However, the weak intramolecular hydrogen bond N-H· · · N is still present in the 7DEF and is replaced by a much stronger N-H· · · O of 2.55 Å and 2.61 Å in active 7AAP and precatalytic 7CTT, respectively. Thanks to the N-H· · · O intramolecular bond, FVP-RTP successfully mimics guanosine. However, this bond strongly affects the intermolecular binding capacity of -NH 2 and =O. In any case, the intramolecular bonds and, consequently, the planar conformation of the FVP moiety, is not only solid-specific but is also retained in solution during the interaction of the molecule with RNA protein. Furthermore, unlike the solid state with the enol-like conformation of the amide, each active/bonded form of FVP is deprived of a proton at -OH, which implies a keto-like conformation of the amide. Thus, the acceptor-donor balance is shifted in favour of the acceptors. The hydrogen bond acceptors (two O, N, and not fulfilling this role, F), which exceed the hydrogen bond donors (-NH 2 ) in number, result in a more polar (thus easily soluble in water) compound. This, in turn, increases secondary electrostatic interactions characterized by no transfer or the sharing of electrons. The latter is an important factor for water-mediated hydrogen bonds that affect protein-RNA recognition. The binding modes described above support the assumption that FVP stops replication through noncovalent interactions in the active site and also suggest that the binding strength, and thus efficiency, may strongly depend on the conformation of the amide and triphosphate groups as well as the acceptor-donor balance. This may explain why the results of clinical trials in the treatment of COVID-19 with FVP have been so conflicting so far [25,63]. Changing the substituents to maintain the preferred conformation of the solid amide and triphosphate groups but also the number of acceptors seem to be good steps in the direction of better efficiency of this drug.
The intrinsic bond strength index [65], a new feature describing the degree of electron density sharing between the atom and its surrounding, DOI, calculated for single FVP and FVP-RTP molecules, partly supports the above observations. The highest DOI values were found for the acceptors -F (37.62%), followed by -C(6) (41.38%), -N(1)= (5.44%), and -N(4)= (1.78%). The lowest DOI were obtained for -NH 2 (0.38%), -OH (0.32%) acting as a donor, and =O (0.13%). The values of DOI for the conformation flipped by 180 • (nonsolid one) are only slightly different. They are the highest for F (37.91%), followed by -C(6) (42.12%), -N(1)= (5.47%), and -N(4)= 1.83%. The donor sites -NH 2 (0.29%) and -OH (from amide; 0.32%) as well as acceptor =O (0.35%) have very low DOIs. The differences caused by the flip of the amide group are, therefore, negligibly small. In contrast, the ribosylation and phosphorylation significantly change the DOI atomic values. The highest values of DOI for FVP-RMP are shifted from the FVP moiety to ribosyl, for C(1) (29.7%), followed by C(2) (12.7 %), C(3) (12.28%), and O(15) (9.5%), which is independent of the RTP conformation. This reflects the loss of the F's ability to share electron density with its surroundings in favour of RTP and the tendency of RTP to form multiple bonds evident in both structures. In RMP, the highest DOI are obtained for phosphoryl moiety atoms (32.4% for O(3), 34.78% for P, 4.9 for O(2), 4.33 for O(3), and 4.34% for O(4)), which reflects its ability to bind phosphate groups. Very low DOIs in the -NH 2 group indicate its strong donor properties, irrespective of the structure of FVP, FVP-RTP, or FVP-RMP. Thus, the -NH 2 group, which links neighbouring molecules in crystals and base pairs in the RNA template, mainly by strong intermolecular N-H· · · O bonds, plays a key role in both the crystal and protein structure. Its involvement in strong intramolecular N-H· · · O hydrogen bonds seems to be an advantage when it comes to planarity and mimicking guanosine, but it may be an obstacle that weakens the bond pairing. Moreover, these results are in very good agreement with the docking results, which indicate that the phosphoryl oxygens bind to the cofactor.

Material
The identity and purity of the active substance sample manufactured in China by Jiangsu Zenji Pharmaceuticals Ltd. have been checked in an independent laboratory with NMR and LC-MS methods, and it does not raise any doubts. The purity of the powdered sample was confirmed to be above 98%. The experimental study was performed using the FPV sample without prior recrystallization or any additional purification. The powdered sample was degassed and sealed in a glass ampoule.

Experimental-1 H-14 N NQDR
The 14 N NQR frequencies in FVP lie in the low-frequency range below 4 MHz, and direct NQR measurement on a small sample (approx. 0.5 g) is not possible due to a very low signal-to-noise ratio. Therefore, the multistep procedure was used. The 14 N NQR frequencies were determined indirectly by measuring the proton NMR signal using three nuclear quadrupole double resonance techniques: 1 H-14 N cross-relaxation (CR) [39], multiple frequency sweeps, and two-frequency irradiation [42,43].

Spectra Simulations-Density Functional Theory (DFT)
Quantum chemical calculations were carried out within the density functional theory (DFT) approach rooted in Kohn-Sham's [66] theorem, generalised by Levy [67]. The Minnesota hybrid meta M062X, high-nonlocality exchange-correlation functional with double the amount of nonlocal exchange (2X), [68] and all-electron split-valence basis set 6 -311 + G(d,p) were applied for the single molecule and cluster calculations. The advantage of M062X with 54% HF exchange is that it gives good results for the system with noncovalent forces, which is one of the biggest deficiencies of standard DFT methods. Because the X-ray crystallography at the usual resolution often fails to directly access the positions of light atoms, the proton positions are poorly identified. Therefore, the corrected (optimized) proton positions are routinely used by us in all analyses based on the electron density distribution. The calculations were performed using Gaussian 16 rev. C01 [69].
For  [76,77] with numerical radial functions basis-double-zeta without and plane-wave basis sets with an energy cutoff of 500 eV, were probed. The GGA (generalised gradient approximation) functionals (often called nonlocal) depend on both dρ/dr and ρ, which provides an increase in accuracy but with an additional cost. The plane-wave basis guarantees a monotonic convergence to the target wavefunction and does not exhibit a basis-set superposition error. A reciprocal space-based technique, useful in the evaluation of long-range interactions, was applied. The sampling of the Brillouin zone was carried out with the scheme of Monkhorst and Pack [78] (1 × 1 × 1 or 3 × 3 × 3 k-point separations) code. Because the components of the electric field gradient (EFG) tensor are highly sensitive to the structure quality, including well-determined positions of hydrogen atoms [28,29], the FPV structure was optimized upon the assumption of a fixed unit cell. All solid-state quantum chemical calculations were carried out within the CASTEP [79] code.
In both above-described approaches, the components of the electric field gradient (EFG) tensor were calculated. The principal components of a second-rank symmetrical electric field gradient (EFG) tensor, q ii = ∂ 2 V(r) ∂x 2 i (i = x, y and z; V(r)-external electrostatic potential, satisfying the relationship |q xx | ≤ q yy ≤ |q zz |), were obtained after the diagonalization of the following tensor, calculated at the selected level of theory: where (r) is the electron density, r i is the projection of the r vector onto x,y, and z-axes, and δ ij is Dirac's delta function. Although the EFG tensor consists of six independent components in the principal axis system, it is fully described by the three principal components and the three eigenvectors describing the orientation of the principal axes with respect to an arbitrary frame. Because it is traceless, only two of the eigenvalues are mutually independent. Therefore, EFG can be described with only two parameters: the quadrupole coupling constant e 2 qQ/h and the asymmetry parameter η. Both are related to the 14 N NQR frequencies, which are observable, according to the following equations: The 14 N NQR resonance frequencies can be derived from these parameters using the set of equations: A nuclear quadrupole moment for 14 N equal to 2.044 fm 2 [80] was assumed.

Quantum Theory of Atoms in Molecules (QTAIM)
A theoretical analysis of the topology of intermolecular interactions was performed within Bader's quantum theory of atoms in molecules (QTAIM) [31] approach. The DFT wavefunction was calculated using the M062X functional and all-electron 6-311 + G(d,p) basis set. The analysis of the stationary points (maxima, saddle points, or minima in the electron density) permits differentiation of the nucleus-, bond-, ring-, and cage critical points, denoted as NAP (Nuclear Critical Point), BCP (Bond Critical Point), RCP (Ring Critical Point), and CCP (Cage Critical Point). The so-called topological descriptors calculated at the BCPs include the electron density at BCP, ρ(r BCP ), three eigenvalues of the principal components of a Hessian matrix composed of second partial derivatives of the electron density describing the curvature of the electron density according to the principal axes λ 1 , λ 2 and λ 3 , Laplacian (the sum of Hessian eigenvalues), ∆ρ, electron density at the bond critical point, BCP, (ρ BCP (r)), its Laplacian ∆ρ BCP (r), the potential electron energy density, V BCP (r), the kinetic electron energy density, G BCP (r), and the total electron energy density, H BCP (r). Some of these descriptors inform about the nature of the interaction (λ 2 and ∆ρ), while the others (λ 3 , V BCP (r), G BCP (r), and H BCP (r)) allow a bond strength estimation. The analysis of the interactions according to QTAIM must be performed attentively because the results can be easily overinterpreted, especially in terms of the number or the nature of bonds. Therefore, the set of quantum indicators characterizing the bonds or sets of rules must be carefully analysed, especially in the case of intramolecular hydrogen bonds [81]. Despite the valid criticisms, experimental binding strengths are often in line with the model's predictions. Thus, the hydrogen bond energy can be estimated using the wellknown formulas derived from using the partitioning of the electron density scheme [82] E HB ≈ E E = 0.5V(r BCP ) by Espinosa−Molins−Lecomte [83], by Matta [84], E HB ≈ E EM = −223.08ρ(r BCP ) + 0.7423 (9) by Emamian (for neutral complexes) [85], by Gili (for O-H···O) [52], and E N = −3.09 + 239ρ(r BCP ) f or OH···O −2.03 + 225ρ(r BCP ) f or NH···O (12) by Nikolaienko [87]. In many cases, the NCI (noncovalent interaction) based on the reduced electron density gradient (RDG) helps to detect interactions which were uncertain or not revealed by QTAIM [32].
The reduced electron density gradient is defined as follows: RDG(r) = |∇ρ(r)| 2(3π) 1/3 (r) 4/3 (13) 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 noncovalent contacts are present in the crystalline 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 by only a slightly negative λ 2 indicates a weakly stabilizing interaction; and a spike associated with λ 2 > 0 indicates the absence of noncovalent interactions. In many cases, spikes do not reach s(r) = 0, indicating that the interactions were missed by QTAIM due to a failure to detect critical points.

Hirshfeld Surfaces (3D HS)
The exploration of intermolecular interaction patterns and packing capacities in the solid was performed within the 3D Hirshfeld surfaces (3D HS) approach [34]. Threedimensional HS is defined as the outer contour of the space that a molecule occupies in a crystalline environment. It is constructed with the use of the molecular weight function (a quotient of the promolecule and procrystal electron density). The descriptors such as d norm , shape index and curvedness of the surface mapped over 3D HS were evaluated [33]. The intermolecular interactions were visualized in the d norm surface mapped over 3D HS in a red-white-blue scheme, where red was used for short contacts such as hydrogen bonds, white for contacts of about van der Waals radii, and blue for the remaining, much longer, contacts. The shape index and curvedness of the surface mapped over 3D HS described its flatness. The decomposition of the 3D HS into a 2D 'molecular fingerprint' (2D FP) map (plot of the distances of each surface point to the nearby interior and exterior atoms d i versus d e ) [34] which summarizes the distribution of interactions of the molecule with its environment, was made. The surface contact data derived from the 2D FP were used to obtain the enrichment ratio, E XY , a descriptor defined as the ratio of the proportion of the actual contacts in a crystal and the theoretical proportion of random contacts. It describes the propensity to form or avoid contacts [88].
A further global quantitative characterization of the noncovalent interactions was delivered by the electrostatic potential (ESP) surface V S (r) (the molecular surface defined by electron density ρ(r) = 0.001 electron bohr 3 ) mapped onto a 3D HS. Using this technique, a complete description of interactions based on the electrostatic complementarity within the so-called Politzer [35,89,90] GIPF (generalised interactions properties function) approach was possible.

Molecular Docking
Molecular docking is a widely used technique in the screening of novel therapeutic agents, which requires not only the knowledge of the ligand structure but also the reliable 3D X-ray crystal structure of the protein. The conversion of the files with receptor and ligand structures to the .pdbqt format was made with MGLTools. The allosteric sites were identified using Allosite-Pro [57] and PASSer [58]. Molecular docking was performed using AutoDock [91] and AutoDock Vina98Vina [92]. Template docking and docking with defined searched space around the active site were probed. The best pose was selected based on the scoring function estimating the protein-ligand binding energy. The score, being a linear combination of steric, van der Waals, hydrogen bonding, electrostatic terms, torsion, and sp2-sp2 terms, was used. In the next step, the residues (side chains) closest to the active ligand were minimised with respect to the pose found, and the ligand was energy-minimised using standard potentials.

Conclusions
In the crystal, the -NH 2 group participates in two intermolecular hydrogen bonds of different strengths N-H···N and N-H···O, -N(4)= participates in a N-H···N intermolecular hydrogen bond, but -N(1)= only in a very weak putative intramolecular hydrogen bond N-H···N (revealed only in NQR spectra, local 2D FP, and RDG). The -NH 2 group plays a dominant role in the intermolecular bonds, irrespective of the environment (crystal or protein). In the crystal, it participates in intermolecular hydrogen bonds N-H···N and N-H···O, in the precatalytic state only in N-H···O, while in the active state in N-H···N and N-H···O, so the same binding motives are kept. Its involvement in intramolecular weak N-H· · · N and strong N-H· · · O hydrogen bonds seems to be an advantage when it comes to planarity and stacking (needed for the mimicking of guanosine), but it is an obstacle that weakens the bond pairing. The -N(1)= nitrogen plays only a supportive role, and -N(4)= is excluded from binding in the active form due to its substitution. This is reflected in the asymmetry parameters-the highest for asymmetrically bonded -NH 2 and the lowest for -N(4)= participating only in one bond coaxial with the Z-axis of the EFG tensor. The ordering of the most important interactions according to the decreasing energy in the crystals is as follows: O-H· · · O > N-H· · · O ≈ π· · · π stacking (parallelly oriented molecules) > N-H· · · N ≈ ≈ π· · · π stacking (flipped molecules) > C-H· · · O >> F· · · F > F-H· · · C, among which N-H· · · O, N-H· · · N, π· · · π stacking, and possibly F-H· · · C, are present also upon binding with the template RNA. The fluorine in this structure plays a minor role and affects only the conformation. Considering the cofactor binding strength, three ligands FVP-RMP and two FVP-RTP can be ordered as FVP-RTP (active state) > FVP-RTP (precatalytic state) > FVP-RMP. The total binding affinity (between target protein and ligand) is the highest for FVP-RTP (−494.30 kJ/mol; active state), followed by FVP-RTP (−471.58 kJ/mol; precatalytic state) and FVP-RMP (−247.02 kJ/mol). The total binding affinity for FVP-RTP in allosteric pocket is as high as −402.96 kJ/mol. Strong direct binding of FVP-RTP to both the active site and cofactor suggests a possible alternative, allosteric mechanism of FVP action, which may explain the scattering of the results of clinical trials or the synergistic effect observed in combined treatment against SARS-CoV-2.