Linear Programming with ScrumPy

First Steps

In order to carry out linear programming investigations on a model, it is first neccessary to generate a specific representation of the model using ScrumPy's LP module:

   1 import LP
   2 my_lp = LP.lp(m)

Where m is an existing ScrumPy model and my_lp is an instance of the lp class, defined by the LP module. Linear programming investigations of the model are then carried out by manipulating my_lp, once the lp object has been generated, no further reference to the LP module is required.

Setting Constraints

The initial lp object has two sets of constaints set: the steady state assumption, and iireversible reactions to carry positive or zero flux. Reversible reactions are handled internally and ther 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 = ScrumPy.Model(m)
   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 lp = ScrumPy.Model(m)
   2 Fluxes = {
   3      "Glucose_tx" : 1,
   4      "EtOH_tx" : -2
   5 }
   6 
   7 lp.SetFluxBounds(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 UnBoundFlux(reac) method:

   1 lp.UnboundFlux("Glucose_tx")

Setting the Objective Function

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

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

If any argument other than the string "Min" or "Max" is given it will be ignored and an exception generated.

Note that here minimisation and maximumisation referes 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 lp = ScrumPy.Model(m)
   2 Fluxes = {
   3      "Glucose_tx" : 1,
   4      "EtOH_tx" : -2
   5 }
   6 
   7 lp.SetFluxBounds(Fluxes)
   8 lp.SetObjDirec()  # we will be minimising
   9 lp.SetObjective(["CO2_tx"]) # attempt to minimise CO2 loss

(More here later, weighted objectives.)

Solving the LP

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

   1 lp.Solve()