**THIS CODE WAS WRITTEN FOR CBE782 J.REED (4/2011)
$onecho > cplex.opt
eprhs 1e-8
epopt 1e-8
epint 1e-8
$offecho
$offlisting
*Read in the appropriate S matrix
$include EcoliCoreTextbookModel.gms
Sets conditions conditions for gene expression datasets /aerobic,anaerobic/
select_condition(conditions) condition to analyze /aerobic/
genes list of genes in the model /aceA,aceB,ackA,acnA,acnB,adhE,adk,aceE,aceF
atpA,atpB,atpC,atpD,atpE,atpF,atpG,atpH,atpI
crr,cydA,cydB,dctA,dcuC,eno,fba,fbp,focA,frdA,frdB,frdC,frdD
fumA,fumB,fumC,gapA,gltA,gnd,icdA,ldhA,lpdA,maeB,mdh
nuoA,nuoB,nuoE,nuoF,nuoG,nuoH,nuoI,nuoJ,nuoK,nuoL,nuoM,nuoN
pckA,pfkA,pfkB,pflA,pflB,pflC,pflD,pgi,pgk,pgl,pgm,pitA,pitB
pntA,pntB,ppc,ppsA,pta,ptsG,ptsH,ptsI,pykA,pykF,rpe,rpi
sdhA,sdhB,sdhC,sdhD,sfcA,sucA,sucB,sucC,sucD,talA,tktA,tktB,tpi,zwf/;
Parameter high_cutoff /2000/, low_cutoff /500/;
*Read in gene expression data
$include GeneExpressionData.gms
*Determine which reactions belong to set Rhigh and Rlow
$include MapData2Rxns.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')=-9.02;
UpperLimits('EX_glc_e')=-9.02;
LowerLimits('EX_o2_e')=-14.93;
UpperLimits('EX_o2_e')=-14.93;
LowerLimits('EX_ac_e')=4.15;
UpperLimits('EX_ac_e')=4.15;
*allow co2,pi,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;
LowerLimits('ATPM')=7.6;
S(i,'Biomass')=1.3*S(i,'Biomass');
Variables
v(j) flux values through reaction in network
Obj this is the value of the objective function for the FBA solutions
growth growth rate;
Binary Variables x(j),y(j);
Parameter epsilon /0.5/;
Equations
massbalance(i) mass balance equations for each metabolite
calcobj calculates the number of fluxes which match the gene expression data
maxbiomass calculate the growth rate (ie. biomass flux)
Rh_lower(j) lower bound for reactions in Rhigh
Rh_upper(j) upper bounds for reaction in Rhigh
Rl_lower(j) lower bounds for reactions in Rlow
Rl_upper(j) upper bounds for reactions in Rlow;
massbalance(i).. sum( j,S(i,j)*v(j) )=e=0;
calcobj.. Obj=e=sum(Rhigh,x(Rhigh)+y(Rhigh))+sum(Rlow,y(Rlow));
Rh_lower(Rhigh).. v(Rhigh) + y(Rhigh)*(LowerLimits(Rhigh) - epsilon)=g=LowerLimits(Rhigh);
Rh_upper(Rhigh).. v(Rhigh) + x(Rhigh)*(UpperLimits(Rhigh) + epsilon)=l=UpperLimits(Rhigh);
Rl_lower(Rlow).. v(Rlow) =g= LowerLimits(Rlow)*(1-y(Rlow));
Rl_upper(Rlow).. v(Rlow) =l= UpperLimits(Rlow)*(1-y(Rlow));
maxbiomass.. growth=e=v('Biomass');
*All reaction between LowerLimits and UpperLimits (this constrains Rmed fluxes)
v.lo(j)=LowerLimits(j);
v.up(j)=UpperLimits(j);
model Shlomi_GeneExp/all/;
Shlomi_GeneExp.optfile=1;
Shlomi_GeneExp.optcr=0;
Parameter maxgrowth(j), mingrowth(j);
file result /Shlomi_GeneExp_results.txt/;
result.pw=500; result.pc=5; result.ps=130;
put result;
*find a flux distribution that best matches gene expression patterns
Solve Shlomi_GeneExp using mip maximizing Obj;
Obj.fx=Obj.l;
put "Model status: "; put Shlomi_GeneExp.modelstat; put /;
put "Solver status: "; put Shlomi_GeneExp.solvestat; put /;
put "Number of Reactions whose flux matches expression patterns"; put Obj.l; put /; put/;
*find the alternate solution with the maximum growth rate
Solve Shlomi_GeneExp using mip maximizing growth;
maxgrowth(j)=v.l(j);
*find the alternate solution with the maximum growth rate
Solve Shlomi_GeneExp using mip minimizing growth;
mingrowth(j)=v.l(j);
put "Reaction"; put "FluxValue"; put "ReactionStatus Values"; put /;
loop(j,put j.tl; put maxgrowth(j):8:4; put mingrowth(j):8:4; put reactionstatus(j); put /);
put /;
put "Highly Expressed Reactions (RH)"; put card(Rhigh); put /;
loop(Rhigh, put Rhigh.tl;); put /; put /;
put "Lowly Expressed Reactions (RL)"; put card(Rlow); put /;
loop(Rlow, put Rlow.tl;); put /; put /;
put "Mediumly Expressed Reactions or Reactions w/o Genes (RM)"; put card(Rmed); put /;
loop(Rmed, put Rmed.tl;); put /; put /;
putclose result;
*This will write the non-zero values for v and Obj in the LST file
display v.l,Obj.l;