We present an analysis of data on single inclusive pion production measured by the LHCf collaboration in high-energy proton-proton and proton-nucleus at ultra-forward rapidities, $8.8 \leq\! y\leq \!10.8$. We also analyse forward RHIC data for calibration purposes. Our analysis relies on the use of a Monte Carlo event generator that combines a perturbative description of the elementary scattering process at partonic level based on the hybrid formalism of the Color Glass Condensate with an implementation of hadronization in the framework of the Lund string fragmentation model. This procedure allows us to reach values of the momenta of the produced particles as low as detected experimentally $p_t\sim0.1$ GeV. We achieve a good description of single inclusive spectra of charged particles and neutral pions at RHIC and the LHC respectively, and nuclear modification factors for proton-lead collisions at the LHC. Our results add evidence to the idea that particle production in the domain of very small Bjorken-$x$ is dominated by the saturation effects encoded in the unintegrated gluon distribution of the target. Being forward particle production of key importance in the development of air showers, we stress that this approach allows for a theoretically controlled extrapolation of our results to the scale of ultra-high energy cosmic rays, thus serving as starting point for future works on this topic.