Notebook

Starting, stopping, pre-pacing, and loops¶

In this example we look at starting, stopping, and restarting simulations, we look at the pre() method, and finally use reset() to run several simulations in a loop. Although we use the myokit.Simulation class in these examples, the interface is shared by other simulation classes in Myokit.

In this notebook, we've used the model of the human ventricular AP by O'Hara et al., you can download a copy here.

Running repeated simulations¶

As before, we start by loading the model and its associated protocol, and creating a simulation:

In [1]:

import myokit
model, protocol, _ = myokit.load('models/c/ohara-2011.mmt')

print('Loaded model: ' + model.name())

Loaded model: ohara-2011

This model has an endo, epi, and mid-myocardial mode. We can switch between the modes using the variable cell.mode:

In [2]:

print(model.get('cell.mode').code())

mode = 1 : The type of cell. Endo = 0, Epi = 1, Mid = 2

Note that this is model-specific: calling the variable cell.mode was a choice made by the person implementing the model.

We can select the mode to use by editing the model code, or by using set_rhs():

In [3]:

# Select epicardial mode
model.get('cell.mode').set_rhs(1)

Next, we run a simulation and plot the results:

In [4]:

sim = myokit.Simulation(model, protocol)
log1 = sim.run(1000, log=['engine.time', 'membrane.V'])

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log1['engine.time'], log1['membrane.V'])
plt.show()

Now we run a second simulation, and see what happens:

In [5]:

log2 = sim.run(1000, log=['engine.time', 'membrane.V'])

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log2['engine.time'], log2['membrane.V'])
plt.show()

It looks very similar, but if look at the x-axis we can see that this second simulation starts at t=1000ms, so after the end of the first simulation.

We can stop and re-start simulations any time we like:

In [6]:

plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
for i in range(40):
    log = sim.run(10, log=['engine.time', 'membrane.V'])
    plt.plot(log['engine.time'], log['membrane.V'], lw=6)
plt.show()

Returning back to the two AP example, we can also overlay the APs to compare their shape. To do this, we subtract 1000ms from each logged time value in the second simulation. This could be done in a for-loop, but we can also use NumPy to create "array" objects from our logged simulation data. A quick way to convert our simulation results to NumPy arrays is to use the npview() method:

In [7]:

# Initially, logs are array.array objects, that can be appended to by the simulation:
print('Before conversion:')
print(type(log2['membrane.V']))
print('After conversion:')
log2 = log2.npview()
print(type(log2['membrane.V']))

Before conversion:
<class 'array.array'>
After conversion:
<class 'numpy.ndarray'>

(Note that we can now no longer use this DataLog for further simulation logging, as it has become a fixed-length object!)

Having the results as a NumPy array means we can now easily manipulate the data, for example:

In [8]:

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log1['engine.time'], log1['membrane.V'])
plt.plot(log2['engine.time'] - 900, log2['membrane.V'] + 10)
plt.show()

Or, more sensibly:

In [9]:

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log1['engine.time'], log1['membrane.V'])
plt.plot(log2['engine.time'] - 1000, log2['membrane.V'])
plt.show()

This shows us that the 2nd simulated AP differs from the first. We can zoom in on the start of the simulations to get a hint why:

In [10]:

plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log1['engine.time'], log1['membrane.V'], label='First AP')
plt.plot(log2['engine.time'] - 1000, log2['membrane.V'], label='Second AP')
plt.xlim(0, 11)
plt.ylim(-90, -80)
plt.legend()
plt.show()

This shows us that the simulation's initial conditions have changed: The first AP was simulated from the model's default initial state, while the second AP was simulated from the state at the end of the first simulation.

Using reset() to return to the default time and state¶

We can return to the original initial conditions using Simulation.reset():

In [11]:

# Reset to original time and initial conditions
sim.reset()

# Run another simulation
log3 = sim.run(1000, log=['engine.time', 'membrane.V'])

