  The forthcoming p+Pb run at the LHC will provide crucial in formation on the
initial state effects of heavy ion collisions and on the gluon saturation
phenomena. In turn, most of the saturation inspired phenomenology in heavy ion
collisions borrows substantial empiric information from the analysis of e+p
data, where abundant high quality data on the small-x kinematic region is
available. Indeed, the very precise combined HERA data provides a testing
ground in which the relevance of novel QCD regimes, other than the successful
linear DGLAP evolution, in small-x inclusive DIS data can be ascertained. We
present a study of the dependence of the AAMQS fits, based on the running
coupling BK non-linear evolution equations (rcBK), on the fitted dataset. This
allows for the identification of the kinematical region where rcBK accurately
describes the data, and thus for the determination of its applicability
boundary. It also set important constraints to the saturation models used to
model the early stages of heavy ion collisions. Finally we compare the rcBK
results with NNLO DGLAP fits, obtained with the NNPDF methodology with
analogous kinematical cuts. Further, we explore the impact on LHC phenomenology
of applying stringent kinematical cuts to the low-x HERA data in a DGLAP fit.
