How to put a poisson stimulation into a fixed point of axon

[thisman]

I read 'the NEURON book' and the following link,
http://www.neuron.yale.edu/phpBB/viewto ... ?f=8&t=416
but I still have no idea of how to put a poisson stimulation into an axon continuously,
let say the poisson stimulation is 5Hz for 5 minutes,
any suggestions?

[ted]

The first step is to create a model of an axon and attach a synaptic mechanism to it. Decide whether a single synaptic activation produces a conductance change that is best approximated by an ExpSyn or an Exp2Syn, and also what the synapse's reversal potential and unitary synaptic conductance should be.

To activate the synapse with a spike train that has "poisson" statistics (i.e. interspike intervals are governed by the negative exponential distribution), create a NetStim object and set its "random" parameter to 1 and its "interval" parameter to 1000/f where f is the mean frequency of activation in Hz. If you want 5 Hz, "interval" should be 200. You'll probably want to reduce "start" to 10 ms or less*. And you'll definitely want to increase "number" to 1e9.

You'll need a NetCon to attach the NetStim to the synaptic mechanism. The NetCon's weight should be set to the peak value of the synaptic conductance elicited by a single synaptic activation.

*--Why wait 50 ms for anything to happen?

Since synaptic activation will occur at random times, there is no way to guarantee that the first event will occur at a particular time START, or that the last event will occur at a particular time START+DURATION. However, you can ensure that the synapse is _not_ activated before START or after START+DURATION.

The following code assumes that (1) the synaptic mechanism is an ExpSyn or Exp2Syn, (2) it has been attached to the axon, and (3) it is called syn.

Code: Select all

WEIGHT = some_value // that you decide. Units are microsiemens.
FREQ = 5 // Hz
START = 10 // earliest possible time of synaptic activation
DURATION = 5000 // ms--length of window in which events may be delivered
  // to the synaptic mechanism

objref ns
ns = new NetStim()
ns.random = 1
ns.start = 0
ns.number = 1e9
ns.interval = 1000/FREQ

objref nc
nc = new NetCon(ns, syn)
nc.delay = 0
nc.weight = WEIGHT

objref fih
fih = new FInitializeHandler("initstim()")

proc initstim() {
  nc.active = 0 // nc is inactive, won't deliver any events
  cvode.event(START, "setstim()")
}

proc setstim() {
//  print "time is ", t
  if (nc.active==0) {
    nc.active = 1 // make netcon active so it will deliver events
    // after netcon has been active for DURATION ms, make it inactive
    cvode.event(t + DURATION, "setstim()")
  } else {
    nc.active = 0 // allow no more events to pass
  }
  // we have changed a parameter abruptly
  // so we really should re-initialize cvode
  if (cvode.active()) {
    cvode.re_init()
  } else {
    fcurrent()
  }
}

[thisman]

Thanks Ted, I did as your suggestion, and below are my .hoc and .mod codes.

In the .hoc part I tried to construct an axon with HH property and Ih current(which described by the .mod part), then I connected this axon with 'syn' that defined by 'ExpSyn'. But there are still some problems when I run it, actually I don't know the meaning of the last two paragraphs of your code(the 'initstim()' and 'setstim()' paragraph).

Here is what I want to do exactly, I want to put 5 minutes poisson stimulation with different frequencies into the beginning of the axon. For each stimuli, it will generate an action potential and transport from the beginning to the end of the axon. If we record the times when the spike appears at the beginning and the middle of the axon separately, then there is a time difference-'delay'(i.e. how long the spike need to transport from the beginning to the middle in the axon). This delay is what I want to get.

Thanks again!

Code: Select all

load_file("nrngui.hoc")

create axon
access axon

axon {
  nseg = 200
  diam = 10
  L = 20000
  Ra = 100
  insert hh
  insert Ihp
  eh = -20
}

objref syn
axon syn = new ExpSyn(0.0025)

