Disagreement between Crank-Nicolson, Backward Euler, and CVODE

[ajpc6d]

I've been continuing to try to understand my inconsistent results (as outlined in previous posts on this forum). My most recent attempt involved taking the stimsphere.py file in stimsphere.zip, hosted here, and adapting it to my scenario.

In the code below, I have tried to retain as much of the original stimsphere code as possible, to illustrate proper functioning, and also added in a highly reduced representation of a complex neuron shape stimulated by electric field.

Three of NEURON's integration methods (Crank-Nicolson, backward Euler, and CVODE) give nearly identical results for the sphere cell scenario. Crank-Nicholson diverges sharply for the scenario of the complex cell.

What is going on here? I know of no mathematical property of the Crank-Nicolson method which prohibits its application for a 'squiggly shape' scenario. (And, since better minds than mine chose to include it in NEURON, I suspect no such property exists.)

The python code:

Code: Select all

from neuron import h, gui
import matplotlib.pyplot as plt
import numpy as np
h.load_file("stdrun.hoc")
h.load_file("rig.ses")

# use interp_xyz at the SECTION level
def interp_xyz(section, origin=[0,0,0]):
    definition = section.nseg
    nn = 0
    sec = section  # h.cas()

    # get the pt3d data for each section
    nn = int(sec.n3d())  # nn is the number of 3-D data points defining each section
    xx = h.Vector(nn)
    yy = h.Vector(nn)
    zz = h.Vector(nn)
    length = h.Vector(nn)
    for ii in range(nn):
        xx[ii] = sec.x3d(ii)
        yy[ii] = sec.y3d(ii)
        zz[ii] = sec.z3d(ii)
        length[ii] = sec.arc3d(ii)

    # h.Vector().interpolate() requires normalized section lengths
    length.div(length[nn - 1])

    # creating the destination vector
    range_var = h.Vector(definition + 2)
    range_var.indgen(1 / definition)
   
    range_var.sub(1 / (2 * definition))
    range_var[0] = 0
    range_var[definition + 1] = 1

    # likewise, the Vectors of interpolated pt3d data should be extended to include endpoints
    xint = h.Vector(definition + 2)
    yint = h.Vector(definition + 2)
    zint = h.Vector(definition + 2)

    # meaning: interpolate the (length, xyz) data onto the range_var scale
    xint.interpolate(range_var, length, xx)
    yint.interpolate(range_var, length, yy)
    zint.interpolate(range_var, length, zz)

    x = [i - origin[0] for i in xint.c(1,definition)]
    y = [j - origin[1] for j in yint.c(1,definition)]
    z = [k - origin[2] for k in zint.c(1,definition)]
    return list(zip(x, y, z))

def makerealdend(nseg):
    data = {'x': [-8.96, -10.42, -11.85, -13.67, -13.41, -14.06, -13.67,
                  -13.54, -13.63, -13.94, -12.34, -11.98, -12.2, -12.14,
                  -10.57, -9.16, -7.46, -3.98, -1.88, -0.19, 2.52, 4.12, 4.08],
            'y': [-2.26, -3.05, -4.68, -6.77, -8, -8.76, -10.25, -11.23,
                  -12.64, -13.63, -14.19, -15.24, -17.89, -19.94, -21.28,
                  -21.9, -21.55, -19.79, -18.59, -18.23, -17.05, -16.1, -16.04],
            'z': [-0.65, -1.1, -1.34, -1.37, -3.4, -3.72, -3.94, -4.3, -6.69,
                  -7.4, -8.15, -8.89, -9.6, -10.52, -12.96, -14.86, -15.57,
                  -15.99, -16.61, -17.56, -19.05, -20.11, -21.63],
            'diam': 0.95       
            }
    
    n3d = len(data['x'])
    dendrite = h.Section(name="dendrite")

    for n in range(n3d):
        xx = data['x'][n]
        yy = data['y'][n]
        zz = data['z'][n]
        dd = data['diam']
        dendrite.pt3dadd(xx,yy,zz,dd)
    for sec in h.allsec():
        h.hh_li_v1.insert(sec)
        sec.ek = -90
        sec.ena = 60
        sec.cm = 0.75
        sec.Ra = 100
        h.extracellular.insert(sec)
        sec.nseg = nseg
    return [dendrite]

def makespherecell(diam, n3dpts, nseg):
    somaR = h.Section(name='somaR')
    somaL = h.Section(name='somaL')
    somaR.connect(somaL(0),0)
    
    for i in range(n3dpts):
        theta = h.PI*0.5*(1 - i / (n3dpts - 1))
        xx = 0.5*diam*h.cos(theta)
        dd = diam*h.sin(theta)
        somaR.pt3dadd(xx,0,0,dd) # remove the 'i'
        somaL.pt3dadd(-xx,0,0,dd)

    for sec in h.allsec():
        h.hh_li_v1.insert(sec)
        sec.ek = -90
        sec.ena = 60
        sec.cm = 0.75
        sec.Ra = 100
        # h.xtra.insert(sec)
        # h.hh.insert(sec)
        h.extracellular.insert(sec)
        sec.nseg = nseg

    return [somaR, somaL]

