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 behaviour.
realtime
This input provides the elapsed system time at each logged point.

No variable labels are required for this simulation type.

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 serialised. When using pickle to serialise 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
An optional myokit.Protocol to use as input for variables bound to pace.
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.

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)

Evaluates and returns the state derivatives.

The state to evaluate for can be given as y. If no state is given the current simulation state is 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 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.

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)”.

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.

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 set 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 indice specifies the dependent variable y, and the second indice 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)

Allows you to manually set the default state.

set_fixed_form_protocol(times=None, values=None)

Configures this simulation to run with a predetermined protocol instead of the usual event-based mechanism.

A 1D time-series should be given as input. During the simulation, the value of the pacing variable will be determined by linearly interpolating between the two nearest points in the series. If the simulation time is outside the bounds of the time-series, the first or last value in the series will be used.

Setting a predetermined protocol clears any previously set (event-based or pre-determined) protocol. To clear all protocols, call this method with times=None. When a simulation is run without any protocol, the value of any variables bound to pace will be set to 0.

Arguments:

times
A non-decreasing array of times. If any times appear more than once, only the value at the highest index will be used.
values
An array of values for the pacing variable. Must have the same size as times.
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=None)

Sets the pacing Protocol used by this simulation.

To run without pacing call this method with protocol = None. In this case, the value of any variables bound to pace will be set to 0.

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, 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.
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.

Legacy code

Myokit is migrating from a CVODE to a CVODES based simulation class. The older code is still available for now as the myokit.LegacySimulation.

myokit.LegacySimulation

alias of myokit._sim.cvodesim.Simulation