Merge pull request #623 from janhahne/improve_comments

Improved comments in gap_junction Pynest examples

extras/userdoc/md/documentation/simulations-with-gap-junctions.md

@@ -8,7 +8,7 @@ due to several improvements in terms of usability.
 
 ## Introduction
 
-Simulations with gap junctions are supported by the Hodgin-Huxley neuron model `hh_psc_alpha_gap`.
+Simulations with gap junctions are supported by the Hodgkin-Huxley neuron model `hh_psc_alpha_gap`.
 The synapse model to create a gap-junction connection is named `gap_junction`. 
 Unlike chemical synapses gap junctions are bidirectional connections. 
 In order to create **one** accurate gap-junction connection 

pynest/examples/gap_junctions_inhibitory_network.py

@@ -22,29 +22,27 @@
 """
 Gap Junctions: Inhibitory network example
 ------------------
-This is the inhibitory network used as test case 2 (see figure 9 and 10) in
-
-    Hahne, J., Helias, M., Kunkel, S., Igarashi, J.,
-    Bolten, M., Frommer, A. and Diesmann, M.,
-    A unified framework for spiking and gap-junction interactions
-    in distributed neuronal network simulations,
-    Front. Neuroinform. 9:22. (2015),
-    doi: 10.3389/fninf.2015.00022
-
-The network contains 500 hh_psc_alpha_gap neurons with random initial
-membrane potentials between -40 and -80 mV. Each neuron receives 50
-inhibitory synaptic inputs that are randomly selected from all other
-neurons, each with synaptic weight JI = -50.0 pA and synaptic delay
-d = 1.0 ms. Each neuron receives an excitatory external Poissonian
-input of 500.0 Hz with synaptic weight JE = 300.0 pA and the same
-delay d. In addition (60*500)/2 gap junctions are added randomly to the
-network resulting in an average of 60 gap-junction connections per neuron.
+This script simulates an inhibitory network of 500 Hodgkin-Huxley neurons.
+Without the gap junctions (meaning for `gap_weight = 0.0`) the network
+shows an asynchronous irregular state that is caused by the external
+excitatory Poissonian drive being balanced by the inhibitory feedback
+within the network. With increasing `gap_weight` the network synchronizes:
+For a lower gap weight of 0.3 nS the network remains in an asynchronous
+state. With a weight of 0.54 nS the network switches randomly between the
+asynchronous to the synchronous state, while for a gap weight of 0.7 nS
+a stable synchronous state is reached.
+
+This example is also used as test case 2 (see figure 9 and 10)
+in Hahne et al. (2015)
+**A unified framework for spiking and gap-junction interactions
+in distributed neuronal network simulations**, *Front. Neuroinform.*
+http://dx.doi.org/10.3389/neuro.11.012.2008
 """
 
-import pylab
 import nest
-import random
+import pylab as pl
 import numpy
+import random
 
 n_neuron = 500
 gap_per_neuron = 60
@@ -55,15 +53,16 @@
 threads = 8
 stepsize = 0.05
 simtime = 501.
+gap_weight = 0.3
 
+nest.ResetKernel()
 
-#Set gap weight here
-gap_weight = 0.32
-
+"""
+First we set the random seed, adjust the kernel settings and create
+`hh_psc_alpha_gap` neurons, `spike_detector` and `poisson_generator`.
+"""
 random.seed(1)
 
-nest.ResetKernel()
-
 nest.SetKernelStatus({'resolution': 0.05,
                       'total_num_virtual_procs': threads,
                       'print_time': True,
@@ -82,6 +81,15 @@
                                            'to_memory': True})
 pg = nest.Create("poisson_generator", params={'rate': 500.0})
 
+"""
+Each neuron shall receive `inh_per_neuron = 50` inhibitory synaptic
+inputs that are randomly selected from all other neurons, each
+with synaptic weight `j_inh = -50.0` pA and a synaptic delay
+of 1.0 ms. Furthermore each neuron shall receive an excitatory
+external Poissonian input of 500.0 Hz with synaptic weight
+`j_exc = 300.0` pA and the same delay.
+The desired connections are created with the following commands:
+"""
 conn_dict = {'rule': 'fixed_indegree',
              'indegree': inh_per_neuron,
              'autapses': False,
@@ -97,29 +105,38 @@
                                                   'weight': j_exc,
                                                   'delay': delay})
 
+"""
+Then the neurons are connected to the `spike_detector` and the initial
+membrane potential of each neuron is set randomly between -40 and -80 mV.
+"""
 nest.Connect(neurons, sd)
 
 for i in range(n_neuron):
     nest.SetStatus([neurons[i]], {'V_m': (-40. - 40. * random.random())})
 
 """
-We must not use the 'fixed_indegree' oder 'fixed_outdegree' functionality of
-nest.Connect() to create the connections, as gap_junction connections are two-
-way connections and we need to make sure that the same neurons are connected
-in both ways (using the "make_symmetric" flag for one-to-one connections).
+Finally gap junctions are added to the network.
+(60*500)/2 `gap_junction` connections are added randomly
+resulting in an average of 60 gap-junction connections per neuron.
+We must not use the `fixed_indegree` oder `fixed_outdegree` functionality of
+`nest.Connect()` to create the connections, as `gap_junction` connections are
+bidirectional connections and we need to make sure that the same neurons
+are connected in both ways. This is achieved by creating the connections on
+the Python level with the `random` module of the Python Standard Library
+and connecting the neurons using the `make_symmetric` flag for
+`one_to_one` connections.
 """
