Prediction of self-assembly of adenosine analogues in solution: a computational approach validated by isothermal titration calorimetry

The recent discovery of the role of adenosine-analogues as neuroprotectants and cognitive enhancers has sparked interest in these molecules as new therapeutic drugs. Understanding the behavior of these molecules in solution and predicting their ability to self-assemble will accelerate new discoveries. We propose a computational approach based on density functional theory, a polarizable continuum solvation description of the aqueous environment, and an eﬃcient search procedure to probe the potential energy surface, to determine the structure and thermodynamic stability of molecular clusters of adenosine analogues in solution, using caﬀeine as a model. The method was validated as a tool for the prediction of the impact of small structural variations on self-assembly using paraxanthine. The computational results were supported by isothermal titration calorimetry experiments. The thermodynamic parameters enabled the quantification of the actual percentage of dimer present in solution as a function of concentration. The data suggest that both caﬀeine and paraxanthine are present at concentrations comparable to the ones found in biological samples.


Introduction
Adenosine is a naturally occurring purine nucleoside present in human cells and known to play a key role in a number of pathologies, including neurological disorders, and cardiovascular and autoimmune diseases. [1][2][3] The role of adenosine and its interactions with different receptors have been studied extensively [4][5][6] with consequent interest in adenosine analogues with agonist and antagonist activity 7 and methods that allow fast screening. 8 Among the different adenosine analogues, caffeine (Fig. 1A) has been the most extensively studied compound. Research carried out on the health implications of caffeine intake has established important beneficial links to a number of illnesses, including depression, 9,10 type 2 diabetes, liver disease, Parkinson's disease 11 and stroke risk. 12 Caffeine is rapidly absorbed in the stomach and it is metabolised to three primary compounds: paraxanthine (82%), theophylline (11%) and theobromine (4%). [13][14][15][16] Studies have demonstrated that the key physiological effects of caffeine 17 are the result of interaction effects on the A 1 and A 2 subtypes of the adenosine receptor. 18 The role of caffeine as a cognitive enhancer 19 and also as a neuroprotectant 20 have also been reported recently. Caffeine used in combination with melatonin was shown to have antiamyloidogenic action, for the treatment of neurodegenerative diseases. 21 Borota et al. have demonstrated recently that caffeine enhances consolidation of long-term memory in humans. 22 There is considerable interest in understanding the mechanism by which caffeine and more generally adenosine analogues can interact with different receptors, as this may subsequently lead to the identification of new compounds with drug-like activity. In this context, information on a new molecule's behaviour in solution and its potential for self-assembly, particularly its degree of complexation (i.e. dimers, trimers and tetramers) and the relative distribution, is a key requirement. The self-association of caffeine has been studied experimentally using a variety of techniques including NMR [23][24][25] and FTIR spectroscopies, 26 dynamic light scattering (DLS), 27,28 osmometry, 29 fluorescence, 30 thermal analysis and differential scanning calorimetry (DSC). 31,32 These methods provide some characterisation although the analysis is often model-dependent. Computational approaches have also been explored to provide a useful tool to understand the behaviour of caffeine. Recent molecular dynamics (MD) simulations of aqueous caffeine solutions by Tavagnacco et al. showed extensive aggregation, with caffeine molecules 'stacking their flat faces against one another like coins'. [33][34][35] These simulations were, however, conducted at the limit of solubility of caffeine in water (about 16 mg mL À1 ) which is approximately five orders of magnitude higher than the concentration range at which caffeine is present in biological systems (0.002 mg mL À1 in plasma). [36][37][38] Such low concentrations would require simulation at the limit of MD techniques, containing hundreds of millions of solvent molecules.
The continuum solvation method is an alternative technique for the determination of the structure and thermodynamic stability of molecular aggregates in solution, 39 where the solute (monomer, dimer, trimer, etc.) is treated at the density functional theory (DFT) level and the solvent is described using a polarizable continuum solvation model. Using this approach, Senthilnithy et al. 40 and Bradley and Hendson 41 report highly positive (+12 kcal mol À1 ) dimerization free energies. However, the B3LYP density functional employed in these studies does not allow dispersion interactions to be described fully, 42 which can affect accuracy. The availability of a computational approach that can accurately predict changes in free energies, because of small structural alterations, would accelerate the identification of new candidate drugs. In this work, we present a quantum mechanical continuum solvation approach as a computational tool that allows prediction of the behaviour in solution of adenosine-analogues (AA) by determining the structure and thermodynamic stability of their molecular clusters (AA) n (n = 2-4), in solution, at concentrations comparable to the ones found in biological samples. Caffeine was initially used as a model molecule; the same tool was then employed to predict the effect of small structural variations in the thermodynamics of molecular aggregation by studying paraxanthine (1,7-dymethylxanthine). The computational approach was then validated experimentally by using isothermal titration calorimetry (ITC), to determine the association enthalpies for both molecules. The thermodynamic data were then used to assess the presence and distribution of potential clusters.

