Exploring Partial Structural Disorder in Anhydrous Paraxanthine through Combined Experiment, Solid-State Computational Modelling, and Molecular Docking

Abstract

Paraxanthine (PX), a major metabolite of caffeine, a protective agent against Alzheimer’s and Parkinson’s disease, and a promising drug for the treatment of post-COVID 2019 anosmia and ageusia, has been studied in the solid state and protein–ligand complex. Partial disorder in PX, caused by the methyl group at the N(7) position, has been modelled and discussed. The relationship between the unusual structural disorder and the propensity to form a specific system of non-covalent bonds was analyzed. Three ¹H-¹⁴N NMR-NQR (nuclear magnetic resonance–nuclear quadrupole resonance) experimental techniques were used, namely multiple frequency sweeps, Larmor frequency scanning, and the two-frequency irradiation, followed by solid-state computational modelling (density functional theory, supplemented by quantum theory of atoms in molecules, 3D Hirshfeld surfaces, and reduced density gradient), and molecular docking approaches. New quantitative methods for estimating changes in the global pattern of interactions under the influence of rotation of the methyl group in N(7) based on the Pompeiu–Hausdorff and Bhattacharayya metrics and the Wasserstein distance have been proposed and applied. A spectrum consisting of 12 lines, indicating the presence of 4 chemically inequivalent nitrogen sites in the PX molecule, was recorded, and the lines’ assignment to particular sites was made. The influence of the methyl rotation on the eigenvalues and eigenvectors of the electric field gradient tensor, NQR parameters, and resonance line positions was modelled in the solid (GGA/RPBE, m-GGA/RSCAN) and cluster (Minnesota M062X hybrid). Three factors have been found to determine structural disorder in PX: larger crystal voids near the methyl at N(7) than at N(1) (opening the path for the disorder), hyperconjugation strongly affecting the density distribution in the five-membered ring, and the involvement of the methyl group at N(7) in many non-covalent bonds that intercept (capture) subsequent jumping protons. The Pompeiu–Hausdorff and Bhattacharayya metrics and the Wasserstein distance confirmed the changes in the distribution and strength of non-covalent interactions throughout the molecule as a result of methyl rotation. This effect is clearly visible regardless of the type of metric, and its order of magnitude is consistent with the modulation effect of the NQR spectra (experimental and calculated). Through molecular docking, it was discovered that the PX moiety in protein–ligand complexes adopt the same methyl group conformation at N(7) as in the solid state. It was found that the cooperation–competition between the C-H⋯O hydrogen bonds and C-H⋯H-C dispersion interactions is the crucial factor that impedes methyl rotation and induces structural disorder, as well as being an important factor in the formation of the protein–ligand complexes.

1. Introduction

Methylxanthines are natural purine alkaloids found in beverages, dietary components, supplements, and pharmacologic agents (drugs and para-pharmaceuticals) that have important clinical applications [1,2]. One of them, paraxanthine (PX; 1,7-dimethylxanthine; 1,7-dimethyl-3H-purine-2,6-dione), Figure 1, is the major metabolite of caffeine (CF; 1,3,7-trimethylxanthine), namely 1,3,7-trimethyl-1H-purine-2,6-(3H,7H)-dione) [3,4].

It is solely the result of demethylation of caffeine at position 3 by the catalytic action of cytochrome P450 subtypes 1A2 (90%) and 1A1, 2E1, 3A4, and 2D6 (10%) in liver parenchymal cells (hepatocytes) [5,6,7,8]. PX exists as a xenobiotic metabolite in all living organisms, from bacteria [9] to mammals [10], but no known plant or food contains PX because it is not accumulated (due to very slow N1-methylation of 7-methylxanthine) [11].

It is believed that the contribution of PX to the pharmacological action attributed to CF is significant. Molecular mechanism of PX action is determined by its potent but nonselective antagonism to adenosine A1, A2A, and A2B receptors, which are involved in such processes as regulation of myocardial, oxygen consumption, and coronary blood flow, as well as in neurotransmitters, such as dopamine and glutamate, and in inflammation and immune responses [12]. Thus, it affects fundamental human processes, such as sleep, arousal, cognition, learning, and memory, and can stimulate the central nervous system, as well as the muscular, respiratory and circularity systems [4,13,14,15,16,17,18,19]. PX treatment may affect short-term memory, reasoning, response time to cognitive challenges, and help to sustain attention [19]. PX, unlike CF, is a potent dopaminergic and protects against neurodegeneration and loss of synaptic function [20]. PX, rather than CF, has recently been reported to reduce the risk of developing Alzheimer’s and Parkinson’s diseases [4,15,21,22]. It was approved by the FDA for the treatment of most reversible obstructive airway diseases (e.g., chronic bronchitis, emphysema, and asthma), but later relegated from the treatment of chronic obstructive pulmonary disease (COPD) due to an imbalance between benefits and side effects [16]. Unapproved by the FDA indications include sleep apnea and infant apnea cardiogenic pulmonary [23], obesity reduction, and treatment of chronic fibrogenic or even carcinogenic liver disorders [24]. PX also attracts some interest in agriculture due to its significant allelopathic potential [25,26]. Recently, during the coronavirus pandemic, PX, which may act as an antagonist of adenosine receptors that modulate taste (A2B) and odor (A1, A2A), has gained interest as a promising drug in the treatment of anosmia (smell blindness) and ageusia (taste blindness), the most common post-COVID 2019 complications [27,28]. Given that anosmia is associated with an increased risk of depression, PX also seems to offer hope for its treatment [4,29].

Although PX has a very broad spectrum of activity and application area, only limited physicochemical data are available on this compound. It is much less studied than other methylxanthines [30] and has not even been classified by the International Agency for Research on Cancer (IARC) [5]. Limiting factors are the lack of natural sources of PX, its fast metabolism, and the high cost of samples. Recently, we explained the role of hydrogen bonding and π⋯π stacking in the context of the disorder and solubility of CF polymorphs [31] and its two minor metabolites, TB and TP [32]. Although PX, unlike CF, TB, and TP, is poorly soluble in water (only 1.0 for PX versus 2160, 330, and 7.36 mg/mL at 25 °C for CF, TB, and TP, respectively [33]), the anhydrous PX crystal was obtained and its structure has been resolved by X-ray diffraction with a final R-factor of 3.81% at 150 K and RT [34]. In contrast to CF, PX easily crystallizes in the anhydrous form, which was attributed exclusively to the formation of lateral hydrogen bonds [34]. CF exhibits typical disorder that arises in the cases when there are different, almost thermodynamically equivalent, variants of the molecule arrangement. Subtle changes in crystallization conditions or phase transitions cause CF molecules to rearrange from symmetrically disordered to totally ordered. However, disordered regions still may be present in the crystal, but do not diffract X-rays due to the effect of structural heterogeneity or dynamic motions. In the PX structure, only one methyl group (the one at N(7)) is disordered, and only at 293 K (at 100 K no disorder is detected); hence, the structural disorder is unique and refers to the multiple conformations of the substituents in the crystal lattice.

This paper attempts to explain the nature of this unusual structural disorder in relation to the propensity to form specific non-covalent bonds. Cooperation–competition of specific non-covalent long-range (electrostatic, inductive, dispersive, and attractive) and short-range (exchange and repulsion) interactions is crucial for the stabilization of solid phase structures. The interplay of various canonical interactions (hydrogen bonds, hydrophobic interactions) and a variety of non-canonical interactions (C–H⋯π, O–H⋯π, N–H⋯π, π⋯π stacking) play an important role in protein stability. So far, the intermolecular hydrogen bonds N–H⋯O have been considered in PX [34]. However, the presence of two kinds of strong acceptors (N, O) and the disorder associated with two methyl groups donating protons suggests that mainly electrostatic C–H⋯O and C–H⋯N hydrogen bonds (intra- or intermolecular) or homopolar C–H⋯H–C dihydrogen bonds governed by dispersion, may also occur and be of importance. Both intramolecular hydrogen bonds and dispersion interactions deliver the energetic contributions, the importance of which has been increasingly noticed in different compounds [35,36,37,38,39]. The key question is what types of non-covalent bonds are actually formed in a disordered crystal structure of PX, how strong they are, and to what extent their individual binding capabilities facilitate or hinder order/disorder in the solid state or formation of the protein–ligand complex.

We applied a unique technique, ¹H–¹⁴N NMR-NQR double resonance spectroscopy, which is highly sensitive to the changes in the local environment and more sensitive than NMR, IR, or powder X-ray diffraction (PXRD) [40]. Although NQR spectra, which are sensitive to all interactions modifying the electron density in the nearest vicinity of the nitrogen site, give some clues about the bonds involving nitrogen or their nature, they cannot provide a complete description of intermolecular interactions in the crystal lattice. However, in combination with the solid-state computational methods (Bader’s quantum theory of atoms in molecules [41], QTAIM; Johnson’s reduced density gradient, RDS [42], and Spackman’s Hirshfeld surfaces [43,44]), it has been discovered to be a valuable tool for this purpose [31,32,38,39]. A combination of the above techniques allows determination of a good quality structure, and, thus, a high-quality wavefunction, which is crucial to discover the role of non-covalent bonds in the disorder of methyl groups. New quantitative methods for estimating changes in the global pattern of interactions under the influence of rotation of the methyl group in N(7) based on the Pompeiu–Hausdorff and Bhattacharayya metrics and the Wasserstein distance have been proposed and applied. This approach is beneficial for evaluating the sources of partial crystal disorder. The specific binding mode, critical to biological activity, can be strengthened or disrupted by the disorder, affecting certain physicochemical properties, such as polymorphism, solubility, and stability. Structurally disordered regions play an important role in protein–protein and protein–ligand interactions, facilitating their mutual matching. It has been shown that a comprehensive study of the interaction pattern in disordered structures can be a predictor of the capability of individual conformers to bind to the target.

In view of the above, we hope that our research will be beneficial in the areas related to drug development.

2. Materials and Methods

2.1. Material

All material was purchased from Sigma-Aldrich, Poznań, Poland, PL, unless otherwise stated. The quality control of raw material was performed by FT IR. The purity of the powdered sample was confirmed to be above 98%, higher than the ACS grade. The PX sample was not subjected to prior recrystallization or any additional purification. The powdered sample was degassed and sealed in a glass ampoule.

2.2. Methods

Experimental Nuclear Quadrupole Resonance (NQR)

Measurements were performed at room temperature (295 K). We used two different techniques, namely multiple frequency sweeps and Larmor frequency scanning, to determine the lowest ¹⁴N NQR frequencies ν₀ at all nitrogen positions in the molecule and to locate the ¹⁴N NQR frequencies ν₊ and ν₋. The two-frequency irradiation technique was finally used to increase the resolution of determination of the NQR frequencies ν₊ and ν₋ to about ±5 kHz. A complete description of the basics of these techniques can be found in our previous papers.

2.3. Computational Density Functional Theory

2.3.1. Spectra Simulation—Cluster Technique

The quantum chemical calculations required for the QTAIM analysis were carried out within the density functional theory (DFT) approach rooted in the Kohn–Sham [45] theorem, generalized by Levy [46]. This approach is an attractive compromise, providing useful accuracy and favorable algorithmic complexity. The Minnesota M062X hybrid meta, a high-nonlocality exchange–correlation functional with a double amount of nonlocal exchange (2X) [47], and an all electron split-valence basis set 6-311 + G(d,p), were used for the single molecule and cluster calculations. M062X, which has a 54% HF exchange, gives high values for the systems with non-covalent forces, which is one of the fundamental disadvantages of classical DFT approaches. X-ray crystallography at the usual resolution often fails to directly access the positions of light atoms; thus, the proton positions were corrected (optimized). A complete description of the basics of the cluster technique can be found in our previous papers [31,32,48]. The calculations were performed using a Gaussian 16 rev. C01 [49].

2.3.2. Spectra Simulation—Solid-State Technique

