Computationally efficient flux variability analysis

where the matrix S is an m × n stoichiometry matrix with m metabolites and n reactions and c is the vector representing the linear objective function. The decision variables v represent fluxes, with V ⊆ ℝ ⁿ and vectors v _l and v _u specify lower and upper bounds, respectively. The constraints Sv = 0 together with the upper and lower bounds specify the feasible region of the problem.

Flux variability analysis (FVA) [3] is used to find the minimum and maximum flux for reactions in the network while maintaining some state of the network, e.g., supporting 90% of maximal possible biomass production rate.

Applications of FVA for molecular systems biology include, but are not limited to, the exploration of alternative optima of (1) [3], studying flux distributions under suboptimal growth [4], investigating network flexibility and network redundancy [5], optimization of process feed formulation for antibiotic production [6], and optimal strain design procedures as a pre-processing step [7, 8].

where Z₀ = w ^T v₀ is an optimal solution to (1), γ is a parameter, which controls whether the analysis is done w.r.t. suboptimal network states (0 ≤ γ < 1) or to the optimal state (γ = 1). Assuming that all n reactions are of interest, FVA requires the solution of 2n LPs. While FVA is clearly an embarrassingly parallel problem and is therefore ideally suited for computer clusters, this note focuses on how FVA can be run efficiently on a single CPU. A multi-CPU implementation of fastFVA can be done in the same way as for FVA, i.e., by distributing subsets of the n reactions to individual CPUs. It is expected to give almost linear speedup for sufficiently large problems.

A direct implementation of FVA iterates through all the n reactions and solves the two optimization problems in (2) from scratch each time by calling a specialized LP solver, such as GLPK [9] or CPLEX (IBM Inc.) At iteration i, i = 1, 2, ..., n all elements of c are zero except c _i = 1. Since the only difference between each iteration is a change in the objective function, i.e., the feasible region does not change, solving the LPs from scratch is wasteful. Each time a LP is solved, the solver has to spend some effort in finding a feasible solution. Once a feasible solution is found, the solver then proceeds to locate the optimum. The small changes in the objective function suggest that, on average, the optimum for iteration i does not lie far away from the optimum for iteration i + 1. With Simplex-type LP algorithms, this property can be exploited by solving problem (1) from scratch and then solving the subsequent 2n problems of (2) by starting from the last optimum solution each time (warm-starts). It should be noted that the default behavior of some Simplex-type solvers is to use warm-starts when a sequence of LPs is solved within the same application call. However, current implementations of FVA do not make use of this option (c.f. [10]). Furthermore, for increased efficiency, model preprocessing (presolving) should be disabled after solving the initial problem P. Given a value of 0 < γ ≤ 1, fastFVA performs the following procedure

Setup problem (1), denote it by P

Solve P from scratch to get v₀ and Z₀

Add the constraint w ^T v ≥ γZ₀ to P

for i = 1 to n

   Let c _i = 1 and c _j = 0, ∀j ≠ i

   Maximize P, starting from v _{i- 1}

                     to get v _i and Z _i

   maxFlux _i = Z _i

Once all maximization problems have been solved, the minimization problems are solved in the same way, starting from v₀ = v _n .

An important difference between the various LP solvers available is their ability to exploit multiple core CPUs or multi-processor CPUs to increase performance. The GLPK solver, for example, is a single threaded application. When running on a quad-core machine with hyperthreading enabled, the CPU load is only at 12-13%. On multi-core machines, a significant speedup can often be achieved by simply running multiple instances of fastFVA, each working on a different subset of the n reactions.

The fastFVA package runs within the Matlab environment, which will facilitate the use of fastFVA by users less experienced in programming. In addition, many biochemical network models can be imported into Matlab using the Systems Biology Markup Language (SBML) and the COBRA toolbox [10].

The fastFVA code is written in C++ and is compiled as a Matlab EXecutable (MEX) file (additional file 1). Matlab's PARFOR command is used to exploit multi-core machines. Two different solvers are supported, the open-source GLPK package [9], and the industrial strength CPLEX solver from IBM.

With this efficient FVA tool in hand, new questions can be addressed to study the flexibility of biochemical reaction networks in different environmental and genetic conditions. It is now possible to design computational experiments requiring hundreds or even thousands of FVAs.

The fastFVA package is freely available at http://notendur.hi.is/ithiele/software/fastfva.html together with pre-compiled binaries for Linux and Microsoft Windows. The fastFVA code runs under Matlab and relies on third-party solvers to solve linear optimization problems. Two such solvers are supported, the open source GLPK [9] and the industrial strength CPLEX (IBM Inc.) The fastFVA code is written in C++ and is compiled as a Matlab EXecutable function (MEX). It is released under GNU LGPL.