There exists an abundance of experimental evidence in a variety of systems. showing that cell populations are heterogeneous systems in the sense that properties such as size, shape, DNA and RNA content are unevenly distributed amongst the cells of the population. The quantitative understanding of heterogeneity is of great significance, since neglecting its effect can lead to false predictions. Cell population balance models are used to address the implications of heterogeneity and can accurately capture the dynamics of heterogeneous cell populations. They are first-order partial-integral differential equations and due to the complexity of formulation, analytical solutions are hard to obtain in the majority of cases. Despite the recent progress, the efficient solution of cell population balance models remains a challenging task. One of the main challenges stems from the fact that the boundaries of the intracellular state space are typically not known a priori and using fixed-boundary algorithms leads to inaccuracies and increased computational time demands. Motivated by this challenge, we formulated a free boundary finite element algorithm, capable of solving cell population balance equations more efficiently than the traditional fixed-boundary algorithms. The implementation of the algorithm is accommodated, in the finite element based software package COMSOL Multiphysics. We demonstrate the efficiency of this algorithm using the lac operon gene regulatory network as our model system and perform transient and asymptotic behavior analysis. In the latter case, the pseudo-arc-length continuation algorithm is incorporated, in order to investigate the existence of a bistability region, also observed at the single-cell level. Our analysis, revealed the existence of a region of bistability when cell heterogeneity is taken into account; however, its extend shrinks comparing to homogeneous cell populations. The free boundary algorithm can be easily extended for problems of higher dimensionality and we present results for a two-dimensional cell population balance model, which can exhibit an oscillatory behavior. It is shown that oscillations do not persist when the intracellular content is unevenly distributed amongst the daughter cells. (C) 2009 Elsevier Ltd. All rights reserved.