In this paper we prove a new combinatorial formula for the modified Macdonald polynomials H λ (X; q, t), motivated by connections to the theory of interacting particle systems from statistical mechanics. The formula involves a new statistic called queue inversions on fillings of tableaux. This statistic is closely related to the multiline queues which were recently used to give a formula for the Macdonald polynomials P λ (X; q, t). In the case q = 1 and X = (1, 1, . . . , 1), that formula had also been shown to compute stationary probabilities for a particle system known as the multispecies ASEP on a ring, and it is natural to ask whether a similar connection exists between the modified Macdonald polynomials and a suitable statistical mechanics model. In a sequel to this work, we demonstrate such a connection, showing that the stationary probabilities of the multispecies totally asymmetric zero-range process (mTAZRP) on a ring can be computed using tableaux formulas with the queue inversion statistic. This connection extends to arbitrary X = (x 1 , . . . , xn); the x i play the role of site-dependent jump rates for the mTAZRP.