**THIS CODE WAS WRITTEN FOR THE CONSTRAINT-BASED MODELING WORKSHOP J.REED (2/2008)

*Read in the appropriate S matrix
sets i /A,B,C,D/
     j /v1*v7/;
Parameter
S(i,j)
/A.v1 1
A.v2 -1
B.v2 1
B.v3 -1
C.v3 1
A.v4 -1
D.v4 1
A.v5 1
D.v5 -1
D.v6 -1
C.v6 1
C.v7 -1/
Vmax;

Vmax=1000;
Sets
irr(j) /v1,v2,v3,v4,v5,v6,v7/;


Parameter
c_irr(irr) used to select which flux to optimize;

Variables
v_irr(irr) irreversible fluxes
Obj this is the value of the objective function for the FBA solutions;


Equations
massbalance(i) mass balance equations for each metabolite
calcobj calculates the dot product of the c vector the flux vector;
massbalance(i).. sum(irr,S(i,irr)*v_irr(irr)) =e=0;
calcobj.. Obj=e=sum(irr,c_irr(irr)*v_irr(irr));

Model FBA /massbalance, calcobj/;

v_irr.lo(irr)=0;
v_irr.up(irr)=Vmax;

******************************************************************
*This Section Finds All Blocked Reactions That Need to Be Removed*
******************************************************************
set unblocked_irr(j);
unblocked_irr(j)=no;
alias(dummyindex,irr);
loop(dummyindex, c_irr(dummyindex)=1;
        solve FBA using lp maximizing Obj;
        if(Obj.l=0,unblocked_irr(dummyindex)=no;
        else unblocked_irr(dummyindex)=yes;);
        c_irr(dummyindex)=0; );

display unblocked_irr;

******************************************************************
***This Section Calculates the Minimum and Maximum Flux Ratios****
******************************************************************
Parameters
d(j);

Variables
v_hat_irr(j) flux values through reaction in network
Obj_hat
t;

v_hat_irr.lo(j)=0;
v_hat_irr.up(j)=inf;
t.lo=0;

Equations
massbalance_hat(i) mass balance equations for each metabolite
calcobj_hat calculates the dot product of the c vector the flux vector
uptakelimit;

massbalance_hat(i).. sum(unblocked_irr,S(i,unblocked_irr)*v_hat_irr(unblocked_irr))=e=0;
calcobj_hat.. Obj_hat=e=sum(unblocked_irr,d(unblocked_irr)*v_hat_irr(unblocked_irr));
uptakelimit.. v_hat_irr('v1')=l=t*10;
Model FluxCoupling /massbalance_hat,calcobj_hat,uptakelimit/;

d(j)=0;
alias(dummyindex2,unblocked_irr);
alias(dummyindex3,unblocked_irr);

Parameter
store_max(j,j)
store_min(j,j)

loop(dummyindex2,
         v_hat_irr.fx(dummyindex2)=1;

         loop(dummyindex3,
                 d(dummyindex3)=1;
                 solve FluxCoupling using lp maximizing Obj_hat;
                 if (FluxCoupling.modelstat=3, store_max(dummyindex2,dummyindex3)=inf;
                    else store_max(dummyindex2,dummyindex3)=Obj_hat.l;);
                 solve FluxCoupling using lp minimizing Obj_hat;
                 store_min(dummyindex2,dummyindex3)=Obj_hat.l;
                 d(unblocked_irr)=0; );

         v_hat_irr.lo(dummyindex2)=0;
         v_hat_irr.up(dummyindex2)=inf;
);


*This section writes out the results to a txt file
file result /FluxCouplingExample.txt/;
result.pw=50000; result.pc=5; result.ps=130;
put result;
put "MAXRATIO";
loop(unblocked_irr,put unblocked_irr.tl;);
put /;
loop(unblocked_irr,put unblocked_irr.tl;
     loop(dummyindex3, put store_max(unblocked_irr,dummyindex3):8:4;);
     put /);

put /;
put "MINRATIO";
loop(unblocked_irr,put unblocked_irr.tl;);
put /;
loop(unblocked_irr,put unblocked_irr.tl;
     loop(dummyindex3, put store_min(unblocked_irr,dummyindex3):8:4;);
     put /);
putclose result;