We cast the metabolism of interacting cells within a statistical mechanics framework considering both the actual phenotypic capacities of each cell and its interaction with its neighbors. Reaction fluxes will be the components of high-dimensional spin vectors, whose values will be constrained by the stochiometry and the energy requirements of the metabolism. Within this picture, finding the phenotypic states of the population turns out to be equivalent to searching for the equilibrium states of a disordered spin model. We provide a general solution of this problem for arbitrary metabolic networks and interactions. We apply this solution to a simplified model of metabolism and to a complex metabolic network, the central core of Escherichia coli, and demonstrate that the combination of selective pressure and interactions defines a complex phenotypic space. We also present numerical results for cells fixed in a grid. These results reproduce the qualitative picture discussed for the mean-field model. Cells may specialize in producing or consuming metabolites complementing each other, and this is described by an equilibrium phase space with multiple minima, like in a spin-glass model.