nmodl + python

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

Post Reply
kadir

nmodl + python

Post by kadir » Thu Mar 24, 2011 8:42 am

I'm using the nmodl mechanism:

Code: Select all

NEURON {
	POINT_PROCESS SpikeOut
	RANGE thresh, refrac, vrefrac, grefrac
	NONSPECIFIC_CURRENT i
}

PARAMETER {
	thresh = -55 (millivolt)
	refrac = 1 (ms)
	vrefrac = -65 (millivolt)
	grefrac = 100 (microsiemens) : clamp to vrefrac
}

ASSIGNED {
	i (nanoamp)
	v (millivolt)
	g (microsiemens)
}

INITIAL {
	net_send(0, 3)
	g = 0
}

: up to here, executed once

BREAKPOINT {
	i = g*(v - vrefrac)
}

NET_RECEIVE(w) {
	if (flag == 1) {
		v = vrefrac
		g = grefrac
		net_event (t) : tells that an event, spike, occured at time t to other neurons.
		net_send(refrac, 2)
	}else if (flag == 2) {
		g = 0
	}else if (flag == 3) {
		WATCH (v > thresh) 1 : if v >= thresh -> net_send(0,1) -> NET_RECEIVE(w=1)  watch is being
                : executed in each time step, you :dont have to call it.
	}	
}
I designed a template:

Code: Select all

begintemplate purkinjecell

public soma, threshold, undershootpotential, spikegenerator

create soma

objref spikegenerator

proc init() {

    // note the NEURON units

    diameter = 25.5 // micrometer
    
    //soma.cm = 0.64 // microfarad / cm^2

    //soma.cm = 5 // microfarad / cm^2. for 1980 paper 

    soma.cm = 0.69 // microfarad / cm^2

    // specificmembraneresistance = 120.2 * 10 ^ 3 // ohm . cm^2

    //specificmembraneresistance = 40 * 10 ^ 3 // ohm . cm^2. for 1980 paper

    specificmembraneresistance = 10 * 10 ^ 3 // ohm . cm^2

    surfacearea = 40821 // micrometer ^ 2

    length = surfacearea / (PI * diameter) // micrometer

    membraneresistance = specificmembraneresistance / (surfacearea * 10 ^ -8)

    soma {

    diam = diameter

    L = length

    // leakage mechanism

    insert pas

    g_pas = (1 / membraneresistance) / (surfacearea * 10 ^ -8) // S/cm^2

    //e_pas = -25 // adjust to get 40 hz

    //e_pas = 25 // adjust to get 40 hz. for 1980 paper

    e_pas = -54.9 // adjust to get 30 hz

    Ra = 141.9 // ohm.cm

    }

    spikegenerator = new SpikeOut(0.5)

    // as thresh is local variable, spikegenerator has the attribute thresh. 

    threshold = spikegenerator.thresh

    undershootpotential = spikegenerator.vrefrac  
   
}

endtemplate purkinjecell
I am using NEURON as a python module. When I generate a single cell from the template and inject current, I get spiking behavior (everything is fine).

However, when I generate 2 cells and stimulate the second cell, it doesn't fire. I checked that it has the nmodl mechanism, but somehow it doesn't work properly. Here's the script:

Code: Select all

from neuron import h

import numpy

from pylab import plot, show

h.load_file ("purkinje cell.hoc")

purkinjecellx = h.purkinjecell ()

purkinjecelly = h.purkinjecell ()

stimulus = h.IClamp (0.5, sec = purkinjecelly.soma)
  
stimulus.delay = 100

stimulus.dur = 300
  
stimulus.amp = 0.1    
  
time = h.Vector()
  
membranepotential = h.Vector()

time.record(h._ref_t)

membranepotential.record (purkinjecelly.soma(0.5)._ref_v)

h.load_file("stdrun.hoc")

h.init()

h.v_init = purkinjecelly.undershootpotential

h.tstop = 500

h.run()

time = numpy.array(time)

membranepotential = numpy.array(membranepotential)

plot(time, membranepotential)

show()
What can it be due to?

Thanks in advance

kadir

Re: nmodl + python

Post by kadir » Thu Mar 24, 2011 1:26 pm

I reproduced the problem in NEURON (without python) using the following script:

Code: Select all

load_file ("purkinje cell.hoc")

objectvar purkinjecellx, purkinjecelly, stimulus

purkinjecellx = new purkinjecell ()

purkinjecelly = new purkinjecell ()

purkinjecelly.soma stimulus = new IClamp (0.5)

access purkinjecelly.soma

stimulus.del = 100

stimulus.dur = 300
  
stimulus.amp = 0.1    

tstop = 500

ted
Site Admin
Posts: 5625
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: nmodl + python

Post by ted » Thu Mar 24, 2011 9:44 pm

The problem(s) are partly in your template and partly in your model cell's parameters. Let's start with the model parameters, which are the easiest to deal with.

The model cell has a very large surface area and a moderately large specific membrane conductance. Consequently its input resistance is quite low, so the stimulus current barely does anything to membrane potential. But it doesn't have to do much, because e_pas is -54.9 and spike threshold is -55, so the cell should fire tonically even in the absence of injected current.

Which is exactly what purkinjecellx does. It's twin, however, does not have its own SpikeOut instance, which is why purkinjecelly is silent. But that doesn't mean there aren't two SpikeOuts in this model--there are, but somehow purkinjecellx got both of them.

Writing templates is a complicated business, at least in hoc, with many details to attend to, and many opportunities to make mistakes. That's where the Network Builder suite of GUI tools comes in handy, because it can export a hoc file that contains properly written templates that define your cell classes.

I just now used it to create a definition of a single cell class, which I gave the name PCell, that has a soma with exactly the same anatomical and biophysical properties as your cell class. (For a couple of tutorials on the Network Builder, see http://www.neuron.yale.edu/neuron/docs)

Then I exported a hoc file and extracted from it a template that I saved to a file called netpcell.hoc. This template required just a bit of editing:
1. I changed the template name from PCell_PCell to PCell, just for convenience.
2. Since the Network Builder tools didn't provide any way for me to attach a SpikeOut to the soma, I had to add the following bits of code to the template. Specifically:
(a) At the top of the template, after all of the other "public" declarations, I inserted

Code: Select all

public spkout
objref spkout
(b) Inside proc synapses(), I inserted

Code: Select all

  soma spkout = new SpikeOut(0.5)
  spkout.thresh=-55
  spkout.refrac=1
  spkout.vrefrac=-65
  spkout.grefrac=100
Yes, I know SpikeOut isn't a synaptic mechanism. proc synapses() was just a convenient place to put some code that I wanted to be executed every time a new instance of the PCell class is created.

Then I revised your hoc file so it reads like so:

Code: Select all

load_file ("netpcell.hoc")

objref pcx, pcy, stim

pcx = new PCell()
pcy = new PCell()

pcy.soma stim = new IClamp (0.5)
stim.del = 100
stim.dur = 300
stim.amp = 0.1

tstop = 500
(short names mean less typing, and less typing means I make fewer mistakes).

So that's how I avoided wrestling with having to write a template. Sure,the PCell template generated by the Network Builder contains some unused stuff, but the code works. It's much easier to start with working code and then weed out the unused bits, than it is start from scratch and build something that works.

Post Reply