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

Set i /A,B,C,D,E,byp,cof/;
Set j /v1,v2,v3,v4,v5,v6,b1,b2,b3/;
Parameter S(i,j)
/A.v1 -1
B.v1 1
B.v2 -2
C.v2 1
byp.v2 1
B.v3 -2
D.v3 1
cof.v3 -1
byp.v3 1
D.v4 -1
E.v4 1
cof.v4 1
C.v5 -1
D.v5 1
cof.v5 -1
C.v6 -1
E.v6 1
A.b1 1
E.b2 -1
byp.b3 -1/
LowerLimits(j)
UpperLimits(j);

LowerLimits(j)=0;
UpperLimits(j)=10;

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  /0.003/
epsilon /0.0001/
wL(j)
wU(j);

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

Variable minnumber;
Binary Variable y(j);

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(j)
l_restriction(j)
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(j, y(j));
u_restriction(j).. v(j)-y(j)*(UpperLimits(j)-wU(j))=l=wU(j);
l_restriction(j).. v(j)-y(j)*(LowerLimits(j)-wL(j))=g=wL(j);
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('b2')=1;
v.fx('v4')=0;
solve FBA using lp maximizing Obj;
wildtype_v(j)=v.l(j);

wU(j)=wildtype_v(j)+delta*abs(wildtype_v(j))+epsilon;
wL(j)=wildtype_v(j)-delta*abs(wildtype_v(j))-epsilon;
v.lo(j)=LowerLimits(j);
v.up(j)=UpperLimits(j);
*This section calculates the MOMA solution for the appropriate knockout
*indicated by the line v.fx('rxnname')=0;
*It also calculates the FBA solution for this same knockout
v.fx('v6')=0;
solve MOMA using nlp minimizing distance;
mutant_moma(j)=v.l(j);
solve ROOM using mip minimizing minnumber;
mutant_room(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 /ShlomiExample_Results.txt/;
result.pw=500; result.pc=5; result.ps=130;
put result;
put "Model status:  ";
put ROOM.modelstat;
put /;
put "Solver status:  ";
put ROOM.solvestat;
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;