www.neuron.yale.edu

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

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

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

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

Code: Select all

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

Code: Select all

cell.L = 200
cell.diam = 6

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_

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