Skip to main content

LooplessFluxSampler: an efficient toolbox for sampling the loopless flux solution space of metabolic models



Uniform random sampling of mass-balanced flux solutions offers an unbiased appraisal of the capabilities of metabolic networks. Unfortunately, it is impossible to avoid thermodynamically infeasible loops in flux samples when using convex samplers on large metabolic models. Current strategies for randomly sampling the non-convex loopless flux space display limited efficiency and lack theoretical guarantees.


Here, we present LooplessFluxSampler, an efficient algorithm for exploring the loopless mass-balanced flux solution space of metabolic models, based on an Adaptive Directions Sampling on a Box (ADSB) algorithm. ADSB is rooted in the general Adaptive Direction Sampling (ADS) framework, specifically the Parallel ADS, for which theoretical convergence and irreducibility results are available for sampling from arbitrary distributions. By sampling directions that adapt to the target distribution, ADSB traverses more efficiently the sample space achieving faster mixing than other methods. Importantly, the presented algorithm is guaranteed to target the uniform distribution over convex regions, and it provably converges on the latter distribution over more general (non-convex) regions provided the sample can have full support.


LooplessFluxSampler enables scalable statistical inference of the loopless mass-balanced solution space of large metabolic models. Grounded in a theoretically sound framework, this toolbox provides not only efficient but also reliable results for exploring the properties of the almost surely non-convex loopless flux space. Finally, LooplessFluxSampler includes a Markov Chain diagnostics suite for assessing the quality of the final sample and the performance of the algorithm.

Peer Review reports


Uniform random sampling has become a useful tool for exploring the solution space of metabolic models in the absence of optimality criteria under different conditions [1]. For a given metabolic model, the solution space is described by a set of equality (mass balances) and inequality (thermodynamic and/or kinetic capacity) constraints. Let us denote by \({{\textbf{S}}} \in \mathbb{R}^{m \times n}\) the stoichiometric matrix of the network, lb and ub \(\in \mathbb{R}^{n}\) the minimum and maximum capacity constraints, respectively, and \({{\textbf{v}}} \in \mathbb{R}^{n}\) a flux solution vector, where n denotes the number of reactions and m the number of balanced metabolites. Then, the mass-balanced solution space \(\varvec{\Omega }\) can be defined as \(\varvec{\Omega } = \{{{\textbf{v}}} \ |\ {{\textbf{S}}} \cdot {{\textbf{v}}} = {{\textbf{0}}},\ {\textbf{lb}} \le {{\textbf{v}}} \le {\textbf{ub}}\}\). Importantly, \(\varvec{\Omega }\) defines a convex body (polytope) on \(\mathbb{R}^{n - \text{rank} \left( {{\textbf{S}}}\right) }\) that can be uniformly sampled and readily explored using Monte Carlo methods, e.g., Hit-and-Run (HR) [2]. Although HR theoretically converges to the uniform distribution on \(\varvec{\Omega }\), its implementation in practice is difficult given the highly heterogeneous scales spanned in metabolic models. To overcome this limitation, different approaches have been proposed for improving its performance [3,4,5]. For instance, the Artificial Centering Hit-and-Run (ACHR) displays accelerated mixing by sampling elongated directions enabling longer steps [5]; unfortunately, the sequence of iterates does not describe a Markov chain, and hence, it is not assured to converge to the target distribution. More recently, the Coordinate Hit-and-Run with Rounding (CHRR) [4] algorithm has been proposed to overcome flux ill-conditioning while maintaining convergence guarantees. Briefly, by computing a maximum volume ellipsoid inscribed in the heterogeneous polytope, the original space can be rounded into a unit ball that can be readily sampled using tractable HR algorithms like the Coordinate Hit-and-Run (CHR) [2]. To this date, this implementation remains the most consistent and efficient method for sampling the mass-balanced flux solution space of genome-scale metabolic models [6], enabling statistical inference of the metabolic phenotype of single- [7, 8] as well as multi-cellular organisms [9, 10], and even microbial communities [11].

