We extend the numerical analysis of the energy and centrality dependence of particle multiplicities at midrapidity in high-energy p+A and A+A collisions from a running coupling $k_T$-factorization formula made in~cite{Dumitru:2018gjm} by considering two unintegrated gluon distributions that were left out. While a good agreement with the experimental data in A+A collisions is achieved, improving the description of those observables in p+A collisions calls for a better understanding of the proton unintegrated gluon distribution at larger values of $x$ and also the use of a realistic impact parameter dependence.