# Plot all three simulations, without any modifications
plt.figure(figsize=(15, 5))
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log1['engine.time'], log1['membrane.V'], label='Simulation 1')
plt.plot(log2['engine.time'], log2['membrane.V'], label='Simulation 2')
plt.plot(log3['engine.time'], log3['membrane.V'], '--', lw=3, label='Simulation 3')
plt.legend()
plt.show()

So reset() has done two things:

It has reset the time variable back to t=0ms
It has returned the model to its default state

The "default state"¶

There are three aspects involved in the example above:

The simulation's current time -- which we can get and set with Simulation.time() and Simulation.set_time().
The simulation's current state -- which we can get and set with Simulation.state() and Simulation.set_state().
The simulation's default state -- which we can get and set with Simulation.default_state() and Simulation.set_default_state().

So a call to reset() is the same as sim.set_time(0) followed by sim.set_state(sim.default_state()).

We've also seen that a call to run() starts from the current time and state, and updates the current time and state.

In the next parts of the example we show how these concepts come into play in everyday simulation work, and introduce a third command: pre().

Converging to a "steady state"¶

Let's have another look at the two APs we simulated earlier:

In [12]:

plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(log1['engine.time'], log1['membrane.V'], label='AP 1')
plt.plot(log2['engine.time'] - 1000, log2['membrane.V'], label='AP 2')
plt.xlim(0, 400)
plt.legend()
plt.show()

By looking closely at the voltage at the start of the AP, we already saw that the initial states for both simulations had changed. We can look in some more detail using Simulation.state() and Model.format_state():

In [13]:

# Reset to t=0 and to the default state
sim.reset()

# Store the initial state
state1 = sim.state()

# Run a simulation
sim.run(1000, log=[])

# Store the state after 1000ms
state2 = sim.state()

# Use model.format_state to compare the two states
print('Variable         state1    state2')
print('-'*51)
print(model.format_state(state1, state2))

Variable         state1    state2
---------------------------------------------------
membrane.V        = -87.0     -8.80076744333505303e+01
sodium.Na_i       = 7.0        7.00067350382111098e+00
sodium.Na_ss      = 7.0        7.00072502364536664e+00
potassium.K_i     = 145.0      1.44993522780129808e+02
potassium.K_ss    = 145.0      1.44993501338738383e+02
calcium.Ca_i      = 0.0001     6.02071865129394192e-05
calcium.Ca_ss     = 0.0001     5.99721909601639493e-05
calcium.Ca_nsr    = 1.2        1.29190735117844402e+00
calcium.Ca_jsr    = 1.2        1.24303497912602223e+00
ina.m             = 0.0        7.33986039713834838e-03
ina.hf            = 1.0        6.98307076731258136e-01
ina.hs            = 1.0        6.98298322628229551e-01
ina.j             = 1.0        6.98249666153142545e-01
ina.hsp           = 1.0        4.55216899109574469e-01
ina.jp            = 1.0        6.98218030111901311e-01
inal.m            = 0.0        1.88055497054156475e-04
inal.h            = 1.0        5.05749429238630421e-01
inal.hp           = 1.0        4.09751355530827310e-01
ito.a             = 0.0        1.00070589506202857e-03
ito.if            = 1.0        9.99554681363342801e-01
ito.is            = 1.0        9.99556687584296943e-01
ito.ap            = 0.0        5.09886566672653331e-04
ito.ifp           = 1.0        9.99554681636059650e-01
ito.isp           = 1.0        9.99559102540451550e-01
ical.d            = 0.0        2.33758500380497421e-09
ical.ff           = 1.0        9.99999990889891555e-01
ical.fs           = 1.0        9.40349032343499758e-01
ical.fcaf         = 1.0        9.99999990889923418e-01
ical.fcas         = 1.0        9.99804558020380041e-01
ical.jca          = 1.0        9.99981793534284180e-01
ical.nca          = 0.0        7.18412993133261154e-04
ical.ffp          = 1.0        9.99999990871068389e-01
ical.fcafp        = 1.0        9.99999990886676571e-01
ikr.xf            = 0.0        7.07409707227045945e-06
ikr.xs            = 0.0        4.48403563353674650e-01
iks.x1            = 0.0        1.30166545730574484e-01
iks.x2            = 0.0        1.92694086685758407e-04
ik1.x             = 1.0        9.96758045172253970e-01
ryr.Jrelnp        = 0.0        9.23485029466150725e-08
ryr.Jrelp         = 0.0        1.15310075899567968e-07
camk.CaMK_trapped = 0.0        2.60827759773024246e-03

