The goal of this work was to investigate and compare the selectivity of three different hydrophilic interaction liquid chromatography (HILIC) charge modulated amide columns, iHILIC®-Fusion, iHILIC®-Fusion(+) and iHILIC®-Fusion(P), for analysis of compounds in hallucinogen mushrooms. An extract of a truffle-like fungus containing psilocin, psilocybin and baeocystin was chosen as test material. Three different modeling methods were applied to describe the retention times of constituents in isocratic separation mode as a function of mobile phase composition, pH and temperature. Two models using DryLab® 2010 assumed quadratic and exponential relationship between the retention time and solvent fraction of aqueous component of the mobile phase, respectively. These models also illustrate the van't Hoff like equation to describe the temperature-dependence of the retention factor and the theory of Snyder et al. to estimate the retention factor as a function of pH of the aqueous mobile phase component. The third model using STATISTICA® multivariate data analysis in a predefined experimental space was able to predict the retention times. All HILIC columns in this comparison were proved to be suitable for separation of the two hallucinogenic alkaloids from each other and from the matrix components. Majority of compounds were separated with satisfactory resolutions required by the comparative analysis despite some of them were not fully baseline separated. It was found the best modeling was obtained by using the quadratic approach to predict chromatograms for predefined chromatographic conditions (volumetric ratio of acetonitrile to buffer, pH of the buffer and temperature), while the exponential model proved to be the worst for prediction. The modeling with multivariate data analysis fell between the other two methods.