Different generalized gradient approximations (GGA) functionals, which depend on both ρ and dρ/dr, were used. They provide an increase in accuracy, but often with an additional computational cost. The GGA functionals, namely RPBE (revised Perdew, Burke, and Ernzerhof) [50], PBE (Perdew, Burke, and Ernzerhof) [51], WC (Wu and Cohen) [52], PBEsol (Perdew et al.) [53], PW91 (Perdew et al.) [54], were probed. The addition of the Tkatchenko–Scheffler [55] (DFT-TS) correction for dispersion to some GGA functionals was found to be useful. Additionally, the usefulness of the modern semi-local density functional, which satisfies all known constraints that an exact density functional must meet, namely m-GGA rSCAN (Sun, Ruzsinsky, and Perdew [56,57], revised by Bartok and Yates) [58] was checked. The gauge including projector-augmented waves (GI)PAW exploiting the full translational symmetry of a crystal and on-the-fly generation ultrasoft (OTFG) potentials were used.

The principal components of a second-rank symmetric electric field gradient (EFG) tensor, $q_{ii}=\frac{{\partial }^{2}V(r)}{\partial x_i^2}$ (i = x, y and z; V(r)—external electrostatic potential, satisfying the relationship $\left|q_{xx}\right|\le \left|q_{yy}\right|\le \left|q_{zz}\right|$), were obtained after diagonalization of the following tensor, calculated at the selected level of theory:

$$\nabla E_{ij}=-\frac{1}{4\pi {\epsilon }_0}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left[\frac{3r_i{r}_j-{\delta }_{ij}{\left|r\right|}^2}{{\left|r\right|}^5}\right]\rho \left(r\right)dr$$

(1)

where ρ(r) is the electron density, r_i is the projection of the vector onto the x, y, and z axes, and δ_ij is Dirac’s delta function. The EFG tensor consists of six independent components in the principal axis system. It is fully described by the three principal components and three eigenvectors, which describe the orientation of the principal axes with respect to an arbitrary frame. The tensor is traceless; therefore, only two eigenvalues are mutually independent. For this reason, two parameters are sufficient to describe the EFG tensor, namely the quadrupole coupling constant e²qQ/h and the asymmetry parameter η. The ¹⁴N NQR frequencies, which are observable, are related to these parameters as follows:

$$\left|\frac{e^{2}Qq}{h}\right|=\frac{2}{3}\left({\nu }_{+}+{\nu }_{-}\right)$$

(2)

$$\eta =\frac{3\left({\nu }_{+}-{\nu }_{-}\right)}{{\nu }_{+}+{\nu }_{-}}$$

(3)

The ¹⁴N NQR resonance frequencies can be derived from these parameters using the following set of equations:

$${\nu }_{+}=\left|\frac{e^{2}Qq}{4h}\right|(3+\eta )$$

(4)

$${\nu }_{-}=\left|\frac{e^{2}Qq}{4h}\right|(3-\eta )$$

(5)

$${\nu }_0={\nu }_{+}-{\nu }_{-}=\left|\frac{e^{2}Qq}{2h}\right|\eta$$

(6)

which are useful in the case of theoretical calculations. The nuclear quadrupole moment of ¹⁴N equal to 2.044 fm² [59] was assumed.

The components of the electric field gradient (EFG) tensor are highly sensitive to the structure quality, including the positions of hydrogen atoms [48]. Therefore, the PX structure was partially and fully optimized. Full optimization was performed upon the assumption of a fixed unit cell.

2.4. Analysis of the Topology of Interactions—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) [41]. The DFT wavefunction calculated using the Minnesota M062X hybrid functional and an all electron 6-311 + G(d,p) basis set was used. The differentiation of the nucleus, bond-, ring-, and cage-critical points, denoted as NAP (nuclear critical point; local maximum of the electron density with all curvatures negative), BCP (bond-critical point; saddle point local maximum of the electron density with two curvatures negative), RCP ring-critical point; saddle point local minimum of the electron density with two curvatures positive), and CCP (cage-critical point; local minimum of the electron density with all curvatures positive) was based on analysis of the stationary points. The so-called topological descriptors, including the electron density at BCP, ρ(r_BCP), three eigenvalues of the principal components of Hessian matrix, namely λ₁, λ₂, and λ₃, 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), were calculated at each BCPs. The λ₂ and Δρ descriptors provide information about the nature of the interactions. Four other descriptors (λ₃, V_BCP(r), G_BCP(r), and H_BCP(r)) allow a bond strength estimation. In general, the hydrogen bond energy can be estimated using the well-known formulas derived using partitioning of the electron density scheme [60] proposed by Espinosa−Molins−Lecomte [61], as follows:

$$E_E=0.5V(r_{BCP})$$

(7)

and by Matta [62], as follows:

$$E_M=-0.429G(r_{BCP})$$

(8)

and by Emamian (for neutral complexes) [63], as follows:

$$E_{EM}=-223.08\mathsf{\rho }\left(r_{BCP}\right)+0.7423$$

(9)

and by Afonin [64], as follows:

$$E_A=-\left(0.277\left|\mathrm{V}\left(r_{BCP}\right)\right|+0.45\right)={-(0.554|E}_{HBE}|+0.45)$$

(10)

and by Nikolaienko [65], as follows:

$$E_N=-2.03+225\mathsf{\rho }\left(r_{BCP}\right)$$

(11)

However, these formulas must be treated attentively because the results can be over-or under-interpreted in the case of intramolecular hydrogen bonds.

The reduced electron density gradient (RDG) is derived from QTAIM parameters and defined as follows:

$$RDG(r)=\frac{\nabla \mathsf{\rho }|(r)|}{{2(3\pi )}^{1/3}{\mathsf{\rho }(r)}^{4/3}}$$

(12)

where ∇ρ(r) is the electron density gradient and ρ(r) is electron density; Equation (12) helps to detect interactions that were uncertain or not revealed by QTAIM [42]. The plot of RDG(r) versus the electron density multiplied by the sign of the λ₂ Hessian eigenvalue, sign(λ₂)ρ(r), reveals characteristic spikes in the low-gradient and low-density regions as long as the non-covalent contacts are present in the crystalline structure. The spikes not reaching RDG(r) indicate interactions missed by QTAIM due to failure to detect critical points. The nature of non-covalent interactions can be inferred from the sign of λ₂. Depending on its value, they are classified as attractive (stabilizing; λ₂ < 0) or repulsive (destabilizing; λ₂ > 0). A spike in the low-gradient, low-density region at λ₂ < 0 suggests a stabilizing interactions, such as a hydrogen bonding; a smaller spike accompanied with only slightly negative λ₂ indicates a weakly stabilizing interaction and a spike associated with λ₂ > 0 indicates the absence of non-covalent interactions. The results of the RDG analysis were visualized using VMD [66].

2.5. Exploration of Intermolecular Interaction Patterns—3D Hirshfeld Surface Analysis

Exploration of intermolecular interaction patterns and packing capacities in the solid was performed within the 3D Hirshfeld surfaces (3D HS) approach [44]. 3D HS is defined as the outer contour of a space, in which a molecule occupies in a crystalline environment. This surface was constructed using the so-called molecular weight function (a quotient of the promolecule and procrystal electron density). Three types of descriptors were evaluated, namely d_norm, shape index, and curvedness of the surface [43]. Intermolecular interactions were visualized on the d_norm surface mapped to the 3D HS in a red–white–blue scheme (red was used for short contacts, white for contacts of about van der Waals radii, and blue for the remaining ones). The flatness of the 3D HS was described using the shape index and curvedness of the surface mapped over the 3D HS being described. The decomposition of the 3D HS into a 2D ‘molecular fingerprint’ (2D FP) map [44] was carried out. The 2D FP (plot of the distances of each surface point to the nearby interior and exterior atoms d_i versus d_e) summarizes the distribution of interactions of the molecule with its environment. The enrichment ratio, E_XY, a descriptor which is defined as the ratio of the proportion of the actual contacts in a crystal and the theoretical proportion of random contacts [67] was calculated based on 2D FP. It describes each atom’s propensity to form or avoid contact.

2.6. Comparison of the Differences in Interactions Patterns Due to the Methyl Group Rotation—Mathematical Metrics

A comparison of the differences/similarity between of the intermolecular interactions pattern due to the methyl group rotation was performed using a few metrics defined as follows:

