Response properties of neocort. neurons to temporally modulated noisy inputs (Koendgen et al. 2008)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:118631
Neocortical neurons are classified by current–frequency relationship. This is a static description and it may be inadequate to interpret neuronal responses to time-varying stimuli. Theoretical studies (Brunel et al., 2001; Fourcaud-Trocmé et al. 2003; Fourcaud-Trocmé and Brunel 2005; Naundorf et al. 2005) suggested that single-cell dynamical response properties are necessary to interpret ensemble responses to fast input transients. Further, it was shown that input-noise linearizes and boosts the response bandwidth, and that the interplay between the barrage of noisy synaptic currents and the spike-initiation mechanisms determine the dynamical properties of the firing rate. In order to allow a reader to explore such simulations, we prepared a simple NEURON implementation of the experiments performed in Köndgen et al., 2008 (see also Fourcaud-Trocmé al. 2003; Fourcaud-Trocmé and Brunel 2005). In addition, we provide sample MATLAB routines for exploring the sandwich model proposed in Köndgen et al., 2008, employing a simple frequdency-domain filtering. The simulations and the MATLAB routines are based on the linear response properties of layer 5 pyramidal cells estimated by injecting a superposition of a small-amplitude sinusoidal wave and a background noise, as in Köndgen et al., 2008.
Reference:
1 . Koendgen H, Geisler C, Wang XJ, Fusi S, Luescher HR, Giugliano M (2004) The dynamical response of single cells to noisy time-varying currents Soc Neurosci Abstr :640
2 . Köndgen H, Geisler C, Fusi S, Wang XJ, Lüscher HR, Giugliano M (2008) The dynamical response properties of neocortical neurons to temporally modulated noisy inputs in vitro. Cereb Cortex 18:2086-97 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Axon;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 L6 pyramidal corticothalamic GLU cell; Abstract Wang-Buzsaki neuron;
Channel(s): I Na,t; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; MATLAB;
Model Concept(s): Parameter Fitting; Methods; Rate-coding model neurons;
Implementer(s): Giugliano, Michele [mgiugliano at gmail.com]; Delattre, Vincent;
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic GLU cell; I Na,t; I K;
//
// DEMO
// by M. Giugliano and V. Delattre(2008)
// Brain Mind Institute, EPFL
//


load_file("nrngui.hoc")                   // Loading of the standard GUI controls...
load_file("mylibs/graphs.hoc")            // Loading some ad-hoc proc for displaying live traces...
load_file("mylibs/Functions.hoc")         // Loading some had-hoc function: mysin, err...

//------------------------------------------------------------------------------------------------------


print "2008 - Michele Giugliano & Vincent Delattre"
print "Brain Mind Institute, EPFL of Lausanne"
print "Demonstrating: 1) noisy stimulation, modulated in time"
print "               2) instantaneous firing rate, estimate and quantitative fit"
print ""


freq   = 0.3			// Oscillation frequency of the input sinusoid
period = 1./freq		// By definition, the period of the oscillation 
T      = 6000.			//[ms], overall duration of the simulation (i.e. should be much larger than (1000*period))
//
//------------------------------------------------------------------------------------------------------------------------
//
print "Let's first create a single-compartmental, conductance-based model neuron.."
create soma          // This creates a single-compartmental neuron..
soma {             
 L     = 65          // [um] length of a cilindrical cable
 diam  = 65          // [um] diameter of a cilindrical cable
 nseg  = 1           // number of segments - see the documentation

 insert pas          // inserts passive mechanisms (i.e. leak currents, voltage-indep.)
 g_pas = 1e-4
 e_pas = -67.        // [mV]      - Leak currents, reversal potential
 cm    = 1.          // [uF/cm^2] - specific membrane capacitance
 Ra    = 35.4        // [ohm cm]  - axial/cytosolic resistivity - useless in a 1-compartmental model

 insert wb			 // inserts active mechanisms (Na, K currents, voltage-dep.)
 gnabar  = .007  //
 gkbar   = .009  //
} // end soma
print "done!"
print ""
//
//------------------------------------------------------------------------------------------------------------------------
//
print "Let's now create and setup a current-clamp stimulation.."
objref fl
fl          = new Isinunoisy(0.5)     // Sinusoidal + Noisy, current-clamp stimulation (modulation on the *mean*); see also the mechanism Isinunoisy2 (modulation on the *standard deviation*)
fl.m        = 0.3					  // [nA], DC offset of the overall current 
fl.s        = 0.2					  // [nA], square root of the steady-state variance of the (noisy) stimulus component
fl.tau      = 2.5					  // [ms], steady-state correlation time-length of the (noisy) stimulus component
fl.amp      = 0.25					  // [nA], amplitude of the (sinusoidal) stimulus component
fl.freq     = freq					  // [Hz], steady-state correlation time-length of the (noisy) stimulus component
print "done!"
print ""
//
//------------------------------------------------------------------------------------------------------------------------
//
print "Let's define a way to count and log the time of each spike that will be fire.."
objectvar apc                   // A.P.C. - Action Potential Count mechanism.
soma apc   = new APCount(0.5)   // Counting is triggered by somatic membrane potential.
apc.thresh = -20.               // [mV], spikes times are recorded as the upward crossing of this threshold
apc.time   = 10000000.          //
print "done!"
print ""
//
//------------------------------------------------------------------------------------------------------------------------
//
print "Time to launch the simulation!"

strdef aptime, linear
objref gg, outspikes, indepvar, myhist, fitsine, inputsine, vec

