The formation of inhomogeneities within fluidized beds, both in terms of the particle configurations and flow structures, have a pronounced effect on the interaction force between the fluid and particles. While recent numerical studies have begun to probe the effects of inhomogeneities on the drag force at the particle scale, the applicability of prior micro-scale constitutive drag relations is still limited to random, homogeneous distributions of particles. Since an accurate model for the drag force is needed to predict the fluidization behavior, the current study utilizes the lessons of prior inhomogeneity studies in order to derive a robust drag relation that is both able to account for the effect of inhomogeneities and applicable as a constitutive closure to larger-scale fluidization simulations. Using fully-resolved lattice Boltzmann simulations of systems composed of fluid and monodisperse spherical particles in the low-Reynolds-number (Re) regime, the fluid-particle drag force, normalized by the ideal Stokes drag force, is found to significantly decrease, over a range of length scales, as the extent of inhomogeneities increases. The extent of inhomogeneities is found to most effectively be quantified through one of two subgrid-scale quantities: the scalar variance of the particle volume fraction or the drift flux, which is the correlation between the particle volume fraction and slip velocity. Scale-similar models are developed to estimate these two sub-grid measures over a wide range of system properties. Two new drag constitutive models are proposed that are not only functions of the particle volume fraction and the Stokes number (St), but also dependent on one of these sub-grid measures for the extent of inhomogeneities. Based on the observed, appreciable effect of inhomogeneities on drag, these new low Re drag models represent a significant advancement over prior constitutive relations.