Skip to content
Permalink
Browse files

Drug binding macro (#387)

* Drug binding macro that mimics experimental settings when a perturbation
is added after a reaction is initiated
  • Loading branch information...
ortega2247 authored and alubbock committed Jan 15, 2019
1 parent b8e036c commit c267f1ec33942a401a6b2d4683f2e195ebd164fd
Showing with 162 additions and 1 deletion.
  1. +162 −1 pysb/macros.py
@@ -47,7 +47,7 @@
'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']
'bind_complex', 'bind_table_complex', 'drug_binding']

# Suppress ModelExistsWarnings in our doctests.
_pysb_doctest_suppress_modelexistswarning = True
@@ -845,6 +845,167 @@ def bind_table_complex(bindtable, row_site, col_site, m1=None, m2=None, kf=None)
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_' (monomers: 0, rules: 0, parameters: 0, expressions: 0, compartments: 0) at ...>
>>> 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


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_' (monomers: 0, rules: 0, parameters: 0, expressions: 0, compartments: 0) at ...>
>>> 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', (t > 10)*kf_drug_substrate),
Expression('kr_expr_drug_substrate', (t > 10)*kr_drug_substrate),
])
"""
_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), (time_obs > t_action) * k1)
kr_expr = Expression('kr_expr_{0}_{1}'.format(drug_monomer_name,
substrate_monomer_name), (time_obs > t_action) * k2)
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
# =========

0 comments on commit c267f1e

Please sign in to comment.
You can’t perform that action at this time.