sinusoid current

[kazmer]

Hello Ted,

I started to simulate subtreshold oscillation of guine-pig cortical neuron. I would like to use sinusoid IClamp current on soma.
equation: Iin=A*(sinB*t^2)+DC
Could I write it in HOV or NMODL?

Thanks,
Kazmer

[kelvin]

Well, I'm not Ted but I figured I could answer this on his behalf. He will probably come along and tell you a more elegant way to do this.

1) Start a new text file in your editor of choice. Call it SinClamp.mod (a suggestion)
2) Copy the code below, paste it in to your file. Then save it and compile it with mknrndll (depending on your platform).
3) I stuck the SinClamp.mod file into a directory where I had an functioning model already, so I started up that model. Then from NEURON Main Menu gui, select Tools>Point Processes>Managers>Point Manager.
4) On the PointProcessManager window that pops up click on SelectPointProcess and down at the bottom you should see SinClamp.
5) I set del = 100, dur = 1000, pkamp = .1, freq = 2 and left phase and bias equal to zero. With this amplitude of current I was still subthreshold in my cell.
6) Get a voltage graph up and ready then hit Init&Run

Hope it works for you. I didn't do too much testing.

Cheers,
Kelvin

Code: Select all

COMMENT
Since this is an electrode current, positive values of i depolarize the cell
and in the presence of the extracellular mechanism there will be a change
in vext since i is not a transmembrane current but a current injected
directly to the inside of the cell.
ENDCOMMENT

NEURON {
        POINT_PROCESS SinClamp
        RANGE del, dur, pkamp, freq, phase, bias
        ELECTRODE_CURRENT i
}

UNITS {
        (nA) = (nanoamp)
             }

PARAMETER {
        del=5	(ms)
        dur=200	(ms)
        pkamp=1 (nA)
        freq=1  (Hz)
        phase=0 (rad)
        bias=0  (nA)
        PI=3.14159265358979323846
}

ASSIGNED {
        i (nA)
}

BREAKPOINT {
       if (t < del) {
	   i=0	
	}else{  
            if (t < del+dur) {
	        i = pkamp*sin(2*PI*freq*(t-del)/1000+phase)+bias
	   }else{  
	        i = 0
}}}

[kazmer]

Thanks, it is fine.
Kazmer

[ted]

Good idea, Kelvin. I have only one suggestion: use at_time() to ensure that the start and
end of the sinusoid are not ignored in simulations that use adaptive integration (cvode).
With this change, the BREAKPOINT block becomes

Code: Select all

BREAKPOINT {
       at_time(del)
       at_time(del + dur)

       if (t < del) {
	   i=0	
	}else{  
            if (t < del+dur) {
	        i = pkamp*sin(2*PI*freq*(t-del)/1000+phase)+bias
	   }else{  
	        i = 0
}}}

[crutchley]

As a point of interest to anyone who may be trying to do the same thing I am, I modified Kelvin's code (thanks, Kelvin and Ted!) to create a swept sine wave to find resonant frequency, in this case from 0-15 Hz linearly over 30 seconds:

Code: Select all

COMMENT 
Since this is an electrode current, positive values of i depolarize the cell 
and in the presence of the extracellular mechanism there will be a change 
in vext since i is not a transmembrane current but a current injected 
directly to the inside of the cell.

Linearly increasing swept sinwave (ZAP function) current injection.
ENDCOMMENT 

NEURON {
        POINT_PROCESS SinSweepClamp 
        RANGE del, dur, pkamp, lofreq, hifreq, phase, bias 
        ELECTRODE_CURRENT i 
} 

UNITS { 
        (nA) = (nanoamp) 
             } 

PARAMETER { 
        del		= 100	(ms) 
        dur    	= 3e4	(ms) 
        pkamp  	= 0.1	(nA) 
        lofreq 	=   0	(Hz)
		hifreq 	=  15	(Hz) 
        phase	=   0	(rad) 
        bias	=   0	(nA) 
        PI		=   3.14159265358979323846 
} 

ASSIGNED { 
        i (nA) 
} 

BREAKPOINT { 
       at_time(del) 
       at_time(del + dur) 

       if (t < del) { 
      i=0    
   }else{  
            if (t < del+dur) { 
           i = pkamp*sin(((2*PI*hifreq - 2*PI*lofreq)/(2*(dur/1000)))*((t-del)/1000)^2+(2*PI*lofreq)*((t-del)/1000)+phase)+bias 
      }else{  
           i = 0 
}}}

I'm very new at this, so be kind, and feel free to point out inefficient code or outright mistakes :-)

[ted]

This implementation of a zap current clamp uses events for compatibility with adaptive integration, guards against nonsense (negative) parameter values, and reports instantaneous frequency. If f0 = f1 it delivers a constant frequency.

Code: Select all

COMMENT
izap.mod
Delivers an oscillating current that starts at t = del >= 0.
The frequency of the oscillation increases linearly with time
from f0 at t == del to f1 at t == del + dur, 
where both del and dur are > 0.
Uses event delivery system to ensure compatibility with adaptive integration.
12/4/2008 NTC
ENDCOMMENT

