Flux Balance Analysis with Mathematica

Flux Balance Analysis with Mathematica

Flux Balance Analysis with Mathematica

You are chief technical officer of a new biomedical research company, LF enterprizes, with a new blockbuster drug, code named X, on your hands - a cure for cancer and a Viagra substitute all in one. You serendipitously discoverd that the drug can be easily manufactured in E. coli by splicing in a gene that is highly expressed in the pituitary gland of the Bonobo. This is represented by the rection Rax in the diagram below. During a test run, the yield of X was quite good, but you need it make it even better, since your chief financial officer just informed you of her plans to do a leveraged buy out of Pfizer.

One thing you can try is to increase the yield of X by adjusting the glucose and oxygen flux rates on the bioreactor, which have the symbols O and G in the diagram. Your options range from 0 to 60M/kgDW-hr for O, and from 0 to 40M/KgDW-hr for G. By inspecting the diagram below, do you think its possible to cut costs without sacrificing yield by giving less than the maximum of glucose and/or oxygen? Might too much of either actually reduce production of X? Dont worry if you can't right now - that's what Flux Balance Analysis is for!

R stands for gRowth, and we assume that E. coli is willing and able to optimize the flux(es) going into R. The willing part means that we have succeeded in represneting E. coli's best interests as the node R in this model of metabolic fluxes. The able part means that E. coli has sufficent regulatory control over the set fluxes in this model to find the model optimum for R.  We can calculate the model optimum using the mass balance equations for each reaction, plus additional constraints that reflect the fact that metabolites must be conserved over all of the reactions. These intra and inter-reaction conservation equations are all linear. The optimization function, simply the magnitude of Rar, is also linear, thus we can use the Mathematica LinearProgramming function, that is, if you have Mathemaica version 4.2. If you dont, bug your TF or figure out how to optimize the equations using the Excel Solver, which is both powerful and end-user friendly.

Its convenient to use a spreadsheet to  set up the equations for the flux model.  The workbook Flux.xls contains the matrix of reactions and metabolites shown below, along with some custom worksheet functions that create matrix definitions that are needed by the Mathematica LinearProgramming function. The matrix represented by the light grey box of cells below is called the Stoichiometric matrix. The reaction mass balance equations are represented by the columns, and the metabolite balance equations are in the rows. The rel, or relation column below can have the values (-1, 0, 1), which is what the LinearProgramming function accepts as codes for (<, =, >). Notice that the thresh colum contains 0 for all of the rows of the Stoichiometric matrix - this is saying that metabolites must be conserved across all reactions. The rows with headers C_G and C_O are the constraint equations for Glucose and Oxygen, and on the last column is our objective function.

The worksheet functions below create the following Matematica matrix definitions from the table above:
ST_matrix - the Stoichiometric matrix, 
OBJ - the objective function, the last row in the excel table and the first argument of the LinearProgramming function. These constraints are set to maximize growth, or as far as the LinearProgramming function is concerned, we are minimizing negative growth, which is why the OBJ_F coefficent for Rar is -1.
LP_matrix - the second argument of the LinearProgramming function. Its all of the columns in the Stoichiometric matrix, intersected with the Stoichiometric matrix rows and the rows containing the metabolite influx constraints, C_G and C_O.
Constraints - the third argument of the LinearProgramming function. Two columns, a threshold level and an equality type code. The row count must match the number of rows in LP_matrix.

Set the matrix definitions to variables, plug the variables into the LinearProgramming function, and turn the crank...
The result is that optimum growth rate of 18 is achieved, and it requires (big suprise) the maximum allowable O and G levels. Notice that we also produce 2 units of our drug in this scenario as well.
The NullSpace function was applied to the Stoichiometric matrix. The result is a set of vectors that span the space of solutions that satisfy all of the constraints. One type of constraint that we did not consider yet is the fact that some of the reaction fluxes may have to be positive. Many reactions are strictly positive, since the reverse reaction may be very slow. Also remember that the true reverse of a given reaction has the same mass balance.

