Cosmological simulations describing both the evolution of supermassive black holes and their host galaxies were performed by using the tree PM-SPH code GADGET-2 (Springel 2005). Physical mechanisms affecting the dynamics and the physical conditions of the gas (ionization and cooling processes, local heating by stars, injection of mechanical energy by supernovae, chemical enrichment) were introduced in the present version of the code (Filloux 2009). Black holes in a state of accretion (AGNs) also inject mechanical energy in the surrounding medium, contributing for quenching the star formation activity. In all simulations a ΛCDM cosmology was adopted (h = 0.7, ΩΛ=0.7, Ωm=0.3, Ωb=0.046 and σ8=0.9). Simulations were performed in a volume with a side of 50h−1 Mpc, starting at z = 50 and through the present time (z = 0). For low and intermediate resolution runs, the initial gas mass particles are respectively 5.35× 108M⊙ and 3.09×108M⊙. Black holes (BHs) are represented by collisionless particles and seeds of 100 M⊙ were introduced in density peaks at z = 15, growing either by accretion or coalescence. The accretion rate from the “disk mode” is based on a turbulent viscous thin disk model whereas in the “spherical mode” the rate is given by the Bondi–Hoyle formula. When accreting matter, jets, modeled by conical regions perpendicular to the disk plane, inject kinetic energy into the surrounding medium. Two models were tested: in the first, the injected energy rate is about 10% of the gravitational energy rate released in the accretion process while in the second, the injected energy rate is based on the Blandford & Znajek (1977) mechanism. All simulations give, at z = 0, similar black hole mass function but they overestimate slightly the BH density for masses above ~ 108M⊙. The resulting BH density in this mass range is affected by feedback processes since they control the amount of gas available for accretion. The present simulations are not able to produce very massive BHs (~109M⊙) at z ~ 6. However the evolution of the BH mass density derived from our simulations are in quite good agreement with that derived from the QSO luminosity function. This indicates that our simulations reproduce quite well the average accretion rate history of BHs. Correlations between the BH mass and properties of the host galaxy (velocity dispersion for bulge systems or the stellar mass or the dark halo mass) are also well reproduced. In conclusion, these exploratory simulations reproduce the data at z = 0 quite well. However, the present adopted recipe for the accretion rate in the “disk mode” seems to be inefficient to produce massive BHs as early as z ~ 6. Higher resolution simulations including a new approach for modeling the “disk mode” are presently under way and that particular difficulty is expected to be solved.