Indirect pharmacodynamic response (IDR) models were developed for agents which alter the generation of cell populations with arbitrary lifespan distributions. These models extend lifespan based IDR models introduced previously [J. Pharmacokinet. Biopharm. 27: 467, 1999] for cell populations with the same lifespan ("delta" distribution). Considered are cell populations exhibiting time-invariant lifespan distributions described by the probability density function l(tau). It is assumed that cell response (R) is produced at a zero-order rate (kin(t)) and is eliminated from the population when the cell lifespan expires. The cell loss rate is calculated as kin*l(t), where '*' denotes the convolution operator. Therapeutic agents can stimulate or inhibit production rates according to the Hill function: 1 +/- H(C(t)) where H(C(t)) contains the capacity (Smax) and potency (SC50) parameters and C(t) is a pharmacokinetic function. The production rate is kin(t) = kin.[ 1+/-H(C(t))]. The operational model is dR/dt = kin(t)-kin*l(t) with the baseline condition R0 = kin.TR, where TR is the mean lifespan. Single populations as well as populations with precursors were examined by simulation to establish the role of lifespan distribution parameters (mean and standard deviation) in controlling the response vs. time profile. Estimability of parameters was assessed. Numerical techniques of solving differential equations with the convolution integral were proposed. In addition, the models were applied to literature data to describe the stimulatory effects of single doses of recombinant human erythropoietin on reticulocytes in blood. The estimates of Smax and SC50 for these agents were obtained along with means and standard deviations for reticulocyte lifespan distributions. The proposed models can be used to analyze the pharmacodynamics of agents which alter natural cell production yielding parameters describing their efficacy and potency as well as means and standard deviations for cell lifespan distributions.