nsdm_hill_tononi_Vp_reset.py

# This tutorial shows you how to implement a simplified version of the
# Hill-Tononi model of the early visual pathway using the NEST Topology
# module.  The model is described in the paper
#
#   S. L. Hill and G. Tononi.
#   Modeling Sleep and Wakefulness in the Thalamocortical System.
#   J Neurophysiology **93**:1671-1698 (2005).
#   Freely available via `doi 10.1152/jn.00915.2004
#   <http://dx.doi.org/10.1152/jn.00915.2004>`_.
#
# And is based on the NEST Topology example of the Hill-Tononi model
# implement by Hans Ekkehard Plesser, which can be found here:
#
# https://github.com/nest/nest-simulator/blob/master/topology/examples/hill_tononi_Vp.py
#
# Authors: Leonardo Barbosa and Keiko Fujii
# Date:    21 September 2016

import os.path

SHOW_FIGURES = False

import pylab
if not SHOW_FIGURES:
    pylab_show = pylab.show
    def nop(s=None): pass
    pylab.show = nop
else:
    pylab.ion()

import matplotlib.pyplot as plt

# Leonardo: I need to import readline for some reason, otherwise nest wont work
import readline

# Load pynest
import nest
import nest.topology as topo
from nest import raster_plot



# Other generic imports
import math
import numpy as np
import scipy.io
import pickle

# --- keiko 9/14/2016
# Import csv to save V_m data
import csv
#f_Vm = open('V_m.csv', 'w')
#writer_Vm = csv.writer(f_Vm)
#f_idx = open('idx.csv', 'w')
#writer_idx = csv.writer(f_idx)

# Configurable Parameters
# =======================


# Number of simulations

Nsim = 1

lambdas = [2.0, -1.0]
#lambdas = [-1.0]

