www.neuron.yale.edu

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:

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()

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?

Posted: **Thu Feb 04, 2016 5:01 pm**

Well, it's a single compartment model, not a cable, but it's a good place to start.

Your code is too clever by half. Python syntax is truly wonderful etc., but notational elegance from the programmer's perspective can obscure the discussion of numerical methods. Can you write something that makes the numerical integration method transparently obvious even to those who aren't familiar with Python? That means not taking advantage of numpy's marvellous methods. Also use a "for" loop. Painful, I know, but then we can focus on numerical methods rather than how to interpret numpy's documentation.

Posted: **Fri Feb 05, 2016 4:23 am**

I think you mean this place where I used numpy magic:

Code: Select all

icap_from_v = np.concatenate((np.array([0]), np.diff(v)/np.diff(t))) * Cm

In a for loop this would be:

Code: Select all

icap_from_v = np.zeros(len(v))
for i in range(1, len(v)):
    icap_from_v[i] = (v[i] - v[i-1]) / (t[i] - t[i-1]) * Cm

But I don't know how to make the numerical integration obvious. As I thought that the numerical integration happens in NEURON's run() function.

Posted: **Mon Feb 08, 2016 5:15 pm**

My apology. Looking back at your original post I see that I misinterpreted your original intent. Having read
"I wanted to check the cable equation implemented in NEURON which says that:
cm * dV/dt - I_inj = - sum(i_ion)"
and
"I computed the left hand side of the equation in two different ways"
I jumped to the conclusion that you intended to integrate the differential equation with your own Python code--since that's what I would have done--and took only the briefest glance at the actual code you provided. In retrospect, I see that you were just using NEURON-generated simulation results to calculate your own estimate of membrane capacitive current, and comparing that to the recording of i_cap.

Now to your questions.

1) Is the implementation correct?

The implementation of the model cell, IClamp, and vector recording are correct.

2) Why do I have to shift the passive current by 1?

You don't. In fact, you shouldn't. The lag between injected current and membrane capacitive current on the one hand, and membrane ionic current on the other, is inevitable, and in a fixed time step simulation it can look like an "off by 1 time step" result. See the discussion in my most recent post to your other thread. This is simply a manifestation of discretization error--an error that is inherent in trying to use a process that is discrete in time as an approximation to a process that is continuous in time. It goes away as you decrease dt (more about this later).

"Well, why can't I just make it vanish by shifting the passive current recording one step earlier?"

That may work in this case, but there's no guarantee that it will work if active currents are present--could make results worse. If you want more accuracy, reduce dt or change to an integration method that has better precision. NEURON's adaptive integrators work very nicely, but you can probably get away by simply cutting dt to 0.005 ms.

"Using adaptive integration would complicate my attempt to verify the value NEURON returns for i_cap."

Yes, you'll have to verify results in a different way. Or you could still use fixed dt but with a smaller time step. Let me suggest what you might try.

First, clean up the simulation by fixing its very poor initialization. Your model's resting potential is -70 mV (same as value of e_pas), so why allow NEURON to initialize it to -65 mV? Insert
h.v_init = -70
before you call h.run(). (By the way, you don't have to call h.init()--h.run() does that for you.)

Now you'll get simulation results that aren't dominated by a bad initialization artifact for the first 10-20 ms or so.

Next, change the stimulus and run parameters so you can see the details of simulation execution at a higher resolution. Make
stim.delay 0.2 ms
stim.dur 0.3 ms
h.tstop 1 ms

Then in a new graph plot the unscaled stimulus current (IClamp[0].i), which is already in nA, plus i_cap and i_pas after scaling them by area(0.5)*1e-2 to convert their units (mA/cm2) to nA.

Finally, in another new graph plot (the unscaled IClamp[0].i - the scaled i_cap), and also plot the scaled i_pas. Run a simulation with dt = 0.025, then change dt to 0.005 and run another simulation. Nicer--look how much closer capacitive current is to the stimulus current when the stimulus turns on, and how much closer i_pas is to the stimulus current minus the capacitive current.

www.neuron.yale.edu

Check the cable equation

Check the cable equation

Re: Check the cable equation

Re: Check the cable equation

Re: Check the cable equation