This tutorial will walk you through the the code (from corm.py) that builds the COX-2 Reaction Model (CORM). It will provide a very basic introduction to the PySB modeling framework that CORM is encoded in as well as the biological interactions that CORM encodes.

The first lines in corm.py import some PySB code that will be needed for building the model:

In [1]:
from pysb import Model, Monomer, Parameter, Initial, Rule, Observable
from pysb.macros import bind, bind_complex, catalyze

The first line imports a number of objects (Python classes, if you're familiar with Python) that will provide building blocks for the model. The second imports several PySB macros. These are Python functions that allow you to quickly define a set of commonly seen biological interactions (like enzyme catalysis).

The first step in building a PySB model is to create a new model object:

In [2]:
Model()
Out[2]:
<Model '_interactive_' (monomers: 0, rules: 0, parameters: 0, expressions: 0, compartments: 0) at 0x10e8eb410>

We now have a new empty model object. The first thing we want to add to it are the monomers we will have in our model. 'Monomers' is a somewhat misleading term as PySB monomers are not necessarily protein monomers (although they can be). A PySB monomer is just the smallest interacting unit you want in your model. This might be anything from a large multi-protein complex (that, in the context of your model, never diassociates) to an actual protein monomer. In the case of CORM, we will define a monomer for the homodimeric enzyme COX-2:

In [3]:
Monomer('COX2', ['allo', 'cat'])
Out[3]:
Monomer('COX2', ['allo', 'cat'])

The first argument to the PySB Monomer class gives the name of the newly created monomer (COX2 in this case). The second argument (in square brackets) is a Python list of binding sites on the monomer. Here we have defined two binding sites for COX2, an allosteric site (named 'allo'), and a catalytic site (named 'cat'). Next we'll define monomers for the two substrates COX2 can bind in our model:

In [4]:
Monomer('AG', ['b'])
Out[4]:
Monomer('AG', ['b'])
In [5]:
Monomer('AA', ['b'])
Out[5]:
Monomer('AA', ['b'])

These are 2-arachidonoylglycerol (called AG in the model) and arachidonic acid (called AA in the model). Each of these has one generic binding site (called 'b'). Lastly, we'll define monomers for the products of COX2 turnover of AA or AG:

In [6]:
Monomer('PG')
Out[6]:
Monomer('PG')
In [7]:
Monomer('PGG')
Out[7]:
Monomer('PGG')

These are prostaglandin (called PG in the model) and prostaglandin glycerol (called PGG in the model). Because we will model the catalysis as a two-step process (E+S <-> E-S --> E + P), we will never have a chemical species in our model in which a product is bound to anything. Therefore, we only need to name PG and PGG but don't have to give them binding sites.

Now that we've defined the interacting components in our model, we need to define parameters for how much of each monomer we begin with. To do this we first use the PySB parameter class to define these parameters:

In [8]:
Parameter('COX2_0', 15e-3)
Out[8]:
Parameter('COX2_0', 0.015)
In [9]:
Parameter('AG_0', 16)
Out[9]:
Parameter('AG_0', 16.0)
In [10]:
Parameter('AA_0', 16)
Out[10]:
Parameter('AA_0', 16.0)
In [11]:
Parameter('PG_0', 0)
Out[11]:
Parameter('PG_0', 0.0)
In [12]:
Parameter('PGG_0', 0)
Out[12]:
Parameter('PGG_0', 0.0)

Each parameter is given a name (COX2_0, for instance) and a value. It is critical to ensure we are consistent with the units we use to define our parameters. In the case of CORM, these initial conditions are given in micromolar. When we define rate constants later for the reactions occuring in our model, we will be sure they are also in terms of micromolar (when concentration dependent).

While we've defined parameters for these initial conditions, we haven't yet told PySB that these parameters are initial conditions, or specified what chemical species each corresponds to. To do that, we use the Initial PySB class:

In [13]:
Initial(COX2(allo=None, cat=None), COX2_0)
In [14]:
Initial(AG(b=None), AG_0)
In [15]:
Initial(AA(b=None), AA_0)
In [16]:
Initial(PG(), PG_0)
In [17]:
Initial(PGG(), PGG_0)

The first argument to the Initial class is the species that will be initialized with that initial condition. For instance, for COX2, that species is COX2(allo=None, cat=None). The none after each of COX2's binding sites indicates that the enzyme should begin with nothing bound at either its allosteric or catalytic sites. The second argument to Initial is the parameter that gives the value of the initial condition (COX2_0, for instance). For COX2, this means we will start the model with a concentration of .015 microM.

Nex, we need to define rate parameters for each of the reactions occuring in our model. In the case of CORM, we will define different rate parameters for each of the reactions that are possible in the model, but note that this is not required. If we believed the same rate parameter accurately described the rate of different reactions, we could have used the same PySB parameter in multiple reactions.

First we will define rates for interactions of AA with COX2 when nothing is present in the allosteric site:

In [18]:
Parameter('kf_AA_cat1', 1000.0)
Out[18]:
Parameter('kf_AA_cat1', 1000.0)
In [19]:
Parameter('kr_AA_cat1', 830.0)
Out[19]:
Parameter('kr_AA_cat1', 830.0)
In [20]:
Parameter('kcat_AA1', 1.3)
Out[20]:
Parameter('kcat_AA1', 1.3)

We have now defined forward and reverse reaction rates for the binding of AA to COX2's catalytic site when COX2 has nothing bound in the allosteric site, and a kcat for COX2 turnover of AA when nothing is present in the allosteric site. Next, we'll define the same type of rates for AA and COX2 interactions when COX2 has an AG molecule bound in the allosteric site:

In [21]:
Parameter('kf_AA_cat2', 1.0e-3)
Out[21]:
Parameter('kf_AA_cat2', 0.001)
In [22]:
Parameter('kr_AA_cat2', 3.3e-6)
Out[22]:
Parameter('kr_AA_cat2', 3.3e-06)
In [23]:
Parameter('kcat_AA2', 2.3)
Out[23]:
Parameter('kcat_AA2', 2.3)

And again for AA and COX2 catalytic interactions when COX2 has a molecule of AA bound in the allosteric site:

In [24]:
Parameter('kf_AA_cat3', 1.0e-3)
Out[24]:
Parameter('kf_AA_cat3', 0.001)
In [25]:
Parameter('kr_AA_cat3', 8.3e-6)
Out[25]:
Parameter('kr_AA_cat3', 8.3e-06)
In [26]:
Parameter('kcat_AA3', 1.3)
Out[26]:
Parameter('kcat_AA3', 1.3)

Next, we'll define the analogous catalytic interactions for binding and turnover of AG by COX2 when COX2 has nothing, AG, or AA bound at the allosteric site:

In [27]:
Parameter('kf_AG_cat1', 1000.0)
Out[27]:
Parameter('kf_AG_cat1', 1000.0)
In [28]:
Parameter('kr_AG_cat1', 760.0)
Out[28]:
Parameter('kr_AG_cat1', 760.0)
In [29]:
Parameter('kcat_AG1', 1.2)
Out[29]:
Parameter('kcat_AG1', 1.2)
In [30]:
Parameter('kf_AG_cat2', 1.0e-3)
Out[30]:
Parameter('kf_AG_cat2', 0.001)
In [31]:
Parameter('kr_AG_cat2', 4.8e-4)
Out[31]:
Parameter('kr_AG_cat2', 0.00048)
In [32]:
Parameter('kf_AG_cat3', 1.0e-3)
Out[32]:
Parameter('kf_AG_cat3', 0.001)
In [33]:
Parameter('kr_AG_cat3', 1.9e-6)
Out[33]:
Parameter('kr_AG_cat3', 1.9e-06)
In [34]:
Parameter('kcat_AG3', 0.21)
Out[34]:
Parameter('kcat_AG3', 0.21)

Note that there is no 'kcat_AG2' parameter because that corresponds to COX2 turnover of AG when AG is also bound in the allosteric site, and since CORM assumes substrate inhibition, this rate would be zero.

Next we need to define reaction rates for binding of COX2 to AA at the allosteric site (with nothing, AA, or AG bound at the catalytic site):

In [35]:
Parameter('kf_AA_allo1', 1000.0)
Out[35]:
Parameter('kf_AA_allo1', 1000.0)
In [36]:
Parameter('kr_AA_allo1', 1.0e5)
Out[36]:
Parameter('kr_AA_allo1', 100000.0)
In [37]:
Parameter('kf_AA_allo2', 1000.0)
Out[37]:
Parameter('kf_AA_allo2', 1000.0)
In [38]:
Parameter('kr_AA_allo2', 1000.0)
Out[38]:
Parameter('kr_AA_allo2', 1000.0)
In [39]:
Parameter('kf_AA_allo3', 1000.0)
Out[39]:
Parameter('kf_AA_allo3', 1000.0)
In [40]:
Parameter('kr_AA_allo3', 250.0)
Out[40]:
Parameter('kr_AA_allo3', 250.0)

And finally we define reaction rates for binding of AG to the allosteric site when nothing, AA, or AG is bound at the catalytic site:

In [41]:
Parameter('kf_AG_allo1', 1000.0)
Out[41]:
Parameter('kf_AG_allo1', 1000.0)
In [42]:
Parameter('kr_AG_allo1', 1.0e5)
Out[42]:
Parameter('kr_AG_allo1', 100000.0)
In [43]:
Parameter('kf_AG_allo2', 1000.0)
Out[43]:
Parameter('kf_AG_allo2', 1000.0)
In [44]:
Parameter('kr_AG_allo2', 400.0)
Out[44]:
Parameter('kr_AG_allo2', 400.0)
In [45]:
Parameter('kf_AG_allo3', 1000.0)
Out[45]:
Parameter('kf_AG_allo3', 1000.0)
In [46]:
Parameter('kr_AG_allo3', 63000.0)
Out[46]:
Parameter('kr_AG_allo3', 63000.0)

We have now defined all the 'pieces' of CORM (its monomers, initial conditions, and rate parameters), but we haven't yet defined what interactions are allowed in our model.

First we will use the PySB catalyze macro to specify that AA can bind COX2 when nothing is bound in its allosteric site and be converted to PG with rates given by the first three rate parameters we defined above:

In [47]:
catalyze(COX2(allo=None), 'cat', AA(), 'b', PG(), [kf_AA_cat1, kr_AA_cat1, kcat_AA1])
Out[47]:
ComponentSet([
 Rule('bind_COX2_AA_to_COX2AA', COX2(allo=None, cat=None) + AA(b=None) <> COX2(allo=None, cat=1) % AA(b=1), kf_AA_cat1, kr_AA_cat1),
 Rule('catalyze_COX2AA_to_COX2_PG', COX2(allo=None, cat=1) % AA(b=1) >> COX2(allo=None, cat=None) + PG(), kcat_AA1),
 ])

The first argument to the catalyze macro is the enzyme (COX2 with nothing bound in the allosteric site), then the site on the enzyme for catalytic binding ('cat'), then the substrate (AA), then the binding site on the substrate ('b'), then the product (PG), and finally, a list of rate parameters (forward and reverse rates for binding of enzyme and substrate and kcat for substrate turnover). The catalyze macro creates two instances of the PySB Rule class each of which defines allowed reaction(s). Each Rule has a name (i.e. 'bind_COX2_AA_to_COX2AA' for the first rule), a defined interaction (i.e. COX2 and AA can reversibly bind when both have nothing previously bound), and rate parameters for the interaction (i.e. kf_AA_cat_1 and kr_AA_cat1).

While in the case of CORM we define rules that apply to only one set of reactants and products, this does not have to be so. The power of PySB and other 'rule-based' modeling languages is that rules can match many different reactions, meaning you do not have to explicitly define every allowed interaction in your model. For CORM we do choose to individually define these reactions so that we can assign different rate parameters to each reaction.

Next, we'll use the bind_complex macro and a PySB rule to state that COX2 with AG bound in the allosteric site can bind AA and create PG:

In [48]:
bind_complex(COX2(allo=1) % AG(b=1), 'cat', AA(), 'b', [kf_AA_cat2, kr_AA_cat2])
Out[48]:
ComponentSet([
 Rule('bind_COX2AG_AA', COX2(allo=1, cat=None) % AG(b=1) + AA(b=None) <> COX2(allo=1, cat=50) % AG(b=1) % AA(b=50), kf_AA_cat2, kr_AA_cat2),
 ])
In [49]:
Rule('kcat_AA_2',
     COX2(allo=1, cat=2) % AG(b=1) % AA(b=2) >> COX2(allo=1, cat=None) % AG(b=1) + PG(),
    kcat_AA2)
Out[49]:
Rule('kcat_AA_2', COX2(allo=1, cat=2) % AG(b=1) % AA(b=2) >> COX2(allo=1, cat=None) % AG(b=1) + PG(), kcat_AA2)

Then we'll do the same to allow COX2 binding and turnover of AA when AA is also bound in the allosteric site:

In [50]:
bind_complex(COX2(allo=1) % AA(b=1), 'cat', AA(), 'b', [kf_AA_cat3, kr_AA_cat3])
Out[50]:
ComponentSet([
 Rule('bind_COX2AA_AA', COX2(allo=1, cat=None) % AA(b=1) + AA(b=None) <> COX2(allo=1, cat=50) % AA(b=1) % AA(b=50), kf_AA_cat3, kr_AA_cat3),
 ])
In [51]:
Rule('kcat_AA_3',
     COX2(allo=1, cat=2) % AA(b=1) % AA(b=2) >> COX2(allo=1, cat=None) % AA(b=1) + PG(),
    kcat_AA3)
Out[51]:
Rule('kcat_AA_3', COX2(allo=1, cat=2) % AA(b=1) % AA(b=2) >> COX2(allo=1, cat=None) % AA(b=1) + PG(), kcat_AA3)

Next we create analagous rules for binding and turnover of AG when nothing, AG, or AA is bound in the allosteric site:

In [52]:
catalyze(COX2(allo=None), 'cat', AG(), 'b', PGG(), [kf_AG_cat1, kr_AG_cat1, kcat_AG1])
Out[52]:
ComponentSet([
 Rule('bind_COX2_AG_to_COX2AG', COX2(allo=None, cat=None) + AG(b=None) <> COX2(allo=None, cat=1) % AG(b=1), kf_AG_cat1, kr_AG_cat1),
 Rule('catalyze_COX2AG_to_COX2_PGG', COX2(allo=None, cat=1) % AG(b=1) >> COX2(allo=None, cat=None) + PGG(), kcat_AG1),
 ])
In [53]:
bind_complex(COX2(allo=1) % AG(b=1), 'cat', AG(), 'b', [kf_AG_cat2, kr_AG_cat2])
Out[53]:
ComponentSet([
 Rule('bind_COX2AG_AG', COX2(allo=1, cat=None) % AG(b=1) + AG(b=None) <> COX2(allo=1, cat=50) % AG(b=1) % AG(b=50), kf_AG_cat2, kr_AG_cat2),
 ])
In [54]:
bind_complex(COX2(allo=1) % AA(b=1), 'cat', AG(), 'b', [kf_AG_cat3, kr_AG_cat3])
Out[54]:
ComponentSet([
 Rule('bind_COX2AA_AG', COX2(allo=1, cat=None) % AA(b=1) + AG(b=None) <> COX2(allo=1, cat=50) % AA(b=1) % AG(b=50), kf_AG_cat3, kr_AG_cat3),
 ])
In [55]:
Rule('kcat_AG_3',
     COX2(allo=1, cat=2) % AA(b=1) % AG(b=2) >> COX2(allo=1, cat=None) % AA(b=1) + PGG(),
    kcat_AG3)
Out[55]:
Rule('kcat_AG_3', COX2(allo=1, cat=2) % AA(b=1) % AG(b=2) >> COX2(allo=1, cat=None) % AA(b=1) + PGG(), kcat_AG3)

The final set of set of rules we need to add will tell PySB that AA and AG can both bind to the allosteric site. First we'll use the bind macro to add the case where AA binds the allosteric site and nothing is in the catalytic site:

In [56]:
bind(COX2(cat=None), 'allo', AA(), 'b', [kf_AA_allo1, kr_AA_allo1])
Out[56]:
ComponentSet([
 Rule('bind_COX2_AA', COX2(allo=None, cat=None) + AA(b=None) <> COX2(allo=1, cat=None) % AA(b=1), kf_AA_allo1, kr_AA_allo1),
 ])

Then we'll use rules to allow the binding of AA to the allosteric site when AA or AG is bound to the catalytic site:

In [57]:
Rule('bind_COX2AA_AA_allo',
     COX2(cat=1, allo=None) % AA(b=1) + AA(b=None) <> COX2(cat=1, allo=2) % AA(b=1) % AA(b=2),
    kf_AA_allo2, kr_AA_allo2)
Out[57]:
Rule('bind_COX2AA_AA_allo', COX2(allo=None, cat=1) % AA(b=1) + AA(b=None) <> COX2(allo=2, cat=1) % AA(b=1) % AA(b=2), kf_AA_allo2, kr_AA_allo2)
In [58]:
Rule('bind_COX2AG_AA_allo',
     COX2(cat=1, allo=None) % AG(b=1) + AA(b=None) <> COX2(cat=1, allo=2) % AG(b=1) % AA(b=2),
    kf_AA_allo3, kr_AA_allo3)
Out[58]:
Rule('bind_COX2AG_AA_allo', COX2(allo=None, cat=1) % AG(b=1) + AA(b=None) <> COX2(allo=2, cat=1) % AG(b=1) % AA(b=2), kf_AA_allo3, kr_AA_allo3)

And finally we'll create the analogous rules for AG binding to the allosteric site with nothing, AA, or AG in the allosteric site:

In [59]:
bind(COX2(cat=None), 'allo', AG(), 'b', [kf_AG_allo1, kr_AG_allo1])
Out[59]:
ComponentSet([
 Rule('bind_COX2_AG', COX2(allo=None, cat=None) + AG(b=None) <> COX2(allo=1, cat=None) % AG(b=1), kf_AG_allo1, kr_AG_allo1),
 ])
In [60]:
Rule('bind_COX2AA_AG_allo',
     COX2(cat=1, allo=None) % AA(b=1) + AG(b=None) <> COX2(cat=1, allo=2) % AA(b=1) % AG(b=2),
    kf_AG_allo2, kr_AG_allo2)
Out[60]:
Rule('bind_COX2AA_AG_allo', COX2(allo=None, cat=1) % AA(b=1) + AG(b=None) <> COX2(allo=2, cat=1) % AA(b=1) % AG(b=2), kf_AG_allo2, kr_AG_allo2)
In [61]:
Rule('bind_COX2AG_AG_allo',
     COX2(cat=1, allo=None) % AG(b=1) + AG(b=None) <> COX2(cat=1, allo=2) % AG(b=1) % AG(b=2),
    kf_AG_allo3, kr_AG_allo3)
Out[61]:
Rule('bind_COX2AG_AG_allo', COX2(allo=None, cat=1) % AG(b=1) + AG(b=None) <> COX2(allo=2, cat=1) % AG(b=1) % AG(b=2), kf_AG_allo3, kr_AG_allo3)

At this point we've defined all the monomers, initial conditions, parameters, and the allowed interactions present in CORM. The model is complete and ready to be simulated.

Before simulating, we can 'tag' particular species we're interested in to make it easy to plot their trajectories after simulation. To do this we use PySB Observables:

In [62]:
Observable('obsPG', PG())
Out[62]:
Observable('obsPG', PG())
In [63]:
Observable('obsPGG', PGG())
Out[63]:
Observable('obsPGG', PGG())
In [64]:
Observable('obsAA', AA())
Out[64]:
Observable('obsAA', AA())
In [65]:
Observable('obsAG', AG())
Out[65]:
Observable('obsAG', AG())

Each observable has a name and a species it matches. The species pattern can match multiple chemical species; for instance, you could match all the COX2 with nothing bound at the allosteric site (but either AA, AG, or nothing bound at the catalytic site) using the pattern COX2(allo=None).

To simulate the model, we import the function PySB odesolve, which will solve the ordinary differential equations (ODEs) that PySB has created under the hood for CORM. We also import the Python package Numpy to allow us to define a trajectory of timepoints:

In [66]:
from pysb.integrate import odesolve
import numpy

Then we define t, the vector of time points we want to simulate at (here 100 evenly spaced time points from 0 to 10 seconds):

In [68]:
t = numpy.linspace(0, 10, num=100)

Finally, we simulate using odesolve and store the matrix of all model species at all requested time points in the variable yout:

In [69]:
yout = odesolve(model, t)
creating /var/folders/z3/m70qnpgj7gx1_bj92txc141h0000gn/T/scipy-Erin-NsKoX4/python27_intermediate/compiler_10a3badd5737329477237726cc2e9ec1

We can pull out the concentrations for a any of the observables we defined using the observable name:

In [70]:
yout['obsPG']
Out[70]:
array([ 0.        ,  0.00140281,  0.00270124,  0.00397287,  0.00523782,
        0.0065011 ,  0.00776396,  0.00902671,  0.01028943,  0.01155214,
        0.01281483,  0.01407753,  0.01534021,  0.0166029 ,  0.01786557,
        0.01912824,  0.02039091,  0.02165357,  0.02291623,  0.02417888,
        0.02544152,  0.02670416,  0.02796679,  0.02922942,  0.03049204,
        0.03175466,  0.03301727,  0.03427988,  0.03554248,  0.03680507,
        0.03806766,  0.03933025,  0.04059283,  0.0418554 ,  0.04311797,
        0.04438053,  0.04564309,  0.04690564,  0.04816818,  0.04943072,
        0.05069326,  0.05195579,  0.05321831,  0.05448083,  0.05574335,
        0.05700585,  0.05826836,  0.05953085,  0.06079334,  0.06205583,
        0.06331831,  0.06458079,  0.06584326,  0.06710572,  0.06836818,
        0.06963063,  0.07089308,  0.07215552,  0.07341796,  0.07468039,
        0.07594282,  0.07720524,  0.07846765,  0.07973006,  0.08099247,
        0.08225487,  0.08351726,  0.08477965,  0.08604203,  0.08730441,
        0.08856678,  0.08982914,  0.0910915 ,  0.09235386,  0.09361621,
        0.09487855,  0.09614089,  0.09740322,  0.09866555,  0.09992787,
        0.10119019,  0.1024525 ,  0.10371481,  0.10497711,  0.1062394 ,
        0.10750169,  0.10876397,  0.11002625,  0.11128853,  0.11255079,
        0.11381306,  0.11507531,  0.11633756,  0.11759981,  0.11886205,
        0.12012428,  0.12138651,  0.12264874,  0.12391096,  0.12517317])

To create plots of our simulated data, we can use the Python package matplotlib:

In [82]:
from matplotlib import pyplot
%matplotlib inline
In [87]:
pyplot.figure()
pyplot.plot(t, yout['obsPG'], label='PG')
pyplot.plot(t, yout['obsPGG'], label='PGG')
pyplot.legend()
Out[87]:
<matplotlib.legend.Legend at 0x1109eb390>

Here, for instance, we see the concentrations of the products (PG and PGG) over time.

This concludes the CORM tutorial. To learn more about other features of PySB models, see the PySB documentation. To see how the CORM model parameters were fitted to experimental data, see the file run_pymc_sampling.py.