Engineering of large-scale metabolic pathways to maximize production of a certain substance or to minimize creation of undesired chemicals require explicit information about the underlying enzyme kinetics. Unfortunately, such information is frequently not available, or when available, it involves uncertainty due to biological variability and to differences in experimental conditions and measurement techniques. In this contribution, we present a computational framework that allows us to characterize the kinetic responses of metabolic network in the presence of uncertainty. The proposed methodology uses Monte Carlo sampling to explore the kinetic data space, and it accounts explicitly for mechanistic properties of enzymes and physico-chemical and thermodynamic constraints. In our bottom-up approach, we start by modeling the kinetic responses of the individual reactions, and then we integrate this mechanistic information with the thermodynamically feasible pairs of substrate and product concentrations to build a population of all biochemically and thermodynamically meaningful models of the metabolic network. In the integration procedure, we exploit the fact that the control coefficients computed for a kinetic mechanism are the elasticities of the enclosing reaction of a metabolic network. The employed sampling method benefits from the specific structure of the underlying constraints and its effectiveness allows implementing the proposed methodology for modeling of large-scale metabolic networks. We demonstrate the properties of the proposed framework through a case study for the identification of rate limiting steps in the central carbon pathway in S. cerevisiae.