DistributedFBA.jl: high-level, high-performance flux balance analysis in Julia.

Motivation: Flux balance analysis and its variants are widely used methods for predicting steady-state reaction rates in biochemical reaction networks. The exploration of high dimensional networks with such methods is currently hampered by software performance limitations.
Results: DistributedFBA.jl is a high-level, high-performance, open-source implementation of flux balance analysis in Julia. It is tailored to solve multiple flux balance analyses on a subset or all the reactions of large and huge-scale networks, on any number of threads or nodes. Availability and Implementation: The code is freely available on github.com/opencobra/COBRA.jl. The documentation can be found at opencobra.github.io/COBRA.jl. Contact: ronan.mt.fleming@gmail.com. Supplementary information: Supplementary data are available at Bioinformatics online.

1 Introduction

Constraint-based reconstruction and analysis (COBRA) (Palsson ) is a widely used approach for modeling genome-scale biochemical networks and integrative analysis of omics data in a network context. All COBRA predictions are derived from optimization problems, typically formulated in the form where represents the rate of each biochemical reaction, is a lower semi-continuous and convex function, is a stoichiometric matrix for m molecular species and n reactions, and b is a vector of known metabolic exchanges. Additional linear inequalities (expressed as a system of equations with matrix C and vector d) may be used to constrain combinations of reaction rates and keep reactions between upper and lower bounds, u and l, respectively.

In flux balance analysis (FBA), one obtains a steady-state by choosing a coefficient vector and letting and . However, the biologically correct coefficient vector is usually not known, so exploration of the set of steady states relies on the embarrassingly parallel problem of solving (1) for many c. Moreover, while is unique for an optimal flux vector , there may be alternate optimal solutions. In flux variability analysis (FVA), one finds the extremes for each optimal reaction rate by choosing a coefficient vector with one nonzero entry, then minimizing and maximizing , subject to the additional constraint for each reaction in turn ().

For kilo-scale models (), the 2n linear optimization problems required for FVA can currently be solved efficiently using existing methods, e.g. FVA of the COBRA Toolbox, fastFVA, or the COBRApy implementation (Schellenberger ; Gudmundsson ; Ebrahim ). However, these implementations perform best when using only one computing node with a few cores, which becomes a temporal limiting factor when exploring the steady state solution space of larger models. Julia is a high-level, high-performance dynamic programming language for technical computing (Bezanson ). Here, we exploit Julia to distribute sets of FBA problems and compare its performance to existing implementations.

2 Overview and implementation

DistributedFBA.jl, part of a novel COBRA.jl package, is implemented in Julia and makes use of the high-level interface MathProgBase.jl (Lubin ; see Supplementary Material). A key feature is the integrated capability of spawning synchronously any number of processes to local and remote workers. Parallelization is primarily achieved through distribution of FBA problems (outer layer), while parallelization of the solution algorithm is solver based (inner layer). COBRA.jl extends the COBRA Toolbox (Schellenberger ) while existing COBRA models (Orth ) can be input.

3 Benchmark results

DistributedFBA.jl and fastFVA (Gudmundsson ) were benchmarked on a set of models of varying dimension (Table 1). All experiments were run on several DELL R630 computing nodes with 2 × 36 threads and 768GB RAM running Linux. As Julia is a just-in-time language, pre-compilation (warm-up) was done on a small-scale model before benchmarking (Orth ). The creation of a parallel pool of workers and the time to spawn the processes are not considered in the reported times.

Table 1

Sizes of S for benchmark models

#	Model name	Metabolites m	Reactions n	References
1	Recon1	2785	3820	Duarte et al. (2007)
2	Recon2	5063	7440	Thiele et al. (2013)
3	Recon3a	7866	12 566
4	Recon2 + 11M	19 714	28 199	Heinken et al. (2015)
5	Multi-organb	47 123	61 230
6	SRS064645	89 756	99 104	Magnusdottir et al. (2016)
7	SRS011061	126 682	139 420	Magnusdottir et al. (2016)
8	SRS012273	186 662	208 714	Magnusdottir et al. (2016)

Brunk, E. et al. (2016) Recon 3d: a three-dimensional view of human metabolism and disease (in revision).

Thiele, I. et al. (2016) Multi-organ model (prototype model) (in preparation).

The serial performance of both implementations is within 10%. The uninodal performance of fastFVA is slightly higher on a few threads, but the performance of distributedFBA.jl is superior for a higher number of threads on a single node (Fig. 1A). The way the FBA problems are distributed among workers (distribution strategy s, see Supplementary Material) yields an additional speedup of 10–20% on a larger number of threads.

Fig. 1

Performance of distributedFBA for selected benchmark models given in Table 1. (A) Speedup factor relative to fastFVA as a function of threads and distribution strategy s (1 node). (B) Multi-nodal speedup in latency and Amdahl’s law (s = 0)

According to Amdahl’s law, the theoretical speedup factor is , where N is the number of threads and p is the fraction of the code (including the model) that can be parallelized. The fraction p increases with an increasing model size (Fig. 1B). The maximum speedup factor for a very large number of threads N is . All reactions of models 6–8 given in Table 1 have been optimized (with full output, s = 0) on 4 nodes/256 threads in only , and , respectively. This demonstrates that for high-dimensional models, it is critical to have a large number of threads on multiple high-memory nodes to accrue a significant speedup.

4 Discussion

The multi-nodal performance of distributedFBA.jl is unparalleled: the scalability of distributedFBA.jl matches theoretical predictions, and resources are optimally used. Key advantages are that the present implementation is open-source, platform independent, and that no pool size limits, memory, or node/thread limitations exist. Its uninodal performance is similar to the performance of fastFVA on a few threads and about 2–3 times higher on a larger number of threads. A key reason is the direct parallelization capabilities of Julia and the wrapper-free interface to the solver. The unilingual and easy-to-use implementation relies on solvers written in other languages, allows the analysis of large and huge-scale biochemical networks in a timely manner, and lifts the analysis possibilities in the COBRA community to another level.

Click here for additional data file.