What does cnexp do?

NMODL and the Channel Builder.
Post Reply
cafischer

What does cnexp do?

Post by cafischer »

cafischer wrote:I could not find a description or reference of the cnexp method in the NEURON documentation and the NEURON book. I would like to know which method it is (equations, reference).
All the source code is at your disposal, but here is a brief summary.

I like to think of cnexp as more of a "directive" for nrnivmodl, the program that translates NMODL code into executable code, but it can also be regarded as the name of a bundle of integrators that are suitable for HH-style (linear ODE) descriptions of channel gating.

In a fixed time step simulation with secondorder == 0, it uses the analytical solutions for the gating variable ODEs to advance their values to t+dt, after which the resulting ionic conductances and currents are used to advance membrane potential to t+dt by numerical integration of the cable equation itself. This is absolutely stable. Error is proportional to dt, and reasonable values of dt produce results that are sufficiently accurate for most purposes. Moore and Ramon did something similar for one of the methods they evaluated in
Moore, J. W. & Ramon, F.
On numerical integration of the Hodgkin and Huxley equations for a membrane action potential.
J Theor Biol, 1974, 45, 249-273

If secondorder==1 or 2, then a staggered Crank-Nicholson method (see Chapter 4 of The NEURON Book) is used to integrate the gating variables. Ionic conductances are second order correct at t+dt/2; ionic currents have only first order accuracy unless secondorder==2, in which case they are second order correct at t-dt/2.

If adaptive integration is used (i.e. if cvode is active), then one of NEURON's adaptive (fully implicit, variable order, variable time step) integrators is used instead of implicit Euler.
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: What does cnexp do?

Post by ted »

Nuts--in replying to your message I seem to have accidentally deleted your original post. Big mistake. Here's what I could salvage:

I create a cell, inserted the transient sodium channel. I stimulated with a step current and measured the current and voltage trace:
. . .
I used the voltage trace to solve the ODEs of the channel myself using the Rosenbrock Euler method and then compared the current trace from the NEURON simulation and the direct integration:

Code: Select all

from neuron import h
import numpy as np
import csv
h.load_file("stdrun.hoc")  # load NEURON libraries
h("""cvode.active(0)""")  # unvariable time step in NEURON

# create simple cell model
cell = h.Section(name='soma')
cell.cm = 1
cell.Ra = 100
cell.diam = 100
cell.L = 100
cell.nseg = 1

# inject current
stim = h.IClamp(cell(0.5))
stim.amp = 1
stim.delay = 100
stim.dur = 500

# insert ionic current
cell.insert('nat')
for seg in cell:
    seg.nat.gbar = 1

# recordings
ina = h.Vector()
ina.record(cell(.5).nat._ref_ina)

v = h.Vector()
v.record(cell(.5)._ref_v)

t = h.Vector()
t.record(h._ref_t)

# run simulation
dt = 0.025
h.tstop = 700
h.steps_per_ms = 1/dt
h.dt = dt
h.run()

# change to arrays
v = np.array(v)
t = np.array(t)
ina = np.array(ina)

# save ina
with open('nat_neuronsim.npy', 'w') as f:
    np.save(f, ina)

# save t, v
header = ['t', 'v']
data_new = np.c_[t, v]
with open('./nat_neuronsim.csv', 'w') as csvoutput:
    writer = csv.writer(csvoutput, lineterminator='\n')
    writer.writerow(header)
    writer.writerows(data_new)
I used the voltage trace to solve the ODEs of the channel myself using the Rosenbrock Euler method and then compared the current trace from the NEURON simulation and the direct integration:

Code: Select all

from __future__ import division
import numpy as np
import pylab as pl
import pandas as pd


# functions
def rates(v, Ra, Rb, qa, tha, celsius, temp):
    a = Ra * qa * efun((tha - v)/qa)
    b = Rb * qa * efun((v - tha)/qa)

    tadj = q10**((celsius - temp)/10)

    mtau = 1/tadj/(a+b)
    minf = a/(a+b)

    a = Rd * qi * efun((thi1 - v)/qi)
    b = Rg * qi * efun((v - thi2)/qi)

    htau = 1/tadj/(a+b)
    hinf = 1/(1+np.exp((v-thinf)/qinf))

    return mtau, minf, htau, hinf


def efun(x):
    if np.abs(x) < 1e-6:
        y = 1 - x/2
    else:
        y = x/(np.exp(x) - 1)
    return y


def slope_intercept_mh(i, vshift, Ra, Rb, qa, tha, celsius, temp):
    mtau, minf, htau, hinf = rates(v[i]+vshift, Ra, Rb, qa, tha, celsius, temp)
    return -1/mtau, -1/htau, minf/mtau, hinf/htau


