librelist archives

« back to archive

Re: Delivery Status Notification (Failure)

Re: Delivery Status Notification (Failure)

From:
John Bachman
Date:
2013-08-11 @ 20:25
Hi Danny,

We're currently working on a feature that lets you write parameters as
functions of time, other parameters, or observables as SymPy
expressions. This would let you, for example, write functional rate
laws for certain reactions. However, this feature is not quite ready
for prime time (though it's close).

However, it sounds like you might not need this function to do what
you're trying to do. I thought I'd give you two examples for typical
use cases for involving "non-autonomous" simulation. In the first, you
simply vary a parameter, like an initial condition or a rate
parameter, and use it to run a series of simulations (e.g., for a
parameter sweep). The second use case handles the case where you want
to the output of one simulation (e.g., a steady-state value) to be fed
in as the input of the second. This can be used for complex
experimental perturbations, step-function inputs, etc. Here's the code
for both:

Example 1: Parameter sweep (plot at http://imgur.com/rLtZ6Mk)

from pylab import *
from pysb import kappa
from pysb.examples.hello_pysb import model

ion()
# Ligand concentrations for the "sweep"
ligand_amounts = linspace(10, 300, 30)
# We'll store the fraction of receptor bound
# for each concentration in here
frac_bound_values = []
# Loop over the ligand concentrations
for ligand_amount in ligand_amounts:
    # Set the initial condition for ligand
    model.parameters['L_0'].value = ligand_amount
    # Run kappa
    simdata = kappa.run_simulation(model, time=500)
    # Save the "steady-state" (final) value for the fraction bound
    frac_bound_values.append(simdata['LR'][-1] /
                             model.parameters['R_0'].value)
# Plot fraction bound vs. ligand concentration
plot(ligand_amounts, frac_bound_values, marker='o')

Example 2: Create an explicit perturbation in an experiment

(see plot at http://imgur.com/Qc6BfpF)

Imagine we want to simulate a type of binding experiment in which we
first allow binding to occur and then we pull out the bound complexes
and watch them dissociate (e.g., to measure the off-rate). To simulate
this we'll first allow L and R to bind and allow LR to reach steady
state, then reset the simulation with no free L and R, and the amount
of LR set as the steady-state amount from the first phase of the
experiment. This two-phase approach is a decent way to look at the
effect of simple step-function inputs and other perturbations, which
may be what you're looking to do.

It's also worth noting that the Kappa team has developed a language
for describing experimental perturbations for Kappa simulations. I
personally have not used this feature but you can read more about it
in the KaSim manual (in the /man subdirectory of your KaSim install).
As a quick-and-dirty way of allowing this feature to be used from
within PySB, there's a "perturbation" keyword argument to the
run_kasim function in pysb.kappa that simply appends any perturbation
language commands at the end of the generated Kappa file.

In the meantime, here's an example of cobbling together a simple
experimental perturbation using PySB:

from pylab import *
from pysb import *
from pysb import kappa
from pysb.util import alias_model_components
from pysb.examples.hello_pysb import model

# Calling this function allows us to refer to components from the hello_pysb
# model as Python variables in this namespace
alias_model_components()

# Reset the off-rate value to get a more interesting plot
kr.value = 1e-2
# Run the binding phase of the experiment
bind_tmax = 200
bind_sim = kappa.run_simulation(model, time=bind_tmax)
# Get the final (steady-state) value for LR
steady_state_LR = bind_sim['LR'][-1]
# Reset the initial conditions of the model for the unbinding phase
L_0.value = 0
R_0.value = 0
# Add a new, nonzero initial condition for LR
Initial(L(s=1) % R(s=1), Parameter('LR_0', steady_state_LR))
# Run the unbinding phase of the experiment
unbind_sim = kappa.run_simulation(model, time=200)
# Link the time and LR vectors together:
# Note that since the last point of the binding phase is identical to the
# first point of the unbinding phase, we make sure we don't double-count it.
# Note also that we offset the time vector by the duration of the first
# experiment:
time_vec = concatenate((bind_sim['time'][0:-1], unbind_sim['time'] + bind_tmax))
LR_vec = concatenate((bind_sim['LR'][0:-1], unbind_sim['LR']))
# Plot
ion()
plot(time_vec, LR_vec)
# Add a vertical line to mark the transition from binding to unbinding
ylim((0, max(LR_vec)))
plot([200, 200], [0, max(LR_vec)], color='r')

Re: Delivery Status Notification (Failure)

From:
John Bachman
Date:
2013-08-11 @ 20:25
Here's the code for the examples (attached)

On Sun, Aug 11, 2013 at 4:25 PM, John Bachman <bachmanjohn@gmail.com> wrote:
> Hi Danny,
>
> We're currently working on a feature that lets you write parameters as
> functions of time, other parameters, or observables as SymPy
> expressions. This would let you, for example, write functional rate
> laws for certain reactions. However, this feature is not quite ready
> for prime time (though it's close).
>
> However, it sounds like you might not need this function to do what
> you're trying to do. I thought I'd give you two examples for typical
> use cases for involving "non-autonomous" simulation. In the first, you
> simply vary a parameter, like an initial condition or a rate
> parameter, and use it to run a series of simulations (e.g., for a
> parameter sweep). The second use case handles the case where you want
> to the output of one simulation (e.g., a steady-state value) to be fed
> in as the input of the second. This can be used for complex
> experimental perturbations, step-function inputs, etc. Here's the code
> for both:
>
> Example 1: Parameter sweep (plot at http://imgur.com/rLtZ6Mk)
>
> from pylab import *
> from pysb import kappa
> from pysb.examples.hello_pysb import model
>
> ion()
> # Ligand concentrations for the "sweep"
> ligand_amounts = linspace(10, 300, 30)
> # We'll store the fraction of receptor bound
> # for each concentration in here
> frac_bound_values = []
> # Loop over the ligand concentrations
> for ligand_amount in ligand_amounts:
>     # Set the initial condition for ligand
>     model.parameters['L_0'].value = ligand_amount
>     # Run kappa
>     simdata = kappa.run_simulation(model, time=500)
>     # Save the "steady-state" (final) value for the fraction bound
>     frac_bound_values.append(simdata['LR'][-1] /
>                              model.parameters['R_0'].value)
> # Plot fraction bound vs. ligand concentration
> plot(ligand_amounts, frac_bound_values, marker='o')
>
> Example 2: Create an explicit perturbation in an experiment
>
> (see plot at http://imgur.com/Qc6BfpF)
>
> Imagine we want to simulate a type of binding experiment in which we
> first allow binding to occur and then we pull out the bound complexes
> and watch them dissociate (e.g., to measure the off-rate). To simulate
> this we'll first allow L and R to bind and allow LR to reach steady
> state, then reset the simulation with no free L and R, and the amount
> of LR set as the steady-state amount from the first phase of the
> experiment. This two-phase approach is a decent way to look at the
> effect of simple step-function inputs and other perturbations, which
> may be what you're looking to do.
>
> It's also worth noting that the Kappa team has developed a language
> for describing experimental perturbations for Kappa simulations. I
> personally have not used this feature but you can read more about it
> in the KaSim manual (in the /man subdirectory of your KaSim install).
> As a quick-and-dirty way of allowing this feature to be used from
> within PySB, there's a "perturbation" keyword argument to the
> run_kasim function in pysb.kappa that simply appends any perturbation
> language commands at the end of the generated Kappa file.
>
> In the meantime, here's an example of cobbling together a simple
> experimental perturbation using PySB:
>
> from pylab import *
> from pysb import *
> from pysb import kappa
> from pysb.util import alias_model_components
> from pysb.examples.hello_pysb import model
>
> # Calling this function allows us to refer to components from the hello_pysb
> # model as Python variables in this namespace
> alias_model_components()
>
> # Reset the off-rate value to get a more interesting plot
> kr.value = 1e-2
> # Run the binding phase of the experiment
> bind_tmax = 200
> bind_sim = kappa.run_simulation(model, time=bind_tmax)
> # Get the final (steady-state) value for LR
> steady_state_LR = bind_sim['LR'][-1]
> # Reset the initial conditions of the model for the unbinding phase
> L_0.value = 0
> R_0.value = 0
> # Add a new, nonzero initial condition for LR
> Initial(L(s=1) % R(s=1), Parameter('LR_0', steady_state_LR))
> # Run the unbinding phase of the experiment
> unbind_sim = kappa.run_simulation(model, time=200)
> # Link the time and LR vectors together:
> # Note that since the last point of the binding phase is identical to the
> # first point of the unbinding phase, we make sure we don't double-count it.
> # Note also that we offset the time vector by the duration of the first
> # experiment:
> time_vec = concatenate((bind_sim['time'][0:-1], unbind_sim['time'] + bind_tmax))
> LR_vec = concatenate((bind_sim['LR'][0:-1], unbind_sim['LR']))
> # Plot
> ion()
> plot(time_vec, LR_vec)
> # Add a vertical line to mark the transition from binding to unbinding
> ylim((0, max(LR_vec)))
> plot([200, 200], [0, max(LR_vec)], color='r')