Using new long-slit spectroscopy obtained with X-Shooter at ESO-VLT, we study, for the first time, radial gradients of optical and near-infrared initial mass function (IMF)-sensitive features in a representative sample of galaxies at the very high mass end of the galaxy population. The sample consists of seven early-type galaxies (ETGs) at z ̃ 0.05, with central velocity dispersion in the range 300 ≲ σ ≲ 350 km s-1. Using state-of-the-art stellar population synthesis models, we fit a number of spectral indices, from different chemical species (including TiO and Na indices), to constrain the IMF slope (i.e. the fraction of low-mass stars), as a function of galactocentric distance, over a radial range out to ̃4 kpc. ETGs in our sample show a significant correlation of IMF slope and surface mass density. The bottom-heavy population (i.e. an excess of low-mass stars in the IMF) is confined to central galaxy regions with surface mass density above ̃ 10^{10} M_☉ kpc^{-2}, or, alternatively, within a characteristic radius of ̃2 kpc. Radial distance, in physical units, and surface mass density are the best correlators to IMF variations, with respect to other dynamical (e.g. velocity dispersion) and stellar population (e.g. metallicity) properties. Our results for the most massive galaxies suggest that there is no single parameter that fully explains variations in the stellar IMF, but IMF radial profiles at z ̃ 0 rather result from the complex formation and mass accretion history of galaxy inner and outer regions.