This is the first paper of a series on the methods and results of the Active Galactic Nuclei In Cosmological Simulations ( AGNICS ) project , which incorporates the physics of AGN into GalICS , a galaxy formation model that combines large cosmological N-body simulations of dark matter hierarchical clustering and a semi-analytic approach to the physics of the baryons . The project explores the quasar-galaxy link in a cosmological perspective , in response to growing observational evidence for a close relation between supermassive black holes ( SMBHs ) and spheroids . The key problems are the quasar fuelling mechanism , the origin of the BH to bulge mass relation , the causal and chronological link between BH growth and galaxy formation , the properties of quasar hosts and the role of AGN feedback in galaxy formation . This first paper has two goals . The first is to describe the general structure and assumptions that provide the framework for the AGNICS series . The second is to apply AGNICS to studying the joint formation of SMBHs and spheroids in galaxy mergers . We investigate under what conditions this scenario can reproduce the local distribution of SMBHs in nearby galaxies and the evolution of the quasar population . AGNICS contains two star formation modes : a quiescent mode in discs and a starburst mode in proto-spheroids , the latter triggered by mergers and disc instabilities . Here we assume that BH growth is linked to the starburst mode . The simplest version of this scenario , in which the black hole accretion rate \dot { M } _ { \bullet } and the star formation rate in the starburst component \dot { M } _ { *burst } are simply related by a constant of proportionality , does not to reproduce the cosmic evolution of the quasar population . A model in which \dot { M } _ { \bullet } \propto \rho _ { burst } ^ { \zeta } \dot { M } _ { *burst } , where \rho _ { burst } is the density of the gas in the starburst and \zeta \simeq 0.5 , can explain the evolution of the quasar luminosity function in B-band and X-rays ( taking into account the presence of obscured AGN inferred from X-ray studies ) . The scatter and the tilt that this model introduces in the BH-to-bulge mass relation are within the observational constraints . The model predicts that the quasar contribution increases with the total bolometric luminosity and that , for a given bulge mass , the most massive black holes are in the bulges with the oldest stars .