A more challenging problem consists of producing a uniform mass-balanced “loopless” flux solution sample from a metabolic model. Flux solutions with active closed loops are not only unrealistic and obscure statistical inference [12], but violate a “loop law” that is analogous to Kirchhoff’s second law for electrical circuits [13, 14]. Any steady-state flux distribution \({{\textbf{v}}}\) can be expressed as \({{\textbf{v}}} = {{\textbf{v}}}_{\text{ll}} + \Delta {{\textbf{v}}}\), where \({{\textbf{v}}}_{\text{ll}}\) denotes a loopless flux distribution and \(\Delta {{\textbf{v}}}\) represents an internal flux distribution (i.e., without active exchanges). The latter flux distributions are both mass-balanced and can be written as [15],

$$\begin{aligned} {{\textbf{v}}}_{\text{ll}}&= \sum _{i} \alpha _i \cdot {{\textbf{e}}}_{i}^{(\text{I})} \nonumber \\ \Delta {{\textbf{v}}}&= \sum _{j} \beta _j \cdot {{\textbf{e}}}_{j}^{(\text{III})} \end{aligned}$$

where \({{\textbf{e}}}_{i}^{(\text{I})}\) and \({{\textbf{e}}}_{j}^{(\text{III})}\) represent type I (pathways with a net conversion of substrates to products) and type III (internal cycles without net conversion) extreme pathways, respectively [16]. In Eq. (1) \({{\textbf{v}}}_{\text{ll}}\) and \(\Delta {{\textbf{v}}}\) are convex combinations of the respective flux modes. Notably, the internal flux distribution fulfills \({{\textbf{S}}}_{\text{int}} \cdot \Delta {{\textbf{v}}}_{\text{int}} = {{\textbf{0}}}\), where “int” refers to the set of internal reactions. Then, the loopless mass-balanced flux solution space can be defined as \(\varvec{\Omega }_{\text{loopless}} = \{{{\textbf{v}}} \ |\ {{\textbf{S}}} \cdot {{\textbf{v}}} = {{\textbf{0}}},\ {\textbf{lb}} \le {{\textbf{v}}} \le {\textbf{ub}},\ \Delta {{\textbf{v}}} = {{\textbf{0}}} \text{ and } {{\textbf{v}}} \ne {{\textbf{0}}}\}\). \(\varvec{\Omega }_{\text{loopless}}\) describes almost surely a non-convex space that is substantially harder to sample from. To tackle this challenge, we previously have proposed an approximate algorithm based on ACHR termed loopless Artificial Centering Hit-and-Run on a Box (ll-ACHRB), which yields a (non-uniform) random sample from \(\varvec{\Omega }_{\text{loopless}}\) in reasonable time [17]. In this work, we present a superior algorithm termed Adaptive Direction Sampling on a Box (ADSB) for generating a provably uniform random flux sample on \(\varvec{\Omega }_{\text{loopless}}\) with substantially higher computational performance and stronger theoretical support.



