Running HH kernel via LLVM backend

[georgemitenkov]

Running hh.mod:

TITLE hh.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 and shifted to reflect a resting potential
    of -65 mV.
    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 hh
    USEION na READ ena WRITE ina
    USEION k READ ek WRITE ik
    NONSPECIFIC_CURRENT il
    RANGE gnabar, gkbar, gl, el, gna, gk
    RANGE minf, hinf, ninf, mtau, htau, ntau
    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
}

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

? currents
BREAKPOINT {
    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
}

? 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+40),10)
    beta =  4 * exp(-(v+65)/18)
    sum = alpha + beta
    mtau = 1/(q10*sum)
    minf = alpha/sum
    :"h" sodium inactivation system
    alpha = .07 * exp(-(v+65)/20)
    beta = 1 / (exp(-(v+35)/10) + 1)
    sum = alpha + beta
    htau = 1/(q10*sum)
    hinf = alpha/sum
    :"n" potassium activation system
    alpha = .01*vtrap(-(v+55),10)
    beta = .125*exp(-(v+65)/80)
    sum = alpha + beta
    ntau = 1/(q10*sum)
    ninf = alpha/sum
}

FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
    if (fabs(x/y) < 1e-6) {
        vtrap = y*(1 - x/y/2)
    } else {
        vtrap = x/(exp(x/y) - 1)
    }
}

UNITSON

currently produces the following NMODL kernel (used passes --inline --const-folding):

VOID nrn_state_hh(INSTANCE_STRUCT *mech){
    INTEGER id
    INTEGER node_id, ena_id, ek_id
    DOUBLE v
    for(id = 0; id<mech->node_count-7; id = id+8) {
        node_id = mech->node_index[id]
        ena_id = mech->ion_ena_index[id]
        ek_id = mech->ion_ek_index[id]
        v = mech->voltage[node_id]
        mech->ena[id] = mech->ion_ena[ena_id]
        mech->ek[id] = mech->ion_ek[ek_id]
        {
            DOUBLE alpha, beta, sum, q10, vtrap_in_0, vtrap_in_1, v_in_1
            v_in_1 = v
            UNITSOFF
            q10 = 3^((mech->celsius-6.3)/10)
            {
                DOUBLE x_in_0, y_in_0
                x_in_0 = -(v_in_1+40)
                y_in_0 = 10
                IF (fabs(x_in_0/y_in_0)<1e-6) {
                    vtrap_in_0 = y_in_0*(1-x_in_0/y_in_0/2)
                } ELSE {
                    vtrap_in_0 = x_in_0/(exp(x_in_0/y_in_0)-1)
                }
            }
            alpha = .1*vtrap_in_0
            beta = 4*exp(-(v_in_1+65)/18)
            sum = alpha+beta
            mech->mtau[id] = 1/(q10*sum)
            mech->minf[id] = alpha/sum
            alpha = .07*exp(-(v_in_1+65)/20)
            beta = 1/(exp(-(v_in_1+35)/10)+1)
            sum = alpha+beta
            mech->htau[id] = 1/(q10*sum)
            mech->hinf[id] = alpha/sum
            {
                DOUBLE x_in_1, y_in_1
                x_in_1 = -(v_in_1+55)
                y_in_1 = 10
                IF (fabs(x_in_1/y_in_1)<1e-6) {
                    vtrap_in_1 = y_in_1*(1-x_in_1/y_in_1/2)
                } ELSE {
                    vtrap_in_1 = x_in_1/(exp(x_in_1/y_in_1)-1)
                }
            }
            alpha = .01*vtrap_in_1
            beta = .125*exp(-(v_in_1+65)/80)
            sum = alpha+beta
            mech->ntau[id] = 1/(q10*sum)
            mech->ninf[id] = alpha/sum
        }
        mech->m[id] = mech->m[id]+(1.0-exp(mech->dt*((((-1.0)))/mech->mtau[id])))*(-(((mech->minf[id]))/mech->mtau[id])/((((-1.0)))/mech->mtau[id])-mech->m[id])
        mech->h[id] = mech->h[id]+(1.0-exp(mech->dt*((((-1.0)))/mech->htau[id])))*(-(((mech->hinf[id]))/mech->htau[id])/((((-1.0)))/mech->htau[id])-mech->h[id])
        mech->n[id] = mech->n[id]+(1.0-exp(mech->dt*((((-1.0)))/mech->ntau[id])))*(-(((mech->ninf[id]))/mech->ntau[id])/((((-1.0)))/mech->ntau[id])-mech->n[id])
    }

  // remainder of the loop
}

