In this notebook we explore the effects of specified conditions on Agents (e.g. bound conditions, modification conditions) and assembly policies on the combinatorial complexity of dynamical models.
First, we import INDRA's TRIPS input API and PySB model assembler.
In [1]:
from indra.sources import trips
from indra.assemblers.pysb import PysbAssembler
# Below is some bookkeeping code needed to display reaction network graphs
%matplotlib inline
from pysb.tools import render_reactions
import pygraphviz, subprocess
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
def draw_reaction_network(model):
pygraphviz.AGraph(render_reactions.run(model)).draw('model.png', prog='dot')
img = mpimg.imread('model.png')
plt.figure(figsize=(50, 50))
plt.imshow(img)
plt.xticks([])
plt.yticks([])
In [2]:
tp = trips.process_text('RAF binds Vemurafenib, RAF phosphorylates MEK and MEK phosphorylates ERK.')
This yields 3 INDRA Statements, as follows. Here empty parentheses after the Agent names indicate that there is no additional context specified on them.
In [3]:
tp.statements
Out[3]:
In [4]:
pa = PysbAssembler()
pa.add_statements(tp.statements)
pa.make_model(policies='one_step')
Out[4]:
In [5]:
model1_one = pa.model
The model has 4 Monomers and 4 Rules.
In [6]:
model1_one.monomers
Out[6]:
In [7]:
model1_one.rules
Out[7]:
Let's examine the last rule which corresponds to MEK phosphorylating ERK. Here, MEK() appears without additional context specified. This means that the rule will apply to any form of MEK, for instance, MEK that is unphosphorylated.
We now generate the rule-based model into a reaction network form using PySB's interface to BioNetGen.
In [8]:
from pysb.bng import generate_equations
generate_equations(model1_one)
We can now plot the reaction network to examine the model. Each colored node of the reaction network is a molecular species, reactions are represented by gray nodes, and arrows show species being consumed and produced by each reaction.
In [9]:
draw_reaction_network(model1_one)
We see from the reaction network that RAF is able to phosphorylate MEK whether or not it is bound to Vemurafenib, and MEK phosphorylates ERK whether or not it is phosphorylated.
In [10]:
pa.make_model(policies='two_step')
Out[10]:
In [11]:
model1_two = pa.model
In [12]:
model1_two.monomers
Out[12]:
In [13]:
for rule in model1_two.rules:
print(rule.rule_expression)
We can now generate the reaction network for model1_two and inspect the reaction network that is created.
In [14]:
generate_equations(model1_two)
In [15]:
model1_two.species
Out[15]:
In [16]:
draw_reaction_network(model1_two)
The two-step policy produced a total of 13 molecular species and 19 reactions. ERK now appears in 6 possible forms:
This means that our initial description allowed for the possibility of ERK, MEK, RAF and Vemurafenib all simultaneously being in a complex. While the existence of such a complex is not impossible, we can introduce additional assumptions to refine the model.
In [17]:
tp = trips.process_text('RAF binds Vemurafenib. '
'RAF not bound to Vemurafenib phosphorylates MEK. '
'Phosphorylated MEK not bound to RAF phosphorylates ERK.')
The INDRA Statements extracted by processing this text are shown below.
In [18]:
tp.statements
Out[18]:
We see that some agents are now subject to additional conditions, for instance, RAF(bound: [VEMURAFENIB, False]) specifies that RAF should not be bound to Vemurafenib in order to phosphorylate MEK.
Let's now assemble the model and generate the reaction network.
In [19]:
pa = PysbAssembler()
pa.add_statements(tp.statements)
pa.make_model(policies='two_step')
Out[19]:
In [20]:
generate_equations(pa.model)
In [21]:
draw_reaction_network(pa.model)
The model is now significantly simpler with a total of 7 reactions down from 19 in the previous model.
As an alternative to the two-step policy, we can assemble the same model using a simpler, one-step policy.
In [22]:
pa.make_model(policies='one_step')
Out[22]:
In [23]:
generate_equations(pa.model)
In [24]:
draw_reaction_network(pa.model)
As we see, this model contains 7 individual species with 3 reactions in total.
One issue with the simple one-step policy is that enzymatic catalysis is modeled with pseudo-first order kinetics. Alternatively, we can use a Michaelis-Menten policy in with case both phosphorylation processes are still effectively modeled as one-step but their kinetic rates will account for enzyme saturation. As seen below, the model still retains the same structure as under the one-step policy, only kinetic rates change.
In [25]:
pa.make_model(policies='michaelis_menten')
Out[25]:
In [26]:
draw_reaction_network(pa.model)