ADS belongs to a family of population-based MCMC methods that involve sampling from an augmented \(k \times \text{dim}(\Theta )\)-dimensional state space over support(\(\Theta\)). At each iteration, ADS maintains a set of k points (vectors) on \(\Theta\) in a current set \({{\textbf{V}}}^{(t)}\) that is used to adaptively construct convenient directions for traversing support(\(\Theta\)) according to the target distribution \(\pi ({{\textbf{v}}})\) [18]. In this case, \(\pi ({{\textbf{v}}})\) is the uniform distribution over \(\varvec{\Omega }_{\text{loopless}}\) (previously \(\Theta\)), which is spanned by k mass-balanced loopless flux vectors contained in \({{\textbf{V}}}^{(t)}\). Particularly in parallel ADS, a new point \({{\textbf{v}}}^{*}\) is proposed by sampling on the line \({\mathcal {L}}\{ \}\) passing through a random point \({{\textbf{v}}}_{c}^{(t)} \in {{\textbf{V}}}^{(t)}\) (current point), and parallel to a direction \({{\textbf{u}}}^{*}\) constructed from two distinct random points (\({{\textbf{v}}}_1^{(t)},{{\textbf{v}}}_2^{(t)}\)) from \({{\textbf{V}}}^{(t)}\) [18] (Fig. 1A). A uniform random step \(\lambda ^{*}\) is then drawn on the cord \({\mathcal {L}}({{\textbf{v}}}_{c}^{(t)} +\lambda {{\textbf{u}}}^{*})\), and \({{\textbf{v}}}_c^{(t)}\) is replaced by \({{\textbf{v}}}^{*}\) in \({{\textbf{V}}}^{(t+1)}\) leaving all other points unchanged [18]. In this way, the parallel ADS algorithm produces a reversible Markov chain \({{\textbf{V}}} = \{{{\textbf{V}}}^{(t)}, t \in \mathbb{Z}^{+}\}\) in the augmented state space according to the uniform distribution on \(\varvec{\Omega }_{\text{loopless}}\) as long as \(\text{dim}({{\textbf{V}}}^{(t)}) = \text{dim}(\varvec{\Omega }_{\text{loopless}})\) [19]. Direct application of ADS to our case is, however, not possible as the actual support \(\varvec{\Omega }_{\text{loopless}}\) is likely non-convex, which renders sampling \(\lambda ^{*}\) hard. To overcome this obstacle, ADSB borrows the “shrinking box” method from slice sampling [20] to efficiently propose new points on \(\varvec{\Omega }\) and retaining those in \(\varvec{\Omega }_{\text{loopless}}\) [21] (Fig. 1B, C). To check whether a flux vector is loopless, the method proposed in [17] for loop detection was employed, which is based on the fact that active closed loops can be topologically detected directly from the sign pattern of a given flux vector [22]. Importantly, as long as \({{\textbf{V}}}^{(t)}\) spans \(\varvec{\Omega }_{\text{loopless}}\), or equivalently, as long as any \({{\textbf{v}}} \in \varvec{\Omega }_{\text{loopless}}\) can be reached from any \({{\textbf{v}}}^{(t)} \in {{\textbf{V}}}^{(t)}\) [23], the asymptotic convergence to the arbitrary target (uniform) distribution is guaranteed [24]. Finally, the computational performance of ADSB is improved by running K non-interacting Markov chains in parallel each with k points in \({{\textbf{V}}}\). Details of the proposed algorithm such as pre-processing, selection of the initial points, number of iterations and the pseudo-code can be found in the Additional file 1.

Fig. 1
figure 1

Illustration of the ADSB algorithm. A Consider a current set V\(^{(t)}\) composed of \(k = 7\) points. At each step, v\(_{c}\) and two other points (v\(_{1}\),v\(_{2}\)) are randomly drawn from V\(^{(t)}\) and used to construct the direction of movement u\(^*\). A uniform random step \(\lambda ^*\) is taken along a cord passing through v\(_{c}\) and parallel to u\(^{*}\). A new point v\(^{*} =\) v\(_{c} + \lambda ^*\) u\(^{*}\) is proposed. If v\(^{*} \notin \varvec{\Omega }_{\text{loopless}}\), then v\(^*\) is rejected and the initial cord is shrunk. B The length of the initial cord is shrunk and a new step is taken. This time the new point v\(^*\) lands on a feasible region and it is therefore accepted. C Once a feasible point is found, v\(^{*}\) is swapped for v\(_{c}\) in the current set V\(^{(t+1)}\) and the algorithm is repeated

Computational implementation

