**THIS CODE WAS WRITTEN FOR THE CONSTRAINT-BASED MODELING WORKSHOP J.REED (2/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')=-18.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')=0;
LowerLimits('EX_pi_e')=-Vmax;
LowerLimits('ATPM')=7.6;
S(i,'Biomass')=1.3*S(i,'Biomass');
*Define the number of steps that you want to take eg. /step1*step25/ will have 25 steps
Sets
steps /step1*step15/
maxmin /maxprod,minprod/;
Parameter
c(j) used to define the objective function for optimization
n_steps number of steps that will be taken and is defined by the elements in steps
range_max maximum flux value through the flux to be varied
range_min minimum flux value through the flux to be varied
flux_value(steps) stores the values for the varied flux
store_obj(steps,maxmin) stores the value of the objective function for each iteration
pick_flux(j) a vector of zeros except for the one flux which will be varying;
*Determine the number of steps based on the number of elements in the set steps
n_steps=card(steps);
*Select the flux that you want to vary
pick_flux('Biomass')=1;
Variables
v(j) flux values through reaction in network
Obj this is the value of the objective function for the individual 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 FBA /massbalance, calcobj/;
v.lo(j)=LowerLimits(j);
v.up(j)=UpperLimits(j);
*calculate the allowable range for the chosen flux
loop(j, if(pick_flux(j)<>0,c(j)=1;); );
solve FBA using lp maximizing Obj;
range_max=Obj.l;
solve FBA using lp minimizing Obj;
range_min=Obj.l;
*reset the objective function to maximize objective of interest
c(j)=0;
c('EX_lacD_e')=1;
loop(steps,
flux_value(steps)=range_min+(ord(steps)-1)*(range_max-range_min)/(n_steps-1);
loop(j, if(pick_flux(j)<>0,v.fx(j)=flux_value(steps); ); );
solve FBA using lp maximizing Obj;
store_obj(steps,'maxprod')=Obj.l;
solve FBA using lp minimizing Obj;
store_obj(steps,'minprod')=Obj.l;
);
file result /Robustness_results.txt/;
result.pw=500; result.pc=5; result.ps=130;
put result;
put "MaxFluxValue";
put range_max:8:4;
put /;
put "MinFluxValue";
put range_min:8:4;
put /;
put "Number_of_Steps";
put n_steps;
put /; put /;
put "FluxValue";
put "MaxProduction";
put "MinProduction";
put /;
loop(steps,put flux_value(steps):8:4; put store_obj(steps,'maxprod'):8:4; put store_obj(steps,'minprod'):8:4; put /);
putclose result;
display flux_value,store_obj;