spike frequency vs current graph

Using the graphical user interface to build and exercise models. Includes customizing the GUI by writing a little bit of hoc or Python
MSN
Posts: 10
Joined: Mon May 23, 2011 5:38 pm

spike frequency vs current graph

Post by MSN »

Dear Seniors

I am a new user of Neuron.
I use a HH neuron and inject constant current into the dendrite (using GUI). For different values of current there is different spike frequency.
Now I want to plot the spike frequency vs current graph.
Can anybody help me, how can I plot it ?

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

Re: spike frequency vs current graph

Post by ted »

In pseudocode

Code: Select all

Create two Vectors called stim and freq.
For each stimulus amplitude
  append the amplitude to stim
  append the corresponding frequency to freq
Create an instance of the Graph class called g
Use the Vector class's plot method to plot stim vs. freq in g
See the Programmer's Reference for documentation of the Vector and Graph classes.

Use the Graph's "View = plot" (see the first few paragraphs at http://www.neuron.yale.edu/neuron/stati ... graph.html) to automatically adjust the axes of the graph.

Use the tricks described in the NEURON Forum's "NEURON hacks" area under
How to use hoc to plot a variable vs. time
to control position of the graph, and for more precise axis scaling.
MSN
Posts: 10
Joined: Mon May 23, 2011 5:38 pm

Re: spike frequency vs current graph

Post by MSN »

dear Ted

i think i did not explain my query clearly.
i again tried to solve the problem but remain unsucessful.
Here i explain it in more detail.

Here is my session file. (i develop it using gui)

Code: Select all

{load_file("nrngui.hoc")}
objectvar save_window_, rvp_
objectvar scene_vector_[5]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{pwman_place(0,0,0)}

