A multi-scale modeling approach is formulated which combines the use of the discrete element method (DEM) and population balance model (PBM) to simulate the evolution of particle size distribution (PSD) in dry milling. As a major novelty, a pseudo-coupled DEM-non-linear PBM approach is proposed, in which the parameters of a non-linear PBM were calibrated using the DEM simulations of the microdynamic environment in a mill by using the DEM intermittently. To demonstrate the application of this approach, breakage dynamics of silica particles in a vibrating cylindrical vessel containing an alumina grinding ball was simulated numerically. DEM simulations were performed to obtain the collision frequency and impact energy spectra, calculate the first-order breakage rate parameters from short milling of mono-sized feeds, and calibrate the parameters of the effectiveness factor of the non-linear PBM. As another novelty, the form of the effectiveness factor chosen accounted for both acceleration effects and deceleration effects simultaneously, which have been observed in various experimental milling studies. By incorporating DEM input, the non-linear PBM with the effectiveness factor was used to simulate the evolution of the PSD. This study demonstrates that a non-linear PBM, whose parameters were calibrated using DEM at early milling times, could predict the PSD evolution during prolonged milling and potentially obviate the need for using computationally expensive DEM simulations for the whole milling duration. (C) 2017 Elsevier B.V. All rights reserved.