1# -*- coding: utf-8 -*- 2# 3# balancedneuron.py 4# 5# This file is part of NEST. 6# 7# Copyright (C) 2004 The NEST Initiative 8# 9# NEST is free software: you can redistribute it and/or modify 10# it under the terms of the GNU General Public License as published by 11# the Free Software Foundation, either version 2 of the License, or 12# (at your option) any later version. 13# 14# NEST is distributed in the hope that it will be useful, 15# but WITHOUT ANY WARRANTY; without even the implied warranty of 16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17# GNU General Public License for more details. 18# 19# You should have received a copy of the GNU General Public License 20# along with NEST. If not, see <http://www.gnu.org/licenses/>. 21 22""" 23Balanced neuron example 24----------------------- 25 26This script simulates a neuron driven by an excitatory and an 27inhibitory population of neurons firing Poisson spike trains. The aim 28is to find a firing rate for the inhibitory population that will make 29the neuron fire at the same rate as the excitatory population. 30 31Optimization is performed using the ``bisection`` method from Scipy, 32simulating the network repeatedly. 33 34This example is also shown in the article [1]_ 35 36References 37~~~~~~~~~~ 38 39.. [1] Eppler JM, Helias M, Mulller E, Diesmann M, Gewaltig MO (2009). PyNEST: A convenient interface to the NEST 40 simulator, Front. Neuroinform. 41 http://dx.doi.org/10.3389/neuro.11.012.2008 42 43""" 44 45############################################################################### 46# First, we import all necessary modules for simulation, analysis and 47# plotting. Scipy should be imported before nest. 48 49from scipy.optimize import bisect 50 51import nest 52import nest.voltage_trace 53import matplotlib.pyplot as plt 54 55############################################################################### 56# Additionally, we set the verbosity using ``set_verbosity`` to 57# suppress info messages. 58 59 60nest.set_verbosity("M_WARNING") 61nest.ResetKernel() 62 63############################################################################### 64# Second, the simulation parameters are assigned to variables. 65 66 67t_sim = 25000.0 # how long we simulate 68n_ex = 16000 # size of the excitatory population 69n_in = 4000 # size of the inhibitory population 70r_ex = 5.0 # mean rate of the excitatory population 71r_in = 20.5 # initial rate of the inhibitory population 72epsc = 45.0 # peak amplitude of excitatory synaptic currents 73ipsc = -45.0 # peak amplitude of inhibitory synaptic currents 74d = 1.0 # synaptic delay 75lower = 15.0 # lower bound of the search interval 76upper = 25.0 # upper bound of the search interval 77prec = 0.01 # how close need the excitatory rates be 78 79############################################################################### 80# Third, the nodes are created using ``Create``. We store the returned 81# handles in variables for later reference. 82 83neuron = nest.Create("iaf_psc_alpha") 84noise = nest.Create("poisson_generator", 2) 85voltmeter = nest.Create("voltmeter") 86spikerecorder = nest.Create("spike_recorder") 87 88################################################################################### 89# Fourth, the ``poisson_generator`` (`noise`) is configured. 90# Note that we need not set parameters for the neuron, the spike recorder, and 91# the voltmeter, since they have satisfactory defaults. 92 93noise.rate = [n_ex * r_ex, n_in * r_in] 94 95############################################################################### 96# Fifth, the ``iaf_psc_alpha`` is connected to the ``spike_recorder`` and the 97# ``voltmeter``, as are the two Poisson generators to the neuron. The command 98# ``Connect`` has different variants. Plain `Connect` just takes the handles of 99# pre- and postsynaptic nodes and uses the default values for weight and 100# delay. It can also be called with a list of weights, as in the connection 101# of the noise below. 102# Note that the connection direction for the ``voltmeter`` is reversed compared 103# to the ``spike_recorder``, because it observes the neuron instead of 104# receiving events from it. Thus, ``Connect`` reflects the direction of signal 105# flow in the simulation kernel rather than the physical process of inserting 106# an electrode into the neuron. The latter semantics is presently not 107# available in NEST. 108 109 110nest.Connect(neuron, spikerecorder) 111nest.Connect(voltmeter, neuron) 112nest.Connect(noise, neuron, syn_spec={'weight': [[epsc, ipsc]], 'delay': 1.0}) 113 114############################################################################### 115# To determine the optimal rate of the neurons in the inhibitory population, 116# the network is simulated several times for different values of the 117# inhibitory rate while measuring the rate of the target neuron. This is done 118# by calling ``Simulate`` until the rate of the target neuron matches the rate 119# of the neurons in the excitatory population with a certain accuracy. The 120# algorithm is implemented in two steps: 121# 122# First, the function ``output_rate`` is defined to measure the firing rate 123# of the target neuron for a given rate of the inhibitory neurons. 124 125 126def output_rate(guess): 127 print("Inhibitory rate estimate: %5.2f Hz" % guess) 128 rate = float(abs(n_in * guess)) 129 noise[1].rate = rate 130 spikerecorder.n_events = 0 131 nest.Simulate(t_sim) 132 out = spikerecorder.n_events * 1000.0 / t_sim 133 print(" -> Neuron rate: %6.2f Hz (goal: %4.2f Hz)" % (out, r_ex)) 134 return out 135 136 137############################################################################### 138# The function takes the firing rate of the inhibitory neurons as an 139# argument. It scales the rate with the size of the inhibitory population and 140# configures the inhibitory Poisson generator (`noise[1]`) accordingly. 141# Then, the spike counter of the ``spike_recorder`` is reset to zero. The 142# network is simulated using ``Simulate``, which takes the desired simulation 143# time in milliseconds and advances the network state by this amount of time. 144# During simulation, the ``spike_recorder`` counts the spikes of the target 145# neuron and the total number is read out at the end of the simulation 146# period. The return value of ``output_rate()`` is the firing rate of the 147# target neuron in Hz. 148# 149# Second, the scipy function ``bisect`` is used to determine the optimal 150# firing rate of the neurons of the inhibitory population. 151 152in_rate = bisect(lambda x: output_rate(x) - r_ex, lower, upper, xtol=prec) 153print("Optimal rate for the inhibitory population: %.2f Hz" % in_rate) 154 155############################################################################### 156# The function ``bisect`` takes four arguments: first a function whose 157# zero crossing is to be determined. Here, the firing rate of the target 158# neuron should equal the firing rate of the neurons of the excitatory 159# population. Thus we define an anonymous function (using `lambda`) that 160# returns the difference between the actual rate of the target neuron and the 161# rate of the excitatory Poisson generator, given a rate for the inhibitory 162# neurons. The next two arguments are the lower and upper bound of the 163# interval in which to search for the zero crossing. The fourth argument of 164# ``bisect`` is the desired relative precision of the zero crossing. 165# 166# Finally, we plot the target neuron's membrane potential as a function of 167# time. 168 169nest.voltage_trace.from_device(voltmeter) 170plt.show() 171