MATLAB code for the NP-complete exhaustive calculation of the minimal cycle basis (MCB) for a biochemical network is given here. A MATLAB script demonstrating calculation of the MCB for a few example networks is also included.
Related references:
-
-
-
-
Our exhaustive algorithm for computing the minimal cycle basis for a stoichiometric system is given in MATLAB code. Using stoichiometric matrix S as an input, function cycles.m computes all minimal cycles, which are output columnwise in matrix N:
function [N] = cycles(S)
N = null(S,'r'); % Compute a basis for null space of S.
% Safeguard wrt numerical error.
epsilon = 1e-10;
N = N.*( abs(N) > epsilon );
[m,n] = size(N); % initial size of N.
i = 1;
while (i <= n)
[i n]
N = [N cycle_i(N,i)];
n = length(N(1,:));
i = i + 1;
end
%Before exiting check for minimality, it is not necessary in most times.
N = test_new_vectors(N,N,2);
Note that in MATLAB syntax N(:,i) denotes the i th column of N and N(i,:) denotes the i th row.
The first step of cycle algorithm uses an internal MATLAB function null to compute an algebraic basis for the null space of S. This basis is stored as the columns of the matrix
.
The function cycle_i.m generates all cycles with respect to i-th column of matrix N. Specifically, all linear combinations of the i-th with i+1 to size(N,2) columns of N are computed by the function linear_combos.m, and combined vectors are tested against the minimal cycle definition.
function [Nnew] = cycle_i(N,i)
n = length(N(1,:));
Nnew = [];
j = i+1;
while j <= n
% Compute new null space vectors from linear combinations of i^th and j^th columns of N.
Nnew = [Nnew linear_combos(N, N(:,i), N(:,j) )];
j = j+1;
end
% Test new null space vectors (Nnew) to see if they form cycles that are minimal against each other and vectors in N.
Nnew = test_new_vectors(Nnew,Nnew,2);
Nnew = test_new_vectors(N,Nnew,1);
function [Nnew] = linear_combos(N,v1,v2);
Nnew = [];
Nnew = [Nnew (v1+v2) (v1-v2)];
% Construct new linear combination so that its i-th entry be zero.
n = length(v1);
for i = 1:n
v1i = v1(i);
v2i = v2(i);
if (v1i ~= 0) & (v2i ~= 0) & ( (abs(v1i)~=1)|(abs(v2i)~=1) )
Nnew = [Nnew (v2i.*v1 - v1i.*v2)];
end
end
% Safeguard wrt numerical error.
epsilon = 1e-10;
Nnew = Nnew.*( abs(Nnew) > epsilon );
% guarantee cycles are minimal against each other and vectors in N.
Nnew = test_new_vectors(Nnew,Nnew,2); %with respect to Nnew itself.
Nnew = test_new_vectors(N,Nnew,1); % with respect to N.
This function constructs a set of putative entries of N from combinations of two vectors
and
, and stores these putative entries in the matrix Nnew. The first two vectors are computed as (
+
) and (
−
). Other vectors are computed as (
−
) where
and
are the i th entries of the vectors
and
. The vector (
−
) has a value of 0 for the ith entry.
The generated set of putative vectors stored in Nnew must next be tested for minimal support using the definition of matroid cycles. Such a test is perfprmed by calling the function test_new_vectors.m:
function [Nadd] = test_new_vectors(N,Nnew,option)
Nadd = [];
n = size(N,2); nnew = size(Nnew,2);
if n == 0 || nnew == 0
return
end
Ns = (N~=0); % Change to binary vector.
Nsnew = abs(sign(Nnew));
% Each cycle should has minimum signed supports.
for j = 1:nnew
if option == 1 % if N =/= Nnew
t = Ns'*Nsnew(:,j) - sum(Ns,1)';
else % if N = Nnew
t = Ns(:,j+1:n)'*Nsnew(:,j) - sum(Ns(:,j+1:n),1)';
end
if sum(t~=0) == length(t)
Nadd = [Nadd Nnew(:,j)];
end
end
In the above function, there are two options—one is to test for minimality in comparison to Nnew itself, the other is to test with respect to N matrix. It is always suggested that new putative vectors are tested for minimality again each other, and then tested with regard to the vectors in N matrix. All non-minimal vectors in Nnew are disqualified from adding to the set N. The complete code, along with text examples is available from the authors.
Few examples have been included in the driver file: cycle_driver.m