Background and Objectives: Marketed formulations ior®EPOCIM, MIRCERA® and two new branched PEGylated erythropoietin formulations (32kDa-PEG2-rHuEPO and 40kDa-PEG2-rHuEPO) were administered a single dose by intravenous bolus in New Zealand rabbits. The aim of this work was to develop a population pharmacokinetic/pharmacodynamic (PK/PD) model describing in a simultaneous and integrated form the time course of erythropoiesis stimulation after erythropoietin administration. Methods: The PK/PD dataset consisted of 266 EPO measurements in the pharmacokinetic analysis and 799 observations of reticulocyte (RET), erythrocyte (RBC) and hemoglobin (HGB) in the pharmacodynamic analysis from 19 and 20 rabbits, respectively. The First Order Conditional Estimation Method with INTERACTION as implemented in NONMEM version 7.3 was used in PK/PD modeling. Results: A semi-mechanistic cell transit model that included a two compartment model with linear elimination and cell proliferation, maturation, and homeostatic regulation provided a good description of the data regardless the type of erythropoietin formulation administered. Inter-subject variability was associated with clearance (50%) and apparent volume of distribution of the central compartment (60%). Residual unexplained variability was estimated to be 20.5 %. The system- and drug-related parameters showed consistency and differed across formulations, respectively. Conclusions: A semi-mechanistic PK/PD model of erythropoiesis provided and adequate description of the observed data of ior®EPOCIM, MIRCERA® and PEG-EPO 32 kDa and PEG-EPO 40 kDa in New Zealand rabbits. The administration of new branched PEG-chains formulations improves PK and PD properties, in terms of increasing elimination half-lives and pharmacological response on RET, RBC and HGB compared to commercially available formulations (ior®EPOCIM and MIRCERA®).