//Begin CellBuild[0]
{
load_file("celbild.hoc", "CellBuild")
}
{ocbox_ = new CellBuild(1)}
{object_push(ocbox_)}
{
version(5.7)
continuous = 1
}
{object_push(topol)}
{
first = 0
slist.remove_all()
sname = "axon"
objref tobj
}
{
tobj = new CellBuildSection("soma",0, 0, tobj, 0) slist.append(tobj)
  tobj.position(0,0,15,0) tobj.lx=7.69231 tobj.ly=-15.2767 tobj.i3d=0
tobj = new CellBuildSection("apical",0, 0, tobj, 1) slist.append(tobj)
tobj.parent=slist.object(0)
  tobj.position(15,0,105,0) tobj.lx=85.4251 tobj.ly=-14.0621 tobj.i3d=0
tobj = new CellBuildSection("axon",0, 0, tobj, 0) slist.append(tobj)
tobj.parent=slist.object(0)
  tobj.position(0,0,-150,0) tobj.lx=-82.1862 tobj.ly=-11.6329 tobj.i3d=0
tobj = new CellBuildSection("basilar",0, 0, tobj, 0) slist.append(tobj)
tobj.parent=slist.object(0)
  tobj.position(0,0,-105,60) tobj.lx=-37.247 tobj.ly=46.6667 tobj.i3d=0
all_init()
}
for i=0, slist.count-1 {slist.object(i).rdses()}
{object_pop()}
{
}
{object_push(subsets)}
{first = 0}
{ tobj = snlist.object(0)}
{tobj = new SNList("has_HH") snlist.append(tobj)}
for i=0,1 tobj.add(bild.topol.slist.object(fscan()))
0
2
{tobj = new SNList("no_HH") snlist.append(tobj)}
for i=0,1 tobj.add(bild.topol.slist.object(fscan()))
1
3
{consist()}
{object_pop()}
{
}
{object_push(geom)}
{
first = 0
tobj = new GeoSpec(7)
tobj.value = 0.1
bild.subsets.snlist.object(0).geo.append(tobj)
tobj = new GeoSpec(2)
tobj.value = 30
bild.topol.slist.object(0).geo.append(tobj)
tobj = new GeoSpec(3)
tobj.value = 30
bild.topol.slist.object(0).geo.append(tobj)
tobj = new GeoSpec(2)
tobj.value = 1000
bild.topol.slist.object(1).geo.append(tobj)
tobj = new GeoSpec(3)
tobj.value = 1
bild.topol.slist.object(1).geo.append(tobj)
tobj = new GeoSpec(2)
tobj.value = 1000
bild.topol.slist.object(2).geo.append(tobj)
tobj = new GeoSpec(3)
tobj.value = 1
bild.topol.slist.object(2).geo.append(tobj)
tobj = new GeoSpec(2)
tobj.value = 1000
bild.topol.slist.object(3).geo.append(tobj)
tobj = new GeoSpec(3)
tobj.value = 2
bild.topol.slist.object(3).geo.append(tobj)
set_default()
}
{object_pop()}
{
}
{object_push(memb)}
{first=0}
{
tobj = new FakeMechStan(1)
tobj.value = 1
tobj.set_default()
tobj = new MStanWrap(tobj, 0)
bild.subsets.snlist.object(0).ml.append(tobj)
}
{
tobj = new FakeMechStan(0)
tobj.value = 100
tobj.set_default()
tobj = new MStanWrap(tobj, 0)
bild.subsets.snlist.object(0).ml.append(tobj)
}
{
tobj = new MechanismStandard("hh")
tobj.set("gnabar_hh", 0.12, 0)
tobj.set("gkbar_hh", 0.036, 0)
tobj.set("gl_hh", 0.0003, 0)
tobj.set("el_hh", -54.3, 0)
tobj = new MStanWrap(tobj, 1)
bild.subsets.snlist.object(1).ml.append(tobj)
}
{
tobj = new MechanismStandard("pas")
tobj.set("g_pas", 0.0002, 0)
tobj.set("e_pas", -65, 0)
tobj = new MStanWrap(tobj, 1)
bild.subsets.snlist.object(2).ml.append(tobj)
}
{object_pop()}
{
}
{object_push(manage)}
{
first = 0
classname = "Cell"
etop=1 esub=1 egeom=1 emem=1
itop=1 isub=0 igeom=0 imem=0
bild.topol.names_off = 0
bild.topol.circles_off = 0
output_index = 0  output_x = 1
thresh = 10
}
{object_pop()}
{
cexport()
}
{object_pop()}
{
save_window_=ocbox_.gtopol
save_window_.size(-200,200,-150,150)
scene_vector_[2] = save_window_
ocbox_.gtopol = save_window_
save_window_.save_name("ocbox_.gtopol")
}
{
ocbox_.map("CellBuild[0]", 978, 42, 774, 755.1)
}
objref ocbox_
//End CellBuild[0]


//Begin PointProcessManager
{
load_file("pointman.hoc")
}
{
basilar ocbox_ = new PointProcessManager(0)
}
{object_push(ocbox_)}
{
mt.select("IClamp") i = mt.selected()
ms[i] = new MechanismStandard("IClamp")
ms[i].set("del", 5, 0)
ms[i].set("dur", 50, 0)
ms[i].set("amp", 1, 0)
mt.select("AlphaSynapse") i = mt.selected()
ms[i] = new MechanismStandard("AlphaSynapse")
ms[i].set("onset", 0.5, 0)
ms[i].set("tau", 0.1, 0)
ms[i].set("gmax", 0.05, 0)
ms[i].set("e", 0, 0)
mt.select("SEClamp") i = mt.selected()
ms[i] = new MechanismStandard("SEClamp")
ms[i].set("rs", 1, 0)
ms[i].set("dur1", 100, 0)
ms[i].set("amp1", 10, 0)
ms[i].set("dur2", 0, 0)
ms[i].set("amp2", 0, 0)
ms[i].set("dur3", 0, 0)
ms[i].set("amp3", 0, 0)
mt.select("IClamp") i = mt.selected() maction(i)
hoc_ac_ = 0.22
sec.sec move() d1.flip_to(0)
}
{object_pop() doNotify()}
{
ocbox_ = ocbox_.v1
ocbox_.map("PointProcessManager", 1224, 564, 279.9, 594)
}
objref ocbox_
//End PointProcessManager

