Cell population balances constitute a broad class of mathematical models which can accurately capture the dynamics of heterogeneous cell populations. They are first-order partial-integral or partial-functional differential equations. Their formulation is fully defined by the single-cell reaction and division rates and the partition probability density function describing the mechanism by which cellular material is distributed amongst the daughter cells at cell division. Due to the complexity of the formulation, analytical solutions are hard to obtain in most realistic cases. Thus, numerical algorithms are necessary for their accurate approximation. Despite recent progress, the efficient simulation of cell population balance models still remains a challenging task. One of the main challenges comes from the fact that the boundaries of the intracellular state space are typically not known a priori. 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 traditional fixed-boundary algorithms. We will demonstrate the efficiency of this algorithm by using the lac operon gene regulatory network as our model system. The relative advantages of the algorithm will be illustrated through comparison with previously developed algorithms as well as the commercial software COMSOL Multiphysics (formerly FEMLab).