Density functional theory methods
Electronic structure calculations were carried out with the NWChem (version 6.3) 43 and Gaussian09 44 codes. The B97 density functional together augmented with the D2 Grimme's dispersion correction (B97-D2) 45 and the 6-31+G(d,p) basis set 46 was chosen to conduct geometry optimization and compute thermal contributions to the gas-phase free energy. The B97-D2 functional was chosen because it is a suitable method to describe molecular complexes with predominant hydrogenbonding and dispersion contribution. [47][48][49][50] The 6-31+G(d,p) basis set was chosen because when evaluating binding energies, the addition of diffuse functions on non-hydrogenic atoms to a doubly polarized valence double-z basis set provides a good compromise between accuracy and computational cost. [51][52][53] The B97-D2/6-31+G(d,p) method, together with the combination of other functionals and basis sets, was also tested in terms in its ability to locate the ''face-to-face'' and ''face-toback'' dimeric motifs found in the crystallographic structure of caffeine (Section S1, ESI †).
The effect of the basis set superposition error (BSSE) was not considered in our simulations. However, the effect of increasing the size of basis set from split-valence double-zeta (6-31+G(d,p)) to split-valence triple-zeta (6-311+G(2df,2p)) on the gas-phase binding energy is only 0.2 kcal mol À1 (Section S2, ESI †). In addition, Di Tommaso and Watson previously demonstrated, using a similar methodology, that the contribution of the BSSE, when evaluated using the counterpoise correction (CPC) method of Boys and Bernardi, 55 to the dimerization constant of carboxylic acids in solutionis less than 0.2 unit. Given that the cost of calculating the CPC is higher than that of the original gas-phase interaction energy, 56 we have not included this term in the evaluation of the gas-phase energy of caffeine and paraxanthine clusters.
The universal solvation model based on solute electron density (SMD) 57 together with the M06-2X density functional 52 and the 6-31+G(d,p) basis set was employed to compute the hydration free energies. 57 Originally tested over a set of 2892 solvation free energies and transfer free energies for both neutral and ionic species in water and non-aqueous solutions, the mean unsigned error obtained with the SMD model over 26 combinations of various basis sets and density functionals was 0.8 kcal mol À1 for neutral solutes. 58 Overall, we expect that the computed Gibbs free energies of a reaction have an uncertainty of AE1.1 kcal mol À1 (Section S2, ESI †). Further calculations to assess the methodology and conduct electron density analysis were also performed using the plane-waves Quantum ESPRESSO 59 software suite with the PBE functional with and without dispersion correction and plane-wave basis set cut offs for the smooth part of the wave functions and the augmented density of 28 and 300 Ry, respectively.