LooplessFluxSampler is implemented in MATLAB (Mathworks, Natick, MA) with an interface to the popular COBRA Toolbox [25], enabling ready application to constraint-based metabolic models. The software runs in Windows, macOS and Linux, on machines with a 64-bit CPU architecture. Vectorized mathematical code has been employed to speed up parallel computations. LooplessFluxSampler makes use of the Parallel Computing Toolbox from MATLAB to speed computations, albeit the latter is not strictly necessary for single-core computations. For sampling large models, this implementation can take advantage of multiple cores with an almost linear scaling when the Parallel Computing Toolbox is available.

LooplessFluxSampler requires the user to provide the metabolic model in COBRA format (.MAT structure), along with the number of samples to be generated. Additional file 1 such as burn-in, thinning, number of cores, sampler (ADSB, ll-ACHRB, and HR), warmup, among others, can be provided if desired. Before sampling, metabolic models are pre-processed to remove blocked reactions using fast-SNP for fast loopless flux optimization [26]. After the sampling, a complete Markov chain diagnostics can be run to assess the quality of the final sample and the performance of the algorithm. Finally, the results presented here were obtained on a 16-CPU 64-GB ram Virtual Machine hosted in the QRIScloud Polaris cell.

Results and discussion

The performance of ADSB is demonstrated in two scenarios assessing different capabilities. First, we evaluated the convergence behavior on the uniform distribution in a core metabolic model of E. coli [27], where flux solutions have full support [28]. Due to the loopless condition imposition, reactions can only carry flux conditional on the fluxes through other reactions. In some networks, all reactions cannot simultaneously carry flux without violating the loopless condition [17], hindering full exploration of \(\varvec{\Omega }_{\text{loopless}}\) by ADSB. In such cases, only a subspace of \(\varvec{\Omega }_{\text{loopless}}\) can be uniformly sampled as some reactions will be necessarily blocked. Thus, we initially assessed the performance of ADSB in metabolic networks where \(\varvec{\Omega }_{\text{loopless}}\) remained conditionally full-dimensional. For this task, four synthetic networks were analyzed using the model with the original bound (capacity) constraints, but forcing a different number of irreversible reactions to be reversible (Additional file 1). As a consequence, the number of potentially active closed loops in the synthetic networks increased, rendering the sampling task more challenging. The resulting models contained 0, 2, 10 and 12 potentially active loops while enabling all their reactions to carry simultaneously flux.

The convergence behavior of ADSB was compared against HR considering that: (1) HR is proven to converge in total variation to a uniform distribution on convex regions [2], and (2) points from HR are in \(\varvec{\Omega }\) and can be later verified to be in \(\varvec{\Omega }_{\text{loopless}}\) by checking if they are loopless. To ensure a fair comparison, both ADSB and HR were initialized using over-dispersed seeds without warmup. Then, if ADSB produces a sample with similar characteristics to HR, there is a strong support for its convergence to the uniform distribution. This was indeed the case (Fig. 2). Overall, the sample means (Fig. 2A) and standard deviations (Fig. 2B) of all the fluxes showed high consistency and no significant statistical difference was detected (adjusted \(p\text{-value}>0.05\), Wilcoxon signed-ranked test, refer to Additional file 1: Table S3 and Figure S1). There was only one instance in model 12 (NADH dehydrogenase reaction, ID: NADH16), where a flux mean differed appreciably between samplers. Subsequent inspection of this case revealed that HR mixed less efficiently in this particular variable, which explained the observed difference (potential scale reduction factor (psrf) \(\approx\) 1.053 vs. 1.00 of ADSB, Additional file 1: Figure S2). Importantly, in terms of sampling performance, ADSB consistently displayed reduced time per effective sample (Fig. 2C) and psrf very close to 1 and with lower variability than HR (Additional file 1: Figure S3). In fact, for models with a higher number of potentially active closed loops, HR displayed higher time per effective sample and hence a lower performance than ADSB (Fig. 2C). Altogether, the above results support the convergence behavior of ADSB and a superior performance for more complex models than the theoretically proven HR.