for lambda_dg in lambdas:

    # Make sure we start with a clean slate, even if we re-run the script
    # in the same Python session.
    nest.ResetKernel()

    # set multi-thread on
    N_vp = 16
    nest.SetKernelStatus({"local_num_threads": N_vp})

    # initialize random seed
    import time
    msd = int(round(time.time() * 1000))
    nest.SetKernelStatus({'grng_seed' : msd})
    nest.SetKernelStatus({'rng_seeds' : range(msd+N_vp+1, msd+2*N_vp+1)})

    #
    # Here we define those parameters that we take to be
    # configurable. The choice of configurable parameters is obviously
    # arbitrary, and in practice one would have far more configurable
    # parameters. We restrict ourselves to:
    #
    # - Network size in neurons ``N``, each layer is ``N x N``.
    # - Network size in subtended visual angle ``visSize``, in degree.
    # - Temporal frequency of drifting grating input ``f_dg``, in Hz.
    # - Spatial wavelength and direction of drifting grating input,
    #   ``lambda_dg`` and ``phi_dg``, in degree/radian.
    # - Background firing rate of retinal nodes and modulation amplitude,
    #   ``retDC`` and ``retAC``, in Hz.
    # - Simulation duration ``simtime``; actual simulation is split into
    #   intervals of ``sim_interval`` length, so that the network state
    #   can be visualized in those intervals. Times are in ms.
    Params = {'N'           :     40,
              'visSize'     :    8.0,
              'f_dg'        :    0.5,
              'lambda_dg'   :    lambda_dg, # spatial structure
              'bar'         :  True,
              #'bar'         :  False,
              'phi_dg'      :    0.0,
              'retDC'       :   10.0,
              'retAC'       :   10.0,
              'simtime'     : 5000.0,
              'sim_interval':    1.0,
              'discardtime' :    0.0
              }


    # Neuron Models
    # =============
    #
    # We declare models in two steps:
    #
    # 1. We define a dictionary specifying the NEST neuron model to use
    #    as well as the parameters for that model.
    # #. We create three copies of this dictionary with parameters
    #    adjusted to the three model variants specified in Table~2 of
    #    Hill & Tononi (2005) (cortical excitatory, cortical inhibitory,
    #    thalamic)
    #
    # In addition, we declare the models for the stimulation and
    # recording devices.
    #
    # The general neuron model
    # ------------------------
    #
    # We use the ``iaf_cond_alpha`` neuron, which is an
    # integrate-and-fire neuron with two conductance-based synapses which
    # have alpha-function time course.  Any input with positive weights
    # will automatically directed to the synapse labeled ``_ex``, any
    # with negative weights to the synapes labeled ``_in``.  We define
    # **all** parameters explicitly here, so that no information is
    # hidden in the model definition in NEST. ``V_m`` is the membrane
    # potential to which the model neurons will be initialized.
    # The model equations and parameters for the Hill-Tononi neuron model
    # are given on pp. 1677f and Tables 2 and 3 in that paper. Note some
    # peculiarities and adjustments:
    #
    # - Hill & Tononi specify their model in terms of the membrane time
    #   constant, while the ``iaf_cond_alpha`` model is based on the
    #   membrane capcitance. Interestingly, conducantces are unitless in
    #   the H&T model. We thus can use the time constant directly as
    #   membrane capacitance.
    # - The model includes sodium and potassium leak conductances. We
    #   combine these into a single one as follows:
    #$   \begin{equation}-g_{NaL}(V-E_{Na}) - g_{KL}(V-E_K)
    #$      =
    #$   -(g_{NaL}+g_{KL})\left(V-\frac{g_{NaL}E_{NaL}+g_{KL}E_K}{g_{NaL}g_{KL}}\right)
    #$   \end{equation}
    # - We write the resulting expressions for g_L and E_L explicitly
    #   below, to avoid errors in copying from our pocket calculator.
    # - The paper gives a range of 1.0-1.85 for g_{KL}, we choose 1.5
    #   here.
    # - The Hill-Tononi model has no explicit reset or refractory
    #   time. We arbitrarily set V_reset and t_ref.
    # - The paper uses double exponential time courses for the synaptic
    #   conductances, with separate time constants for the rising and
    #   fallings flanks. Alpha functions have only a single time
    #   constant: we use twice the rising time constant given by Hill and
    #   Tononi.
    # - In the general model below, we use the values for the cortical
    #   excitatory cells as defaults. Values will then be adapted below.
    #
    nest.CopyModel('iaf_cond_alpha', 'NeuronModel',
                   params = {'C_m'       :  16.0,
                             'E_L'       : (0.2 * 30.0 + 1.5 * -90.0)/(0.2 + 1.5),
                             'g_L'       : 0.2 + 1.5,
                             #'E_L'       : 30.0,
                             #'g_L'       : 0.2,
                             'E_ex'      :   0.0,
                             'E_in'      : -70.0,
                             'V_reset'   : -60.0,
                             'V_th'      : -51.0,
                             't_ref'     :   2.0,
                             'tau_syn_ex':   1.0,
                             'tau_syn_in':   2.0,
                             'I_e'       :   0.0,
                             'V_m'       : -70.0})

    # Adaptation of models for different populations
    # ----------------------------------------------

    # We must copy the `NeuronModel` dictionary explicitly, otherwise
    # Python would just create a reference.

    # Cortical excitatory cells
    # .........................
    # Parameters are the same as above, so we need not adapt anything
    nest.CopyModel('NeuronModel', 'CtxExNeuron')

    # Cortical inhibitory cells
    # .........................
    nest.CopyModel('NeuronModel', 'CtxInNeuron',
                   params = {'C_m'  :   8.0,
                             'V_th' : -53.0,
                             't_ref':   1.0})

    # Thalamic cells
    # ..............
    nest.CopyModel('NeuronModel', 'ThalamicNeuron',
                   params = {'C_m'  :   8.0,
                             'V_th' : -53.0,
                             't_ref':   1.0,
                             'E_in' : -80.0})

    # Input generating nodes
    # ----------------------

    # Input is generated by sinusoidally modulate Poisson generators,
    # organized in a square layer of retina nodes. These nodes require a
    # slightly more complicated initialization than all other elements of
    # the network:
    #
    # - Average firing rate ``rate``, firing rate modulation depth ``amplitude``, and
    #   temporal modulation frequency ``frequency`` are the same for all retinal
    #   nodes and are set directly below.
    # - The temporal phase ``phase`` of each node depends on its position in
    #   the grating and can only be assigned after the retinal layer has
    #   been created. We therefore specify a function for initalizing the
    #   ``phase``. This function will be called for each node.
    def phaseInit(pos, lam, alpha):
        '''Initializer function for phase of drifting grating nodes.

           pos  : position (x,y) of node, in degree
           lam  : wavelength of grating, in degree
           alpha: angle of grating in radian, zero is horizontal

           Returns number to be used as phase of sinusoidal Poisson generator.
        '''
        return 360.0 / lam * (math.cos(alpha) * pos[0] + math.sin(alpha) * pos[1])

    def barInit(pos, w, a, s):
        '''Initializer function for the rate of a fixed bar

           pos  : position (x,y) of node, in degree
           lam  : width of the bar, in degree
           alpha: angle of the bar in radian, zero is horizontal
           vsize: size of the visual field, bar will be in the center

           Returns 0 if outside the bar, 1 otherwise
        '''

        max = w * (s*math.cos(a) + s*math.sin(a))
        inibar = -(w/s/2)
        endbar = +(w/s/2)

        x = pos[0]
        y = pos[1]
        A = w * (x*math.cos(a) + y*math.sin(a)) / max
        ok = (A > inibar) & (A < endbar)

        return 1*ok

    def phaseBarInit(pos, w, a, s):

        max = w * (s*math.cos(a) + s*math.sin(a))
        inibar = -(w/s/2)
        endbar = +(w/s/2)

        x = pos[0]
        y = pos[1]
        phase = w * (x*math.cos(a) + y*math.sin(a))
        A = phase / max
        ok = (A > inibar) & (A < endbar)

        if ok:
            return 360.0 / phase
        else:
            return 702.0

    if Params['bar']:
        nest.CopyModel('poisson_generator', 'RetinaNode',
                       params = {'rate'     : 0.0})
    else:
        nest.CopyModel('sinusoidal_poisson_generator', 'RetinaNode',
                       params = {'amplitude': Params['retAC'],
                                 'rate'     : Params['retDC'],
                                 'frequency': Params['f_dg'],
                                 'phase'    : 0.0,
                                 'individual_spike_trains': True})

    #nest.CopyModel('sinusoidal_poisson_generator', 'RetinaNode',
    #               params = {'amplitude': Params['retAC'],
    #                         'rate'     : Params['retDC'],
    #                         'frequency': Params['f_dg'],
    #                         'phase'    : 0.0,
    #                         'individual_spike_trains': True})

    # Recording nodes
    # ---------------

    # We use the new ``multimeter`` device for recording from the model
    # neurons. At present, ``iaf_cond_alpha`` is one of few models
    # supporting ``multimeter`` recording.  Support for more models will
    # be added soon; until then, you need to use ``voltmeter`` to record
    # from other models.
    #
    # We configure multimeter to record membrane potential to membrane
    # potential at certain intervals to memory only. We record the GID of
    # the recorded neurons, but not the time.
    nest.CopyModel('multimeter', 'RecordingNode',
                   params = {'interval'   : Params['sim_interval'],
                             'record_from': ['V_m'],
                             'record_to'  : ['memory'],
                             'withgid'    : True,
                             'withtime'   : False})


    # Populations
    # ===========

    # We now create the neuron populations in the model, again in the
    # form of Python dictionaries. We define them in order from eye via
    # thalamus to cortex.
    #
    # We first define a dictionary defining common properties for all
    # populations
    layerProps = {'rows'     : Params['N'],
                  'columns'  : Params['N'],
                  'extent'   : [Params['visSize'], Params['visSize']],
                  'edge_wrap': True}
    # This dictionary does not yet specify the elements to put into the
    # layer, since they will differ from layer to layer. We will add them
    # below by updating the ``'elements'`` dictionary entry for each
    # population.

    # Retina
    # ------


    layerProps.update({'elements': 'RetinaNode'})
    retina = topo.CreateLayer(layerProps)

    # Original: Gabor retina input

    # GetPosition => [column, row] !
    # GetLeave => [2,3,4,5,6,7] =>
    #  2 4 6
    #  3 5 7

    retina_positions_x = np.array([topo.GetPosition([n])[0][1] for n in nest.GetLeaves(retina)[0]]).reshape(Params['N'], Params['N']).T
    retina_positions_y = np.array([topo.GetPosition([n])[0][0] for n in nest.GetLeaves(retina)[0]]).reshape(Params['N'], Params['N']).T

    if Params['bar']:
        if Params['lambda_dg'] >= 0:
            ## Now set phases of retinal oscillators; we use a list comprehension instead
            ## of a loop.
            tmp = [Params['retDC'] * barInit(topo.GetPosition([n])[0],Params["lambda_dg"], Params["phi_dg"],Params["visSize"])
                   for n in nest.GetLeaves(retina)[0]]
            screen = np.array(tmp).reshape(Params['N'], Params['N']).T
            plt.imshow(screen)
            plt.show()

            [nest.SetStatus([n], {"rate": Params['retDC'] * barInit(topo.GetPosition([n])[0], Params["lambda_dg"], Params["phi_dg"],Params["visSize"])})
                    for n in nest.GetLeaves(retina)[0]]

            #tmp = [phaseInit(topo.GetPosition([n])[0],Params["lambda_dg"], Params["phi_dg"])
            #       for n in nest.GetLeaves(retina)[0]]
            #screen = np.array(tmp).reshape(Params['N'], Params['N'])
            #plt.imshow(screen)
            #plt.colorbar()
            #plt.show()

            #[nest.SetStatus([n], {"phase": phaseInit(topo.GetPosition([n])[0],
            #                                     Params["lambda_dg"],
            #                                     Params["phi_dg"])})
            #for n in nest.GetLeaves(retina)[0]]
            #[nest.SetStatus([n], {"amplitude": 0.0 * barInit(topo.GetPosition([n])[0],
            #                                       Params["lambda_dg"],
            #                                       Params["phi_dg"],
            #                                       Params["visSize"])})
            #for n in nest.GetLeaves(retina)[0]]

            #[nest.SetStatus([n], {"phase": phaseBarInit(topo.GetPosition([n])[0],
            #                                       Params["lambda_dg"],
            #                                       Params["phi_dg"],
            #                                       Params["visSize"])})
            #for n in nest.GetLeaves(retina)[0]]
    else:
        if Params['lambda_dg'] >= 0:
            # Now set phases of retinal oscillators; we use a list comprehension instead
            # of a loop.
            [nest.SetStatus([n], {"phase": phaseInit(topo.GetPosition([n])[0],
                                                 Params["lambda_dg"],
                                                 Params["phi_dg"])})
            for n in nest.GetLeaves(retina)[0]]
        else:
            # Leonardo: Random retina input

            [nest.SetStatus([n], {"phase": phaseInit(topo.GetPosition([n])[0],
                                                 np.pi * np.random.rand(),
                                                 np.pi * np.random.rand())})
            for n in nest.GetLeaves(retina)[0]]

    # Thalamus
    # --------

    # We first introduce specific neuron models for the thalamic relay
    # cells and interneurons. These have identical properties, but by
    # treating them as different models, we can address them specifically
    # when building connections.
    #
    # We use a list comprehension to do the model copies.
    [nest.CopyModel('ThalamicNeuron', SpecificModel) for SpecificModel in ('TpRelay', 'TpInter')]

    # Now we can create the layer, with one relay cell and one
    # interneuron per location:
    layerProps.update({'elements': ['TpRelay', 'TpInter']})
    Tp = topo.CreateLayer(layerProps)

    # Reticular nucleus
    # -----------------
    # We follow the same approach as above, even though we have only a
    # single neuron in each location.
    [nest.CopyModel('ThalamicNeuron', SpecificModel) for SpecificModel in ('RpNeuron',)]
    layerProps.update({'elements': 'RpNeuron'})
    Rp = topo.CreateLayer(layerProps)

    # Primary visual cortex
    # ---------------------

    # We follow again the same approach. We differentiate neuron types
    # between layers and between pyramidal cells and interneurons. At
    # each location, there are two pyramidal cells and one interneuron in
    # each of layers 2-3, 4, and 5-6. Finally, we need to differentiate
    # between vertically and horizontally tuned populations. When creating
    # the populations, we create the vertically and the horizontally
    # tuned populations as separate populations.

    # We use list comprehesions to create all neuron types:
    [nest.CopyModel('CtxExNeuron', layer+'pyr') for layer in ('L23','L4','L56')]
    [nest.CopyModel('CtxInNeuron', layer+'in' ) for layer in ('L23','L4','L56')]

    # Now we can create the populations, suffixes h and v indicate tuning
    layerProps.update({'elements': ['L23pyr', 2, 'L23in', 1,
                                    'L4pyr' , 2, 'L4in' , 1,
                                    'L56pyr', 2, 'L56in', 1]})
    Vp_h = topo.CreateLayer(layerProps)
    Vp_v = topo.CreateLayer(layerProps)

    # Collect all populations
    # -----------------------

    # For reference purposes, e.g., printing, we collect all populations
    # in a tuple:
    populations = (retina, Tp, Rp, Vp_h, Vp_v)

    # Inspection
    # ----------

    # We can now look at the network using `PrintNetwork`:
    nest.PrintNetwork()

    # We can also try to plot a single layer in a network. For
    # simplicity, we use Rp, which has only a single neuron per position.
    #topo.PlotLayer(Rp)
    #pylab.title('Layer Rp')
    #pylab.show()


    # Synapse models
    # ==============

    # Actual synapse dynamics, e.g., properties such as the synaptic time
    # course, time constants, reversal potentials, are properties of
    # neuron models in NEST and we set them in section `Neuron models`_
    # above. When we refer to *synapse models* in NEST, we actually mean
    # connectors which store information about connection weights and
    # delays, as well as port numbers at the target neuron (``rport``)
    # and implement synaptic plasticity. The latter two aspects are not
    # relevant here.
    #
    # We just use NEST's ``static_synapse`` connector but copy it to
    # synapse models ``AMPA`` and ``GABA_A`` for the sake of
    # explicitness. Weights and delays are set as needed in section
    # `Connections`_ below, as they are different from projection to
    # projection. De facto, the sign of the synaptic weight decides
    # whether input via a connection is handle by the ``_ex`` or the
    # ``_in`` synapse.
    nest.CopyModel('static_synapse', 'AMPA')
    nest.CopyModel('static_synapse', 'GABA_A')

    # Connections
    # ====================

    # Building connections is the most complex part of network
    # construction. Connections are specified in Table 1 in the
    # Hill-Tononi paper. As pointed out above, we only consider AMPA and
    # GABA_A synapses here.  Adding other synapses is tedious work, but
    # should pose no new principal challenges. We also use a uniform in
    # stead of a Gaussian distribution for the weights.
    #
    # The model has two identical primary visual cortex populations,
    # ``Vp_v`` and ``Vp_h``, tuned to vertical and horizonal gratings,
    # respectively. The *only* difference in the connection patterns
    # between the two populations is the thalamocortical input to layers
    # L4 and L5-6 is from a population of 8x2 and 2x8 grid locations,
    # respectively. Furthermore, inhibitory connection in cortex go to
    # the opposing orientation population as to the own.
    #
    # To save us a lot of code doubling, we thus defined properties
    # dictionaries for all connections first and then use this to connect
    # both populations. We follow the subdivision of connections as in
    # the Hill & Tononi paper.
    #
    # **Note:** Hill & Tononi state that their model spans 8 degrees of
    # visual angle and stimuli are specified according to this. On the
    # other hand, all connection patterns are defined in terms of cell
    # grid positions. Since the NEST Topology Module defines connection
    # patterns in terms of the extent given in degrees, we need to apply
    # the following scaling factor to all lengths in connections:
    dpc = Params['visSize'] / (Params['N'] - 1)

    # We will collect all same-orientation cortico-cortical connections in
    ccConnections = []
    # the cross-orientation cortico-cortical connections in
    ccxConnections = []
    # and all cortico-thalamic connections in
    ctConnections = []

    # spread of connections
    horIntraCtxExRadius = 12.0
    horIntraCtxExRadiusL4 = 7.0
    verIntraCtxExRadius = 2.0
    allIntraInhRadius = 7.0

    allCtxTalRadius = 5.0
    horCtxTalRecBig = 4.0
    horCtxTalRecSmall = 1.0
    difTalCtxRadius = 5.0

    # Horizontal intralaminar
    # -----------------------
    # *Note:* "Horizontal" means "within the same cortical layer" in this
    # case.
    #
    # We first define a dictionary with the (most) common properties for
    # horizontal intralaminar connection. We then create copies in which
    # we adapt those values that need adapting, and
    horIntraBase = {"connection_type": "divergent",
                    "synapse_model": "AMPA",
                    "mask": {"circular": {"radius": horIntraCtxExRadius * dpc}},
                    "kernel": {"gaussian": {"p_center": 0.05, "sigma": 7.5 * dpc}},
                    "weights": 1.0,
                    "delays": {"uniform": {"min": 1.75, "max": 2.25}}}

    # We use a loop to do the for for us. The loop runs over a list of
    # dictionaries with all values that need updating
    for conn in [{"sources": {"model": "L23pyr"}, "targets": {"model": "L23pyr"}},
                 {"sources": {"model": "L23pyr"}, "targets": {"model": "L23in" }},
                 {"sources": {"model": "L4pyr" }, "targets": {"model": "L4pyr" },
                  "mask"   : {"circular": {"radius": horIntraCtxExRadiusL4 * dpc}}},
                 {"sources": {"model": "L4pyr" }, "targets": {"model": "L4in"  },
                  "mask"   : {"circular": {"radius": horIntraCtxExRadiusL4 * dpc}}},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "L56pyr" }},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "L56in"  }}]:
        ndict = horIntraBase.copy()
        ndict.update(conn)
        ccConnections.append(ndict)

    # Vertical intralaminar
    # -----------------------
    # *Note:* "Vertical" means "between cortical layers" in this
    # case.
    #
    # We proceed as above.
    verIntraBase = {"connection_type": "divergent",
                    "synapse_model": "AMPA",
                    "mask": {"circular": {"radius": verIntraCtxExRadius * dpc}},
                    "kernel": {"gaussian": {"p_center": 1.0, "sigma": 7.5 * dpc}},
                    "weights": 2.0,
                    "delays": {"uniform": {"min": 1.75, "max": 2.25}}}

    for conn in [{"sources": {"model": "L23pyr"}, "targets": {"model": "L56pyr"}, "weights": 1.0},
                 {"sources": {"model": "L23pyr"}, "targets": {"model": "L56in"}, "weights": 1.0}, #---keiko
                 #{"sources": {"model": "L23pyr"}, "targets": {"model": "L23in"}, "weights": 1.0}, #---original
                 {"sources": {"model": "L4pyr" }, "targets": {"model": "L23pyr"}},
                 {"sources": {"model": "L4pyr" }, "targets": {"model": "L23in" }},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "L23pyr"}},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "L23in" }},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "L4pyr" }},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "L4in"  }}]:
        ndict = verIntraBase.copy()
        ndict.update(conn)
        ccConnections.append(ndict)

    # Intracortical inhibitory
    # ------------------------
    #
    # We proceed as above, with the following difference: each connection
    # is added to the same-orientation and the cross-orientation list of
    # connections.
    #
    # **Note:** Weights increased from -1.0 to -2.0, to make up for missing GabaB
    #
    # Note that we have to specify the **weight with negative sign** to make
    # the connections inhibitory.
    intraInhBase = {"connection_type": "divergent",
                    "synapse_model": "GABA_A",
                    "mask": {"circular": {"radius": allIntraInhRadius * dpc}},
                    "kernel": {"gaussian": {"p_center": 0.25, "sigma": 7.5 * dpc}},
                    "weights": -2.0,
                    "delays": {"uniform": {"min": 1.75, "max": 2.25}}}

    # We use a loop to do the for for us. The loop runs over a list of
    # dictionaries with all values that need updating
    for conn in [{"sources": {"model": "L23in"}, "targets": {"model": "L23pyr"}},
                 {"sources": {"model": "L23in"}, "targets": {"model": "L23in" }},
                 {"sources": {"model": "L4in" }, "targets": {"model": "L4pyr" }},
                 {"sources": {"model": "L4in" }, "targets": {"model": "L4in"  }},
                 {"sources": {"model": "L56in"}, "targets": {"model": "L56pyr"}},
                 {"sources": {"model": "L56in"}, "targets": {"model": "L56in" }}]:
        ndict = intraInhBase.copy()
        ndict.update(conn)
        ccConnections.append(ndict)
        ccxConnections.append(ndict)

    # Corticothalamic
    # ---------------
    corThalBase = {"connection_type": "divergent",
                   "synapse_model": "AMPA",
                   "mask": {"circular": {"radius": allCtxTalRadius * dpc}},
                   "kernel": {"gaussian": {"p_center": 0.5, "sigma": 7.5 * dpc}},
                   "weights": 1.0,
                   "delays": {"uniform": {"min": 7.5, "max": 8.5}}}

    # We use a loop to do the for for us. The loop runs over a list of
    # dictionaries with all values that need updating
    for conn in [{"sources": {"model": "L56pyr"}, "targets": {"model": "TpRelay" }},
                 {"sources": {"model": "L56pyr"}, "targets": {"model": "TpInter" }}]:
        #ndict = corThalBase.copy() #---keiko
        ndict = intraInhBase.copy() #---original
        ndict.update(conn)
        ctConnections.append(ndict)

    # Corticoreticular
    # ----------------

    # In this case, there is only a single connection, so we write the
    # dictionary itself; it is very similar to the corThalBase, and to
    # show that, we copy first, then update. We need no ``targets`` entry,
    # since Rp has only one neuron per location.
    corRet = corThalBase.copy()
    corRet.update({"sources": {"model": "L56pyr"}, "weights": 2.5})


    # Build all connections beginning in cortex
    # -----------------------------------------

    # Cortico-cortical, same orientation
    print("Connecting: cortico-cortical, same orientation")
    [topo.ConnectLayers(Vp_h, Vp_h, conn) for conn in ccConnections]
    [topo.ConnectLayers(Vp_v, Vp_v, conn) for conn in ccConnections]

    # Cortico-cortical, cross-orientation
    print("Connecting: cortico-cortical, other orientation")
    [topo.ConnectLayers(Vp_h, Vp_v, conn) for conn in ccxConnections]
    [topo.ConnectLayers(Vp_v, Vp_h, conn) for conn in ccxConnections]

    # Cortico-thalamic connections
    print("Connecting: cortico-thalamic")
    [topo.ConnectLayers(Vp_h, Tp, conn) for conn in ctConnections]
    [topo.ConnectLayers(Vp_v, Tp, conn) for conn in ctConnections]
    topo.ConnectLayers(Vp_h, Rp, corRet)
    topo.ConnectLayers(Vp_v, Rp, corRet)

    # Thalamo-cortical connections
    # ----------------------------

    # **Note:** According to the text on p. 1674, bottom right, of
    # the Hill & Tononi paper, thalamocortical connections are
    # created by selecting from the thalamic population for each
    # L4 pyramidal cell, ie, are *convergent* connections.
    #
    # We first handle the rectangular thalamocortical connections.
    thalCorRect = {"connection_type": "convergent",
                   "sources": {"model": "TpRelay"},
                   "synapse_model": "AMPA",
                   "weights": 5.0,
                   "delays": {"uniform": {"min": 2.75, "max": 3.25}}}

    print("Connecting: thalamo-cortical")

    # Horizontally tuned
    thalCorRect.update({"mask": {"rectangular": {"lower_left" : [-horCtxTalRecBig*dpc, -horCtxTalRecSmall*dpc],
                                                 "upper_right": [ horCtxTalRecBig*dpc,  horCtxTalRecSmall*dpc]}}})
    for conn in [{"targets": {"model": "L4pyr" }, "kernel": 0.5},
                 {"targets": {"model": "L56pyr"}, "kernel": 0.3}]:
        thalCorRect.update(conn)
        topo.ConnectLayers(Tp, Vp_h, thalCorRect)

    # Vertically tuned
    thalCorRect.update({"mask": {"rectangular": {"lower_left" : [-horCtxTalRecSmall*dpc, -horCtxTalRecBig*dpc],
                                                 "upper_right": [ horCtxTalRecSmall*dpc,  horCtxTalRecBig*dpc]}}})
    for conn in [{"targets": {"model": "L4pyr" }, "kernel": 0.5},
                 {"targets": {"model": "L56pyr"}, "kernel": 0.3}]:
        thalCorRect.update(conn)
        topo.ConnectLayers(Tp, Vp_v, thalCorRect)

    # Diffuse connections
    thalCorDiff = {"connection_type": "convergent",
                   "sources": {"model": "TpRelay"},
                   "synapse_model": "AMPA",
                   "weights": 5.0,
                   "mask": {"circular": {"radius": difTalCtxRadius * dpc}},
                   "kernel": {"gaussian": {"p_center": 0.1, "sigma": 7.5 * dpc}},
                   "delays": {"uniform": {"min": 2.75, "max": 3.25}}}

    #---keiko
    for conn in [{"targets": {"model": "L4in"}},
                 {"targets": {"model": "L56in"}}]:
    #---original
    #for conn in [{"targets": {"model": "L4pyr"}},
    #             {"targets": {"model": "L56pyr"}}]:
        thalCorDiff.update(conn)
        topo.ConnectLayers(Tp, Vp_h, thalCorDiff)
        topo.ConnectLayers(Tp, Vp_v, thalCorDiff)

    # Thalamic connections
    # --------------------

    # Connections inside thalamus, including Rp
    #
    # *Note:* In Hill & Tononi, the inhibition between Rp cells is mediated by
    # GABA_B receptors. We use GABA_A receptors here to provide some self-dampening
    # of Rp.
    #
    # **Note:** The following code had a serious bug in v. 0.1: During the first
    # iteration of the loop, "synapse_model" and "weights" were set to "AMPA" and "0.1",
    # respectively and remained unchanged, so that all connections were created as
    # excitatory connections, even though they should have been inhibitory. We now
    # specify synapse_model and weight explicitly for each connection to avoid this.

    thalBase = {"connection_type": "divergent",
                "delays": {"uniform": {"min": 1.75, "max": 2.25}}}

    print("Connecting: intra-thalamic")

    allTalRetEx = 2.0
    allTalTalIntRelInh = 2.0
    allTalTalIntIntInh = 2.0

    allRpTp = 12.0

    for src, tgt, conn in [(Tp, Rp, {"sources": {"model": "TpRelay"},
                                     "synapse_model": "AMPA",
                                     "mask": {"circular": {"radius": allTalRetEx * dpc}},
                                     "kernel": {"gaussian": {"p_center": 1.0, "sigma": 7.5 * dpc}},
                                     "weights": 2.0}),
                           (Tp, Tp, {"sources": {"model": "TpInter"},
                                     "targets": {"model": "TpRelay"},
                                     "synapse_model": "GABA_A",
                                     "weights": -1.0,
                                     "mask": {"circular": {"radius":  allTalTalIntRelInh * dpc}},
                                     "kernel": {"gaussian": {"p_center": 0.25, "sigma": 7.5 * dpc}}}),
                           (Tp, Tp, {"sources": {"model": "TpInter"},
                                     "targets": {"model": "TpInter"},
                                     "synapse_model": "GABA_A",
                                     "weights": -1.0,
                                     "mask": {"circular": {"radius": allTalTalIntIntInh * dpc}},
                                     "kernel": {"gaussian": {"p_center": 0.25, "sigma": 7.5 * dpc}}}),
                           (Rp, Tp, {"targets": {"model": "TpRelay"},
                                     "synapse_model": "GABA_A",
                                     "weights": -1.0,
                                     "mask": {"circular": {"radius": allRpTp * dpc}},
                                     "kernel": {"gaussian": {"p_center": 0.15, "sigma": 7.5 * dpc}}}),
                           (Rp, Tp, {"targets": {"model": "TpInter"},
                                     "synapse_model": "GABA_A",
                                     "weights": -1.0,
                                     "mask": {"circular": {"radius": allRpTp * dpc}},
                                     "kernel": {"gaussian": {"p_center": 0.15, "sigma": 7.5 * dpc}}}),
                           (Rp, Rp, {"targets": {"model": "RpNeuron"},
                                     "synapse_model": "GABA_A",
                                     "weights": -1.0,
                                     "mask": {"circular": {"radius": allRpTp * dpc}},
                                     "kernel": {"gaussian": {"p_center": 0.5, "sigma": 7.5 * dpc}}})]:
        thalBase.update(conn)
        topo.ConnectLayers(src, tgt, thalBase)

    # Thalamic input
    # --------------

    # Input to the thalamus from the retina.
    #
    # **Note:** Hill & Tononi specify a delay of 0 ms for this connection.
    # We use 1 ms here.
    retThal = {"connection_type": "divergent",
               "synapse_model": "AMPA",
               "mask": {"circular": {"radius": 1.0 * dpc}},
               "kernel": {"gaussian": {"p_center": 0.75, "sigma": 2.5 * dpc}},
               "weights": 10.0,
               "delays": 1.0}

    print("Connecting: retino-thalamic")

    for conn in [{"targets": {"model": "TpRelay"}},
                 {"targets": {"model": "TpInter"}}]:
        retThal.update(conn)
        topo.ConnectLayers(retina, Tp, retThal)


    # Checks on connections
    # ---------------------

    # As a very simple check on the connections created, we inspect
    # the connections from the central node of various layers.

    # Connections from Retina to TpRelay
    #topo.PlotTargets(topo.FindCenterElement(retina), Tp, 'TpRelay', 'AMPA')
    #pylab.title('Connections Retina -> TpRelay')
    #pylab.show()

    # Connections from TpRelay to L4pyr in Vp (horizontally tuned)
    #topo.PlotTargets(topo.FindCenterElement(Tp), Vp_h, 'L4pyr', 'AMPA')
    #pylab.title('Connections TpRelay -> Vp(h) L4pyr')
    #pylab.show()

    # Connections from TpRelay to L4pyr in Vp (vertically tuned)
    #topo.PlotTargets(topo.FindCenterElement(Tp), Vp_v, 'L4pyr', 'AMPA')
    #pylab.title('Connections TpRelay -> Vp(v) L4pyr')
    #pylab.show()



    # [n for n in nest.GetLeaves(retina)[0]]

    tini = time.time()

    for isim in range(1,Nsim+1,1):

        nest.ResetNetwork()
        msd = int(round(time.time() * 1000))
        nest.SetKernelStatus({'grng_seed' : msd})
        nest.SetKernelStatus({'rng_seeds' : range(msd+N_vp+1, msd+2*N_vp+1)})
        nest.SetStatus([0], {"time": 0.0})

        # Set folder name from parameter names
        folder_name = "sim_%d_" % isim
        for p,v in Params.items():
            folder_name += p + "_%.2f_" % v

        # Example simulation
        # ====================

        # This simulation is set up to create a step-wise visualization of
        # the membrane potential. To do so, we simulate ``sim_interval``
        # milliseconds at a time, then read out data from the multimeters,
        # clear data from the multimeters and plot the data as pseudocolor
        # plots.

        # show time during simulation
        nest.SetStatus([0],{'print_time': True})

        # first run the network for a while until it stabilises
        if Params['discardtime'] > 0.0:
            nest.Simulate(Params['discardtime'])

        # then connect the recorders
        print("Connecting: Recording devices")
        recorders = {}
        for name, loc, population, model in [('TpRelay', 1, Tp, 'TpRelay'),
                                             ('Rp', 2, Rp, 'RpNeuron'),
                                             ('Vp_v L4pyr', 3, Vp_v, 'L4pyr'),
                                             ('Vp_h L4pyr', 4, Vp_h, 'L4pyr')]:
            recorders[name] = (nest.Create('RecordingNode'), loc)
            tgts = [nd for nd in nest.GetLeaves(population)[0]
                    if nest.GetStatus([nd], 'model')[0] == model]
            nest.Connect(recorders[name][0], tgts)  # one recorder to all targets

        # create the spike detectors
        detectors = {}
        for name, population, model in [('Vp_v L4pyr', Vp_v, 'L4pyr'),
                                        ('Vp_h L4pyr', Vp_h, 'L4pyr'),
                                        ('Vp_v L23pyr', Vp_v, 'L23pyr'),
                                        ('Vp_h L23pyr', Vp_h, 'L23pyr'),
                                        ('Vp_v L56pyr', Vp_v, 'L56pyr'),
                                        ('Vp_h L56pyr', Vp_h, 'L56pyr'),
                                        ('TpRelay', Tp, 'TpRelay'),
                                        ('Rp', Rp, 'RpNeuron'),
                                        ('Retina', retina, 'RetinaNode')]:
            tgts = [nd for nd in nest.GetLeaves(population)[0]
                    if nest.GetStatus([nd], 'model')[0] == model]
            detectors[name] = (nest.Create('spike_detector', params={"withgid": True, "withtime": True}), loc)
            print(name + ' : %d' % max(tgts))

            # TODO For some reason, some retina nodes have different synapse types. Connect one by one for now
            if name == 'Retina':
                for t in tgts:
                    try:
                        nest.Connect([t], detectors[name][0])
                    except:
                        print('%d did not work' % t)
            else:
                nest.Connect(tgts, detectors[name][0])

        #TODO check difference with respect to running each step at a time
        # and run the simulation
        nest.Simulate(Params['simtime'])

        data_folder = './data/' + folder_name + 'detectors'
        if not os.path.isdir(data_folder):
            os.makedirs(data_folder)

        figure_folder = './figures/' + folder_name + 'detectors'
        if not os.path.isdir(figure_folder):
            os.makedirs(figure_folder)

        # save raster plot and data
        for name, r in detectors.items():
            rec = r[0]

            # raster plot
            p = raster_plot.from_device(rec, hist=True)

            pylab.title(name)
            f = p[0].figure
            f.set_size_inches(15, 9)
            f.savefig(figure_folder + '/spikes_' + name + '.png', dpi=100)

            # data
            spikes = nest.GetStatus(rec, "events")[0]
            print ('Saving spikes for ' + name)
            with open(data_folder + '/spikes_' + name + '.pickle', 'w') as f:
                pickle.dump(spikes, f)

            scipy.io.savemat(data_folder + '/spikes_' + name + '.mat', mdict={'senders': spikes['senders'], 'times': spikes['times']})

            # clear data from spike detectors
            nest.SetStatus(rec, {'n_events': 0})

            plt.close()

        # just for some information at the end
        #print(nest.GetKernelStatus())

        #--- keiko 9/14/2016
        #f_Vm.close()
        #f_idx.close()

    elapsed = time.time() - tini
    m = "Elapsed time : %.2f s" % elapsed
    print(m)