## sinusoid current

Anything that doesn't fit elsewhere.
kazmer
Posts: 10
Joined: Thu May 17, 2007 2:37 pm

### sinusoid current

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
Posts: 9
Joined: Tue Jun 07, 2005 10:00 pm
Contact:

### Sinusoidal current injection

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)
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
Posts: 10
Joined: Thu May 17, 2007 2:37 pm
Thanks, it is fine.
Kazmer

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

### Re: Sinusoidal current injection

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
Posts: 1
Joined: Thu Jun 14, 2007 1:11 pm
Contact:

### Modification -- swept sine wave (ZAP function)

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)
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
Posts: 5727
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: swept sine (ZAP function)

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

: 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
Posts: 5727
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: swept sine

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
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

### Negative t values on Current Axis using CVODE

Hi,

I was using the Izap.mod file on my model and was getting a few unexpected negative t values when plotting Izap.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.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
Posts: 5727
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: sinusoid current

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
Posts: 5727
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: sinusoid current

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

### Re: sinusoid current

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.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
Posts: 5727
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: sinusoid current

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

### Re: sinusoid current

It seems to work now. May I have done something wrong with it?
Anyway, thanks a lot Ted.
Javier
{"ft_lang_label":"__label__en","ft_lang_prob":0.83020186,"math_prob":0.86898386,"size":17578,"snap":"2020-24-2020-29","text_gpt3_token_len":5077,"char_repetition_ratio":0.093604185,"word_repetition_ratio":0.54083663,"special_character_ratio":0.28985095,"punctuation_ratio":0.1370768,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9780199,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2020-05-31T17:23:00Z\",\"WARC-Record-ID\":\"<urn:uuid:0734a5f6-83db-46a8-a42f-585e9a31d9f4>\",\"Content-Length\":\"72554\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:a97813ae-d41d-4b80-94c6-297c0faf494b>\",\"WARC-Concurrent-To\":\"<urn:uuid:85d5586b-637e-4f14-849e-a23bdd88f3ea>\",\"WARC-IP-Address\":\"128.36.111.212\",\"WARC-Target-URI\":\"https://www.neuron.yale.edu/phpBB/viewtopic.php?f=8&t=897\",\"WARC-Payload-Digest\":\"sha1:CGT2EQXCQ4OXFCZPJDH67U26L2IHMHFL\",\"WARC-Block-Digest\":\"sha1:4OW7D3JFSTACMKC2Y6NNKLFVW476CN57\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2020/CC-MAIN-2020-24/CC-MAIN-2020-24_segments_1590347413551.52_warc_CC-MAIN-20200531151414-20200531181414-00468.warc.gz\"}"}