**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;