def spherecalc(sec, seg, dedx=1):
    xloc = (sec.diam3d(0)/2)*seg.x # distance of seg's node from 0 end
    if sec.x3d(1) < 0: 
        xloc = -xloc # left hemisphere
    rx = 1e-3 * (dedx * xloc)
    return rx

def realcalc(pts, j):
    rx = pts[j][1] * 1e-3
    if abs(rx) < 1e-9:
        rx = 1e-9 * np.sign(rx)
    return rx


##########
# This now works with hh/xtra, hh_sphere, and hh_li_v1
# WHY doesn't it work in the big script?
##########

# cell = makespherecell(5, 25, 9)
cell = makerealdend(9)
amp = 1e3 * h.PI * 2.8 # near threshold for sphere cell
amp = 8e3 # near threshold (maybe?) for real dendrite
h.celsius = 37

t_sim = 20
h.tstop = t_sim
dpts = 2000
stimpts = 150
padpts = (dpts - stimpts) // 4
rempts = dpts - padpts - stimpts
wave = np.ones(stimpts)
stim = np.pad(wave, (padpts,rempts), constant_values=0)
svec = h.Vector(stim).mul(amp)
tvec = h.Vector(np.linspace(0,t_sim,dpts))

svec.play(h._ref_ex_hh_li_v1, tvec, True)

for sec in h.allsec():
    pt3ds = interp_xyz(sec)
    j = -1
    for seg in sec:
        j+=1
        seg.hh_li_v1._ref_Vs = seg._ref_e_extracellular
        
        seg.hh_li_v1.rx = spherecalc(sec, seg)
        seg.hh_li_v1.rx = realcalc(pt3ds, j)
       

# h.dt = 1e-3
rdt = 0.025
t = h.Vector().record(h._ref_t, rdt)
v = h.Vector().record(cell[0](0.5)._ref_v, rdt)

h.secondorder = 0
h.finitialize(-65)
h.continuerun(t_sim)
print(f"the chosen solver is {h.CVode().current_method()}")
lns = plt.plot(t,v, label='backward Euler, soma')

h.secondorder = 2
h.finitialize(-65)
h.continuerun(t_sim)
print(f"the chosen solver is {h.CVode().current_method()}")
lns += plt.plot(t,v, label='Crank-Nicolson, soma')

cvode = h.CVode()
cvode.active(True)
cvode.use_daspk(True)
cvode.dae_init_dteps(1e-9)
h.finitialize(-65)
h.continuerun(t_sim)
print(f"the chosen solver is {h.CVode().current_method()}")
lns += plt.plot(t,v, label='CVODE, soma')

lns += plt.twinx().plot(tvec, svec, label='RF', color='r', linestyle=':')

lbls = [x.get_label() for x in lns]

plt.legend(lns, lbls)
plt.show()

The NMODL code:

Code: Select all

: A custom Hodgkin-Huxley-style ion channel built on the ODEs found in Li et al 2015

NEURON {
    SUFFIX hh_li_v1
    USEION k READ ek WRITE ik
    USEION na READ ena WRITE ina
    NONSPECIFIC_CURRENT il:, i
    RANGE gnabar, ena, gna
    RANGE gkbar, ek, gk
    RANGE glbar, el, gl
    RANGE rx
    GLOBAL ex
    POINTER Vs
}

UNITS {
    (mV) = (millivolt)
    (mA) = (milliamp)
    (mS) = (millisiemens)
    (S) = (siemens)
    (nW) = (nanowatts)
    (uA) = (microamp)
    (uF) = (microfarad)
    (um) = (micrometer)
}

PARAMETER {
    gkbar = 40 (mS/cm2) <0,1e9>
    gnabar = 150 (mS/cm2) <0,1e9>
    glbar = 0.033 (mS/cm2) <0,1e9>

    el = -70 (mV)
    :ek = -90 (mV)
    :ena = 60 (mV)
}

ASSIGNED {
    celsius (degC) 
    phi
    dt (ms)

    v (mV)
    Vs (mV)
    ex (mV/um)
    rx (um)

    ena (mV)
    ek (mV)

    gna (mS/cm2)
    gk (mS/cm2)
    gl (mS/cm2)
    
    ina (mA/cm2)
    ik (mA/cm2)
    il (mA/cm2)
}

STATE { n m h}

