1# -*- coding: utf-8 -*-
2#
3# gap_junctions_inhibitory_network.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"""
23Gap Junctions: Inhibitory network example
24-----------------------------------------
25
26This script simulates an inhibitory network of 500 Hodgkin-Huxley neurons.
27Without the gap junctions (meaning for ``gap_weight = 0.0``) the network shows
28an asynchronous irregular state that is caused by the external excitatory
29Poissonian drive being balanced by the inhibitory feedback within the
30network. With increasing `gap_weight` the network synchronizes:
31
32For a lower gap weight of 0.3 nS the network remains in an asynchronous
33state. With a weight of 0.54 nS the network switches randomly between the
34asynchronous to the synchronous state, while for a gap weight of 0.7 nS a
35stable synchronous state is reached.
36
37This example is also used as test case 2 (see Figure 9 and 10)
38in [1]_.
39
40References
41~~~~~~~~~~
42
43.. [1] Hahne et al. (2015) A unified framework for spiking and gap-junction
44       interactions in distributed neuronal network simulations, Front.
45       Neuroinform. http://dx.doi.org/10.3389/neuro.11.012.2008
46"""
47
48import nest
49import matplotlib.pyplot as plt
50import numpy
51
52n_neuron = 500
53gap_per_neuron = 60
54inh_per_neuron = 50
55delay = 1.0
56j_exc = 300.
57j_inh = -50.
58threads = 8
59stepsize = 0.05
60simtime = 501.
61gap_weight = 0.3
62
63nest.ResetKernel()
64
65###############################################################################
66# First we set the random seed, adjust the kernel settings and create
67# ``hh_psc_alpha_gap`` neurons, ``spike_recorder`` and ``poisson_generator``.
68
69numpy.random.seed(1)
70
71nest.resolution = 0.05
72nest.total_num_virtual_procs = threads
73nest.print_time = True
74
75# Settings for waveform relaxation. If 'use_wfr' is set to False,
76# communication takes place in every step instead of using an
77# iterative solution
78nest.use_wfr = True
79nest.wfr_comm_interval = 1.0
80nest.wfr_tol = 0.0001
81nest.wfr_max_iterations = 15
82nest.wfr_interpolation_order = 3
83
84neurons = nest.Create('hh_psc_alpha_gap', n_neuron)
85
86sr = nest.Create("spike_recorder")
87pg = nest.Create("poisson_generator", params={'rate': 500.0})
88
89###############################################################################
90# Each neuron shall receive ``inh_per_neuron = 50`` inhibitory synaptic inputs
91# that are randomly selected from all other neurons, each with synaptic
92# weight ``j_inh = -50.0`` pA and a synaptic delay of 1.0 ms. Furthermore each
93# neuron shall receive an excitatory external Poissonian input of 500.0 Hz
94# with synaptic weight ``j_exc = 300.0`` pA and the same delay.
95# The desired connections are created with the following commands:
96
97conn_dict = {'rule': 'fixed_indegree',
98             'indegree': inh_per_neuron,
99             'allow_autapses': False,
100             'allow_multapses': True}
101
102syn_dict = {'synapse_model': 'static_synapse',
103            'weight': j_inh,
104            'delay': delay}
105
106nest.Connect(neurons, neurons, conn_dict, syn_dict)
107
108nest.Connect(pg, neurons, 'all_to_all',
109             syn_spec={'synapse_model': 'static_synapse',
110                       'weight': j_exc,
111                       'delay': delay})
112
113###############################################################################
114# Then the neurons are connected to the ``spike_recorder`` and the initial
115# membrane potential of each neuron is set randomly between -40 and -80 mV.
116
117nest.Connect(neurons, sr)
118
119neurons.V_m = nest.random.uniform(min=-80., max=-40.)
120
121#######################################################################################
122# Finally gap junctions are added to the network. :math:`(60*500)/2` ``gap_junction``
123# connections are added randomly resulting in an average of 60 gap-junction
124# connections per neuron. We must not use the ``fixed_indegree`` oder
125# ``fixed_outdegree`` functionality of ``nest.Connect()`` to create the
126# connections, as ``gap_junction`` connections are bidirectional connections
127# and we need to make sure that the same neurons are connected in both ways.
128# This is achieved by creating the connections on the Python level with the
129# `random` module of the Python Standard Library and connecting the neurons
130# using the ``make_symmetric`` flag for ``one_to_one`` connections.
131
132n_connection = int(n_neuron * gap_per_neuron / 2)
133neuron_list = neurons.tolist()
134connections = numpy.random.choice(neuron_list, [n_connection, 2])
135
136for source_node_id, target_node_id in connections:
137    nest.Connect(nest.NodeCollection([source_node_id]),
138                 nest.NodeCollection([target_node_id]),
139                 {'rule': 'one_to_one', 'make_symmetric': True},
140                 {'synapse_model': 'gap_junction', 'weight': gap_weight})
141
142###############################################################################
143# In the end we start the simulation and plot the spike pattern.
144
145nest.Simulate(simtime)
146
147events = sr.events
148times = events['times']
149spikes = events['senders']
150n_spikes = sr.n_events
151
152hz_rate = (1000.0 * n_spikes / simtime) / n_neuron
153
154plt.figure(1)
155plt.plot(times, spikes, 'o')
156plt.title(f'Average spike rate (Hz): {hz_rate:.2f}')
157plt.xlabel('time (ms)')
158plt.ylabel('neuron no')
159plt.show()
160