Free energies of association in solution
The free energy of formation of caffeine self-associates, (CAF) n (n = 2-4) was computed according to this equation: where G X * is the total Gibbs free energy of caffeine species X (X = monomer, dimer, trimer or tetramer) in the liquid at 298 K. For example, the dimerization free energy (DG D *) was computed as: The term G* was determined as the sum of the following terms in eqn (3), E e,gas is the gas-phase total electronic energy of the solute, dG VRT;gas is the vibrational-rotational-translational contribution to the gas-phase Gibbs free energy at T = 298 K under a standard-state partial pressure of 1 atm, DG solv * is the solvation free energy of the solute corresponding to transfer from an ideal gas at a concentration of 1 mol L À1 to an ideal solution at a liquid-phase concentration of 1 mol L À1 , and RT ln[RT] is the free energy change of 1 mol of an ideal gas from 1 atm to 1 mol L À1 (RT ln[RT] = 1.89 kcal mol À1 at T = 298.15 K and R = 0.082052 K À1 ).
The terms E e,gas and dG VRT;gas were computed at the B97-D2/ 6-31+G(d,p) level of theory, whereas the hydration free energy component, DG solv *, was determined using the SMD/M06-2X/ 6-31+G(d,p) solvation model using the gas-phase B97-D2/ 6-31+G(d,p) optimised geometries. Although an approach that includes nuclear relaxation in solution of the solute is supposed to improve a method that uses gas-phase geometries, the hydration free energy (DG solv *) of caffeine calculated as the difference between the solute molecule's energy in solution (SMD/M06-2X/6-31+G(d,p)) at the geometry optimized in gas-phase (B97-D2/6-31+G(d,p)) and the gas phase energy (M06-2X/6-31+G(d,p)) at the geometry optimized in the gas phase (DG solv * = À12.8 kcal mol À1 ) is in excellent agreement with the experimental value of À12.1 kcal mol À1 (Table S3, ESI †). The value of DG solv * of caffeine changes only slightly if the hydration free energy is calculated as a difference between the caffeine's energy in solution at the geometry optimized in solution and the gas phase free energy at the geometry optimized in the gas phase (DG solv * = À13.2 kcal mol À1 ) (Table S3, ESI †). This approach also follows the recommendation by Ho et al. that free energies in solution should be obtained from separate gas-and solution-phase calculations, 60 because well-established equations based on the statistical thermodynamics of an ideal gas, under the harmonic oscillator rigid rotor approximation, have been developed to calculate gas-phase partition functions and associated thermodynamic functions (entropy, enthalpy and free energy) from the geometries, frequencies, and energies derived from standard ab initio calculations. If a structure is optimised in solution, then the thermal contribution is obtained by applying the same ideal gas partition functions to the frequencies calculated in a dielectric continuum. This approach is questionable because the translational and rotational contributions are unlikely to be valid in solution, and it should only be applied when stationary points in the solution do not correspond to stationary points in the gas-phase. 61

Localisation of energy minima on the PES of dimers
The potential energy surfaces (PESs) of weakly interacting molecular complexes are characterised by a large number of minima. The most stable isomers of the caffeine dimer, (CAF) 2 , were located using modified version of the computational protocol developed by Di Tommaso and Watson: 54 (i) hundreds of thousands of caffeine dimeric configurations were initially generated using the Granada software. 62 (ii) A first selection of the caffeine dimers was based on geometrical considerations: only those configurations with caffeine molecules at a distance of less than 4.0 Å were considered for subsequent calculations. (iii) The energies of the configurations generated at step (ii) were computed at the B97-D2/6-31+G(d,p) level of theory and used to calculate the Boltzmann factor ( f i ): where E i is the energy of the ith candidate dimeric structure, E 0 is the energy of the most stable generated dimer, the index j runs over all selected configurations, R is the universal gas constant and T is the absolute temperature, which was set at 300 K. (iv) The 110 candidate structures with a Boltzmann factor f i Z 0.01 were subject to full geometry optimization, and the thermochemical properties and solvation energies were calculated. (v) The Boltzmann averaged free energy of association was then calculated according to the following equation: where f i is the Boltzmann factor calculated according to eqn (4) and DG Di * is the association free energy for the ith dimeric configuration. The sum runs over all the selected dimers. The procedure to localise minima on the PES of (CAF) 2 involved 13 000 single point energy calculations (B154k CPU hours using NWChem and 288 cores), followed by 95 structure optimizations (B13k CPU hours using Gaussian09 and 12 cores), frequency calculation (B8k CPU hours using Gaussian and 12 cores) and solvation energy evaluation (B300 CPU hours using Gaussian09 and 12 cores). The total number of CPU hours for the caffeine dimer at the B97-D2/6-31+G(d,p) level of theory was more than 175k CPU hours.

X-ray diffraction analysis
Paraxanthine was recrystallized from methanol by slow evaporation at room temperature. A translucent colourless blocklike crystal with approximate dimensions of 0.080 Â 0.150 Â 0.300 mm was used for X-ray crystallographic analysis. A total of 1088 frames were collected. The total exposure time was 3.02 hours. The frames were integrated with the Bruker SAINT software package using a narrow-frame algorithm. Crystal data were collected at 100 (2)