• Pompeiu–Hausdorff metrics [68], which measures how far two subsets of a metric space are from each other, as in the following Equation (13):

  $$d_{PH}(p,q)=\max \{{\mathrm{s}\mathrm{u}\mathrm{p}}_p(d\left(p_i,p\right)),{\mathrm{s}\mathrm{u}\mathrm{p}}_q(d\left(q_i,q\right))$$

  (13)

  where d can be calculated as the Manhattan distance, as follows:

  $$d_{ML}(p,q)=\sum _{i=1}^n|p_i-q_i|$$

  (14)

  or as the Euclidian distance, as follows:

  $$d_E(p,q)=\sum _{i=1}^n{(p_i-q_i)}^2$$

  (15)
• The Bhattacharayya coefficient [69], which is a distance between two probability distributions, as follows:

  $$d_{BC}(p,q)=\sum _{i=1}^n{p}_i{q}_i$$

  (16)
• The Wasserstein distance (Kantorovich–Rubinstein metric) [70], which is a function defined between probability distributions on a given metric space, as follows:

  $$d_W\left(p,q\right)={\mathrm{i}\mathrm{n}\mathrm{f}}_{dist}\sum _{i=1}^n|p_i-q_i|$$

  (17)

  where sup represents the supremum, inf represents the infimum, p_i, q_i are interactions in each structure, and dist are permutations of all p_i, q_i , and sum runs overall.

We previously successfully used the Manhattan and Euclidian metrics and Bhattacharayya distance to compare interaction patterns across polymorphic forms [71].

2.7. Exploration of Binding—Molecular Docking

Molecular docking (MD) is a widely used technique in screening of novel therapeutic agents. This technique requires both the knowledge of the ligand structure and a reliable 3D X-ray crystal structure of protein. The MGLTools software (ver 1.5.7) was used to convert the files with these structures to the .pdbqt format. MD was performed using AutoDock [72] and AutoDock Vina98Vina [73]. Two techniques were used to dock a template docking and docking with defined searched space around the active site. The score function being a linear combination of different terms (steric, van der Waals, hydrogen bonding, electrostatic, torsion and sp²-sp²) was used. Based on the scoring function’s estimation of the protein–ligand binding energy, the optimum pose was chosen. The approach for docking the ligands in the protein pocket was identical to that described earlier [38]. Exploration of protein–ligand hydrophobic interactions was performed using LigPlot+ [74,75].

3. Results and Discussion

3.1. ¹H-¹⁴N NQDR Spectrum

As the first technique for the location of ¹⁴N NQR frequencies, we used multiple frequency sweeps [76] and Larmor frequency scanning [77]. The proton polarization period was set to 30 s and the relaxation period was set to 0.5 s. We used linear frequency sweeps with the frequency limits of 1.6 MHz and 3.5 MHz. The duration of a single sweep was set to 10 ms. The average amplitude of the rf magnetic field during the sweep was approximately equal to 2 mT. The proton Larmor frequency scan was performed in steps of 10 kHz. At each value of the proton Larmor frequency, we took four measurements and made an average of the four proton NMR signals. The measurements were performed at room temperature (295 K).

Using this technique, we observed a broad dip in the frequency range between 220 kHz and 270 kHz and two narrower dips centered at 530 kHz and 730 kHz. The dips are shown in Figure 2 as the low-frequency part of the ¹⁴N NQR spectrum of PX. Since there are four nitrogen positions in a PX molecule, the broad dip represents an overlap of two dips.

In the next step, we located the higher ¹⁴N NQR frequencies ν₊ and ν₋ by varying the frequency limits of the sweeps. First, we set the proton Larmor frequency to 220 kHz and, by varying the upper frequency limit, we observed a step in the proton NMR signal at a frequency of 2.690 MHz. By varying the lower frequency limit, we observed a step at a frequency of 2.460 MHz, so the three ¹⁴N NQR frequencies corresponding to this nitrogen position, which we label as N(1), are approximately equal to 2.690 MHz, 2.460 MHz, and 230 kHz, respectively. Then, we fixed the proton Larmor frequency at 270 kHz and repeated the measurement. We observed the steps in the proton NMR signal at the frequencies 1.910 MHz and 1.650 MHz, so the three ¹⁴N NQR frequencies from this nitrogen position, which we label as N(2), are approximately equal to 1.910 MHz, 1.650 MHz, and 260 kHz. At the third nitrogen position, N(3), the ¹⁴N NQR frequencies were located at 2.870 MHz, 2.340 MHz, and 530 kHz, whereas at the fourth nitrogen position, N(4), they were located at 2.780 MHz, 2.050 MHz, and 730 kHz.

In order to increase the resolution to about ±5 kHz, we used the two-frequency irradiation technique as the final step [76]. The two rf magnetic fields were applied as a series of repetitive pulses with a pulse length of 1 ms. The amplitudes of the two rf magnetic fields were set to 0.2 mT. The experimentally determined ¹⁴N NQR spectrum is shown in Figure 2, and the frequencies e2qQ/h and η, calculated using Equation (3), are listed in

Table 1. ¹⁴N NQR frequencies, quadrupole coupling constants e²qQ/h, and the asymmetry parameters η in paraxanthine at T = 295 K.

Nitrogen Site	ν+ [MHz]	ν− [MHz]	ν0 [MHz]	e2qQ/h [MHz]	η	Final Assignment *
N(1)	2.685	2.457	0.228	−3.428	0.133	>N(1)CH3
N(4)	2.780	2.047	0.733	−3.218	0.456	>N(3)H
N(2)	1.910	1.650	0.260	−2.373	0.219	>N(7)CH3
N(3)	2.870	2.340	0.530	−3.473	0.305	-N(9)=

* The final assignment made on the basis of quantum chemical calculations.

.

Figure 2 is in fact a merge of five double resonance spectra, with the spectrum below 800 kHz obtained by multiple frequency sweeps and proton Larmor frequency scanning, and the other four spectra obtained by the two-frequency irradiation technique. The low-frequency dips are broader than the high-frequency ones, which is a result of the idiosyncrasy of the techniques used. The width of the low-frequency dip is mainly determined by the proton NMR linewidth, whereas the width of the high-frequency dip is determined by magnetic field broadening of the ¹⁴N NQR line and by the inhomogeneity of the low magnetic field.

The spectrum of PX, as in Figure 2, similarly to that of TP [32], is “clean”, i.e., the lines are not anomalously broadened, with a similar full width at half maximum (FWHM), which suggests no evidence of structural disorder. The lack of multiplicity suggests the presence of equivalent molecules per unit cell, which is in good agreement with the X-ray data [34]. The electronic environment of all four nitrogen nuclei in the PX molecule is different, so they are chemically inequivalent. All four nitrogen sites have distinctly different e2Qq/h and η values, although two of them, N(1) and N(7), are equally substituted with an electron donating and hyperconjugative methyl group, whose σ-bond is delocalized toward the adjacent π orbital, but N(1) and N(7) belong to different heterocyclic rings, six- and five-membered, respectively. The five-membered ring is more strained, which affects the values of the components of the EFG tensor and the orientation of its principal axes. In addition, the methyl group at N(1) is located between two oxygen atom species, whose lone pairs are adjacent to the π system. The oxygen atom exerts an electron-withdrawing (inductive) effect, but its lone pairs cause an electron-donating effect (resonance effect), which is the opposite and much stronger. Thus, N(1) should be electron-rich compared to N(7), and its η should be smaller. Other clues to help assign the resonance lines come directly from a comparison of NQR parameters with those previously obtained for TB or TP [31,32]. Based on these data, it seems that the two sets of resonance lines yielding the lowest η should be assigned to the sites N(1) and N(7), while the other two should be assigned to the remaining N(3) and N(9) ones, but it is by no means certain. Several overlapping factors, in addition to simple substituent effects, such as the poorly differentiated environment of quadrupolar sites (especially methyl ones), their high sensitivity to proton positions, and the disorder of protons in the methyl group at N(7), make it difficult to reliably assign the signals to nitrogen sites in PX. The quantitative analysis of NQR parameters in solid systems may be further complicated by influences from the surroundings, including various non-covalent bonds, which disturb the bonding situation of the nitrogen atoms.

As a result, an unambiguous assignment of nitrogen signals based solely on the NQR spectra is not possible. The simulation of the ¹⁴N NQR spectra of the other heterocyclic compounds has revealed that NQR parameters were extremely sensitive, even to the influence of distant nitrogen atoms [31,48]. Therefore, we performed quantum chemistry calculations in two variants (cluster and solid-state). They make it possible to take into account the effects of substituents, hybridization, and intermolecular interactions, including those that are distant, when modeling NQR spectra. Both strategies require a consistent, high-quality crystalline structure as the input. The Cambridge Structural Database (CSD) contains only two PX crystalline structures. While the low-temperature structure of PX (CCDC FIDSAN [34]) is fully ordered, in the RT structure (CCDC FIDSAN1 [34]) the methyl group substituted at N(7) equally occupied two orientations. These orientations, named ABC and DEF, differ by an angle of 60° (which corresponds to the torsion angle 120°), Figure 1. Rotational disorder, per se, is not surprising, since the disorder was observed in CF [31,34], and the classic Arrhenius-type rotation of methyl groups is expected near RT.

Therefore, various models have been proposed to study and describe the structural disorder that can affect the assignment of lines in the resonance spectrum. In the first model, a low-temperature ordered structure (with X-ray protons, with optimized protons, and a fully optimized structure with a fixed unit cell) was used as the input for the calculations. Since in the classical, Arrhenius type of methyl group rotation, the classical jumps of methyl groups at low temperatures are considered to be frozen, this model with a good approximation describes a completely ordered system and can be treated as the so-called reference point (corresponding to a rigid structure). Then, a simple model of disorder (the so-called direct disorder), assuming a disordered PX structure at RT, with proton occupation determined directly from X-ray data, was analyzed. This structure was probed in three variants (with X-ray protons, with protons optimized, and a fully optimized structure with a fixed unit cell). Finally, the averaging model, that assumes the average of the two substructures RT ABC and RT DEF retracted from disordered structure at RT and differing only in the conformation of the methyl group at N(7), was proposed. The calculations were performed for this model in three variants (with X-ray protons, with protons optimized, and a fully optimized structure with a fixed unit cell). The results are listed in

Table 2. The ¹⁴N NQR parameters for paraxanthine calculated theoretically at the GGA/RPBE level with TS DFT-D correction.

Structure	Site	ν+ [MHz]	ν− [MHz]	ν0 [MHz]	e2qQ/h [MHz]	η	s, r2, Intercept, Slope
LT 150 K protons optimized	>N(1)CH3	2.696	2.556	0.140	−3.501	0.08	0.052
	>N(3)H	2.799	2.027	0.772	−3.217	0.48	0.995
	>N(7)CH3	2.019	1.767	0.252	−2.524	0.20	0.016
	-N(9)=	2.785	2.388	0.397	−3.449	0.23	1.002
RT ABC protons optimized	>N(1)CH3	2.675	2.553	0.122	−3.485	0.07	0.073
	>N(3)H	2.809	2.021	0.789	−3.22	0.49	0.993
	>N(7)CH3	2.017	1.777	0.240	−2.529	0.19	−0.041
	-N(9)=	2.776	2.381	0.395	−3.438	0.23	1.027
RT DEF protons optimized	>N(1)CH3	2.696	2.556	0.140	−3.501	0.08	0.063
	>N(3)H	2.799	2.027	0.772	−3.217	0.48	0.994
	>N(7)CH3	2.019	1.767	0.252	−2.524	0.20	−0.038
	-N(9)=	2.785	2.388	0.397	−3.449	0.23	1.027
RT ABCDEF averaged protons optimized	>N(1)CH3	2.685	2.554	0.131	−3.493	0.08	−0.039
	>N(3)H	2.804	2.024	0.780	−3.219	0.49	1.027
	>N(7)CH3	2.018	1.772	0.246	−2.527	0.20	0.068
	-N(9)=	2.781	2.385	0.396	−3.444	0.23	0.994
direct disorder protons optimized	>N(1)CH3	2.884	2.734	0.150	−3.745	0.08	9.365
	>N(3)H	2.989	1.806	1.183	−3.197	0.74	0.391
	>N(7)CH3	4.709	2.866	1.843	−5.050	0.73	0.781
	-N(9)=	2.518	2.102	0.416	−3.080	0.27	0.822

, Tables S1 and S2, and visualized in Figure 3 and Figure S1.

The assumption of the low-temperature (150 K) PX structure and the crystallographic proton positions resulted in a significant overestimation of η for >N(3)H and unreliable set of the resonance lines, which suggests that its proton is involved in a strong hydrogen bond and shifted towards the acceptor (=O). This unsatisfactory result was found to be independent of the choice of the functional or basis set. Full optimization of the geometry, with a fixed unit cell allowed us to obtain NQR parameters that were quite well correlated with the experimental ones (high Pearson’s coefficient and small standard deviation) and underestimated by only about 4%. Unfortunately, such a procedure requires costly and time-consuming calculations, regardless of the level of theory. Moreover, the conformation of the methyl groups after optimization is far from the original one.

Calculations performed at the theoretical level of GGA/RBPE, assuming a low-temperature structure and proton positions optimized at the same level, allow reasonable NQR parameters to be obtained at a relatively low computational cost. The Pearson correlation coefficient is high, the dispersion is low and the NQR frequencies obtained from these parameters are reproduced very well. The modern m-GGA/RSCAN functional, which is much more computationally demanding, produces results similarly well correlated with experimental data, but the NQR parameters are reproduced with slightly lower accuracy. The use of the other GGA functionals (e.g., BLYP, WC, PBSOL, PW91, M11L, M06) gave the same signal assignments, although their use results in a set of frequencies less consistent with the experimental ones (the discrepancy concerns mainly nitrogen >N(3)H and, to a lesser extent, -N(9)=).

An attempt to model disorder defined directly i.e., by assigning it to particular protons based on the X-ray data (RT direct disorder), failed. The electron densities at both nitrogen atoms, N(7) and, to a much lesser extent, N(1), are poorly reproduced using this model and have led to the unreliable EFG tensor components and resonance frequencies, as in Table 2 and Figure S1c,d. It is, therefore, obvious that the wavefunction obtained by this method cannot be used for further analysis of the non-covalent interactions.

On the contrary, the model that averages two structures RT ABC and RT DEF (named as RT ABCDEF averaged) is quite reliable, as seen in Table 2 and Figure S1a,b. Both structures, RT ABC and RT DEF, differ in the conformation of the methyl group at N(7), as in Figure 1, in which only in the RT DEF structure is close to the low-temperature one. The underestimation of e²Qq/h and overestimation of η only for the >N(3)H site suggests its participation in the strong hydrogen bond N–H⋯O and, consequently, the shift in the proton from nitrogen towards the acceptor (=O). In the PX structure, the N–H⋯O bonds do not occur singly, but in pairs i.e., two N(3)H⋯O(2) hydrogen bond link two PX molecules into a dimer. At higher temperatures, the probability that the concerted proton jumps in such a double bond increase. These jumps average the electron density at the nitrogen site >N(3)H and the EFG tensor, which results in the lower-than-expected values of e²Qq/h and η.

The cluster calculations confirmed the assignment of the signals obtained based on the solid-state calculations, but also provided a high-quality wavefunction for further QTAIM analysis. The cluster technique, which takes into account all intermolecular interactions only in the immediate vicinity of the molecule and modern hybrid M062X functional with all-electron basis set 6-311 + G(d,p), yielded NQR parameters that are very well correlated with the experimental data, as in

Table 3. The ¹⁴N NQR parameters for paraxanthine, calculated theoretically at M062X level, assuming clusters of 15 molecules.

	Site	ν+ [MHz]	ν− [MHz]	ν0 [MHz]	e2qQ/h [MHz]	η	s, r2, Intercept, Slope
LT 150 K	>N(1)CH3	2.835	2.511	0.324	−3.564	0.182	0.104
	>N(3)H	2.878	2.113	0.765	−3.327	0.460	0.989
	>N(7)CH3	1.929	1.723	0.206	−2.435	0.169	0.025
	-N(9)=	2.707	2.189	0.517	−3.264	0.317	0.995
RT ABC	>N(1)CH3	2.765	2.517	0.248	−3.521	0.141	0.064
	>N(3)H	2.907	2.225	0.682	−3.421	0.399	0.994
	>N(7)CH3	2.033	1.844	0.189	−2.585	0.146	−0.041
	-N(9)=	2.876	2.440	0.436	−3.544	0.246	1.057
RT DEF	>N(1)CH3	2.774	2.524	0.251	−3.532	0.142	0.058
	>N(3)H	2.910	2.227	0.683	−3.425	0.399	0.994
	>N(7)CH3	2.027	1.840	0.187	−2.578	0.145	−0.033
	-N(9)=	2.873	2.403	0.470	−3.517	0.267	1.052
RT ABCDEF averaged	>N(1)CH3	2.769	2.520	0.250	−3.522	0.142	0.046
	>N(3)H	2.908	2.226	0.683	−3.423	0.415	0.995
	>N(7)CH3	2.030	1.842	0.188	−2.541	0.144	−0.036
	-N(9)=	2.874	2.422	0.453	−3.491	0.265	1.054

and Figures S2 and S3.

The results obtained at the M062X level, assuming the RT ABCDEF averaged model, are very close to those obtained using the solid-state technique at the GGA/RPBE level (calculations for the direct model with weighted proton positions are not feasible using the cluster technique). The bouncing of points off the linear regression ν_M062X(ν_NQR) at low temperature, as in Figure S2, suggests that for =N(9)–, the eQqh⁻¹ should be higher but the asymmetry parameter should be lower. This effect may result from the averaging of the EFG tensor components caused by the rotation of the methyl group at N(7), proton transfer in the N(9) hydrogen bond, or other actors modulating non-covalent bonds.

The rotation of methyl group is not free, and different factors including steric and electronic effects modify the barrier and bring contributions to the structural disorder. Moreover, methyl group surroundings, especially nearby groups, respond to or cooperate with the methyl rotation and modify the barrier. The simulation of the rotation of the methyl group at N(7) permits testing of the influence of these factors. The estimated barrier for the rotation along the threefold axis, as in Figure 4, is 53.95 kJ/mol per cluster, i.e., about 3.60 kJ/mol per single molecule.

The shape of the E(θ) curve is typical of threefold rotational potential. Thus, if an excitation takes place, the staggered conformer would be favored over the eclipsed. The barrier for the methyl group rotation is higher than RT at 296 K (2.479 kJ/mol) and 77 K (0.640 kJ/mol); thus, at room temperature it should exhibit quantum tunnelling enabled by the function overlap. In this temperature range, however, it is more likely that the rotation of the methyl group will be hindered.

Crystal voids in RT ABC (overall volume 14.50 ų, area 192.86 Ų,) are much smaller than those in RT DEF (overall volume 26.93 ų, area 208.16 Ų); their asphericity is higher (0.237 versus 0.192) but their globularity is lower (0.138 versus 0.192). The shape of the crystal voids around the methyl group at N(7) reveals more space around it than around N(1) and opens a path for structural disorder at N(7). It supports the predominance of the θ = 60° conformation, which, in the absence of activation, is more stable. Simulation proves that the rotation of the methyl group at N(7) modulates the EFG tensor, and, consequently, the NQR parameters and resonance frequencies, Figure 5.

Unexpectedly, it affects the eigenvalues and, to a lesser extent, eigenvectors of EFG tensors, at all nitrogens, not only N(7). For nitrogen atoms in six-membered and five-membered rings, two different trends in the e²Qq/h(θ) relationship are evident, as in Figure 5.

As the angle increases, the e²Qq/h value for the N(1) and N(3) sites decreases, but for N(7) and N(9), it increases. In addition, e²Qq/h(θ) of N(9) shows the greatest variability. The e²Qq/h(θ) variation can be relatively easily explained by the pure hyperconjugation, which causes the delocalization of the σ electrons of the C-H bond of the methyl group to the π orbital of the sp²-hybridized nitrogen atom (due to the rotation of the methyl group around the N-C bond, all three C-H bonds of the methyl group are involved in this hyperconjugation). Although e²Qq/h is modulated for all nitrogen sites, the asymmetry parameter changes significantly only for nitrogens N(7) and N(9), which is surprising. This suggests a strong susceptibility of N(9) to electron delocalization associated with hyperconjugation and the influence of additional factors elsewhere. The shape of the e²Qq/h(θ) and η(θ) curves for N(9), similar to E(θ), confirms this conclusion. The orientation of the X, Y, and Z axes of the EFG tensor is different for –N(1)CH₃, >N(3)H, and –N(7)CH₃ than it is for =N(9)–, as in Figure 6.

While in the case of the former orientation, the ZZ-axis is perpendicular to the plane of the molecule, in the case of the latter one it lies in its plane. The XX and YY axes of the EFG tensor for –N(1)CH₃ and –N(7)CH₃ are oriented differently (rotated by an angle of about 30° around the ZZ-axis) and its YY-axis points outside the ring for the former orientation and inside for the latter one. It seems that the rotation of the methyl group should disturb the axes of the EFG tensor at N(9) because its ZZ-axis lies in the plane of the molecule and is directed to the inside of the ring. However, the plot of eQq_xx/h, eQq_yy/h versus eQq_zz/h, as in Figure S4, reveals a completely different nature of the EFG tensor at N(7) compared to that of the other nitrogen sites. This effect is visible both in experimental results and modelling, regardless of the calculation method.

Analysis of the orientation of the main EFG tensor axes upon methyl rotation reveals how they are really affected by the reorientation. The reorientation angles of the EFG tensor axes are small; their values do not exceed 2.4°, as in Figure 7.

The largest axes deviations occur for nitrogen N(7) and concern the XX and YY axes, while the smallest ones occur for N(3). The differences between the orientations of all three EFG axes for the RT ABC and RT DEF conformations are negligible for N(9). Analysis of the resonant frequencies allows assessment of the impact of rotation on the NQR spectrum. The largest frequency shifts are predicted for N(9), but their magnitude is small (5%) and within the FWHM of the resonance line, as in Figure S5. The lack of differences between the individual FWHMs of the experimental NQR resonance lines indicates that the effect of the rotation at each site should be similar or that the shape of the resonance lines does not reveal the presence of rotation. However, the rotation of the methyl group modulates values of the EFG tensor and the positions of the axes, which results in the shifts in all the spectral lines, especially for N(9). This is clearly visible from the calculation results.

There is one more factor to consider, namely the role of the non-covalent bonds in the methyl rotation in the crystal structure. Analysis of these interactions should shed more light on the very different nature of N(7) and N(9) and should be useful in describing the nature of the disorder of the N(7) methyl group.

3.2. Intermolecular Interactions Pattern in Crystal

The strongest non-covalent interactions in crystals are usually hydrogen bonds. According to the Etter [78] rule, the molecule tends to form as many hydrogen bonds as possible. PX can theoretically realize up to twenty-one different types of intermolecular hydrogen bonds (via three donor hydrogen atoms of each methyl group and one hydrogen atom of >N(3)H, and three acceptor atoms N(9) and two O). However, two methyl groups actually limit the ability of nitrogen atoms to participate in strong hydrogen bonds, and, thus, only one strong hydrogen bond (through the >N(3)H group as donor and any of the three acceptor atoms) seem feasible. The essential question is what types of bonds are formed in the crystal structure and what their strengths are. It should help clarify which of them involve nitrogen atoms and to what extent their specific bonding abilities facilitate methyl group disorder.

A combined analysis of 3D Hirshfeld surfaces (3D HS) and 2D molecular fingerprints (2D FP) derived from 3D HS, provided insight into the ten types of the homo- and hetero-nuclear intermolecular interactions in PX crystals, as in

Table 4. Percentage contributions to the 3D Hirshfeld surface area calculated for each pair of species.

Structure	Homonuclear				Heteronuclear
	CC	HH	NN	OO	CH	CN	CO	NH	OH	NO
LT 150 K	3.9	35.5	1.1	0	7.1	6.8	2.1	14.4	26.6	2.6
RT ABC	3.9	36.4	1.6	0	7	6.2	2.4	13.6	26.2	2.6
RT DEF	4	36	1.2	0	6.7	6.3	2.1	14.1	26.7	2.8
RT ABCDEF averaged	3.95	36.2	1.4	0	6.85	6.25	2.25	13.85	26.45	2.7
RT disorder	4	37.8	1.2	0	7.1	6.0	2.2	13.7	25.5	2.6

and Figure 8.

Approximately 83% of all interactions of the complex surface consist of four contact types, namely H⋯H, N⋯H/H⋯N, O⋯H/H⋯O, and C⋯H/H⋯C, all of which involve protons. The distribution of these contacts is much more influenced by temperature (6.92%) than by disorder (0.54%). As the temperature increases, the percentage contributions in PX shift from H⋯C/C⋯H, O⋯H/H⋯O, and N⋯H/H⋯N to H⋯H. In general, for all heavy atoms, differences are within error (0.1%), but those for protons are only slightly higher (0.5%). Thus, the effect of temperature and disorder is weakly manifested (the surface percentage only changes by 0.407% for RT DEF, 0.284% for RT ABC, and 0.296% for the RT disordered structure, respectively).

As much as approximately 60% of the PX molecular surface is generated by H atoms, while the other elements (C, N, and O) contribute to much smaller, but almost equivalent percentages, from about 12 to 16%, as in

Table 5. Enrichment ratios (E_XY) characterizing the various contacts in PX.

Structure	Atom	C	H	N	O
LT 150 K	Surface%	11.9	59.55	13.0	15.65
	C	2.75	-	-	-
	H	0.50	1.00	-	-
	N	2.20	0.93	0.65	-
	O	0.52	1.43	0.64	0
RT ABC	Surface%	11.8	59.8	12.8	15.6
	C	2.87
	H	0.50	1.02
	N	2.05	0.89	0.98
	O	0.65	1.40	0.65	0
RT DEF	Surface%	11.55	59.75	12.8	15.8
	C	3.00
	H	0.49	1.01
	N	2.13	0.92	0.73
	O	0.58	1.41	0.69	0
RT ABCDEF averaged	Surface%	11.68	59.78	12.80	15.70
	C	2.94	-	-	-
	H	0.50	1.02	-	-
	N	2.09	0.91	0.86	-
	O	0.62	1.41	0.67	0
RT disordered	Surface%	11.65	60.95	12.35	15.15
	C	2.95	-	-	-
	H	0.50	1.02	-	-
	N	2.09	0.91	0.79	-
	O	0.62	1.38	0.69	0

. The values of the enrichment ratio of the intermolecular contacts, E_XY, which reveal privileged (E_XY > 1) and disfavored (E_XY < 1) contacts between each of the two atomic species, X and Y, are collected in Table 5 and shown in Figure 9. They increase with temperature by only 0.218%, which mainly reflects the reorientation of the -CH₃ group.

The five C and four N atoms in the two heterocyclic rings are exposed to the molecular surface because they each have three (C) and two or three (N) bonds (trigonal planar sp² hybridization) and their p orbitals participate in π bonds above and below the ring. The C⋯C contacts in PX are significantly enriched (E_CN ≈ 2.75); therefore, C species seem the preferred interaction partners for C, but are less privileged compared to O (E_CO ≈ 1.27) or N (E_CN ≈ 1.34). These three types of contacts have enrichment ratios much higher than unity, suggesting their high contribution to π⋯π stacking. The decrease in temperature (i.e., change from RT DEF to LT) reduces the E_CC enrichment factor by only 0.19%, which, while maintaining the E_CO and E_CN values, suggests a weakening of π⋯π stacking between carbon atoms. However, the contribution of C atoms to the surface is small (11%), and the percentage of C⋯C to the 3D Hirshfeld surface area is only 4%, so the role of C⋯C is overestimated. C⋯H contacts are not favored, their enrichment ratio is low, and E_CH ≈ 0.50, because H atoms participate in the formation of hydrogen bonds. The H⋯H contacts appear with an enrichment ratio equal to unity (E_HH ≈ 1.0), despite their strong competition with O and N hydrogen acceptors (H⋯O, E_OH ≈ 1.22 and H⋯N, E_NH ≈ 1.18). They are not impoverished due to the presence of high hydrogen content. The O⋯H/H⋯O contacts exhibit high enrichment values (E_OH ≈ 1.22), because they are favored in crystal packing. Two oxygen atoms form intermolecular bonds of the type N–H⋯O (–NH from –NH) and C–H⋯O. For the same reason, the O⋯O contacts are impoverished (E_OO ≈ 0.57). The contacts N⋯H/H⋯N show enrichment values higher than unity (E_NH ≈ 1.18). Thus, the existence of the C–H⋯N hydrogen bond is raised, and the N⋯N contacts are impoverished (E_NN ≈ 0.70). Overall, the favored interaction partners for N, whose contribution to the 3D HS molecular surface is of about 14.1%, are only C and H species, while N and O are disfavored.

The percentage contributions to local 3D HS were calculated for a single species, and four nitrogen atoms, two oxygen atoms, two carbons of methyl groups, and two methyl groups are listed in

Table 6. Percentage contributions to the local 3D Hirshfeld surface area calculated for a single species.

Structure	Contact Type	–N(1)=	–N(3)=	–N(7)	N(9)	Contact Type	=O(4)	=O(2)	Contact Type	C (N(1) CH3)	C N(7) CH3	Contact Type	CH3 (N(1))	CH3 (N(7))
LT 150 K	N⋯C	88.2	48.7	79.6	50.4	O⋯C	30.5	24.9	C⋯C	0	0	H⋯C	6.8	0.7
	N⋯H	11.8	49.5	20.0	48.0	O⋯H	65.7	75.1	C⋯H	88.9	87.9	H⋯H	42.9	57.1
	N⋯O	0	1.8	0	1.6	O⋯O	0	0	C⋯O	0	0.5	H⋯O	21.4	21.4
	N⋯N	0	0	0.4	0	O⋯N	3.8	0	C⋯N	11.1	11.7	H⋯N	21.5	13.5
RT ABC	N⋯C	87.3	48.6	79.9	50.8	O⋯C	29.4	24.7	C⋯C	0	0	H⋯C	6.3	0.2
	N⋯H	12.7	50	19.9	47.8	O⋯H	67	75	C⋯H	89.2	88.6	H⋯H	44.8	55.4
	N⋯O	0	1.3	0	1.4	O⋯O	0	0	C⋯O	0	0.2	H⋯O	20.5	22.4
	N⋯N	0	0.1	0.1	0	O⋯N	3.6	0.3	C⋯N	10.8	11.2	H⋯N	21.1	13.6
RT DEF	N⋯C	87.3	48.6	79.6	50.3	O⋯C	30.4	24.8	C⋯C	0	0	H⋯C	6.4	0.9
	N⋯H	12.7	49.7	20.4	48.3	O⋯H	65.6	74.8	C⋯H	89.1	88.1	H⋯H	43.7	56.5
	N⋯O	0	1.7	0	1.4	O⋯O	0	0.4	C⋯O	0	0.3	H⋯O	20.8	22.1
	N⋯N	0	0	0	0	O⋯N	4	0	C⋯N	10.9	11.7	H⋯N	21.8	13.1
RT ABCDEF averaged	N⋯C	87.3	48.6	79.8	50.6	O⋯C	29.9	24.8	C⋯C	0.0	0.0	H⋯C	6.4	0.6
	N⋯H	12.7	49.9	20.2	48.1	O⋯H	66.3	74.9	C⋯H	89.2	88.4	H⋯H	44.3	56.0
	N⋯O	0.0	1.5	0.0	1.4	O⋯O	0.0	0.2	C⋯O	0.0	0.3	H⋯O	20.7	22.3
	N⋯N	0.0	0.1	0.1	0.0	O⋯N	3.8	0.2	C⋯N	10.9	11.5	H⋯N	21.5	13.4

.

Temperature has little effect on these percentages. The difference between the percentage contributions to the local 3D Hirshfeld surface area from nitrogen atoms does not exceed 1.404%. At both oxygen sites, these percentages at LT and RT differ by only 0.917%, while at the carbons from methyl groups at N(1) and N(7) they differ by 0.678%. Surprisingly, the difference in the percentage contributions to the local 3D Hirshfeld surface calculated for N(9) and N(3) is insignificant, even though N(3) is covalently bound to H, while N(9) only accepts a proton and forms a N–H⋯O hydrogen bond. This confirms that the N–H⋯O hydrogen bond is strong and thermally stable. The contributions of the methyl groups substituted at N(1) and N(7) differ by 2.067% and temperature influence them in 2.358%. Thus, the main changes in the interactions pattern are rather minor, involving the methyl groups and reflecting their disorder.

A little more information about the types of interactions and their location is provided by 3D HS results, which, with the normalized contact distance d_norm, shape index, and curvedness mapped over this surface, are shown in Figure 10.

The 3D HS surfaces are sufficient to identify most of the forces stabilizing the crystal structure. However, the shape index and curvedness maps are much more useful for identifying the presence of π–π stacking or intramolecular interactions. The red/blue triangular patches in the shape index and the flat patches in the curvature are essentially the only reliable methods of identifying the presence of these interactions, which are very faint on 3D HS normal surfaces.

Intermolecular hydrogen bonds are indicated in the 3D HS surface by intense red to white areas, as in Figure 10. The relevant parameters, namely d_norm, shape index, and curvedness, which describe these hydrogen bonds, are listed in

Table 7. The parameters d_norm, shape index, and curvedness characterizing the intermolecular hydrogen bonds in paraxanthine.

Hydrogen Bond	Parameter	LT 150 K	RT ABC	RT DEF
N(3)–H⋯O(2)	dnorm	−0.6156	−0.6012	−0.5994
	Shape index	−0.9788	−0.9781	−0.9800
	Curvedness	−1.8587	−1.6237	−1.6267
C(8)–H⋯O(6)	dnorm	−0.1792	−0.1258	−0.1200
	Shape index	−0.9536	−0.8976	−0.9133
	Curvedness	−1.8293	−1.8489	−1.8332
C–H⋯O(2)	dnorm	0.1752	0.3293	0.2277
	Shape index	−0.1031	−0.2422	−0.6734
	Curvedness	−2.6470	−2.5477	−2.7226
C–H⋯N(9)	dnorm	−0.019	−0.052	−0.098
	Shape index	−0.8637	−0.8101	−0.4796
	Curvedness	−1.8697	−1.8412	−2.0382
C–H⋯O(6)	dnorm	−0.0753	−0.0569	−0.0612
	Shape index	−0.9330	−0.7467	−0.9141
	Curvedness	−1.4493	−1.3704	−1.3276

.

Intense red areas in the 3D HS located near =O(2) and H (from –N(3)H) represent the lowest value of d_norm = −0.6156 a.u. (−0.6012 a.u. for RT ABC; −0.5994 a.u. for RT DEF), which confirms the presence of a pair of strong N(3)–H⋯O(2) hydrogen bonds forming a dimer. These bonds are short R_O⋯N = 2.809 Å (2.820 Å for RT ABC and RT DEF) and only slightly nonlinear <NHO = 171.1° (167° for RT ABC and RT DEF). The brighter red area near –O(6)= and H from –C(8)H moiety represents the d_norm value of −0.1792 a.u. (−0.1258 a.u. for RT ABC; −0.120 a.u. for RT DEF). It indicates the presence of a weak C(8)–H⋯O(6) hydrogen bond, that is long, R_C⋯O = 3.439 Å (3.511 Å for RT ABC and RT DEF), and highly nonlinear, <CHO = 175.53° (176.35° for RT ABC and RT DEF). The light pink area near –O(2)= and H (from –N(1)CH₃ moiety) represents the value of d_norm 0.1752 a.u. (0.3293 a.u. for RT ABC and 0.2277 a.u. for RT DEF) and indicates the presence of a weak C–H⋯O(2) hydrogen bond, that is long, RC⋯O = 3.406 Å (3.504 Å for RT ABC and RT DEF), and highly nonlinear <CHO = 119.56° (122.37° for RT ABC and RT DEF).

The almost white area near –N(9)= and H (from –N(1)CH₃) moiety represents the value of d_norm 0.019 a.u. (0.052 a.u. for RT ABC; 0.098 a.u. for RT DEF). It indicates the presence of a weak C–H⋯N(9)) hydrogen bond, which is long RC⋯N = 3.390 Å (3.357 Å for RT ABC and RT DEF) and highly nonlinear <CHO = 138.26° (127.3° for RT ABC and RT DEF). The light pink areas near O(6) and H (from –N(7)CH₃) represent the value of d_norm −0.0753 a.u. (−0.0569 a.u. for RT ABC; −0.0612 a.u. for RT DEF). They imply the presence of a pair of weak C–H⋯O(6) hydrogen bonds (C–H from –N(7)CH₃), which are long RC⋯O = 3.0916 Å (3.116 Å RT for ABC and RT DEF) and very highly nonlinear <CHO = 90.03° (91.38° for RT ABC and RT DEF).

