When it is taken into account that hydrophilic interaction liquid chromatography (HILIC) as an analytical method is relatively young compared with the other techniques, retention modeling could still bring scientifically valuable data to the field. Therefore, in this paper, olanzapine and its 8 impurities were selected as a test mixture, considering that they have never been analyzed in HILIC before. Their investigation on 4 different HILIC columns (bare silica, cyanopropyl, diol and zwitterionic) has been performed. The mixture of 9 structurally similar substances allows the examination of complex HILIC retention behavior depending on the chemical properties of the analytes, as well as of the stationary phase. To describe the nature of the relationship between the retention and the stronger eluent content in the mobile phase, we fitted experimentally obtained data to several theoretical (localized adsorption, nonlocalized partition, quadratic, and mixed) models. Results show that the best fit is the quadratic model with the highest R 2 and cross-validated coefficient of determination (Q 2 ) values, but its usage has some drawbacks. With the aim to improve the possibility to predict retention behavior in HILIC, a new empirical model was proposed. For that purpose, a spline interpolation technique was performed, by dividing the experimental range into several subdivisions. This type of interpolation was performed for the first time in the chromatographic field.The estimation of the polynomial equations was performed using Q 2 values. Obtained Q 2 values pointed out the goodness of fit of the model, as well as its good predictive capabilities. In the end, the prediction capabilities were experimentally verified, under randomly chosen conditions from the experimental range. The errors in prediction were all under 10%, which is satisfying for HILIC.