Order of accuracy for an active cylinder solved with CN

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

Order of accuracy for an active cylinder solved with CN

Post by danarwed »

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)) {
        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
}
ted
Site Admin
Posts: 6394
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

Post by ted »

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

Post by danarwed »

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
        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
}
 
? 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
Site Admin
Posts: 6394
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

Post by ted »

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.
Post Reply