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?