{
xpanel("RunControl", 0)
v_init = -65
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 50
xvalue("t","t", 2 )
tstop = 50
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 0.05
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(733,582)
}
{
save_window_ = new Graph(0)
save_window_.size(0,50,-80,40)
scene_vector_[4] = save_window_
{save_window_.view(0, -80, 50, 120, 1260, 151, 300.6, 200.8)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)
}
objectvar scene_vector_[1]
{doNotify()}
when i run it, it gives a graph of soma voltage vs time.
there are 4 spikes/peaks in this graph.

Now if we increase the Iclamp amplitude to 2nA, there will be 5 spikes.

similarly for 4nA there will be 6 spikes/peaks.

we can say that the number of peaks (spike frequency) increases with the
increase in current amplitude.

I want to plot this graph ie on x-axis it will be Iclamp amplitude and on
Y-axis spike frequency corresponding to that current.

for this, first we have to find the spike frequency for a particular current.
for this i looked the following post at the Neuron forum
http://www.neuron.yale.edu/phpBB/viewto ... ?f=8&t=468

and i used the following code:

Code: Select all

objref rec, nil
soma rec =  NetCon(&v(.5), nil)
when i put these line in the .ses file. (with in the Runcontrol braces)
it gives error.

Looking forward for your reply

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

Re: spike frequency vs current graph

Post by ted »

Don't do anything to the ses file. Until you know more about NEURON and hoc, leave ses files completely alone, exactly as NEURON generated them.

Create a plain text file (don't use Word or Notepad or any MicroSoft program--use a plain text editor, like NoteTab or one of the other editors recommended in the NEURON Forum for whatever OS you're using). In that plain text file put the following statements:

Code: Select all

load_file("nrngui.hoc")
load_file("foo.ses")
where you replace foo.ses with whatever you're calling your ses file.

Save this to a file called init.hoc. Double click on init.hoc, and NEURON will start and load your ses file. You should now be able to run your simulation just as if you had started NEURON, then used
File / load session
to load the ses file.

Exit NEURON.

Now you can add whatever statements you like to your init.hoc file,
_after_ the
load_file("foo.ses")
statement.

For example, you should be able to add these two statements

Code: Select all

objref rec, nil
soma rec =  NetCon(&v(.5), nil)
to the end of your init.hoc file. These statements will create a NetCon that watches soma.v(0.5) for spikes. You should be able to double click on init.hoc, then launch simulations by clicking on the Init & Run button.

But you'll want to record spike times to a Vector. Then, after each simulation, you'll want to use those spike times to calculate firing frequency. And you'll also want to keep track of the current amplitude that you're using to stimulate the cell. In pseudocode,

Code: Select all

create Vectors called ivec and stimvec

for each stimulus amplitude of interest {
  run a simulation
  analyze the recorded spike times to find the firing frequency
  append the stimulus current to ivec, and the firing frequency to fvec
}

create a Graph called g
use Vector plot() to plot fvec vs. ivec
rescale g's axes by using Graph class's exec_menu() to execute the command "View = plot"
Analyzing the recorded spike times to find the firing frequency should be done by a function that takes one argument: the Vector of recorded spike times. It should test to make sure that there are at least two spikes--otherwise you can't calculate frequency. Also, you might want to make sure that the cell has had sufficient time for its firing frequency to stabilize. How long to wait depends on your model. Models that involve hh usually take about 4 or 5 spikes to settle down to a steady firing frequency. You will want to make sure your program runs long enough to record a sufficent number of spikes.
MSN
Posts: 10
Joined: Mon May 23, 2011 5:38 pm

Re: spike frequency vs current graph

Post by MSN »

I just finished the task and here is the code i wrote. Kindly have a look. May the way i did is not efficient but i got the graph.

Code: Select all

load_file("nrngui.hoc")
load_file("constant_source freq.ses")


objref avec, mvec[21], indx[21], numofspikes, num, fvec, g 
avec = new Vector ()
numofspikes = new Vector ()
num = new Vector ()
fvec = new Vector ()
g = new Graph ()


