Further analysis of Arbor vs. Neuron discrepancy for BBP mechanisms (K_Tst and K_Pst)

[lukasgd]

In the BluePyOpt integration, a validation test between Arbor and Neuron has been developed and evaluated for the BBP catalogue that is used in the layer-5 pyramidal cell model. While showing good agreement for most mechanisms, there is a discrepancy for K_Tst and K_Pst that was discussed here (cf. latest results) and cannot be explained by differences in the NMODL definition (driving Neuron with Arbor's NMODL gives the same discrepancy).

The purpose of this issue is to identify the source of this error, possibly by first time-integrating the ODEs for a single-compartimental setup with an individual mechanism using scipy to identify the correct solution and then narrowing down the error in the affected simulator.

[thorstenhater]

Interestingly, just running a Ball-and-Stick cell with any of the BBP channels --- incl K_P and K_T, but excl Nap --- doesn't show measurable differences. Nap is different, there NRN 8 diverges with t, while ARB just show a beginning of unbounded growth.

I'd like a more condensed version of the test case that demonstrates the issue, excl all the surrounding noise and with a significantly simpler morphology.

[thorstenhater]

I found the issue. As is New Year's tradition I am trawling the backlog for issues to close and prioritize. Since I had some problems with the Allen mechanisms I decided to check for the same problem here, but now luck. Then I decided to take one look at the generated C++. With SIMD it's ironically easier to see:

DERIVATIVE states {
    LOCAL qt, mInf, hInf, mTau, hTau

    qt   = 2.3^((34 - 21)/10)
    mInf = m_inf(v)
    # vvvvv
    mTau = 0.34 + 0.92*exp(-((v + 81)/59)^2)
    # ^^^^^
    hInf = h_inf(v)
    hTau = 8 + 49*exp(-((v + 83)/23)^2)

    m' = qt*(mInf - m)/mTau
    h' = qt*(hInf - h)/hTau
}

The marked line becomes in the the C++ code:

assign(mTau, S::add(0.34,
                     S::mul(0.92,
                             S::exp(S::mul(S::neg(S::mul(S::add(v, 81.0), 0.016949152542372881)),
                                           S::neg(S::mul(S::add(v, 81.0), 0.016949152542372881)))))));

(modulo clean-up) and one can basically see how precedence is mishandled. For reference, -(...)^2 is misinterpreted as (-(...))^2.

[thorstenhater]

after adding some special handling for power when parsing unary (negation)

        assign(mTau,
               S::add(0.34,
                      S::mul(0.92,
                             S::exp(S::neg(S::mul(S::mul(S::add(v, 81.0), 0.016949152542372881),
                                                  S::mul(S::add(v, 81.0), 0.016949152542372881)))))));

which seems correct.

[thorstenhater]

For reference; before #2432 And afterwards

[thorstenhater]

Case closed