From 69c1561f8f05f3e78299484e5309c4ed0a04cbda Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Thu, 13 Oct 2022 14:37:13 +0200
Subject: [PATCH] Discuss q10 pattern in NMODL docs. (#1995)
* Discuss q10 pattern in NMODL docs.
---
doc/fileformat/nmodl.rst | 45 +++++++++++++++++++++++++++++++++++++++
mechanisms/default/hh.mod | 26 +++++++++-------------
2 files changed, 55 insertions(+), 16 deletions(-)
diff --git a/doc/fileformat/nmodl.rst b/doc/fileformat/nmodl.rst
index 3ebf3d3d..b975fc59 100644
--- a/doc/fileformat/nmodl.rst
+++ b/doc/fileformat/nmodl.rst
@@ -374,6 +374,51 @@ Note that we do not lose accuracy here, since Arbor does not support
higher-order ODEs and thus will treat ``g`` as a constant across
a single timestep even if ``g`` actually depends on ``v``.
+Using Memory versus Computation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Commonly ion channels need to correct for temperature differences, which yields
+a term similar to
+
+.. code::
+
+ q = 3^(0.1*celsius - 0.63)
+
+Here, we find that the cost of the exponential when computing ``q`` in the
+``DERIVATIVE`` block is high enough to make pre-computing ``q`` in ``INITIAL``
+and loading the value later an optimisation. Shown below is a simplified version
+of this pattern from ``hh.mod`` in the Arbor sources
+
+.. code::
+
+ NEURON {
+ ...
+ RANGE ..., q
+ }
+
+ ASSIGNED { q }
+
+ PARAMETER {
+ ...
+ celsius (degC)
+ }
+
+ STATE { ... }
+
+ BREAKPOINT {
+ SOLVE dS METHOD cnexp
+ ...
+ }
+
+ INITIAL {
+ q = 3^(0.1*celsius - 0.63)
+ ...
+ }
+
+ DERIVATIVE states {
+ ... : uses q
+ }
+
Specialised Functions
~~~~~~~~~~~~~~~~~~~~~
diff --git a/mechanisms/default/hh.mod b/mechanisms/default/hh.mod
index 9c64142f..88c28346 100644
--- a/mechanisms/default/hh.mod
+++ b/mechanisms/default/hh.mod
@@ -38,7 +38,7 @@ BREAKPOINT {
INITIAL {
LOCAL alpha, beta
- q10 = 3^((celsius - 6.3)/10.0)
+ q10 = 3^(0.1*celsius - 0.63)
: sodium activation system
alpha = m_alpha(v)
@@ -53,38 +53,32 @@ INITIAL {
: potassium activation system
alpha = n_alpha(v)
beta = n_beta(v)
-
n = alpha/(alpha + beta)
}
DERIVATIVE states {
- LOCAL alpha, beta, sum
+ LOCAL alpha, beta
: sodium activation system
alpha = m_alpha(v)
beta = m_beta(v)
- sum = alpha + beta
- m' = (alpha - m*sum)*q10
+ m' = (alpha - m*(alpha + beta))*q10
: sodium inactivation system
alpha = h_alpha(v)
beta = h_beta(v)
- sum = alpha + beta
- h' = (alpha - h*sum)*q10
+ h' = (alpha - h*(alpha + beta))*q10
: potassium activation system
alpha = n_alpha(v)
beta = n_beta(v)
- sum = alpha + beta
- n' = (alpha - n*sum)*q10
+ n' = (alpha - n*(alpha + beta))*q10
}
-FUNCTION vtrap(x,y) { vtrap = y*exprelr(x/y) }
-
-FUNCTION m_alpha(v) { m_alpha = 0.1*vtrap(-(v + 40.0), 10.0) }
-FUNCTION h_alpha(v) { h_alpha = 0.07*exp(-(v + 65.0)/20.0) }
-FUNCTION n_alpha(v) { n_alpha = 0.01*vtrap(-(v + 55.0), 10.0) }
+FUNCTION m_alpha(v) { m_alpha = exprelr(-0.1*v - 4.0) }
+FUNCTION h_alpha(v) { h_alpha = 0.07*exp(-0.05*v - 3.25) }
+FUNCTION n_alpha(v) { n_alpha = 0.1*exprelr(-0.1*v - 5.5) }
FUNCTION m_beta(v) { m_beta = 4.0*exp(-(v + 65.0)/18.0) }
-FUNCTION h_beta(v) { h_beta = 1.0/(exp(-(v + 35.0)/10.0) + 1.0) }
-FUNCTION n_beta(v) { n_beta = 0.125*exp(-(v + 65.0)/80.0) }
+FUNCTION h_beta(v) { h_beta = 1.0/(exp(-0.1*v - 3.5) + 1.0) }
+FUNCTION n_beta(v) { n_beta = 0.125*exp(-0.0125*v - 0.8125) }
--
GitLab