Disagreement between Crank-Nicolson, Backward Euler, and CVODE

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

ajpc6d
Posts: 36
Joined: Mon Oct 24, 2022 6:33 pm

Disagreement between Crank-Nicolson, Backward Euler, and CVODE

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

# 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']
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)

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)
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)
(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
Posts: 36
Joined: Mon Oct 24, 2022 6:33 pm

Re: Disagreement between Crank-Nicolson, Backward Euler, and CVODE

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
Posts: 36
Joined: Mon Oct 24, 2022 6:33 pm

Re: Disagreement between Crank-Nicolson, Backward Euler, and CVODE

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..!