WEIGHT = 1 // that you decide. Units are microsiemens.
FREQ = 5 // Hz
START = 10 // earliest possible time of synaptic activation
DURATION = 5000 // ms--length of window in which events may be delivered
                // to the synaptic mechanism

objref ns
ns = new NetStim()
ns.noise = 1
ns.start = 0
ns.number = 1e9
ns.interval = 1000/FREQ

objref nc
nc = new NetCon(ns, syn)
nc.delay = 0
nc.weight = WEIGHT

objref fih
fih = new FInitializeHandler("initstim()")

proc initstim() {
  nc.active = 0 // nc is inactive, won't deliver any events
  cvode.event(START, "setstim()")
}

proc setstim() {
//  print "time is ", t
  if (nc.active==0) {
    nc.active = 1 // make netcon active so it will deliver events
    // after netcon has been active for DURATION ms, make it inactive
    cvode.event(t + DURATION, "setstim()")
  } else {
    nc.active = 0 // allow no more events to pass
  }
  // we have changed a parameter abruptly
  // so we really should re-initialize cvode
  if (cvode.active()) {
    cvode.re_init()
  } else {
    fcurrent()
  }
}

Code: Select all

TITLE Hyperpolarization activated inward current for axon

UNITS {
  (mV) = (millivolt)
  (mA) = (milliamp)
}

NEURON {
    SUFFIX Ihp
    NONSPECIFIC_CURRENT ih
    RANGE gmax
}

PARAMETER {
    gmax  = 0.0001 (mho/cm2)
}

ASSIGNED { 
    v (mV)
    eh (mV)
    ih (mA/cm2)
    taur (/ms)
    rinf (/ms)
}

STATE {
    r
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    ih  = gmax*r*(v-eh)
}

INITIAL {
    settables(v)
    r = rinf
    
}

DERIVATIVE states {  
    settables(v)      
    r' = (rinf-r)/taur

}

UNITSOFF

PROCEDURE settables(v (mV)) {
    TABLE taur, rinf 
    FROM -100 TO 100 WITH 200

    taur = 500+400/(1+exp(-(v+70)/5))
    rinf = 1/(1+exp((v+65)/4))
}

UNITSON

[ted]

there are still some problems when I run it

If you don't describe them, I will never know what they are.

actually I don't know the meaning of the last two paragraphs of your code(the 'initstim()' and 'setstim()' paragraph).

You want to start stimulating at a particular time START, and you want the stimulation to recur for a specified amount of time DURATION. This can be done by deactivating the NetCon until t == START, when the NetCon is activated, then deactivating the NetCon when t == START+DURATION.
procs initstim() and setstim() take care of this. To understand how, read the Programmer's Reference documentation of the FInitializeHandler class and the CVode class's event() method.

In the meantime, here's a brief explanation:
proc initstim() is executed during simulation initiation. It deactivates the NetCon at t == 0 and executes cvode.event(start, "setstim()").
cvode.event(start, "setstim()") makes proc setstim() be executed when t == START.
When setstim() is executed at t == START, the NetCon is not active, so setstim() activates it and executes cvode.event(START + DURATION, "setstim()").
cvode.event(START + DURATION, "setstim()") makes proc setstim() be executed when t == START + DURATION.
When setstim() is executed at t == START + DURATION, the NetCon is active, so setstim() inactivates it.

I want to put 5 minutes poisson stimulation with different frequencies into the beginning of the axon. For each stimuli, it will generate an action potential

Not if the stimulus arrives during the axon's refractory period. Remember that for interstimulus intervals governed by the negative exponential distribution, the most likely interstimulus interval is 0.

how long the spike need to transport from the beginning to the middle in the axon

Then you will want to use an odd number for nseg--see
Why should I use an odd value for nseg?
in NEURON's FAQ list
http://www.neuron.yale.edu/neuron/faq/general-questions

[thisman]

