Standard Eulerian-Lagrangian (EL) methods generally employ drag force models that only represent the mean hydrodynamic force acting upon a particle-laden suspension. Consequently, higherorder drag force statistics, arising from neighbor-induced flow perturbations, are not accounted for; with implications on predictions for particle velocity variance and dispersion. We develop a force Langevin (FL) model that treats neighbor-induced drag fluctuations as a stochastic force within an EL framework. The stochastic drag force follows an Ornstein-Uhlenbeck process and requires closure of the integral time scale for the fluctuating hydrodynamic force and the standard deviation in drag. The former is closed using the mean-free time between successive collisions, derived from the kinetic theory of non-uniform gases. For the latter, particle-resolved direct numerical simulation (PR-DNS) of fixed particle assemblies is utilized to develop a correlation. The stochastic EL framework specifies unresolved drag force statistics, leading to the correct evolution and sustainment of particle velocity variance over a wide range of Reynolds numbers and solids volume fractions when compared to PR-DNS of freely-evolving homogeneous suspensions. By contrast, standard EL infers drag statistics from variations in the resolved flow and thus under-predicts the growth and steady particle velocity variance in homogeneous suspensions. Velocity statistics from standard EL approaches are found to depend on the bandwidth of the projection function used for two-way momentum coupling, while results obtained from the stochastic EL approach are insensitive to the projection bandwidth.