1# -*- coding: utf-8 -*-
2#
3# glif_psc_neuron.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"""
23Current-based generalized leaky integrate and fire (GLIF) neuron example
24------------------------------------------------------------------------
25
26Simple example of how to use the ``glif_psc`` neuron model for
27five different levels of GLIF neurons.
28
29Four stimulation paradigms are illustrated for the GLIF model
30with externally applied current and spikes impinging
31
32Voltage traces, current traces, threshold traces, and spikes are shown.
33
34KEYWORDS: glif_psc
35"""
36
37##############################################################################
38# First, we import all necessary modules to simulate, analyze and plot this
39# example.
40
41import nest
42import matplotlib.pyplot as plt
43import matplotlib.gridspec as gridspec
44
45##############################################################################
46# We initialize NEST and set the simulation resolution.
47
48nest.ResetKernel()
49resolution = 0.05
50nest.resolution = resolution
51
52##############################################################################
53# We also pre-define the synapse time constant array, [2.0, 1.0] ms for
54# the two desired synaptic ports of the GLIF neurons. Note that the default
55# synapse time constant is [2.0] ms, which is for neuron with one port.
56
57syn_tau = [2.0, 1.0]
58
59###############################################################################
60# We create the five levels of GLIF model to be tested, i.e.,
61# ``lif``, ``lif_r``, ``lif_asc``, ``lif_r_asc``, ``lif_r_asc_a``.
62# For each level of GLIF model, we create  a ``glif_psc`` node. The node is
63# created by setting relative model mechanism parameters and the time constant
64# of the 2 synaptic ports as mentioned above. Other neuron parameters are set
65# as default. The five ``glif_psc`` node handles were combined as a list.
66
67n_lif = nest.Create("glif_psc",
68                    params={"spike_dependent_threshold": False,
69                            "after_spike_currents": False,
70                            "adapting_threshold": False,
71                            "tau_syn": syn_tau})
72n_lif_r = nest.Create("glif_psc",
73                      params={"spike_dependent_threshold": True,
74                              "after_spike_currents": False,
75                              "adapting_threshold": False,
76                              "tau_syn": syn_tau})
77n_lif_asc = nest.Create("glif_psc",
78                        params={"spike_dependent_threshold": False,
79                                "after_spike_currents": True,
80                                "adapting_threshold": False,
81                                "tau_syn": syn_tau})
82n_lif_r_asc = nest.Create("glif_psc",
83                          params={"spike_dependent_threshold": True,
84                                  "after_spike_currents": True,
85                                  "adapting_threshold": False,
86                                  "tau_syn": syn_tau})
87n_lif_r_asc_a = nest.Create("glif_psc",
88                            params={"spike_dependent_threshold": True,
89                                    "after_spike_currents": True,
90                                    "adapting_threshold": True,
91                                    "tau_syn": syn_tau})
92
93neurons = n_lif + n_lif_r + n_lif_asc + n_lif_r_asc + n_lif_r_asc_a
94
95###############################################################################
96# For the stimulation input to the glif_psc neurons, we create one excitation
97# spike generator and one inhibition spike generator, each of which generates
98# three spikes; we also create one step current generator and a Poisson
99# generator, a parrot neuron (to be paired with the Poisson generator).
100# The three different injections are spread to three different time periods,
101# i.e., 0 ms ~ 200 ms, 200 ms ~ 500 ms, 600 ms ~ 900 ms.
102# Each of the excitation and inhibition spike generators generates three spikes
103# at different time points. Configuration of the current generator includes the
104# definition of the start and stop times and the amplitude of the injected
105# current. Configuration of the Poisson generator includes the definition of
106# the start and stop times and the rate of the injected spike train.
107
108espikes = nest.Create("spike_generator",
109                      params={"spike_times": [10., 100., 150.],
110                              "spike_weights": [20.] * 3})
111ispikes = nest.Create("spike_generator",
112                      params={"spike_times": [15., 99., 150.],
113                              "spike_weights": [-20.] * 3})
114cg = nest.Create("step_current_generator",
115                 params={"amplitude_values": [400., ],
116                         "amplitude_times": [200., ],
117                         "start": 200., "stop": 500.})
118pg = nest.Create("poisson_generator",
119                 params={"rate": 150000., "start": 600., "stop": 900.})
120pn = nest.Create("parrot_neuron")
121
122###############################################################################
123# The generators are then connected to the neurons. Specification of
124# the ``receptor_type`` uniquely defines the target receptor.
125# We connect current generator, the spike generators, Poisson generator (via
126# parrot neuron) to receptor 0, 1, and 2 of the GLIF neurons, respectively.
127# Note that Poisson generator is connected to parrot neuron to transit the
128# spikes to the glif_psc neuron.
129
130nest.Connect(cg, neurons, syn_spec={"delay": resolution})
131nest.Connect(espikes, neurons,
132             syn_spec={"delay": resolution, "receptor_type": 1})
133nest.Connect(ispikes, neurons,
134             syn_spec={"delay": resolution, "receptor_type": 1})
135nest.Connect(pg, pn, syn_spec={"delay": resolution})
136nest.Connect(pn, neurons, syn_spec={"delay": resolution, "receptor_type": 2})
137
138###############################################################################
139# A ``multimeter`` is created and connected to the neurons. The parameters
140# specified for the multimeter include the list of quantities that should be
141# recorded and the time interval at which quantities are measured.
142
143mm = nest.Create("multimeter",
144                 params={"interval": resolution,
145                         "record_from": ["V_m", "I", "I_syn", "threshold",
146                                         "threshold_spike",
147                                         "threshold_voltage",
148                                         "ASCurrents_sum"]})
149nest.Connect(mm, neurons)
150
151###############################################################################
152# A ``spike_recorder`` is created and connected to the neurons record the
153# spikes generated by the glif_psc neurons.
154
155sr = nest.Create("spike_recorder")
156nest.Connect(neurons, sr)
157
158###############################################################################
159# Run the simulation for 1000 ms and retrieve recorded data from
160# the multimeter and spike recorder.
161
162nest.Simulate(1000.)
163
164data = mm.events
165senders = data["senders"]
166
167spike_data = sr.events
168spike_senders = spike_data["senders"]
169spikes = spike_data["times"]
170
171###############################################################################
172# We plot the time traces of the membrane potential (in blue) and
173# the overall threshold (in green), and the spikes (as red dots) in one panel;
174# the spike component of threshold (in yellow) and the voltage component of
175# threshold (in black) in another panel; the injected currents (in strong blue),
176# the sum of after spike currents (in cyan), and the synaptic currents (in
177# magenta) in responding to the spike inputs to the neurons in the third panel.
178# We plot all these three panels for each level of GLIF model in a separated
179# figure.
180
181glif_models = ["lif", "lif_r", "lif_asc", "lif_r_asc", "lif_r_asc_a"]
182for i in range(len(glif_models)):
183
184    glif_model = glif_models[i]
185    node_id = neurons[i].global_id
186    plt.figure(glif_model)
187    gs = gridspec.GridSpec(3, 1, height_ratios=[2, 1, 1])
188    t = data["times"][senders == 1]
189
190    ax1 = plt.subplot(gs[0])
191    plt.plot(t, data["V_m"][senders == node_id], "b")
192    plt.plot(t, data["threshold"][senders == node_id], "g--")
193    plt.plot(spikes[spike_senders == node_id],
194             [max(data["threshold"][senders == node_id]) * 0.95] *
195             len(spikes[spike_senders == node_id]), "r.")
196    plt.legend(["V_m", "threshold", "spike"])
197    plt.ylabel("V (mV)")
198    plt.title("Simulation of glif_psc neuron of " + glif_model)
199
200    ax2 = plt.subplot(gs[1])
201    plt.plot(t, data["threshold_spike"][senders == node_id], "y")
202    plt.plot(t, data["threshold_voltage"][senders == node_id], "k--")
203    plt.legend(["threshold_spike", "threshold_voltage"])
204    plt.ylabel("V (mV)")
205
206    ax3 = plt.subplot(gs[2])
207    plt.plot(t, data["I"][senders == node_id], "--")
208    plt.plot(t, data["ASCurrents_sum"][senders == node_id], "c-.")
209    plt.plot(t, data["I_syn"][senders == node_id], "m")
210    plt.legend(["I_e", "ASCurrents_sum", "I_syn"])
211    plt.ylabel("I (pA)")
212    plt.xlabel("t (ms)")
213
214plt.show()
215