diff --git a/ft-example/.gitignore b/ft-example/.gitignore
new file mode 100644
index 0000000..7f47fdc
--- /dev/null
+++ b/ft-example/.gitignore
@@ -0,0 +1 @@
+combined.zip
diff --git a/ft-example/ft-example.py b/ft-example/ft-example.py
new file mode 100755
index 0000000..48b47cb
--- /dev/null
+++ b/ft-example/ft-example.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+#
+# FiberTissueSimulation example
+#
+import sys
+import myokit
+
+debug = 'debug' in sys.argv
+
+# Load models
+mf = myokit.load_model('stewart-2009.mmt')
+mt = myokit.load_model('tentusscher-2006-lhopital.mmt')
+
+# Create a protocol
+p = myokit.pacing.blocktrain(period=1000, duration=2)
+
+# Set up a simulation
+s = myokit.FiberTissueSimulation(
+    mf,
+    mt,
+    p,
+    ncells_fiber=(32, 2),
+    ncells_tissue=(32, 32),
+    nx_paced=5,
+    g_fiber=(3, 1),
+    g_tissue=(3, 1),
+    g_fiber_tissue=3,
+)
+s.set_step_size(0.0005)
+
+# Debugging? Then log all variables
+if debug:
+    logf = logt = myokit.LOG_STATE + myokit.LOG_BOUND
+else:
+    logf = ['engine.time', 'membrane.V', 'calcium.Ca_i']
+    logt = ['engine.time', 'membrane.V', 'calcium.Cai']
+
+# Run (using a progress printer to see how far we got)
+pr = myokit.ProgressPrinter()
+df, dt = s.run(50, logf=logf, logt=logt, log_interval=0.5, progress=pr)
+if debug:
+    print('Sim OK!')
+    sys.exit(1)
+
+# Convert the dictionary logs to structured 2d arrays
+bf = df.block2d()
+bt = dt.block2d()
+
+# Create a combined data block
+# First define a map that indicates equivalent variables. Format:
+# new_name : (modelf name, modelt name, padding value)
+varmap = {
+    'V' : ('membrane.V', 'membrane.V', -100),
+    'Cai': ('calcium.Ca_i', 'calcium.Cai', 0),
+}
+b = myokit.DataBlock2d.combine(bf, bt, map2d=varmap, pos1=(0, 15))
+
+# Store the data block
+fname = 'combined.zip'
+print(f'Storing in {fname}')
+b.save(fname)
+
+print('Done!')
+print('View with: python3 -m myokit block combined.zip')
+
diff --git a/ft-example/stewart-2009.mmt b/ft-example/stewart-2009.mmt
new file mode 100644
index 0000000..fa18924
--- /dev/null
+++ b/ft-example/stewart-2009.mmt
@@ -0,0 +1,577 @@
+[[model]]
+name: stewart-2009
+display_name: Stewart et al., 2009
+version: 20211027
+mmt_authors: Michael Clerx
+desc: """
+    Model of the human Purinkje action potential by Stewart et al. [1], based
+    on the ventricular endocardial model by Ten Tusscher et al. [2].
+
+    This implementation was created from the CellML code made available with
+    the original publication [2]. Several "cosmetic" changes were made without
+    affecting the equations, which were tested against the CellML from PMR and
+    matched to within rounding error precision.
+    - The state variables were ordered as in the publication
+    - The components were ordered as in the publication, and were annotated 
+      with links to the relevant document section.
+    In addition, incorrect units in F, R, Cm, gCaL, Vc, Vss, and Vsr were
+    corrected.
+
+    [1] Mathematical model of the electrical action potential of Purkinje fibre
+        cells.
+        Philip Stewart, Oleg V. Aslanidi, Denis Noble, Penelope J. Noble, Mark
+        R. Boyett, and Henggui Zhang, 2009
+        Philosophical Transactions of the Royal Society, 367, 2225-2255.
+        https://doi.org/10.1098/rsta.2008.0283
+
+    [2] Alternans and spiral breakup in a human ventricular tissue model.
+        K.H.W.J. ten Tusscher, A.V. Panfilov, Sep 2006,
+        American Journal of Physiology, Heart and Circulatory Physiology,
+        291 3, H1088-1100.
+        https://doi.org/10.1152/ajpheart.00109.2006
+
+    [3] https://models.physiomeproject.org/exposure/38cf8387b0707f0ef6947f009710aeb5
+        Retrieved on 2021-10-27
+
+    An extract from the CellML file's meta data follows below:
+    ---------------------------------------------------------------------------
+    Penny Noble
+    DPAG, University of Oxford
+
+    Model Status
+    This CellML ersion of the model has been checked in COR and OpenCell. The
+    units are consistent and the model runs to recreate the published results.
+
+    Model Structure
+    In the paper described here, Philip Stewart and colleagues present a
+    details of their newly developed model for the human Purkinje cell
+    including validation against experimental data. Ionic mechanisms underlying
+    the heterogeneity between the Purkinje fibre and ventricular action
+    potentials in humans and other species were analysed. The newly developed
+    Purkinje fibre cell model adds a new member to the family of human cardiac
+    cell models developed previously for the sino-atrial node, atrial and
+    ventricular cells, which can be incorporated into an anatomical model of
+    human heart with details of its electrophysiological heterogeneity and
+    anatomical complexity.
+    """
+# Initial values
+membrane.V    = -7.12864384994752527e+01
+calcium.Ca_i  =  1.02363913704157998e-04
+calcium.Ca_SR =  3.14149868138687083e+00
+calcium.Ca_ss =  3.81250245527617196e-04
+sodium.Na_i   =  8.80505373054986329e+00
+potassium.K_i =  1.36773426842998674e+02
+ina.m         =  2.83293473203283658e-02
+ina.h         =  2.51110676365191021e-01
+ina.j         =  2.65680422809985550e-01
+ikr.Xr1       =  1.34628009070351228e-02
+ikr.Xr2       =  3.32697379577155727e-01
+iks.Xs        =  8.40127007165385360e-03
+ito.r         =  8.84290083740064937e-04
+ito.s         =  9.68749754333663282e-01
+ical.d        =  2.16290842830038621e-04
+ical.f1       =  9.68969771502977917e-01
+ical.f2       =  9.96328488888527430e-01
+ical.fCass    =  9.99935288766206076e-01
+if.y          =  3.84052487230444398e-02
+irel.R_prime  =  9.84086493577585375e-01
+
+#
+# Simulation variables
+#
+[engine]
+time = 0 [ms]
+    in [ms]
+    bind time
+pace = 0
+    bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+dot(V) = -(i_ion + i_stim + i_diff)
+    in [mV]
+    label membrane_potential
+i_ion = sodium.INa_tot + potassium.IK_tot + calcium.ICa_tot
+    in [A/F]
+i_diff = 0 [A/F]
+    in [A/F]
+    bind diffusion_current
+i_stim = engine.pace * amplitude
+    in [A/F]
+amplitude = -30 [A/F]
+    in [A/F]
+
+#
+# (a) Inward rectifier current, IK1
+# Appendix page 2242
+#
+[ik1]
+use membrane.V
+inf = 1 / (1 + exp(0.1 [1/mV] * (V + 75.44 [mV])))
+g_K1 = 0.065 [mS/uF]
+    in [mS/uF]
+i_K1 = g_K1 * inf * (V - 8 [mV] - nernst.E_K)
+    in [A/F]
+
+#
+# (b) Transient outward current, ITo
+# Appendix page 2242
+#
+[ito]
+use membrane.V
+dot(r) = (r_inf - r) / tau_r
+    r_inf = 1 / (1 + exp((20 [mV] - V) / 13 [mV]))
+    tau_r = 10.45 [ms] * exp(-(V + 40 [mV])^2 / 1800 [mV^2]) + 7.3 [ms]
+        in [ms]
+dot(s) = (s_inf - s) / tau_s
+    s_inf = 1 / (1 + exp((V + 27 [mV]) / 13 [mV]))
+    tau_s = 85 [ms] * exp(-(V + 25 [mV])^2 / 320 [mV^2]) + 5 [ms] / (1 + exp((V - 40 [mV]) / 5 [mV])) + 42 [ms]
+        in [ms]
+g_to = 0.08184 [mS/uF]
+    in [mS/uF]
+i_to = g_to * r * s * (V - nernst.E_K)
+    in [A/F]
+
+#
+# (c) Sustained current, Isus
+# Appendix page 2243
+#
+[isus]
+use membrane.V
+a = 1 / (1 + exp((5 [mV] - V) / 17 [mV]))
+g_sus = 0.0227 [mS/uF]
+    in [mS/uF]
+i_sus = g_sus * a * (V - nernst.E_K)
+    in [A/F]
+
+#
+# (d) Hyperpolarization-activated current, If
+# Appendix page 2243
+# Described in detail on page 2229
+#
+[if]
+use membrane.V
+dot(y) = (y_inf - y) / tau_y
+    alpha_y = 1 [1/ms] * exp(-2.9 - 0.04 [1/mV] * V)
+        in [1/ms]
+    beta_y = 1 [1/ms] * exp(3.6 + 0.11 [1/mV] * V)
+        in [1/ms]
+    tau_y = 4000 / (alpha_y + beta_y)
+        in [ms]
+    y_inf = 1 / (1 + exp((V + 80.6 [mV]) / 6.8 [mV]))
+g_f_K = 0.0234346 [mS/uF]
+    in [mS/uF]
+g_f_Na = 0.0145654 [mS/uF]
+    in [mS/uF]
+i_f_K = y * g_f_K * (V - nernst.E_K)
+    in [A/F]
+i_f_Na = y * g_f_Na * (V - nernst.E_Na)
+    in [A/F]
+i_f = i_f_Na + i_f_K
+    in [A/F]
+
+#
+# (e) Fast sodium current, INa
+# Appendix page 2243
+#
+[ina]
+use membrane.V
+dot(m) = (inf - m) / tau
+    alpha = 1 / (1 + exp((-60 [mV] - V) / 5 [mV]))
+    beta = 0.1 / (1 + exp((V + 35 [mV]) / 5 [mV])) + 0.1 / (1 + exp((V - 50 [mV]) / 200 [mV]))
+    tau = 1 [ms] * alpha * beta
+        in [ms]
+    inf = 1 / (1 + exp((-56.86 [mV] - V) / 9.03 [mV]))^2
+dot(h) = (inf - h) / tau
+    alpha = if(V < -40 [mV], 0.057 [1/ms] * exp(-(V + 80 [mV]) / 6.8 [mV]), 0 [1/ms])
+        in [1/ms]
+    beta = if(V < -40 [mV],
+              2.7 [1/ms] * exp(0.079 [1/mV] * V) + 310000 [1/ms] * exp(0.3485 [1/mV] * V),
+              0.77 [1/ms] / (0.13 * (1 + exp((V + 10.66 [mV]) / -11.1 [mV]))))
+        in [1/ms]
+    tau = 1 / (alpha + beta)
+        in [ms]
+    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
+dot(j) = (inf - j) / tau
+    alpha = if(V < -40 [mV],
+               (-25428 [1/ms] * exp(0.2444 [1/mV] * V) - 6.948e-6 [1/ms] * exp(-0.04391 [1/mV] * V)) * (V + 37.78 [mV]) / 1 [mV] / (1 + exp(0.311 [1/mV] * (V + 79.23 [mV]))),
+               0 [1/ms])
+        in [1/ms]
+    beta = if(V < -40 [mV],
+              0.02424 [1/ms] * exp(-0.01052 [1/mV] * V) / (1 + exp(-0.1378 [1/mV] * (V + 40.14 [mV]))),
+              0.6 [1/ms] * exp(0.057 [1/mV] * V) / (1 + exp(-0.1 [1/mV] * (V + 32 [mV]))))
+        in [1/ms]
+    tau = 1 / (alpha + beta)
+        in [ms]
+    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
+g_Na = 130.5744 [mS/uF]
+    in [mS/uF]
+i_Na = g_Na * m^3 * h * j * (V - nernst.E_Na)
+    in [A/F]
+
+#
+# (f) L-type calcium current, ICaL
+# Appendix page 2244
+#
+[ical]
+use membrane.V
+use phys.FRT, phys.FFRT
+use calcium.Ca_ss, extra.Ca_o
+dot(d) = (d_inf - d) / tau_d
+    alpha_d = 1.4 / (1 + exp((-35 [mV] - V) / 13 [mV])) + 0.25
+    beta_d = 1.4 / (1 + exp((V + 5 [mV]) / 5 [mV]))
+    gamma_d = 1 [ms] / (1 + exp((50 [mV] - V) / 20 [mV]))
+        in [ms]
+    tau_d = 1 [ms] * alpha_d * beta_d + gamma_d
+        in [ms]
+    d_inf = 1 / (1 + exp((-8 [mV] - V) / 7.5 [mV]))
+dot(f1) = (f_inf - f1) / tau_f
+    f_inf = 1 / (1 + exp((V + 20 [mV]) / 7 [mV]))
+    tau_f = 1102.5 [ms] * exp(-(V + 27 [mV])^2 / 225 [mV^2]) + 200 [ms] / (1 + exp((13 [mV] - V) / 10 [mV])) + 180 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV])) + 20 [ms]
+        in [ms]
+dot(f2) = (f2_inf - f2) / tau_f2
+    f2_inf = 0.67 / (1 + exp((V + 35 [mV]) / 7 [mV])) + 0.33
+    tau_f2 = 562 [ms] * exp(-(V + 27 [mV])^2 / 240 [mV^2]) + 31 [ms] / (1 + exp((25 [mV] - V) / 10 [mV])) + 80 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV]))
+        in [ms]
+dot(fCass) = (fCass_inf - fCass) / tau_fCass
+    fCass_inf = 0.6 / (1 + (calcium.Ca_ss / 0.05 [mM])^2) + 0.4
+    tau_fCass = 80 [ms] / (1 + (calcium.Ca_ss / 0.05 [mM])^2) + 2 [ms]
+        in [ms]
+g_CaL = 0.0398 [L/F/s]
+    in [L/F/s]
+i_CaL = if(abs(V - 15 [mV]) < 1e-6 [mV], a * (b - Ca_o) / (2 * FRT), numer / denom)
+    in [A/F]
+    a = g_CaL * d * f1 * f2 * fCass * 4 * FFRT
+        in [1/mM/ms]
+    b = 0.25 * Ca_ss
+        in [mM]
+    numer = a * (V - 15 [mV]) * (b * exp((V - 15 [mV]) * 2 * FRT) - Ca_o)
+        in [A/F]
+    denom = (exp(2 * (V - 15 [mV]) * FRT) - 1)
+
+#
+# (g) Slow delayed rectifier current, IKs
+# Appendix page 2245
+# Also called the "slow time-dependent potassium current"
+#
+[iks]
+use membrane.V
+dot(Xs) = (inf - Xs) / tau
+    alpha = 1400 / sqrt(1 + exp((5 [mV] - V) / 6 [mV]))
+    beta = 1 / (1 + exp((V - 35 [mV]) / 15 [mV]))
+    tau = 1 [ms] * alpha * beta + 80 [ms]
+        in [ms]
+    inf = 1 / (1 + exp((-5 [mV] - V) / 14 [mV]))
+g_Ks = 0.2352 [mS/uF]
+    in [mS/uF]
+i_Ks = g_Ks * Xs^2 * (V - nernst.E_Ks)
+    in [A/F]
+
+#
+# (h) Rapid delayed rectifier current, IKr
+# Appendix page 2246
+# Also called the "rapid time-dependent potassium current"
+#
+[ikr]
+use membrane.V
+dot(Xr1) = (xr1_inf - Xr1) / tau_xr1
+    alpha_xr1 = 450 / (1 + exp((-45 [mV] - V) / 10 [mV]))
+    beta_xr1 = 6 / (1 + exp((V + 30 [mV]) / 11.5 [mV]))
+    tau_xr1 = 1 [ms] * alpha_xr1 * beta_xr1
+        in [ms]
+    xr1_inf = 1 / (1 + exp((-26 [mV] - V) / 7 [mV]))
+dot(Xr2) = (xr2_inf - Xr2) / tau_xr2
+    alpha_xr2 = 3 / (1 + exp((-60 [mV] - V) / 20 [mV]))
+    beta_xr2 = 1.12 / (1 + exp((V - 60 [mV]) / 20 [mV]))
+    tau_xr2 = 1 [ms] * alpha_xr2 * beta_xr2
+        in [ms]
+    xr2_inf = 1 / (1 + exp((V + 88 [mV]) / 24 [mV]))
+g_Kr = 0.0918 [mS/uF]
+    in [mS/uF]
+i_Kr = g_Kr * sqrt(extra.K_o / 5.4 [mM]) * Xr1 * Xr2 * (V - nernst.E_K)
+    in [A/F]
+
+#
+# (i) Na/Ca exchange current, INaCa
+# Appendix page 2246
+#
+[inaca]
+use membrane.V, phys.FRT
+use extra.Na_o, extra.Ca_o
+use sodium.Na_i, calcium.Ca_i
+alpha = 2.5
+gamma = 0.35
+K_sat = 0.1
+K_NaCa = 1000 [A/F]
+    in [A/F]
+Km_Nai = 87.5 [mM]
+    in [mM]
+Km_Ca = 1.38 [mM]
+    in [mM]
+i_NaCa = K_NaCa * (exp(gamma * V * FRT) * Na_i^3 * Ca_o - exp((gamma - 1) * V * FRT) * Na_o^3 * Ca_i * alpha) / ((Na_o^3 + Km_Nai^3) * (Km_Ca + Ca_o) * (1 + K_sat * exp((gamma - 1) * V * FRT)))
+    in [A/F]
+
+#
+# (j) Na/K pump current, INaK
+# Appendix page 2247
+#
+[inak]
+use membrane.V, phys.FRT
+use extra.K_o, sodium.Na_i
+K_mNa = 40 [mM]
+    in [mM]
+K_mk = 1 [mM]
+    in [mM]
+P_NaK = 2.724 [A/F]
+    in [A/F]
+i_NaK = P_NaK * K_o / (K_o + K_mk) * Na_i / (Na_i + K_mNa) / (1 + 0.1245 * exp(-0.1 * V * FRT) + 0.0353 * exp(-V * FRT))
+    in [A/F]
+
+#
+# (j) Calcium pump current, IpCa
+# Appendix page 2247
+#
+[ipca]
+K_pCa = 0.0005 [mM]
+    in [mM]
+g_pCa = 0.1238 [A/F]
+    in [A/F]
+i_p_Ca = g_pCa * calcium.Ca_i / (calcium.Ca_i + K_pCa)
+    in [A/F]
+
+#
+# (j) Potassium pump current, IpK
+# Appendix page 2247
+#
+[ipk]
+use membrane.V
+g_pK = 0.0146 [mS/uF]
+    in [mS/uF]
+i_p_K = g_pK * (V - nernst.E_K) / (1 + exp((25 [mV] - V) / 5.98 [mV]))
+    in [A/F]
+
+#
+# (k) Sodium background current, IbNa
+# Appendix page 2247
+#
+[ibna]
+g_bna = 0.00029 [mS/uF]
+    in [mS/uF]
+i_b_Na = g_bna * (membrane.V - nernst.E_Na)
+    in [A/F]
+
+#
+# (k) Calcium background current, IbCa
+# Appendix page 2247
+#
+[ibca]
+g_bca = 0.000592 [mS/uF]
+    in [mS/uF]
+i_b_Ca = g_bca * (membrane.V - nernst.E_Ca)
+    in [A/F]
+
+#
+# (l) Calcium uptake into the SR (SERCA)
+#
+[iup]
+use calcium.Ca_i
+Vmax_up = 0.006375 [mM/ms]
+    in [mM/ms]
+K_up = 0.00025 [mM]
+    in [mM]
+i_up = Vmax_up / (1 + K_up^2 / Ca_i^2)
+    in [mM/ms]
+
+#
+# (l) Calcium release from the SR (RyR)
+#
+[irel]
+use calcium.Ca_ss, calcium.Ca_SR
+k1_prime = 0.15 [1/mM^2/ms]
+    in [1/mM^2/ms]
+k2_prime = 0.045 [1/mM/ms]
+    in [1/mM/ms]
+kcasr = max_sr - (max_sr - min_sr) / (1 + (EC / Ca_SR)^2)
+    max_sr = 2.5
+    min_sr = 1
+    EC = 1.5 [mM]
+        in [mM]
+k1 = k1_prime / kcasr
+    in [1/mM^2/ms]
+k2 = k2_prime * kcasr
+    in [1/mM/ms]
+k3 = 0.06 [1/ms]
+    in [1/ms]
+k4 = 0.005 [1/ms]
+    in [1/ms]
+dot(R_prime) = -k2 * Ca_ss * R_prime + k4 * (1 - R_prime)
+O = k1 * Ca_ss^2 * R_prime / (k3 + k1 * Ca_ss^2)
+V_rel = 0.102 [mS/uF]
+    in [mS/uF]
+i_rel = V_rel * O * (Ca_SR - Ca_ss)
+    in [mM/ms]
+
+#
+# (l) Calcium leak from SR
+#
+[ileak]
+use calcium.Ca_SR, calcium.Ca_i
+V_leak = 0.00036 [mS/uF]
+    in [mS/uF]
+i_leak = V_leak * (Ca_SR - Ca_i)
+    in [mM/ms]
+
+#
+# (l) Diffusion from SS to Bulk
+#
+[ixfer]
+use calcium.Ca_i, calcium.Ca_ss
+V_xfer = 0.0038 [mS/uF]
+    in [mS/uF]
+i_xfer = V_xfer * (Ca_ss - Ca_i)
+    in [mM/ms]
+
+#
+# (l) Calcium dynamics
+# Appendix page 2247
+#
+[calcium]
+use membrane.V
+use cell.V_c, cell.V_sr, cell.V_ss, cell.Cm
+use ileak.i_leak, iup.i_up, irel.i_rel, ixfer.i_xfer
+# Buffering constants
+Buf_c = 0.2 [mM]
+    in [mM]
+Buf_ss = 0.4 [mM]
+    in [mM]
+Buf_sr = 10 [mM]
+    in [mM]
+K_buf_c = 0.001 [mM]
+    in [mM]
+K_buf_ss = 0.00025 [mM]
+    in [mM]
+K_buf_sr = 0.3 [mM]
+    in [mM]
+# Buffered concentrations
+Ca_i_bufc = 1 / (1 + Buf_c * K_buf_c / (Ca_i + K_buf_c)^2)
+Ca_ss_bufss = 1 / (1 + Buf_ss * K_buf_ss / (Ca_ss + K_buf_ss)^2)
+Ca_sr_bufsr = 1 / (1 + Buf_sr * K_buf_sr / (Ca_SR + K_buf_sr)^2)
+# Free calcium in cytosol, SS, and SR
+dot(Ca_i) = Ca_i_bufc * ((i_leak - i_up) * V_sr / V_c + i_xfer - ICa_cyt * Cm / (2 * V_c * phys.F))
+    in [mM]
+dot(Ca_ss) = Ca_ss_bufss * (-ical.i_CaL * Cm / (2 * V_ss * phys.F) + i_rel * V_sr / V_ss - i_xfer * V_c / V_ss)
+    in [mM]
+dot(Ca_SR) = Ca_sr_bufsr * (i_up - (i_rel + i_leak))
+    in [mM]
+# Total calcium currents
+ICa_cyt = ibca.i_b_Ca + ipca.i_p_Ca - 2 * inaca.i_NaCa
+    in [A/F]
+ICa_tot = ICa_cyt + ical.i_CaL
+    in [A/F]
+
+#
+# (m) Sodium dynamics
+# Appendix page 2248
+#
+[sodium]
+INa_tot = ina.i_Na + ibna.i_b_Na + if.i_f_Na + 3 * inak.i_NaK + 3 * inaca.i_NaCa
+    in [A/F]
+dot(Na_i) = -INa_tot / (cell.V_c * phys.F) * cell.Cm
+    in [mM]
+
+#
+# (m) Potassium dynamics
+# Appendix page 2248
+#
+[potassium]
+IK_tot = ik1.i_K1 + ito.i_to + if.i_f_K + isus.i_sus + ikr.i_Kr + iks.i_Ks + ipk.i_p_K - 2 * inak.i_NaK
+    in [A/F]
+dot(K_i) = -IK_tot / (cell.V_c * phys.F) * cell.Cm  # cell.Cm is not in the appendix!
+    in [mM]
+
+#
+# Physical constants
+#
+[phys]
+R = 8.314472 [J/mol/K]
+    in [J/mol/K]
+T = 310 [K]
+    in [K]
+F = 9.64853414999999950e+01 [C/mmol]
+    in [C/mmol]
+RTF = R * T / F
+    in [mV]
+FRT = F / R / T
+    in [1/mV]
+FFRT = F * FRT
+    in [C/mmol/mV]
+
+#
+# Cell geometry
+#
+[cell]
+Cm = 185 [pF]
+    in [pF]
+    desc: Cell capacitance
+V_c = 16404 [um^3]
+    in [um^3]
+    desc: Bulk cytoplasm volume
+V_ss = 54.68 [um^3]
+    in [um^3]
+    desc: Dyadic (junctional) subspace volume
+V_sr = 1094 [um^3]
+    in [um^3]
+    desc: Sarcoplasmic reticulum volume
+
+#
+# Extracellular concentrations
+#
+[extra]
+K_o = 5.4 [mM]
+    in [mM]
+Na_o = 140 [mM]
+    in [mM]
+Ca_o = 2 [mM]
+    in [mM]
+
+#
+# Reversal potentials
+#
+[nernst]
+use extra.K_o, extra.Na_o, extra.Ca_o
+use potassium.K_i, sodium.Na_i, calcium.Ca_i
+E_Ca = 0.5 * phys.RTF * log(Ca_o / Ca_i)
+    in [mV]
+E_K = phys.RTF * log(K_o / K_i)
+    in [mV]
+E_Ks = phys.RTF * log((K_o + P_kna * Na_o) / (K_i + P_kna * Na_i))
+    in [mV]
+E_Na = phys.RTF * log(Na_o / Na_i)
+    in [mV]
+P_kna = 0.03
+
+[[protocol]]
+# Level  Start    Length   Period   Multiplier
+1        100      0.5      1000     0
+
+[[script]]
+import matplotlib.pyplot as plt
+import myokit
+
+# Get model and protocol, create simulation
+m = get_model()
+#p = get_protocol()
+s = myokit.Simulation(m)
+
+# Run simulation
+d = s.run(3000)
+
+# Display the result
+plt.figure()
+plt.xlabel('Time (ms)')
+plt.ylabel('Membrane potential (mV)')
+plt.plot(d['engine.time'], d['membrane.V'])
+plt.show()
+
diff --git a/ft-example/tentusscher-2006-lhopital.mmt b/ft-example/tentusscher-2006-lhopital.mmt
new file mode 100644
index 0000000..276c39f
--- /dev/null
+++ b/ft-example/tentusscher-2006-lhopital.mmt
@@ -0,0 +1,618 @@
+[[model]]
+name: tentusscher-2006-lhopital
+version: 20220502
+mmt_authors: Michael Clerx
+desc: """
+    Human ventricular myocyte model by ten Tusscher et al., 2006.
+
+    The model contains a switch for endo, epi and mid-myocardial modes.
+
+    References:
+
+    [1] The CellML versions of the model (http://models.cellml.org)
+    [2] Alternans and spiral breakup in a human ventricular tissue model,
+        ten Tusscher, Panfilov (2006) AJPHeart
+    [3] A model for human ventricular tissue,
+        ten Tusscher, Noble, Noble, Panfilov (2003) AJPHeart
+
+    Adapted from CellML (October 2020 version). The original CellML model meta
+    data follows below.
+
+    ---------------------------------------------------------------------------
+
+    Alternans and spiral breakup in a human ventricular tissue model
+
+    Penny Noble
+    Oxford University Cardiac Electrophysiology Group
+
+    Model Status
+
+    This is the EPICARDIAL CELL VARIANT of the model. This model was created by
+    Penny Noble of Oxford University and is known to read in COR and PCEnv. A
+    stimulus protocol has been added that allows the model to simulate multiple
+    action potentials at 1Hz.
+
+    Model Structure
+
+    ABSTRACT: Ventricular fibrillation (VF) is one of the main causes of death
+    in the Western world. According to one hypothesis, the chaotic excitation
+    dynamics during VF are the result of dynamical instabilities in action
+    potential duration (APD) the occurrence of which requires that the slope of
+    the APD restitution curve exceeds 1. Other factors such as electrotonic
+    coupling and cardiac memory also determine whether these instabilities can
+    develop. In this paper we study the conditions for alternans and spiral
+    breakup in human cardiac tissue. Therefore, we develop a new version of our
+    human ventricular cell model, which is based on recent experimental
+    measurements of human APD restitution and includes a more extensive
+    description of intracellular calcium dynamics. We apply this model to study
+    the conditions for electrical instability in single cells, for reentrant
+    waves in a ring of cells, and for reentry in two-dimensional sheets of
+    ventricular tissue. We show that an important determinant for the onset of
+    instability is the recovery dynamics of the fast sodium current. Slower
+    sodium current recovery leads to longer periods of spiral wave rotation and
+    more gradual conduction velocity restitution, both of which suppress
+    restitution-mediated instability. As a result, maximum restitution slopes
+    considerably exceeding 1 (up to 1.5) may be necessary for electrical
+    instability to occur. Although slopes necessary for the onset of
+    instabilities found in our study exceed 1, they are within the range of
+    experimentally measured slopes. Therefore, we conclude that steep APD
+    restitution-mediated instability is a potential mechanism for VF in the
+    human heart.
+
+    The original paper reference is cited below:
+
+    Alternans and spiral breakup in a human ventricular tissue model, K.H.W.J.
+    ten Tusscher, A.V. Panfilov, Sep 2006,
+    American Journal of Physiology, Heart and Circulatory Physiology , 291 3,
+    H1088-1100. PubMed ID: 16565318
+    """
+# Initial values (ordered as in the original code)
+membrane.V   = -85.23
+calcium.Cai  = 0.000126
+calcium.CaSR = 3.64
+calcium.CaSS = 0.00036
+sodium.Nai   = 8.604
+potassium.Ki = 136.89
+ina.m        = 0.00172
+ina.h        = 0.7444
+ina.j        = 0.7045
+ikr.xr1      = 0.00621
+ikr.xr2      = 0.4712
+iks.xs       = 0.0095
+ito.r        = 2.42e-8
+ito.s        = 0.999998
+ical.d       = 3.373e-5
+ical.f       = 0.7888
+ical.f2      = 0.9755
+ical.fCaSS   = 0.9953
+jrel.R_prime = 0.9073
+
+#
+# Simulator variables
+#
+[engine]
+time = 0 [ms]
+    in [ms]
+    bind time
+pace = 0
+    bind pace
+
+#
+# Membrane potential
+#
+# Appendix to [3]
+#
+[membrane]
+use stimulus.i_stim
+dot(V) = -(i_ion + i_stim + i_diff)
+    in [mV]
+    label membrane_potential
+i_ion = (+ ina.INa
+         + ik1.IK1
+         + ikr.IKr
+         + iks.IKs
+         + ito.Ito
+         + ical.ICaL
+         + inak.INaK
+         + inaca.INaCa
+         + ipca.IpCa
+         + ipk.IpK
+         + icab.ICab
+         + inab.INab
+        )
+    in [A/F]
+i_diff = 0 [A/F]
+    in [A/F]
+    bind diffusion_current
+
+#
+# Stimulus current
+#
+# Appendix to [3]
+#
+[stimulus]
+i_stim = engine.pace * amplitude
+    in [A/F]
+amplitude = -52 [A/F]
+    in [A/F]
+
+#
+# Cell parameters
+#
+[cell]
+type = 1
+    desc: The type of cell: endocardial = 0, epicardial = 1, mid-myocardial = 2
+Vc = 16404 [um^3]
+    in [um^3]
+    desc: Bulk cytoplasm volume
+Vss = 54.68 [um^3]
+    in [um^3]
+    desc: Dyadic (junctional) subspace volume
+Vsr = 1094 [um^3]
+    in [um^3]
+    desc: Sarcoplasmic reticulum volume
+Cm = 185 [pF]
+    in [pF]
+    desc: Cell capacitance
+
+#
+# Physical constants
+#
+# Appendix to [3]
+#
+[phys]
+F = 96.485 [C/mmol]
+    in [C/mmol]
+R = 8.314 [J/mol/K]
+    in [J/mol/K]
+T = 310 [K]
+    in [K]
+RTF = R * T / F
+    in [mV]
+FRT = F / (R * T)
+    in [1/mV]
+FFRT = F * FRT
+    in [C/mmol/mV]
+
+#
+# Fast sodium current
+#
+# Appendix to [3]
+#
+[ina]
+use membrane.V, nernst.ENa
+INa = gNa * m ^ 3 * h * j * (V - ENa)
+    in [A/F]
+gNa = 14.838 [mS/uF]
+    in [mS/uF]
+dot(m) = (inf - m) / tau
+    inf = 1 / (1 + exp((-56.86 [mV] - V) / 9.03 [mV])) ^ 2
+    tau = 1 [ms] * alpha * beta
+        in [ms]
+    alpha = 1 / (1 + exp((-60 [mV] - V) / 5 [mV]))
+    beta = 0.1 / (1 + exp((V + 35 [mV]) / 5 [mV])) + 0.1 / (1 + exp((V - 50 [mV]) / 200 [mV]))
+dot(h) = (inf - h) / tau
+    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV])) ^ 2
+    tau = 1 / (alpha + beta)
+        in [ms]
+    alpha = if(V < -40 [mV],
+                 0.057 [mS/uF] * exp(-(V + 80 [mV]) / 6.8 [mV]),
+                 0 [mS/uF])
+        in [1/ms]
+    beta = if(V < -40 [mV],
+                2.7 [mS/uF] * exp(0.079 [1/mV] * V) + 310000 [mS/uF] * exp(0.3485 [1/mV] * V),
+                0.77 [mS/uF] / (0.13 * (1 + exp((V + 10.66 [mV]) / -11.1 [mV]))))
+        in [1/ms]
+dot(j) = (inf - j) / tau
+    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV])) ^ 2
+    tau = 1 / (alpha + beta)
+        in [ms]
+    alpha = if(V < -40 [mV],
+               (-25428 [mS/uF] * exp(0.2444 [1/mV] * V) - 6.948e-6 [mS/uF] * exp(-0.04391 [1/mV] * V)) * (V + 37.78 [mV]) / 1 [mV] / (1 + exp(0.311 [1/mV] * (V + 79.23 [mV]))),
+               0 [mS/uF])
+        in [1/ms]
+    beta = if(V < -40 [mV],
+              0.02424 [mS/uF] * exp(-0.01052 [1/mV] * V) / (1 + exp(-0.1378 [1/mV] * (V + 40.14 [mV]))),
+              0.6 [mS/uF] * exp(0.057 [1/mV] * V) / (1 + exp(-0.1 [1/mV] * (V + 32 [mV]))))
+        in [1/ms]
+
+#
+# Inward rectifier potassium current
+#
+# Appendix to [3]
+#
+[ik1]
+use membrane.V, nernst.EK
+IK1 = gK1 * inf * (V - EK)
+    in [A/F]
+gK1 = 5.405 [mS/uF] * sqrt(extra.Ko / 5.4 [mM])
+    in [mS/uF]
+inf = alpha / (alpha + beta)
+    alpha = 0.1 [1/ms] / (1 + exp(0.06 [1/mV] * (V - EK - 200 [mV])))
+        in [1/ms]
+    beta = (3 [1/ms] * exp(0.0002 [1/mV] * (V - EK + 100 [mV])) + 1 [1/ms] * exp(0.1 [1/mV] * (V - EK - 10 [mV]))) / (1 + exp(-0.5 [1/mV] * (V - EK)))
+        in [1/ms]
+
+#
+# Rapid time-dependent potassium current
+#
+# Appendix to [3]
+#
+[ikr]
+use membrane.V, nernst.EK
+IKr = gKr * xr1 * xr2 * (V - EK)
+    in [A/F]
+gKr = 0.153 [mS/uF] * sqrt(extra.Ko / 5.4 [mM])
+    in [mS/uF]
+dot(xr1) = (inf - xr1) / tau
+    inf = 1 / (1 + exp((-26 [mV] - V) / 7 [mV]))
+    tau = 1 [ms] * alpha * beta
+        in [ms]
+    alpha = 450 / (1 + exp((-45 [mV] - V) / 10 [mV]))
+    beta = 6 / (1 + exp((V + 30 [mV]) / 11.5 [mV]))
+dot(xr2) = (inf - xr2) / tau
+    inf = 1 / (1 + exp((V + 88 [mV]) / 24 [mV]))
+    tau = 1 [ms] * alpha * beta
+        in [ms]
+    alpha = 3 / (1 + exp((-60 [mV] - V) / 20 [mV]))
+    beta = 1.12 / (1 + exp((V - 60 [mV]) / 20 [mV]))
+
+#
+# Slow time-dependent potassium current
+#
+# Appendix to [2]
+#
+[iks]
+use membrane.V, nernst.EKs
+IKs = gKs * xs ^ 2 * (V - EKs)
+    in [A/F]
+gKs = if(cell.type != 2,
+         0.392 [mS/uF],      # Endo and epicardial
+         0.098 [mS/uF]       # Mid-myocardial
+        )
+    in [mS/uF]
+dot(xs) = (inf - xs) / tau
+    inf = 1 / (1 + exp((-5 [mV] - V) / 14 [mV]))
+    tau = 1 [ms] * alpha * beta + 80 [ms]
+        in [ms]
+    alpha = 1400 / sqrt(1 + exp((5 [mV] - V) / 6 [mV]))
+    beta = 1 / (1 + exp((V - 35 [mV]) / 15 [mV]))
+
+#
+# Transient outward current
+#
+# Appendix to [2]
+#
+[ito]
+use membrane.V, nernst.EK
+Ito = gto * r * s * (V - EK)
+    in [A/F]
+gto = if(cell.type == 0,
+         0.073 [mS/uF],      # Endocardial
+         0.294 [mS/uF]       # Epicardial and mid-myocardial
+        )
+    in [mS/uF]
+dot(r) = (inf - r) / tau
+    inf = 1 / (1 + exp((20 [mV] - V) / 6 [mV]))
+    tau = 9.5 [ms] * exp(-(V + 40 [mV]) ^ 2 / 1800 [mV^2]) + 0.8 [ms]
+        in [ms]
+dot(s) = (inf - s) / tau
+    inf = if(cell.type == 0,
+             1 / (1 + exp((V + 28 [mV]) / 5 [mV])),     # Endocardial
+             1 / (1 + exp((V + 20 [mV]) / 5 [mV]))      # Epicardial and mid-myocardial
+            )
+    tau = if(cell.type == 0,
+             # Endocardial:
+             1000 [ms] * exp(-(V + 67 [mV]) ^ 2 / 1000 [mV^2]) + 8 [ms],
+             # Epicardial and mid-myocardial:
+             85 [ms] * exp(-(V + 45 [mV]) ^ 2 / 320 [mV^2]) + 5 [ms] / (1 + exp((V - 20 [mV]) / 5 [mV])) + 3 [ms]
+            )
+        in [ms]
+
+#
+# L-type calcium current
+#
+# Appendix to [2]
+#
+[ical]
+use membrane.V
+use phys.FRT, phys.FFRT
+use extra.Cao, calcium.Cai, calcium.CaSS
+ICaL = gCaL * d * f * f2 * fCaSS * 4  * FFRT * if(abs(V - 15) < 1e-5,
+           (0.25 * CaSS - Cao) / (2 * FRT),
+           (V - 15 [mV]) * (0.25 * CaSS * exp(2 * (V - 15 [mV]) * FRT) - Cao) / (exp(2 * (V - 15 [mV]) * FRT) - 1))
+    in [A/F]
+gCaL = 0.0398 [L/F/s]
+    in [L/F/s]
+dot(d) = (inf - d) / tau
+    inf = 1 / (1 + exp((-8 [mV] - V) / 7.5 [mV]))
+    tau = 1 [ms] * (alpha * beta + gamma)
+        in [ms]
+    alpha = 1.4 / (1 + exp((-35 [mV] - V) / 13 [mV])) + 0.25
+    beta = 1.4 / (1 + exp((V + 5 [mV]) / 5 [mV]))
+    gamma = 1 / (1 + exp((50 [mV] - V) / 20 [mV]))
+dot(f) = (inf - f) / tau
+    inf = 1 / (1 + exp((V + 20 [mV]) / 7 [mV]))
+    tau = 1102.5 [ms] * exp(-(V + 27 [mV]) ^ 2 / 225 [mV^2]) + 200 [ms] / (1 + exp((13 [mV] - V) / 10 [mV])) + 180 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV])) + 20 [ms]
+        in [ms]
+dot(f2) = (inf - f2) / tau
+    inf = 0.67 / (1 + exp((V + 35 [mV]) / 7 [mV])) + 0.33
+    tau = 562 [ms] * exp(-(V + 27 [mV]) ^ 2 / 240 [mV^2]) + 31 [ms] / (1 + exp((25 [mV] - V) / 10 [mV])) + 80 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV]))
+        in [ms]
+dot(fCaSS) = (inf - fCaSS) / tau
+    inf = 0.6 / (1 + (CaSS / 0.05 [mM]) ^ 2) + 0.4
+    tau = 80 [ms] / (1 + (CaSS / 0.05 [mM]) ^ 2) + 2 [ms]
+        in [ms]
+
+#
+# Sodium-potassium pump
+#
+# Appendix to [3]
+#
+[inak]
+use membrane.V, phys.FRT
+use extra.Ko, sodium.Nai
+INaK = P_NaK * Ko / (Ko + K_mk) * Nai / (Nai + K_mNa) / (1 + 0.1245 * exp(-0.1 * V * FRT) + 0.0353 * exp(-V * FRT))
+    in [A/F]
+P_NaK = 2.724 [A/F]
+    in [A/F]
+K_mNa = 40 [mM]
+    in [mM]
+K_mk = 1 [mM]
+    in [mM]
+
+#
+# Sodium-calcium exchanger
+#
+# Appendix to [3]
+#
+[inaca]
+use membrane.V, phys.FRT
+use extra.Nao, extra.Cao
+use sodium.Nai, calcium.Cai
+INaCa = (K_NaCa
+         * (exp(gamma * V * FRT) * Nai ^ 3 * Cao - exp((gamma - 1) * V * FRT) * Nao ^ 3 * Cai * alpha)
+         / ((Km_Nai ^ 3 + Nao ^ 3) * (Km_Ca + Cao) * (1 + K_sat * exp((gamma - 1) * V * FRT))))
+    in [A/F]
+K_NaCa = 1000 [A/F]
+    in [A/F]
+Km_Ca = 1.38 [mM]
+    in [mM]
+Km_Nai = 87.5 [mM]
+    in [mM]
+K_sat = 0.1
+alpha = 2.5
+gamma = 0.35
+
+#
+# Calcium pump
+#
+# Appendix to [3]
+#
+[ipca]
+use calcium.Cai
+IpCa = gpCa * Cai / (Cai + KpCa)
+    in [A/F]
+gpCa = 0.1238 [A/F]
+    in [A/F]
+KpCa = 0.0005 [mM]
+    in [mM]
+
+#
+# Potassium pump
+#
+# Appendix to [3]
+#
+[ipk]
+use membrane.V, nernst.EK
+IpK = gpK * (V - EK) / (1 + exp((25 [mV] - V) / 5.98 [mV]))
+    in [A/F]
+gpK = 0.0146 [mS/uF]
+    in [mS/uF]
+
+#
+# Background calcium current
+#
+# Appendix to [3]
+#
+[icab]
+use membrane.V, nernst.ECa
+ICab = gCab * (V - ECa)
+    in [A/F]
+gCab = 0.000592 [mS/uF]
+    in [mS/uF]
+
+#
+# Background sodium current
+#
+# Appendix to [3]
+#
+[inab]
+use membrane.V, nernst.ENa
+INab = gNab * (V - ENa)
+    in [A/F]
+gNab = 0.00029 [mS/uF]
+    in [mS/uF]
+
+#
+# External concentrations
+#
+# Appendix to [3]
+#
+[extra]
+Cao = 2 [mM]
+    in [mM]
+Nao = 140 [mM]
+    in [mM]
+Ko = 5.4 [mM]
+    in [mM]
+
+#
+# Nernst/Reversal potentials
+#
+# Appendix to [3]
+#
+[nernst]
+use extra.Cao, calcium.Cai
+use extra.Nao, sodium.Nai
+use extra.Ko, potassium.Ki
+use phys.RTF
+ECa = RTF * log(Cao / Cai) * 0.5
+    in [mV]
+ENa = RTF * log(Nao / Nai)
+    in [mV]
+EK = RTF * log(Ko / Ki)
+    in [mV]
+EKs = RTF * log((Ko + P_kna * Nao) / (Ki + P_kna * Nai))
+    in [mV]
+P_kna = 0.03
+
+#
+# Calcium release from the SR (RyR)
+#
+[jrel]
+use membrane.V
+use calcium.CaSS, calcium.CaSR
+Jrel = Vrel * O * (CaSR - CaSS)
+    in [mM/ms]
+Vrel = 0.102 [mS/uF]
+    in [mS/uF]
+O = k1 * CaSS ^ 2 * R_prime / (k3 + k1 * CaSS ^ 2)
+kcasr = max_sr - (max_sr - min_sr) / (1 + (EC / CaSR) ^ 2)
+    max_sr = 2.5
+    min_sr = 1
+    EC = 1.5 [mM]
+        in [mM]
+k1 = k1_prime / kcasr
+    in [1/mM^2/ms]
+k1_prime = 0.15 [1/mM^2/ms]
+    in [1/mM^2/ms]
+k2 = k2_prime * kcasr
+    in [1/mM/ms]
+k2_prime = 0.045 [1/mM/ms]
+    in [1/mM/ms]
+k3 = 0.06 [mS/uF]
+    in [mS/uF]
+k4 = 0.005 [mS/uF]
+    in [mS/uF]
+dot(R_prime) = -k2 * CaSS * R_prime + k4 * (1 - R_prime)
+
+#
+# Leak from the SR
+#
+[jleak]
+use calcium.Cai, calcium.CaSR
+Jleak = Vleak * (CaSR - Cai)
+    in [mM/ms]
+Vleak = 0.00036 [mS/uF]
+    in [mS/uF]
+
+#
+# Calcium uptake into the SR (SERCA)
+#
+[jup]
+use calcium.Cai
+Jup = Vmax_up / (1 + K_up ^ 2 / Cai ^ 2)
+    in [mM/ms]
+Vmax_up = 0.006375 [mM/ms]
+    in [mM/ms]
+K_up = 0.00025 [mM]
+    in [mM]
+
+#
+# Diffusion from SS to Bulk
+#
+[jxfer]
+use calcium.Cai, calcium.CaSS
+Jxfer = Vxfer * (CaSS - Cai)
+    in [mM/ms]
+Vxfer = 0.0038 [mS/uF]
+    in [mS/uF]
+
+#
+# Calcium dynamics
+#
+# Appendix to [2]
+#
+[calcium]
+use jup.Jup, jrel.Jrel, jleak.Jleak, jxfer.Jxfer
+use phys.F, cell.Cm, cell.Vc, cell.Vss, cell.Vsr
+# Free calcium in cytosol, SS, and SR
+dot(Cai) = ddt_Cai_total * f_JCai_free
+    in [mM]
+dot(CaSS) = ddt_CaSS_total * f_JCaSS_free
+    in [mM]
+dot(CaSR) = ddt_CaSR_total * f_JCaSR_free
+    in [mM]
+# Derivative of total calcium concentrations (free + buffered)
+ddt_Cai_total = -(icab.ICab + ipca.IpCa - 2 * inaca.INaCa) * Cm / (2 * Vc * F) + (Jleak - Jup) * Vsr / Vc + Jxfer
+    in [mM/ms]
+ddt_CaSS_total = -ical.ICaL * Cm / (2 * Vss * F) + Jrel * Vsr / Vss - Jxfer * Vc / Vss
+    in [mM/ms]
+ddt_CaSR_total = Jup - (Jrel + Jleak)
+    in [mM/ms]
+# Conversion factors from d/dt total to d/dt free
+f_JCai_free = 1 / (1 + Buf_c * K_buf_c / (Cai + K_buf_c) ^ 2)
+f_JCaSS_free = 1 / (1 + Buf_SS * K_buf_SS / (CaSS + K_buf_SS) ^ 2)
+f_JCaSR_free = 1 / (1 + Buf_SR * K_buf_SR / (CaSR + K_buf_SR) ^ 2)
+# Buffering: Ca_buffered = Ca_i * Buf_c / (Ca_i + K_bufc)
+Buf_c = 0.2 [mM]
+    in [mM]
+Buf_SS = 0.4 [mM]
+    in [mM]
+Buf_SR = 10 [mM]
+    in [mM]
+K_buf_c = 0.001 [mM]
+    in [mM]
+K_buf_SS = 0.00025 [mM]
+    in [mM]
+K_buf_SR = 0.3 [mM]
+    in [mM]
+
+#
+# Sodium dynamics
+#
+# Appendix to [3]
+#
+[sodium]
+use phys.F, cell.Cm, cell.Vc
+dot(Nai) = -INa_total * Cm / (Vc * F)
+    in [mM]
+INa_total = ina.INa + inab.INab + 3 * inak.INaK + 3 * inaca.INaCa
+    in [A/F]
+
+#
+# Potassium dynamics
+#
+# Appendix to [3]
+#
+[potassium]
+use phys.F, cell.Cm, cell.Vc
+dot(Ki) = -IK_total * Cm / (Vc * F)
+    in [mM]
+IK_total = ik1.IK1 + ito.Ito + ikr.IKr + iks.IKs + ipk.IpK + stimulus.i_stim - 2 * inak.INaK
+    in [A/F]
+
+[[protocol]]
+# Level  Start    Length   Period   Multiplier
+1.0      10.0     1.0      1000.0   0
+
+[[script]]
+import matplotlib.pyplot as plt
+import myokit
+
+# Get model and protocol, create simulation
+m = get_model()
+p = get_protocol()
+s = myokit.Simulation(m, p)
+
+# Run simulation
+d = s.run(1000)
+
+# Display the results
+plt.figure()
+plt.xlabel('Time (ms)')
+plt.ylabel('Membrane potential (mV)')
+plt.plot(d.time(), d['membrane.V'])
+plt.show()
+