A Simple and Scalable Sampling Algorithm for Exploration of the Feasible Solution Space of Metabolic Networks

A diversity of computational methods, commonly referred to as constrained-based methods, has been developed to study the characteristics and capabilities of Genome-scale Models (GEMs). Among these methods, uniform random sampling is a popular choice for exploring the full solution space and has, for example, been used to evaluate the impact of physiological constraints on achievable network states and to identify correlated reaction sets. For this type of analysis to be valid however, a key – and often overlooked – assumption is thermodynamic feasibility of the samples, i.e., the samples must satisfy the loopless condition. Imposition of this constraint is challenging as it renders the feasible space non-convex and non-convex spaces are difficult to sample. A simple solution is to generate samples for the convex space and subsequently reject those that do not satisfy the loopless condition; this strategy however lowers the sampling efficiency. Here, we develop a loopless Artificially-Centered Hit-And-Run sampling algorithm (ll-ACHR) capable of generating thermodynamically feasible flux samples. Our algorithm is motivated by a recent variation of the Hit-And-Run algorithm for uniform sampling from general regions termed ‘Hit-And-Run on a Box’, which we have adapted to the problem at hand. Random selection of directions is carried out following the conventional ‘Artificially-Centered Hit-And-Run’, but now using only feasible warmup points. Resolution of a simple LP problem during sampling ensures thermodynamic realizability of the sampled flux distributions, greatly improving speed and sampling efficiency. We illustrate the performance of our algorithm in several metabolic models ranging from small-size models to large GEMs. Here, we show that generation of feasible directions is critical for efficient mixing of the sampler, especially for medium- and large-size models (#reactions>100). In these cases, traditional selection of directions results in practically frozen Hit-And-Run dynamics. We further compared the samples generated by this new algorithm to those generated using ACHR without considering loops. For large models, almost every sample generated using conventional ACHR contains loops and it would be difficult to extract a loopless sample; ignoring loops, however, greatly affects the sample statistics. The current algorithm overcomes these problems and is scalable to large GEMs.