Linear Programming with ScrumPy

1. First Steps

In order to carry out linear programming investigations on a model, it is first neccessary to generate a linear program associated with the model:

   1 lp = m.GetLP()

Where m is an existing ScrumPy model and lp is an instance an object representing a linear program set to minimise a set of fluxes subject to the steady-state and irreversibility constraints. No fluxes are are added to the objective function, and no other constraints are set at this point (other than steady-state). Linear programming investigations of the model are then carried out by manipulating the lp object.

2. Setting Constraints

The initial lp object has two sets of constraints set: the steady state assumption, and irreversible reactions to carry positive or zero flux. Reversible reactions are handled internally and there is no need to explicitly split reversible reactions into their forward and reverse components.

However, no flux constraints are set, nor is an objective function, both of which must done to carry out any useful work.

The most general way in which flux bounds may be set is with the lp.setFluxBounds(BDict) method, where BDict is a python dictionary in which the keys are reaction names and the values are a tuple describing the lower and upper flux bound of that reaction, eg:

   1 lp = m.GetLP()
   2 
   3 Bounds = {
   4      "Glucose_tx" : (0, 1),
   5      "CO2_tx": (-1,0),
   6      "EtOH_tx" : (-2,0)
   7 }
   8 
   9 lp.SetFluxBounds(Bounds)

In this example we are using the convention (NOT enforced by ScrumPy) that _tx appended to a metabolite name represents a transport reaction for that metabolite and that positive flux through a transporter represents import of the metabolite, and conversly negative flux represents export.

Thus in the above example we are sepcifying that the maximum rate of glucose uptake is one flux unit, the maximum rate of CO2 production is one flux unit and the maximum rate of ethanol production is 2 flux units. Note that these limits cannot be reached simultaneously (unless there is an error in the model !).

Of course any reaction, not just transporters can have bounds set in a similar fashion. Flux bounds may be changed or removed at any time. Instead of a numerical value, a bound may specified with the python constant None, this has the effect of removing the upper or lower bound (reversibility considerations still being taken into account). Setting a negative bound on an irreversible reaction will generate and exception (and otherwise ignored).

It often desired to set one or more fluxes to a fixed value ,and although this could be achieved using SetFluxBounds() with lower bound equal to the upper bound there is also the convenience function SetFixedFlux(FDict), e.g.

   1 Fluxes = {
   2      "Glucose_tx" : 1,
   3      "EtOH_tx" : -2
   4 }
   5 
   6 lp.SetFixedFlux(Fluxes)

This also allows negative flux to be set on an irreversible reaction, allthough a warning will be given. This should be used only very sparingly and with care, as it may lead to numerical instability later.

A constraint may be conveniently removed from a reaction with the ClearFluxConstraint(reac) method:

   1 lp.ClearFluxConstraint("Glucose_tx")

3. Setting the Objective Function

Firstly it must be specified if the optimisation is a minimisation or maximisation:

   1 lp.SetObjDirec("Min") # minimisation
   2 lp.SetObjDirec("Max") # maximisation
   3 lp.SetObjDirec()      # default - minimisation

If any argument other than the string "Min" or "Max" is given it will be ignored and an exception generated. If the direction is not specified, minimisation will be used by default.

Note that here minimisation and maximisation refers to the ABSOLUTE value of the flux in the case of reversible reactions.

Next one must specify which reactions are to be minimised or maximised, this achieved with the SetObjective(ReactionList) method, e.g.

   1 Fluxes = {
   2      "Glucose_tx" : 1,
   3      "EtOH_tx" : -2
   4 }
   5 
   6 lp.SetFixedFlux(Fluxes)
   7 lp.SetObjective(["CO2_tx"]) # minimise CO2 loss

Note that by convention "_tx" is used to denote a transporter defined such that negative flux indicates loss of the metabolite from the system.

4. Solving the LP

Once lp has been set up, this is as simple as

   1 lp.Solve()
   2 sol = lp.GetPrimSol()
   3 len(sol)
   4 print (sol)

Where sol is a dictionary of reaction name and fluxes.

5. Filtering solutions

For example to examine only the flux in transport reactions in the solution:

   1 for r in sol:
   2     if r.endswith("_tx"):
   3         print (r, sol[r])
   4 
   5 # or using python's buitin filter and lambda functions:
   6 for r in filter(lambda s: s.endswith("_tx"), sol):
   7     print (r, sol[r])
   8 
   9 NEED PYTHON COMPREHENSION