Software & Databases

Based on the network stoichiometry and sign constraints imposed on the boundary fluxes, this Matlab-based package is used to predict thermodynamically feasible reaction directions. With the computed feasible flux directions for a network, the thermodynamic constraint on fluxes may be imposed as a set of linear inequality constraints. In the examples outlined below, this computed set of linear inequality constraints is sufficient to ensure thermodynamic feasibility for certain cases. The methodology is detailed in the following reference.

Components of package:

Step 1. Define network stoichiometry and boundary fluxes. A small reaction networks is illustrated in the following figure. In this network, there are 4 metabolites (labeled A, B, C, and D) and 5 internal reactions (labeled J1 through J5) with boundary fluxes, labeled J6, J7, J8 . The arrows in the figure indicating the direction defined as positive. Boundary fluxes are assumed to have sign pattern sign { J6, J7, J8} = {+, +, +}.

The stoichiometric matrix for this network is:

Step 2. With stoichiometric matrix as an input, cycles.m calculate the minimal cycle basis. The syntax for calling this script is illustrated below:

>> C = cycles(S); % all basis including internal and external cycles.

For the above network, the resulted cycle basis is:

Step 3. Extract the internal cycle basis and impose sign constraints. For example, the first three columns of the above cycle matrix are the basis of internal cycles and they are stored in C_in. The others have to be imposed with boundary constraints and resulted cycle matrix are stored in C_ex. This task was performed by calling extract_basis.m function.

% C_in is the internal cycles, while
% C_ex is cycles with imposed boundary constraints and excluding cycles in C_in.

[C_in, C_ex] = extract_basis(C,ext_flux);

Here is the output from the function extract_basis.m:

Step 4. We compute the feasible signs of the reaction fluxes, given knowledge of the directions of the boundary fluxes, as follows.

Optimization method

The allowable range for the jth flux is computed by finding two optimum values of :

Here, T_constraint.m is to check if a flux violates the thermodynamic constraint in terms of the matroid formalism, and flux.m are objective function. Note that you have to choose different initial points or the optimization process may fail.

Enumeration method

As an alternative to the optimization method described above, it is possible to apply an enumeration method to determine the feasible reaction directions. This method is outlined by the following three steps:

1. Construct a set of feasible sign vectors from the columns of C by deleting the columns that belong to C_in and the columns that violate the imposed boundary sign constraints. This step is already implemented when we call extract_basis.m:

[C_in, C_ex] = extract_basis(C,ext_flux);

2. Now scan the sign vectors obtained from Step 1 row by row. If all entries of a row of matrix C_ex are {+1} and/or {0} , then it is a forward reaction. If entries {+1} and {─} appear in a row, then the sign of this flux is not determined. The predicted flux direction is stored in the vector Sign_Pattern. The function sign_pattern.m performs this step.

Step 5: Check the thermodynamic feasibility. For each internal cycle basis , compare to the predicted flux direction to see whether ‘Sign_Pattern’ possibly contains the same pattern as . If so, the computed linear inequality constraints are not sufficient to ensure thermodynamic feasibility. The codes for this step are as follows.

This output means that the predicted Sign_Pattern contains the first cycle of C_in, which violates the thermodynamic constraints. For this example, the linear constraints {J2, J3, J5} > 0, and sign {J1, J4} {−, 0, +}, may be appended with the non-linear constraint

if

to complete a necessary and sufficient set of thermodynamic constraints.