Fig. 2
figure 2

Convergence comparison of ADSB and HR in E. coli core metabolic models with varying number of potentially active closed loops. Correlation between flux means (A) and standard deviations (B) between both samplers in all models. Synthetic models R0, R2, R10 and R12 contain 0, 2, 10 and 12 potentially active closed loops, respectively. The reactions forced to be reversible to generate these models can be found in Section 6 of the Additional file 1. C Comparison of the convergence behavior of ADSB and HR. Each model was simulated on minimal medium growing aerobically. Both algorithms were run under the same conditions: \(2 \cdot 10^{5}\) final points, 100 steps per point (thinning), and first \(2 \cdot 10^{4}\) points discarded (burn-in)

Next, we evaluated the computational performance and scalability of ADSB by comparing its performance against ll-ACHRB. Again, to ensure a fair comparison between samplers, both ll-ACHRB and ADSB were initialized using over-dispersed seeds determined through optimization without warmup. For this task, large metabolic models spanning different sizes were chosen (Fig. 3). In these models, all reactions could not carry flux simultaneously without violating the loopless condition, and thus, there is no guarantee of convergence to the uniform distribution for the full model. We note that this guarantee could be recovered for a submodel made of only flux-carrying reactions.

Previous results have shown that in such complex models, i.e., with various potentially active loop laws, it is infeasible to generate a loopless flux sample without blocking specific reactions [17]. In these cases, a strategy like the previously employed in the E. coli core model will be extremely inefficient at producing loopless flux samples. An alternative is to formulate and solve an optimization problem where infeasible loops are removed from flux samples after uniform random sampling [29]. However, even if efficient cycle removal formulations are used (see e.g., CycleFreeFlux [30]), the computation time will be substantial as various thousands of flux samples are typically drawn for sampling analyses. More importantly, this random flux sample has no guarantee of targeting the uniform distribution. This background further supports the need of loopless flux samplers for analysing large metabolic models.

ADSB displayed from an approximately 4-fold increased average computational performance in the small E. coli core model, up to an approx. 1000-fold increased performance in the large iMM904 model when compared to ll-ACHRB (Fig. 3A, Additional file 1: Table S4). The computational performance was measured as the running time divided by the average (over the fluxes) effective sample size \(N_{\text{eff}}\) (the number of independent samples the final correlated sample is worth), which is a fair indicator of the sample quality. Overall, ADSB was substantially superior to ll-ACHRB achieving two orders of magnitude increased computational performance in the smallest genome-scale model (see iIT341, Fig. 3A). Moreover, when comparing mixing properties, only ADSB showed a consistent behavior across all the tested models (\(psrf \approx 1\), Fig. 3B). This was not the case for ll-ACHRB in the larger models iYO844 and iMM904 (Fig. 3B). Importantly, ADSB displayed low psrf variability, which points to a consistent and scalable sampling performance in larger models.

Fig. 3
figure 3

Performance comparison of ADSB and ll-ACHRB in different metabolic models. A Sampling performance and B mixing behavior were evaluated in ADSB and ll-ACHRB in large metabolic models. Each model was simulated on aerobic minimal medium (default simulation conditions). Both algorithms were run using the settings: \(2 \cdot 10^{5}\) final points, 100 steps per point (thinning), and first \(2 \cdot 10^{4}\) points discarded (burn-in)


LooplessFluxSampler provides an efficient tool for randomly sampling the loopless flux solution space of large-scale metabolic models, displaying both superior computational performance and convergence behavior than the alternatives. Under conditions where the loopless flux solution space can be described by flux vectors with full support, this tool provably converges on the uniform distribution ensuring consistent results. Finally, the LooplessFluxSampler is fully compatible with the popular COBRA Toolbox, enabling rapid adoption by modelers and practitioners in the field. We expect to deploy this software on other open-source platforms such as Python and Julia in the near future to expand its use by the scientific community.

