We discuss open heavy flavor ($D$-meson) and quarkonium ($J/\psi$) production in high multiplicity proton-proton ($p+p$) and proton-nucleus ($p+A$) collisions at the LHC in the impact parameter dependent saturation model. A key ingredient in our calculations is the fluctuations of the gluon saturation scale at large-$x$ at each impact parameter, reflecting rare parton configurations inside the projectile and target. Such initial state fluctuations have been used to describe the multiplicity distribution of $D$-meson and $J/\psi$ production yield at mid rapidity in small collision systems at the LHC. In these proceedings, we examine the rapidity dependence of $J/\psi$ production and emphasize its importance in constraining initial state models of high multiplicity events in small systems.