Dissociation experiments by isothermal titration calorimetry
Isothermal titration calorimetry (ITC) measurements were carried out on a VP-ITC Microcalorimeter (MicroCal, Malvern Instruments) at 25 1C. A solution of caffeine (82.4 mM) or paraxanthine (5.55 mM) was loaded in the syringe and stepwise injected in the cell, in degassed H 2 O. The dilution effect when gradually injecting the concentrated solution of caffeine or paraxanthine in H 2 O made the aggregates to dissociate and the heat generated was measured. The reference cell was also filled with distilled H 2 O. The reference power was set to 4 mcal s À1 ; initial delay 60 s; stirring speed 400 rpm; and equilibration options were set to fast equilibration; feedback mode high. For both experiments, the initial injection was 5 mL and the rest aliquots were of 15 mL with 300 s time spacing between each injection to ensures the signal returns to the baseline. The filter period was 2 s. The duration of the first injection was 10 s and 30 s for the following ones. The total number of injections was 30. The enthalpy (DH diss ) and equilibrium constant (K diss ) of the dissociation in solutions containing either caffeine or paraxanthine were obtained by integrating each heat peak in the raw calorimetric data ( Fig. S7 and S9, ESI †) using a dissociation curve fitting model. The integrated heat plots were acquired ( Fig. S8 and S10, ESI †) and the values of DH diss and K diss were obtained together with their associated errors. The association constant (K ass ), the Gibbs free energy (DG diss ) and the entropy (DH diss ) changes were calculated according to the equation K ass = 1/K diss , DG diss = ÀRT ln K diss , and DG diss = DH diss À TDS diss , where R = 1.987 cal mol À1 K À1 (gas constant) and T = 298 K.

Caffeine dimers
Dimers are the first aggregates that can form during the process of self-association and their structure often corresponds to the building unit of the crystal that grows from solution. [63][64][65][66] The chemical nature of the intermolecular interactions in the caffeine dimer (CAF) 2 can therefore be expected to control its aggregation. The rearrangement of the electronic charge of the face-to-face and face-to-back dimeric configurations found in the X-ray crystalline structure of caffeine 67,68 (Fig. S4b and S4b 0 , ESI †) confirms the presence of six C-HÁ Á ÁO and C-HÁ Á ÁN intermolecular hydrogen bonds between the two caffeine molecules. [69][70][71][72] Although the difference in electronegativity between C and H is small, the weak hydrogen bonds between caffeine molecules play an important role, together with the p-p stacking of the aromatic rings, in controlling the selfassociation process. 73 To evaluate all possible interactions, a computational approach based on probing the potential energy surface of caffeine dimers was applied. Nine low energy configurations for the caffeine dimer (D1-D9) were identified using this approach (Fig. S5, ESI †), five of which had never been reported before. 40,69 All these structures present the CH 3 Á Á ÁX (X = N, O) interactions between the two caffeine molecules, therefore confirming their key role. Overall, the gas-phases energy (DE e,gas ) and Gibbs free energy (DG D 1) of dimerization were found to be negative. However, in aqueous solution, the formation of dimers D1-D8 is thermodynamically unfavourable ( Table 1). The exception is represented by dimer D9 (Fig. 1B), whose liquid-phase reaction is slightly exergonic (DG D * = À0.4 kcal mol À1 ). This therefore suggests that dimer D9 can exist in aqueous solutions. There is much debate on the structural relationship between the most populated dimers in solution and the structural motifs found in solid state forms. 53,64,72 Fig. 1B and 1C highlight how the structure of the most stable caffeine dimer in water (D9), obtained computationally, correlates very well with the structural synthon found in the X-ray crystal structure (Fig. 1C). On this basis, it can be confirmed that the aggregation of caffeine obeys the so-called ''link hypothesis'', 72 according to which the transfer of structural information from the liquid-to the solid state-phase is linked to the presence in the solution of stable ''face-to-back'' dimers.
For dimer D9 to have the lowest free energy in solution, Table 1 suggests that it has a contribution of approximately À2 kcal mol À1 from the hydration free energy. This contrasts with dimers D1-D8, whose solvation corrections are all positive, with a range of 1.6 to 4.5 kcal mol À1 . Assuming a continuum representation of the solvent, the solvation free energy (DG solv *) is given by two components: 57 the electronic polarization term (DG EP *), associated with the electrostatic interaction between the solute and the solvent; the free energy change contribution (G CDS *) from solvent cavitation (C), changes in dispersion (D) energy, and possible changes in the local solvent structure (S). For all dimers, the G CDS * contribution to the solvation energy is between 8.7 and 9.2 kcal mol À1 . Consequently, the extrastabilization in the solution of D9 is due to the electronic polarization term. In this configuration, the monomers are almost superimposed (Fig. 1A) and this configuration does not maximize the attractive interaction between electropositive and electronegative parts of the two caffeine molecules. The gas-phase interaction between the two caffeine molecules (DE e,gas = À12.87 kcal mol À1 ) is lower than in dimers D1-D8 (À17 kcal mol À1 o DE e,gas o À19 kcal mol À1 ) and the extra-stabilization is due to the higher electrostatic interaction between D9 and the polarizable solvation model.
The effect of changes to the solvation environment on caffeine dimerization was also investigated by computing the formation of (CAF) 2 in three different systems: (i) methanol, a polar protic solvent with a static dielectric constant (e) lower than water; (ii) acetone, a polar aprotic solvent, and (iii) chloroform, a non polar aprotic solvent. The Boltzmann averaged Gibbs free energy of dimerization obtained for all three systems (Table S5, ESI †) shows an inverse correlation with the dielectric constants.