The white areas in 3D HS describe the dispersive contacts C–H⋯H–C mainly between the protons from neighboring N(7)CH₃ groups. Their two pairs (R_HH = 2.300 Å and R_CC = 3.543 Å at LT, and R_HH = 2.324 Å and R_CC = 3.572 Å RT DEF) form a dimer. The remaining proton of N(7)CH₃ participates in dispersive contact C–H⋯H–C (R_HH = 2.135 Å and R_CC = 3.572 Å) with a proton from N(1)CH₃. There are fewer dispersive contacts in the RT ABC structure—only one proton of N(7)CH₃ participates in the H⋯H of only 1.655 Å. The H⋯H contacts in each case are shorter than twice the van der Waals radius of hydrogen (2.40 Å), which suggests that dispersion forces involve both H⋯H and C⋯H. For some compounds, dispersion forces involving C⋯H may be even stronger than those involving H⋯H, but such an arrangement in PX is precluded by the E_XY enrichment ratios, as in Table 5.

Intramolecular interactions are also present but difficult to find. They can only be detected on the local 3D HS surface with the curvature mapped on this surface generated for the acceptor (O(6)). In RT ABC, both methyl groups at N(1) and N(7) seem to be involved in the intramolecular hydrogen bonds CH⋯O, while in RT DEF only methyl substituted at N(7) seems to participate in such a bond.

