Initial-condition sensitivity

class myokit.ICSimulation(model, protocol=None)

Runs a forward-Euler based simulation and calculates the partial derivatives of the state vector with respect to the initial conditions.

The simulation is based on automatic differentiation implemented using a C++ data type that replaces a single scalar float with a float and a list of partial derivatives. Any operations on this pair update both the float and the set of derivatives. A normal simulation starts with a state y(tmin) and a right-hand side function (RHS) f(y) = dy/dt. It then integrates f(y) from tmin to tmax resulting in an output state y(tmax). In this simulation the data type of f is replaced by a (f, df/dy), where df/dy is the matrix of partial derivatives of f with respect to y. By integrating f from tmin to tmax we obtain the state at tmax. This can be seen as a function F(y(tmin)), that gives the state at tmax given y(tmin). By integrating df/dy the derivative of F to y(tmin) is obtained. This result allows the sensitivity of the system to its initial conditions to be evaluated.

N.B. The partial derivatives can not be calculated for the following functions: floor, ceil, abs, quotients and remainders. If these are encountered the resulting derivatives will be yielded as NaN. However, in many cases, these functions will only occur as part of a condition in an if statement, so the NaN‘s won’t propagate to the final result.

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 derivatives of the current state with respect to the initial state

When a simulation is created, the simulation time is set to 0 and the state is obtained from the given model. The initial derivatives matrix is an identity matrix of size (n, n), where n is the number of states in the model. After each call to run() the time, state and derivative variables are updated so that each successive call to run continues where the previous one left off. A reset() method is provided that will set the time back to 0, revert the current state to the default state and set the derivatives back to I.

The simulation provides two inputs a variable can bind to:

time
This variable contains the simulation time.
pace
This variable contains the current value of the pacing variable as given by the protocol passed to the Simulation.

No labeled variables are required.

block(log, derivatives)

Takes the output of a simulation (a simulation log and a list of derivatives) and combines it into a single DataBlock2d object.

Each entry in the log is converted to a 0d entry in the log. The calculated derivatives are stored as the 2d field derivatives.

default_state()

Returns the default state.

derivatives()

Return the partial derivatives of the current state with respect to the initial state.

reset()

Resets the simulation:

  • The time variable is set to 0
  • The state is set back to the default state
run(duration, log=None, log_interval=5, progress=None, msg='Running ICSimulation')

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 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_INTER, myokit.LOG_BOUND.
  • A list of qnames or variable objects
  • A myokit.DataLog obtained from a previous simulation. In this case, the newly logged data will be appended to the existing log.

For more details on the log argument, see the function myokit.prepare_log().

The method returns a myokit.DataLog and a 3d numpy array. In the returned array, the first axis represents the time, the second axis is a state x and the third is a state y such that the point (t, x, y) represents dx/dy(0) at time t. For example, if p is the array of derivatives, to get the derivative of state 0 with respect to the initial value of state 1, use p[:,0,1].

A log entry is created every time at least log_interval time units have passed.

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.

set_protocol(protocol=None)

Changes the pacing protocol used by this simulation.

set_step_size(dt=0.01)

Sets the step size used in the forward Euler solving routine.

state()

Returns the current state.

time()

Returns the current simulation time.