NEURON {
  POINT_PROCESS Izap
  RANGE del, dur, f0, f1, amp, f, i
  ELECTRODE_CURRENT i
}

UNITS {
  (nA) = (nanoamp)
  PI = (pi) (1)
}

PARAMETER {
  del (ms)
  dur (ms)
  f0 (1/s)  : frequency is in Hz
  f1 (1/s)
  amp (nA)
}

ASSIGNED {
  f (1/s)
  i (nA)
  on (1)
}

INITIAL {
  f = 0
  i = 0
  on = 0

  if (del<0) { del=0 }
  if (dur<0) { dur=0 }
  if (f0<=0) { f0=0 (1/s) }
  if (f1<=0) { f1=0 (1/s) }

  : do nothing if dur == 0
  if (dur>0) {
    net_send(del, 1)  : to turn it on and start frequency ramp
    net_send(del+dur, 1)  : to stop frequency ramp, freezing frequency at f1
  }
}

COMMENT
The angular velocity in radians/sec is w = 2*PI*f, 
where f is the instantaneous frequency in Hz.

Assume for the moment that the frequency ramp starts at t = 0.
f = f0 + (f1 - f0)*t/dur

Then the angular displacement is
theta = 2*PI * ( f0*t + (f1 - f0)*(t^2)/(2*dur) ) 
      = 2*PI * t * (f0 + (f1 - f0)*t/(2*dur))
But the ramp starts at t = del, so just substitute t-del for every occurrence of t
in the formula for theta.
ENDCOMMENT

BREAKPOINT {
  if (on==0) {
    f = 0
    i = 0
  } else {
    f = f0 + (f1 - f0)*(t-del)/dur
    i = amp * sin( 2*PI * (t-del) * (f0 + (f1 - f0)*(t-del)/(2*dur)) * (0.001) )
  }
}

NET_RECEIVE (w) {
  : respond only to self-events with flag > 0
  if (flag == 1) {
    if (on==0) {
      on = 1  : turn it on
    } else {
      on = 0  : turn it off
    }
  }
}

[ted]

I have updated izap.mod so that it now uses the event delivery system to control start and stop time of the sinusoid. The revised file, plus a hoc file and a ses file that demonstrate it, is available here
http://www.neuron.yale.edu/ftp/ted/neuron/izap.zip

Usage:
1. Extract the files from izap.zip
2. Compile the mod file (under MSWin or OS X use mknrndll; for UNIX/Linux use nrnivmodl)
3. Use NEURON to execute initizap.hoc

[shailesh]

Hi,

I was using the Izap.mod file on my model and was getting a few unexpected negative t values when plotting Izap[0].i when CVODE was active. It didn't give any problems with CVODE turned off.

To track down the problem, I reduced my code to its very basic (without making any changes to Izap.mod):

Code: Select all

load_file("nrngui.hoc")
create cell
access cell

The following time sequence was generated:

-0.0125
-71.6102
5
5.0559
...
...

It shows a few initial negative values and then proceeds normally.

On updating a couple of parameters:

Code: Select all

cell.L = 200
cell.diam = 6

The time sequence became:

-0.0125
1.71613
-17156.3
5
5.00866
...
...

It shows a larger negative value initially and then proceeds normally.

I repeated the above a few times and found that this error was observed only on plotting on a current-axis graph and the same isn't seen on a voltage-axis graph. Also the same happens only for CVODE active.

Just to be sure, I checked it on the sample hoc file that was zipped alonged with Izap.mod. I added a current axis graph to plot Izap[0].i in addition to the existing voltage-axis graph for the same. There too the same thing was observed on the current axis graph, but only when pas and hh both were excluded, with CVODE active. The following time sequence was generated:

-0.0125
1.71613
-17151.3
10
10.1316
...
...

Note sure whether it means anything but here as well the exact same initial sequence is generated (-0.0125, 1.71613, -17151.3, ... then normal)

So was wondering what is the cause of the above observation and how does a voltage axis graph differ from a current axis graph in this respect. I assume its something to do with the way CVODE operates, but cant figure out why it would differ for the two graphs, and also how to avoid it.

Thanks,
Shailesh Appukuttan

[ted]

Very strange. I don't see this at all with adaptive integration. With fixed time step integration, the first point plotted in current axis graphs is at -dt/2 by design (because ionic currents are second order correct at t-dt/2 when the Crank-Nicholson method is used--see the Programmer's Reference entry on secondorder, and read about the "staggered time step" algorithm http://www.neuron.yale.edu/neuron/stati ... /nc3p3.htm. However, when adaptive integration is used, all variables are computed and plotted with identical accuracy (3rd order or better) at the same time.

So let me ask what version of NEURON you are using (what does hoc print to its xterm when NEURON starts) and could you please zip up and email to me
ted dot carnevale at yale dot edu
the code you used to generate your observations, so that I can reproduce them for myself.