Multiple interatomic contacts contribute to the interlayered π⋯π stacking, such as N(7)⋯N(3) of 3.415 Å (3.463 Å for RT ABC and RT DEF); C(8)⋯N(1) of 3.361 Å (3.396 Å for RT ABC and RT DEF); C(6)⋯N(9) of 3.352 Å (3.399 Å for RT ABC and RT DEF); C(4)⋯C(5) of 3.378 Å (3.420 Å for RT ABC and RT DEF) and C(2)⋯N(7) of 3.441 Å (3.439 Å for RT ABC and RT DEF), with the distance between layers of 3.351 Å. The displacement of the molecules relative to each other in adjacent layers weakens the π⋯π stacking. Therefore, there are only two types of interatomic contacts that contribute to the π⋯π stacking: N(1)⋯C(6) of 3.419 Å (3.497 Å RT ABC) and N(7)⋯C (C from N(1)CH₃) of 3.440 Å (3.481 Å RT ABC) with an interlayer distance of 3.436 Å. All these π⋯π stacking interactions are visible only for the shape index mapped over 3D HS surface as characteristic mixed red/blue triangular patches.

In order to verify the presence of interactions, mainly directional, the characteristic features of the 2D fingerprint (2D FP), obtained from the 3D HS, were further analyzed, as in Figure 11.

