Markov Models¶
The module myokit.lib.markov
contains functions for working with
Markov models of ion channel currents.
- The class
LinearModel
can be used to extract a Markov model from a Myokit model.
- By fixing the voltage and all other states that aren’t part of the
Markov model a linear system is created.
- Fast simulations can then be performed using the
AnalyticalSimulation
class, which is particularly useful when trying to estimate a Markov model’s parameters from a piecewise-linear voltage protocol (e.g. a normal step protocol).- Discrete, stochastic simulations can be performed using the
DiscreteSimulation
class.
Analytical and discrete simulation¶
- class myokit.lib.markov.LinearModel(model, states, parameters=None, current=None, vm=None)¶
Represents a linear Markov model of an ion channel extracted from a
myokit.Model
.The class assumes the markov model can be written as:
dot(x) = A(V,p) * x I = B(V,p) * x
where
V
is the membrane potential,p
is a set of parameters andA
andB
are the matrices that relate the statex
to its derivativedot(x)
and a currentI
.A
andB
can contain non-linear functions, but should be simple scalar matrices when evaluated for a fixedV
andp
. For example, the current equationI = g * O * (V - E)
would havep = [g, E]
, so that for fixedp
andV
this would resolve toI = (g * (V - E)) * O
, such thatg * (V - E)
is a constant that can be included inB
.The model variables to treat as parameter are specified by the user when the model is created. Any other variables, for example state variables such as intercellular calcium or constants such as temperature, are fixed when the markov model is created and can no longer be changed. Initial values written as expressions are evaluated when the model is made.
To create a
Markov
, pass in amyokit.Model
and select a list of states. All other states will be fixed at their current value and an attempt will be made to write the remaining state equations as linear combinations of the states. If this is not possible, aValueError
is raised. The membrane potential must be indicated using the labelmembrane_potential
or by passing it in asvm
.The current variable is optional, if no current is specified by the user the relation
I = B * x
is dropped and noB
is calculated.Example:
import myokit import myokit.lib.markov as markov # Load a model from disk model = myokit.load_model('some-model.mmt') # Select the relevant states and parameters states = [ 'ina.C1', 'ina.C2', 'ina.O', ... ] parameters = [ 'ina.p1', 'ina.p2', ... ] current = 'ina.INa' # Extract a markov model mm = markov.LinearModel(model, states, parameters, current) # Get the matrices A and B such that dot(x) = A * x and I = B * x # where ``x`` is the state vector and ``I`` is the current. A, B = mm.matrices(membrane_potential=-40) print(A)
Alternatively, a LinearModel can be constructed from a single component using the method
from_component()
:import myokit import myokit.lib.markov as markov # Load a model from disk model = myokit.load_model('some-model.mmt') # Extract a markov model mm = markov.LinearModel.from_component(model.get('ina'))
Arguments:
model
The model to work with.
states
An ordered list of state variables (or state variable names) from
model
. All remaining state variables will be frozen in place. Each state’s derivative must be a linear combination of the other states.parameters
A list of parameters to maintain in their symbolic form.
current
The markov model’s current variable. The current must be a linear combination of the states (for example
g * (V - E) * (O1 + O2)
) whereO1
andO2
are states. If no current variable is specifiedNone
can be used instead.vm
The variable indicating membrane potential. If set to
None
(default) the method will search for a variable with the labelmembrane_potential
.
- current()¶
Returns the name of the current variable used by this model, or None if no current variable was specified.
- default_membrane_potential()¶
Returns this markov model’s default membrane potential value.
- default_parameters()¶
Returns this markov model’s default parameter values
- default_state()¶
Returns this markov model’s default state values.
- static from_component(component, states=None, parameters=None, current=None, vm=None)¶
Creates a Markov model from a component, using the following rules:
Every state in the component is a state in the Markov model
Every unnested constant in the component is a parameter
The component should contain exactly one unnested intermediary variable whose value depends on the model states, this will be used as the current variable.
The model contains a variable labeled “membrane_potential”.
Any of the automatically set variables can be overridden using the keyword arguments
states
,parameters
,current
andvm
.The parameters, if determined automatically, will be specified in alphabetical order (using a natural sort).
- matrices(membrane_potential=None, parameters=None)¶
For a given value of the
membrane_potential
and a list of values for theparameters
, this method calculates and returns the matricesA
andB
such that:dot(x) = A * x I = B * x
where
x
is the state vector andI
is the current.Arguments:
membrane_potential
The value to use for the membrane potential, or
None
to use the value from the originalmyokit.Model
.parameters
The values to use for the parameters, given in the order they were originally specified in (if the model was created using
from_component()
, this will be alphabetical order).
- membrane_potential()¶
Returns the name of the membrane potential variable used by this model.
- parameters()¶
Returns the names of the parameter variables used by this model.
- rates(membrane_potential=None, parameters=None)¶
For a given value of the
membrane_potential
and a list of values for theparameters
, this method calculates and returns an ordered list of tuples(i, j, rij)
such thatrij
is a non-zero transition rate from thei
-th state to thej
-th state.Arguments:
membrane_potential
The value to use for the membrane potential, or
None
to use the value from the originalmyokit.Model
.parameters
The values to use for the parameters, given in the order they were originally specified in (if the model was created using
from_component()
, this will be alphabetical order).
- states()¶
Returns the names of the state variables used by this model.
- steady_state(membrane_potential=None, parameters=None)¶
Analytically determines a steady state solution for this Markov model.
membrane_potential
The value to use for the membrane potential, or
None
to use the value from the originalmyokit.Model
.parameters
The values to use for the parameters, given in the order they were originally specified in (if the model was created using
from_component()
, this will be alphabetical order).
- class myokit.lib.markov.LinearModelError(message)¶
Raised for issues with constructing or using a
LinearModel
.
- class myokit.lib.markov.AnalyticalSimulation(model, protocol=None)¶
Analytically evaluates a
LinearModel
’s state over a given set of points in time.Solutions are calculated for the “law of large numbers” case, i.e. without stochastic behavior. The solution algorithm is based on eigenvalue decomposition.
Each simulation object maintains an internal state consisting of
The current simulation time
The current state
The default state
When a simulation is created, the simulation time is set to zero and both the current and default state are initialized using the
LinearModel
. After each call torun()
the time and current state are updated, so that each successive call to run continues where the previous simulation left off.A
protocol
can be used to set the membrane potential during the simulation, or the membrane potential can be adjusted manually between runs.Example:
import myokit import myokit.lib.markov as markov # Create a linear markov model m = myokit.load_model('clancy-1999.mmt') m = markov.LinearModel.from_component(m.get('ina')) # Create an analytical simulation object s = markov.AnalyticalSimulation(m) # Run a simulation s.set_membrane_potential(-30) d = s.run(10) # Show the results import matplotlib.pyplot as plt plt.figure() plt.subplot(211) for state in m.states(): plt.plot(d.time(), d[state], label=state) plt.legend(loc='center right') plt.subplot(212) plt.plot(d.time(), d[m.current()]) plt.show()
- current(state)¶
Calculates the current for a given state.
- default_state()¶
Returns the default state used by this simulation.
- membrane_potential()¶
Returns the currently set membrane potential.
- parameters()¶
Returns the currently set parameter values.
- pre(duration)¶
Performs an unlogged simulation for
duration
time units and uses the final state as the new default state.After the simulation:
The simulation time is not affected
The current state and the default state are updated to the final state reached in the simulation.
Calls to
reset()
after usingpre()
will set the current state to this new default state.
- reset()¶
Resets the simulation:
The time variable is set to zero.
The state is set to the default state.
- run(duration, log=None, log_interval=0.01, log_times=None)¶
Runs a simulation for
duration
time units.After the simulation:
The simulation time will be increased by
duration
time units.The simulation state will be updated to the last reached state.
Arguments:
duration
The number of time units to simulate.
log
A log from a previous run can be passed in, in which case the results will be appended to this one.
log_interval
The time between logged points.
log_times
A pre-defined sequence of times to log at. If set,
log_interval
will be ignored.
Returns a
myokit.DataLog
with the simulation results.
- set_constant(variable, value)¶
Updates a single parameter to a new value.
- set_default_state(state)¶
Changes this simulation’s default state.
- set_membrane_potential(v)¶
Changes the membrane potential used in this simulation.
- set_parameters(parameters)¶
Changes the parameter values used in this simulation.
- set_state(state)¶
Changes the initial state used by in this simulation.
- solve(times)¶
Evaluates and returns the states at the given times.
In contrast to
run()
, this method simply evaluates the states (and current) at the given times, using the last known settings for the state and membrane potential. It does not use a protocol and does not take into account the simulation time. After running this method, the state and simulation time are not updated.Arguments:
times
A series of times, where each time must be some
t >= 0
.
For models with a current variable, this method returns a tuple
(state, current)
wherestate
is a matrix of shape(len(states), len(times))
andcurrent
is a vector of lengthlen(times)
.For models without a current variable, only
state
is returned.
- state()¶
Returns the initial state used by this simulation.
- class myokit.lib.markov.DiscreteSimulation(model, protocol=None, nchannels=100)¶
Performs stochastic simulations of a
LinearModel
’s behavior for a finite number of channels.Simulations are run using the “Direct method” proposed by Gillespie [1].
Each simulation object maintains an internal state consisting of
The current simulation time
The current state
The default state
When a simulation is created, the simulation time is set to zero and both the current and default state are initialized using the
LinearModel
. After each call torun()
the time and current state are updated, so that each successive call to run continues where the previous simulation left off.A
protocol
can be used to set the membrane potential during the simulation, or the membrane potential can be adjusted manually between runs.Example:
import myokit import myokit.lib.markov as markov # Create linear markov model m = myokit.load_model('clancy-1999.mmt') m = markov.LinearModel.from_component(m.get('ina')) # Run discrete simulation s = markov.DiscreteSimulation(m, nchannels=1000) s.set_membrane_potential(-30) d = s.run(10) import matplotlib.pyplot as plt plt.figure() for state in m.states(): plt.step(d.time(), d[state], label=state) plt.legend() plt.show()
References
- [1] Gillespie (1976) A General Method for Numerically Simulating the
stochastic time evolution of coupled chemical reactions The Journal of Computational Physics, 22, 403-434.
Arguments:
model
A
LinearModel
.nchannels
The number of channels to simulate.
- default_state()¶
Returns the default simulation state.
- discretize_state(x)¶
Converts a list of fractional state occupancies to a list of channel counts.
Arguments:
x
A fractional state where
sum(x) == 1
.
Returns a discretized state
y
wheresum(y) = nchannels
.
- membrane_potential()¶
Returns the current membrane potential.
- number_of_channels()¶
Returns the number of channels used in this simulation.
- parameters()¶
Returns the current parameter values.
- pre(duration)¶
Performs an unlogged simulation for
duration
time units and uses the final state as the new default state.After the simulation:
The simulation time is not affected
The current state and the default state are updated to the final state reached in the simulation.
Calls to
reset()
after usingpre()
will set the current state to this new default state.
- reset()¶
Resets the simulation:
The time variable is set to zero.
The state is set to the default state.
- run(duration, log=None)¶
Runs a simulation for
duration
time units.After the simulation:
The simulation time will be increased by
duration
time units.The simulation state will be updated to the last reached state.
Arguments:
duration
The number of time units to simulate.
log
A log from a previous run can be passed in, in which case the results will be appended to this one.
Returns a
myokit.DataLog
with the simulation results.
- set_constant(variable, value)¶
Updates a single parameter to a new value.
- set_default_state(state)¶
Changes the default state used in the simulation.
- set_membrane_potential(v)¶
Changes the membrane potential used in this simulation.
- set_parameters(parameters)¶
Changes the parameter values used in this simulation.
- set_state(state)¶
Changes the current state used in the simulation (i.e. the number of channels in every markov model state).
- state()¶
Returns the current simulation state.
Finding Markov models¶
- myokit.lib.markov.convert_markov_models_to_compact_form(model)¶
Scans a
myokit.Model
for Markov models, and ensures they contain one state that’s not evaluated as an ODE, but as1 - sum(x[i])
, where the sum is over all other statesx[i]
.Arguments:
model
The
myokit.Model
to scan.
Returns an updated
myokit.Model
.
- myokit.lib.markov.convert_markov_models_to_full_ode_form(model)¶
Scans a
myokit.Model
for Markov models, and ensures they are written in a form where every Markov state is evaluated as an ODE.Arguments:
model
The
myokit.Model
to scan.
Returns an updated
myokit.Model
.
- myokit.lib.markov.find_markov_models(model)¶
Searches a
myokit.Model
for groups of states that constitute a Markov model.Returns a list of lists, where the inner lists are groups of variables that form a Markov model together.
Note that this method performs a shallow check of the equation shapes, and does not perform any simplification or rewriting to see if the expressions can be made to fit a Markov form.
Arguments:
model
The
myokit.Model
to search.
Deprecated¶
The following class was used in previous versions of Myokit (before 1.22.0). It now exists only as an interface to the newer classes. The MarkovModel class will be removed in future versions.
- class myokit.lib.markov.MarkovModel(model, states, parameters=None, current=None, vm=None)¶
Deprecated: This class has been replaced by the classes
LinearModel
andAnalyticalSimulation
. Please update your code to use these classes instead. This class will be removed in future versions of Myokit.- static from_component(component, states=None, parameters=None, current=None, vm=None)¶
Creates and returns an
AnalyticalSimulation
using aLinearModel
based on a Myokit model component.