A statistical description of multiphase flows is inevitable due to the inherent variability observed in such systems. The theoretical foundation for the Eulerian-Eulerian (EE) statistical representation of two-phase flows, which is the primary focus of the current study, is established using a probability density function formalism. It is shown that this probabilistic formalism leads naturally to the widely-used ensemble-averaged equations in the EE statistical representation. The relationship between the Lagrangian-Eulerian (LE) statistical representation and the EE formalism is clearly established. In particular, it is shown that the EE and the LE representations bear an exact relationship to each other only under restrictive conditions of local homogeneity of the two-phase flow. The correspondence between unclosed terms in the governing equations that are derived in the two statistical representations is presented. This correspondence allows one to transfer information from one representation to the other at the level of the means. A comparison of the two approaches reveals that the information content in the two representations is indeed different. The interchangeability between Lagrangian and Eulerian descriptions of the carrier phase is investigated. This exercise leads to the formulation of a new statistical representation, namely the Lagrangian-Lagrangian (LL) representation. In the LL formalism, it is shown that the only meaningful way to describe the carrier phase in a Lagrangian frame is through "surrogate" fluid particles. Together, the EE, LE and LL statistical representations presented in this study form a complete framework for the consistent single-point description of two-phase flows. Extension of the EE and LE statistical representations to systems with three and more co-existing phases is outlined. A clearly established theoretical foundation is indeed necessary; also essential is a concomitant improvement in the capability to model unclosed terms in the governing equations of a two-phase flow. Particle dispersion and modulation of turbulent kinetic energy (TKE) by the dispersing particles are two important coupled phenomena that are observed in two-phase flows. Direct numerical simulations (DNS) of canonical homogeneous two-phase flows reveal that the timescales that govern these phenomena behave differently with Stokes number, which is an important non-dimensional quantity that characterizes the relative ease with which the dispersed phase responds to the disturbances in the carrier phase. A new dual-timescale Langevin model (DLM), which is essentially an LL model, is proposed. This model has the unique feature of simultaneously capturing the disparate timescale trends of particle dispersion and interphase TKE transfer with Stokes number. An important ingredient of DLM is a multiscale interaction timescale which is proposed to capture the multiscale nature of particle-turbulence interaction. The behavior of DLM in three canonical homogeneous particle-laden flows, namely freely-decaying turbulence, artificially-forced stationary turbulence and homogeneous shear, is investigated. The versatility of DLM is illustrated by its ability to capture the trends of important statistics that are observed in DNS of the aforementioned canonical two-phase flows with varying Stokes number and mass loading, which is another important non-dimensional quantity that characterizes the relative mass of each phase in a two-phase system. DLM can be extended to inhomogeneous flows with the help of the sound theoretical foundation for multiphase flows that has been established in this work.