If you don't describe them, I will never know what they are.

After I compile the .mod file as below

Code: Select all

TITLE Hyperpolarization activated inward current for axon

UNITS {
  (mV) = (millivolt)
  (mA) = (milliamp)
}

NEURON {
    SUFFIX Ihp
    NONSPECIFIC_CURRENT ih
    RANGE gmax
}

PARAMETER {
    gmax  = 0.0001 (mho/cm2)
}

ASSIGNED { 
    v (mV)
    eh (mV)
    ih (mA/cm2)
    taur (/ms)
    rinf (/ms)
}

STATE {
    r
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    ih  = gmax*r*(v-eh)
}

INITIAL {
    settables(v)
    r = rinf
    
}

DERIVATIVE states {  
    settables(v)      
    r' = (rinf-r)/taur

}

UNITSOFF

PROCEDURE settables(v (mV)) {
    TABLE taur, rinf 
    FROM -100 TO 100 WITH 200

    taur = 500+400/(1+exp(-(v+70)/5))
    rinf = 1/(1+exp((v+65)/4))
}

UNITSON

I run my .hoc file as below

Code: Select all

load_file("nrngui.hoc")

create axon
access axon

axon {
  nseg = 201
  diam = 10
  L = 20000
  Ra = 100
  insert hh
  insert Ihp
  eh = -20
}

objref syn
axon syn = new ExpSyn(0.0025)

WEIGHT = 0.5 // that you decide. Units are microsiemens.
FREQ = 5 // Hz
START = 10 // earliest possible time of synaptic activation
DURATION = 5000 // ms--length of window in which events may be delivered
                // to the synaptic mechanism

objref ns
ns = new NetStim()
ns.noise = 1
ns.start = 0
ns.number = 1e9
ns.interval = 1000/FREQ

objref nc
nc = new NetCon(ns, syn)
nc.delay = 0
nc.weight = WEIGHT

objref fih
fih = new FInitializeHandler("initstim()")

proc initstim() {
  nc.active = 0 // nc is inactive, won't deliver any events
  cvode.event(START, "setstim()")
}

proc setstim() {
  //print "time is ", t
  if (nc.active==0) {
    nc.active = 1 // make netcon active so it will deliver events
                  // after netcon has been active for DURATION ms, make it inactive
    cvode.event(t + DURATION, "setstim()")
  } else {
    nc.active = 0 // allow no more events to pass
  }
  // we have changed a parameter abruptly
  // so we really should re-initialize cvode
  if (cvode.active()) {
    cvode.re_init()
  } else {
    fcurrent()
  }
}

then I have this error

oc>nrniv: Cannot assign to left hand side
near line 72
{initstim()}
^
initstim( )
finitialize(-65 )
init( )
stdinit( )
and others

So what's mean of "cannot assign to left hand side"? Could you please fix this .hoc file?
What's mean of the function "fih = new FInitializeHandler("initstim()")" in your code?
If this code works, how can I stop the simulation? and how can I control the amplitude of the stimulation?
Did I miss some function in my .hoc code?

Thank you very much Ted, again!

[ted]

Could you please fix this .hoc file?

If I did, what would you learn? Here's how to figure it out for yourself.

So what's mean of "cannot assign to left hand side"?

Seems pretty clear: it's complaining about an assignment statement, i.e. a statment of the form
a = b
So the problem is in one of the assignment statements in the hoc file.

But which one? This is a run time error message (it occurs during simulation), and hoc often gives rather vague clues about the causes of run time errors. But this one seems rather specific: there's a problem with initstim(). So take a close look at the contents of proc initstim().

A good strategy for localizing errors is to look for statements that use unfamiliar keywords. If you're not sure about the usage of a particular keyword, read about it in the Programmer's Reference. proc initstim() doesn't contain a lot of keywords, so you won't have much reading to do before you discover the problem. Once you find the error, look for other instances of it and fix them too.

