The rate of long-duration gamma ray bursts (GRBs) has been identified as a potential proxy for the star formation rate (SFR) across redshift, but the exact relationship depends on GRB progenitor models (single versus binary). The single-progenitor collapsar model accounts for the preference towards low-metallicity GRB progenitors, but is in apparent tension with some high-metallicity GRB host galaxy measurements. As a possible solution, we consider the scenario where high-metallicity GRB hosts harbour low metallicity regions in which GRB progenitors form. For this, we use the IllustrisTNG cosmological hydrodynamical simulation to investigate the internal metallicity distribution of GRB hosts, implementing in post-processing different GRB formation models. Predictions (GRB rate, host metallicities and stellar masses) are compared to the high-completeness GRB legacy surveys BAT6 and SHOALS and a sample of high-redshift GRB-DLA metallicities, allowing us to compute their relative likelihoods. When the internal metallicity distribution of galaxies is ignored, the best-fitting model requires a metallicity-independent channel, as previously proposed by Trenti, Perna & Jimenez. However, when the internal metallicity distribution is considered, a basic metallicity bias model with a cutoff at $Z_{max}=0.35Z_\odot$ is the best fitting one. Current data are insufficient to discriminate among more detailed metallicity bias models, such as weak metallicity dependence of massive binaries vs stronger metallicity bias of collapsars. An increased sample of objects, and direct measurements of host stellar masses at redshift $z>2$ would allow to further constrain the origin of long GRBs.