We employ a binary population synthesis, combined with a detailed treatment of mass transfer onto neutron stars and black holes from the optical component with the help of the MESA stellar evolution code, to model populations of binary X-ray sources in galaxies. Wind-fed accreting neutron stars are calculated with taking into account different regimes of quasi-spherical accretion onto magnetised neutron stars (Bondi-Hoyle or settling subsonic, depending on the plasma cooling near the neutron star magnetosphere). Features of the supercritical disk accretion onto black holes and magnetised neutron stars are also taken into account. Our model enables us to reproduce the observed properties of specific X-ray binary types, including symbiotic X-ray binaries, Be-X-ray binaries and ultraluminous X-ray sources. We also model the X-ray luminosity function from different X-ray binary populations in galaxies with an account of the star formation rate history and metallicity evolution. The X-ray luminosity function can be used as a tool to study properties of X-ray binary populations observed in X-ray sky surveys.