for i=0, 20 {
mvec[i]=new Vector ()
indx[i]=new Vector ()
}


proc loop () {
for j=0, 20 {IClamp[0].amp = j*0.1
avec.record(&v(0.5))
init ()
run ()
mvec[j].copy(avec)
}
}


proc go() {
loop ()


for k=0, 20 {
for m=0, tstop/dt-1 {
if (mvec[k].x[m]*mvec[k].x[m+1] < 0) {    
indx[k].append(m)                                   
}
}
numofspikes.append(int(indx[k].size()))       /* number of spike which are above the x-axis */
}

num.copy(numofspikes)
fvec=num.div(tstop)                          /* frequency in KHz */
fvec.mul(1000)                               /* frequency in Hz */

g.size(0, 2, 0, 200)
fvec.plot(g, .1)
}
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: spike frequency vs current graph

Post by ted »

MSN wrote:I just finished the task and here is the code i wrote. Kindly have a look. May the way i did is not efficient but i got the graph.
Well, it should run, but does it give the right answer?

This is a clever way to detect zero crossings

Code: Select all

for m=0, tstop/dt-1 {
  if (mvec[k].x[m]*mvec[k].x[m+1] < 0) {    
    indx[k].append(m)                                   
  }
}
but
(1) won't it miss a zero crossing if one datum equals 0 (that would be a false negative)
(2) for zero crossings that do not involve a zero valued datum, won't it count both positive- and negative-going crossings (that would be a double count, or a false positive)
?

So false negatives and double counts would make it highly inaccurate.

NetCon only detects positive going threshold crossings, and the NetCon class's record() method will capture spike times to a vector in the course of a run so there's no need to store the detailed time course of membrane potential or go back and look for spikes after the runs are complete.
MSN
Posts: 10
Joined: Mon May 23, 2011 5:38 pm

Re: spike frequency vs current graph

Post by MSN »

you are right. it made following changes to the code:

Code: Select all

if (mvec[k].x[m]*mvec[k].x[m+1] <= 0) {    
indx[k].append(m)                                   
}
so it will not miss if the datum is 0 and

Code: Select all

numofspikes.append(int(indx[k].size()/2))       /* number of spike which are above the x-axis */
so this will cancel the effect of double count.

I also try to use the NetCon class record() function. But i did not get sucess.
i also try it in a simple example. the code is as follows:

Code: Select all

load_file("nrngui.hoc")
load_file("example.hoc")

objref nc, nil, veca, vecb
veca = new Vector()
vecb = new Vector()

proc algo () {
nc = new NetCon(&v(0.5), nil) 
nc.threshold = 0
veca.record(&v(0.5))
nc.record(veca)
init ()
run ()
vecb = nc.get_recordvec()
}

proc go () {
algo ()
}
I was hoping that vecb will be the collection of indexes of zero crossing. but it is just the duplicate the vector veca.

Can you point out what mistake i am making ?

Regards
MSN
Posts: 10
Joined: Mon May 23, 2011 5:38 pm

Re: spike frequency vs current graph

Post by MSN »

