A wide range of cellular processes are inherently stochastic. While stochasticity of gene transcription and translation or cellular growth profiles is well-understood, little is known about the stochastic properties of metabolism. Recent experimental findings strongly suggest that metabolism may indeed be subject to stochastic phenomena, which has questioned the traditional deterministic view of metabolism and casts crucial doubt on the general validity of this modelling paradigm. In this thesis, we examine stochastic aspects of metabolic reactions in detail. We focus on stochastic versions of classic deterministic models for metabolic reactions coupled with well-established stochastic models for gene expression. We incorporate experimental measurements of kinetic parameters in the study, which results in a specific multiscale structure of the presented class of models. In the course of this thesis, we present numerous models with increasing complexity, focussing on three key contributions. Firstly, we present the derivation of an analytical tool to approximate stationary metabolite distributions in closed-form by exploiting the multiple scales. As a result, we propose a strikingly-accurate analytical tool for exploring the parameter space. Secondly, we reveal which parameters have strong impacts on the stationary metabolite distributions and identify conditions for increased coefficients of variation and highly-complex bimodal and multimodal patterns. Finally, we propose a general strategy to obtain closed-form approximations in more complex models, such as multi-step pathways and regulatory processes commonly found in metabolism, such as allostery or end-product inhibition. The results in this thesis lay the groundwork for future studies of metabolic heterogeneity and offer numerous biological hypotheses that could soon be tested in light of recent progress in single-cell measurements of cellular metabolites.