"Aha! ted made a mistake!"

Yes, and you can figure out how to fix it. There are several lessons here. First, trust noone's code, not even your own (I know this one particularly well). Second, the only way to stop making mistakes is to stop doing things. Third, with continued practice it is possible to develop skill in debugging, especially when there is no "integrated development environment." Fourth, except for the most trivial cases, program development always proceeds in a cycle of incremental revision and testing. Incremental revison means that, even though one may--and should--have an overall plan in mind, and may even have constructed a detailed outline of what the final program will look like, only the most trivial programs can be implemented at one pass. Anything with reasonable complexity will require multiple iterations.

What's mean of the function "fih = new FInitializeHandler("initstim()")" in your code?

Have you read the Programmer's Reference documentation of the FInitializeHandler class? Have you searched the NEURON Forum for discussions that involve this keyword?

If this code works, how can I stop the simulation?

If I were you, I'd use NEURON's standard run system to control simulation flow ("simulation flow" includes starting and stopping simulations). The easiest way to do this during program development is with the GUI. To learn how to use the GUI, work through this tutorial
Construction and Use of Models: Part 1. Elementary Tools
which is listed on the Documentation page of NEURON's WWW site http://www.neuron.yale.edu/.

After you have everything working the way you want it, you may want to use hoc statements to control program execution. But believe me, you have a lot to discover about this simple model you have implemented, and it will be much easier if you start by using the NEURON Main Menu toolbar to construct a custom user interface with a RunControl panel and a "Voltage axis" graph that shows axon membrane potential at both 0 and 0.5. For example, you need to discover the axon's resting potential so you can specify the proper initial value for membrane potential (otherwise your simulations will contain an initialization transient that may take many tens of milliseconds to decay). You also need to know how long a simulation must run in order to have a reasonable chance of including at least one synaptic event. And you may want to adjust the strength of the synaptic input, or the synaptic time constant.

how can I control the amplitude of the stimulation?

By changing the NetCon's weight; read in the Programmer's Reference about NetStim and ExpSyn.


Comments about the NMODL code:

I'd suggest changing the suffix in the NMODL block to lower case, specifically hp. Why?
1. First, about lower case: Capitalized keywords are generally reserved for the names of object classes. Point processes are objects, so the NMODL code that defines a class of point process uses a capitalized word for the name of the object class. Your h current is a density mechanism, not an object class. NMODL code that defines a density mechanism therefore uses a lower case word for the name of the density mechanism.
2. Next, about whether or not to use an "i" or an "I" (upper case i) as the first letter in the name of a mechanism. Most mechanisms that generate a current, whether they are point processes or density mechanisms, declare a variable that starts with "i" in the NEURON block. This variable can be used by user-written hoc statements to discover the current that is generated by the mechanism.
The complete hoc name of this variable will be of the form
i..._suffix for a density mechanism,
or, in the case of point processes,
ClassName.i... or objrefname.i...
If the suffix starts with "i", or the ClassName or objrefname starts with i, then the user-written hoc code will have variable names with redundant i's, like this
i_ihp
or
IGABA.i
or
igaba.i
If you were writing a paper, you might have a variable called i_hp or GABAi or igaba, but nobody would deliberately choose names like i_ihp, IGABAi, or igabai (or even worse, ih_ih, ik_ik etc.). So why write NMODL code in a way that forces similar hiccups in hoc?

The ASSIGNED block declares
taur (/ms)
rinf (/ms)
which is silly. Won't cause an error in your simulation, but will gag modlunit (have you checked this mod file with modlunit?) and will make experienced NEURON users lift an eyebrow. Time constants in NMODL use (ms), and steady state values of dimensionless state variables are dimensionless (the corresponding units specification is (1)).

Comments about the hoc code:

