A quantitative structure–property relation study was performed on the solubility of C60 and C70 fullerene derivatives. Topological and geometrical as well as quantum mechanical energy-related and charge distribution-related descriptors, generated from CODESSA, were calculated to define the molecule structures requirement for measuring their correlations with solubility. The best four variables among the other subsets were selected by the genetic algorithm variable subset selection procedure. Modeling of the relationship between selected molecular descriptors and solubility data was achieved by multiple linear regression method (R2train=0.801, R2test=0.792, Q2LOO=0.716, Q2BOOT=0.674). The robustness and the predictive performance of the proposed model were verified using both internal (cross-validation by leave one out, bootstrap, Y-scrambling) and external statistical validations (external test set by splitting the original data set into training and test sets by k-nearest neighbor (kNN) classification method). Further, the external predictive power of the developed model was examined by considering modified r2 and concordance correlation coefficient values. The reactivity, the polar interactions, the electron–electron repulsion energy, the electronuclear attraction energy, the nuclear–nuclear repulsion energy, and the rotational–vibrational energies were the main independent factors contributing to the solubility of the fullerenes.