[georgemitenkov]

Problems and observations

1. Inlined function still remain in AST, and the visitor is dispatched for them. Hence, if the function uses instance variable field (e.g. rates uses celsius) LLVM backend throws an error that the variable is not defined.

2. LOCAL allocations are not pulled out (fine for now)

3. There are unsupported operators (e.g. ^)

[georgemitenkov]

Working hh.mod

TITLE hh.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 and shifted to reflect a resting potential
    of -65 mV.
    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)
}

NEURON {
    SUFFIX hh
    USEION na READ ena WRITE ina
    USEION k READ ek WRITE ik
    NONSPECIFIC_CURRENT il
    RANGE gnabar, gkbar, gl, el, gna, gk
    RANGE minf, hinf, ninf, mtau, htau, ntau
    THREADSAFE
}

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
}

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

BREAKPOINT {
    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 {
    {
        : inlined rates
        LOCAL alpha, beta, sum, q10, vtrap_in_0, v_in_0
        v_in_0 = v
        q10 = 3*((celsius-6.3)/10)
        alpha = .07*exp(-(v_in_0+65)/20)
        beta = 1/(exp(-(v_in_0+35)/10)+1)
        sum = alpha+beta
        htau = 1/(q10*sum)
        hinf = alpha/sum
        {
            : inlined vtrap
            LOCAL x_in_0, y_in_0
            x_in_0 = alpha
            y_in_0 = alpha

            : no control flow
            vtrap_in_0 = y_in_0*(1-x_in_0/y_in_0/2)
        }
        hinf = vtrap_in_0
    }
    m = minf
    h = hinf
    n = ninf
}

DERIVATIVE states {
    {
        : inlined rates
        LOCAL alpha, beta, sum, q10, vtrap_in_0, v_in_1
        v_in_1 = v
        q10 = 3*((celsius-6.3)/10)
        alpha = .07*exp(-(v_in_1+65)/20)
        beta = 1/(exp(-(v_in_1+35)/10)+1)
        sum = alpha+beta
        htau = 1/(q10*sum)
        hinf = alpha/sum
        {
           : inlined vtrap
            LOCAL x_in_0, y_in_0
            x_in_0 = alpha
            y_in_0 = alpha

            : no control flow
            vtrap_in_0 = y_in_0*(1-x_in_0/y_in_0/2)  
        }
        hinf = vtrap_in_0
    }
    m = m+(1.0-exp(dt*((((-1.0)))/mtau)))*(-(((minf))/mtau)/((((-1.0)))/mtau)-m)
    h = h+(1.0-exp(dt*((((-1.0)))/htau)))*(-(((hinf))/htau)/((((-1.0)))/htau)-h)
    n = n+(1.0-exp(dt*((((-1.0)))/ntau)))*(-(((ninf))/ntau)/((((-1.0)))/ntau)-n)
}

UNITSON

[pramodk]

@georgemitenkov :

Inlined function still remain in AST,

this we can remove easily. Will get this done tomorrow.

There are unsupported operators (e.g. ^)

this would be straightforward / minor I assume (pow function).

LOCAL allocations are not pulled out (fine for now)

we can try to get this done, shouldn't be much efforts.

[pramodk]

@iomaganaris : we can discuss this tomorrow morning and you can help to get 1st and 3rd done.

[georgemitenkov]

2nd is addressed here #661

Regarding the 3rd one: this is not urgent, as we will optimize IR for benchmarking anyway (so all allocas will go away). But nice to have that!

[georgemitenkov]

@iomaganaris @pramodk LOCAL allocations are now pulled out (see #663), so its only the duplicate AST nodes that stop us from running hh.mod 😊

[pramodk]

Perfect! Ioannis has WIP changes already in pipeline to handle duplicate AST nodes. So before tomorrow's sync, we should have everything ready!