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