BREAKPOINT {
    SOLVE states METHOD cnexp
    : The ionic conductances w.r.t time
    gna = gnabar*m^3*h
    gk = gkbar*n
    gl = glbar

    : The ionic currents
    ik = gk*(v - ek)*(1e-3)
    ina = gna*(v - ena)*(1e-3)
    il = gl*(v - el)*(1e-3)

    : The voltage which will be pointed into _e_extracellular
    Vs = ex*rx:*(1e6)

}

DERIVATIVE states {
    n' = alpha_n(v, phi)*(1 - n) - beta_n(v, phi)*n
    m' = alpha_m(v, phi)*(1 - m) - beta_m(v, phi)*m
    h' = (1 / (1 + exp((v + 60) / 6.2 (mV))) - h)*(alpha_h(v, phi) + beta_h(v, phi))
}

INITIAL {
    phi = 2.3^((celsius - 23)/10 (degC))
    n = alpha_n(v, phi) / (alpha_n(v,phi) + beta_n(v,phi))
    m = alpha_m(v, phi) / (alpha_m(v,phi) + beta_m(v,phi))
    h = alpha_h(v, phi) / (alpha_h(v,phi) + beta_h(v,phi))
    Vs = ex*rx:*(1e6)
}

FUNCTION alpha_n(v (mV), phi) (/ms) {
    LOCAL x
    UNITSOFF
    x = 1 - exp(-(v - 30)/9 (mV))
    if (fabs(x) < 1e-6) {
        alpha_n = phi*0.01*x / (1 - exp(-x))
    }
    else {
        alpha_n = phi*(0.01*(v - 30)) / x
    }
    UNITSON
}

FUNCTION beta_n(v (mV), phi) (/ms) {
    LOCAL x
    UNITSOFF
    x = 1 - exp((v - 30)/9 (mV))
    if (fabs(x) < 1e-6) {
        beta_n = -phi*0.002*x / (1 - exp(-x))
    }
    else {
        beta_n = -phi*(0.002*(v - 30)) / x
    }
    :beta_n = phi*0.125*exp(-(v + 65) / 80)
    UNITSON
}

FUNCTION alpha_m(v (mV), phi) (/ms) {
    LOCAL x
    UNITSOFF
    x = 1 - exp(-(v + 30)/8 (mV))
    if (fabs(x) < 1e-6) {
        alpha_m = phi*0.182*x / (1 - exp(-x))
    }
    else {
        alpha_m = phi*(0.182*(v + 30)) / x
    }
    UNITSON
}

FUNCTION beta_m(v (mV), phi) (/ms) {
    LOCAL x
    UNITSOFF
    x = 1 - exp((v + 30)/8 (mV))
    if (fabs(x) < 1e-6) {
        beta_m = -phi*0.124*x / (1 - exp(-x))
    }
    else {
        beta_m = -phi*(0.124*(v + 30)) / x
    }
    :beta_m =  phi * 4 * exp(-(v+65)/18)
    UNITSON
}

FUNCTION alpha_h(v (mV), phi) (/ms) {
    LOCAL x
    UNITSOFF
    x = 1 - exp(-(v + 45)/6 (mV))
    if (fabs(x) < 1e-6) {
        alpha_h = phi*0.028*x / (1 - exp(-x))
    }
    else {
        alpha_h = phi*(0.028*(v + 45)) / x
    }
    :alpha_h = phi*.07 * exp(-(v+65)/20)
    UNITSON
}

FUNCTION beta_h(v (mV), phi) (/ms) {
    LOCAL x
    UNITSOFF
    x = 1 - exp((v + 70)/6 (mV))
    if (fabs(x) < 1e-6) {
        beta_h = -phi*0.0091*x / (1 - exp(-x))
    }
    else {
        beta_h = -phi*(0.0091*(v + 70)) / x
    }
    :beta_h = phi / (exp(-(v+35)/10) + 1)
    UNITSON
}

[ajpc6d]

I found this post which indicates CN ought not to be used when using voltage clamps. Is it possible that this same logic applies to extracellular stimulation?

[ajpc6d]

From The NEURON Book (pp.87, emphasis added):

Although the Crank–Nicholson method is formally stable, models with stiff equations require small del_t to avoid numerical oscillations (Fig. 4.8). It is unusable in the presence of voltage clamps, extracellular mechanisms, or linear circuits, since the solution of algebraic equations gives results with large numerical oscillations.

So, I guess this answers the question. I don't see this tidbit anywhere in the 'readthedocs' documentation. Even The NEURON Book doesn't make it very clear, in my opinion; pp.87 is the penultimate page of the summary of the chapter, after the topic has already been discussed in detail.
Then again, since the question hasn't been asked before, maybe I'm the problem..!