Quantitative characterization of intracellular metabolite concentrations is central to predictive models of cellular functions. Current constraint-based models of metabolism and macromolecular expression provide genotype-phenotype predictions without accounting for intracellular concentrations using simplifying assumptions and optimality principles. However, optimality principles do not generally apply under stress conditions, rendering these predictions unreliable. Incorporation of intracellular concentrations into constraint-based or kinetic models of cellular functions introduces nonlinearities that are computationally challenging to handle, limiting their applicability to small-scale biological networks. Here, we introduce tractable computational techniques to characterize intracellular metabolite concentrations within a constraint-based modeling framework. This model provides a feasible concentration set, which can generally be nonconvex and disconnected. We examine three approaches based on polynomial optimization, random sampling, and global optimization. We leverage the sparsity and algebraic structure of the underlying biophysical models to enhance the computational efficiency of these techniques. We then compare their performance in two case studies, showing that the global-optimization formulation exhibits more desirable scaling properties than the random-sampling and polynomial-optimization formulation, and, thus, is a promising candidate for handling large-scale metabolic networks.