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

*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_pi_e')=-Vmax;
UpperLimits('ATPM')=7.6;
LowerLimits('ATPM')=7.6;
UpperLimits('EX_glc_e')=0;

*Scaling Factor Introduced By Varma and Palsson AEM 1994
S(i,'Biomass')=1.3*S(i,'Biomass');

sets timesteps /time1*time121/
media /biomass,glucose,acetate,form,lactate,ethanol,succinate/;

Parameter
c(j) used to define the objective function for FBA
time(timesteps)
oxygenuptake  mmol per gDW per hour (max is 15) /0/
glucoseuptake mmol per gDW per hour (max is 10.5 aerobic and 18.5 anaerobic) /18.5/
acetateuptake mmol per gDW per hour /3.1/
initialglucose mmol per liter /11/
initialacetate mmol per liter /0/
initialbiomass gDW per liter /0.02/
dt units of hours /0.1/
concentrations(timesteps,media) mmol per liter;

c('Biomass')=1;
LowerLimits('EX_o2_e')=-oxygenuptake;

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 FBA /massbalance, calcobj/;

v.lo(j)=LowerLimits(j);
v.up(j)=UpperLimits(j);
concentrations('time1','biomass')=initialbiomass;
concentrations('time1','glucose')=initialglucose;
concentrations('time1','acetate')=initialacetate;

loop(timesteps,
    if(ord(timesteps)>1,time(timesteps)=(ord(timesteps)-1)*dt;
         v.lo('EX_glc_e')=max(-concentrations(timesteps-1,'glucose')/dt/concentrations(timesteps-1,'biomass'),-glucoseuptake);
         v.lo('EX_ac_e')=max(-concentrations(timesteps-1,'acetate')/dt/concentrations(timesteps-1,'biomass'),-acetateuptake);
         solve FBA using lp maximizing Obj;
         if(FBA.modelstat=1,
          concentrations(timesteps,'biomass')=concentrations(timesteps-1,'biomass')*exp(v.l('Biomass')*dt);
          concentrations(timesteps,'glucose')=max(concentrations(timesteps-1,'glucose')-v.l('EX_glc_e')/v.l('Biomass')*concentrations(timesteps-1,'biomass')*(1-exp(v.l('Biomass')*dt)),0);
          concentrations(timesteps,'acetate')=max(concentrations(timesteps-1,'acetate')-v.l('EX_ac_e')/v.l('Biomass')*concentrations(timesteps-1,'biomass')*(1-exp(v.l('Biomass')*dt)),0);
          concentrations(timesteps,'form')=concentrations(timesteps-1,'form')-v.l('EX_for_e')/v.l('Biomass')*concentrations(timesteps-1,'biomass')*(1-exp(v.l('Biomass')*dt));
          concentrations(timesteps,'lactate')=concentrations(timesteps-1,'lactate')-v.l('EX_lacD_e')/v.l('Biomass')*concentrations(timesteps-1,'biomass')*(1-exp(v.l('Biomass')*dt));
          concentrations(timesteps,'ethanol')=concentrations(timesteps-1,'ethanol')-v.l('EX_etoh_e')/v.l('Biomass')*concentrations(timesteps-1,'biomass')*(1-exp(v.l('Biomass')*dt));
          concentrations(timesteps,'succinate')=concentrations(timesteps-1,'succinate')-v.l('EX_succ_e')/v.l('Biomass')*concentrations(timesteps-1,'biomass')*(1-exp(v.l('Biomass')*dt));
         else
          concentrations(timesteps,'biomass')=concentrations(timesteps-1,'biomass');
          concentrations(timesteps,'glucose')=concentrations(timesteps-1,'glucose');
          concentrations(timesteps,'acetate')=concentrations(timesteps-1,'acetate');
          concentrations(timesteps,'form')=concentrations(timesteps-1,'form');
          concentrations(timesteps,'lactate')=concentrations(timesteps-1,'lactate');
          concentrations(timesteps,'ethanol')=concentrations(timesteps-1,'ethanol');
          concentrations(timesteps,'succinate')=concentrations(timesteps-1,'succinate');
          );
    );
);
*This section writes out the results to a txt file
file result /BatchResultsAnaerobic.txt/;
result.pw=500; result.pc=5; result.ps=130;
put result;
put "Time(hours)";
put "Biomass(gDW/L)";
put "Glucose (mmol/L)";
put "Acetate (mmol/L)";
put "Formate (mmol/L)";
put "Lactate (mmol/L)";
put "Ethanol (mmol/L)";
put "Succinate (mmol/L)";
put /;
loop(timesteps,
    put time(timesteps):8:4;
    loop(media, put concentrations(timesteps,media):8:4;);
    put /;
);
putclose result;

*This will write the non-zero values for v and Obj in the LST file
display concentrations;