CONTRIBUTORS: Written: Regan ODONGO, Edited: Tunahan ÇAKIR

Procedure for using Gurobi solver in MATLAB for FBA

Gurobi solver can used to solve most programming problems. In the linear system we defined above, we seek to solve the following problem:

$$min/max \sum_{i=0}^n f(v_i),$$

such that

$$S×v \leq b,$$

$$Seq×v = beq,$$

$$lb \leq v \leq ub.$$

The output of this optimization problem is the vector v, which optimizes the function f. The general implementation steps are as follows:

$$sol=gurobi(model).$$

Here, the object ‘model’ is a MATLAB structure object with the following main fields:

• model.A → An mxn matrix. the S matrix (stoichiometric matrix) must be renamed as model.A and stored as a sparse matrix. Using sparse(S) returns a sparse matrix for the S matrix.

• model.obj → An nx1 vector. The objective function vector, f, explained above must be renamed as model.obj.

• model.lb → An nx1 vector. The lower bound vector.

• model.ub → An nx1 vector. The upper bound vector.

• model.rhs → A zero m×1 vector. The right-hand-side vector, b, explained above must be renamed as model.rhs.

• model.vtype → An nx1 vector. It is a cell array that can be of ‘B’, ‘C’ or ‘I’ for binary, continuous and integer respectively based on the problem being solved. In FBA, all the unknowns (rates) are continuous. Therefore, all elements of model.vtype is set to ‘C’.

• model.sense → An mx1 vector. It is a cell array that can be of ‘=’, ‘<’ or ‘>’ for equality or inequality. It defines the model equality/inequality constraint. In FBA, all the linear set of equations derived from the differential mass balances are equality equations. Therefore, all elements of model.sense is set to ‘=’.

For simplicity, the above object can be defined as follows in MATLAB:

        model.A = sparse(S);

        model.obj = zeros(1,size(A,1));

        model.lb = lb;

        model.ub = ub;

        model.rhs = zeros(size(A,2),1);

        model.vtype = 'C';

        model.sense = '=';

Finally, running

$$sol=gurobi(model);$$

will generate an output, ‘sol’, which is a structure object with several fields. The most important fields in the structure object are sol.x, sol.status, and sol.objval, which stores the predicted v vector (the flux vector), whether the obtained solution is optimal (sometimes it might be ‘INFEASIBLE’, ‘SUB-OPTIMAL’ or ‘NUMERIC’. Care should be taken especially for ‘SUB-OPTIMAL’ since it also returns a flux vector which is not possible when the status is ‘NUMERIC’ or ‘INFEASIBLE’) and the value of the objective function, respectively.

Under circumstances where no optimal solution could be found by the solver, finding the root cause of the problem should begin by checking out numerical consistencies in the defined system constraints. When this does not solve the problem, getting into more advanced troubleshooting approaches using parameter control methods should follow. More detailed information on Gurobi solver parameter control can be found here. What should guide the choice of a new parameter is taking into consideration the solver’s default parameters and weighing the cost of each change on the accuracy of the results to be obtained.