This article reviews the most critical issues in the simulation of turbulent polydisperse gas-liquid systems. Here the discussion is limited to bubbly flows, where the gas appears in the form of separate individual bubbles. First, the governing equations are presented with particular focus on the generalized population balance equation (GPBE). Then, the mesoscale models defining the evolution of the gas-liquid system (e.g., interface forces, mass transfer, coalescence, and breakup) are introduced and critically discussed. Particular attention is devoted to the choice of the drag model to properly simulate dense gas-liquid systems in the presence of microscale turbulence. Finally, the different solution methods, namely, Lagrangian and Eulerian, are presented and discussed. The link between mixture, two- and multi-fluid models, and the GPBE is also analyzed. Eventually, the different methodologies to account for polydispersity, with focus on Lagrangian or direct simulation Monte Carlo methods and Eulerian quadrature-based moment methods, are also presented. A number of practical examples are discussed and the review is concluded by presenting the advantages and disadvantages of the different methods and the corresponding computational costs.