cm * dV/dt - I_inj = - sum(i_ion)
(ignoring axial currents)
The following code is a minimal example of doing this whereby I computed the left hand side of the equation in two different ways:
Code: Select all
from neuron import h
import numpy as np
import pylab as pl
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 currents
cell.insert('pas')
# recordings
ipas = h.Vector()
ipas.record(cell(.5).pas._ref_i)
v = h.Vector()
v.record(cell(.5)._ref_v)
icap = h.Vector()
icap.record(cell(.5)._ref_i_cap)
i_inj = h.Vector()
i_inj.record(stim._ref_i)
# run simulation
h.tstop = 700
h.dt = 0.025
h.init()
h.run()
# change to arrays
v = np.array(v)
i_inj = np.array(i_inj)
icap = np.array(icap)
ipas = np.array(ipas)
t = np.arange(0, h.tstop+h.dt, h.dt)
# change units
Cm = cell.cm * cell(.5).area() * 1e-8
icap_from_v = np.concatenate((np.array([0]), np.diff(v)/np.diff(t))) * Cm
icap = 1e3 * icap * cell(.5).area() * 1e-8
ipas = ipas * 1e3 * cell(.5).area() * 1e-8
i_inj *= 1e-3
# plot
# the three traces should be the same
pl.figure()
pl.plot(t, icap_from_v - i_inj, 'k')
pl.plot(t, (icap) - i_inj, 'r')
pl.plot(t[:-1], -1 * ipas[1:] , 'y')
pl.show()
1) Is the implementation correct?
2) Why do I have to shift the passive current by 1?
3) Why is the manual computation of icap (icap_from_v) a little different from the other two traces?
Furthermore if I include another current and add it to ipas, then the deviation of the right hand side of the equation to the left hand side gets stronger. Why does this happen?