Metabolism plays a key role in a multitude of different biological processes ranging from food production and biofuel production to human health. Predicting the metabolism of a living organism, however, can be a challenging task. Genome-scale models (GEMs) can provide this predictive power by accounting for all metabolic reactions in an organism's genome. So far, GEMs have been used to model metabolism through optimization approaches, but this approach shows limitations. We propose a new approach based on a combination of Markov Chain Monte Carlo and Bayesian inference that provides all metabolic states compatible with the available experimental data. We discuss efficient sampling techniques which can leverage high performance computing to efficiently handle the associated computational burden. These techniques are based on Hamiltonian Monte Carlo methods that leverage artificial neural networks for efficient gradient calculation. The corresponding numerical results for case studies related to predictive modeling of metabolism are presented and analyzed. This technique represents a first step towards modeling microbial communities in the future.