sToik  =  {{0, 1, 0, 0, -1, 0, 0, 0}, {0, 0, 0, 0, 
              1, -1, -1, 0}, {0, 0, -1, 0, 0, 1, 0, 0}, {0, 0, 
              0, -1, -2, 0, 5, 0}, {0, 0, 0, 1, 2, 0, -5, 0}, {0, 
              0, 0, -1, 0, 0, 4, -1}, {0, 0, 0, 1, 0, 0, -4, 
              1}, {1, 0, 0, -1, 0, 0, 0, 0}};
lP1 =  {0, 0, 0, 0, 0, 0, -1, 0};
lP2  = {{0, 1, 0, 0, -1, 0, 0, 0}, {0, 0, 0, 0, 
              1, -1, -1, 0}, {0, 0, -1, 0, 0, 1, 0, 0}, {0, 0, 
              0, -1, -2, 0, 5, 0}, {0, 0, 0, 1, 2, 0, -5, 0}, {0, 
              0, 0, -1, 0, 0, 4, -1}, {0, 0, 0, 1, 0, 0, -4, 
              1}, {1, 0, 0, -1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 
              0}, {1, 0, 0, 0, 0, 0, 0, 0}};
lP3  =  {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 
              0}, {0, 0}, {20, -1}, {50, -1}};
NullSpace[sToik]
LinearProgramming[lP1, lP2, lP3]

In the cell below, we determine the combinations of our two nullspace spanning vectors that produce all positive fluxes.

Clear[AllPlus]
AllPlus[
          a_]   :=  
Block[
{br}, 
br  =  1; 
For[i  =  1,  
            i  ≤  Length[a],  i++, 
If[
              a[[i]]  <  0,  
              br = 0, , ]
];  
br
]
a  =  {-2, 1, 1, -2, 1, 1, 0, 2};
b  =  {8, 1, -1, 8, 1, -1, 2, 0};
VVs  = 
          Table[AllPlus[x a  +  y b],  {x,  -1,  5,  
              0.2},  {y,  -1,  5,  0.2}];
ListDensityPlot[VVs,  PlotLabel -> "(a,b)", 
          MeshRange -> {{-1,  5},  {-1, 5}}];

Set the matrix definitions to variables, plug the variables into the LinearProgramming function, and turn the crank...
The NullSpace function was applied to the Stoi matrix. The result is a set of vectors that span the space of solutions that satisfy all of the constraints. One type of constraint that we did not consider yet is the fact that some of the reaction fluxes may have to be positive. Many reactions are strictly positive, since the reverse reaction rate may be very slow and have a different mass balance.

Dont forget that our objective was to maximize drug production, so we will assume that our model is valid for glucose and oxygen levels between (0, 40) and (0, 70) respectively, and we will see how much drug is produced assuming that E. coli will adjust it's fluxes to maximize growth in each input flux scenario. This is what we do with the FluxScan function below

lP1 =  {0, 0, 0, 0, 0, 0, -1, 0};
lP2  = {{0, 1, 0, 0, -1, 0, 0, 0}, {0, 0, 0, 0, 
              1, -1, -1, 0}, {0, 0, -1, 0, 0, 1, 0, 0}, {0, 0, 
              0, -1, -2, 0, 5, 0}, {0, 0, 0, 1, 2, 0, -5, 0}, {0, 
              0, 0, -1, 0, 0, 4, -1}, {0, 0, 0, 1, 0, 0, -4, 
              1}, {1, 0, 0, -1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 
              0}, {1, 0, 0, 0, 0, 0, 0, 0}};
FluxScan[mG_,  
            mO_]  :=  
Table[LinearProgramming[lP1,  
                lP2,  {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0},   {0, 
                    0}, {0, 0}, {0, 
                    0}, {i, -1}, {j, -1}}][[3]],  {i,  1,  
              mG},  {j,  1,  mO}];
ListPlot3D[FluxScan[40,  70]]

And the results are both unquestionable as well as exciting -  LF enterprizes will become very profitable if they grow E. coli anerobicaly! Time to hurry along now and collect my consulting fee....