{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# NEST spatial example: A case-based tutorial\n\n:Author: Hans Ekkehard Plesser\n:Institution: Norwegian University of Life Sciences\n:Version: 0.4\n:Date: 21 November 2012\n:Copyright: The NEST Initiative (2004)\n:License: Creative Commons Attribution License\n\n**NOTE:** The network generated by this script does generate\ndynamics in which the activity of the entire system, especially\nRp and Vp oscillates with approx 5 Hz. This is different from\nthe full model. Deviations are due to the different model type\nand the elimination of a number of connections, with no changes\nto the weights.\n\n## Introduction\n\nThis tutorial shows you how to implement a simplified version of the\nHill-Tononi model of the early visual pathway using NEST. The model\nis described in the paper\n\n S. L. Hill and G. Tononi.\n Modeling Sleep and Wakefulness in the Thalamocortical System.\n J Neurophysiology **93**:1671-1698 (2005).\n Freely available via `doi 10.1152/jn.00915.2004\n `_.\n\nWe simplify the model somewhat both to keep this tutorial a bit\nshorter, and because some details of the Hill-Tononi model are not\ncurrently supported by NEST. Simplifications include:\n\n1. We use the ``iaf_cond_alpha`` neuron model, which is\n simpler than the Hill-Tononi model.\n\n#. As the ``iaf_cond_alpha`` neuron model only supports two\n synapses (labeled \"ex\" and \"in\"), we only include AMPA and\n GABA_A synapses.\n\n#. We ignore the secondary pathway (Ts, Rs, Vs), since it adds just\n more of the same from a technical point of view.\n\n#. Synaptic delays follow a Gaussian distribution in the HT\n model. This implies actually a Gaussian distributions clipped at\n some small, non-zero delay, since delays must be\n positive. Currently, there is a bug in the module when using clipped\n Gaussian distribution. We therefore draw delays from a\n uniform distribution.\n\n#. Some further adaptations are given at the appropriate locations in\n the script.\n\nThis tutorial is divided in the following sections:\n\nPhilosophy_\n Discusses the philosophy applied to model implementation in this\n tutorial\n\nPreparations_\n Neccessary steps to use NEST\n\n`Configurable Parameters`_\n Define adjustable network parameters\n\n`Neuron Models`_\n Define the neuron models needed by the network model\n\nPopulations_\n Create Populations\n\n`Synapse models`_\n Define the synapse models used in the network model\n\nConnections_\n Create Connections\n\n`Example simulation`_\n Perform a small simulation for illustration. This\n section also discusses the setup for recording.\n\n## Philosophy\n\nA network models has two essential components: *populations* and\n*projections*. We first use NEST's ``CopyModel()`` mechanism to\ncreate specific models for all populations and subpopulations in\nthe network, and then create the populations using the\n``Create()`` function.\n\nWe use a two-stage process to create the connections, mainly\nbecause the same configurations are required for a number of\nprojections: we first define dictionaries specifying the\nconnections, then apply these dictionaries later.\n\nThe way in which we declare the network model here is an\nexample. You should not consider it the last word: we expect to see\na significant development in strategies and tools for network\ndescriptions in the future. The following contributions to CNS*09\nseem particularly interesting\n\n- Ralf Ansorg & Lars Schwabe. Declarative model description and\n code generation for hybrid individual- and population-based\n simulations of the early visual system (P57);\n- Sharon Crook, R. Angus Silver, & Padraig Gleeson. Describing\n and exchanging models of neurons and neuronal networks with\n NeuroML (F1);\n\nas well as the following paper which will apply in PLoS\nComputational Biology shortly:\n\n- Eilen Nordlie, Marc-Oliver Gewaltig, & Hans Ekkehard Plesser.\n Towards reproducible descriptions of neuronal network models.\n\n## Preparations\n\nPlease make sure that your ``PYTHONPATH`` is set correctly, so\nthat Python can find the NEST Python module.\n\n**Note:** By default, the script does not show any graphics.\nSet ``SHOW_FIGURES`` to ``True`` to activate graphics.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pprint import pprint\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nSHOW_FIGURES = False\n\nif SHOW_FIGURES:\n plt.ion()\nelse:\n plt_show = plt.show\n\n def nop(s=None, block=None):\n pass\n\n plt.show = nop" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This tutorial gives a brief introduction to the ConnPlotter\ntoolbox. It is by no means complete.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Load pynest\nimport nest\n\n# Make sure we start with a clean slate, even if we re-run the script\n# in the same Python session.\nnest.ResetKernel()\n\n# Import math, we need Pi\nimport math" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configurable Parameters\n\nHere we define those parameters that we take to be\nconfigurable. The choice of configurable parameters is obviously\narbitrary, and in practice one would have far more configurable\nparameters. We restrict ourselves to:\n\n- Network size in neurons ``N``, each layer is ``N x N``.\n- Network size in subtended visual angle ``visSize``, in degree.\n- Temporal frequency of drifting grating input ``f_dg``, in Hz.\n- Spatial wavelength and direction of drifting grating input,\n ``lambda_dg`` and ``phi_dg``, in degree/radian.\n- Background firing rate of retinal nodes and modulation amplitude,\n ``retDC`` and ``retAC``, in Hz.\n- Simulation duration ``simtime``; actual simulation is split into\n intervals of ``sim_interval`` length, so that the network state\n can be visualized in those intervals. Times are in ms.\n- Periodic boundary conditions, ``edge_wrap``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Params = {'N': 40,\n 'visSize': 8.0,\n 'f_dg': 2.0,\n 'lambda_dg': 2.0,\n 'phi_dg': 0.0,\n 'retDC': 30.0,\n 'retAC': 30.0,\n 'simtime': 100.0,\n 'sim_interval': 1.0,\n 'edge_wrap': True\n }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Neuron Models\n\nWe declare models in two steps:\n\n1. We define a dictionary specifying the NEST neuron model to use\n as well as the parameters for that model.\n#. We create three copies of this dictionary with parameters\n adjusted to the three model variants specified in Table~2 of\n Hill & Tononi (2005) (cortical excitatory, cortical inhibitory,\n thalamic)\n\nIn addition, we declare the models for the stimulation and\nrecording devices.\n\n## The general neuron model\n\nWe use the ``iaf_cond_alpha`` neuron, which is an\nintegrate-and-fire neuron with two conductance-based synapses which\nhave alpha-function time course. Any input with positive weights\nwill automatically directed to the synapse labeled ``_ex``, any\nwith negative weights to the synapes labeled ``_in``. We define\n**all** parameters explicitly here, so that no information is\nhidden in the model definition in NEST. ``V_m`` is the membrane\npotential to which the model neurons will be initialized.\nThe model equations and parameters for the Hill-Tononi neuron model\nare given on pp. 1677f and Tables 2 and 3 in that paper. Note some\npeculiarities and adjustments:\n\n- Hill & Tononi specify their model in terms of the membrane time\n constant, while the ``iaf_cond_alpha`` model is based on the\n membrane capcitance. Interestingly, conductances are unitless in\n the H&T model. We thus can use the time constant directly as\n membrane capacitance.\n- The model includes sodium and potassium leak conductances. We\n combine these into a single one as follows:\n\n\\begin{align}-g_{NaL}(V-E_{Na}) - g_{KL}(V-E_K)\n = -(g_{NaL}+g_{KL})\n \\left(V-\\frac{g_{NaL}E_{NaL}+g_{KL}E_K}{g_{NaL}g_{KL}}\\right)\\end{align}\n\n- We write the resulting expressions for g_L and E_L explicitly\n below, to avoid errors in copying from our pocket calculator.\n- The paper gives a range of 1.0-1.85 for g_{KL}, we choose 1.5\n here.\n- The Hill-Tononi model has no explicit reset or refractory\n time. We arbitrarily set V_reset and t_ref.\n- The paper uses double exponential time courses for the synaptic\n conductances, with separate time constants for the rising and\n fallings flanks. Alpha functions have only a single time\n constant: we use twice the rising time constant given by Hill and\n Tononi.\n- In the general model below, we use the values for the cortical\n excitatory cells as defaults. Values will then be adapted below.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.CopyModel('iaf_cond_alpha', 'NeuronModel',\n params={'C_m': 16.0,\n 'E_L': (0.2 * 30.0 + 1.5 * -90.0) / (0.2 + 1.5),\n 'g_L': 0.2 + 1.5,\n 'E_ex': 0.0,\n 'E_in': -70.0,\n 'V_reset': -60.0,\n 'V_th': -51.0,\n 't_ref': 2.0,\n 'tau_syn_ex': 1.0,\n 'tau_syn_in': 2.0,\n 'I_e': 0.0,\n 'V_m': -70.0})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adaptation of models for different populations\n\nWe must copy the `NeuronModel` dictionary explicitly, otherwise\nPython would just create a reference.\n\n### Cortical excitatory cells\nParameters are the same as above, so we need not adapt anything\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.CopyModel('NeuronModel', 'CtxExNeuron')\n\n# Cortical inhibitory cells\n# .........................\nnest.CopyModel('NeuronModel', 'CtxInNeuron',\n params={'C_m': 8.0,\n 'V_th': -53.0,\n 't_ref': 1.0})\n\n# Thalamic cells\n# ..............\nnest.CopyModel('NeuronModel', 'ThalamicNeuron',\n params={'C_m': 8.0,\n 'V_th': -53.0,\n 't_ref': 1.0,\n 'E_in': -80.0})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Input generating nodes\n\nInput is generated by sinusoidally modulate Poisson generators,\norganized in a square layer of retina nodes. These nodes require a\nslightly more complicated initialization than all other elements of\nthe network:\n\n- Average firing rate ``rate``, firing rate modulation depth ``amplitude``,\n and temporal modulation frequency ``frequency`` are the same for all\n retinal nodes and are set directly below.\n- The temporal phase ``phase`` of each node depends on its position in\n the grating and can only be assigned after the retinal layer has\n been created.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.CopyModel('sinusoidal_poisson_generator', 'RetinaNode',\n params={'amplitude': Params['retAC'],\n 'rate': Params['retDC'],\n 'frequency': Params['f_dg'],\n 'phase': 0.0,\n 'individual_spike_trains': False})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Recording nodes\n\nWe use the ``multimeter`` device for recording from the model\nneurons. At present, ``iaf_cond_alpha`` is one of few models\nsupporting ``multimeter`` recording. Support for more models will\nbe added soon; until then, you need to use ``voltmeter`` to record\nfrom other models.\n\nWe configure multimeter to record membrane potential to membrane\npotential at certain intervals to memory only. We record the node ID of\nthe recorded neurons, but not the time.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.CopyModel('multimeter', 'RecordingNode',\n params={'interval': Params['sim_interval'],\n 'record_from': ['V_m'],\n 'record_to': 'memory'})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Populations\n\nWe now create the neuron populations in the model. We define\nthem in order from eye via thalamus to cortex.\n\nWe first define a spatial grid defining common positions and\nparameters for all populations\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "layerGrid = nest.spatial.grid(shape=[Params['N'], Params['N']],\n extent=[Params['visSize'], Params['visSize']],\n edge_wrap=Params['edge_wrap'])\n# We can pass this object to the ``positions`` argument in ``Create``\n# to define the positions of the neurons." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Retina\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "retina = nest.Create('RetinaNode', positions=layerGrid)\n\n# Now set phases of retinal oscillators; we create a Parameter\n# which represents the phase based on the spatial properties of\n# the neuron.\n\nretina_phase = 360.0 / Params['lambda_dg'] * (math.cos(Params['phi_dg']) * nest.spatial.pos.x +\n math.sin(Params['phi_dg']) * nest.spatial.pos.y)\nretina.phase = retina_phase" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thalamus\n\nWe first introduce specific neuron models for the thalamic relay\ncells and interneurons. These have identical properties, but by\ntreating them as different populations, we can address them specifically\nwhen building connections.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for model_name in ('TpRelay', 'TpInter'):\n nest.CopyModel('ThalamicNeuron', model_name)\n\n# Now we can create the layers, one with relay cells,\n# and one with interneurons:\nTpRelay = nest.Create('TpRelay', positions=layerGrid)\nTpInter = nest.Create('TpInter', positions=layerGrid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reticular nucleus\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.CopyModel('ThalamicNeuron', 'RpNeuron')\nRp = nest.Create('RpNeuron', positions=layerGrid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Primary visual cortex\n\nWe follow again the same approach as with Thalamus. We differentiate\nneuron types between layers and between pyramidal cells and\ninterneurons. We have two layers for pyramidal cells, and two layers for\ninterneurons for each of layers 2-3, 4, and 5-6. Finally, we need to\ndifferentiate between vertically and horizontally tuned populations.\nWhen creating the populations, we create the vertically and the\nhorizontally tuned populations as separate dictionaries holding the\nlayers.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for layer in ('L23', 'L4', 'L56'):\n nest.CopyModel('CtxExNeuron', layer + 'pyr')\nfor layer in ('L23', 'L4', 'L56'):\n nest.CopyModel('CtxInNeuron', layer + 'in')\n\nname_dict = {'L23pyr': 2, 'L23in': 1,\n 'L4pyr': 2, 'L4in': 1,\n 'L56pyr': 2, 'L56in': 1}\n\n# Now we can create the populations, suffixes h and v indicate tuning\nVp_h_layers = {}\nVp_v_layers = {}\nfor layer_name, num_layers in name_dict.items():\n for i in range(num_layers):\n Vp_h_layers['{}_{}'.format(layer_name, i)] = nest.Create(layer_name, positions=layerGrid)\n Vp_v_layers['{}_{}'.format(layer_name, i)] = nest.Create(layer_name, positions=layerGrid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Collect all populations\n\nFor reference purposes, e.g., printing, we collect all populations\nin a tuple:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "populations = (retina, TpRelay, TpInter, Rp) + tuple(Vp_h_layers.values()) + tuple(Vp_v_layers.values())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inspection\n\nWe can now look at the network using `PrintNodes`:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.PrintNodes()\n\n# We can also try to plot a single layer in a network. All layers have\n# equal positions of the nodes.\nnest.PlotLayer(Rp)\nplt.title('Layer Rp')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Synapse models\n\nActual synapse dynamics, e.g., properties such as the synaptic time\ncourse, time constants, reversal potentials, are properties of\nneuron models in NEST and we set them in section `Neuron models`_\nabove. When we refer to *synapse models* in NEST, we actually mean\nconnectors which store information about connection weights and\ndelays, as well as port numbers at the target neuron (``rport``)\nand implement synaptic plasticity. The latter two aspects are not\nrelevant here.\n\nWe just use NEST's ``static_synapse`` connector but copy it to\nsynapse models ``AMPA`` and ``GABA_A`` for the sake of\nexplicitness. Weights and delays are set as needed in section\n`Connections`_ below, as they are different from projection to\nprojection. De facto, the sign of the synaptic weight decides\nwhether input via a connection is handle by the ``_ex`` or the\n``_in`` synapse.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.CopyModel('static_synapse', 'AMPA')\nnest.CopyModel('static_synapse', 'GABA_A')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Connections\n\nBuilding connections is the most complex part of network\nconstruction. Connections are specified in Table 1 in the\nHill-Tononi paper. As pointed out above, we only consider AMPA and\nGABA_A synapses here. Adding other synapses is tedious work, but\nshould pose no new principal challenges. We also use a uniform in\nstead of a Gaussian distribution for the weights.\n\nThe model has two identical primary visual cortex populations,\n``Vp_v`` and ``Vp_h``, tuned to vertical and horizonal gratings,\nrespectively. The *only* difference in the connection patterns\nbetween the two populations is the thalamocortical input to layers\nL4 and L5-6 is from a population of 8x2 and 2x8 grid locations,\nrespectively. Furthermore, inhibitory connection in cortex go to\nthe opposing orientation population as to the own.\n\nTo save us a lot of code doubling, we thus defined properties\ndictionaries for all connections first and then use this to connect\nboth populations. We follow the subdivision of connections as in\nthe Hill & Tononi paper.\n\nTODO: Rewrite this note.\n**Note:** Hill & Tononi state that their model spans 8 degrees of\nvisual angle and stimuli are specified according to this. On the\nother hand, all connection patterns are defined in terms of cell\ngrid positions. Since the NEST defines connection patterns in terms\nof the extent given in degrees, we need to apply the following\nscaling factor to all lengths in connections:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dpc = Params['visSize'] / (Params['N'] - 1)\n\n# We will collect all same-orientation cortico-cortical connections in\nccConnections = []\n# the cross-orientation cortico-cortical connections in\nccxConnections = []\n# and all cortico-thalamic connections in\nctConnections = []" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Horizontal intralaminar\n\n*Note:* \"Horizontal\" means \"within the same cortical layer\" in this\ncase.\n\nWe first define dictionaries with the (most) common properties for\nhorizontal intralaminar connection. We then create copies in which\nwe adapt those values that need adapting, and\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "horIntra_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"p\": 0.05*nest.spatial_distributions.gaussian(nest.spatial.distance, std=7.5 * dpc)}\n\nhorIntra_syn_spec = {\"synapse_model\": \"AMPA\",\n \"weight\": 1.0,\n \"delay\": nest.random.uniform(min=1.75, max=2.25)}\n\n# In a loop, we run over the sources and targets and the corresponding\n# dictionaries with values that needs updating.\nfor conn in [{\"sources\": \"L23pyr\", \"targets\": \"L23pyr\", \"conn_spec\": {}},\n {\"sources\": \"L23pyr\", \"targets\": \"L23in\", \"conn_spec\": {}},\n {\"sources\": \"L4pyr\", \"targets\": \"L4pyr\", \"conn_spec\": {\"mask\": {\"circular\": {\"radius\": 7.0 * dpc}}}},\n {\"sources\": \"L4pyr\", \"targets\": \"L4in\", \"conn_spec\": {\"mask\": {\"circular\": {\"radius\": 7.0 * dpc}}}},\n {\"sources\": \"L56pyr\", \"targets\": \"L56pyr\", \"conn_spec\": {}},\n {\"sources\": \"L56pyr\", \"targets\": \"L56in\", \"conn_spec\": {}}]:\n conn_spec = horIntra_conn_spec.copy()\n conn_spec.update(conn['conn_spec'])\n ccConnections.append([conn['sources'], conn['targets'], conn_spec, horIntra_syn_spec])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vertical intralaminar\n*Note:* \"Vertical\" means \"between cortical layers\" in this\ncase.\n\nWe proceed as above.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "verIntra_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"p\": nest.spatial_distributions.gaussian(nest.spatial.distance, std=7.5 * dpc)}\n\nverIntra_syn_spec = {\"synapse_model\": \"AMPA\",\n \"weight\": 2.0,\n \"delay\": nest.random.uniform(min=1.75, max=2.25)}\n\nfor conn in [{\"sources\": \"L23pyr\", \"targets\": \"L56pyr\",\n \"syn_spec\": {\"weight\": 1.0}},\n {\"sources\": \"L23pyr\", \"targets\": \"L23in\",\n \"syn_spec\": {\"weight\": 1.0}},\n {\"sources\": \"L4pyr\", \"targets\": \"L23pyr\", \"syn_spec\": {}},\n {\"sources\": \"L4pyr\", \"targets\": \"L23in\", \"syn_spec\": {}},\n {\"sources\": \"L56pyr\", \"targets\": \"L23pyr\", \"syn_spec\": {}},\n {\"sources\": \"L56pyr\", \"targets\": \"L23in\", \"syn_spec\": {}},\n {\"sources\": \"L56pyr\", \"targets\": \"L4pyr\", \"syn_spec\": {}},\n {\"sources\": \"L56pyr\", \"targets\": \"L4in\", \"syn_spec\": {}}]:\n syn_spec = verIntra_syn_spec.copy()\n syn_spec.update(conn['syn_spec'])\n ccConnections.append([conn['sources'], conn['targets'], verIntra_conn_spec, syn_spec])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Intracortical inhibitory\n\nWe proceed as above, with the following difference: each connection\nis added to both the same-orientation and the cross-orientation list of\nconnections.\n\n**Note:** Weights increased from -1.0 to -2.0, to make up for missing GabaB\n\nNote that we have to specify the **weight with negative sign** to make\nthe connections inhibitory.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "intraInh_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}},\n \"p\": 0.25*nest.spatial_distributions.gaussian(nest.spatial.distance, std=7.5 * dpc)}\n\nintraInh_syn_spec = {\"synapse_model\": \"GABA_A\",\n \"weight\": -2.0,\n \"delay\": nest.random.uniform(min=1.75, max=2.25)}\n\nfor conn in [{\"sources\": \"L23in\", \"targets\": \"L23pyr\", \"conn_spec\": {}},\n {\"sources\": \"L23in\", \"targets\": \"L23in\", \"conn_spec\": {}},\n {\"sources\": \"L4in\", \"targets\": \"L4pyr\", \"conn_spec\": {}},\n {\"sources\": \"L4in\", \"targets\": \"L4in\", \"conn_spec\": {}},\n {\"sources\": \"L56in\", \"targets\": \"L56pyr\", \"conn_spec\": {}},\n {\"sources\": \"L56in\", \"targets\": \"L56in\", \"conn_spec\": {}}]:\n conn_spec = intraInh_conn_spec.copy()\n conn_spec.update(conn['conn_spec'])\n ccConnections.append([conn['sources'], conn['targets'], conn_spec, intraInh_syn_spec])\n ccxConnections.append([conn['sources'], conn['targets'], conn_spec, intraInh_syn_spec])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Corticothalamic\nWe proceed as above.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "corThal_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"mask\": {\"circular\": {\"radius\": 5.0 * dpc}},\n \"p\": 0.5*nest.spatial_distributions.gaussian(nest.spatial.distance, std=7.5 * dpc)}\n\ncorThal_syn_spec = {\"synapse_model\": \"AMPA\",\n \"weight\": 1.0,\n \"delay\": nest.random.uniform(min=7.5, max=8.5)}\n\nfor conn in [{\"sources\": \"L56pyr\", \"conn_spec\": {}}]:\n conn_spec = intraInh_conn_spec.copy()\n conn_spec.update(conn['conn_spec'])\n syn_spec = intraInh_syn_spec.copy()\n ctConnections.append([conn['sources'], conn_spec, syn_spec])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Corticoreticular\n\nIn this case, there is only a single connection, so we define the\ndictionaries directly; it is very similar to corThal, and to show that,\nwe copy first, then update.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "corRet = corThal_conn_spec.copy()\ncorRet_syn_spec = corThal_syn_spec.copy()\ncorRet_syn_spec.update({\"weight\": 2.5})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build all connections beginning in cortex\n\nCortico-cortical, same orientation\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"Connecting: cortico-cortical, same orientation\")\nfor source, target, conn_spec, syn_spec in ccConnections:\n for src_i in range(name_dict[source]):\n for tgt_i in range(name_dict[target]):\n source_name = '{}_{}'.format(source, src_i)\n target_name = '{}_{}'.format(target, tgt_i)\n nest.Connect(Vp_h_layers[source_name], Vp_h_layers[target_name], conn_spec, syn_spec)\n nest.Connect(Vp_v_layers[source_name], Vp_v_layers[target_name], conn_spec, syn_spec)\n\n# Cortico-cortical, cross-orientation\nprint(\"Connecting: cortico-cortical, other orientation\")\nfor source, target, conn_spec, syn_spec in ccxConnections:\n for src_i in range(name_dict[source]):\n for tgt_i in range(name_dict[target]):\n source_name = '{}_{}'.format(source, src_i)\n target_name = '{}_{}'.format(target, tgt_i)\n nest.Connect(Vp_h_layers[source_name], Vp_v_layers[target_name], conn_spec, syn_spec)\n nest.Connect(Vp_v_layers[source_name], Vp_h_layers[target_name], conn_spec, syn_spec)\n\n# Cortico-thalamic connections\nprint(\"Connecting: cortico-thalamic\")\nfor source, conn_spec, syn_spec in ctConnections:\n for src_i in range(name_dict[source]):\n source_name = '{}_{}'.format(source, src_i)\n nest.Connect(Vp_h_layers[source_name], TpRelay, conn_spec, syn_spec)\n nest.Connect(Vp_h_layers[source_name], TpInter, conn_spec, syn_spec)\n nest.Connect(Vp_v_layers[source_name], TpRelay, conn_spec, syn_spec)\n nest.Connect(Vp_v_layers[source_name], TpInter, conn_spec, syn_spec)\n\nfor src_i in range(name_dict['L56pyr']):\n source_name = 'L56pyr_{}'.format(src_i)\n nest.Connect(Vp_h_layers[source_name], Rp, corRet, corRet_syn_spec)\n nest.Connect(Vp_v_layers[source_name], Rp, corRet, corRet_syn_spec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thalamo-cortical connections\n\n**Note:** According to the text on p. 1674, bottom right, of the Hill &\nTononi paper, thalamocortical connections are created by selecting from\nthe thalamic population for each L4 pyramidal cell. We must therefore\nspecify that we want to select from the source neurons.\n\nWe first handle the rectangular thalamocortical connections.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "thalCorRect_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"use_on_source\": True}\n\nthalCorRect_syn_spec = {\"synapse_model\": \"AMPA\",\n \"weight\": 5.0,\n \"delay\": nest.random.uniform(min=2.75, max=3.25)}\n\nprint(\"Connecting: thalamo-cortical\")\n\n# Horizontally tuned\nthalCorRect_conn_spec.update(\n {\"mask\": {\"rectangular\": {\"lower_left\": [-4.0 * dpc, -1.0 * dpc],\n \"upper_right\": [4.0 * dpc, 1.0 * dpc]}}})\n\nfor conn in [{\"targets\": \"L4pyr\", \"conn_spec\": {\"p\": 0.5}},\n {\"targets\": \"L56pyr\", \"conn_spec\": {\"p\": 0.3}}]:\n conn_spec = thalCorRect_conn_spec.copy()\n conn_spec.update(conn['conn_spec'])\n for trg_i in range(name_dict[conn['targets']]):\n target_name = '{}_{}'.format(conn['targets'], trg_i)\n nest.Connect(\n TpRelay, Vp_h_layers[target_name], conn_spec, thalCorRect_syn_spec)\n\n# Vertically tuned\nthalCorRect_conn_spec.update(\n {\"mask\": {\"rectangular\": {\"lower_left\": [-1.0 * dpc, -4.0 * dpc],\n \"upper_right\": [1.0 * dpc, 4.0 * dpc]}}})\n\nfor conn in [{\"targets\": \"L4pyr\", \"conn_spec\": {\"p\": 0.5}},\n {\"targets\": \"L56pyr\", \"conn_spec\": {\"p\": 0.3}}]:\n conn_spec = thalCorRect_conn_spec.copy()\n conn_spec.update(conn['conn_spec'])\n for trg_i in range(name_dict[conn['targets']]):\n target_name = '{}_{}'.format(conn['targets'], trg_i)\n nest.Connect(\n TpRelay, Vp_v_layers[target_name], conn_spec, thalCorRect_syn_spec)\n\n# Diffuse connections\nthalCorDiff_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"use_on_source\": True,\n \"mask\": {\"circular\": {\"radius\": 5.0 * dpc}},\n \"p\": 0.1*nest.spatial_distributions.gaussian(nest.spatial.distance, std=7.5*dpc)}\n\nthalCorDiff_syn_spec = {\"synapse_model\": \"AMPA\",\n \"weight\": 5.0,\n \"delay\": nest.random.uniform(min=2.75, max=3.25)}\n\nfor conn in [{\"targets\": \"L4pyr\"},\n {\"targets\": \"L56pyr\"}]:\n for trg_i in range(name_dict[conn['targets']]):\n target_name = '{}_{}'.format(conn['targets'], trg_i)\n nest.Connect(TpRelay, Vp_h_layers[target_name], thalCorDiff_conn_spec, thalCorDiff_syn_spec)\n nest.Connect(TpRelay, Vp_v_layers[target_name], thalCorDiff_conn_spec, thalCorDiff_syn_spec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thalamic connections\n\nConnections inside thalamus, including Rp.\n\n*Note:* In Hill & Tononi, the inhibition between Rp cells is mediated by\nGABA_B receptors. We use GABA_A receptors here to provide some\nself-dampening of Rp.\n\n**Note 1:** The following code had a serious bug in v. 0.1: During the first\niteration of the loop, \"synapse_model\" and \"weights\" were set to \"AMPA\" and\n\"0.1\", respectively and remained unchanged, so that all connections were\ncreated as excitatory connections, even though they should have been\ninhibitory. We now specify synapse_model and weight explicitly for each\nconnection to avoid this.\n\n**Note 2:** The following code also had a serious bug in v. 0.4: In the\nloop the connection dictionary would be updated directly, i.e. without\nmaking a copy. This lead to the entry ``'sources': 'TpInter'`` being\nleft in the dictionary when connecting with ``Rp`` sources. Therefore no\nconnections for the connections with ``Rp`` as source would be created\nhere.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "thal_conn_spec = {\"rule\": \"pairwise_bernoulli\"}\nthal_syn_spec = {\"delay\": nest.random.uniform(min=1.75, max=2.25)}\n\nprint(\"Connecting: intra-thalamic\")\n\nfor src, tgt, conn, syn in [(TpRelay, Rp,\n {\"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"p\": nest.spatial_distributions.gaussian(\n nest.spatial.distance, std=7.5*dpc)},\n {\"synapse_model\": \"AMPA\",\n \"weight\": 2.0}),\n (TpInter, TpRelay,\n {\"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"p\": 0.25*nest.spatial_distributions.gaussian(\n nest.spatial.distance, std=7.5*dpc)},\n {\"synapse_model\": \"GABA_A\",\n \"weight\": -1.0}),\n (TpInter, TpInter,\n {\"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"p\": 0.25*nest.spatial_distributions.gaussian(\n nest.spatial.distance, std=7.5*dpc)},\n {\"synapse_model\": \"GABA_A\", \"weight\": -1.0}),\n (Rp, TpRelay, {\"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"p\": 0.15*nest.spatial_distributions.gaussian(\n nest.spatial.distance, std=7.5*dpc)},\n {\"synapse_model\": \"GABA_A\", \"weight\": -1.0}),\n (Rp, TpInter, {\"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"p\": 0.15*nest.spatial_distributions.gaussian(\n nest.spatial.distance, std=7.5*dpc)},\n {\"synapse_model\": \"GABA_A\", \"weight\": -1.0}),\n (Rp, Rp, {\"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"p\": 0.5*nest.spatial_distributions.gaussian(\n nest.spatial.distance, std=7.5*dpc)},\n {\"synapse_model\": \"GABA_A\", \"weight\": -1.0})\n ]:\n conn_spec = thal_conn_spec.copy()\n conn_spec.update(conn)\n syn_spec = thal_syn_spec.copy()\n syn_spec.update(syn)\n nest.Connect(src, tgt, conn_spec, syn_spec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thalamic input\n\nInput to the thalamus from the retina.\n\n**Note:** Hill & Tononi specify a delay of 0 ms for this connection.\nWe use 1 ms here.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "retThal_conn_spec = {\"rule\": \"pairwise_bernoulli\",\n \"mask\": {\"circular\": {\"radius\": 1.0 * dpc}},\n \"p\": 0.75*nest.spatial_distributions.gaussian(nest.spatial.distance, std=2.5*dpc)}\n\nretThal_syn_spec = {\"weight\": 10.0,\n \"delay\": 1.0,\n \"synapse_model\": \"AMPA\"}\n\nprint(\"Connecting: retino-thalamic\")\n\nnest.Connect(retina, TpRelay, retThal_conn_spec, retThal_syn_spec)\nnest.Connect(retina, TpInter, retThal_conn_spec, retThal_syn_spec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checks on connections\n\nAs a very simple check on the connections created, we inspect\nthe connections from the central node of various layers.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Connections from Retina to TpRelay\nretina_ctr_node_id = nest.FindCenterElement(retina)\nretina_ctr_index = retina.index(retina_ctr_node_id.global_id)\nconns = nest.GetConnections(retina[retina_ctr_index], TpRelay)\nnest.PlotTargets(retina[retina_ctr_index], TpRelay, 'AMPA')\nplt.title('Connections Retina -> TpRelay')\n\n# Connections from TpRelay to L4pyr in Vp (horizontally tuned)\nTpRelay_ctr_node_id = nest.FindCenterElement(TpRelay)\nTpRelay_ctr_index = TpRelay.index(TpRelay_ctr_node_id.global_id)\nnest.PlotTargets(TpRelay[TpRelay_ctr_index], Vp_h_layers['L4pyr_0'], 'AMPA')\nplt.title('Connections TpRelay -> Vp(h) L4pyr')\n\n# Connections from TpRelay to L4pyr in Vp (vertically tuned)\nnest.PlotTargets(TpRelay[TpRelay_ctr_index], Vp_v_layers['L4pyr_0'], 'AMPA')\nplt.title('Connections TpRelay -> Vp(v) L4pyr')\n\n# Block until the figures are closed before we continue.\nplt.show(block=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Recording devices\n\nThis recording device setup is a bit makeshift. For each population\nwe want to record from, we create one ``multimeter``, then select\nall nodes of the right model from the target population and\nconnect. ``loc`` is the subplot location for the layer.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"Connecting: Recording devices\")\nrecorders = {}\n\nfor name, loc, population in [('TpRelay', 0, TpRelay),\n ('Rp', 1, Rp),\n ('Vp_v L4pyr 1', 2, Vp_v_layers['L4pyr_0']),\n ('Vp_v L4pyr 2', 3, Vp_v_layers['L4pyr_1']),\n ('Vp_h L4pyr 1', 4, Vp_h_layers['L4pyr_0']),\n ('Vp_h L4pyr 2', 5, Vp_h_layers['L4pyr_1'])]:\n recorders[name] = (nest.Create('RecordingNode'), loc)\n # one recorder to all targets\n nest.Connect(recorders[name][0], population)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example simulation\n\nThis simulation is set up to create a step-wise visualization of\nthe membrane potential. To do so, we simulate ``sim_interval``\nmilliseconds at a time, then read out data from the multimeters,\nclear data from the multimeters and plot the data as pseudocolor\nplots.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# show time during simulation\nnest.SetKernelStatus({'print_time': True})\n\n# lower and upper limits for color scale, for each of the\n# populations recorded.\nvmn = [-80, -80, -80, -80, -80, -80]\nvmx = [-50, -50, -50, -50, -50, -50]\n\nnest.Simulate(Params['sim_interval'])\n\n# Set up the figure, assume six recorders.\nfig, axes = plt.subplots(2, 3)\nimages = []\n\nfor i, ax in enumerate(axes.flat):\n # We initialize with an empty image\n images.append(ax.imshow([[0.]], aspect='equal', interpolation='nearest',\n extent=(0, Params['N'] + 1, 0, Params['N'] + 1),\n vmin=vmn[i], vmax=vmx[i], cmap='plasma'))\n fig.colorbar(images[-1], ax=ax)\n\n# loop over simulation intervals\nfor t in np.arange(0, Params['simtime'], Params['sim_interval']):\n\n # do the simulation\n nest.Simulate(Params['sim_interval'])\n\n # now plot data from each recorder in turn\n for name, rec_item in recorders.items():\n recorder, subplot_pos = rec_item\n ax = axes.flat[subplot_pos]\n im = images[subplot_pos]\n\n d = recorder.get('events', 'V_m')\n # clear data from multimeter\n recorder.n_events = 0\n\n # update image data and title\n im.set_data(np.reshape(d, (Params['N'], Params['N'])))\n ax.set_title(name + ', t = %6.1f ms' % nest.GetKernelStatus()['time'])\n\n # We need to pause because drawing of the figure happens while the main code is sleeping\n plt.pause(0.0001)\n\n# just for some information at the end\npprint(nest.GetKernelStatus())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.12" } }, "nbformat": 4, "nbformat_minor": 0 }