{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Structural Plasticity example\n\nThis example shows a simple network of two populations where structural\nplasticity is used. The network has 1000 neurons, 80% excitatory and\n20% inhibitory. The simulation starts without any connectivity. A set of\nhomeostatic rules are defined, according to which structural plasticity will\ncreate and delete synapses dynamically during the simulation until a desired\nlevel of electrical activity is reached. The model of structural plasticity\nused here corresponds to the formulation presented in [1]_.\n\nAt the end of the simulation, a plot of the evolution of the connectivity\nin the network and the average calcium concentration in the neurons is created.\n\n## References\n\n.. [1] Butz, M., and van Ooyen, A. (2013). A simple rule for dendritic spine and axonal bouton formation can\n       account for cortical reorganization after focal retinal lesions. PLoS Comput. Biol. 9 (10), e1003259.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we have import all necessary modules.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport numpy\nimport matplotlib.pyplot as plt\nimport sys"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We define general simulation parameters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class StructralPlasticityExample:\n\n    def __init__(self):\n        # simulated time (ms)\n        self.t_sim = 200000.0\n        # simulation step (ms).\n        self.dt = 0.1\n        self.number_excitatory_neurons = 800\n        self.number_inhibitory_neurons = 200\n\n        # Structural_plasticity properties\n        self.update_interval = 10000.0\n        self.record_interval = 1000.0\n        # rate of background Poisson input\n        self.bg_rate = 10000.0\n        self.neuron_model = 'iaf_psc_exp'"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this implementation of structural plasticity, neurons grow\nconnection points called synaptic elements. Synapses can be created\nbetween compatible synaptic elements. The growth of these elements is\nguided by homeostatic rules, defined as growth curves.\nHere we specify the growth curves for synaptic elements of excitatory\nand inhibitory neurons.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Excitatory synaptic elements of excitatory neurons\n        self.growth_curve_e_e = {\n            'growth_curve': \"gaussian\",\n            'growth_rate': 0.0001,  # (elements/ms)\n            'continuous': False,\n            'eta': 0.0,  # Ca2+\n            'eps': 0.05,  # Ca2+\n        }\n\n        # Inhibitory synaptic elements of excitatory neurons\n        self.growth_curve_e_i = {\n            'growth_curve': \"gaussian\",\n            'growth_rate': 0.0001,  # (elements/ms)\n            'continuous': False,\n            'eta': 0.0,  # Ca2+\n            'eps': self.growth_curve_e_e['eps'],  # Ca2+\n        }\n\n        # Excitatory synaptic elements of inhibitory neurons\n        self.growth_curve_i_e = {\n            'growth_curve': \"gaussian\",\n            'growth_rate': 0.0004,  # (elements/ms)\n            'continuous': False,\n            'eta': 0.0,  # Ca2+\n            'eps': 0.2,  # Ca2+\n        }\n\n        # Inhibitory synaptic elements of inhibitory neurons\n        self.growth_curve_i_i = {\n            'growth_curve': \"gaussian\",\n            'growth_rate': 0.0001,  # (elements/ms)\n            'continuous': False,\n            'eta': 0.0,  # Ca2+\n            'eps': self.growth_curve_i_e['eps']  # Ca2+\n        }\n\n        # Now we specify the neuron model.\n\n        self.model_params = {'tau_m': 10.0,  # membrane time constant (ms)\n                             # excitatory synaptic time constant (ms)\n                             'tau_syn_ex': 0.5,\n                             # inhibitory synaptic time constant (ms)\n                             'tau_syn_in': 0.5,\n                             't_ref': 2.0,  # absolute refractory period (ms)\n                             'E_L': -65.0,  # resting membrane potential (mV)\n                             'V_th': -50.0,  # spike threshold (mV)\n                             'C_m': 250.0,  # membrane capacitance (pF)\n                             'V_reset': -65.0  # reset potential (mV)\n                             }\n\n        self.nodes_e = None\n        self.nodes_i = None\n        self.mean_ca_e = []\n        self.mean_ca_i = []\n        self.total_connections_e = []\n        self.total_connections_i = []"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We initialize variables for the postsynaptic currents of the\nexcitatory, inhibitory, and external synapses. These values were\ncalculated from a PSP amplitude of 1 for excitatory synapses,\n-1 for inhibitory synapses and 0.11 for external synapses.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "self.psc_e = 585.0\n        self.psc_i = -585.0\n        self.psc_ext = 6.2\n\n    def prepare_simulation(self):\n        nest.ResetKernel()\n        nest.set_verbosity('M_ERROR')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We set global kernel parameters. Here we define the resolution\nfor the simulation, which is also the time resolution for the update\nof the synaptic elements.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.resolution = self.dt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set Structural Plasticity synaptic update interval which is how often\nthe connectivity will be updated inside the network. It is important\nto notice that synaptic elements and connections change on different\ntime scales.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.structural_plasticity_update_interval = self.update_interval"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we define Structural Plasticity synapses. In this example we create\ntwo synapse models, one for excitatory and one for inhibitory synapses.\nThen we define that excitatory synapses can only be created between a\npre-synaptic element called `Axon_ex` and a postsynaptic element\ncalled `Den_ex`. In a similar manner, synaptic elements for inhibitory\nsynapses are defined.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.CopyModel('static_synapse', 'synapse_ex')\n        nest.SetDefaults('synapse_ex', {'weight': self.psc_e, 'delay': 1.0})\n        nest.CopyModel('static_synapse', 'synapse_in')\n        nest.SetDefaults('synapse_in', {'weight': self.psc_i, 'delay': 1.0})\n        nest.structural_plasticity_synapses = {\n            'synapse_ex': {\n                'synapse_model': 'synapse_ex',\n                'post_synaptic_element': 'Den_ex',\n                'pre_synaptic_element': 'Axon_ex'\n            },\n            'synapse_in': {\n                'synapse_model': 'synapse_in',\n                'post_synaptic_element': 'Den_in',\n                'pre_synaptic_element': 'Axon_in'\n            }\n        }\n\n    def create_nodes(self):\n        \"\"\"\n        Assign growth curves to synaptic elements\n        \"\"\"\n\n        synaptic_elements = {\n            'Den_ex': self.growth_curve_e_e,\n            'Den_in': self.growth_curve_e_i,\n            'Axon_ex': self.growth_curve_e_e,\n        }\n\n        synaptic_elements_i = {\n            'Den_ex': self.growth_curve_i_e,\n            'Den_in': self.growth_curve_i_i,\n            'Axon_in': self.growth_curve_i_i,\n        }"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then it is time to create a population with 80% of the total network\nsize excitatory neurons and another one with 20% of the total network\nsize of inhibitory neurons.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "self.nodes_e = nest.Create('iaf_psc_alpha',\n                                   self.number_excitatory_neurons,\n                                   {'synaptic_elements': synaptic_elements})\n\n        self.nodes_i = nest.Create('iaf_psc_alpha',\n                                   self.number_inhibitory_neurons,\n                                   {'synaptic_elements': synaptic_elements_i})\n        self.nodes_e.synaptic_elements = synaptic_elements\n        self.nodes_i.synaptic_elements = synaptic_elements_i\n\n    def connect_external_input(self):\n        \"\"\"\n        We create and connect the Poisson generator for external input\n        \"\"\"\n        noise = nest.Create('poisson_generator')\n        noise.rate = self.bg_rate\n        nest.Connect(noise, self.nodes_e, 'all_to_all',\n                     {'weight': self.psc_ext, 'delay': 1.0})\n        nest.Connect(noise, self.nodes_i, 'all_to_all',\n                     {'weight': self.psc_ext, 'delay': 1.0})\n\n    ####################################################################################\n    # In order to save the amount of average calcium concentration in each\n    # population through time we create the function ``record_ca``. Here we use\n    # the value of `Ca` for every neuron in the network and then\n    # store the average.\n    def record_ca(self):\n        ca_e = self.nodes_e.Ca,  # Calcium concentration\n        self.mean_ca_e.append(numpy.mean(ca_e))\n\n        ca_i = self.nodes_i.Ca,  # Calcium concentration\n        self.mean_ca_i.append(numpy.mean(ca_i))\n\n    ####################################################################################\n    # In order to save the state of the connectivity in the network through time\n    # we create the function ``record_connectivity``. Here we retrieve the number\n    # of connected pre-synaptic elements of each neuron. The total amount of\n    # excitatory connections is equal to the total amount of connected excitatory\n    # pre-synaptic elements. The same applies for inhibitory connections.\n    def record_connectivity(self):\n        syn_elems_e = self.nodes_e.synaptic_elements\n        syn_elems_i = self.nodes_i.synaptic_elements\n        self.total_connections_e.append(sum(neuron['Axon_ex']['z_connected']\n                                            for neuron in syn_elems_e))\n        self.total_connections_i.append(sum(neuron['Axon_in']['z_connected']\n                                            for neuron in syn_elems_i))\n\n    ####################################################################################\n    # We define a function to plot the recorded values\n    # at the end of the simulation.\n    def plot_data(self):\n        fig, ax1 = plt.subplots()\n        ax1.axhline(self.growth_curve_e_e['eps'],\n                    linewidth=4.0, color='#9999FF')\n        ax1.plot(self.mean_ca_e, 'b',\n                 label='Ca Concentration Excitatory Neurons', linewidth=2.0)\n        ax1.axhline(self.growth_curve_i_e['eps'],\n                    linewidth=4.0, color='#FF9999')\n        ax1.plot(self.mean_ca_i, 'r',\n                 label='Ca Concentration Inhibitory Neurons', linewidth=2.0)\n        ax1.set_ylim([0, 0.275])\n        ax1.set_xlabel(\"Time in [s]\")\n        ax1.set_ylabel(\"Ca concentration\")\n        ax2 = ax1.twinx()\n        ax2.plot(self.total_connections_e, 'm',\n                 label='Excitatory connections', linewidth=2.0, linestyle='--')\n        ax2.plot(self.total_connections_i, 'k',\n                 label='Inhibitory connections', linewidth=2.0, linestyle='--')\n        ax2.set_ylim([0, 2500])\n        ax2.set_ylabel(\"Connections\")\n        ax1.legend(loc=1)\n        ax2.legend(loc=4)\n        plt.savefig('StructuralPlasticityExample.eps', format='eps')\n\n    ####################################################################################\n    # It is time to specify how we want to perform the simulation. In this\n    # function we first enable structural plasticity in the network and then we\n    # simulate in steps. On each step we record the calcium concentration and the\n    # connectivity. At the end of the simulation, the plot of connections and\n    # calcium concentration through time is generated.\n    def simulate(self):\n        if nest.NumProcesses() > 1:\n            sys.exit(\"For simplicity, this example only works \" +\n                     \"for a single process.\")\n        nest.EnableStructuralPlasticity()\n        print(\"Starting simulation\")\n        sim_steps = numpy.arange(0, self.t_sim, self.record_interval)\n        for i, step in enumerate(sim_steps):\n            nest.Simulate(self.record_interval)\n            self.record_ca()\n            self.record_connectivity()\n            if i % 20 == 0:\n                print(\"Progress: \" + str(i / 2) + \"%\")\n        print(\"Simulation finished successfully\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we take all the functions that we have defined and create the sequence\nfor our example. We prepare the simulation, create the nodes for the network,\nconnect the external input and then simulate. Please note that as we are\nsimulating 200 biological seconds in this example, it will take a few minutes\nto complete.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "if __name__ == '__main__':\n    example = StructralPlasticityExample()\n    # Prepare simulation\n    example.prepare_simulation()\n    example.create_nodes()\n    example.connect_external_input()\n    # Start simulation\n    example.simulate()\n    example.plot_data()"
      ]
    }
  ],
  "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}