Vrest              = -67          // [mV]
tstart             = 0            // [ms] 
t                  = tstart
tstop              = T            // [ms] Total simulation duration [ms]
dt                 = 0.025       // ms - integration time step

outspikes = new Vector()
apc.record(outspikes)   // This statement attaches to the spike counter 'myapc' a vector to collect spike times..

//
//----------- GRAPH DISPLAYING ETC. ETC. ----------------------------------------------------------------------------
//
objectvar save_window_, rvp_
objectvar scene_vector_[4]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
objectvar g0, g1

finitialize(Vrest)

proc graph_and_run() {
  {ocbox_list_ = new List()  scene_list_ = new List()}
  {pwman_place(0,0,0)}
  {
  save_window_ = new Graph(0)
  save_window_.size(0,tstop,-80,40)
  scene_vector_[2] = save_window_
  {save_window_.view(0, -80, tstop, 120, 707, 22, 625, 395)}
  graphList[0].append(save_window_)
  save_window_.save_name("graphList[0].")
  save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)
  g0=save_window_
  }
  {
  save_window_ = new Graph(0)
  save_window_.size(0,tstop,-0.8,1.2)
  scene_vector_[3] = save_window_
  {save_window_.view(0, 0, tstop, 1, 708, 445, 626, 394)}
  graphList[2].append(save_window_)
  save_window_.save_name("graphList[2].")
  save_window_.addvar("fl.i", 1, 1, 0.8, 0.9, 2)
  g1=save_window_
  }
//
//------------------------------------------------------------------------------------------------------------------------
//


  finitialize(Vrest)
  run()
  N      = apc.n          // N spikes have been counted.
  R0     = 1000 * N / T   // [Hz], Mean Firing rate
  print "done!"
  print ""

//
//------------------------------------------------------------------------------------------------------------------------
//

  i = 0
  while (i < N) {                               // I go through each spike time, one by one..
          spk = outspikes.x[i] / 1000
          tmp = int(spk * freq)                 // tmp now contains the integer part of the division spk /period..
          outspikes.x[i] = spk - tmp / freq     // Then this makes it as the reminder of the division spk / period (by definition)..
          i = i + 1                             // Next spike incrementing the current index 'I'...
  } // end while()


////////////////////////////BUILDING HISTOGRAM///////////////////////////////////////////////////////////////////

  myhist = new Vector()                     // contains the histogram
  R0 = 1000 * N / T                         // This is the mean firing rate [spikes/s]   
  Ncycles = 0.001 * freq * T                // This will be used for the normalization
  binW = 1 / freq * 1 / 33                  // Bin width
  binL = 0                                  // Low born of the histogram

  myhist = myhist.hist(outspikes,binL,33,binW)                  // We want to histogram into an estimate of the instantaneous firing rate, so we need to make a normalization..
  myhist.div(Ncycles * binW)                                    // 'div' is a fast way to divide each member of the vector by the same normalization factor: Ncycles * width.
  Mbins = myhist.size()                                         // Number of bins

  // Now elements are measured in [Hz]..

  indepvar = new Vector()                                       // Let's now create an index vector (our x-axis) with appropriately
  indepvar.indgen(binW * .5, 1 / freq - binW * 0.5, binW)       // spaced points..

///////////////////////////////////FITTING A FUNCTION////////////////////////////////////////////////////////////
  inputsine = new Vector(Mbins)                                                // The input sinusoid
  fitsine = new Vector(Mbins)                                                  // The best fitted sinusoid
  vec     = new Vector(3)                                                      // Contains the parameters to be fitted
  vec.x[0]= (myhist.max()-myhist.min())                                        // The amplitude r1         
  vec.x[1]= indepvar.x[myhist.max_ind()] * (360. * freq) - 90.                 // The phase ph
  vec.x[2]= R0                                                                 // The offset

  attr_praxis(0.01, 1000, 0)            // Set the fitting attributes: tolerance...
  fit_praxis(3,"err",&vec.x[0])         // Do the fitting

            r1 = vec.x[0]
            ph = vec.x[1]
            ro = vec.x[2]

  if (r1 < 0) {                         //
   r1 = -r1                          // You never know what the fitting procedure will choose
   ph = ph - 180.                    // and you want to refer to 'r1' as a positive number..
  }

   i = 0
   while (i < Mbins) {                               // I go through each spike time, one by one..
        fitsine.x[i] = mysin(i*binW, r1, ph, R0)
		inputsine.x[i] = mysin(i*binW, r1, 0, R0)
		i = i + 1                             // Next spike incrementing the current index 'I'...
   } // end while()


  {
  gg = new Graph(0)
  gg.size(0,period,0,R0+r1)
  gg.view(0, 0, period, R0+r1, 1, 82, 700, 555)}

  fitsine.plot(gg, indepvar, 2,1)
  myhist.plot(gg, indepvar, 1,1)
  inputsine.plot(gg, indepvar, 3,1)

  print "In black, the instantaneous firing rate [Hz] is indicated (i.e. the PSTH)"
  print "In red, the best fit sinusoid (mean: ", R0, "[Hz], amplitude: ", r1, "[Hz], phase: ", ph
  print "In blue the sinusoid subsequently modified by noise to produce the input."
  g0.exec_menu("View = plot")  // zoom out to show everything within the graphs with these commands
  g1.exec_menu("View = plot")
  gg.exec_menu("View = plot")
}

xpanel("Launch panel")
  xbutton("Click here to run","graph_and_run()")
xpanel()