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 efficient search procedure to probe the potential energy surface, to determine the structure and thermodynamic stability of molecular clusters of adenosine analogues in solution, using caffeine 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 caffeine and paraxanthine are present at concentrations comparable to the ones found in biological samples.