I try to solve an active model on a uniform cylinder and investigate the order of convergence one can achieve with the Crank-Nicolson method (CN). The channels in the active model are the classical Hodgkin-Huxley channels from the giant quid axon, but rescaled so as to have the resting potential v_r = 0. In the first listing below you can see the hoc-file, in the second my local myhh.mod (a slightly changed hh.mod) file containing the definition of the channels.
To check convergence, I start with a given spatial resolution (nseg=) and time step size (dt=) and refine them uniformly obeying Neuron's odd segment and multiply by 3 rule for nseg and, as I would expect second order for CN, I divide the time step size by 3 eveytime I multiply nseg by 3. Here is an example of my results:
// tend nseg dt soma.v(0.5) at t=10ms
10 1 0.01 -9.15806317753224
10 3 0.0033333333 -9.148834936342411
10 9 0.0011111111 -9.145007547480585
10 27 0.00037037 -9.143599849408666
10 81 0.00012345679 -9.143115464857985
10 243 4.11523e-05 -9.14295231432898
This gives a reduction of the error (from nseg=3 on, I compute the difference between the solution value with the current dt and nseg and the solution between the solution on one level coarser):
0.00922824, 0.0038273, 0.0014076, 0.00048438, 0.00016315
From the documentation I understand that even in the active case one would expect second order convergence/reduction in the potential, but, as one can see here, the error reduction is of third order only.
My question was if my setupt is wrong (I get third order convergence for all combinations of nseg/dt I tried) or do I misunderstood something or miss a thing? Any help will be greatly aprreciated!
Thanks & best regards
Dan
Code: Select all
create soma
secondorder = 2
/* specify anatomical and biophysical properties */
soma {
nseg=9
diam = 4.0 // [um], soma diameter tapers along its length
L = 100. 0 // [um], length
}
objectvar syn1
soma syn1 = new AlphaSynapse(0.5)
syn1.tau = 1.0 // [ms] t_peak
syn1.onset = 1.0 // [ms] t_start
syn1.e = 50 // [mV] synaptical reversal potential
syn1.gmax = 0.005 // [uS] maximal synaptical conductivity
proc init() {
celsius = 6.3
dt=0.01
tstop = 10. // [ms], end time
finitialize(0.0) // [mv], init membrane potential
forall {
nseg=9
// channels
insert myhh
Ra = 200
el_myhh = 0
ek = -12
ena = 115
gl_myhh = 0.0005
cm = 1. // [uf/cm^2]
// initialize gating particles
m_myhh = 0.0529551
n_myhh = 0.317732
h_myhh = 0.595994
}
if (cvode.active()) {
cvode.re_init()
} else {
fcurrent()
}
frecord_init()
}
proc integrate() {
printf("%.8g %8d %12.8g %20.16g %20.16g %20.16g %20.16g\n", t, nseg, dt, soma.v(0.5), soma.m_myhh, soma.n_myhh, soma.h_myhh)
while(t < tstop-(dt*1e-5)) {
fadvance()
printf("%.8g %8d %12.8g %20.16g %20.16g %20.16g %20.16g\n", t, nseg, dt, soma.v(0.5), soma.m_myhh, soma.n_myhh, soma.h_myhh)
}
}
/* run */
init()
integrate()
quit()
Code: Select all
TITLE myhh.mod squid sodium, potassium, and leak channels
COMMENT
This is the original Hodgkin-Huxley treatment for the set of sodium,
potassium, and leakage channels found in the squid giant axon membrane.
("A quantitative description of membrane current and its application
conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).)
Membrane voltage is in absolute mV and has been reversed in polarity
from the original HH convention.
Remember to set celsius=6.3 (or whatever) in your HOC file.
See squid.hoc for an example of a simulation using this model.
SW Jaslove 6 March, 1992
ENDCOMMENT
UNITS {
(mA) = (milliamp)
(mV) = (millivolt)
(S) = (siemens)
}
? interface
NEURON {
SUFFIX myhh
USEION na READ ena WRITE ina
USEION k READ ek WRITE ik
NONSPECIFIC_CURRENT il
RANGE gnabar, gkbar, gl, el, gna, gk
GLOBAL minf, hinf, ninf, mtau, htau, ntau, alpha_m, beta_m, alpha_n, beta_n, alpha_h, beta_h
THREADSAFE : assigned GLOBALs will be per thread
}
PARAMETER {
gnabar = .12 (S/cm2) <0,1e9>
gkbar = .036 (S/cm2) <0,1e9>
gl = .0003 (S/cm2) <0,1e9>
el = -54.3 (mV)
}
STATE {
m h n alpham betam alphan betan alphah betah
}
ASSIGNED {
v (mV)
celsius (degC)
ena (mV)
ek (mV)
gna (S/cm2)
gk (S/cm2)
ina (mA/cm2)
ik (mA/cm2)
il (mA/cm2)
minf hinf ninf
mtau (ms) htau (ms) ntau (ms)
alpha_m beta_m
alpha_n beta_n
alpha_h beta_h
}