The statement
eh = -20
does nothing useful because the name on the left hand side is not used anywhere else in the program. You probably meant eh_hp.
By the way, to eliminate this hiccup you might want to change the declaration in the ASSIGNED block of the mod file to
e (mV)
Don't forget to revise the "ih =" statement in the BREAKPOINT block.

[thisman]

Hi Ted, thanks for your patient reply.
After read the Programmer's Reference, I used

boolean = nc.active()

instead of nc.active = 0, and

boolean = nc.active(boolean)

instead of nc.active = 1.
Also I fixed my .mod as your suggestion.

For each stimuli, if it generates an action potential in the beginning of the axon, let say v(0.0025),
then we also can find the same action potential in the middle of the axon, let say v(0.5).
Suppose we fix a threshold for whole stimulation process, then we can record the times when action potentials
arrive at this threshold for all spikes, then I can find the time difference (delay) for each action potential.
For now, the

print "time is", t

in your code just give me one record for the whole stimulation process, which equal to START.
How can I get the times that all spikes arrive at the same threshold in the whole process?
(Well, this is my idea that how can we measure the delay of the same action potential transfers from the beginning to
the middle of the axon, you may have better idea to get the same result.)

[ted]

After read the Programmer's Reference, I used
boolean = nc.active()
instead of nc.active = 0, and
boolean = nc.active(boolean)
instead of nc.active = 1.

Not quite. The syntax examples are concise to the point of unintentional obfuscation. Let me explain.

boolean = netcon.active()
means that, when netcon.active() is called without an argument, it returns a boolean value (a 0 or a 1) that reflects "the current state" of the NetCon.
boolean = netcon.active(boolean)
means that if netcon.active() is called with an argument, two things happen:
1. the value of that argument is interpreted as a boolean value that specifies whether the NetCon is to be turned on or off (depending on whether the argument is 1 or 0)
and
2. a boolean value is returned that reports the previous state of the NetCon (1 if on and 0 if off).

Example: start NEURON, and at the oc prompt enter the following commands

Code: Select all

oc>create soma
oc>access soma
oc>objref syn
oc>syn = new ExpSyn(0.5)
oc>objref ns, nc
oc>ns = new NetStim()
oc>nc = new NetCon(ns, syn)
oc>nc.active()
        1            << means the NetCon will allow events to pass
oc>nc.active(0)  << this turns the NetCon off
        1            << and reports the previous state of the NetCon
oc>nc.active()  << reports the current state of the NetCon
        0              << sure enough, it's off
oc>nc.active(1)  << turns the NetCon on
        0             << and reports the previous state of the NetCon
oc>

In the original code that I suggested, I was treating nc.active as if it were a constant, to which I could assign a value. Every instance where I suggested
nc.active = 0
or
nc.active = 1
should be replaced by
nc.active(0)
or
nc.active(1)
respectively.

For each stimuli, if it generates an action potential in the beginning of the axon, let say v(0.0025)

