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

*Read in the appropriate S matrix
$include CoreTextbookModel.gms

*Place limits on the exchange fluxes based on the minimal media
*for a negative flux through the exchange reactions implies that
*the metabolites are being taken up or consumed by the cell.
*By default the upperlimits of the exchange fluxes are all set to
*the Vmax, indicating that the cell can secrete any of the extracellular
*metabolites
UpperLimits(j)=Vmax;
*CARBON SOURCE: select upper and lower limits for exchange flux
LowerLimits('EX_glc_e')=-5;
UpperLimits('EX_glc_e')=-5;
*allow co2,pi,o2,h,h2o to be taken up by the cell
LowerLimits('EX_co2_e')=-Vmax;
LowerLimits('EX_h2o_e')=-Vmax;
LowerLimits('EX_h_e')=-Vmax;
LowerLimits('EX_o2_e')=0;
LowerLimits('EX_pi_e')=-Vmax;

Parameter
c(j) used to define the objective function for FBA
store_maxs(j) stores the maximum value for each flux
store_mins(j) stores the minimum value for each flux;

Variables
v(j) flux values through reaction in network
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( j,S(i,j)*v(j) )=e=0;
calcobj.. Obj=e=sum( j,c(j)*v(j) );

Model fluxvariability /massbalance, calcobj/;

alias (j,duplicate_j);
v.lo(j)=LowerLimits(j);
v.up(j)=UpperLimits(j);

*optimize for a particular flux and fix level
*NOTE: if you want to look at variability across
*all solutions then comment out the next three lines
*so that the optimized flux isn't fixed to the optimal level
c('Biomass')=1;
solve fluxvariability using lp maximizing Obj;
v.fx('Biomass')=Obj.l;
c(j)=0;

loop(duplicate_j,c(duplicate_j)=1;
         solve fluxvariability using lp maximizing Obj;
         store_maxs(duplicate_j)=Obj.l;
         solve fluxvariability using lp minimizing Obj;
         store_mins(duplicate_j)=Obj.l;
         c(duplicate_j)=0;  );

*This section writes out the results to a txt file
file result /maxs_mins_fluxes.txt/;
result.pw=500; result.pc=5; result.ps=130;
put result;
loop(j,put j.tl put store_maxs(j):8:4;put store_mins(j):8:4; put /);
putclose result;

display store_maxs;
display store_mins;