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