Single cell simulation¶
- class myokit.Simulation(model, protocol=None, sensitivities=None, path=None)¶
Runs single cell simulations using the CVODES solver (see [1]); CVODES uses an implicit multi-step method to achieve high accuracy and stability with adaptive step sizes.
The model passed to the simulation is cloned and stored internally, so changes to the original model object will not affect the simulation. A protocol can be passed in as
protocol
or set later usingset_protocol()
.Simulations maintain an internal state consisting of
the current simulation time
the current state
the default state
(optional) the current and default state’s sensitivities with respect to the selected independent variables.
When a simulation is created, the simulation time is set to 0 and both the current and the default state are copied from the model. After each call to
Simulation.run()
the time variable and current state are updated, so that each successive call to run continues where the previous simulation left off. Areset()
method is provided that will set the time back to 0 and revert the current state to the default state. To change the time or state manually, useset_time()
andset_state()
.A pre-pacing method
pre()
is provided that doesn’t affect the simulation time but will update the current and the default state. This allows you to pre-pace, run a simulation, reset to the pre-paced state, run another simulation etc.Sensitivities
Sensitivities of model variables with respect to parameters (any unbound model variable defined as a literal constant) or initial values of states can be calculated using the
sensitivities
argument. If set, this should be a tuple(ys, xs)
containing the “dependent” and “independent” expressions used to make up the calculated sensitivities (partial derivatives)dy/dx
.The “dependent variables” must refer to state variables, time-derivatives of state variables, or intermediary variables. These can be specified as
myokit.Variable
objects, asmyokit.Name
ormyokit.Derivative
expressions, or as strings e.g."ina.INa"
or"dot(membrane.V)"
.The “independent variables” in sensitivities must refer to either literal variables (variables with no dependencies) or initial values of state variables. These can be specified as
myokit.Variable
objects, asmyokit.Name
ormyokit.InitialValue
expressions, or as strings e.g."ikr.gKr"
or"init(membrane.V)"
.Bound variables and labels
The simulation provides four inputs a model variable can be bound to:
time
This input provides the simulation time.
pace
This input provides the current value of the pacing variable. This is determined using the protocol passed into the Simulation.
evaluations
This input provides the number of rhs evaluations used at each point in time and can be used to gain some insight into the solver’s behavior.
realtime
This input provides the elapsed system time at each logged point.
No variable labels are required for this simulation type.
Multiple protocols or no protocol
This simulation supports pacing with more than one protocol. To this end, pass in a dictionary mapping pacing labels (bindings) to
Protocol
orTimeSeriesProtocol
objects, e.g.protocol={'pace_1': protocol_1, 'pace_2': protocol_2}
.For backwards compatibility, if no protocol is set and
protocol=None
, then the pacing labelpace
is still registered (allowing later calls to add a protocol withset_protocol
. Alternatively, ifprotocol={}
then no pacing labels will be registered, and any subsequent calls toset_protocol
will fail.Storing and loading simulation objects
There are two ways to store Simulation objects to the file system: 1. using serialisation (the
pickle
module), and 2. using thepath
constructor argument.Each time a simulation is created, a C module is compiled, imported, and deleted in the background. This means that part of a
Simulation
object is a reference to an imported C extension, and cannot be serialized. When usingpickle
to serialize the simulation, this compiled part is not stored. Instead, when the pickled simulation is read from disk the compilation is repeated. Unpickling simulations also restores their state: variables such as model state, simulation time, and tolerance are preserved when pickling.As an alternative to serialisation, Simulations can be created with a
path
variable that specifies a location where the generated module can be stored. For example, callingSimulation(model, path='my-sim.zip')
will create a file calledmy-sim.zip
in which the C extension is stored. To load, useSimulation.from_path('my-sim.zip')
. Unlike pickling, simulations stored and loaded this way do not maintain state: variables such as model state, simulation time, and tolerance are not preserved with this method. Finally, note that (again unlike pickled simulations), the generated zip files are highly platform dependent: a zip file generated on one machine may not work on another.Arguments
model
The model to simulate
protocol
A
myokit.Protocol
ormyokit.TimeSeriesProtocol
to use for the variable with bindingpace
. Atlernatively, a dictionary mapping binding labels tomyokit.Protocol
objects can be used to run with multiple protocols. Finally, can beNone
to run without a protocol.sensitivities
An optional tuple
(dependents, independents)
wheredependents
is a list of variables or expressions to take derivatives of (y
indy/dx
) andindependents
is a list of variables or expressions to calculate derivatives to (x
indy/dx
). Each entry independents
must be amyokit.Variable
, amyokit.Name
, amyokit.Derivative
, or a string with either a fully qualified variable name or adot()
expression. Each entry inindependents
must be amyokit.Variable
, amyokit.Name
, amyokit.InitialValue
, or a string with either a fully qualified variable name or aninit()
expression.path
An optional path used to load or store compiled simulation objects. See “Storing and loading simulation objects”, above.
References
[1] SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. Hindmarsh, Brown, Woodward, et al. (2005) ACM Transactions on Mathematical Software.
- crash_inputs()¶
If the last call to
Simulation.pre()
orSimulation.run()
resulted in an error, this will return the last “inputs” (time, pace, etc) reached during that simulation.Will return
None
if no simulation was run or the simulation did not result in an error.
- crash_state()¶
If the last call to
Simulation.pre()
orSimulation.run()
resulted in an error, this will return the last state reached during that simulation.Will return
None
if no simulation was run or the simulation did not result in an error.
- default_state()¶
Returns the default state.
- default_state_sensitivities()¶
Returns the default sensitivities with respect to state variables, or
None
if sensitivities are not enabled.
- eval_derivatives(y=None, pacing=None)¶
Deprecated alias of
evaluate_derivatives()
.Note that while
eval_derivatives()
used thestate()
as default, the new method usesdefault_state()
.
- evaluate_derivatives(state=None, inputs=None)¶
Evaluates and returns the state derivatives for the given
state
andinputs
.If no
state
is given, the value returned bydefault_state()
is used.An optional dict
inputs
can be passed in to set the values of time and other inputs (e.g. pacing values). If “time” is not set ininputs
, the value 0 will be used, regardless of the current simulation. Similarly, if pacing values are not set ininputs
, 0 will be used regardless of the protocols or the current time.If any variables have been changed with set_constant their new value will be used.
- static from_path(path)¶
Loads a simulation from the zip file in
path
.The file at
path
must have been created bySimulation
and on the same platform (e.g. the same operating system, processor architecture, Myokit version etc.).Note that the simulation state (e.g. time, model state, solver tolerance etc. is not restored).
Returns a
Simulation
object.
- last_number_of_evaluations()¶
Returns the number of rhs evaluations performed by the solver during the last simulation.
- last_number_of_steps()¶
Returns the number of steps taken by the solver during the last simulation.
- last_state()¶
If the last call to
Simulation.pre()
orSimulation.run()
resulted in an error, this will return the last state reached during that simulation.Will return
None
if no simulation was run or the simulation did not result in an error.
- pre(duration, progress=None, msg='Pre-pacing simulation')¶
This method can be used to perform an unlogged simulation, typically to pre-pace to a (semi-)stable orbit.
After running this method
The simulation time is not affected
The current state and the default state are updated to the final state reached in the simulation.
The current state-sensitivities and default state-sensitivities are updated to the final state reached in the simulation, except for state-sensitivities with respect to initial conditions, which are reset (see below).
Calls to
reset()
after usingpre()
will set the current state to this new default state.When calculating sensitivities with respect to the initial value of some state variable with index
i
, the sensitivity of that state w.r.t. the initial value will be1
at time zero. Similarly, the sensitivity of any other state to that initial value at this time will be0
. Becausepre()
is used to set a new state for time zero, this method will set the current and default state sensitivities for initial value sensitivities to 0 or 1 as above, instead of using the values reached in the simulation.The number of time units to simulate can be set with
duration
. Theduration
argument cannot be negative, and special care needs to be taken when very small (positive) values are used. For some very short (but non-zero) durations, amyokit.SimulationError
may be raised.To obtain feedback on the simulation progress, an object implementing the
myokit.ProgressReporter
interface can be passed in. passed in asprogress
. An optional description of the current simulation to use in the ProgressReporter can be passed in as msg.
- reset()¶
Resets the simulation:
The time variable is set to 0
The state is set to the default state
The state sensitivities are set to the default state sensitivities
- run(duration, log=None, log_interval=None, log_times=None, sensitivities=None, apd_variable=None, apd_threshold=None, progress=None, msg='Running simulation')¶
Runs a simulation and returns the logged results. Running a simulation has the following effects:
The internal state is updated to the last state in the simulation.
- The simulation’s time variable is updated to reflect the time
elapsed during the simulation.
The number of time units to simulate can be set with
duration
. Theduration
argument cannot be negative, and special care needs to be taken when very small (positive) values are used. For some very short (but non-zero) durations, amyokit.SimulationError
may be raised.The method returns a
myokit.DataLog
dictionary that maps variable names to lists of logged values. The variables to log can be indicated using thelog
argument. There are several options for its value:None
(default), to log all states.An integer flag or a combination of flags. Options:
myokit.LOG_NONE
,myokit.LOG_STATE
,myokit.LOG_BOUND
,myokit.LOG_INTER
,myokit.LOG_DERIV
ormyokit.LOG_ALL
.A sequence of variable names. To log derivatives, use “dot(membrane.V)”.
A
myokit.DataLog
object. In this case, the new data will be appended to the existing log.
For detailed information about the
log
argument, see the functionmyokit.prepare_log()
.By default, every step the solver takes is logged. This is usually advantageous, since more points are added exactly at the times the system gets more interesting. However, if equidistant points are required a
log_interval
can be set. Alternatively, thelog_times
argument can be used to specify logging times directly.To get action potential duration (APD) measurements, the simulation can be run with threshold crossing detection. To enable this, pass in a state variable as
apd_variable
and a threshold value asapd_threshold
. Please note the APD calculation implemented here uses a fixed threshold, and so differs from the often used dynamical thresholds such as “90% of max(V) - min(V)”.To obtain feedback on the simulation progress, an object implementing the
myokit.ProgressReporter
interface can be passed in. passed in asprogress
. An optional description of the current simulation to use in the ProgressReporter can be passed in asmsg
.Arguments:
duration
The time to simulate.
log
The variables to log (see detailed description above).
log_interval
An optional fixed size log interval. Must be
None
iflog_times
is used. If both areNone
every step is logged.log_times
An optional sequence (e.g. a list or a numpy array) of pre-determined logging times. Must be
None
iflog_interval
is used. If both areNone
every step is logged.sensitivities
An optional list-of-lists to append the calculated sensitivities to.
apd_variable
An optional
myokit.Variable
or fully-qualified variable name to use in APD calculations. If set, anapd_threshold
must also be specified.apd_threshold
An optional (fixed) threshold to use in APD calculations. Must be set if and
apd_variable
is set, andNone
if not.progress
An optional
myokit.ProgressReporter
used to obtain feedback about simulation progress.msg
An optional message to pass to any progress reporter.
By default, this method returns a
myokit.DataLog
containing the logged variables.However, if sensitivity calculations are enabled a tuple is returned, where the first entry is the
myokit.DataLog
and the second entry is a matrix of sensitivitiesdy/dx
, where the first index specifies the dependent variabley
, and the second index specifies the independent variablex
.Finally, if APD calculation is enabled, the method returns a tuple
(log, apds)
or(log, sensitivities, apds)
whereapds
is amyokit.DataLog
with entriesstart
andduration
, representing the start and duration of all measured APDs.
- set_constant(var, value)¶
Changes a model constant. Only literal constants (constants not dependent on any other variable) can be changed.
The constant
var
can be given as aVariable
or a string containing a variable qname. Thevalue
should be given as a float.
- set_default_state(state)¶
Change the default state to
state
.
- set_fixed_form_protocol(times=None, values=None)¶
Sets a
TimeSeriesProtocol
specified bytimes
andvalues
for the labelpace
.This method is provided for backwards compatibility with older versions, please use
set_protocol()
and theTimeSeriesProtocol
class instead.
- set_max_step_size(dtmax=None)¶
Sets a maximum step size. To let the solver pick any step size it likes use
dtmax = None
.
- set_min_step_size(dtmin=None)¶
Sets a minimum step size. To let the solver pick any step size it likes use
dtmin = None
.
- set_protocol(protocol, label='pace')¶
Set an event-based pacing
Protocol
or aTimeSeriesProtocol
for the givenlabel
.To remove a previously set binding call this method with
protocol = None
. In this case, the value of any variables bound tolabel
will be set to 0.The label must be one of the pacing labels set in the constructor.
- set_state(state)¶
Sets the current state.
- set_time(time=0)¶
Sets the current simulation time.
- set_tolerance(abs_tol=1e-06, rel_tol=0.0001)¶
Sets the solver tolerances. Absolute tolerance is set using
abs_tol
, relative tolerance usingrel_tol
. For more information on these values, see the Sundials CVODES documentation.
- state()¶
Returns the current state.
- time()¶
Returns the current simulation time.
Sundials utility classes¶
C Model code¶
The standard simulation uses an ansi-C model template, generated by the
myokit.CModel
class.
- class myokit.CModel(model, pacing_labels, sensitivities)¶
Generates ansi-C code for a model.
Example:
model = myokit.parse_model('example') cmodel = CModel(model) print(cmodel.code)
Arguments:
model
A
myokit.Model
.pacing_labels
A list of
str
variable labels, each corresponding to a paced variable.sensitivities
Either
None
or a tuple(dependents, independents)
. Seemyokit.Simulation
for details.
The following properties are all public for easy access. But note that they do not interact with the compiled header so changing them will have little effect.
model
The myokit.Model this CModel was compiled for.
code
The generated model code. Can be placed in a header file or incorporated directly.
Sensitivities:
has_sensitivities
True if this model was created with sensitivity calculations enabled.
dependents
A list of “dependent variables” for sensitivity calculations, all specified as expressions (
myokit.Name
ormyokit.Derivative
).independents
A list of “independent variables” for sensitivity calculations, all specified as expressions (
myokit.Name
ormyokit.InitialValue
).
Constants (and parameters) are all stored inside ordered dicts mapping variable objects onto equations. Each dict is stored in a solvable order.
parameters
Constants used as “independent variable” in sensitivity calculations.
parameter_derived
Constants that depend on parameters (and possibly on literals too).
literals
Constants that do not depend on any other constants (and are not parameters).
literal_derived
Constants that depend on literals, but not on parameters.