Availability of data and materials

The LooplessFluxSampler function library along with example codes are freely available via the GitHub repository for this project:

• Project name: LoopessFluxSampler.

• Project home page (including examples of use):

• Operating system: Platform independent.

• Programming language: MATLAB environment.

• Other requirements: COBRA Toolbox v3.0.

• License: MIT.

• Any restrictions to use by non-academics: None.


  1. Herrmann HA, Dyson BC, Vass L, Johnson GN, Schwartz J-M. Flux sampling is a powerful tool to study metabolism under changing environmental conditions. NPJ Syst Biol Appl. 2019;5:32.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Smith RL. Efficient Monte Carlo procedures for generating points uniformly distributed over bounded regions. Oper Res. 1984;32(6):1296–308.

    Article  Google Scholar 

  3. De Martino D, Mori M, Parisi V. Uniform sampling of steady states in metabolic networks: heterogeneous scales and rounding. PLoS ONE. 2015;10(4):0122670.

    Article  Google Scholar 

  4. Haraldsdóttir HS, Cousins B, Thiele I, Fleming RM, Vempala S. Chrr: coordinate hit-and-run with rounding for uniform sampling of constraint-based models. Bioinformatics. 2017;33(11):1741–3.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Kaufman DE, Smith RL. Direction choice for accelerated convergence in hit-and-run sampling. Oper Res. 1998;46(1):84–95.

    Article  Google Scholar 

  6. Fallahi S, Skaug HJ, Alendal G. A comparison of Monte Carlo sampling methods for metabolic network models. PLoS ONE. 2020;15(7):1–24.

    Article  CAS  Google Scholar 

  7. Price ND, Schellenberger J, Palsson BO. Uniform sampling of steady-state flux spaces: means to design experiments and to interpret enzymopathies. Biophys J. 2004;87(4):2172–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Bordel S, Agren R, Nielsen J. Sampling the solution space in genome-scale metabolic networks reveals transcriptional regulation in key enzymes. PLoS Comput Biol. 2013;6(7):e1000859.

    Article  Google Scholar 

  9. Gomes de Oliveira Dal’Molin C, Quek L-E, Saa PA, Nielsen LK. A multi-tissue genome-scale metabolic modeling framework for the analysis of whole plant systems. Front Plant Sci. 2015;6:4.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Gomes de Oliveira Dal’Molin C, Quek L-E, Saa PA, Palfreyman R, Nielsen LK. From reconstruction to c4 metabolic engineering: a case study for overproduction of polyhydroxybutyrate in bioenergy grasses. Plant Sci. 2018;273:50–60.

    Article  CAS  PubMed  Google Scholar 

  11. Hirmas B, Gasaly N, Orellana G, Vega-Sagardía M, Saa P, Gotteland M, Garrido D. Metabolic modeling and bidirectional culturing of two gut microbes reveal cross-feeding interactions and protective effects on intestinal cells. mSystems. 2022.

    Article  PubMed  PubMed Central  Google Scholar 

  12. De Martino D. Scales and multimodal flux distributions in stationary metabolic network models via thermodynamics. Phys Rev E. 2017;95(6): 062419.

    Article  PubMed  Google Scholar 

  13. Beard DA, Liang S-D, Qian H. Energy balance for analysis of complex metabolic networks. Biophys J. 2002;83(1):79–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Noor E, Lewis N, Milo R. A proof for loop-law constraints in stoichiometric metabolic networks. BMC Syst Biol. 2012;6:140.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Price N, Famili I, Beard D, Palsson B. Extreme pathways and Kirchhoff’s second law. Biophys J. 2002;83(5):2879–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Price ND, Famili I, Beard DA, Palsson B. Extreme pathways and Kirchhoff’s second law. Biophys J. 2002;83(5):2879–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Saa PA, Nielsen LK. ll-achrb: a scalable algorithm for sampling the feasible solution space of metabolic networks. Bioinformatics. 2016;32(15):2330–7.

    Article  CAS  PubMed  Google Scholar 

  18. Gilks WR, Roberts GO, George EI. Adaptive direction sampling. J R Stat Soc Ser D (Stat). 1994;43(1):179–89.

    Google Scholar 

  19. Roberts G, Gilks W. Convergence of adaptive direction sampling. J Multivar Anal. 1994;49(2):287–98.

    Article  Google Scholar 

  20. Neal RM. Slice sampling. Ann Stat. 2003;31(3):705–67.

    Article  Google Scholar 

  21. Kiatsupaibul S, Smith RL, Zabinsky ZB. An analysis of a variation of hit-and-run for uniform sampling from general regions. ACM Trans Model Comput Simul. 2011;21(3):1–11.

    Article  Google Scholar 

  22. Nigam R, Liang S. Algorithm for perturbing thermodynamically infeasible metabolic networks. Comput Biol Med. 2007;37(2):126–33.

    Article  CAS  PubMed  Google Scholar 

  23. Gollub MG, Kaltenbach H-M, Stelling J. Probabilistic thermodynamic analysis of metabolic networks. Bioinformatics. 2021;37(18):2938–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Bélisle CJ, Romeijn HE, Smith RL. Hit-and-run algorithms for generating multivariate distributions. Math Oper Res. 1993;18(2):255–66.

    Article  Google Scholar 

  25. Heirendt L, Arreckx S, Pfau T, Mendoza SN, Richelle A, Heinken A, Haraldsdóttir HS, Wachowiak J, Keating SM, Vlasov V, et al. Creation and analysis of biochemical constraint-based models using the cobra toolbox v. 3.0. Nat Protoc. 2019;14(3):639–702.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Saa PA, Nielsen LK. Fast-snp: a fast matrix pre-processing algorithm for efficient loopless flux optimization of metabolic models. Bioinformatics. 2016;32(24):3807–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Orth JD, Fleming RMT, Palsson B. Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guide. EcoSal Plus. 2010;4(1):10.

    Article  CAS  Google Scholar 

  28. Gerstl MP, Müller S, Regensburger G, Zanghellini J. Flux tope analysis: studying the coordination of reaction directions in metabolic networks. Bioinformatics. 2018;35(2):266–73.

    Article  CAS  PubMed Central  Google Scholar 

  29. Schellenberger J, Lewis NE, Palsson B. Elimination of thermodynamically infeasible loops in steady-state metabolic models. Biophys J. 2011;100(3):544–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Desouki AA, Jarre F, Gelius-Dietrich G, Lercher MJ. CycleFreeFlux: efficient removal of thermodynamically infeasible loops from flux distributions. Bioinformatics. 2015;31(13):2159–65.

    Article  CAS  PubMed  Google Scholar 

Download references


This research was conducted using computational resources from the Queensland Cyber Infrastructure Foundation (


PAS acknowledges the financial support from the Open Seed Fund UC-UQ Grant OSF2021_00193 and from the National Center for Artificial Intelligence CENIA FB210017, Basal ANID. CD acknowledges financial support from an Australian Research Council Future Fellowship (FT210100260). LKN acknowledges support from the Novo Nordisk Foundation (NNF20CC0035580, NNF14OC0009473).

Author information

Authors and Affiliations



Conceptualization: PAS, LKN. Methodology: PAS, CD. Software: PAS. Formal analysis: PAS, SZ. Writing-Original Draft: PAS, SZ. Visualization: PAS, SZ. Writing-Review and Editing: PAS, CD, LKN. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Lars K. Nielsen.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Additional information describing algorithmic details and results of the application of the sampler in benchmark cases.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Saa, P.A., Zapararte, S., Drovandi, C.C. et al. LooplessFluxSampler: an efficient toolbox for sampling the loopless flux solution space of metabolic models. BMC Bioinformatics 25, 3 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: