"""
A collection of generally useful modeling macros.
These macros are written to be as generic and reusable as possible, serving as a
collection of best practices and implementation ideas. They conform to the
following general guidelines:
* All components created by the macro are implicitly added to the current model
and explicitly returned in a ComponentSet.
* Parameters may be passed as Parameter or Expression objects, or as plain
numbers for which Parameter objects will be automatically created using an
appropriate naming convention.
* Arguments which accept a MonomerPattern should also accept Monomers, which are
to be interpreted as MonomerPatterns on that Monomer with an empty condition
list. This is typically implemented by having the macro apply the "call"
(parentheses) operator to the argument with an empty argument list and using
the resulting value instead of the original argument when creating Rules, e.g.
``arg = arg()``. Calling a Monomer will return a MonomerPattern, and calling a
MonomerPattern will return a copy of itself, so calling either is guaranteed
to return a MonomerPattern.
The _macro_rule helper function contains much of the logic needed to follow
these guidelines. Every macro in this module either uses _macro_rule directly or
calls another macro which does.
Another useful function is _verify_sites which will raise an exception if a
Monomer or MonomerPattern does not possess every one of a given list of sites.
This can be used to trigger such errors up front rather than letting an
exception occur at the point where the macro tries to use the invalid site in a
pattern, which can be harder for the caller to debug.
"""
from pysb import *
from pysb.core import ComponentSet, as_complex_pattern, MonomerPattern, ComplexPattern
import numbers
import functools
import itertools
from sympy import Piecewise
__all__ = ['equilibrate',
'bind', 'bind_table',
'catalyze', 'catalyze_state', 'catalyze_complex',
'catalyze_one_step', 'catalyze_one_step_reversible',
'synthesize', 'degrade', 'synthesize_degrade_table',
'assemble_pore_sequential', 'pore_transport', 'pore_bind', 'assemble_chain_sequential_base',
'bind_complex', 'bind_table_complex', 'drug_binding']
# Suppress ModelExistsWarnings in our doctests.
_pysb_doctest_suppress_modelexistswarning = True
# Internal helper functions
# =========================
def _complex_pattern_label(cp):
"""Return a string label for a ComplexPattern."""
if cp is None:
return ''
mp_labels = [_monomer_pattern_label(mp) for mp in cp.monomer_patterns]
return ''.join(mp_labels)
def _monomer_pattern_label(mp):
"""Return a string label for a MonomerPattern."""
site_values = [str(x) for x in mp.site_conditions.values()
if x is not None
and not isinstance(x, list)
and not isinstance(x, tuple)
and not isinstance(x, numbers.Real)]
return mp.monomer.name + ''.join(site_values)
def _rule_name_generic(rule_expression):
"""Return a generic string label for a RuleExpression."""
# Get ReactionPatterns
react_p = rule_expression.reactant_pattern
prod_p = rule_expression.product_pattern
# Build the label components
lhs_label = [_complex_pattern_label(cp) for cp in react_p.complex_patterns]
lhs_label = '_'.join(lhs_label)
rhs_label = [_complex_pattern_label(cp) for cp in prod_p.complex_patterns]
rhs_label = '_'.join(rhs_label)
return '%s_to_%s' % (lhs_label, rhs_label)
def _macro_rule(rule_prefix, rule_expression, klist, ksuffixes,
name_func=_rule_name_generic):
"""
A helper function for writing macros that generates a single rule.
Parameters
----------
rule_prefix : string
The prefix that is prepended to the (automatically generated) name for
the rule.
rule_expression : RuleExpression
An expression specifying the form of the rule; gets passed directly
to the Rule constructor.
klist : list of Parameters or Expressions, or list of numbers
If the rule is unidirectional, the list must contain one element
(either a Parameter/Expression or number); if the rule is reversible,
it must contain two elements. If the rule is reversible, the first
element in the list is taken to be the forward rate, and the second
element is taken as the reverse rate.
ksuffixes : list of strings
If klist contains numbers rather than Parameters or Expressions, the
strings in ksuffixes are used to automatically generate the necessary
Parameter objects. The suffixes are appended to the rule name to
generate the associated parameter name. ksuffixes must contain one
element if the rule is unidirectional, two if it is reversible.
name_func : function, optional
A function which takes a RuleExpression and returns a string label for
it, to be called as part of the automatic rule name generation. If not
provided, a built-in default naming function will be used.
Returns
-------
components : ComponentSet
The generated components. Contains the generated Rule and up to two
generated Parameter objects (if klist was given as numbers).
Notes
-----
The default naming scheme (if `name_func` is not passed) follows the form::
'%s_%s_to_%s' % (rule_prefix, lhs_label, rhs_label)
where lhs_label and rhs_label are each concatenations of the Monomer names
and specified sites in the ComplexPatterns on each side of the
RuleExpression. The actual implementation is in the function
_rule_name_generic, which in turn calls _complex_pattern_label and
_monomer_pattern_label. For some specialized reactions it may be helpful to
devise a custom naming scheme rather than rely on this default.
Examples
--------
Using distinct Monomers for substrate and product::
>>> from pysb import *
>>> from pysb.macros import _macro_rule
>>>
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A', ['s'])
Monomer('A', ['s'])
>>> Monomer('B', ['s'])
Monomer('B', ['s'])
>>>
>>> _macro_rule('bind', A(s=None) + B(s=None) | A(s=1) % B(s=1),
... [1e6, 1e-1], ['kf', 'kr']) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_A_B_to_AB', A(s=None) + B(s=None) | A(s=1) % B(s=1),
bind_A_B_to_AB_kf, bind_A_B_to_AB_kr),
Parameter('bind_A_B_to_AB_kf', 1000000.0),
Parameter('bind_A_B_to_AB_kr', 0.1),
])
"""
r_name = '%s_%s' % (rule_prefix, name_func(rule_expression))
# If rule is unidirectional, make sure we only have one parameter
if (not rule_expression.is_reversible):
if len(klist) != 1 or len(ksuffixes) != 1:
raise ValueError("A unidirectional rule must have one parameter.")
# If rule is bidirectional, make sure we have two parameters
else:
if len(klist) != 2 or len(ksuffixes) != 2:
raise ValueError("A bidirectional rule must have two parameters.")
if all(isinstance(x, (Parameter, Expression)) for x in klist):
k1 = klist[0]
if rule_expression.is_reversible:
k2 = klist[1]
params_created = ComponentSet()
# if klist is numbers, generate the Parameters
elif all(isinstance(x, numbers.Real) for x in klist):
k1 = Parameter('%s_%s' % (r_name, ksuffixes[0]), klist[0])
params_created = ComponentSet([k1])
if rule_expression.is_reversible:
k2 = Parameter('%s_%s' % (r_name, ksuffixes[1]),
klist[1])
params_created.add(k2)
else:
raise ValueError("klist must contain Parameters, Expressions, or numbers.")
if rule_expression.is_reversible:
r = Rule(r_name, rule_expression, k1, k2)
else:
r = Rule(r_name, rule_expression, k1)
# Build a set of components that were created
return ComponentSet([r]) | params_created
def _verify_sites(m, *site_list):
"""
Checks that the monomer m contains all of the sites in site_list.
Parameters
----------
m : Monomer or MonomerPattern
The monomer to check.
site1, site2, ... : string
One or more site names to check on m
Returns
-------
True if m contains all sites; raises a ValueError otherwise.
Raises
------
ValueError
If any of the sites are not found.
"""
if isinstance(m, ComplexPattern):
return _verify_sites_complex(m, *site_list)
else:
for site in site_list:
if site not in m().monomer.sites:
raise ValueError("Monomer '%s' must contain the site '%s'" %
(m().monomer.name, site))
return True
def _verify_sites_complex(c, *site_list):
"""
Checks that the complex c contains all of the sites in site_list.
Parameters
----------
c : ComplexPattern
The complex to check.
site1, site2, ... : string
One or more site names to check on c
Returns
-------
If all sites are found within the complex, a dictionary of monomers and the sites within site_list they contain. Raises a ValueError if one or more sites not in the complex.
Raises
------
ValueError
If any of the sites are not found within the complex.
"""
allsitesdict = {}
for mon in c.monomer_patterns:
allsitesdict[mon] = mon.monomer.sites
for site in site_list:
specsitesdict = {}
for monomer, li in allsitesdict.items():
for s in li:
if site in li:
specsitesdict[monomer] = site
if len(specsitesdict) == 0:
raise ValueError("Site '%s' not found in complex '%s'" % (site, c))
return specsitesdict
# Unimolecular patterns
# =====================
[docs]def equilibrate(s1, s2, klist):
"""
Generate the unimolecular reversible equilibrium reaction S1 <-> S2.
Parameters
----------
s1, s2 : Monomer or MonomerPattern
S1 and S2 in the above reaction.
klist : list of 2 Parameters or list of 2 numbers
Forward (S1 -> S2) and reverse rate constants (in that order). If
Parameters are passed, they will be used directly in the generated
Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the names and states of S1 and S2
and these parameters will be included at the end of the returned
component list.
Returns
-------
components : ComponentSet
The generated components. Contains one reversible Rule and optionally
two Parameters if klist was given as plain numbers.
Examples
--------
Simple two-state equilibrium between A and B::
Model()
Monomer('A')
Monomer('B')
equilibrate(A(), B(), [1, 1])
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A')
Monomer('A')
>>> Monomer('B')
Monomer('B')
>>> equilibrate(A(), B(), [1, 1]) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('equilibrate_A_to_B', A() | B(), equilibrate_A_to_B_kf, equilibrate_A_to_B_kr),
Parameter('equilibrate_A_to_B_kf', 1.0),
Parameter('equilibrate_A_to_B_kr', 1.0),
])
"""
# turn any Monomers into MonomerPatterns
return _macro_rule('equilibrate', s1 | s2, klist, ['kf', 'kr'])
# Binding
# =======
[docs]def bind(s1, site1, s2, site2, klist):
"""
Generate the reversible binding reaction S1 + S2 | S1:S2.
Parameters
----------
s1, s2 : Monomer or MonomerPattern
Monomers participating in the binding reaction.
site1, site2 : string
The names of the sites on s1 and s2 used for binding.
klist : list of 2 Parameters or list of 2 numbers
Forward and reverse rate constants (in that order). If Parameters are
passed, they will be used directly in the generated Rules. If numbers
are passed, Parameters will be created with automatically generated
names based on the names and states of S1 and S2 and these parameters
will be included at the end of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains the bidirectional binding Rule
and optionally two Parameters if klist was given as numbers.
Examples
--------
Binding between A and B::
Model()
Monomer('A', ['x'])
Monomer('B', ['y'])
bind(A, 'x', B, 'y', [1e-4, 1e-1])
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A', ['x'])
Monomer('A', ['x'])
>>> Monomer('B', ['y'])
Monomer('B', ['y'])
>>> bind(A, 'x', B, 'y', [1e-4, 1e-1]) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_A_B', A(x=None) + B(y=None) | A(x=1) % B(y=1), bind_A_B_kf, bind_A_B_kr),
Parameter('bind_A_B_kf', 0.0001),
Parameter('bind_A_B_kr', 0.1),
])
"""
_verify_sites(s1, site1)
_verify_sites(s2, site2)
return _macro_rule('bind',
s1(**{site1: None}) + s2(**{site2: None}) |
s1(**{site1: 1}) % s2(**{site2: 1}),
klist, ['kf', 'kr'], name_func=bind_name_func)
def bind_name_func(rule_expression):
# Get ComplexPatterns
react_cps = rule_expression.reactant_pattern.complex_patterns
# Build the label components
return '_'.join(_complex_pattern_label(cp) for cp in react_cps)
[docs]def bind_complex(s1, site1, s2, site2, klist, m1=None, m2=None):
"""
Generate the reversible binding reaction ``S1 + S2 | S1:S2``,
with optional complexes attached to either
``S1`` (``C1:S1 + S2 | C1:S1:S2``), ``S2`` (``S1 + C2:S2 | C2:S2:S1``),
or both (``C1:S1 + C2:S2 | C1:S1:S2:C2``).
Parameters
----------
s1, s2 : Monomer, MonomerPattern, or ComplexPattern
Monomers or complexes participating in the binding reaction.
site1, site2 : string
The names of the sites on s1 and s2 used for binding.
klist : list of 2 Parameters or list of 2 numbers
Forward and reverse rate constants (in that order). If Parameters are
passed, they will be used directly in the generated Rules. If numbers
are passed, Parameters will be created with automatically generated
names based on the names and states of S1 and S2 and these parameters
will be included at the end of the returned component list.
m1, m2 : Monomer or MonomerPattern
If s1 or s2 binding site is present in multiple monomers
within a complex, the specific monomer desired for binding must be specified.
Returns
-------
components : ComponentSet
The generated components. Contains the bidirectional binding Rule
and optionally two Parameters if klist was given as numbers.
Examples
--------
Binding between ``A:B`` and ``C:D``:
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A', ['a', 'b'])
Monomer('A', ['a', 'b'])
>>> Monomer('B', ['c', 'd'])
Monomer('B', ['c', 'd'])
>>> Monomer('C', ['e', 'f'])
Monomer('C', ['e', 'f'])
>>> Monomer('D', ['g', 'h'])
Monomer('D', ['g', 'h'])
>>> bind_complex(A(a=1) % B(c=1), 'b', C(e=2) % D(g=2), 'h', [1e-4, \
1e-1]) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_AB_DC', A(a=1, b=None) % B(c=1) + D(g=3, h=None) % C(e=3)
| A(a=1, b=50) % B(c=1) % D(g=3, h=50) % C(e=3), bind_AB_DC_kf,
bind_AB_DC_kr),
Parameter('bind_AB_DC_kf', 0.0001),
Parameter('bind_AB_DC_kr', 0.1),
])
Execution:
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A', ['a', 'b'])
Monomer('A', ['a', 'b'])
>>> Monomer('B', ['c', 'd'])
Monomer('B', ['c', 'd'])
>>> Monomer('C', ['e', 'f'])
Monomer('C', ['e', 'f'])
>>> Monomer('D', ['g', 'h'])
Monomer('D', ['g', 'h'])
>>> bind(A, 'a', B, 'c', [1e4, 1e-1]) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_A_B',
A(a=None) + B(c=None) | A(a=1) % B(c=1),
bind_A_B_kf, bind_A_B_kr),
Parameter('bind_A_B_kf', 10000.0),
Parameter('bind_A_B_kr', 0.1),
])
>>> bind(C, 'e', D, 'g', [1e4, 1e-1]) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_C_D',
C(e=None) + D(g=None) | C(e=1) % D(g=1),
bind_C_D_kf, bind_C_D_kr),
Parameter('bind_C_D_kf', 10000.0),
Parameter('bind_C_D_kr', 0.1),
])
>>> bind_complex(A(a=1) % B(c=1), 'b', C(e=2) % D(g=2), 'h', [1e-4, \
1e-1]) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_AB_DC',
A(a=1, b=None) % B(c=1) + D(g=3, h=None) % C(e=3) | A(a=1,
b=50) % B(c=1) % D(g=3, h=50) % C(e=3),
bind_AB_DC_kf, bind_AB_DC_kr),
Parameter('bind_AB_DC_kf', 0.0001),
Parameter('bind_AB_DC_kr', 0.1),
])
"""
if isinstance(m1, Monomer):
m1 = m1()
if isinstance(m2, Monomer):
m2 = m2()
#Define some functions for checking complex sites, building complexes up from monomers, and creating rules.
def comp_mono_func(s1, site1, s2, site2, m1):
_verify_sites(s2, site2)
#Retrieve a dictionary specifying the MonomerPattern within the complex that contains the given binding site.
specsites = list(_verify_sites_complex(s1, site1))
s1complexpatub, s1complexpatb = check_sites_comp_build(s1, site1, m1, specsites)
return create_rule(s1complexpatub, s1complexpatb, s2({site2:None}), s2({site2: 50}))
def check_sites_comp_build(s1, site1, m1, specsites):
#Return error if binding site exists on multiple monomers and a monomer for binding (m1) hasn't been specified.
if len(specsites) > 1 and m1==None:
raise ValueError("Binding site '%s' present in more than one monomer in complex '%s'. Specify variable m1, the monomer used for binding within the complex." % (site1, s1))
if not s1.is_concrete:
raise ValueError("Complex '%s' must be concrete." % (s1))
#If the given binding site is only present in one monomer in the complex:
if m1==None:
#Build up ComplexPattern for use in rule (with state of given binding site specified).
s1complexpatub = specsites[0]({site1:None})
s1complexpatb = specsites[0]({site1:50})
for monomer in s1.monomer_patterns:
if monomer not in specsites:
s1complexpatub %= monomer
s1complexpatb %= monomer
#If the binding site is present on more than one monomer in the complex, the monomer must be specified by the user. Use specified m1 to build ComplexPattern.
else:
#Make sure binding states of MonomerPattern m1 match those of the monomer within the ComplexPattern s1 (ComplexPattern monomer takes precedence if not).
i = 0
identical_monomers = []
other_monomers = []
for mon in s1.monomer_patterns:
#Only change the binding site for the first monomer that matches. Keep any others unchanged to add to final complex that is returned.
if mon.monomer.name == m1.monomer.name and mon.site_conditions==m1.site_conditions:
i += 1
if i == 1:
s1complexpatub = mon({site1:None})
s1complexpatb = mon({site1:50})
else:
identical_monomers.append(mon)
else:
other_monomers.append(mon)
#Throw an error if no monomer pattern in the complex matched the pattern given for m1
if i == 0:
raise ValueError("No monomer pattern in complex '%s' matches the pattern given for m1, '%s'." % (s1, m1))
#Build up ComplexPattern for use in rule (with state of given binding site on m1 specified).
for mon in other_monomers:
s1complexpatub %= mon
s1complexpatb %= mon
if identical_monomers:
for i in range(len(identical_monomers)):
s1complexpatub %= identical_monomers[i]
s1complexpatb %= identical_monomers[i]
return s1complexpatub, s1complexpatb
#Create rules.
def create_rule(s1ub, s1b, s2ub, s2b):
return _macro_rule('bind',
s1ub + s2ub |
s1b % s2b,
klist, ['kf', 'kr'], name_func=bind_name_func)
#If no complexes given, revert to normal bind macro.
if (isinstance(s1, MonomerPattern) or isinstance(s1, Monomer)) and (isinstance(s2, MonomerPattern) or isinstance(s2, Monomer)):
_verify_sites(s1, site1)
_verify_sites(s2, site2)
return bind(s1, site1, s2, site2, klist)
#Create rules if only one complex or the other is present.
elif isinstance(s1, ComplexPattern) and (isinstance(s2, MonomerPattern) or isinstance(s2, Monomer)):
return comp_mono_func(s1, site1, s2, site2, m1)
elif (isinstance(s1, MonomerPattern) or isinstance(s1, Monomer)) and isinstance(s2, ComplexPattern):
return comp_mono_func(s2, site2, s1, site1, m2)
#Create rule when both s1 and s2 are complexes.
else:
#Retrieve a dictionary specifiying the MonomerPattern within
#the complex that contains the given binding site. Convert to list.
specsites1 = list(_verify_sites_complex(s1, site1))
specsites2 = list(_verify_sites_complex(s2, site2))
#Return error if binding site exists on multiple monomers and a monomer for binding (m1/m2) hasn't been specified.
if len(specsites1) > 1 and m1==None:
raise ValueError("Binding site '%s' present in more than one monomer in complex '%s'. Specify variable m1, the monomer used for binding within the complex." % (site1, s1))
if len(specsites2) > 1 and m2==None:
raise ValueError("Binding site '%s' present in more than one monomer in complex '%s'. Specify variable m2, the monomer used for binding within the complex." % (site2, s2))
if not s1.is_concrete:
raise ValueError("Complex '%s' must be concrete." % (s1))
if not s2.is_concrete:
raise ValueError("Complex '%s' must be concrete." % (s2))
#To avoid creating rules with multiple bonds to the same site when combining the two complexes, check for the maximum bond integer in s1 and add to all s2 bond integers.
maxint = 0
for monomer in s1.monomer_patterns:
for stateint in monomer.site_conditions.values():
if isinstance(stateint, int):
if stateint > maxint:
maxint = stateint
match = 'N'
for monomer in s2.monomer_patterns:
if m2 is not None:
if m2.site_conditions == monomer.site_conditions and m2.monomer.name == monomer.monomer.name:
match = 'Y'
for site, stateint in monomer.site_conditions.items():
if isinstance(stateint, int):
monomer.site_conditions[site] += maxint
if match == 'Y':
m2.site_conditions = monomer.site_conditions
match = 'N'
#Actually create rules
s1complexpatub, s1complexpatb = check_sites_comp_build(s1, site1, m1, specsites1)
s2complexpatub, s2complexpatb = check_sites_comp_build(s2, site2, m2, specsites2)
return create_rule(s1complexpatub, s1complexpatb, s2complexpatub, s2complexpatb)
[docs]def bind_table(bindtable, row_site, col_site, kf=None):
"""
Generate a table of reversible binding reactions.
Given two lists of species R and C, calls the `bind` macro on each pairwise
combination (R[i], C[j]). The species lists and the parameter values are
passed as a list of lists (i.e. a table) with elements of R passed as the
"row headers", elements of C as the "column headers", and forward / reverse
rate pairs (in that order) as tuples in the "cells". For example with two
elements in each of R and C, the table would appear as follows (note that
the first row has one fewer element than the subsequent rows)::
[[ C1, C2],
[R1, (1e-4, 1e-1), (2e-4, 2e-1)],
[R2, (3e-4, 3e-1), (4e-4, 4e-1)]]
Each parameter tuple may contain Parameters or numbers. If Parameters are
passed, they will be used directly in the generated Rules. If numbers are
passed, Parameters will be created with automatically generated names based
on the names and states of the relevant species and these parameters will be
included at the end of the returned component list. To omit any individual
reaction, pass None in place of the corresponding parameter tuple.
Alternately, single kd values (dissociation constant, kr/kf) may be
specified instead of (kf, kr) tuples. If kds are used, a single shared kf
Parameter or number must be passed as an extra `kf` argument. kr values for
each binding reaction will be calculated as kd*kf. It is important to
remember that the forward rate constant is a single parameter shared across
the entire bind table, as this may have implications for parameter fitting.
Parameters
----------
bindtable : list of lists
Table of reactants and rates, as described above.
row_site, col_site : string
The names of the sites on the elements of R and C, respectively, used
for binding.
kf : Parameter or number, optional
If the "cells" in bindtable are given as single kd values, this is the
shared kf used to calculate the kr values.
Returns
-------
components : ComponentSet
The generated components. Contains the bidirectional binding Rules and
optionally the Parameters for any parameters given as numbers.
Examples
--------
Binding table for two species types (R and C), each with two members::
Model()
Monomer('R1', ['x'])
Monomer('R2', ['x'])
Monomer('C1', ['y'])
Monomer('C2', ['y'])
bind_table([[ C1, C2],
[R1, (1e-4, 1e-1), (2e-4, 2e-1)],
[R2, (3e-4, 3e-1), None]],
'x', 'y')
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('R1', ['x'])
Monomer('R1', ['x'])
>>> Monomer('R2', ['x'])
Monomer('R2', ['x'])
>>> Monomer('C1', ['y'])
Monomer('C1', ['y'])
>>> Monomer('C2', ['y'])
Monomer('C2', ['y'])
>>> bind_table([[ C1, C2],
... [R1, (1e-4, 1e-1), (2e-4, 2e-1)],
... [R2, (3e-4, 3e-1), None]],
... 'x', 'y') # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_R1_C1', R1(x=None) + C1(y=None) | R1(x=1) % C1(y=1),
bind_R1_C1_kf, bind_R1_C1_kr),
Parameter('bind_R1_C1_kf', 0.0001),
Parameter('bind_R1_C1_kr', 0.1),
Rule('bind_R1_C2', R1(x=None) + C2(y=None) | R1(x=1) % C2(y=1),
bind_R1_C2_kf, bind_R1_C2_kr),
Parameter('bind_R1_C2_kf', 0.0002),
Parameter('bind_R1_C2_kr', 0.2),
Rule('bind_R2_C1', R2(x=None) + C1(y=None) | R2(x=1) % C1(y=1),
bind_R2_C1_kf, bind_R2_C1_kr),
Parameter('bind_R2_C1_kf', 0.0003),
Parameter('bind_R2_C1_kr', 0.3),
])
"""
# extract species lists and matrix of rates
s_rows = [row[0] for row in bindtable[1:]]
s_cols = bindtable[0]
kmatrix = [row[1:] for row in bindtable[1:]]
# ensure kf is passed when necessary
kiter = itertools.chain.from_iterable(kmatrix)
if any(isinstance(x, numbers.Real) for x in kiter) and kf is None:
raise ValueError("must specify kf when using single kd values")
# loop over interactions
components = ComponentSet()
for r, s_row in enumerate(s_rows):
for c, s_col in enumerate(s_cols):
klist = kmatrix[r][c]
if klist is not None:
# if user gave a single kd, calculate kr
if isinstance(klist, numbers.Real):
kd = klist
klist = (kf, kd*kf)
components |= bind(s_row(), row_site, s_col(), col_site, klist)
return components
[docs]def bind_table_complex(bindtable, row_site, col_site, m1=None, m2=None, kf=None):
"""
Generate a table of reversible binding reactions when either the row or column species (or both) have a complex bound to them.
Given two lists of species R and C (which can be complexes or monomers),
calls the `bind_complex` macro on each pairwise
combination (R[i], C[j]). The species lists and the parameter values are
passed as a list of lists (i.e. a table) with elements of R passed as the
"row headers", elements of C as the "column headers", and forward / reverse
rate pairs (in that order) as tuples in the "cells". For example with two
elements in each of R and C, the table would appear as follows (note that
the first row has one fewer element than the subsequent rows)::
[[ C1, C2],
[R1, (1e-4, 1e-1), (2e-4, 2e-1)],
[R2, (3e-4, 3e-1), (4e-4, 4e-1)]]
Each parameter tuple may contain Parameters or numbers. If Parameters are
passed, they will be used directly in the generated Rules. If numbers are
passed, Parameters will be created with automatically generated names based
on the names and states of the relevant species and these parameters will be
included at the end of the returned component list. To omit any individual
reaction, pass None in place of the corresponding parameter tuple.
Alternately, single kd values (dissociation constant, kr/kf) may be
specified instead of (kf, kr) tuples. If kds are used, a single shared kf
Parameter or number must be passed as an extra `kf` argument. kr values for
each binding reaction will be calculated as kd*kf. It is important to
remember that the forward rate constant is a single parameter shared across
the entire bind table, as this may have implications for parameter fitting.
Parameters
----------
bindtable : list of lists
Table of reactants and rates, as described above.
row_site, col_site : string
The names of the sites on the elements of R and C, respectively, used
for binding.
m1 : Monomer or MonomerPattern, optional
Monomer in row complex for binding. Must be specified if there are multiple monomers
that have the row_site within a complex.
m2 : Monomer or MonomerPattern, optional
Monomer in column complex for binding. Must be specified if there are multiple monomers
that have the col_site within a complex.
kf : Parameter or number, optional
If the "cells" in bindtable are given as single kd values, this is the
shared kf used to calculate the kr values.
Returns
-------
components : ComponentSet
The generated components. Contains the bidirectional binding Rules and
optionally the Parameters for any parameters given as numbers.
Examples
--------
Binding table for two species types (R and C, which can be complexes or monomers)::
Model()
Monomer('R1', ['x', 'c1'])
Monomer('R2', ['x', 'c1'])
Monomer('C1', ['y', 'c2'])
Monomer('C2', ['y', 'c2'])
bind(C1(y=None), 'c2', C1(y=None), 'c2', (1e-3, 1e-2))
bind(R1(x=None), 'c1', R2(x=None), 'c1', (1e-3, 1e-2))
bind_table_complex([[ C1(c2=1, y=None)%C1(c2=1), C2],
[R1()%R2(), (1e-4, 1e-1), (2e-4, 2e-1)],
[R2, (3e-4, 3e-1), None]],
'x', 'y', m1=R1(), m2=C1(y=None, c2=1))
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('R1', ['x', 'c1'])
Monomer('R1', ['x', 'c1'])
>>> Monomer('R2', ['x', 'c1'])
Monomer('R2', ['x', 'c1'])
>>> Monomer('C1', ['y', 'c2'])
Monomer('C1', ['y', 'c2'])
>>> Monomer('C2', ['y', 'c2'])
Monomer('C2', ['y', 'c2'])
>>> bind(C1(y=None), 'c2', C1(y=None), 'c2', (1e-3, 1e-2)) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_C1_C1', C1(y=None, c2=None) + C1(y=None, c2=None) | C1(y=None, c2=1) % C1(y=None, c2=1), bind_C1_C1_kf, bind_C1_C1_kr),
Parameter('bind_C1_C1_kf', 0.001),
Parameter('bind_C1_C1_kr', 0.01),
])
>>> bind(R1(x=None), 'c1', R2(x=None), 'c1', (1e-3, 1e-2)) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_R1_R2', R1(x=None, c1=None) + R2(x=None, c1=None) | R1(x=None, c1=1) % R2(x=None, c1=1), bind_R1_R2_kf, bind_R1_R2_kr),
Parameter('bind_R1_R2_kf', 0.001),
Parameter('bind_R1_R2_kr', 0.01),
])
>>> bind_table_complex([[ C1(c2=1, y=None)%C1(c2=1), C2],
... [R1()%R2(), (1e-4, 1e-1), (2e-4, 2e-1)],
... [R2, (3e-4, 3e-1), None]],
... 'x', 'y', m1=R1(), m2=C1(y=None, c2=1)) #doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_R1R2_C1C1', R1(x=None) % R2() + C1(y=None, c2=1) % C1(c2=1) | R1(x=50) % R2() % C1(y=50, c2=1) % C1(c2=1), bind_R1R2_C1C1_kf, bind_R1R2_C1C1_kr),
Parameter('bind_R1R2_C1C1_kf', 0.0001),
Parameter('bind_R1R2_C1C1_kr', 0.1),
Rule('bind_R1R2_C2', R1(x=None) % R2() + C2(y=None) | R1(x=50) % R2() % C2(y=50), bind_R1R2_C2_kf, bind_R1R2_C2_kr),
Parameter('bind_R1R2_C2_kf', 0.0002),
Parameter('bind_R1R2_C2_kr', 0.2),
Rule('bind_C1C1_R2', C1(y=None, c2=1) % C1(c2=1) + R2(x=None) | C1(y=50, c2=1) % C1(c2=1) % R2(x=50), bind_C1C1_R2_kf, bind_C1C1_R2_kr),
Parameter('bind_C1C1_R2_kf', 0.0003),
Parameter('bind_C1C1_R2_kr', 0.3),
])
"""
# extract species lists and matrix of rates
s_rows = [row[0] for row in bindtable[1:]]
s_cols = bindtable[0]
kmatrix = [row[1:] for row in bindtable[1:]]
# ensure kf is passed when necessary
kiter = itertools.chain.from_iterable(kmatrix)
if any(isinstance(x, numbers.Real) for x in kiter) and kf is None:
raise ValueError("must specify kf when using single kd values")
# loop over interactions
components = ComponentSet()
for r, s_row in enumerate(s_rows):
for c, s_col in enumerate(s_cols):
klist = kmatrix[r][c]
if klist is not None:
# if user gave a single kd, calculate kr
if isinstance(klist, numbers.Real):
kd = klist
klist = (kf, kd*kf)
components |= bind_complex(s_row, row_site, s_col, col_site, klist, m1, m2)
return components
def create_t_obs():
"""
Generate a rule to simulate passing of time and create a time observable
that can be used in complex Expression rates.
.. warning::
This macro is usually used to create rate laws that depend on time.
Time tracking rate laws using this macro only work for deterministic simulations.
Returns
-------
components : ComponentSet
The generated components. Contains the time monomer, Parameter rate of time creation,
Rule to simulate passing of time, ime Observable.
Examples
--------
Create rule to simulate passing of time and time observable::
Model()
create_t_obs()
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> create_t_obs()
ComponentSet([
Rule('synthesize___t', None >> __t(), __k_t),
Monomer('__t'),
Parameter('__k_t', 1.0),
Observable('t', __t()),
])
"""
# Add a time monomer and reaction to be able to create an observable to
# track the time within the simulation
time = Monomer('__t')
k_time = Parameter('__k_t', 1)
time_obs = Observable('t', time())
components = synthesize(time(), k_time)
components |= [time, k_time, time_obs]
return components
[docs]def drug_binding(drug, d_site, substrate, s_site, t_action, klist):
"""
Generate the reversible binding reaction DRUG + SUBSTRATE | DRUG:SUBSTRATE
that only gets triggered when the simulation reaches the time point t_action.
The idea of this macro is to mimic experimental settings when a reaction is
started and later on some kind of perturbation is added to the system.
.. warning::
This macro only works when a model is simulated using a deterministic simulator.
Parameters
----------
drug, substrate: Monomer or MonomerPattern
Monomers participating in the binding reaction.
d_site, s_site: string
The names of the sites on s1 and s2 used for binding.
t_action: float
Time of the simulation at which the drug is added
klist: list of 2 Parameters or list of 2 numbers
Forward and reverse rate constants (in that order). If Parameters are
passed, they will be used directly in the generated Rules. If numbers
are passed, Parameters will be created with automatically generated
names based on the names and states of S1 and S2 and these parameters
will be included at the end of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains the bidirectional binding Rule,
the time monomer, Parameter rate of time creation, Rule to simulate passing of time,
time Observable, two Expression rates that take into account when the interaction
between the drug and that substrate start to occur and optionally two Parameters
if klist was given as numbers
as numbers
Examples
--------
Binding between drug and substrate::
Model()
Monomer('drug', ['b'])
Monomer('substrate', ['b'])
drug_binding(drug(), 'b', substrate(), 'b', 10, [2,4])
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('drug', ['b'])
Monomer('drug', ['b'])
>>> Monomer('substrate', ['b'])
Monomer('substrate', ['b'])
>>> drug_binding(drug(), 'b', substrate(), 'b', 10, [0.1, 0.01])
ComponentSet([
Rule('bind_drug_substrate_to_drugsubstrate', drug(b=None) + substrate(b=None) | drug(b=1) % substrate(b=1), kf_expr_drug_substrate, kr_expr_drug_substrate),
Parameter('kf_drug_substrate', 0.1),
Parameter('kr_drug_substrate', 0.01),
Rule('synthesize___t', None >> __t(), __k_t),
Monomer('__t'),
Parameter('__k_t', 1.0),
Observable('t', __t()),
Expression('kf_expr_drug_substrate', Piecewise((kf_drug_substrate, t > 10), (0, True))),
Expression('kr_expr_drug_substrate', Piecewise((kr_drug_substrate, t > 10), (0, True))),
])
"""
_verify_sites(drug, d_site)
_verify_sites(substrate, s_site)
# Create a time observable using the create_t_obs macro
components_time_obs = create_t_obs()
time_obs = components_time_obs.t
# Set up some aliases to the patterns we'll use in the rules
drug_free = drug({d_site: None})
# retain any existing state for substrate's s_site, otherwise set it to None
if s_site in substrate.site_conditions:
substrate_free = substrate()
s_state = (substrate.site_conditions[s_site], 1)
else:
substrate_free = substrate({s_site: None})
s_state = 1
ds_complex = drug({d_site: 1}) % substrate({s_site: s_state})
substrate_monomer_name = substrate.monomer.name
drug_monomer_name = drug.monomer.name
if all(isinstance(x, (Parameter, Expression)) for x in klist):
k1 = klist[0]
k2 = klist[1]
params_created = ComponentSet()
elif all(isinstance(x, numbers.Real) for x in klist):
k1 = Parameter('kf_{0}_{1}'.format(drug_monomer_name, substrate_monomer_name), klist[0])
params_created = ComponentSet([k1])
k2 = Parameter('kr_{0}_{1}'.format(drug_monomer_name, substrate_monomer_name), klist[1])
params_created.add(k2)
else:
raise ValueError("klist must contain Parameters, Expressions, or numbers.")
kf_expr = Expression('kf_expr_{0}_{1}'.format(drug_monomer_name,
substrate_monomer_name), Piecewise((k1, time_obs > t_action), (0, True)))
kr_expr = Expression('kr_expr_{0}_{1}'.format(drug_monomer_name,
substrate_monomer_name), Piecewise((k2, time_obs > t_action), (0, True)))
bind_kpars = [kf_expr, kr_expr]
components_added_macro = components_time_obs
components_added_macro |= bind_kpars
components = _macro_rule('bind', drug_free + substrate_free | ds_complex, bind_kpars, ['kf', 'kr'])
components |= params_created
components |= components_added_macro
return components
# Catalysis
# =========
[docs]def catalyze(enzyme, e_site, substrate, s_site, product, klist):
"""
Generate the two-step catalytic reaction E + S | E:S >> E + P.
Parameters
----------
enzyme, substrate, product : Monomer or MonomerPattern
E, S and P in the above reaction.
e_site, s_site : string
The names of the sites on `enzyme` and `substrate` (respectively) where
they bind each other to form the E:S complex.
klist : list of 3 Parameters or list of 3 numbers
Forward, reverse and catalytic rate constants (in that order). If
Parameters are passed, they will be used directly in the generated
Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the names and states of enzyme,
substrate and product and these parameters will be included at the end
of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains two Rules (bidirectional complex
formation and unidirectional product dissociation), and optionally three
Parameters if klist was given as plain numbers.
Notes
-----
When passing a MonomerPattern for `enzyme` or `substrate`, do not include
`e_site` or `s_site` in the respective patterns. The macro will handle this.
Examples
--------
Using distinct Monomers for substrate and product::
Model()
Monomer('E', ['b'])
Monomer('S', ['b'])
Monomer('P')
catalyze(E(), 'b', S(), 'b', P(), (1e-4, 1e-1, 1))
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('E', ['b'])
Monomer('E', ['b'])
>>> Monomer('S', ['b'])
Monomer('S', ['b'])
>>> Monomer('P')
Monomer('P')
>>> catalyze(E(), 'b', S(), 'b', P(), (1e-4, 1e-1, 1)) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_E_S_to_ES', E(b=None) + S(b=None) | E(b=1) % S(b=1),
bind_E_S_to_ES_kf, bind_E_S_to_ES_kr),
Parameter('bind_E_S_to_ES_kf', 0.0001),
Parameter('bind_E_S_to_ES_kr', 0.1),
Rule('catalyze_ES_to_E_P', E(b=1) % S(b=1) >> E(b=None) + P(),
catalyze_ES_to_E_P_kc),
Parameter('catalyze_ES_to_E_P_kc', 1.0),
])
Using a single Monomer for substrate and product with a state change::
Monomer('Kinase', ['b'])
Monomer('Substrate', ['b', 'y'], {'y': ('U', 'P')})
catalyze(Kinase(), 'b', Substrate(y='U'), 'b', Substrate(y='P'),
(1e-4, 1e-1, 1))
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Kinase', ['b'])
Monomer('Kinase', ['b'])
>>> Monomer('Substrate', ['b', 'y'], {'y': ('U', 'P')})
Monomer('Substrate', ['b', 'y'], {'y': ('U', 'P')})
>>> catalyze(Kinase(), 'b', Substrate(y='U'), 'b', Substrate(y='P'), (1e-4, 1e-1, 1)) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_Kinase_SubstrateU_to_KinaseSubstrateU',
Kinase(b=None) + Substrate(b=None, y='U') | Kinase(b=1) % Substrate(b=1, y='U'),
bind_Kinase_SubstrateU_to_KinaseSubstrateU_kf,
bind_Kinase_SubstrateU_to_KinaseSubstrateU_kr),
Parameter('bind_Kinase_SubstrateU_to_KinaseSubstrateU_kf', 0.0001),
Parameter('bind_Kinase_SubstrateU_to_KinaseSubstrateU_kr', 0.1),
Rule('catalyze_KinaseSubstrateU_to_Kinase_SubstrateP',
Kinase(b=1) % Substrate(b=1, y='U') >> Kinase(b=None) + Substrate(b=None, y='P'),
catalyze_KinaseSubstrateU_to_Kinase_SubstrateP_kc),
Parameter('catalyze_KinaseSubstrateU_to_Kinase_SubstrateP_kc', 1.0),
])
"""
_verify_sites(enzyme, e_site)
_verify_sites(substrate, s_site)
# Set up some aliases to the patterns we'll use in the rules
enzyme_free = enzyme({e_site: None})
# retain any existing state for substrate's s_site, otherwise set it to None
if s_site in substrate.site_conditions:
substrate_free = substrate()
s_state = (substrate.site_conditions[s_site], 1)
else:
substrate_free = substrate({s_site: None})
s_state = 1
es_complex = enzyme({e_site: 1}) % substrate({s_site: s_state})
# If product is actually a variant of substrate, we need to explicitly say
# that it is no longer bound to enzyme, unless product already specifies a
# state for s_site.
if product().monomer is substrate().monomer \
and s_site not in product.site_conditions:
product = product({s_site: None})
# create the rules
components = _macro_rule('bind',
enzyme_free + substrate_free | es_complex,
klist[0:2], ['kf', 'kr'])
components |= _macro_rule('catalyze',
es_complex >> enzyme_free + product,
[klist[2]], ['kc'])
return components
[docs]def catalyze_complex(enzyme, e_site, substrate, s_site, product, klist, m1=None, m2=None):
""" Generate the two-step catalytic reaction E + S | E:S >> E + P, while allowing complexes to serve as enzyme, substrate and/or product.
E:S1 + S:S2 | E:S1:S:S2 >> E:S1 + P:S2
Parameters
----------
enzyme, substrate, product : Monomer, MonomerPattern, or ComplexPattern
Monomers or complexes participating in the binding reaction.
e_site, s_site : string
The names of the sites on `enzyme` and `substrate` (respectively) where
they bind each other to form the E:S complex.
klist : list of 3 Parameters or list of 3 numbers
Forward, reverse and catalytic rate constants (in that order). If
Parameters are passed, they will be used directly in the generated
Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the names and states of enzyme,
substrate and product and these parameters will be included at the end
of the returned component list.
m1, m2 : Monomer or MonomerPattern
If enzyme or substrate binding site is present in multiple monomers
within a complex, the specific monomer desired for binding must be specified.
Returns
-------
components : ComponentSet
The generated components. Contains the bidirectional binding Rule
and optionally three Parameters if klist was given as numbers.
"""
if isinstance(m1, Monomer):
m1 = m1()
if isinstance(m2, Monomer):
m2 = m2()
def build_complex(s1, site1, m1):
_verify_sites_complex(s1, site1)
#Retrieve a dictionary specifying the MonomerPattern within the complex that contains the given binding site.
specsitesdict = _verify_sites_complex(s1, site1)
s1complexpatub, s1complexpatb = check_sites_comp_build(s1, site1, m1, specsitesdict)
return s1complexpatb, s1complexpatub
def check_sites_comp_build(s1, site1, m1, specsitesdict):
#Return error if binding site exists on multiple monomers and a monomer for binding (m1) hasn't been specified.
if len(specsitesdict) > 1 and m1==None:
raise ValueError("Binding site '%s' present in more than one monomer in complex '%s'. Specify variable m1, the monomer used for binding within the complex." % (site1, s1))
if not s1.is_concrete:
raise ValueError("Complex '%s' must be concrete." % (s1))
#If the given binding site is only present in one monomer in the complex:
if m1==None:
#Build up ComplexPattern for use in rule (with state of given binding site specified).
s1complexpatub = list(specsitesdict.keys())[0]({site1:None})
s1complexpatb = list(specsitesdict.keys())[0]({site1:50})
for monomer in s1.monomer_patterns:
if monomer not in specsitesdict.keys():
s1complexpatub %= monomer
s1complexpatb %= monomer
#If the binding site is present on more than one monomer in the complex, the monomer must be specified by the user. Use specified m1 to build ComplexPattern.
else:
#Make sure binding states of MonomerPattern m1 match those of the monomer within the ComplexPattern s1 (ComplexPattern monomer takes precedence if not).
i = 0
identical_monomers = []
for mon in s1.monomer_patterns:
#Only change the binding site for the first monomer that matches. Keep any others unchanged to add to final complex that is returned.
if mon.monomer.name == m1.monomer.name:
i += 1
if i == 1:
s1complexpatub = m1({site1:None})
s1complexpatb = m1({site1:50})
else:
identical_monomers.append(mon)
#Build up ComplexPattern for use in rule (with state of given binding site on m1 specified).
for mon in s1.monomer_patterns:
if mon.monomer.name != m1.monomer.name:
s1complexpatub %= mon
s1complexpatb %= mon
if identical_monomers:
for i in range(len(identical_monomers)):
s1complexpatub %= identical_monomers[i]
s1complexpatb %= identical_monomers[i]
return s1complexpatub, s1complexpatb
#If no complexes exist in the reaction, revert to catalyze().
if (isinstance(enzyme, MonomerPattern) or isinstance(enzyme, Monomer)) and (isinstance(substrate, MonomerPattern) or isinstance(substrate, Monomer)):
_verify_sites(enzyme, e_site)
_verify_sites(substrate, s_site)
return catalyze(enzyme, e_site, substrate, s_site, product, klist,)
# Build E:S
if isinstance(enzyme, ComplexPattern):
enzymepatb, enzyme_free = build_complex(enzyme, e_site, m1)
else:
enzymepatb, enzyme_free = enzyme({e_site: 1}), enzyme({e_site: None})
if isinstance(substrate, ComplexPattern):
substratepatb, substratepatub = build_complex(substrate, s_site, m2)
else:
substratepatb = substrate({s_site: 50})
"""if s_site in substrate.site_conditions:
substrate_free = substrate()
s_state = (substrate.site_conditions[s_site], 1)
else:
substrate_free = substrate({s_site: None})
s_state = 1
substratepatb = substrate({s_site: s_state})
"""
es_complex = enzymepatb % substratepatb
# Use bind complex to binding rule.
components = bind_complex(enzyme, e_site, substrate, s_site, klist[0:2], m1, m2)
components |= _macro_rule('catalyze',
es_complex >> enzyme_free + product,
[klist[2]], ['kc'])
return components
[docs]def catalyze_state(enzyme, e_site, substrate, s_site, mod_site,
state1, state2, klist):
"""
Generate the two-step catalytic reaction E + S | E:S >> E + P. A wrapper
around catalyze() with a signature specifying the state change of the
substrate resulting from catalysis.
Parameters
----------
enzyme : Monomer or MonomerPattern
E in the above reaction.
substrate : Monomer or MonomerPattern
S and P in the above reaction. The product species is assumed to be
identical to the substrate species in all respects except the state
of the modification site. The state of the modification site should
not be specified in the MonomerPattern for the substrate.
e_site, s_site : string
The names of the sites on `enzyme` and `substrate` (respectively) where
they bind each other to form the E:S complex.
mod_site : string
The name of the site on the substrate that is modified by catalysis.
state1, state2 : strings
The states of the modification site (mod_site) on the substrate before
(state1) and after (state2) catalysis.
klist : list of 3 Parameters or list of 3 numbers
Forward, reverse and catalytic rate constants (in that order). If
Parameters are passed, they will be used directly in the generated
Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the names and states of enzyme,
substrate and product and these parameters will be included at the end
of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains two Rules (bidirectional complex
formation and unidirectional product dissociation), and optionally three
Parameters if klist was given as plain numbers.
Notes
-----
When passing a MonomerPattern for `enzyme` or `substrate`, do not include
`e_site` or `s_site` in the respective patterns. In addition, do not
include the state of the modification site on the substrate. The macro
will handle this.
Examples
--------
Using a single Monomer for substrate and product with a state change::
Monomer('Kinase', ['b'])
Monomer('Substrate', ['b', 'y'], {'y': ('U', 'P')})
catalyze_state(Kinase, 'b', Substrate, 'b', 'y', 'U', 'P',
(1e-4, 1e-1, 1))
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Kinase', ['b'])
Monomer('Kinase', ['b'])
>>> Monomer('Substrate', ['b', 'y'], {'y': ('U', 'P')})
Monomer('Substrate', ['b', 'y'], {'y': ('U', 'P')})
>>> catalyze_state(Kinase, 'b', Substrate, 'b', 'y', 'U', 'P', (1e-4, 1e-1, 1)) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('bind_Kinase_SubstrateU_to_KinaseSubstrateU',
Kinase(b=None) + Substrate(b=None, y='U') | Kinase(b=1) % Substrate(b=1, y='U'),
bind_Kinase_SubstrateU_to_KinaseSubstrateU_kf,
bind_Kinase_SubstrateU_to_KinaseSubstrateU_kr),
Parameter('bind_Kinase_SubstrateU_to_KinaseSubstrateU_kf', 0.0001),
Parameter('bind_Kinase_SubstrateU_to_KinaseSubstrateU_kr', 0.1),
Rule('catalyze_KinaseSubstrateU_to_Kinase_SubstrateP',
Kinase(b=1) % Substrate(b=1, y='U') >> Kinase(b=None) + Substrate(b=None, y='P'),
catalyze_KinaseSubstrateU_to_Kinase_SubstrateP_kc),
Parameter('catalyze_KinaseSubstrateU_to_Kinase_SubstrateP_kc', 1.0),
])
"""
return catalyze(enzyme, e_site, substrate({mod_site: state1}),
s_site, substrate({mod_site: state2}), klist)
[docs]def catalyze_one_step(enzyme, substrate, product, kf):
"""
Generate the one-step catalytic reaction E + S >> E + P.
Parameters
----------
enzyme, substrate, product : Monomer or MonomerPattern
E, S and P in the above reaction.
kf : a Parameter or a number
Forward rate constant for the reaction. If a
Parameter is passed, it will be used directly in the generated
Rules. If a number is passed, a Parameter will be created with an
automatically generated name based on the names and states of the
enzyme, substrate and product and this parameter will be included
at the end of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains the unidirectional reaction Rule
and optionally the forward rate Parameter if klist was given as a
number.
Notes
-----
In this macro, there is no direct binding between enzyme and substrate,
so binding sites do not have to be specified. This represents an
approximation for the case when the enzyme is operating in its linear
range. However, if catalysis is nevertheless contingent on the enzyme or
substrate being unbound on some site, then that information must be encoded
in the MonomerPattern for the enzyme or substrate. See the examples, below.
Examples
--------
Convert S to P by E::
Model()
Monomer('E', ['b'])
Monomer('S', ['b'])
Monomer('P')
catalyze_one_step(E, S, P, 1e-4)
If the ability of the enzyme E to catalyze this reaction is dependent
on the site 'b' of E being unbound, then this macro must be called as
catalyze_one_step(E(b=None), S, P, 1e-4)
and similarly if the substrate or product must be unbound.
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('E', ['b'])
Monomer('E', ['b'])
>>> Monomer('S', ['b'])
Monomer('S', ['b'])
>>> Monomer('P')
Monomer('P')
>>> catalyze_one_step(E, S, P, 1e-4) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('one_step_E_S_to_E_P', E() + S() >> E() + P(), one_step_E_S_to_E_P_kf),
Parameter('one_step_E_S_to_E_P_kf', 0.0001),
])
"""
if isinstance(enzyme, Monomer):
enzyme = enzyme()
if isinstance(substrate, Monomer):
substrate = substrate()
if isinstance(product, Monomer):
product = product()
return _macro_rule('one_step',
enzyme + substrate >> enzyme + product,
[kf], ['kf'])
[docs]def catalyze_one_step_reversible(enzyme, substrate, product, klist):
"""
Create fwd and reverse rules for catalysis of the form::
E + S -> E + P
P -> S
Parameters
----------
enzyme, substrate, product : Monomer or MonomerPattern
E, S and P in the above reactions.
klist : list of 2 Parameters or list of 2 numbers
A list containing the rate constant for catalysis and the rate constant
for the conversion of product back to substrate (in that order). If
Parameters are passed, they will be used directly in the generated
Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the names and states of S1 and
S2 and these parameters will be included at the end of the returned
component list.
Returns
-------
components : ComponentSet
The generated components. Contains two rules (the single-step catalysis
rule and the product reversion rule) and optionally the two generated
Parameter objects if klist was given as numbers.
Notes
-----
Calls the macro catalyze_one_step to generate the catalysis rule.
Examples
--------
One-step, pseudo-first order conversion of S to P by E::
Model()
Monomer('E', ['b'])
Monomer('S', ['b'])
Monomer('P')
catalyze_one_step_reversible(E, S, P, [1e-1, 1e-4])
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('E', ['b'])
Monomer('E', ['b'])
>>> Monomer('S', ['b'])
Monomer('S', ['b'])
>>> Monomer('P')
Monomer('P')
>>> catalyze_one_step_reversible(E, S, P, [1e-1, 1e-4]) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('one_step_E_S_to_E_P', E() + S() >> E() + P(), one_step_E_S_to_E_P_kf),
Parameter('one_step_E_S_to_E_P_kf', 0.1),
Rule('reverse_P_to_S', P() >> S(), reverse_P_to_S_kr),
Parameter('reverse_P_to_S_kr', 0.0001),
])
"""
if isinstance(enzyme, Monomer):
enzyme = enzyme()
if isinstance(substrate, Monomer):
substrate = substrate()
if isinstance(product, Monomer):
product = product()
components = catalyze_one_step(enzyme, substrate, product, klist[0])
components |= _macro_rule('reverse', product >> substrate,
[klist[1]], ['kr'])
return components
# Synthesis and degradation
# =========================
[docs]def synthesize(species, ksynth):
"""
Generate a reaction which synthesizes a species.
Note that `species` must be "concrete", i.e. the state of all
sites in all of its monomers must be specified. No site may be
left unmentioned.
Parameters
----------
species : Monomer, MonomerPattern or ComplexPattern
The species to synthesize. If a Monomer, sites are considered
as unbound and in their default state. If a pattern, must be
concrete.
ksynth : Parameters or number
Synthesis rate. If a Parameter is passed, it will be used directly in
the generated Rule. If a number is passed, a Parameter will be created
with an automatically generated name based on the names and site states
of the components of `species` and this parameter will be included at
the end of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains the unidirectional synthesis Rule and
optionally a Parameter if ksynth was given as a number.
Examples
--------
Synthesize A with site x unbound and site y in state 'e'::
Model()
Monomer('A', ['x', 'y'], {'y': ['e', 'f']})
synthesize(A(x=None, y='e'), 1e-4)
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A', ['x', 'y'], {'y': ['e', 'f']})
Monomer('A', ['x', 'y'], {'y': ['e', 'f']})
>>> synthesize(A(x=None, y='e'), 1e-4) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('synthesize_Ae', None >> A(x=None, y='e'), synthesize_Ae_k),
Parameter('synthesize_Ae_k', 0.0001),
])
"""
def synthesize_name_func(rule_expression):
cps = rule_expression.product_pattern.complex_patterns
return '_'.join(_complex_pattern_label(cp) for cp in cps)
if isinstance(species, Monomer):
species = species()
species = as_complex_pattern(species)
if not species.is_concrete():
raise ValueError("species must be concrete")
return _macro_rule('synthesize', None >> species, [ksynth], ['k'],
name_func=synthesize_name_func)
[docs]def degrade(species, kdeg):
"""
Generate a reaction which degrades a species.
Note that `species` is not required to be "concrete".
Parameters
----------
species : Monomer, MonomerPattern or ComplexPattern
The species to synthesize. If a Monomer, sites are considered
as unbound and in their default state. If a pattern, must be
concrete.
kdeg : Parameters or number
Degradation rate. If a Parameter is passed, it will be used directly in
the generated Rule. If a number is passed, a Parameter will be created
with an automatically generated name based on the names and site states
of the components of `species` and this parameter will be included at
the end of the returned component list.
Returns
-------
components : ComponentSet
The generated components. Contains the unidirectional degradation Rule
and optionally a Parameter if ksynth was given as a number.
Examples
--------
Degrade all B, even bound species::
Model()
Monomer('B', ['x'])
degrade(B(), 1e-6)
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('B', ['x'])
Monomer('B', ['x'])
>>> degrade(B(), 1e-6) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('degrade_B', B() >> None, degrade_B_k),
Parameter('degrade_B_k', 1e-06),
])
"""
def degrade_name_func(rule_expression):
cps = rule_expression.reactant_pattern.complex_patterns
return '_'.join(_complex_pattern_label(cp) for cp in cps)
if isinstance(species, Monomer):
species = species()
species = as_complex_pattern(species)
return _macro_rule('degrade', species >> None, [kdeg], ['k'],
name_func=degrade_name_func)
[docs]def synthesize_degrade_table(table):
"""
Generate a table of synthesis and degradation reactions.
Given a list of species, calls the `synthesize` and `degrade` macros on each
one. The species and the parameter values are passed as a list of lists
(i.e. a table) with each inner list consisting of the species, forward and
reverse rates (in that order).
Each species' associated pair of rates may be either Parameters or
numbers. If Parameters are passed, they will be used directly in the
generated Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the names and states of the relevant
species and these parameters will be included in the returned component
list. To omit any individual reaction, pass None in place of the
corresponding parameter.
Note that any `species` with a non-None synthesis rate must be "concrete".
Parameters
----------
table : list of lists
Table of species and rates, as described above.
Returns
-------
components : ComponentSet
The generated components. Contains the unidirectional synthesis and
degradation Rules and optionally the Parameters for any rates given as
numbers.
Examples
--------
Specify synthesis and degradation reactions for A and B in a table::
Model()
Monomer('A', ['x', 'y'], {'y': ['e', 'f']})
Monomer('B', ['x'])
synthesize_degrade_table([[A(x=None, y='e'), 1e-4, 1e-6],
[B(), None, 1e-7]])
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('A', ['x', 'y'], {'y': ['e', 'f']})
Monomer('A', ['x', 'y'], {'y': ['e', 'f']})
>>> Monomer('B', ['x'])
Monomer('B', ['x'])
>>> synthesize_degrade_table([[A(x=None, y='e'), 1e-4, 1e-6],
... [B(), None, 1e-7]]) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('synthesize_Ae', None >> A(x=None, y='e'), synthesize_Ae_k),
Parameter('synthesize_Ae_k', 0.0001),
Rule('degrade_Ae', A(x=None, y='e') >> None, degrade_Ae_k),
Parameter('degrade_Ae_k', 1e-06),
Rule('degrade_B', B() >> None, degrade_B_k),
Parameter('degrade_B_k', 1e-07),
])
"""
# loop over interactions
components = ComponentSet()
for row in table:
species, ksynth, kdeg = row
if ksynth is not None:
components |= synthesize(species, ksynth)
if kdeg is not None:
components |= degrade(species, kdeg)
return components
# Polymer assembly (pores/rings and chains)
# =========================================
def polymer_species(subunit, site1, site2, size, closed=False):
"""
Return a ComplexPattern representing a linear or closed circular polymer.
Parameters
----------
subunit : Monomer or MonomerPattern
The subunit of which the polymer is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
size : integer
The number of subunits in the polymer.
closed : boolean
If False (default), the polymer is linear, with unbound sites at each
end. If True, the polymer is a closed circle, like a ring or pore.
Returns
-------
A ComplexPattern corresponding to the polymer.
Notes
-----
Used by both chain_species and pore_species.
"""
_verify_sites(subunit, site1, site2)
if size <= 0:
raise ValueError("size must be an integer greater than 0")
if size == 1:
polymer = subunit({site1: None, site2: None})
elif size == 2:
polymer = subunit({site1: None, site2: 1}) % \
subunit({site1: 1, site2: None})
else:
# If a closed circle, use 0 as the bond number for the "seam";
# if linear, use None for the unbound ends
seam_site_num = size if closed else None
# First subunit
polymer = subunit({site1: seam_site_num, site2: 1})
# Build up the ComplexPattern for the polymer, starting with the first
# subunit
for i in range(1, size-1):
polymer %= subunit({site1: i, site2: i + 1})
# Attach the last subunit
polymer %= subunit({site1: size-1, site2: seam_site_num})
# Set ComplexPattern to MatchOnce
polymer.match_once = True
return polymer
def assemble_polymer_sequential(subunit, site1, site2, max_size, ktable,
closed=False):
"""Generate rules to assemble a polymer by sequential subunit addition.
The polymer species are created by sequential addition of `subunit` monomers,
i.e. larger oligomeric species never fuse together. The polymer structure is
defined by the `polymer_species` macro.
Parameters
----------
subunit : Monomer or MonomerPattern
The subunit of which the polymer is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
max_size : integer
The maximum number of subunits in the polymer.
ktable : list of lists of Parameters or numbers
Table of forward and reverse rate constants for the assembly steps. The
outer list must be of length `max_size` - 1, and the inner lists must
all be of length 2. In the outer list, the first element corresponds to
the first assembly step in which two monomeric subunits bind to form a
2-subunit complex, and the last element corresponds to the final step in
which the `max_size`th subunit is added. Each inner list contains the
forward and reverse rate constants (in that order) for the corresponding
assembly reaction, and each of these pairs must comprise solely
Parameter objects or solely numbers (never one of each). If Parameters
are passed, they will be used directly in the generated Rules. If
numbers are passed, Parameters will be created with automatically
generated names based on `subunit`, `site1`, `site2` and the polymer sizes
and these parameters will be included at the end of the returned
component list.
closed : boolean
If False (default), assembles a linear (non-circular) polymer. If True,
assembles a circular ring/pore polymer.
Notes
-----
See documentation for :py:func:`assemble_chain_sequential` and
:py:func:`assemble_pore_sequential` for examples.
"""
if len(ktable) != max_size - 1:
raise ValueError("len(ktable) must be equal to max_size - 1")
def polymer_rule_name(rule_expression, size):
react_p = rule_expression.reactant_pattern
monomer = react_p.complex_patterns[0].monomer_patterns[0].monomer
return '%s_%d' % (monomer.name, size)
components = ComponentSet()
s = polymer_species(subunit, site1, site2, 1, closed=closed)
for size, klist in zip(range(2, max_size + 1), ktable):
polymer_prev = polymer_species(subunit, site1, site2, size - 1,
closed=closed)
polymer_next = polymer_species(subunit, site1, site2, size,
closed=closed)
name_func = functools.partial(polymer_rule_name, size=size)
rule_name_base = 'assemble_%s_sequential' % \
('pore' if closed else 'chain')
components |= _macro_rule(rule_name_base,
s + polymer_prev | polymer_next,
klist, ['kf', 'kr'],
name_func=name_func)
return components
# Pore assembly
# =============
def pore_species(subunit, site1, site2, size):
"""
Return a ComplexPattern representing a circular homomeric pore.
Parameters
----------
subunit : Monomer or MonomerPattern
The subunit of which the pore is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
size : integer
The number of subunits in the pore.
Returns
-------
A ComplexPattern corresponding to the pore.
Notes
-----
At sizes 1 and 2 the ring is not closed, i.e. there is one site1 and one
site2 which remain unbound. At size 3 and up the ring is closed and all
site1 sites are bound to a site2.
Examples
--------
Get the ComplexPattern object representing a pore of size 4::
Model()
Monomer('Unit', ['p1', 'p2'])
pore_tetramer = pore_species(Unit, 'p1', 'p2', 4)
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2'])
Monomer('Unit', ['p1', 'p2'])
>>> pore_species(Unit, 'p1', 'p2', 4)
MatchOnce(Unit(p1=4, p2=1) % Unit(p1=1, p2=2) % Unit(p1=2, p2=3) % Unit(p1=3, p2=4))
"""
return polymer_species(subunit, site1, site2, size, closed=True)
[docs]def assemble_pore_sequential(subunit, site1, site2, max_size, ktable):
"""Generate rules to assemble a circular homomeric pore sequentially.
The pore species are created by sequential addition of `subunit` monomers,
i.e. larger oligomeric species never fuse together. The pore structure is
defined by the `pore_species` macro.
Parameters
----------
subunit : Monomer or MonomerPattern
The subunit of which the pore is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
max_size : integer
The maximum number of subunits in the pore.
ktable : list of lists of Parameters or numbers
Table of forward and reverse rate constants for the assembly steps. The
outer list must be of length `max_size` - 1, and the inner lists must
all be of length 2. In the outer list, the first element corresponds to
the first assembly step in which two monomeric subunits bind to form a
2-subunit complex, and the last element corresponds to the final step in
which the `max_size`th subunit is added. Each inner list contains the
forward and reverse rate constants (in that order) for the corresponding
assembly reaction, and each of these pairs must comprise solely
Parameter objects or solely numbers (never one of each). If Parameters
are passed, they will be used directly in the generated Rules. If
numbers are passed, Parameters will be created with automatically
generated names based on `subunit`, `site1`, `site2` and the pore sizes
and these parameters will be included at the end of the returned
component list.
Examples
--------
Assemble a three-membered pore by sequential addition of monomers,
with the same forward/reverse rates for monomer-monomer and monomer-dimer
interactions::
Model()
Monomer('Unit', ['p1', 'p2'])
assemble_pore_sequential(Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2)
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2'])
Monomer('Unit', ['p1', 'p2'])
>>> assemble_pore_sequential(Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('assemble_pore_sequential_Unit_2',
Unit(p1=None, p2=None) + Unit(p1=None, p2=None) |
Unit(p1=None, p2=1) % Unit(p1=1, p2=None),
assemble_pore_sequential_Unit_2_kf,
assemble_pore_sequential_Unit_2_kr),
Parameter('assemble_pore_sequential_Unit_2_kf', 0.0001),
Parameter('assemble_pore_sequential_Unit_2_kr', 0.1),
Rule('assemble_pore_sequential_Unit_3',
Unit(p1=None, p2=None) + Unit(p1=None, p2=1) % Unit(p1=1, p2=None) |
MatchOnce(Unit(p1=3, p2=1) % Unit(p1=1, p2=2) % Unit(p1=2, p2=3)),
assemble_pore_sequential_Unit_3_kf,
assemble_pore_sequential_Unit_3_kr),
Parameter('assemble_pore_sequential_Unit_3_kf', 0.0001),
Parameter('assemble_pore_sequential_Unit_3_kr', 0.1),
])
"""
return assemble_polymer_sequential(subunit, site1, site2, max_size, ktable,
closed=True)
[docs]def pore_transport(subunit, sp_site1, sp_site2, sc_site, min_size, max_size,
csource, c_site, cdest, ktable):
"""
Generate rules to transport cargo through a circular homomeric pore.
The pore structure is defined by the `pore_species` macro -- `subunit`
monomers bind to each other from `sp_site1` to `sp_site2` to form a closed
ring. The transport reaction is modeled as a catalytic process of the form
pore + csource | pore:csource >> pore + cdest
Parameters
----------
subunit : Monomer or MonomerPattern
Subunit of which the pore is composed.
sp_site1, sp_site2 : string
Names of the sites where one copy of `subunit` binds to the next.
sc_site : string
Name of the site on `subunit` where it binds to the cargo `csource`.
min_size, max_size : integer
Minimum and maximum number of subunits in the pore at which transport
will occur.
csource : Monomer or MonomerPattern
Cargo "source", i.e. the entity to be transported.
c_site : string
Name of the site on `csource` where it binds to `subunit`.
cdest : Monomer or MonomerPattern
Cargo "destination", i.e. the resulting state after the transport event.
ktable : list of lists of Parameters or numbers
Table of forward, reverse and catalytic rate constants for the transport
reactions. The outer list must be of length `max_size` - `min_size` + 1,
and the inner lists must all be of length 3. In the outer list, the
first element corresponds to the transport through the pore of size
`min_size` and the last element to that of size `max_size`. Each inner
list contains the forward, reverse and catalytic rate constants (in that
order) for the corresponding transport reaction, and each of these pairs
must comprise solely Parameter objects or solely numbers (never some of
each). If Parameters are passed, they will be used directly in the
generated Rules. If numbers are passed, Parameters will be created with
automatically generated names based on the subunit, the pore size and
the cargo, and these parameters will be included at the end of the
returned component list.
Examples
--------
Specify that a three-membered pore is capable of
transporting cargo from the mitochondria to the cytoplasm::
Model()
Monomer('Unit', ['p1', 'p2', 'sc_site'])
Monomer('Cargo', ['c_site', 'loc'], {'loc':['mito', 'cyto']})
pore_transport(Unit, 'p1', 'p2', 'sc_site', 3, 3,
Cargo(loc='mito'), 'c_site', Cargo(loc='cyto'),
[[1e-4, 1e-1, 1]])
Generates two rules--one (reversible) binding rule and one transport
rule--and the three associated parameters.
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2', 'sc_site'])
Monomer('Unit', ['p1', 'p2', 'sc_site'])
>>> Monomer('Cargo', ['c_site', 'loc'], {'loc':['mito', 'cyto']})
Monomer('Cargo', ['c_site', 'loc'], {'loc': ['mito', 'cyto']})
>>> pore_transport(Unit, 'p1', 'p2', 'sc_site', 3, 3,
... Cargo(loc='mito'), 'c_site', Cargo(loc='cyto'),
... [[1e-4, 1e-1, 1]]) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('pore_transport_complex_Unit_3_Cargomito',
MatchOnce(Unit(p1=3, p2=1, sc_site=None) %
Unit(p1=1, p2=2, sc_site=None) %
Unit(p1=2, p2=3, sc_site=None)) +
Cargo(c_site=None, loc='mito') |
MatchOnce(Unit(p1=3, p2=1, sc_site=4) %
Unit(p1=1, p2=2, sc_site=None) %
Unit(p1=2, p2=3, sc_site=None) %
Cargo(c_site=4, loc='mito')),
pore_transport_complex_Unit_3_Cargomito_kf,
pore_transport_complex_Unit_3_Cargomito_kr),
Parameter('pore_transport_complex_Unit_3_Cargomito_kf', 0.0001),
Parameter('pore_transport_complex_Unit_3_Cargomito_kr', 0.1),
Rule('pore_transport_dissociate_Unit_3_Cargocyto',
MatchOnce(Unit(p1=3, p2=1, sc_site=4) %
Unit(p1=1, p2=2, sc_site=None) %
Unit(p1=2, p2=3, sc_site=None) %
Cargo(c_site=4, loc='mito')) >>
MatchOnce(Unit(p1=3, p2=1, sc_site=None) %
Unit(p1=1, p2=2, sc_site=None) %
Unit(p1=2, p2=3, sc_site=None)) +
Cargo(c_site=None, loc='cyto'),
pore_transport_dissociate_Unit_3_Cargocyto_kc),
Parameter('pore_transport_dissociate_Unit_3_Cargocyto_kc', 1.0),
])
"""
_verify_sites(subunit, sc_site)
_verify_sites(csource, c_site)
if len(ktable) != max_size - min_size + 1:
raise ValueError("len(ktable) must be equal to max_size - min_size + 1")
def pore_transport_rule_name(rule_expression, size):
# Get ReactionPatterns
react_p = rule_expression.reactant_pattern
prod_p = rule_expression.product_pattern
# Build the label components
# Pore is always first complex of LHS due to how we build the rules
subunit = react_p.complex_patterns[0].monomer_patterns[0]
if len(react_p.complex_patterns) == 2:
# This is the complexation reaction
cargo = react_p.complex_patterns[1].monomer_patterns[0]
else:
# This is the dissociation reaction
cargo = prod_p.complex_patterns[1].monomer_patterns[0]
return '%s_%d_%s' % (_monomer_pattern_label(subunit), size,
_monomer_pattern_label(cargo))
components = ComponentSet()
# Set up some aliases that are invariant with pore size
subunit_free = subunit({sc_site: None})
csource_free = csource({c_site: None})
# If cdest is actually a variant of csource, we need to explicitly say that
# it is no longer bound to the pore
if cdest().monomer is csource().monomer:
cdest = cdest({c_site: None})
for size, klist in zip(range(min_size, max_size + 1), ktable):
# More aliases which do depend on pore size
pore_free = pore_species(subunit_free, sp_site1, sp_site2, size)
# This one is a bit tricky. The pore:csource complex must only introduce
# one additional bond even though there are multiple subunits in the
# pore. We create partial patterns for bound pore and csource, using a
# bond number that is high enough not to conflict with the bonds within
# the pore ring itself.
# Start by copying pore_free, which has all cargo binding sites empty
pore_bound = pore_free.copy()
# Get the next bond number not yet used in the pore structure itself
cargo_bond_num = size + 1
# Assign that bond to the first subunit in the pore
pore_bound.monomer_patterns[0].site_conditions[sc_site] = cargo_bond_num
# Create a cargo source pattern with that same bond
csource_bound = csource({c_site: cargo_bond_num})
# Finally we can define the complex trivially; the bond numbers are
# already present in the patterns
pc_complex = pore_bound % csource_bound
# Create the rules (just like catalyze)
name_func = functools.partial(pore_transport_rule_name, size=size)
components |= _macro_rule('pore_transport_complex',
pore_free + csource_free | pc_complex,
klist[0:2], ['kf', 'kr'],
name_func=name_func)
components |= _macro_rule('pore_transport_dissociate',
pc_complex >> pore_free + cdest,
[klist[2]], ['kc'],
name_func=name_func)
return components
[docs]def pore_bind(subunit, sp_site1, sp_site2, sc_site, size, cargo, c_site,
klist):
"""
Generate rules to bind a monomer to a circular homomeric pore.
The pore structure is defined by the `pore_species` macro -- `subunit`
monomers bind to each other from `sp_site1` to `sp_site2` to form a closed
ring. The binding reaction takes the form pore + cargo | pore:cargo.
Parameters
----------
subunit : Monomer or MonomerPattern
Subunit of which the pore is composed.
sp_site1, sp_site2 : string
Names of the sites where one copy of `subunit` binds to the next.
sc_site : string
Name of the site on `subunit` where it binds to the cargo `cargo`.
size : integer
Number of subunits in the pore at which binding will occur.
cargo : Monomer or MonomerPattern
Cargo that binds to the pore complex.
c_site : string
Name of the site on `cargo` where it binds to `subunit`.
klist : list of Parameters or numbers
List containing forward and reverse rate constants for the binding
reaction (in that order). Rate constants should either be both Parameter
objects or both numbers. If Parameters are passed, they will be used
directly in the generated Rules. If numbers are passed, Parameters
will be created with automatically generated names based on the
subunit, the pore size and the cargo, and these parameters will be
included at the end of the returned component list.
Examples
--------
Specify that a cargo molecule can bind reversibly to a 3-membered
pore::
Model()
Monomer('Unit', ['p1', 'p2', 'sc_site'])
Monomer('Cargo', ['c_site'])
pore_bind(Unit, 'p1', 'p2', 'sc_site', 3,
Cargo(), 'c_site', [1e-4, 1e-1, 1])
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2', 'sc_site'])
Monomer('Unit', ['p1', 'p2', 'sc_site'])
>>> Monomer('Cargo', ['c_site'])
Monomer('Cargo', ['c_site'])
>>> pore_bind(Unit, 'p1', 'p2', 'sc_site', 3,
... Cargo(), 'c_site', [1e-4, 1e-1, 1]) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('pore_bind_Unit_3_Cargo',
MatchOnce(Unit(p1=3, p2=1, sc_site=None) %
Unit(p1=1, p2=2, sc_site=None) %
Unit(p1=2, p2=3, sc_site=None)) +
Cargo(c_site=None) |
MatchOnce(Unit(p1=3, p2=1, sc_site=4) %
Unit(p1=1, p2=2, sc_site=None) %
Unit(p1=2, p2=3, sc_site=None) %
Cargo(c_site=4)),
pore_bind_Unit_3_Cargo_kf, pore_bind_Unit_3_Cargo_kr),
Parameter('pore_bind_Unit_3_Cargo_kf', 0.0001),
Parameter('pore_bind_Unit_3_Cargo_kr', 0.1),
])
"""
_verify_sites(subunit, sc_site)
_verify_sites(cargo, c_site)
def pore_bind_rule_name(rule_expression, size):
# Get ReactionPatterns
react_p = rule_expression.reactant_pattern
prod_p = rule_expression.product_pattern
# Build the label components
# Pore is always first complex of LHS due to how we build the rules
subunit = react_p.complex_patterns[0].monomer_patterns[0].monomer
if len(react_p.complex_patterns) == 2:
# This is the complexation reaction
cargo = react_p.complex_patterns[1].monomer_patterns[0]
else:
# This is the dissociation reaction
cargo = prod_p.complex_patterns[1].monomer_patterns[0]
return '%s_%d_%s' % (subunit.name, size,
_monomer_pattern_label(cargo))
components = ComponentSet()
# Set up some aliases that are invariant with pore size
subunit_free = subunit({sc_site: None})
cargo_free = cargo({c_site: None})
#for size, klist in zip(range(min_size, max_size + 1), ktable):
# More aliases which do depend on pore size
pore_free = pore_species(subunit_free, sp_site1, sp_site2, size)
# This one is a bit tricky. The pore:cargo complex must only introduce
# one additional bond even though there are multiple subunits in the
# pore. We create partial patterns for bound pore and cargo, using a
# bond number that is high enough not to conflict with the bonds within
# the pore ring itself.
# Start by copying pore_free, which has all cargo binding sites empty
pore_bound = pore_free.copy()
# Get the next bond number not yet used in the pore structure itself
cargo_bond_num = size + 1
# Assign that bond to the first subunit in the pore
pore_bound.monomer_patterns[0].site_conditions[sc_site] = cargo_bond_num
# Create a cargo source pattern with that same bond
cargo_bound = cargo({c_site: cargo_bond_num})
# Finally we can define the complex trivially; the bond numbers are
# already present in the patterns
pc_complex = pore_bound % cargo_bound
# Create the rules
name_func = functools.partial(pore_bind_rule_name, size=size)
components |= _macro_rule('pore_bind',
pore_free + cargo_free | pc_complex,
klist[0:2], ['kf', 'kr'],
name_func=name_func)
return components
# Chain assembly
# =============
def chain_species(subunit, site1, site2, size):
"""
Return a ComplexPattern representing a linear, chained polymer.
Parameters
----------
subunit : Monomer or MonomerPattern
The subunit of which the chain is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
size : integer
The number of subunits in the chain.
Returns
-------
A ComplexPattern corresponding to the chain.
Notes
-----
Similar to pore_species, but never closes the chain.
Examples
--------
Get the ComplexPattern object representing a chain of length 4::
Model()
Monomer('Unit', ['p1', 'p2'])
chain_tetramer = chain_species(Unit, 'p1', 'p2', 4)
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2'])
Monomer('Unit', ['p1', 'p2'])
>>> chain_species(Unit, 'p1', 'p2', 4)
MatchOnce(Unit(p1=None, p2=1) % Unit(p1=1, p2=2) % Unit(p1=2, p2=3) % Unit(p1=3, p2=None))
"""
return polymer_species(subunit, site1, site2, size, closed=False)
def assemble_chain_sequential(subunit, site1, site2, max_size, ktable):
"""
Generate rules to assemble a homomeric chain sequentially.
The chain species are created by sequential addition of `subunit` monomers.
The chain structure is defined by the `chain_species` macro.
Parameters
----------
subunit : Monomer or MonomerPattern
The subunit of which the chain is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
max_size : integer
The maximum number of subunits in the chain.
ktable : list of lists of Parameters or numbers
Table of forward and reverse rate constants for the assembly steps. The
outer list must be of length `max_size` - 1, and the inner lists must
all be of length 2. In the outer list, the first element corresponds to
the first assembly step in which two monomeric subunits bind to form a
2-subunit complex, and the last element corresponds to the final step in
which the `max_size`th subunit is added. Each inner list contains the
forward and reverse rate constants (in that order) for the corresponding
assembly reaction, and each of these pairs must comprise solely
Parameter objects or solely numbers (never one of each). If Parameters
are passed, they will be used directly in the generated Rules. If
numbers are passed, Parameters will be created with automatically
generated names based on `subunit`, `site1`, `site2` and the chain sizes
and these parameters will be included at the end of the returned
component list.
Examples
--------
Assemble a three-membered chain by sequential addition of monomers,
with the same forward/reverse rates for monomer-monomer and monomer-dimer
interactions::
Model()
Monomer('Unit', ['p1', 'p2'])
assemble_chain_sequential(Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2)
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2'])
Monomer('Unit', ['p1', 'p2'])
>>> assemble_chain_sequential(Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('assemble_chain_sequential_Unit_2', Unit(p1=None, p2=None) + Unit(p1=None, p2=None) | Unit(p1=None, p2=1) % Unit(p1=1, p2=None), assemble_chain_sequential_Unit_2_kf, assemble_chain_sequential_Unit_2_kr),
Parameter('assemble_chain_sequential_Unit_2_kf', 0.0001),
Parameter('assemble_chain_sequential_Unit_2_kr', 0.1),
Rule('assemble_chain_sequential_Unit_3', Unit(p1=None, p2=None) + Unit(p1=None, p2=1) % Unit(p1=1, p2=None) | MatchOnce(Unit(p1=None, p2=1) % Unit(p1=1, p2=2) % Unit(p1=2, p2=None)), assemble_chain_sequential_Unit_3_kf, assemble_chain_sequential_Unit_3_kr),
Parameter('assemble_chain_sequential_Unit_3_kf', 0.0001),
Parameter('assemble_chain_sequential_Unit_3_kr', 0.1),
])
"""
return assemble_polymer_sequential(subunit, site1, site2, max_size, ktable,
closed=False)
def chain_species_base(base, basesite, subunit, site1, site2, size, comp=1):
"""
Return a MonomerPattern representing a chained species, chained to a base complex.
Parameters
----------
base : Monomer or MonomerPattern
The base complex to which the growing chain will be attached.
basesite : string
Name of the site on complex where first subunit binds.
subunit : Monomer or MonomerPattern
The subunit of which the chain is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next.
size : integer
The number of subunits in the chain.
comp : optional; a ComplexPattern to which the base molecule is attached.
Returns
-------
A ComplexPattern corresponding to the chain.
Notes
-----
Similar to pore_species, but never closes the chain.
Examples
--------
Get the ComplexPattern object representing a chain of size 4 bound to a base, which is itself bound to a complex:
Model()
Monomer('Base', ['b1', 'b2'])
Monomer('Unit', ['p1', 'p2'])
Monomer('Complex1', ['s1'])
Monomer('Complex2', ['s1', 's2'])
chain_tetramer = chain_species_base(Base(b1=1, b2=ANY), 'b1', Unit, 'p1', 'p2', 4, Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY))
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Unit', ['p1', 'p2'])
Monomer('Unit', ['p1', 'p2'])
>>> Monomer('Base', ['b1', 'b2'])
Monomer('Base', ['b1', 'b2'])
>>> Monomer('Complex1', ['s1'])
Monomer('Complex1', ['s1'])
>>> Monomer('Complex2', ['s1', 's2'])
Monomer('Complex2', ['s1', 's2'])
>>> chain_species_base(Base(b2=ANY), 'b1', Unit, 'p1', 'p2', 4, Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY))
MatchOnce(Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY) % Base(b1=1, b2=ANY) % Unit(p1=1, p2=2) % Unit(p1=2, p2=3) % Unit(p1=3, p2=4) % Unit(p1=4, p2=None))
"""
_verify_sites(base, basesite)
_verify_sites(subunit, site1, site2)
if size <= 0:
raise ValueError("size must be an integer greater than 0")
if comp == 1:
compbase = base({basesite: 1})
else:
compbase = comp % base({basesite: 1})
if size == 1:
chainlink = compbase % subunit({site1: 1, site2: None})
elif size == 2:
chainlink = compbase % subunit({site1: 1, site2: 2}) % \
subunit({site1: 2, site2: None})
else:
# build up a ComplexPattern, starting with a single subunit
chainbase = compbase
chainlink = chainbase % subunit({site1: 1, site2: 2})
for i in range(2, size):
chainlink %= subunit({site1: i, site2: i+1})
chainlink %= subunit({site1: size, site2: None})
chainlink.match_once = True
return chainlink
[docs]def assemble_chain_sequential_base(base, basesite, subunit, site1, site2, max_size, ktable, comp=1):
"""
Generate rules to assemble a homomeric chain sequentially onto a base complex (only the subunit creates repeating chain, not the base).
The chain species are created by sequential addition of `subunit` monomers.
The chain structure is defined by the `pore_species_base` macro.
Parameters
----------
base : Monomer or MonomerPattern
The base complex to which the chain is attached.
basesite : string
The name of the site on the complex to which chain attaches.
subunit : Monomer or MonomerPattern
The subunit of which the chain is composed.
site1, site2 : string
The names of the sites where one copy of `subunit` binds to the next; the first will also be the site where the first subunit binds the base.
max_size : integer
The maximum number of subunits in the chain.
ktable : list of lists of Parameters or numbers
Table of forward and reverse rate constants for the assembly steps. The
outer list must be of length `max_size` + 1, and the inner lists must
all be of length 2. In the outer list, the first element corresponds to
the first assembly step in which the complex binds the first subunit. The next corresponds to a bound subunit binding to form a
2-subunit complex, and the last element corresponds to the final step in
which the `max_size`th subunit is added. Each inner list contains the
forward and reverse rate constants (in that order) for the corresponding
assembly reaction, and each of these pairs must comprise solely
Parameter objects or solely numbers (never one of each). If Parameters
are passed, they will be used directly in the generated Rules. If
numbers are passed, Parameters will be created with automatically
generated names based on `subunit`, `site1`, `site2` and the chain sizes
and these parameters will be included at the end of the returned
component list.
comp : optional; a ComplexPattern to which the base molecule is attached.
Examples
--------
Assemble a three-membered chain by sequential addition of monomers to a base, which is in turn attached to a complex,
with the same forward/reverse rates for monomer-monomer and monomer-dimer
interactions::
Model()
Monomer('Base', ['b1', 'b2'])
Monomer('Unit', ['p1', 'p2'])
Monomer('Complex1', ['s1'])
Monomer('Complex2', ['s1', s2'])
assemble_chain_sequential(Base(b2=ANY), 'b1', Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2, Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY))
Execution::
>>> Model() # doctest:+ELLIPSIS
<Model '_interactive_' ...>
>>> Monomer('Base', ['b1', 'b2'])
Monomer('Base', ['b1', 'b2'])
>>> Monomer('Unit', ['p1', 'p2'])
Monomer('Unit', ['p1', 'p2'])
>>> Monomer('Complex1', ['s1'])
Monomer('Complex1', ['s1'])
>>> Monomer('Complex2', ['s1', 's2'])
Monomer('Complex2', ['s1', 's2'])
>>> assemble_chain_sequential_base(Base(b2=ANY), 'b1', Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2, Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY)) # doctest:+NORMALIZE_WHITESPACE
ComponentSet([
Rule('assemble_chain_sequential_base_Unit_2', Unit(p1=None, p2=None) + Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY) % Base(b1=1, b2=ANY) % Unit(p1=1, p2=None) | Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY) % Base(b1=1, b2=ANY) % Unit(p1=1, p2=2) % Unit(p1=2, p2=None), assemble_chain_sequential_base_Unit_2_kf, assemble_chain_sequential_base_Unit_2_kr),
Parameter('assemble_chain_sequential_base_Unit_2_kf', 0.0001),
Parameter('assemble_chain_sequential_base_Unit_2_kr', 0.1),
Rule('assemble_chain_sequential_base_Unit_3', Unit(p1=None, p2=None) + Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY) % Base(b1=1, b2=ANY) % Unit(p1=1, p2=2) % Unit(p1=2, p2=None) | MatchOnce(Complex1(s1=ANY) % Complex2(s1=ANY, s2=ANY) % Base(b1=1, b2=ANY) % Unit(p1=1, p2=2) % Unit(p1=2, p2=3) % Unit(p1=3, p2=None)), assemble_chain_sequential_base_Unit_3_kf, assemble_chain_sequential_base_Unit_3_kr),
Parameter('assemble_chain_sequential_base_Unit_3_kf', 0.0001),
Parameter('assemble_chain_sequential_base_Unit_3_kr', 0.1),
])
"""
if len(ktable) != max_size-1:
raise ValueError("len(ktable) must be equal to max_size-1")
def chain_rule_name(rule_expression, size):
react_p = rule_expression.reactant_pattern
monomer = react_p.complex_patterns[0].monomer_patterns[0].monomer
return '%s_%d' % (monomer.name, size)
components = ComponentSet()
s = subunit({site1:None, site2:None})
for size, klist in zip(range(2, max_size + 1), ktable):
chain_prev = chain_species_base(base, basesite, subunit, site1, site2, size - 1, comp)
chain_next = chain_species_base(base, basesite, subunit, site1, site2, size, comp)
name_func = functools.partial(chain_rule_name, size=size)
components |= _macro_rule('assemble_chain_sequential_base',
s + chain_prev | chain_next,
klist, ['kf', 'kr'],
name_func=name_func)
return components
if __name__ == "__main__":
import doctest
doctest.testmod()