In general, the 2D FP results for LT and RT DEF are in 98.9% similar, while those for RT ABC and the disordered structure differ by as much as 27.7%, mainly in the central part of the fingerprint and along the axis of symmetry i.e., d_e ≅ d_i. In the 2D FP plot in Figure 11b, the O⋯H/H⋯O contacts describing intramolecular C–H⋯O and intermolecular N–H⋯O hydrogen bonds are represented by symmetric spikes (‘‘wings’’), which are very sharp and long at d_e + d_i ~1.8 Å (d_e + d_i ~1.75 Å for the RT ABC and DEF). These interactions cover an area of 26.6% (26.2% for the RT ABC, 26.7% for RF DEF, 26.45% for the RT averaged, and only 25.5% for the disordered structures) of the total 3D HS. The N⋯H/H⋯N contacts are represented by the most external symmetric spikes consisting of two components, which are wide at d_e + d_i ~2.7 Å and sharp at d_e + d_i ~2.4 Å (d_e + d_i ~2.0 Å for the optimised geometry), as in Figure 11c. These interactions cover a smaller area of 14.4% (13.6% for the RT ABC, 14.1% for RF DEF, 13.85% for the RT averaged, and 13.7% for the disordered structures) of the total 3D HS. The individual 2D FPs differ significantly in a wide spike in the central region that describes the dispersion interactions: a single H⋯H at d_e + d_i ~1.6 A in RT ABC and disordered structure and a pair of H⋯H at d_e + d_i ~2.3 Å in RT DEF and LT, as in Figure 11d. This indicates significantly stronger H⋯H interactions in RT ABC compared to LT and RT DEF.

The C⋯H/H⋯C contacts are represented in 2D FP by the widest spikes at d_e + d_i ~3.00 Å, as in Figure 11e. They cover 7.1% (7% for the RT ABC, 6.7% for RF DEF, 6.85% for the RT averaged, and 7.1% for the disordered structures) of the total 3D HS. These contacts are significantly weaker than those involving C–H⋯H–C.

The C⋯O/O⋯C contacts bring a very small contribution of 2.1% (2.4% for the RT ABC, 2.1% for RF DEF, 2.25% for the RT averaged, and 2.2% for the disordered structures) and are represented by wide wings d_e + d_i ~3.40 Å (3.10 Å for the optimized geometry) located in the middle area of the entire fingerprint. The O⋯N/N⋯O contacts bring a small contribution of 2.6 (2.6% for the RT ABC, 2.8% for RF DEF, 2.7% for the RT averaged, and 2.6% for the disordered structures). The high contribution is brought by weak H⋯H interactions, which cover 35.5% (36.4% for the RT ABC, 36% for RF DEF, 36.2% for the RT averaged, and 37.8% for the disordered structures) of the total 3D HS. In the 2D FP, they are reflected by the cloud of scattered points and wide joined spikes at d_e + d_i ~2.0 Å. The contribution brought by C⋯N/N⋯C and N⋯N contacts, which covers 6.8% (6.2% for the RT ABC, 6.3% for RF DEF, 6.25% for the RT averaged, and 6% for the disordered structures) and 1.1% (1.6% for the RT ABC, 1.2% for RF DEF, 1.4% for the RT averaged, and 1.2% for the disordered structures), as in Figure 11f, mainly describes the interlayer π⋯π stacking interactions. The contribution brought by C⋯C contacts covers 3.9% (4.0% for the remaining structures), as in Figure 11g. Carbon, being less electronegative than nitrogen, has a lower capability to stabilize valence electron density. Hence, it exhibits greater reactivity and plays a key role in π-stacking.

Local 2D molecular fingerprints of –N(1)CH₃, >N(3)H, –N(7)CH₃, and –N(9)=, generated on 3D HS for individual nitrogen atoms, as in Figure S6, vary from one nitrogen site to another, both in 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 and type of non-covalent interactions.

The local 2D FPs suggests that N(1) and N(7) participate only in covalent NC bonds. The >N(3)H also participates in such bonds but additionally to a covalent NH bond, which is a donor for strong hydrogen bond N–H⋯O. Nitrogen –N(9)=, apart from NC, participates in the intermolecular C–H⋯N hydrogen bond. These two nitrogen atoms, N(3) and N(9), are crucial to the bonding pattern in the PX structure. A deeper insight into local FPs reveals a small, easily overlooked wing that indicates the involvement of N(7) in the C–H⋯N bonding between adjacent molecules. It supports the π-π stacking, which can, therefore, be stronger than average.

The local 2D FP for the oxygen atom O(2) reveals two sharp wings, with the upper one attributed to the OC covalent bond at d_e + d_i = 1.3 Å and the bottom one attributed to two hydrogen bonds: N–H⋯O at d_e + d_i = 1.85 Å and C–H⋯O at d_e + d_i = 2.7 Å. The analysis of the local 2D FP for the oxygen atom O(6) 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 (weaker intermolecular C–H⋯O at d_e + d_i = 2.5 Å and stronger intramolecular C–H⋯O at d_e + d_i = 2.2 Å). The wing describing the intramolecular hydrogen bond C–H⋯O differs significantly from those describing the intermolecular hydrogen bonds. It is specific, and very wide and scattered. The analysis of local, atomic, 2D molecular fingerprints helps to find the very weak intramolecular hydrogen bonds C–H⋯O, which meet the Gilli criterion [79].

The examination of changes in global 2P FP caused by the rotation of the methyl group in N(7) reveals changes in the distribution and strength of non-covalent interactions in the entire molecule. They can be quantified using the Pompeiu–Hausdorff and Bhattacharayya metrics and the Wasserstein distance, as in Figure 12.

The Pompeiu–Hausdorff distance between 2D FP upon methyl rotation does not exceed 0.5% (the global minimum and maximum, at θ = 0° and θ = 60°, respectively). The Wasserstein distance, which allows identification of the dissimilarity between the 2D FP distributions, does not exceed 2% (the global minimum and maximum, at θ = 0° and θ = 60°, respectively). Both methods describe, quite well, the differences in the distribution of interactions, which occupy mainly the central region of the fingerprint (d_i = d_e axis). The Bhattacharayya distance, which allows identification of differences in the strength of the interactions, does not exceed 6% (the global minimum and maximum, at θ = 0° and θ = 60°, respectively). The distribution of C-H⋯H-C contacts is clearly different for individual conformers, as in Figure 10d. The 2D FP changes due to methyl rotation are for all contacts taken together and only C-H⋯H-C contacts of comparable magnitude. The differences between them do not exceed 1%, and the global extremes occur at the same abscissa values. Thus, among the intermolecular interactions, the C-H⋯H-C contacts are mainly affected during the rotation of the methyl group at N(7). Moreover, the 2D FP(θ) dependence suggests that the change in the number of interactions is smaller than the change in the strength of interactions. The magnitude of the percentages, as in Figure 12, is in good agreement with the modulation of the NQR parameters due to rotation, as in Figure 5 and Figure 7, as well as the changes of the total energy upon the rotation of the methyl group, as in Figure 4.

3.3. Characterization of the Strength of the Interactions

The packing similarity between the LT and RT DEF crystal structures is high (99.75%), while that between the LT and RT ABC structures is lower (99.68%). 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, as in

Table 8. Separation of electrostatic, E_e, polarization. E_p, dispersion, E_d, and repulsion, and E_r, terms for the total lattice energy, E_t, calculated with the Spackman scaled approach (optimized proton positions).

Structure	Ee [kJ/mol]	Ep [kJ/mol]	Ed [kJ/mol]	Er [kJ/mol]	Etotal [kJ/mol]
LT 150 K	−105.00	−27.00	−177.80	160.30	−186.65
RT ABC	−102.15	−25.75	−168.45	154.75	−178.35
RT DEF	−99.20	−25.45	−169.65	141.65	−183.95
RT averaged	−100.68	−25.60	−169.05	148.20	−181.15

, has revealed that irrespective of the structure (ordered or disordered) and temperature (LT or RT), the attractive interactions contribute about 72% of the total energy and cannot be cancelled by the much smaller repulsion of about 38%.

At 150 K, the total lattice energy (−105 kJ/mol) is dominated by the dispersion (−177.8 kJ/mol) and electrostatic forces (−105.82 kJ/mol), which are stronger than the polarization (−27 kJ/mol) or repulsion (56.79 kJ/mol). With increasing temperature, all the contributions decrease in absolute value, resulting in a slightly lower total energy, as in Figure 13.

The particular energy terms in all these structures are strongly correlated (the Pearson coefficient exceeds 0.999), Figure S7, but the scattering is moderate for LT versus RT DEF (6.64%), high for LT versus RT ABC (20.63%), and marginal for the LT versus RT averaged (2.79%).

This suggests that LT and RT DEF are “energetically” more similar than LT and RT ABC. The two types of leading interactions, electrostatic (mainly hydrogen bonds; shown in red) and dispersion (mainly π⋯π; shown in green), form networks of non-covalent interactions in two orthogonal directions, as in Figure 14. An increase in temperature only slightly weakens the total energy of the lattice (by 4.06 kJ/mol), and the proportions between individual terms are within 1%, so the temperature slightly affects the shape of this energy skeleton.

Closer analysis of the data reveals that the major contributions to the total lattice energy, as in

Table 9. The strength of the particular interactions in paraxanthine (optimized proton positions), estimated using the pairs of molecules (electrostatic, E_e, polarization. E_p, dispersion, E_d, and repulsion, E_r, terms for the total lattice energy, E_t).

Structure	Interaction Type and Their Number		Ee [kJ/mol]	Ep [kJ/mol]	Ed [kJ/mol]	Er [kJ/mol]	Et [kJ/mol]
LT 150 K	N(3)H⋯O(2)	2	−87.7	−19.6	−14.3	92.6	−62.4
	C(8)–H⋯O(6)	1	−10.8	−3.0	−11.6	14.5	−14.8
	(N7)C–H⋯O(6) + N(7)C–H⋯H–CN(7)	2	−17.3	−5.4	−16.7	16.0	−26.9
	N(9)⋯HC(8) v N(9)⋯H–CN(7)	1	−4.7	−2.2	−11.3	7.2	−11.9
	N(3)CH⋯N(9)	1	−6.3	−1.8	−10.8	11.5	−10.2
	N(1)CH⋯O(6) + N(1)CH⋯HCN(1)	1 + 1	−4.4	−1.8	−7.2	3.9	−9.8
	π⋯π stacking	1	−8.4	−2.8	−63.2	35.2	−44.3
	π⋯π stacking	1	−17.9	−2.9	−58.2	33.7	−51.0
RT ABC	N(3)H⋯O(2)	2	−84.5	−18.8	−13.9	86.0	−62.3
	C(8)–H⋯O(6)	1	−8.4	−2.2	−10.1	8.3	−14.2
	(N7)C–H⋯O(6) +N(7)C–H⋯H–CN(7)	2	−17.6	−5.7	−15.8	33.3	−16
	N(9)⋯HC(8) v N(9)⋯H–CN(7)	1	−7.3	−2.8	−11.5	14.4	−10.9
	N(3)CH⋯N(9)	1	−4.5	−1.5	−9.9	7.9	−9.6
	N(1)CH⋯O(6) + N(1)CH⋯HCN(1)	1 + 1	−4.2	−1.6	−6.5	3.1	−9.4
	π⋯π stacking	1	−12.5	−3.0	−61.2	32.4	−48.8
	π⋯π stacking	1	−14.2	−2.4	−54.4	29.0	−46.3