I discover my mistake.
wow, it was so easy.
I think, if it is explicitly mentioned in the help (http://www.neuron.yale.edu/neuron/stati ... tml#record) that the resultant indexes will be saved in the vector (which is given in the argument) then it might be easier for the new users to understand.

Any how my new code is

Code: Select all

load_file("nrngui.hoc")
load_file("constant_source freq.ses")

objref nc, nil, indx[201], numofspikes, num, fvec, g 
numofspikes = new Vector ()
num = new Vector ()
fvec = new Vector ()
g = new Graph ()


for i=0, 200 {
indx[i]=new Vector ()
}


proc loop () {
for j=0, 200 {IClamp[0].amp = j*0.01
nc = new NetCon(&v(0.5), nil) 
nc.threshold = 0
nc.record(indx[j])
init ()
run ()
}
}


proc go() {
loop ()

for k=0, 200 {
numofspikes.append(indx[k].size())       /* number of spike which are above the x-axis */
}

num.copy(numofspikes)
fvec=num.div(tstop-50)                        /* frequency in KHz */
fvec.mul(1000)                               /* frequency in Hz */

g.size(0, 2, 0, 100)
fvec.plot(g, .01)
} 
and TED thank you very much for your time.

I did not use finitialize() in my code, but i noticed that it is is mentioned number of times in the help.
If i want to use it in my code where should i use it, i mean in the loop or in the end of the code.

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

Re: spike frequency vs current graph

Post by ted »

Glad to hear that you worked it out.
MSN wrote:I think, if it is explicitly mentioned in the help (http://www.neuron.yale.edu/neuron/stati ... tml#record) that the resultant indexes will be saved in the vector (which is given in the argument) then it might be easier for the new users to understand.
There is always a balance to be struck between conciseness and verbosity. The documentation of the Vector class's record method begins with
record
. . .

SYNTAX
netcon.record(Vector)
. . .

DESCRIPTION
Records the event times at the source the netcon connects to.
. . .
which seems fairly explicit.

With regard to the new version of the program, here are some tips for program organization and efficiency--

Creation of a NetCon, the Vector to which it records, and the record() statement itself, are what we call "instrumentation code" (to distinguish it from "model setup code" which specifies the biological properties that are represented, and "control code" which controls simulation execution (initialization and "runs")). In cases where instrumentation remains the same from run to run, there is no need to execute the instrumentation code more than once. This is such a case. Your program will be more efficient and easier to debug if you take the instrumentation code (the stuff that creates the NetCon, specifies its threshold, and the record() statement) out of your "control code" (the "for" loop that runs a bunch of simulations), and execute it _before_ you get to the for loop.

"But I won't be able to record to each of the 201 indx Vectors!"

Right, and you don't have to either. If all you want is the firing frequency, create just one indx vector (why not call it a self-documenting name like spiketimes?) change your proc loop() to this (again, choosing a self-documenting name for the procedure)

Code: Select all

proc batchrun () { local j
  for j=0,$1 {
    IClamp[0].amp = j*0.01
    run ()
    postprocess()
  }
}
where postprocess() is a procedure that calculates the frequency and appends it to the freq vector (and it is called postprocess because it contains the code that analyzes simulation results after the simulation ends). That would scale much better because a separate indx vector is not required for each run.

The calculation of firing frequency leaves something to be desired. Many cells (and many models) fire just once and then remain silent, or fire a few spikes early on and then remain silent for a long period. It is meaningless to speak of firing frequency unless there is more than one spike. It is better to base frequency calculations on interspike intervals (ISIs). Given a Vector that contains firing times, it is easy to determine whether there are enough spikes (i.e. more than one) to bother going further, and also easy to calculate the ISIs (read about the Vector class's deriv method, in particular vdest.deriv(vsrc) or perhaps better vsrcdest.deriv() which does the calculation "in place" so does not need to eat up extra storage space).

Of course this opens up the question of what to do about cells that fire bursts of spikes . . .
I did not use finitialize() in my code, but i noticed that it is is mentioned number of times in the help.
If i want to use it in my code where should i use it, i mean in the loop or in the end of the code.
run() calls a procedure that takes care of initialization. This is all part of NEURON's standard run system, which controls the sequence of calculations in a simulation. You might want to examine its source code--it's written in hoc and it's quite readable. See
Secrets of NEURON: the hoc library
viewtopic.php?f=28&t=506
MSN
Posts: 10
Joined: Mon May 23, 2011 5:38 pm

Re: spike frequency vs current graph

Post by MSN »

nice suggestion, now the code is looking organized and is efficent . here it is

Code: Select all

load_file("nrngui.hoc")
load_file("constant_source freq.ses")


// processing //
proc batchrun () { local j
  for j=0, 200 {
    IClamp[0].amp = j*0.01
    run ()
    postprocess()
  }
}


// postprocessing //
objref nc, nil, indx, numofspikes, num, fvec
indx = new Vector ()
numofspikes = new Vector ()
num = new Vector ()
fvec = new Vector ()

proc postprocess () {
nc = new NetCon(&v(0.5), nil) 
nc.threshold = 0
nc.record(indx)                              /* record indexes of the positive zero crossings */
numofspikes.append(indx.size())              /* number of spikes for each value of IClamp */
}

proc integrate () {
num.copy(numofspikes)
fvec=num.div(tstop-50)                        /* frequency in KHz */
fvec.mul(1000)                               /* frequency in Hz */
}



// graph display //
objref g
g = new Graph ()

proc gph () {
g.size(0, 2, 0, 100)
fvec.plot(g, .01)
}


// simulation //
proc go() {
init()
batchrun ()
integrate ()
gph ()
}
But i donot understand how to use vsrcdest.deriv() to calculate frequency or ISI
what should be dx
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: spike frequency vs current graph

Post by ted »

Code: Select all

proc postprocess () {
nc = new NetCon(&v(0.5), nil) 
nc.threshold = 0
nc.record(indx)                              /* record indexes of the positive zero crossings */
numofspikes.append(indx.size())              /* number of spikes for each value of IClamp */
}
This will spawn a new NetCon on each run. No need to do that. Create one NetCon immediately after model setup is complete and reuse that same NetCon, and the Vector to which spike times are recorded, for each run.

Also, it doesn't do anything beyond append the number of spikes to indx. Specifically, it doesn't actually complete the entire calculation required to find the firing frequency. You only know how many spikes there were. You don't know when the first one occurred, or when the last one occurred. You only know the total duration of a simulation. That's insufficient to find firing frequency.

Assuming your program is organized like so
model setup code (creates the model cell itself)
instrumentation code (attaches stimulus source, and sets up NetCon recording of spik e times to a Vector)
and that the Vector is called spiketimes
then
objref isivec // read as "vector of ISIs" i.e. of interspike intervals
isivec = spiketimes.c // so operations on isivec don't destroy the original spik e times
isivec.deriv // in-place calculation of ISIs
At this point the elements of isis are the interspike intervals of the spike train generated by the model cell during the simulation. From these you can determine instantaneous firing frequency, and whether the cell's firing rate has stabilized or is still changing.
santana_loren9
Posts: 12
Joined: Thu Oct 19, 2017 6:42 pm

Re: spike frequency vs current graph

Post by santana_loren9 »

First off, I am so sorry for asking so many questions. I do try to do some research on the forum before asking. I unfortunately chose to learn NEURON in the last few weeks of the semester. I promise I'll get better!
With that said, I am having issues figuring out the spike frequency. I am trying to do two things:
1) Get the spike frequency without considering inputs. I expect to get a number out of this.
2) Getting spike frequency considering inputs. I expect to get a graph out of this.

I have a hoc file with the following:

Code: Select all

load_file("nrngui.hoc")

create soma
access soma

soma {
    L = 15
    diam = 15
    Ra = 70
    insert nahh
    	gbar_nahh = 0.01
    insert kv1_gp
    	gbar_kv1_gp = 0
    insert kv2_hh
    	gbar_kv2_hh = 0
    insert kv2_simplehh
    	gbar_kv2_simplehh = 0.00015
    insert kv4hh
    	gbar_kv4hh = 0
    insert kcnq_hh
    	gbar_kcnq_hh = 7e-07
    insert kir2_dop
    	gbar_kir2_dop = 0
    insert sk_dop
	gbar_sk_dop = 5e-06
    insert hcn_siegelbaum
	gbar_hcn_siegelbaum = 7e-05
	ehcn_hcn_siegelbaum = -10
    insert cal_dop
	gbar_cal_dop = 5e-05
    insert cap_dop
	gbar_cap_dop = 2.5e-05
    insert leak_dop
    	gbar_leak_dop = 5e-06
    	e_leak_dop = -45
    insert CaDiffuse3
}

v_init = -68
celsius = 35
Pmax_CaDiffuse3 = 125
TotalBuffer_CaDiffuse3 = 0.0003
ek = -90
cai_init = 5e-05

objectvar stim
stim = new IClamp(.5)
stim.del = 500
stim.dur = 200
stim.amp = 0

cvode_active(0)
proc init() {
	finitialize(v_init)
	soma.cai = cai_init
}

tstop = 10000
run()
In attempt to follow the conversation above, for the first part, I have in a hoc file titled init.hoc:

Code: Select all

load_file("nrngui.hoc")
load_file("build.hoc")

objref rec, nil, ind, spikenum, fvec
ind = new Vector ()
spikenum = new Vector ()
fvec = new Vector ()

soma rec = new NetCon(&v(0.5), nil) 		      

proc sim () {
run () 							
recordspikes()
compfreq()
}

proc recordspikes () {
rec.threshold = 0
rec.record(ind)                              
spikenum.append(ind.size())              
}

proc compfreq () {
fvec=spikenum.div(tstop)                  
fvec.mul(1000)                               
}
I'm not sure how to get a number out of this. When I type sim(), nothing is produced.

When I run init.hoc, this is the message I get:
-e
NEURON -- VERSION 7.5 master (6b4c19f) 2017-09-25
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

loading membrane mechanisms from /Users/lorensantana/Desktop/snc_compensation/x86_64/.libs/libnrnmech.so
Additional mechanisms from files
Ca.mod CaDiffuse2.mod CaDiffuse3.mod CaHandler2.mod bk_dop.mod cal_dop.mod cap_dop.mod hcn_siegelbaum.mod kcnq_hh.mod kir2_dop.mod kv1_gp.mod kv2_hh.mod kv2_simplehh.mod kv4hh.mod leak_dop.mod nahh.mod sk_dop.mod vramp.mod
1
1
0
1

I wasn't sure if loading a hoc rather than ses file was a problem so I ran build.hoc, checked that the voltage graph looked as I expected and saved the session as build.ses within the same folder as the hoc and mod files. Now, when I run init.hoc, I get:
-e
NEURON -- VERSION 7.5 master (6b4c19f) 2017-09-25
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

loading membrane mechanisms from /Users/lorensantana/Desktop/snc_compensation/x86_64/.libs/libnrnmech.so
Additional mechanisms from files
Ca.mod CaDiffuse2.mod CaDiffuse3.mod CaHandler2.mod bk_dop.mod cal_dop.mod cap_dop.mod hcn_siegelbaum.mod kcnq_hh.mod kir2_dop.mod kv1_gp.mod kv2_hh.mod kv2_simplehh.mod kv4hh.mod leak_dop.mod nahh.mod sk_dop.mod vramp.mod
1
1
/Applications/NEURON-7.5/nrn/x86_64/bin/nrniv: syntax error
in init.hoc near line 9
soma rec = new NetCon(&v(0.5), nil)
^
oc>
My session file looks like:
-e
NEURON -- VERSION 7.5 master (6b4c19f) 2017-09-25
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

loading membrane mechanisms from /Users/lorensantana/Desktop/snc_compensation/x86_64/.libs/libnrnmech.so
Additional mechanisms from files
Ca.mod CaDiffuse2.mod CaDiffuse3.mod CaHandler2.mod bk_dop.mod cal_dop.mod cap_dop.mod hcn_siegelbaum.mod kcnq_hh.mod kir2_dop.mod kv1_gp.mod kv2_hh.mod kv2_simplehh.mod kv4hh.mod leak_dop.mod nahh.mod sk_dop.mod vramp.mod
oc>
For the second part, I figured the same code that MSN (user above) used should work for me as long as I use "stim" instead of "IClamp[0]". The only difference between this code and the one that MSN posted last is that 1) I'm loading a hoc file rather than a session file (is this bad?), 2) I moved up the netcon object to avoid having it be reproduce each time through the loop, 3) I eliminated the num variable since it was a replica of numofspikes, and 4) I changed IClamp[0] to stim.

