### Check the cable equation

Posted:

**Thu Feb 04, 2016 12:20 pm**I wanted to check the cable equation implemented in NEURON which says that:

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:
My question is:

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?

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?