RT DEF	N(3)H⋯O(2)	2	−84.3	−18.8	−13.9	85.6	−62.3
	C(8)–H⋯O(6)	1	−9.2	−2.6	−11.0	10.6	−14.6
	(N7)C–H⋯O(6) + N(7)C–H⋯H–CN(7)	2	−17.7	−5.3	−16.2	14.7	−27.8
	N(9)⋯HC(8) v N(9)⋯H–CN(7)	1	−4.0	−2.1	−10.6	6.6	−10.9
	N(3)CH⋯N(9)	1	−4.5	−1.5	−9.8	7.9	−9.5
	N(1)CH⋯O(6) + N(1)CH⋯HCN(1)	2 + 1	−4.1	−1.6	−6.5	3.1	−9.3
	π⋯π stacking	1	−8.6	−2.7	−61.1	30.9	−45.2
	π⋯π stacking	1	−17.8	−2.9	−55.6	32.4	−49.4

and Figure S8, come from a pair of strongly concerted N–H⋯O hydrogen bonds (about 31.2 kJ/mol) and a pair of (N7)C–H⋯O(4) bonds (about 13 kJ/mol for LT, RT DEF structures, and 8 kJ/mol for RT ABC structure).

Weak hydrogen bonds C(8)–H⋯O(6) of 14.2 kJ/mol, (N9)C–H⋯C(8), and (N3)C–H⋯(N9) of about 10 kJ/mol and N(1)CH⋯O(6) of about 9 kJ/mol bind the PX molecules together. In addition, non-polar N(1)C–H⋯H–CN(1) dispersive contacts of about 9 kJ/mol bind protons of methyl groups from neighboring molecules. The entire structure is stabilized by the interlayered π⋯π stackings of about 40–50 kJ/mol. The ordering of the most important interactions according to their decreasing energy is as follows: a pair of N–H⋯O > π⋯π stacking > C–H⋯O > C–H⋯N > C–H⋯H–C. Neither temperature nor disorder change it. The N(3)–H⋯O(2) hydrogen bonds and π⋯π stacking are the main driving forces forming the “energetical” skeleton in PX, as shown in Figure 14. However, a significant difference occurs in the (N7)C–H⋯O(6) hydrogen bond and N(7)C–H⋯H–CN(7) dispersive contacts, which are much weaker in RT DEF than in LT or RT ABC. These dihydrogen contacts are subtle but definitely not weak. Apparently, the change in the conformation of the methyl group at N(7) affects the strength of the (N7)C–H⋯O(6) hydrogen bond and dispersive contact N(7)C–H⋯H–CN(7). The change from syn to anti causes their weakening by up to 5.5 kJ/mol, Table 8, but only the cumulative effect of (N7)C–H⋯O(6) and N(7)C–H⋯H–CN(7)) can be assessed. Conformational disorder is clearly visible in Figure S8. The points scattering from E_RT(E_LT) correspond to repulsion and, to a lesser extent, to electrostatic terms for (N7)C–H⋯O(6) and N(7)C–H⋯H–CN(7), and, in a minor degree, to N(9)⋯H–CN(7). Thus, conformational disorder significantly modulates these two contributions to the total energy.