This is the code I have:

Code: Select all

load_file("nrngui.hoc")
load_file("build.hoc")

// postprocessing //
objref nc, nil, indx, numofspikes, fvec
indx = new Vector ()
numofspikes = new Vector ()
fvec = new Vector ()
soma nc = new NetCon(&v(0.5), nil) 

// processing //
proc batchrun () { local j
  for j=0, 200 {
    stim.amp = j*0.01
    run ()
    postprocess()
  }
}

proc postprocess () {
nc.threshold = 0
nc.record(indx)                              /* record indexes of the positive zero crossings */
numofspikes.append(indx.size())              /* number of spikes for each value of IClamp */
}

proc integrate () {
fvec=numofspikes.div(tstop-50)                        /* frequency in KHz */
fvec.mul(1000)                               /* frequency in Hz */
}



// graph display //
objref g
g = new Graph ()

proc gph () {
g.size(0, 2, 0, 100)
fvec.plot(g, .01)
}


// simulation //
proc go() {
init()
batchrun ()
integrate ()
gph ()
}
When I run this code, I get:
-e
NEURON -- VERSION 7.5 master (6b4c19f) 2017-09-25
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

loading membrane mechanisms from /Users/lorensantana/Desktop/snc_compensation/x86_64/.libs/libnrnmech.so
Additional mechanisms from files
Ca.mod CaDiffuse2.mod CaDiffuse3.mod CaHandler2.mod bk_dop.mod cal_dop.mod cap_dop.mod hcn_siegelbaum.mod kcnq_hh.mod kir2_dop.mod kv1_gp.mod kv2_hh.mod kv2_simplehh.mod kv4hh.mod leak_dop.mod nahh.mod sk_dop.mod vramp.mod
1
1
0
1
oc>
In the terminal I type go() and I get a graph but I do not believe it is correct. For example, it shows that when the stimulus amplitude is 0, the frequency is 0, which I know that's not true. The model is of a pacemaker so it spikes even without external input.

