  Empirical "backward" galaxy evolution models for infrared bright galaxies are
constrained using multi-band infrared surveys. We developed a new Monte-Carlo
algorithm for this task, implementing luminosity dependent distribution
functions for the galaxies' infrared spectral energy distributions (SEDs) and
for the AGN contribution, allowing for evolution of these quantities. The
adopted SEDs take into account the contributions of both starbursts and AGN to
the infrared emission, for the first time in a coherent treatment rather than
invoking separate AGN and star-forming populations. In the first part of the
paper we consider the quantification of the AGN contribution for local universe
galaxies, as a function of total infrared luminosity. It is made using a large
sample of LIRGs and ULIRGs for which mid-infrared spectra are available in the
Spitzer archive. In the second part we present the model. Our best-fit model
adopts very strong luminosity evolution, $L=L_0(1+z)^{3.4}$, up to $z=2.3$, and
density evolution, $\rho=\rho_0(1+z)^2$, up to $z=1$, for the population of
infrared galaxies. At higher $z$, the evolution rates drop as $(1+z)^{-1}$ and
$(1+z)^{-1.5}$ respectively. To reproduce mid-infrared to submillimeter number
counts and redshift distributions, it is necessary to introduce both an
evolution in the AGN contribution and an evolution in the
luminosity-temperature relation. Our models are in plausible agreement with
current photometry-based estimates of the typical AGN contribution as a
function of mid-infrared flux, and well placed to be compared to upcoming
Spitzer spectroscopic results. As an example of future applications, we use our
best-fitting model to make predictions for surveys with Herschel.