Although the above analysis provides a preliminary and rough insight into the distribution of the interactions, it is still incomplete. Therefore, the quantum theory of atoms in molecules (QTAIM) was used to verify the presence and estimate the strength of the inter- and intra-molecular interactions. The topological descriptors of electron density (ρ(r), Δρ(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 S3. The molecular graphs based on the crystalline structure are shown in Figure 15. Overall, the QTAIM analysis confirmed the presence of the intermolecular hydrogen bonds and H⋯H contacts which contribute to the stabilization of the crystalline structure, as detected by the total 2D FP, as in

Table 10. The energies of the intermolecular interactions calculated according to the Espinosa (E_E), Matta (E_M), Emamian (E_EM), and Afonin (E_A) equations at the M062X/6-311 + G** level.

Structure	Type	Proton Donor Group	Nature	Contact	EE [kJ/mol]	EM [kJ/mol]	EEM [kJ/mol]	EA [kJ/mol]	EN [kJ/mol]
	intra	N(1)CH3	HB, electrostatic	CH⋯O(6)	−20.9	−20.54	−15.02	−12.03	-
LT 150 K	inter	N(3)H	HB, electrostatic	N(3)H⋯O(2)	−33.71	−30.20	−27.45	−19.12	−22.32
	inter	N(3)H	HB, electrostatic	N(3)H⋯O(2)	−33.47	−29.80	−27.70	−18.99	−22.58
	inter	N(1)CH3	HB, electrostatic	CH⋯N(9)	−9.15	−8.73	−7.29	−5.52	-
	inter	N(1)CH3	HB, electrostatic	CH⋯N(9)	−9.16	−8.73	−7.29	−5.52	-
	contact	-	van der Waals	N(7)⋯O(6)	−5.71	−6.16	−2.73	−3.61	-
	inter	N(7)CH3	HB, electrostatic	CH⋯O(6)	−5.63	−6.12	−2.69	−3.57	-
	inter	N(7)CH3	Dispersive	CH⋯HC	−3.74	−4.84	−2.43	−2.52	-
	inter	N(7)CH3	Dispersive	CH⋯HC	−3.73	−4.83	−2.42	−2.52	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(2)	−4.56	−5.04	−2.11	−2.97	-
	inter	N(7)CH3	HB, electrostatic	CH⋯N(9)	−4.29	−4.71	−2.34	−2.83	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(6)	−2.27	−2.83	0.08	−1.71	-
	inter	N(1)CH3	Dispersive	CH⋯HC	−1.09	−1.60	1.06	−1.05	-
	inter	C(8)H	HB, electrostatic	C(8)H⋯O(6)	−10.60	−9.32	−7.12	−6.32	-
RT ABC	intra	N(1)CH3	HB, electrostatic	CH⋯O(6)	−21.46	−21.04	−15.61	−12.34	-
	intra	N(7)CH3	HB, electrostatic	CH⋯O(6)	−12.41	−11.20	−9.00	−7.32	-
	inter	N(3)H	HB, electrostatic	N(3)H⋯O(2)	−32.23	−28.79	−25.94	−18.31	−20.80
	inter	N(3)H	HB, electrostatic	N(3)H⋯O(2)	−32.00	−28.43	−26.17	−18.18	−21.04
	inter	N(7)CH3	Dispersive	CH⋯HC	−20.54	−17.24	−15.74	−11.83	-
	inter	N(7)CH3	HB, electrostatic	CH⋯N(9)	−9.65	−9.18	−7.30	−5.80	-
	inter	N(1)CH3	HB, electrostatic	CH⋯N(9)	−6.66	−6.76	−4.98	−4.14	-
	inter	N(1)CH3	HB, electrostatic	CH⋯N(9)	−6.66	−6.76	−4.98	−4.14	-
	inter	N(7)CH3	HB, electrostatic	CH⋯O(6)	−5.52	−6.02	−2.57	−3.51	-
	inter	N(7)CH3	HB, electrostatic	CH⋯O(6)	−5.43	−5.97	−2.51	−3.46	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(2)	−3.73	−4.27	−1.32	−2.52	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(6)	−1.88	−2.45	0.48	−1.49	-
	inter	N(1)CH3	Dispersive	CH⋯HC	−1.16	−1.74	0.90	−1.09	-
	inter	C(8)H	HB, electrostatic	C(8)H⋯O(6)	−8.·90	−8.06	−5.75	−5.38	-
RT DEF	intra	N(1)CH3	HB, electrostatic	CH⋯O(6)	−21.52	−21.07	−15.64	−12.37	-
	inter	N(3)H	HB, electrostatic	N(3)H⋯O(2)	−32.15	−28.71	−25.84	−18.26	−20.70
	inter	N(3)H	HB, electrostatic	N(3)H⋯O(2)	−31.92	−28.36	−26.07	−18.13	−20.93
	inter	N(1)CH3	HB, electrostatic	CH⋯N(9)	−6.66	−6.77	−4.98	−4.14	-
	inter	N(1)CH3	HB, electrostatic	CH⋯N(9)	−6.66	−6.77	−4.98	−4.14	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(6)	−5.42	−5.91	−2.48	−3.45	-
	inter	N(7)CH3	HB, electrostatic	CH⋯O(6)	−5.35	−5.87	−2.43	−3.41	-
	inter	N(7)CH3	Dispersive	CH⋯HC	−3.50	−4.56	−2.15	−2.39	-
	inter	N(7)CH3	Dispersive	CH⋯HC	−3.51	−4.56	−2.15	−2.40	-
	inter	N(7)CH3	HB, electrostatic	CH⋯N(9)	−4.40	−4.81	−2.45	−2.89	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(2)	−3.73	−4.27	−1.32	−2.52	-
	inter	N(1)CH3	HB, electrostatic	CH⋯O(2)	−1.90	−2.47	0.46	−1.50	-
	inter	C(8)H	HB, electrostatic	C(8)H⋯O(6)	−8.90	−8.06	−5.75	−5.38	-

. For each of the interactions listed in Table S3, the values of ρ(r) in the BCP do not exceed 0.1 a.u., while the values of the Laplacian, Δρ_BCP(r), are positive. Thus, according to the Koch and Poplier [80] criterion, all these interactions can be defined as non-covalent. Moreover, in each case the –G(r_BCP)/V(r_BCP) ratio exceeds 1; thus, all of them are truly non-covalent (i.e., none are partially covalent). The energies of the N(3)H⋯O(2) and C(8)H⋯O(6) hydrogen bonds, estimated using QTAIM and Equations (7)–(11), as in Table 10 and Table S3, are in good agreement with those calculated using the classical energy framework approach, but the energies of the remaining non-covalent bonds are much lower. Some contributions are combined, and their energies were misestimated, while others were missed. However, the energies of the C-H⋯H-C dispersive interactions may be underestimated because in QTAIM these bonds are treated individually and the coupling effect between them is ignored. Regardless, these results confirm the conclusion that LT and RT DEF are more similar than LT and RT ABC.

In contrast to the 2D FP in which intramolecular hydrogen bonds are indistinguishable from intermolecular ones, or simply undetectable, QTAIM detected intramolecular N(1)C–H⋯O(6) hydrogen bond in all LT and RT structures, but N(7)C–H⋯O(6) only in the RT ABC structure. However, in some cases, the standard QTAIM approach is also insufficient to detect all intramolecular interactions, so the isosurface of the reduced density gradient, RDG(r), was investigated. The reduced density gradient isosurfaces with sign(λ₂)ρ_BCP mapped over the surface are shown in Figure 15.

The light cyan slabs between N(3)H and O(2) indicate the presence of stabilizing intermolecular N(3)H⋯O(2) hydrogen bonds. Smaller and green slabs between C(8)H and O(6) suggests the existence of much weaker intermolecular C(8)H⋯O(6) hydrogen bonds. The disc-shaped green surfaces of the reduced density gradient located between the protons from two –N(7)CH₃ groups belonging to the neighboring molecules confirm the presence of the dispersive C–H⋯H–C interactions, as in Figure 15a–c. The small green slab that is disc-shaped, located between the proton from –N(1)CH₃ and O(6), indicates the presence of the weak attractive intramolecular C–H⋯O hydrogen bond, as in Figure 15d,e. Another, slightly larger green slab, that is also disc-shaped, is located between the proton from –N(7)CH₃ and O(6), and confirms our assumptions about the presence of the weak attractive intramolecular C–H⋯O hydrogen bond in RT ABC structure, as in Figure 15e, even though no critical point was detected. The tiny green area between the proton of –N(1)CH₃ and O(2) suggests the existence of a very weak attractive intramolecular C–H⋯O hydrogen bond, but only in the LT structure Figure 15d. Both interactions affect the barrier height to the rotation of the methyl group at the N(7) position, which leads to hindered rotation and structural disorder. As illustrated in Figure 15c–e, one or two dark red elliptical slabs are present at the center of the six-membered ring and five-membered ring, respectively. Both indicate a significant classical steric effect (repulsion) that results from the spatial arrangement of N and C atoms in the rings. This steric crowding is stronger in the six-membered ring than in the five-membered ring.

In the RDG(r) versus sign(λ₂)ρ(r) plot, as in Figure 16, the non-covalent interactions, i.e., hydrogen bonds (strong and attractive), van der Waals contacts (weak), and steric interactions (ring closures, repulsive) are visible as spikes (on the left, λ₂ > 0, strong/weak repulsive, while on the right, λ₂ < 0, weak/strong attractive). The spikes on the left correspond to negative sign(λ₂)ρ and the most negative value −0.04 a.u. The spikes on the right correspond to positive sign(λ₂)ρ, where the least positive value is 0.055 a.u. Such a wide range of sign(λ₂)ρ suggests a combination of strong hydrogen bonding, van der Waals, π⋯π stacking, and steric interactions contributing to the stability of the crystal. Moreover, the plot indicates that the RT ABC structure is much less stable than RT DEF or LT. Low RDS spikes at large negative sign(λ₂)ρ(r) indicate stabilizing intermolecular N(3)H⋯O(2) and C(8)H⋯O(6) hydrogen bonds. Numerous narrow spikes located in the central part of the plot (i.e., near-zero sign(λ₂)ρ(r)) depict dispersive C–H⋯H–C and π⋯π stacking interactions. Their positions for the RT ABC structure differ significantly from those for LT and RT DEF. The RDS spikes at about −0.4 a.u, −0.2 a.u., and −0.15 a.u., as in Figure 16 (top), describe intramolecular N(1)CH⋯O(6), N(1)CH⋯O(2), and N(7)CH⋯O(2) hydrogen bonds. Two spikes at large positive sign(λ₂)ρ(r) values for sterically crowded rings indicate the lack of bonding in their central area.

The spike positions, assigned to nitrogen interactions, as in Figure 17, vary from one structure to another, but that for RT ABC differs significantly from those for LT and RT DEF.

The most stable peak positions are those describing the N(3)H⋯O hydrogen bond, while the most shifted ones are those describing non-covalent interactions in which N(7) and N(9) are involved. This confirms the greater sensitivity of the five-membered ring to the methyl rotation. Moreover, the positions of the peaks describing non-covalent interactions in which N(1) is involved are relatively stable, which confirms the conclusion about the completely different nature of the two methyl groups in PX. The positions of the spikes in the plot of RDG(r) versus sign(λ₂)ρ(r) for N(3) and N(9), both of which have a direct involvement in hydrogen bonds, are quite well correlated with e²qQ/h.

3.4. Binding Mode of Paraxanthine with A2A Receptor

PX has a similar structure to adenosine and, like caffeine, acts as a competitive antagonist to adenosine receptors (A1 and A2A). Unfortunately, the structure of the PX complex with A2A has not yet been experimentally investigated. Therefore, to determine which conformation of the methyl group at N(7) is preferred in a protein–ligand complex, we used molecular docking. The crystal structure of stabilized A2A adenosine receptor A2AR-StaR2-bRIL in complex with caffeine, 5MZP [81], was retrieved from the PDB database. The receptor was prepared by removing caffeine and water molecules. However, two water molecules binding the actual caffeine ligand have been retained and treated as part of the receptor environment. The potential binding site (cavity) of 124.416 ų in volume was detected. This cavity is formed by the following receptor residues: Phe163, Leu249, Ile274, Asn253, Met177, Ile66, Val84, and Glu169 (all except polar Asn253 are hydrophobic). Its shape provides the largest contact area between the ligand and receptor surfaces. The search space of the simulation exploited in the docking studies was defined as a subset region of 9.0–15.0 Å around the active site cleft. To evaluate the docking process, we performed a redocking evaluation task. The ligand was removed from the crystal structure and redocked in the initial ligand-binding pocket of the receptor. The redocking protocol was successful, as the root-mean-square deviation (RMSD) was only 0.11 Å. PX was then docked to a rigid protein structure in exactly the same way. The best pose of conformer that leads to the stabilization of the complex with the highest binding/docking score was selected after docking (the approach for docking the ligands in the protein pocket was the same as reported previously [64]).

The stability of the binding of a complex resulting from molecular docking is often verified by means of molecular dynamics simulations. In this case, random jumps of functional groups are observed, making the study of the binding modes of individual conformers problematic. Therefore, we prepared and docked twelve conformers of the methyl group at N(7) differing in the θ angle 0° ≤ θ ≤ 120° (rotation about its threefold symmetry axis). The strength of the binding between the protein and ligand, measured by the total binding energy and docking score, are strongly influenced by θ. However, both show global minima at 60°, as in Figure 18.

Regardless of the conformer, oxygen O(6) formed two O–H⋯O hydrogen bonds with nearby water molecules, and oxygen O(2) of PX formed one N–H⋯O hydrogen bond with the Asn253 residue. However, the length and strength of N–H⋯O bond, which seems to be a key factor in the A2A receptor–PX interaction, depends on the θ, as in Figure 19.

The best score conformers (global minimum at θ = 60°), as depicted in Figure 20, participate in the strongest hydrogen bond (of 2.902 Å), which is only marginally weaker than that in the actual caffeine ligand (of 2.894 Å). The best score conformer of the methyl group at N(7) in PX is very close to that of the actual ligand (caffeine) in the 5MZP (angle deviation does not exceed 1°). Moreover, the PX moiety in the protein–ligand complex has adopted the same methyl group conformation at N(7) as in the solid state (i.e., the RT DEF conformer).

Detailed insight into the binding mode with LigPlot+ [74,75] confirms the above observations and reveals a set of hydrophobic interactions binding PX with Phe163, Ile66, Ile274, Leu249, Met177, and Val84, as in Figure 21, which further stabilizes the protein–ligand complex.

It should be noted that LigPlot+ suggests a similar set of hydrophobic interactions for the caffeine binding with the same protein. The hydrophobic interaction energies obtained from molecular docking are listed in

Table 11. The energies of the hydrophobic protein–ligand interactions.

Residue Id	EABC [kcal/mol]	EDEF [kcal/mol]	Ecaffeine [kcal/mol]
Phe163	−29.28	−29.35	−30.95
Leu249	−9.68	−9.70	−10.56
Ile274	−8.27	−8.28	−8.25
As253	−5.32	−5.33	−6.49
Met177	−2.66	−2.66	−2.85
Val84	−2.05	−2.06	−1.69
Met270	−1.11	−1.12	−1.40
Ile66	−1.09	−1.08	−1.16
Ala63	−0.98	−0.98	−1.02
Glu169	−0.96	−0.96	−2.40

. All of them, apart from the binding between the ligand and Val84, are stronger when the ligand is caffeine.

The sign[λ₂(r)]ρ(r) surface mapped to the RDG(r) isosurface in the red–green–blue scheme for the protein–ligand complex is shown in Figure 22. This surface visualizes the binding mode and reveals the nature of the non-covalent interactions between the A2A protein and PX ligand. In general, the protein environment makes the electronic density slightly less localized and smooth. Weak van der Waals-type interactions clearly dominate, which is indicated by the green colour of the surface. Most of them are C-H⋯H-C interactions involving the methyl group at the ligand’s N(7) nitrogen and the methyl groups of Ile66, Ile274, and Ala63, and the methyl group at the ligands’s N(1) nitrogen and the methyl groups of Leu249 and Met177 residues. The large and nearly flat green area above the six-membered ring of the PX ring proves π⋯π-stacking between the ligand and Phe163 (average π-stack contact distance: 3.4 Å; center-to-center axis deviation: 2°). The N–H⋯O hydrogen bond is clearly visible, and is depicted by a small light cyan disc-shaped area near O(2).

Conformation changes have a significant effect on the protein–ligand binding, as in Figure 23. The mutual shift in the red and green areas and the red and green spikes, respectively, shows the differences in the binding mode caused by the conformational change from θ = 0° to θ = 60°. In general, the distance from the ligand and the shape of the surface of the binding regions differ among conformers. However, the differences in non-covalent intermolecular interactions are relatively small compared to intramolecular interactions.

Although the methyl group is not directly involved in any strong non-covalent protein–ligand bond, it does contribute to the overall strength of the ligand-to-protein binding. This indirect effect is related to the hydrophobic interactions and the modulation of the electron density throughout the ligand due to the rotation of the methyl group. However, methyl group jumps, if they occur, will cause a subtle effect that will modify but not destroy the key protein–ligand bond. Involvement of the methyl group at N(7) in a strong intramolecular C–H⋯O hydrogen bond supported by dispersive C–H⋯H–C interactions, appears to be a source of disorder as well as an important factor in the formation of protein–ligand complexes.

4. Conclusions

Three factors determine the structural disorder in PX: larger crystal voids near methyl at N(7) than N(1), opening the path for disorder, hyperconjugation that strongly affects the density distribution in the five-membered ring, and the involvement of the methyl group at N(7) in many non-covalent bonds that intercept subsequent jumping protons. By the modulation of the strong hyperconjugation effect, the rotation of the methyl group at the N(7) position affects the electron density in the entire molecule. Therefore, the disorder, that is considered to be local and restricted to the N(7) site, actually affects the entire molecule. As a consequence, methyl rotation significantly modulates the components of the EFG tensor and slightly alters the orientation of its axes for each nitrogen.

Changes in the global 2P FP caused by the rotation of the methyl group at N(7), quantified using the Pompeiu–Hausdorff and Bhattacharayya metrics and the Wasserstein distance, confirm the changes in the distribution and strength of non-covalent interactions occurring throughout the molecule. The influence of methyl group rotation on the pattern of interactions is clearly visible regardless of the type of metric. The magnitude of these percentages is consistent with the changes in the total energy and modulation of the NQR spectra due to the methyl group rotation.

Non-covalent interactions of four types, namely intermolecular hydrogen bonds N−H⋯O, C−H⋯N and C−H⋯O, weak C−H⋯O intramolecular hydrogen bonds (closing 6-member rings and stiffening the structure), numerous H⋯H dispersive interactions, and π⋯π stacking, form the bonding network in PX. The rotation of methyl groups slightly affects the shape of the PX crystal’s energy skeleton, but in-depth analysis reveals differences in the strength of individual interactions and their total distribution. In the LT and RT DEF structures, N(1) participates in two intermolecular hydrogen bonds of different strengths, i.e., C–H⋯O and C–H⋯N, with one dispersive H⋯H contact and one intramolecular hydrogen bond C–H⋯O, while –N(7) participates in one C–H⋯N intermolecular hydrogen bond, two dispersive H⋯H contacts, and a van der Waals N⋯O contact. Nitrogen N(3) acts as a strong proton donor and is involved only in the N–H⋯O bond, while –N(9)= is involved only in the weak C–H⋯N intermolecular hydrogen bond. Thus, only N(3) and N(9) are involved in hydrogen bonds directly, while N(1) and N(7) contribute as the donors of methyl groups, indirectly. Two additional interactions play an important role in the structure of RT ABC: strong dispersive contacts H⋯H and relatively weak intramolecular bonding of the nitrogen methyl group at N(7), C–H⋯O. In the RT ABC structure both methyl groups at N(1) and N(7) seem to be involved in the intramolecular hydrogen bonds C–H⋯O, while in LT and RT DEF only the methyl substituted at N(7) is involved in such a bond. Thus, in RT ABC, the two oxygen atoms strongly bind N(1) with intramolecular hydrogen bonds, resulting in a high increase in the rotation barrier of the methyl group at N(1). Moreover, in the structure of RT ABC, the nitrogen N(7) binds oxygen almost as strongly as N(1) and a few times stronger than in LT or RT DEF. During the rotation of the methyl group, the strength of two intramolecular hydrogen bonds changes significantly. In addition, during rotation, the protons of the methyl group are captured and bound by the CH⋯O bond, which temporarily stabilizes the structure.

The PX moiety in the A2A adenosine receptor–ligand complex has adopted the same methyl group conformation at N(7) as in the solid state (i.e., RT DEF conformer). Thus, the interplay between the intramolecular C−H⋯O hydrogen bonds and C−H⋯H−C dispersive interactions is a key factor in hindering methyl rotation and introducing structural disorder, and an important factor in the formation of protein–ligand complexes. In addition, dispersion interactions play an important role in molecular associations, regardless of the state (solid or protein–ligand complex).