%
% wiech3.m
%
% setup and solve the cumomer equations for the
% example net in    
%	Bidirectional Reaction Steps in
%	Metabolic Networks: III
%

syms F1 F2 F3 F3 F4 F5 F6 F7 B2 B3	% the flux terms
syms ax1 a1x a11			% the input terms

% setup the 1-cumomer system  : equation (15)

A1 = [F2+F3+F5 0 0 -B3 0 0 0 0 0 0 -B2;
      0 F2+F3+F5 -B3 0 0 0 0 0 0 -B2 0;
      0 -F3 B3+F6 0 0 0 0 0 0 0 0;
      -F3 0 0 B3+F6 0 0 0 0 0 0 0;
      0 0 0 0 B3+F6 0 0 0 0 -F3 0;
      0 0 0 0 0 B3+F6 0 0 0 0 -F3
      0 0 0 -F6 0 0 F7 0 0 0 0;
      0 0 0 0 -F6 0 0 F7 0 0 0;
      0 0 0 0 0 -F6 0 0 F7 0 0;
      0 -F2 0 0 -B3 0 -F7 0 0 B2+F3+F4 0;
      -F2 0 0 0 0 -B3 0 -F7 0 0 B2+F3+F4];

r1 = [ax1*F1 a1x*F1 0 0 0 0 0 0 0 0 0].';

c1 = A1 \ r1;		% solve for the 1-cumomers

% setup the 2-cumomer system  : equation (16)

A2 = [F2+F3+F5 -B3 0 0 0 0 0 0 0 0 -B2;
      -F3 B3+F6 0 0 0 0 0 0 0 0 0;
      0 0 B3+F6 0 0 0 0 0 0 0 0;
      0 0 0 B3+F6 0 0 0 0 0 0 0;
      0 0 0 0 B3+F6 0 0 0 0 0 0;
      0 0 0 0 0 B3+F6 0 0 0 0 0;
      0 0 0 0 0 0 B3+F6 0 0 0 -F3;
      0 0 0 0 -F6 0 0 F7 0 0 0;
      0 0 0 0 0 -F6 0 0 F7 0 0;
      0 0 0 0 0 0 -F6 0 0 F7 0;
      -F2 0 0 0 0 0 -B3 -F7 0 0 B2+F3+F4];

r2 = [a11*F1;
      0;
      c1(2)*c1(10)*F3;
      c1(2)*c1(11)*F3;
      c1(1)*c1(10)*F3;
      c1(1)*c1(11)*F3;
      0;
      0;
      0;
      0;
      0];

c2 = inv(A2)*r2; 	% solve for the 2-cumomers (beware of backslash)

% setup the 3 cumomer system  :  equation (17)

A3 = [B3+F6 0 0 0 0;
      0 B3+F6 0 0 0;
      0 0 B3+F6 0 0;
      0 0 0 B3+F6 0;
      0 0 0 -F6 F7];

r3 = [c2(1)*c1(10)*F3;
      c2(1)*c1(11)*F3;
      c1(2)*c2(11)*F3;
      c1(1)*c2(11)*F3;
      0];

c3 = inv(A3)*r3;

% setup the 4 cumomer system  :  equation (18)

c4 = c2(1)*c2(11)*F3/(B3+F6);