-
-# create gap_junction connections
 n_connection = n_neuron * gap_per_neuron / 2
 connections = numpy.transpose(
     [random.sample(neurons, 2) for _ in range(n_connection)])
 
-# Connect sources -> targets and targets -> sources with
-# one call to nest.Connect using the "make_symmetric" flag
 nest.Connect(connections[0], connections[1],
              {'rule': 'one_to_one', 'make_symmetric': True},
              {'model': 'gap_junction', 'weight': gap_weight})
 
+"""
+In the end we start the simulation and plot the spike pattern.
+"""
 nest.Simulate(simtime)
 
 times = nest.GetStatus(sd, 'events')[0]['times']
@@ -128,9 +145,9 @@
 
 hz_rate = (1000.0 * n_spikes / simtime) / n_neuron
 
-pylab.figure(1)
-pylab.plot(times, spikes, 'o')
-pylab.title('Average spike rate (Hz): %.2f' % hz_rate)
-pylab.xlabel('time (ms)')
-pylab.ylabel('neuron no')
-pylab.show()
+pl.figure(1)
+pl.plot(times, spikes, 'o')
+pl.title('Average spike rate (Hz): %.2f' % hz_rate)
+pl.xlabel('time (ms)')
+pl.ylabel('neuron no')
+pl.show()

pynest/examples/gap_junctions_two_neurons.py

@@ -20,29 +20,25 @@
 # along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 
 """
-Gap Junctions: Two neurons example
+Gap Junctions: Two neuron example
 ------------------
-This is a simple example of two hh_psc_alpha_gap neurons connected
-by a gap-junction. Please note that gap junctions are two-way connections:
-In order to create an accurate gap-junction connection between two
-neurons i and j the use of the flag "make_symmetric" is required.
+This script simulates two Hodgkin-Huxley neurons of type `hh_psc_alpha_gap`
+connected by a gap junction. Both neurons receive a constant current of
+100.0 pA. The neurons are initialized with different membrane potentials and
+synchronize over time due to the gap-junction connection.
 """
 
 import nest
-import pylab
+import pylab as pl
 import numpy
 
 nest.ResetKernel()
 
-nest.SetKernelStatus({'resolution': 0.05,
-                      # Settings for waveform relaxation
-                      # 'use_wfr': False uses communication in every step
-                      # instead of an iterative solution
-                      'use_wfr': True,
-                      'wfr_comm_interval': 1.0,
-                      'wfr_tol': 0.0001,
-                      'wfr_max_iterations': 15,
-                      'wfr_interpolation_order': 3})
+"""
+First we set the resolution of the simulation, create two neurons and
+create a `voltmeter` for recording.
+"""
+nest.SetKernelStatus({'resolution': 0.05})
 
 neuron = nest.Create('hh_psc_alpha_gap', 2)
 
@@ -51,32 +47,39 @@
                                       'withtime': True,
                                       'interval': 0.1})
 
+"""
+Then we set the constant current input, modify the inital membrane
+potential of one of the neurons and connect the neurons to the `voltmeter`.
+"""
 nest.SetStatus(neuron, {'I_e': 100.})
 nest.SetStatus([neuron[0]], {'V_m': -10.})
 
 nest.Connect(vm, neuron, 'all_to_all')
 
 """
-Use the 'make_symmetric' flag to connect
-neuron[0] -> neuron[1] and neuron[1] -> neuron[0]
-with a single call to nest.Connect in order to create
-an accurate gap junction between both neurons
+In order to create the `gap_junction` connection we employ the `all_to_all`
+connection rule: Gap junctions are bidirectional connections, therefore we
+need to connect `neuron[0]` to `neuron[1]` and `neuron[1]` to `neuron[0]`:
 """
-nest.Connect([neuron[0]], [neuron[1]],
-             {'rule': 'one_to_one', 'make_symmetric': True},
+nest.Connect(neuron, neuron,
+             {'rule': 'all_to_all', 'autapses': False},
              {'model': 'gap_junction', 'weight': 0.5})
 
+"""
+Finally we start the simulation and plot the membrane potentials of
+both neurons.
+"""
 nest.Simulate(351.)
 
 senders = nest.GetStatus(vm, 'events')[0]['senders']
 times = nest.GetStatus(vm, 'events')[0]['times']
 V = nest.GetStatus(vm, 'events')[0]['V_m']
 
-pylab.figure(1)
-pylab.plot(times[numpy.where(senders == 1)],
-           V[numpy.where(senders == 1)], 'r-')
-pylab.plot(times[numpy.where(senders == 2)],
-           V[numpy.where(senders == 2)], 'g-')
-pylab.xlabel('time (ms)')
-pylab.ylabel('membrane potential (mV)')
-pylab.show()
+pl.figure(1)
+pl.plot(times[numpy.where(senders == 1)],
+        V[numpy.where(senders == 1)], 'r-')
+pl.plot(times[numpy.where(senders == 2)],
+        V[numpy.where(senders == 2)], 'g-')
+pl.xlabel('time (ms)')
+pl.ylabel('membrane potential (mV)')
+pl.show()