Does anyone know what I'm doing incorrectly in any of these?
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: spike frequency vs current graph

Post by ted »

Before discussing determining the f-I relationship, let's talk about your model. Does it matter to you whether the model cell is initialized to a steady resting state? If yes, how do you know that

Code: Select all

proc init() {
	finitialize(v_init)
	soma.cai = cai_init
}
initializes the model cell to steady state? Besides, this initialization breaks vector recording of soma.cai, and it will not work with adaptive integration.
santana_loren9
Posts: 12
Joined: Thu Oct 19, 2017 6:42 pm

Re: spike frequency vs current graph

Post by santana_loren9 »

Hi Ted,

Thanks for the question. To be honest, I don't know the answer as to what the authors wanted. I am trying to replicate a paper and this was the model that they put up. The f-I curve is something I'm checking for something else, but I intended to use the base model as given.

I normally use Matlab when simulating and I do start at steady state.
The mod files contain infinity curves obtained with the initial voltage and within the mod files there are inital "blocks" that set the steady state gating variables. Do you not think that having finitialize(v_init) does set the cell in steady state?

If this initialization does break the vector recording, is there a better way to go about it?
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: spike frequency vs current graph

Post by ted »

finitialize(v_init) initializes all voltage-dependent gating variables to their steady state values at v = v_init. It also executes whatever initialization code exists in ion transport and accumulation mechanisms. But it does nothing--and it can't do anything--that ensures ionic concentrations have reached their steady state values in all compartments. Why is this important? If ionic concentrations are not at their steady state at t = 0, concentrations and reversal potentials will begin changing from the start of the simulation. By definition, the model is not at its steady state.

Ensuring that concentrations are at their steady states requires that this condition be met:

Code: Select all

for each ionic species x that has an ion accumulation mechanism
  for each compartment that has that ion accumulation mechanism
    ix equals 0
This condition means that dxo/dt = dxi/dt = 0 at t = 0 (i.e. rates of change of extracellular and intracellular concentrations of x are zero at the start of the simulation). It can be satisfied by custom code that is either located in proc init() or executed by an FInitializeHandler. Does such code accompany the model you're working with? Assigning a particular value to xi, e.g. an assignment statement like
soma.cai = cai_init
is not the answer.
Post Reply