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