We investigate a spatially explicit metapopulation model consisting of one predator and two hierarchically competing prey species on a discrete lattice. Each local population follows stochastic rules for extinction, colonisation, competition, and predation. From the master equation of this individual-based model, we rigorously derive the corresponding mean-field equations. The analysis of these first-principles mean-field equations reveals the existence of a rich phase diagram with different coexistence regimes depending on parameter values. We identify both stable and spiral nodes, in agreement with the damped oscillations observed in Monte Carlo simulations in the latter case. We find qualitative agreement between the mean-field results and Monte Carlo simulations for the three-species coexistence and for the coexistence of the two prey species. However, the mean-field description fails to reproduce the coexistence between the predator and the inferior prey at high predator coverages. We argue that this discrepancy arises from spatial prey aggregation, which the mean-field approach cannot capture since it neglects correlations. In the stochastic model, spatial clustering acts as a crucial protective mechanism against predation, particularly for the best coloniser. Our findings suggest that prey aggregation contributes to system stability when colonisation and predation operate at comparable spatial scales. The combination of first-principles mean-field equations and stochastic simulations constitutes a powerful framework for clarifying the roles of hierarchical interactions, predation, colonisation, spatial organisation and stochasticity in multi-species coexistence.