you mean soma(0.0025). soma(0.0025) is the location in soma that is 0.25% of the distance from soma's 0 end to soma's 1 end. v(0.00250 is the membrane potential at the location 0.0025 in the currently accessed section, i.e. at soma(0.0025) if the soma is the currently accessed section.

Suppose we fix a threshold for whole stimulation process, then we can record the times when action potentials
arrive at this threshold for all spikes, then I can find the time difference (delay) for each action potential.

OK, you want to determine the conduction latency between two points.

For now, the

print "time is", t

in your code just give me one record for the whole stimulation process, which equal to START.

That was just a diagnostic statement, to report the time at which the NetCon would allow events to pass from NetStim to synaptic mechanism. It does not mean that the axon is stimulated at that time. Remember that the NetStim's output occurs at random times. Every so often you may get a simulation in which NO stimulus occurs at all during the 5 seconds that the NetCon is active.

How can I get the times that all spikes arrive at the same threshold in the whole process?
(Well, this is my idea that how can we measure the delay of the same action potential transfers from the beginning to the middle of the axon, you may have better idea to get the same result.)

Easily done. But the first question to ask is whether you really want to measure at those two locations. There are confounding factors that you may wish to avoid:

1. You probably don't want to measure near the spike initiation zone. Spike amplitude and time course (and conduction velocity) near the point of spike initiation are not the same as they are after the spike has had the opportunity to propagate a few tens or hundreds of microns.
2. You also probably don't want to measure near the distal end of the axon, where sealed end effects will change spike amplitude and time course (and conduction velocity).

A reasonable rule of thumb would be to measure at points that are at 1/3 and 2/3 of the length of the axon. If you make nseg a multiple of 5, there will be nodes at 0.3 and 0.7, which are close enough to 1/3 and 2/3. These would be the sites at which to observe the spike waveform.

For greatest accuracy, you will want to place your "threshold" for measuring spike latency at the point on the waveform where the slope is steepest. As a rule of thumb, this will be at the voltage that lies halfway between resting poential vrest and the depolarized peak vmax of the spike. You'll have to determine vrest and vmax yourself.

You can use a NetCon for each spike detector; just set its threshold to (vrest + vmax)/2. Use the NetCon class's record() method to record the spike times detected by the two NetCons to a pair of Vectors (read about this in the Programmer's Reference).

Code: Select all

objref nc1, nc2, tvec1, tvec2

// set up detection and spike recording at location 1
axon nc1 = new NetCon(&v(0.3), calclatency())
nc1.threshold(whatever_you_choose)
tvec1 = new Vector()
nc1.record(tvec1)

// set up detection and spike recording at location 2
axon nc2 = new NetCon(&v(0.5), calclatency())
nc2.threshold(whatever_you_choose)
etc.

You can automate the process of recording and calculating spike times by creating a procedure like this:

Code: Select all

proc doall() { local i
  run()
  // use the recorded spike pairs to calculate the conduction latencies.
  // there will be as many pairs as there are elements in tvec2.
  if (tvec2.size > 0) {
    print "conduction latencies:"
    for i=0,tvec2.size-1 print tvec2.x[i] - tvec1.x[i]
  } else {
    print "sorry, no conduction latencies this time!"
  }
  print " "
}

You can then type
doit()
at the oc> prompt, and presto! a simulation executes, and you get your results.

The latencies will only have first order accuracy if you are simulating with fixed time steps. If you want second order accurate results, insert the following statements before the definition of proc doall()--

Code: Select all

cvode.active(1) // use adaptive integration
cvode.condition_order(2) // use linear interpolation to determine time of threshold crossing

Be sure to read about the CVode class, cvode.active, and cvode.condition_order in the Programmer's Reference.

[thisman]

Hi Ted, thanks again!
I read the information about "nc.record(Vector)", but I don't know where should I put it(the following code) in the original code

Code: Select all

objref nc1, nc2, tvec1, tvec2, calclatency

axon nc1 = new NetCon(&v(0.3), calclatency)
nc1.threshold = -10
tvec1 = new Vector()
nc1.record(tvec1)

axon nc2 = new NetCon(&v(0.7), calclatency)
nc2.threshold = -10
tvec2 = new Vector()
nc2.record(tvec2)

Should I put it before/after/into the "proc setstim()" function in the original code? (I tried all of them, seems not work.)
What's mean of "calclatency"? I defined it in the "objref" sentence and removed the "()" in the "new NetCon(&v(0.3),calclatency)".
I put the

Code: Select all

proc doall() { local i
  run()
  if (tvec2.size > 0) {
    print "conduction latencies:"
    for i=0,tvec2.size-1 print tvec2.x[i] - tvec1.x[i]
  } else {
    print "sorry, no conduction latencies this time!"
  }
  print " "
}

in the end of my code, but I still did not get the print out result. I think this may be caused by the wrong place of the detection code.
Since I use GUI to control my simulation time, do I need "run()" here?