  Quite generally, constraint-based metabolic flux analysis describes the space
of viable flux configurations for a metabolic network as a high-dimensional
polytope defined by the linear constraints that enforce the balancing of
production and consumption fluxes for each chemical species in the system. In
some cases, the complexity of the solution space can be reduced by performing
an additional optimization, while in other cases, knowing the range of
variability of fluxes over the polytope provides a sufficient characterization
of the allowed configurations. There are cases, however, in which the thorough
information encoded in the individual distributions of viable fluxes over the
polytope is required. Obtaining such distributions is known to be a highly
challenging computational task when the dimensionality of the polytope is
sufficiently large, and the problem of developing cost-effective {\it ad hoc}
algorithms has recently seen a major surge of interest. Here, we propose a
method that allows us to perform the required computation heuristically in a
time scaling {\it linearly} with the number of reactions in the network,
overcoming some limitations of similar techniques employed in recent years. As
a case study, we apply it to the analysis of the human red blood cell metabolic
network, whose solution space can be sampled by different exact techniques,
like Hit-and-Run Monte Carlo (scaling roughly like the third power of the
system size). Remarkably accurate estimates for the true distributions of
viable reaction fluxes are obtained, suggesting that, although further
improvements are desirable, our method enhances our ability to analyze the
space of allowed configurations for large biochemical reaction networks.
