Numerical schemes#
gotranx can generate more than the right-hand side \(f(t, y)\) of an ODE. It can also generate a function that takes one time step with a given numerical scheme. Such a function has the signature
def scheme_name(states, t, dt, parameters): ...
and returns the state values at time \(t + \Delta t\). This is useful when the ODE solver must be very cheap per step, for example when the ODE is solved at every node of a mesh in a coupled PDE–ODE simulation (see Vectorized computations), or when the generated code is used from a language without a good ODE solver library.
Throughout this section we consider a system of \(N\) ordinary differential equations
discretized with a constant time step \(\Delta t\) such that \(t_n = t_0 + n \Delta t\) and \(y^n \approx y(t_n)\).
Available schemes#
Scheme |
|
Order |
Treats exponentially |
Page |
|---|---|---|---|---|
Forward (explicit) Euler |
|
1 |
no state |
|
Rush–Larsen |
|
1 |
the gating variables |
|
Generalized Rush–Larsen |
|
1 |
every state |
|
Hybrid Rush–Larsen |
|
1 |
the states you choose |
You can list the schemes from the command line with
gotranx list-schemes
and generate code for one or more of them with the --scheme flag, e.g.
gotranx ode2py model.ode --scheme explicit_euler --scheme generalized_rush_larsen -o model.py
The same options can be set in a configuration file (see Configuration), and from Python with the scheme argument of the code generators, e.g.
gotranx.cli.gotran2py.get_code().
Which scheme should I use?#
All four schemes are explicit, single-step and first-order accurate, so they cost roughly one evaluation of the right-hand side per step. They differ in stability:
Forward Euler is the simplest and cheapest scheme, but it is only stable when the time step is small compared to the fastest time scale of the model.
The Rush–Larsen family integrates some or all states with an exponential update that is exact for linear equations. This removes the time step restriction coming from those states, at the cost of one exponential per state per step.
Generalized Rush–Larsen applies the exponential update to every state and is the most robust choice. It is a good default for cardiac cell models.
Hybrid Rush–Larsen applies it only to the states you list, so you only pay for the exponentials you need.
A larger stable time step does not mean a more accurate solution. All schemes are first order, and for a given \(\Delta t\) forward Euler is sometimes more accurate than the exponential schemes when it is stable. Check the accuracy of your simulation by refining the time step.
The running example#
To compare the schemes we use the classical model of the squid giant axon by Hodgkin and Huxley (1952). It has a membrane potential \(V\) and three gating variables \(m\), \(h\) and \(n\). You can download the model in .ode format here.
# Hodgkin-Huxley model of the squid giant axon
#
# Hodgkin AL, Huxley AF. A quantitative description of membrane current and its
# application to conduction and excitation in nerve. J Physiol. 1952;117(4):500-544.
# https://doi.org/10.1113/jphysiol.1952.sp004764
#
# Written with the modern sign convention (resting potential at -65 mV).
parameters("membrane",
Cm = ScalarParam(1.0, unit="uF/cm**2", description="Membrane capacitance"),
stim_amp = ScalarParam(10.0, unit="uA/cm**2", description="Stimulus amplitude"),
stim_start = ScalarParam(1.0, unit="ms", description="Stimulus start"),
stim_duration = ScalarParam(1.0, unit="ms", description="Stimulus duration")
)
parameters("sodium_channel",
g_Na = ScalarParam(120.0, unit="mS/cm**2", description="Maximal sodium conductance"),
E_Na = ScalarParam(50.0, unit="mV", description="Sodium reversal potential")
)
parameters("potassium_channel",
g_K = ScalarParam(36.0, unit="mS/cm**2", description="Maximal potassium conductance"),
E_K = ScalarParam(-77.0, unit="mV", description="Potassium reversal potential")
)
parameters("leakage_current",
g_L = ScalarParam(0.3, unit="mS/cm**2", description="Leak conductance"),
E_L = ScalarParam(-54.387, unit="mV", description="Leak reversal potential")
)
states("membrane",
V = ScalarParam(-65.0, unit="mV", description="Membrane potential")
)
states("sodium_channel_m_gate",
m = ScalarParam(0.05, description="Sodium activation gate")
)
states("sodium_channel_h_gate",
h = ScalarParam(0.6, description="Sodium inactivation gate")
)
states("potassium_channel_n_gate",
n = ScalarParam(0.32, description="Potassium activation gate")
)
expressions("sodium_channel_m_gate")
alpha_m = 0.1 * (V + 40.0) / (1.0 - exp(-(V + 40.0) / 10.0)) # 1/ms
beta_m = 4.0 * exp(-(V + 65.0) / 18.0) # 1/ms
dm_dt = alpha_m * (1 - m) - beta_m * m
expressions("sodium_channel_h_gate")
alpha_h = 0.07 * exp(-(V + 65.0) / 20.0) # 1/ms
beta_h = 1.0 / (1.0 + exp(-(V + 35.0) / 10.0)) # 1/ms
dh_dt = alpha_h * (1 - h) - beta_h * h
expressions("potassium_channel_n_gate")
alpha_n = 0.01 * (V + 55.0) / (1.0 - exp(-(V + 55.0) / 10.0)) # 1/ms
beta_n = 0.125 * exp(-(V + 65.0) / 80.0) # 1/ms
dn_dt = alpha_n * (1 - n) - beta_n * n
expressions("membrane")
i_Na = g_Na * m**3 * h * (V - E_Na) # uA/cm**2
i_K = g_K * n**4 * (V - E_K) # uA/cm**2
i_L = g_L * (V - E_L) # uA/cm**2
i_stim = Conditional(And(Ge(t, stim_start), Le(t, stim_start + stim_duration)), stim_amp, 0.0) # uA/cm**2
dV_dt = (i_stim - i_Na - i_K - i_L) / Cm # mV/ms