Deriving tractable reduced equations of biological neural networks capturing the macroscopic dynamics of sub-populations of neurons has been a longstanding problem in computational neuroscience. In this paper, we propose a reduction of large-scale multi-population stochastic networks based on the mean-field theory. We derive, for a wide class of spiking neuron models, a system of differential equations of the type of the usual Wilson-Cowan systems describing the macroscopic activity of populations, under the assumption that synaptic integration is linear with random coefficients. Our reduction involves one unknown function, the effective non-linearity of the network of populations, which can be analytically determined in simple cases, and numerically computed in general. This function depends on the underlying properties of the cells, and in particular the noise level. Appropriate parameters and functions involved in the reduction are given for different models of neurons: McKean, Fitzhugh-Nagumo and Hodgkin-Huxley models. Simulations of the reduced model show a precise agreement with the macroscopic dynamics of the networks for the first two models.