The critical size for stable single-domain (SD) behavior has been calculated as a function of grain elongation for magnetite grains using a numerical micromagnetic algorithm. Importantly, for the first time, we consider the contribution of intergrain magnetostatic interactions on the SD/multidomain (MD) critical size (d0). For individual grains our numerical estimates for d0 for elongated grains are lower than that determined by previous analytical and numerical calculations. Nevertheless, the inclusion of magnetostatic interactions into the model was found to increase d0 to values significantly higher than any previously published estimates of d0 for individual grains. Therefore the model calculations show that there is a relatively wide range of grain sizes within which depending on the degree of magnetostatic interactions and elongation, a grain can be either SD or MD. The model results are compared to observations of magnetosomes found in magnetotactic bacteria. The newly calculated upper d0 limit for the interacting grains now accommodates the largest magnetosomes reported in the literature. These large magnetosomes were previously thought to be MD, suggesting that evolutionary processes are highly efficient at optimizing magnetosome grain size and spatial distribution.