Direct Numerical Simulations (DNS) of multiphase flows have progressed rapidly over the last decade and it is now possible to simulate motions of hundreds of deformable bubbles in turbulent flows. The availability of different statistics calculated from such DNS data could help advance the development of new reduced order models of the average or large-scale flows. DNS simulation of a laminar system with nearly spherical bubbles have been used to examine the full transient motion and shows a non-monotonic evolution where all the bubbles first move toward the walls and then the liquid slows down, eventually allowing some bubbles to return to the center of the channel. DNS simulation of a turbulent system with bubbles of different sizes at a friction Reynolds number of 250 shows that small bubbles quickly migrate to the wall, but the bulk flow takes much longer to adjust to the new bubble distribution. The DNS results are then used to help develop the averaged model equations with unknown closures accounting for the effect of the unresolved scales. A database is generated by averaging the DNS results over planes parallel to the stream-wise direction of the flow. The most important turbulent quantities are selected first with the feature selection algorithm. Then a Neural Network (NN) is used for the a priori test, to examine the relationships between unknown closure terms in averaged models for the flow and quantities that are available through the DNS results. For laminar cases, the closure relations are then tested, by following the evolution of different initial conditions, and it is found that the model predictions are in reasonably good agreement with DNS results. For turbulent cases, the preliminary results for the feature selection and a priori test are promising and robust, the future work is to validate the predictive performance of the turbulent modeling.