A number of agent-based models (ABMs) for biological wastewater treatment processes have been developed, but their skill in predicting heterogeneity of intracellular storage states has not been tested against observations due to the lack of analytical methods for measuring single-cell intracellular properties. Further, several mechanisms can produce and maintain heterogeneity (e.g., different histories, uneven division) and their relative importance has not been explored. This article presents an ABM for the enhanced biological phosphorus removal (EBPR) treatment process that resolves heterogeneity in three intracellular polymer storage compounds (i.e., polyphosphate, polyhydroxybutyrate, and glycogen) in three functional microbial populations (i.e., polyphosphate-accumulating, glycogen-accumulating, and ordinary heterotrophic organisms). Model predicted distributions were compared to those based on single-cell estimates obtained using a Raman microscopy method for a laboratory-scale sequencing batch reactor (SBR) system. The model can reproduce many features of the observed heterogeneity. Two methods for introducing heterogeneity were evaluated. First, biological variability in individual cell behavior was simulated by randomizing model parameters (e.g., maximum acetate uptake rate) at division. This method produced the best fit to the data. An optimization algorithm was used to determine the best variability (i.e., coefficient of variance) for each parameter, which suggests large variability in acetate uptake. Second, biological variability in individual cell states was simulated by randomizing state variables (e.g., internal nutrient) at division, which was not able to maintain heterogeneity because the memory in the internal states is too short. These results demonstrate the ability of ABM to predict heterogeneity and provide insights into the factors that contribute to it. Comparison of the ABM with an equivalent population-level model illustrates the effect of accounting for the heterogeneity in models.