## Order of accuracy for an active cylinder solved with CN

Anything that doesn't fit elsewhere.
danarwed
Posts: 6
Joined: Fri Jan 08, 2010 8:05 am

### Order of accuracy for an active cylinder solved with CN

Dear Neuroners,

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)) {
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
}

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
}
``````
ted
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Order of accuracy for an active cylinder solved with CN

Time is discretized and space is discretized. secondorder specifies the order of accuracy in time. What happens to error if you change dt?
danarwed
Posts: 6
Joined: Fri Jan 08, 2010 8:05 am

### Order of accuracy for an active cylinder solved with CN

Hi Ted,

thanks for your fast reply (as always!). Let me figure out that I started with nseg=3 and dt=0.01 ms and then successively refined the time steps size, dt -> dt/3 while increasing nseg by 3, nseg -> nseg*3, and these are the results I got:

Code: Select all

``````//tend   nseg       dt                    soma.v(0.5) at t=10ms
a) 10        1      0.01                  -9.15806317753224
b) 10        3      0.0033333333          -9.148834936342411
c) 10        9      0.0011111111          -9.145007547480585
d) 10       27      0.00037037            -9.143599849408666
e) 10       81      0.00012345679         -9.143115464857985
f) 10      243      4.11523e-05           -9.14295231432898
``````
Here I think I observe the 3rd order convergence, regarding the difference between the voltage values b) -a), c)-b), d)-c) and so on:

0.00922824, 0.0038273, 0.0014076, 0.00048438, 0.00016315

Let me further note that I get second order convergence in the case of a passive mechanism inserted instead of the active ones, such that I assume that my mistake has to deal with my 'myhh.mod'-file [Ooops, it was somehow truncated in the last post, hopefully it'll be better now]:

Code: Select all

``````TITLE myhh.mod   squid sodium, potassium, and leak channels

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
}
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
}

? currents
BREAKPOINT {
SOLVE states METHOD derivimplicit
: SOLVE states METHOD cnexp
gna = gnabar*m*m*m*h
ina = gna*(v - ena)
gk = gkbar*n*n*n*n
ik = gk*(v - ek)
il = gl*(v - el)
}

INITIAL {
rates(v)
m = minf
h = hinf
n = ninf
alpham = alpha_m
betam = beta_m
alphan = alpha_n
betan = beta_n
alphah = alpha_h
betah = beta_h
}

? states
DERIVATIVE states {
rates(v)
m' =  (minf-m)/mtau
h' = (hinf-h)/htau
n' = (ninf-n)/ntau
}

:LOCAL q10

? rates
PROCEDURE rates(v(mV)) {  :Computes rate and other constants at current v.
:Call once from HOC to initialize inf at resting v.
LOCAL  alpha, beta, sum, q10

UNITSOFF
q10 = 3^((celsius - 6.3)/10)
:"m" sodium activation system
alpha = .1 * vtrap(-(v-25),10)
beta =  4 * exp(-(v)/18)

alpha_m = alpha
beta_m = beta

sum = alpha + beta
mtau = 1/(q10*sum)
minf = alpha/sum
:"h" sodium inactivation system
alpha = .07 * exp(-v/20)
beta = 1 / (exp(-(v-30)/10) + 1)

alpha_h = alpha
beta_h = beta

sum = alpha + beta
htau = 1/(q10*sum)
hinf = alpha/sum
:"n" potassium activation system
alpha = .01*vtrap(-(v-10),10)
beta = .125*exp(-v/80)

alpha_n = alpha
beta_n = beta

sum = alpha + beta
ntau = 1/(q10*sum)
ninf = alpha/sum
}

FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
: Derivation of this approximation:
: see http://www.neuron.yale.edu/phpBB/viewtopic.php?f=15&t=1075
if (fabs(x/y) < 1e-6) {
vtrap = y*(1 - x/y/2)
}else{
vtrap = x/(exp(x/y) - 1)
}
}

UNITSON
``````
Is the integration method

Code: Select all

``````      SOLVE states METHOD derivimplicit
: SOLVE states METHOD cnexp
``````
the proper choice for integrating the channels?

Best,
Dan
ted
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Order of accuracy for an active cylinder solved with CN

Good questions. It's important to think clearly about these issues, and take nothing for granted.
danarwed wrote:I started with nseg=3 and dt=0.01 ms and then successively refined the time steps size, dt -> dt/3 while increasing nseg by 3, nseg -> nseg*3
You're changing two things that affect accuracy: the spatial grid, and the temporal grid. This makes it difficult to determine which change did what to the solution. Keep it simple. To investigate the relationship between dt (time discretization) and error, eliminate space by using a model with just one compartment, and change dt. To study the relationship between dx (space discretization) and error, eliminate time by examining how nseg affects the steady state response to injection of current at one point in a model of a "long" cylinder with uniform diameter and electrical properties.

In either case (time course of v when current is injected into a single compartment model; steady state v along a finite cable elicited by injecting a constant current at one point) the analytical solution is known. The error is the difference between the analytical solution and the solution generated with a particular dt (or, for the cable model, for a particular L/nseg).
I think I observe the 3rd order convergence
Well, no, but more about that later. First a definition of terms: this has nothing to do with convergence. It's all about "order of accuracy," also known as "discretization order" or "truncation order". For a very brief discussion see "Order of Accuracy" http://www.cmth.ph.ic.ac.uk/people/a.ma ... node5.html, and for a very clear discussion see http://en.wikipedia.org/wiki/Finite_dif ... _and_order. I have changed the "Subject" lines for this thread accordingly.

If you want to determine the order of accuracy of a numerical method with regard to some independent variable z, you need to determine how the error of the solutions generated with that method are affected by changes in the discretization of z. This implies that you start with a _standard_ solution, i.e. an _accurate_ solution, against which you compare inaccurate solutions generated with two or more discretizations of z. Call the accurate solution sa. Then
dz1 produces solution s1 = sa + error1
dz2 produces solution s2 = sa + error2
If it turns out that error1/error2 is close to dz1/dz2, then error is proportional to dz and the method has "first order accuracy". If error1/error2 is closer to (dz1/dz2)^2, error is proportional to dz^2 and the method has "second order accuracy".

But what if you try to compare inaccurate solutions with each other by finding
s1 - s2 = error1 - error2
s2 - s3 = error2 - error3
?
You can't just divide s1-s2 by s2-s3 and call that the "order of accuracy."

For practical examples using NEURON, see chapter 4 of The NEURON Book, especially
4.2.3.1 Effcient handling of nonlinearity and figures 4.9 and 4.10 on pages 70-72.
Is the integration method derivimplicit the proper choice for integrating the channels?
It is no more accurate than cnexp, but it chews up more cpu time. See Integration methods for SOLVE statements in the Hot tips area of the NEURON Forum.