This shows us all the model's state variables, their original values, and their values after the first simulation. The first thing to notice is that original values shown here are all nice rounded numbers: it seems these were all set manually by the person who coded up the model!

Looking a bit closer at the state variables (and perhaps opening up the model definition file, and looking at the original paper!), we can see that the top rows all contain concentrations, which seem to have changed after the first beat.

We can run a longer simulation and see what's going on:

In [14]:

# Select a few concentrations for plotting:
concentrations = [
    'sodium.Na_i',
    #'sodium.Na_ss',
    'potassium.K_i',
    #'potassium.K_ss',
    'calcium.Ca_i',
    #'calcium.Ca_ss',
    #'calcium.Ca_nsr',
    #'calcium.Ca_jsr',
]

# Simulate 5 beats (each 1000ms long)
sim.reset()
log = sim.run(5 * 1000, log=['engine.time'] + concentrations)

# Show how the concentrations change over time
n = len(concentrations)
plt.figure(figsize=(15, 9))
for i, variable in enumerate(concentrations):
    plt.subplot(n, 1, 1 + i)
    plt.plot(log['engine.time'], log[variable], label=variable)
    plt.legend()
plt.xlabel('Time (ms)')
plt.show()

The traces for the concentrations are almost, but not quite periodical. They make roughly the same pattern at each beat, but they don't quite end up where they started, and so the system changes subtly from beat to beat.

We can run an even longer simulation to see what happens over time:

In [15]:

# Simulate 50 beats (each 1000ms long)
sim.reset()
log = sim.run(50 * 1000, log=['engine.time'] + concentrations)

# Show how the concentrations change over time
n = len(concentrations)
plt.figure(figsize=(15, 9))
for i, variable in enumerate(concentrations):
    plt.subplot(n, 1, 1 + i)
    plt.plot(log['engine.time'], log[variable], label=variable)
    plt.legend()
plt.xlabel('Time (ms)')
plt.show()

Uh oh! The sodium concentration, which had been going down in the first 10 beats, has now started going up. And the potassium concentration is moving steadily downwards. Let's run an even longer simulation...

In [16]:

# Simulate 1000 beats (each 1000ms long)
sim.reset()
log = sim.run(1000 * 1000, log=['engine.time'] + concentrations)

# Show how the concentrations change over time
n = len(concentrations)
plt.figure(figsize=(10, 9))
for i, variable in enumerate(concentrations):
    plt.subplot(n, 1, 1 + i)
    plt.plot(log['engine.time'], log[variable], label=variable)
    plt.legend()
plt.xlabel('Time (ms)')
plt.show()

This looks a bit better! It seems the concentrations for sodium and potassium have stabilised, although for calcium it's a bit harder to see what's going on.

At this scale, the periodic behaviour is very hard to see in the plot, so we can consider omitting it altogether by only logging a single point in each beat:

In [17]:

# Simulate 1000 beats (each 1000ms long), but only log one beat
sim.reset()
log = sim.run(
    1000 * 1000, 
    log=['engine.time'] + concentrations, 
    log_interval=1000)

# Show how the concentrations change over time
n = len(concentrations)
plt.figure(figsize=(10, 9))
for i, variable in enumerate(concentrations):
    plt.subplot(n, 1, 1 + i)
    plt.plot(log['engine.time'], log[variable], label=variable)
    plt.legend()
plt.xlabel('Time (ms)')
plt.show()

So it looks like the concentrations have stabilised at this point: they still vary during each beat, but they return to the same value at the end of every beat.

