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 using set_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. A reset() 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, use set_time() and set_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, as myokit.Name or myokit.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, as myokit.Name or myokit.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 or TimeSeriesProtocol 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 label pace is still registered (allowing later calls to add a protocol with set_protocol. Alternatively, if protocol={} then no pacing labels will be registered, and any subsequent calls to set_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 the path 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 using pickle 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, calling Simulation(model, path='my-sim.zip') will create a file called my-sim.zip in which the C extension is stored. To load, use Simulation.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 or myokit.TimeSeriesProtocol to use for the variable with binding pace. Atlernatively, a dictionary mapping binding labels to myokit.Protocol objects can be used to run with multiple protocols. Finally, can be None to run without a protocol.
sensitivities: An optional tuple (dependents, independents) where dependents is a list of variables or expressions to take derivatives of (y in dy/dx) and independents is a list of variables or expressions to calculate derivatives to (x in dy/dx). Each entry in dependents must be a myokit.Variable, a myokit.Name, a myokit.Derivative, or a string with either a fully qualified variable name or a dot() expression. Each entry in independents must be a myokit.Variable, a myokit.Name, a myokit.InitialValue, or a string with either a fully qualified variable name or an init() 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() or Simulation.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.

A typical use case is to inspect the calculation of derivatives after a crash, with:

sim.evaluate_derivatives(sim.crash_state(), sim.crash_inputs())

crash_log()¶

If the last call to Simulation.pre() or Simulation.run() resulted in an error, this will return the myokit.DataLog containing the data logged up until the crash.

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() or Simulation.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.

A typical use case is to inspect the calculation of derivatives after a crash, with:

sim.evaluate_derivatives(sim.crash_state(), sim.crash_inputs())

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 the state() as default, the new method uses default_state().

evaluate_derivatives(state=None, inputs=None)¶

Evaluates and returns the state derivatives for the given state and inputs.

If no state is given, the value returned by default_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 in inputs, the value 0 will be used, regardless of the current simulation. Similarly, if pacing values are not set in inputs, 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 by Simulation 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() or Simulation.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 using pre() 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 be 1 at time zero. Similarly, the sensitivity of any other state to that initial value at this time will be 0. Because pre() 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. The duration 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, a myokit.SimulationError may be raised.

To obtain feedback on the simulation progress, an object implementing the myokit.ProgressReporter interface can be passed in. passed in as progress. 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 method returns a myokit.DataLog dictionary that maps variable names to lists of logged values. The variables to log can be indicated using the log 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 or myokit.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 function myokit.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, the log_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 as apd_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)”.

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 if log_times is used. If both are None 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 if log_interval is used. If both are None 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, an apd_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, and None 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 sensitivities dy/dx, where the first index specifies the dependent variable y, and the second index specifies the independent variable x.

Finally, if APD calculation is enabled, the method returns a tuple (log, apds) or (log, sensitivities, apds) where apds is a myokit.DataLog with entries start and duration, 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 a Variable or a string containing a variable qname. The value 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 by times and values for the label pace.

This method is provided for backwards compatibility with older versions, please use set_protocol() and the TimeSeriesProtocol 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 a TimeSeriesProtocol for the given label.

To remove a previously set binding call this method with protocol = None. In this case, the value of any variables bound to label 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 using rel_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¶

class myokit.Sundials¶

Tests for Sundials/Sundials support.

static version()¶: Returns the detected Sundials version on this system, or None if no version of Sundials was found.

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). See myokit.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 or myokit.Derivative).
independents: A list of “independent variables” for sensitivity calculations, all specified as expressions (myokit.Name or myokit.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.

Myokit :: API and User Guide

Search

Single cell simulation¶

Sundials utility classes¶

C Model code¶