**THIS CODE WAS WRITTEN FOR THE CONSTRAINT-BASED MODELING WORKSHOP BY J.REED (1/2008)
**Due to GAMS demo limits not all reactions are included in the ROOM objective function metric
**it (excludes biomass, ATPM, and transport/exchange reactions)

*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')=0;
*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')=-Vmax;
LowerLimits('EX_pi_e')=-Vmax;

*Define reactions that are used in the ROOM objective function
Set subj(j) /ACKr,ACONT,ADHEr,ADK1,ATPS4r,CS,CYTBD,ENO,FBA,FBP,FRD,FUM,G6PDH2r,GAPD,GND,ICDHyr,ICL,LDH_D,MALS,MDH,ME1,ME2,NADH11,PDH,PFK,PFL,PGI,PGK,PGL,PGM,PPC,PPCK,PPS,PTAr,PYK,PYRt2r,RPE,RPI,SUCD1i,SUCD4,SUCOAS,TALA,AKGDH,NADTRHD,THD2,TKT1,TKT2,TPI/;

Parameter
c(j)           used to define the objective function for FBA
wildtype_v(j)  used to store wildtype FBA fluxes
mutant_room(j) used to store mutant MOMA fluxes
mutant_fba(j)  used to store mutant FBA fluxes
mutant_moma(j) used to store mutant MOMA fluxes
delta  used to indicate what flux changes are significant (ROOM)
epsilon used to indicate what flux changes are significant (ROOM)
wL(subj) used to indicate what flux changes are significant (ROOM)
wU(subj) used to indicate what flux changes are significant (ROOM);

delta = 0.03;
epsilon = 0.001;

Variables
v(j) flux values through reaction in network
Obj  objective function value for the FBA solution
distance objective function value for for the MOMA solution;

Variable minnumber objective function value for ROOM solution ie. number of significantly altered fluxes;
Binary Variable y(subj) binary variables 1 means that the flux changes significantly;

Equations
massbalance(i) mass balance equations for each metabolite
calcobj calculates the dot product of the c vector the flux vector
sum_y calculates the euclidean distance between mutant fluxes and wt fluxes
u_restriction(subj) upper flux bounds for ROOM
l_restriction(subj) lower flux bounds for ROOM
eucl_distance calculates the euclidean distance between mutant fluxes and wt fluxes;

massbalance(i).. sum( j,S(i,j)*v(j) )=e=0;
calcobj.. Obj=e=sum( j,c(j)*v(j) );
sum_y.. minnumber=e=sum(subj, y(subj));
u_restriction(subj).. v(subj)-y(subj)*(UpperLimits(subj)-wU(subj))=l=wU(subj);
l_restriction(subj).. v(subj)-y(subj)*(LowerLimits(subj)-wL(subj))=g=wL(subj);
eucl_distance.. distance=e=sum(j, sqr(wildtype_v(j)-v(j)));

Model FBA /massbalance, calcobj/;
Model ROOM /massbalance,sum_y,u_restriction,l_restriction/;
Model MOMA /massbalance,eucl_distance/;

******************************************************************************************************************
*This section calculates the FBA solution for maximizing biomass
*for the wildtype strain and stores the fluxes in the wildtype_v parameter
v.lo(j)=LowerLimits(j);
v.up(j)=UpperLimits(j);
c('Biomass')=1;
solve FBA using lp maximizing Obj;
wildtype_v(j)=v.l(j);

*Defines allowable variation before becoming significant for ROOM calculations
wU(subj)=wildtype_v(subj)+delta*abs(wildtype_v(subj))+epsilon;
wL(subj)=wildtype_v(subj)-delta*abs(wildtype_v(subj))-epsilon;

v.fx('ACONT')=0;
**********************************************************************************************************************
*This section calculates the ROOM and MOMA solutions for the appropriate knockout
*indicated by the line v.fx('rxnname')=0;
*It also calculates the FBA solution for this same knockout
solve ROOM using mip minimizing minnumber;
mutant_room(j)=v.l(j);

solve MOMA using nlp minimizing distance;
mutant_moma(j)=v.l(j);
solve FBA using lp maximizing Obj;
mutant_fba(j)=v.l(j);


*This section writes out the results to a txt file
file result /GeneDeletionPredictions.txt/;
result.pw=500; result.pc=5; result.ps=130;
put result;
put "MOMA Solver status:  ";
put "ROOM Solver status:  ";
put "FBA Solver status:  ";
put /;
put MOMA.solvestat;
put ROOM.solvestat;
put FBA.solvestat;
put /;
put "MOMA Objective:  ";
put "ROOM Objective:  ";
put "FBA Objective:  ";
put /;
put distance.l;
put minnumber.l;
put Obj.l;
put /;
put "Flux_Name"; put "FBA_Wildtype"; put "MOMA_Mutant";  put "ROOM_Mutant"; put "FBA_Mutant"; put /;
loop(j,put j.tl; put wildtype_v(j):8:4; put mutant_moma(j):8:4; put mutant_room(j):8:4; put mutant_fba(j):8:4;put /);
putclose result;