[ted]

Your code does indeed work as you described. This puzzled the heck out of me, so I asked Michael Hines about this, and here's his answer:

The problem is the direct use of cvode.active(1) instead of the wrapper cvode_active(1) which would, among a number of things, set the variable using_cvode_. So when the NEURONMainMenu/Graph/Currentaxis is used, instead of plotting at times t-dt/2, things would be plotted at t. (note that on the second SingleStep, dt goes to 34323).

This would also be important for NEURONMainMenu/Graph/Stateaxis which plots things at time t+dt/2 for the fixed step method.

Which is news to me. So I cd'd to /usr/local/nrn/share/nrn/lib/hoc (c:\nrn-x.x\lib\hoc under MSWin), then looked for occurrences of using_cvode_

Code: Select all

[ted@loki hoc]$ grep using_cvode_ *hoc
stdrun.hoc:using_cvode_ = 0 // now controlled in NumericalMethodPanel
stdrun.hoc:     if (using_cvode_) return
stdrun.hoc:     if (using_cvode_) {
stdrun.hoc:     if (using_cvode_ && stoprun == 0) { // handle the "tstop" event
stdrun.hoc:     if (using_cvode_) {
stdrun.hoc:     if (using_cvode_) {
varmeth1.hoc:using_cvode_ = 0
varmeth1.hoc:external using_cvode_
varmeth1.hoc:   using_cvode_ = isvar_

Examining stdrun.hoc I find code that indicates that using_cvode_ is part of simulation flow control in NEURON's run time system, and affects the plotting of intermediate results.

Bottom line is that it is best to turn adaptive integration on or off either by using the VariableStepControl GUI tool, or by calling cvode_active() with an argument of 1 or 0, respectively.

[javiZ]

Dear all,
Thanks for the previous discussion, reading this thread helped me to clarify several things.
I've been using the Izap.mod ( taken from http://www.neuron.yale.edu/ftp/ted/neuron/izap.zip) in a simple cell model and it works, however I've noticed that I get some issues with the duration and the frequency range covered.

When I set del = 50; dur = 50
then the current injection starts at 50 ms and lasts until 150 ms (real duration is always del+dur)

On the other hand, if I set
f0 = 0; f1 = 100
then the maximal frequency reached is 199.975 and there is always a byas in the highest frequency.

So far, I couldn't find the problem, I'd appreciate if somebody has a clue.
I'm copying the code I'm using to test it so you can test it as well and see the errors in the graphs.
Cheers!
Javier

Code: Select all

import numpy as np
import matplotlib as plt
from pylab import *

from neuron import h
Section = h.Section

# initiate topology ############
soma = Section()
axon = Section()

axon.connect(soma , 0, 0)

# define geometry                   #############
soma.L    = 10
soma.diam = 10
soma.nseg = 1

axon.L    = 200
axon.nseg = 200
axon.diam = 0.6

# biophysics                        #############
# soma 
soma.insert('pas')
soma.g_pas = 0.00002
soma.e_pas = -70

axon.Ra   = 100

axon.insert('pas')
axon.g_pas = 0.00002
axon.e_pas = -70.

axon.insert('hh')
axon.gnabar_hh = 0.05
axon.gkbar_hh  = 0.4
axon.ena       = 50
axon.ek        = -70
axon.el_hh     = -70
#
axon.insert('nap')
axon.gbar_nap  = 0.006
#
axon.insert('Ih')
axon.gkhbar_Ih = 0.001

# Izap                              #############
stim = h.Izap(0.5, sec=soma)
setattr(stim, 'amp', 0.03)
setattr(stim, 'del', 50)
setattr(stim, 'dur', 50)
setattr(stim, 'f0', 0.)
setattr(stim, 'f1', 100)

# record variables                  #############
v_soma = h.Vector()
tt = h.Vector()
injCurr = h.Vector()
freq = h.Vector()

v_soma.record(soma(0.5)._ref_v)
tt.record(h._ref_t)
injCurr.record(stim._ref_i)
freq.record(stim._ref_f)

# simulation control                ###########
h.load_file("stdrun.hoc")
h.tstop  = 200
v_init = -65

h.init()
h.finitialize(v_init)
h.fcurrent()
h.run()


print 'del: '+str(getattr(stim, 'del'))
print 'dur: '+str(getattr(stim, 'dur'))
print 'f0: '+str(getattr(stim, 'f0'))
print 'f1: '+str(getattr(stim, 'f1'))
print 'Final freq recorded from electrode: '+ str(max(freq))
plt.figure(0)
plt.plot(tt, v_soma)
plt.figure(1)
plt.plot(tt, injCurr)
plt.figure(2)
plt.plot(tt,freq)
plt.show()

[ted]

Strange--thought I uploaded the fix for that almost a year ago. Thanks for bringing it to my attention. Download the latest izap.zip from the same URL and let me know if that takes care of the problem.

[javiZ]

It seems to work now. May I have done something wrong with it?
Anyway, thanks a lot Ted.
Javier