def phi(z):
    return (np.exp(z) - 1) / z

# parameters from simulation
celsius = 36
ena = 83

# load data
data = pd.read_csv('./nat_neuronsim.csv').convert_objects(convert_numeric=True)
v = np.array(data.v)
t = np.array(data.t)
dt = t[1] - t[0]

# channel parameter
vshift = 0  # (mV)
gbar = 1  # (pS/um2)
tha = -35  # (mV)
qa = 9  # (mV)
Ra = 0.182  # (/ms)
Rb = 0.124  # (/ms)
thi1 = -50  # (mV)
thi2 = -75  # (mV)
qi = 5  # (mV)
thinf = -65  # (mV)
qinf = 6.2  # (mV)
Rg = 0.0091  # (/ms)
Rd = 0.024  # (/ms)
temp = 23  # (degC)
q10 = 3
vin = -120  # (mV)
vax = 100  # (mV)

# initialization
m = np.zeros(len(t))
h = np.zeros(len(t))
ina = np.zeros(len(t))
mtau, minf, htau, hinf = rates(v[0]+vshift, Ra, Rb, qa, tha, celsius, temp)
m[0] = minf
h[0] = hinf
ina[0] = 1e-4 * gbar*m[0]**3*h[0] * (v[0] - ena)

m_slope = np.zeros(len(t))
h_slope = np.zeros(len(t))
m_intercept = np.zeros(len(t))
h_intercept = np.zeros(len(t))
m_slope[0], h_slope[0], m_intercept[0], h_intercept[0] = slope_intercept_mh(0, vshift, Ra, Rb, qa, tha, celsius, temp)
for i in range(1, len(t)): 
    m_slope[i], h_slope[i], m_intercept[i], h_intercept[i] = slope_intercept_mh(i-1, vshift, Ra, Rb, qa, tha, celsius, temp)  # current was flowing in response to v from the previous time step

# integration (exponential Euler)
for i in range(1, len(t)):
    m[i] = dt * phi(m_slope[i]*dt) * m_intercept[i-1] + m[i-1] * np.exp(m_slope[i]*dt)
    h[i] = dt * phi(h_slope[i]*dt) * h_intercept[i-1] + h[i-1] * np.exp(h_slope[i]*dt)

    # compute current
    ina[i] = 1e-4 * gbar*m[i]**3*h[i] * (v[i-1] - ena)   # current was flowing in response to v from the previous time step


# plot
with open('nat_neuronsim.npy', 'r') as f:
    ina_neuron = np.load(f)

pl.figure()
pl.plot(t, ina_neuron, 'r', linewidth=1.5, label='NEURON \nsimulation')
pl.plot(t, ina, 'k', linewidth=1.5, label='Direct \nIntegration')
pl.xlabel('$Time (ms)$', fontsize=18)
pl.ylabel('$Current (pS/\mu m^2)$', fontsize=18)
pl.legend(loc='lower right', fontsize=18)
pl.show()
I found that in the very first step there is a difference between the current traces. What could be the cause?
Also the rest of the current traces a quite different. Did I use the wrong integration method or is there another reason?
cafischer

Re: What does cnexp do?

Post by cafischer »

I would like to learn more about the the integtrators that are used within cnexp in a fixed time step simulation.
So if the user does not set manually secondorder = 0, then gating variables of ion channels are always computed analytically (using symbolic manipulation software)?
When won't they be computed analytically? And what will be the method (does it then fall back to the staggered Crank-Nicholson method)?
How can the user determine whether the gating variables were computed analytically or by numerical integration?
Could you point out to me the file in the src code doing this?
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: What does cnexp do?

Post by ted »

cafischer wrote:if the user does not set manually secondorder = 0, then gating variables of ion channels are always computed analytically
True for fixed step simulation of a channel described by NMODL code that uses hh-style linear ODEs to describe gating variable dynamics.
using symbolic manipulation software?
Not as such. Those ODEs are so stereotyped that parsing them is relatively straightforward.
When won't they be computed analytically?
If any of the following are true:
--the ODEs are nonlinear
--the channel description is expressed in kinetic scheme notation as a family of state transitions
--adaptive integration is used
does it then fall back to the staggered Crank-Nicholson method?
Staggered Crank-Nicholson has nothing to do with it.
How can the user determine whether the gating variables were computed analytically or by numerical integration?
See the reply to "when won't they be computed analytically"
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: What does cnexp do?

Post by hines »

Can you send me a zip file with the necessary code. The fragments above are missing nat.mod file.
michael dot hines at yale dot edu
Post Reply