(More complicated things can also happen. For example you can arrive in a situation where the system alternates between two similar beats. This is called "alternans" and has been shown to occur in real cells as well as models. It's usually related to variations in the Calcium concentration.)

In cardiac electrophysiology, we say that the model is now in steady state. If a mathematician is listening, it's better to say the system has converged to a periodic orbit: it still varies all the time, but the system's orbit through state space has converged to a periodic one. We saw above that the system has quite a lot of state variables, so this orbit is hard to visualise, but we can at least plot "slices" of it, to see that this is the case:

In [18]:

# Simulate a few beats, starting at the end of the previous simulation
plt.figure(figsize=(12, 12))
plt.xlabel('[Na]i')
plt.ylabel('[K]i')
for i in range(5):
    log = sim.run(
        1000,
        log=['engine.time'] + concentrations,
        log_interval=5
    )
    plt.plot(log['sodium.Na_i'], log['potassium.K_i'], 'x-', alpha=0.5)
plt.show()

The first thing you might notice is that the "orbit" shown here is quite complex! Action potential models are not very simple systems.

Another thing you can clearly see in this graph is the straight line at the top: this corresponds to the rhythmic pacing that is applied to this model, and the subsequent upstroke.

Finally, note that the lines from the 10 simulations are almost identical, making it hard to distinguish them. However, the line still looks somewhat blurry, showing that there is still some adaptation going on from beat to beat.

There's a lot of questions you might want to ask at this point:

Has the system really reached a "steady state"? Or will it start diverging from the periodic orbit at some point?
Does every model have such a steady state?
And if so, is there a single steady state or will we end up in different places depending on the initial conditions?
Are there clever ways to find these "steady states" or "periodic orbits", or do we just need to simulate for a long time?
Most importantly, is this something that happens only in models, or is it a feature of real cells?

Some notable studies that touch on the 2nd question include Endresen et al., 2000 and Hund et al., 2001, which both focus on calculating the membrane potential from the concentrations directly, and Livshitz & Rudy, 2009, which looks at the role of the stimulus current and conservation of charge.

A paper by Jacquemet, 2007 shows that, in some models at least, there is more than one "steady state".

However, in the vast majority of published simulation studies, modellers have taken a practical approach and simply "pre-paced" the model to something approximating a steady state, before running further experiments.

Pre-pacing AP models¶

Because pre-pacing is so common, Myokit's simulations have a method pre() for pre-pacing.

Like the run() method, pre() takes a duration as its input, for example:

In [19]:

# Create a simulation
sim = myokit.Simulation(model, protocol)

In [20]:

# Pre-pace for 500s
sim.pre(500 * 1000)

There are three major differences from run():

The pre() method does not do any logging (making this a lot faster), and does not return a DataLog.
The pre() does not affect the simulation time.
The pre() method updates both the state and the default state.

We can check the second two differences: The time is still 0 (ms):

In [21]:

print(sim.time())

The default state has been updated, and is no longer set to the user-defined original values:

In [22]:

print(model.format_state(sim.default_state(), sim.state()))

membrane.V        = -8.79366724794197410e+01    -8.79366724794197410e+01
sodium.Na_i       =  7.86126498244010019e+00     7.86126498244010019e+00
sodium.Na_ss      =  7.86134716986204030e+00     7.86134716986204030e+00
potassium.K_i     =  1.43994952092648504e+02     1.43994952092648504e+02
potassium.K_ss    =  1.43994919690563904e+02     1.43994919690563904e+02
calcium.Ca_i      =  7.64959399791582224e-05     7.64959399791582224e-05
calcium.Ca_ss     =  7.54965657440926968e-05     7.54965657440926968e-05
calcium.Ca_nsr    =  2.00009369544246196e+00     2.00009369544246196e+00
calcium.Ca_jsr    =  1.94170772144389936e+00     1.94170772144389936e+00
ina.m             =  7.39245449482547282e-03     7.39245449482547282e-03
ina.hf            =  6.95843566016024928e-01     6.95843566016024928e-01
ina.hs            =  6.95832944414410348e-01     6.95832944414410348e-01
ina.j             =  6.95771580881882401e-01     6.95771580881882401e-01
ina.hsp           =  4.52317914455469816e-01     4.52317914455469816e-01
ina.jp            =  6.95728616753768514e-01     6.95728616753768514e-01
inal.m            =  1.90608712351722818e-04     1.90608712351722818e-04
inal.h            =  5.01629103365055284e-01     5.01629103365055284e-01
inal.hp           =  2.75411135383717920e-01     2.75411135383717920e-01
ito.a             =  1.00550739918049215e-03     1.00550739918049215e-03
ito.if            =  9.99549112477026491e-01     9.99549112477026491e-01
ito.is            =  9.99542939129281938e-01     9.99542939129281938e-01
ito.ap            =  5.12334269909103779e-04     5.12334269909103779e-04
ito.ifp           =  9.99549112809742790e-01     9.99549112809742790e-01
ito.isp           =  9.99547953162419178e-01     9.99547953162419178e-01
ical.d            =  2.37715649270646627e-09     2.37715649270646627e-09
ical.ff           =  9.99999990712996500e-01     9.99999990712996500e-01
ical.fs           =  9.21983529780828093e-01     9.21983529780828093e-01
ical.fcaf         =  9.99999990713036691e-01     9.99999990713036691e-01
ical.fcas         =  9.99885710699513597e-01     9.99885710699513597e-01
ical.jca          =  9.99984993142850942e-01     9.99984993142850942e-01
ical.nca          =  1.74881882071773148e-03     1.74881882071773148e-03
ical.ffp          =  9.99999990711129660e-01     9.99999990711129660e-01
ical.fcafp        =  9.99999990711243125e-01     9.99999990711243125e-01
ikr.xf            =  8.09453698088356724e-06     8.09453698088356724e-06
ikr.xs            =  4.27668970548617744e-01     4.27668970548617744e-01
iks.x1            =  2.42746952623205281e-01     2.42746952623205281e-01
iks.x2            =  1.94239320459573233e-04     1.94239320459573233e-04
ik1.x             =  9.96776710047558079e-01     9.96776710047558079e-01
ryr.Jrelnp        =  4.59571517646687238e-07     4.59571517646687238e-07
ryr.Jrelp         =  5.74361409907950119e-07     5.74361409907950119e-07
camk.CaMK_trapped =  1.60382047012359863e-02     1.60382047012359863e-02

Because the default state has changed, subsequent calls to reset() will reset the simulation to the last state reached by pre().

Running simulations in a loop¶

Using pre(), reset(), and run() we now have the tools to explore the effects of varying model parameters.

In the next example we:

pre-pace the model to a quasi-steady state
modify the model's IKr conductance, set by the variable ikr.gKr
run another series of unlogged beats, this time with run()
plot an AP
reset the simulation to the sitatuation at the end of step 1, and repeat starting from step 2

In [23]:

# Create a simulation
sim = myokit.Simulation(model, protocol)

# Pre-pace 500 beats
sim.pre(500 * 1000)

# Original gKr value
gKr_original = 0.046

# Create a figure
plt.figure(figsize=(12, 6))
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')

# Explore parameter variation
levels = [1.0, 0.9, 0.8, 0.7, 0.6, 0.5]
for level in levels:
    # Reset to state after initial pre-pacing
    sim.reset()
    
    # Modify IKr conductance
    sim.set_constant('ikr.gKr', level * gKr_original)
   
    # Pre-pace for the adapted level
    sim.run(50 * 1000, log=[])
    sim.set_time(0)

    # Calculate percentage block, for display
    block = round(100 * (1 - level))
    
    # Plot AP
    log = sim.run(450, log=['engine.time', 'membrane.V'])
    plt.plot(log['engine.time'], log['membrane.V'], label=str(block) + '% IKr block')

# Show results
plt.legend()
plt.show()

Note that, while we're pre-pacing now, in this example we don't actually check whether the pre-pacing has had the desired effect. To do this, we could simulate another beat (after sim.run(50 * 1000, log=[])) and e.g.

compare the result of sim.state() before and after the beat, or
log and compare two action potentials.

If you write a good method to do this, you could consider adding it to Myokit!