Caffeine clusters
Given the existing literature data on higher order clusters for caffeine, the next step focused on the stability of these systems. 28,35 The configuration of the tetramer was obtained by stacking two dimers of type D9 (the most stable in water) on top of each other. The configuration of the trimer was then obtained by removing one caffeine molecule from the (CAF) 4 . The optimised structures of the higher order clusters (CAF) 3 , and (CAF) 4 are reported in Fig. 2.
Two reaction schemes were used to compute the thermodynamic stability of the trimer in solution: (i) the condensation of D9 with an additional caffeine molecule, (CAF) D9 + CAF -(CAF) 3 ; (ii) the aggregation of three isolated CAF molecules, 3CAF -(CAF) 3 . For the tetrameter (CAF 4 ), we considered the condensation of two D9 dimers, 2 CAF D9 -(CAF) 4 , as well as the association of one caffeine molecule to the trimer structure (CAF) 3 + CAF -(CAF) 4 . The energetics of formation of the caffeine trimer and tetramer are listed in Table 2 and show that the free energies of formation of these two species are approximately +2 kcal mol À1 . Trimers and tetramers are unstable in water, compared to both the monomer and dimer D9, and unlikely to be present in caffeine aqueous solutions at the concentrations found in biological samples.

Paraxanthine dimers
Understanding the molecular structure of the analogues and their self-association behaviour is an important requirement for the evaluation of interactions with biological receptors. 74 Paraxanthine (1,7-dimethylxanthine) is the major primary metabolite of caffeine, and although the two molecules are very similar in structure, the removal of one methyl group results in Fig. 2 Optimized structures of trimer (A) and tetramer (B) of caffeine. For the trimer, a top and a side view are reported, which show that the reciprocal orientation of the molecules is retained. In the case of the tetramer, the first two neighbor molecules have an identical orientation while the second two caffeine are rotated, leading to a final geometry of the XX-YY type. The interaction geometry between X and Y is reported. a tenfold decrease in the limit of solubility in water. 75 Compared with the nine low energy dimeric configurations of caffeine, the DFT tool identified 20 stable paraxanthine dimers (Fig. S6, ESI †), showing the significant effect of small structural changes in caffeine analogues on the process of self-association. Two types of structures can be identified: (i) dimers D1, D2, D4 and D13, where the two paraxanthine units are linked by two lateral H-bonds; (ii) all remaining (PX) 2 structures, where the interaction is controlled by p-p stacking. Dimers D4 and D7 correspond to the synthons in the crystallographic structure of paraxanthine (Fig. 3). The agreement between the simulated building units of paraxanthine in solution and its solid-state structure confirms the efficiency of the computational tool here adopted to identify the possible dimeric structures of caffeine analogues. In the gas phase, most dimers are thermodynamically stable but in aqueous solution, only the configurations forming p-p stacking interactions have a negative free energy of formation (À1.0 { DG D * o À0.1 kcal mol À1 ) (Table S8, ESI †). In water, the Boltzmann averaged value for the dimerization of paraxanthine is À0.11 kcal mol À1 . The comparison of the computed thermodynamic stability of the dimers of caffeine and paraxanthine shows that, despite their very small structural difference, the free energy of formation of (PX) 2 is 1.1 kcal mol À1 lower than (CAF) 2 .

Experimental ITC verification and solution speciation
ITC experiments were conducted to obtain the enthalpy and equilibrium constant of dissociation of caffeine and paraxanthine in aqueous solution ( Table 3). The raw calorimetric data are characterised by positive heat peaks ( Fig. S8 and S10, ESI †) indicating that the dissociation for both caffeine and paraxanthine is an endothermic process (DH diss 4 0). Fitting of these data to a dimer dissociation model provides constants for both caffeine (K diss = 0.112 AE 3.3 Â 10 À3 M) and paraxanthine (K diss = 0.03 AE 1.8 m Â 10 À3 M). In both cases, the Gibbs free energy of association (DG ass o 0) is negative. The DG ass value obtained using ITC measurements (À1.3 kcal mol À1 ) is in excellent agreement with a previously reported DSC study (DG ass = À1.5 kcal mol À1 at 25 1C). 76 The free energies of dimerization obtained computationally for caffeine (DG dim = +0.92 kcal mol À1 ) and for paraxanthine (DG dim = À0.11 kcal mol À1 ) are approximately 2 kcal mol À1 higher than the experimental values obtained from the ITC experiments (Table 3). However, if the errors associated with the use of the B97-D2/6-31+G(d,p) method together with the SMD/M06-2X/ 6-31+G(d,p) solvation model [e(DG dim *) =e(DG gas 1) +e(DDG solv *) = 1.1 kcal mol À1 ] are taken into account (Section S2, ESI †) there is a good agreement between computational and experimental results. The small values obtained for the free energy of dimerization of both caffeine and paraxanthine (o2 kcal mol À1 ) also confirm that the monomer and dimer can coexist in solution. The results from the computational approach are considerably closer to the experimental values, and represent a significant improvement compared to previously published data, which reported free energies of dimerization (+12 kcal mol À1 ) considerably higher. 40,41 Furthermore, the difference in the Gibbs free energy of association between caffeine and paraxanthine (DDG ass = DG CAF ass À DG PX ass ) obtained by ITC (0.6 kcal mol À1 ) and DFT (1.1 kcal mol À1 ) is also in excellent accordance, providing strong evidence for the use of our computational approach to predict the self-association behaviour in solution of adenosine-analogues.  Having demonstrated that at low concentrations both caffeine and paraxanthine are potentially present in solution only as a combination of monomers and dimers, the experimental and computational values of DG ass were used in the monomerdimer model to estimate the percentage of dimer in solution as a function of total concentration: 77  Fig. 4 shows the computational and experimental profiles of the dimer percentage in solution for caffeine and paraxanthine in the 5-50 mM concentration range, which is close to the one found in biological samples. The computational data consistently suggest a lower percentage of the dimer in solution compared to ITC for both paraxanthine and caffeine. This is due to the value of DG D * obtained from DFT calculations, for both molecules, being about 2 kcal mol À1 more positive than DG ass obtained from ITC measurements. For both approaches, computational and ITC, paraxanthine is found to be present as a dimer in higher concentrations compared to caffeine, which may be related to the differences in solubility of the two molecules (PX = 5 mM, CAF 82 mM). This behaviour is observed not only at very low concentrations, but also towards the limit of solubility (Fig. S11, ESI †). The results therefore indicate that both paraxanthine and caffeine, when evaluated at the concentrations found in biological samples (10 mM 36 ), are predominantly present as monomers, with less than 0.05% dimer in solution.

Conclusions
We presented a DFT-based computational tool that allowed the study of the self-assembly behaviour in solution of two adenosine-analogues like caffeine and paraxanthine, and we demonstrated that for both molecules, the formation of higher order clusters is thermodynamically disfavoured. The results were validated experimentally by ITC measurements, which provided values for the Gibbs free energy of association consistent with the DFT results. Significantly, we showed how the computational approach allows prediction of the effect of small structural changes on the chemical structure and thermodynamic stability of clusters resulting from self-assembly. In low concentration solutions, both caffeine and paraxanthine are mainly present as monomers, with less than 0.05% dimer present in 10 mM solutions. The computational method used here is a valid tool to quantify the speciation of adenosineanalogues and predict accurately their behaviour in solution, which can provide a standardised and faster approach to the evaluation of potential new drug leads.

Conflicts of interest
There are no conflicts to declare. Fig. 4 Comparison between the computational and experimental profiles of the population of dimers in solution at different total concentrations of caffeine and paraxanthine (mM). The computational points were determined using the Gibbs free energy of dimerization in water of the most stable dimers: caffeine (black circles) and paraxanthine (red circles). The experimental points were calculated using the Gibbs free energy of association determined using ITC experiments: caffeine: (black star) and paraxanthine (red star).