Ab Initio Prediction of Thermodynamically Feasible Reaction Directions from Network Stoichiometric Matrix
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 :
= max | = 0, , for each , and
sign { J6, J7, J8} = {+, 0, +}, {+, +, 0}, or {+, +, +} , and
, for i = 1, …, n
= min | = 0, , for each , and
sign { J6, J7, J8} = {+, 0, +}, {+, +, 0}, or {+, +, +} , and
, for i = 1, …, n
This optimization process can be performed by following codes:
LB = zeros(nc_ex,1);
UB = +10.*ones(nc_ex,1); %UB is an array of arbitrary finite positive real numbers
% A feasible initial guess for calculations below:
xo = ones(nc_ex,1);
for i = 1:n_intmax_min= -1;
[x,F,xflag(i,1)] = fmincon(@flux,xo,-C_ex(n_int+1:n_all,:),…
zeros(n_ext,1),[],[],LB,UB,@T_constraint,[],C_in,C_ex,max_min,i);
J_min(i) = F;max_min= 1;
[x,F,xflag(i,2)] = fmincon(@flux,xo,-C_ex(n_int+1:n_all,:),…
zeros(n_ext,1),[],[],LB,UB,@T_constraint,[],C_in,C_ex,max_min,i);
J_max(i) = -F;end
if sum(sum( xflag<0)) ~= 0
disp(‘error in optimization’);
xflagend
[(1:n_int)’ J_min’ J_max’ ]
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.
You should get a result that looks like this:
Index
|
J_min
|
J_max
|
1
|
-10.0000
|
20.0000
|
2
|
0 |
40.0000
|
3
|
0
|
13.9055
|
4
|
-20.0000
|
10.0000
|
5
|
0
|
30.0000
|
The feasible signs for the jth flux are easily deduced from the above results. With the boundary fluxes are assumed to have signs: sign {J6, J7, J8} = {+, +, +}. , the computed inequality constraints are { J2, J3, J5} > 0, and sign { J1, J4} {−, 0, +}.
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.
Sign_Pattern = sign_pattern(S, C_ex, ext_flux);
n = length(Sign_Pattern);
[(1:n)’ Sign_Pattern]
This is the output result from ‘sign_pattern.m’. Here, ‘-11’ denotes reversible reaction; ‘+1’ denote forward direction and ‘-1’ denotes backward direction.
Index
|
Flux direction
|
1
|
-11
|
2
|
1
|
3
|
1
|
4
|
-11
|
5
|
1
|
6
|
1
|
7
|
1
|
8
|
1
|
Therefore, flux direction are: sign{J2, J3, J5} > 0, and sign {J1, J4} {−, 0, +}.
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.
%have to check T-constraint for generated sign pattern.
k = 1;
match = [];for j = 1:nc_in
for i = 1:n_all
if C_in(i,j) ~= 0 && Sign_Pattern(i) ~= -11
if C_in(i,j)*Sign_Pattern(i) < 0
match(k) = 1;
else
match(k) = -1;
end
k = k + 1;
end
end
if length(match) > 0 && abs(sum(match)) ==
sum(ones(length(match),1))
cycle(j) = 1;
disp(‘found a possible cycle in sign pattern’);
else
cycle(j) = -1;
end
end
cycle
For this example network, it will output:
“found a possible cycle in sign pattern!” and
cycle = [ 1 -1 -1]
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.