{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# PyGeNN implementation of local cortical microcircuit model\nThis example model is a reimplementation of the microcircuit model developed by Tobias C. Potjans and Markus Diesmann [Potjans2014]_.\nIt is a full-scale spiking network model of the local cortical microcircuit.\nThe simulated spontaneous activity is asynchronous irregular and cell-type specific firing rates are in agreement\nwith in vivo recordings in awake animals, including the low rate of layer 2/3 excitatory cells.\nThis example can be used as follows:\n\n.. argparse::\n :filename: ../userproject/potjans_microcircuit.py\n :func: get_parser\n :prog: potjans_microcircuit\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np \n\nfrom argparse import ArgumentParser\nfrom pygenn import (GeNNModel, VarLocation, init_postsynaptic,\n init_sparse_connectivity, init_weight_update, init_var)\nfrom scipy.stats import norm\nfrom time import perf_counter\n\n# ----------------------------------------------------------------------------\n# Parameters\n# ----------------------------------------------------------------------------\n# Layer names\nLAYER_NAMES = [\"23\", \"4\", \"5\", \"6\"]\n\n# Population names\nPOPULATION_NAMES = [\"E\", \"I\"]\n\n# Simulation timestep [ms]\nDT_MS = 0.1\n\n# How many threads to use per spike for procedural connectivity?\nNUM_THREADS_PER_SPIKE = 8\n\n# Background rate per synapse\nBACKGROUND_RATE = 8.0 # spikes/s\n\n# Relative inhibitory synaptic weight\nG = -4.0\n\n# Mean synaptic weight for all excitatory projections except L4e->L2/3e\nMEAN_W = 87.8e-3 # nA\nEXTERNAL_W = 87.8e-3 # nA\n\n# Mean synaptic weight for L4e->L2/3e connections\n# See p. 801 of the paper, second paragraph under 'Model Parameterization',\n# and the caption to Supplementary Fig. 7\nLAYER_23_4_W = 2.0 * MEAN_W # nA\n\n# Standard deviation of weight distribution relative to mean for\n# all projections except L4e->L2/3e\nREL_W = 0.1\n\n# Standard deviation of weight distribution relative to mean for L4e->L2/3e\n# This value is not mentioned in the paper, but is chosen to match the\n# original code by Tobias Potjans\nLAYER_23_4_RELW = 0.05\n\n# Numbers of neurons in full-scale model\nNUM_NEURONS = {\n \"23\": {\"E\":20683, \"I\": 5834},\n \"4\": {\"E\":21915, \"I\": 5479},\n \"5\": {\"E\":4850, \"I\": 1065},\n \"6\": {\"E\":14395, \"I\": 2948}}\n\n# Probabilities for >=1 connection between neurons in the given populations.\n# The first index is for the target population; the second for the source population\nCONNECTION_PROBABILTIES = {\n \"23E\": {\"23E\": 0.1009, \"23I\": 0.1689, \"4E\": 0.0437, \"4I\": 0.0818, \"5E\": 0.0323, \"5I\": 0.0, \"6E\": 0.0076, \"6I\": 0.0},\n \"23I\": {\"23E\": 0.1346, \"23I\": 0.1371, \"4E\": 0.0316, \"4I\": 0.0515, \"5E\": 0.0755, \"5I\": 0.0, \"6E\": 0.0042, \"6I\": 0.0},\n \"4E\": {\"23E\": 0.0077, \"23I\": 0.0059, \"4E\": 0.0497, \"4I\": 0.135, \"5E\": 0.0067, \"5I\": 0.0003, \"6E\": 0.0453, \"6I\": 0.0},\n \"4I\": {\"23E\": 0.0691, \"23I\": 0.0029, \"4E\": 0.0794, \"4I\": 0.1597, \"5E\": 0.0033, \"5I\": 0.0, \"6E\": 0.1057, \"6I\": 0.0},\n \"5E\": {\"23E\": 0.1004, \"23I\": 0.0622, \"4E\": 0.0505, \"4I\": 0.0057, \"5E\": 0.0831, \"5I\": 0.3726, \"6E\": 0.0204, \"6I\": 0.0},\n \"5I\": {\"23E\": 0.0548, \"23I\": 0.0269, \"4E\": 0.0257, \"4I\": 0.0022, \"5E\": 0.06, \"5I\": 0.3158, \"6E\": 0.0086, \"6I\": 0.0},\n \"6E\": {\"23E\": 0.0156, \"23I\": 0.0066, \"4E\": 0.0211, \"4I\": 0.0166, \"5E\": 0.0572, \"5I\": 0.0197, \"6E\": 0.0396, \"6I\": 0.2252},\n \"6I\": {\"23E\": 0.0364, \"23I\": 0.001, \"4E\": 0.0034, \"4I\": 0.0005, \"5E\": 0.0277, \"5I\": 0.008, \"6E\": 0.0658, \"6I\": 0.1443}}\n \n\n# In-degrees for external inputs\nNUM_EXTERNAL_INPUTS = {\n \"23\": {\"E\": 1600, \"I\": 1500},\n \"4\": {\"E\": 2100, \"I\": 1900},\n \"5\": {\"E\": 2000, \"I\": 1900},\n \"6\": {\"E\": 2900, \"I\": 2100}}\n\n# Mean rates in the full-scale model, necessary for scaling\n# Precise values differ somewhat between network realizations\nMEAN_FIRING_RATES = {\n \"23\": {\"E\": 0.971, \"I\": 2.868},\n \"4\": {\"E\": 4.746, \"I\": 5.396},\n \"5\": {\"E\": 8.142, \"I\": 9.078},\n \"6\": {\"E\": 0.991, \"I\": 7.523}}\n\n# Means and standard deviations of delays from given source populations (ms)\nMEAN_DELAY = {\"E\": 1.5, \"I\": 0.75}\n\nDELAY_SD = {\"E\": 0.75, \"I\": 0.375}\n\n# ----------------------------------------------------------------------------\n# Helper functions\n# ----------------------------------------------------------------------------\ndef get_scaled_num_neurons(layer, pop, neuron_scale):\n return int(round(neuron_scale * NUM_NEURONS[layer][pop]))\n\ndef get_full_num_inputs(src_layer, src_pop, trg_layer, trg_pop):\n num_src = NUM_NEURONS[src_layer][src_pop]\n num_trg = NUM_NEURONS[trg_layer][trg_pop]\n connection_prob = CONNECTION_PROBABILTIES[trg_layer + trg_pop][src_layer + src_pop]\n\n return int(round(np.log(1.0 - connection_prob) / np.log(float(num_trg * num_src - 1) / float(num_trg * num_src))) / num_trg)\n\ndef get_mean_weight(src_layer, src_pop, trg_layer, trg_pop):\n # Determine mean weight\n if src_pop == \"E\":\n if src_layer == \"4\" and trg_layer == \"23\" and trg_pop == \"E\":\n return LAYER_23_4_W\n else:\n return MEAN_W\n else:\n return G * MEAN_W\n\ndef get_scaled_num_connections(src_layer, src_pop, trg_layer, trg_pop,\n neuron_scale, connectivity_scale):\n # Scale full number of inputs by scaling factor\n num_inputs = get_full_num_inputs(src_layer, src_pop, trg_layer, trg_pop) * connectivity_scale\n assert num_inputs >= 0.0\n\n # Multiply this by number of postsynaptic neurons\n num_neurons = get_scaled_num_neurons(trg_layer, trg_pop, neuron_scale)\n return int(round(num_inputs * float(num_neurons)))\n\ndef get_full_mean_input_current(layer, pop):\n # Loop through source populations\n mean_input_current = 0.0\n for src_layer in LAYER_NAMES:\n for src_pop in POPULATION_NAMES:\n mean_input_current += (get_mean_weight(src_layer, src_pop, layer, pop) *\n get_full_num_inputs(src_layer, src_pop, layer, pop) *\n MEAN_FIRING_RATES[src_layer][src_pop])\n\n # Add mean external input current\n mean_input_current += EXTERNAL_W * NUM_EXTERNAL_INPUTS[layer][pop] * BACKGROUND_RATE\n assert mean_input_current >= 0.0\n return mean_input_current\n\ndef get_parser():\n parser = ArgumentParser()\n parser.add_argument(\"--duration\", type=float, default=1000.0, help=\"Duration to simulate (ms)\")\n parser.add_argument(\"--neuron-scale\", type=float, default=1.0, help=\"Scaling factor to apply to number of neurons\")\n parser.add_argument(\"--connectivity-scale\", type=float, default=1.0, help=\"Scaling factor to apply to number of neurons\")\n parser.add_argument(\"--kernel-profiling\", action=\"store_true\", help=\"Output kernel profiling data\")\n parser.add_argument(\"--procedural-connectivity\", action=\"store_true\", help=\"Use procedural connectivity\")\n parser.add_argument(\"--save-data\", action=\"store_true\", help=\"Save spike data (rather than plotting it)\")\n return parser\n\n# ----------------------------------------------------------------------------\n# Entry point\n# ----------------------------------------------------------------------------\nif __name__ == \"__main__\":\n args = get_parser().parse_args()\n\n # ----------------------------------------------------------------------------\n # Network creation\n # ----------------------------------------------------------------------------\n model = GeNNModel(\"float\", \"potjans_microcircuit\")\n model.dt = DT_MS\n model.fuse_postsynaptic_models = True\n model.default_narrow_sparse_ind_enabled = True\n model.timing_enabled = args.kernel_profiling\n model.default_var_location = VarLocation.DEVICE\n model.default_sparse_connectivity_location = VarLocation.DEVICE\n\n lif_init = {\"V\": init_var(\"Normal\", {\"mean\": -58.0, \"sd\": 5.0}), \"RefracTime\": 0.0}\n poisson_init = {\"current\": 0.0}\n\n exp_curr_init = init_postsynaptic(\"ExpCurr\", {\"tau\": 0.5})\n\n quantile = 0.9999\n normal_quantile_cdf = norm.ppf(quantile)\n max_delay = {pop: MEAN_DELAY[pop] + (DELAY_SD[pop] * normal_quantile_cdf)\n for pop in POPULATION_NAMES}\n print(\"Max excitatory delay:%fms , max inhibitory delay:%fms\" % (max_delay[\"E\"], max_delay[\"I\"]))\n\n # Calculate maximum dendritic delay slots\n # **NOTE** it seems inefficient using maximum for all but this allows more aggressive merging of postsynaptic models\n max_dendritic_delay_slots = int(round(max(max_delay.values()) / DT_MS))\n print(\"Max dendritic delay slots:%d\" % max_dendritic_delay_slots)\n\n print(\"Creating neuron populations:\")\n total_neurons = 0\n neuron_populations = {}\n for layer in LAYER_NAMES:\n for pop in POPULATION_NAMES:\n pop_name = layer + pop\n\n # Calculate external input rate, weight and current\n ext_input_rate = NUM_EXTERNAL_INPUTS[layer][pop] * args.connectivity_scale * BACKGROUND_RATE\n ext_weight = EXTERNAL_W / np.sqrt(args.connectivity_scale)\n ext_input_current = 0.001 * 0.5 * (1.0 - np.sqrt(args.connectivity_scale)) * get_full_mean_input_current(layer, pop)\n assert ext_input_current >= 0.0\n\n lif_params = {\"C\": 0.25, \"TauM\": 10.0, \"Vrest\": -65.0, \"Vreset\": -65.0, \"Vthresh\" : -50.0,\n \"Ioffset\": ext_input_current, \"TauRefrac\": 2.0}\n poisson_params = {\"weight\": ext_weight, \"tauSyn\": 0.5, \"rate\": ext_input_rate}\n\n pop_size = get_scaled_num_neurons(layer, pop, args.neuron_scale)\n neuron_pop = model.add_neuron_population(pop_name, pop_size, \"LIF\", lif_params, lif_init)\n model.add_current_source(pop_name + \"_poisson\", \"PoissonExp\", neuron_pop, poisson_params, poisson_init)\n\n # Enable spike recording\n neuron_pop.spike_recording_enabled = True\n\n print(\"\\tPopulation %s: num neurons:%u, external input rate:%f, external weight:%f, external DC offset:%f\" % (pop_name, pop_size, ext_input_rate, ext_weight, ext_input_current))\n\n # Add number of neurons to total\n total_neurons += pop_size\n\n # Add neuron population to dictionary\n neuron_populations[pop_name] = neuron_pop\n\n # Loop through target populations and layers\n print(\"Creating synapse populations:\")\n total_synapses = 0\n num_sub_rows = NUM_THREADS_PER_SPIKE if args.procedural_connectivity else 1\n for trg_layer in LAYER_NAMES:\n for trg_pop in POPULATION_NAMES:\n trg_name = trg_layer + trg_pop\n\n # Loop through source populations and layers\n for src_layer in LAYER_NAMES:\n for src_pop in POPULATION_NAMES:\n src_name = src_layer + src_pop\n\n # Determine mean weight\n mean_weight = get_mean_weight(src_layer, src_pop, trg_layer, trg_pop) / np.sqrt(args.connectivity_scale)\n\n # Determine weight standard deviation\n if src_pop == \"E\" and src_layer == \"4\" and trg_layer == \"23\" and trg_pop == \"E\":\n weight_sd = mean_weight * LAYER_23_4_RELW\n else:\n weight_sd = abs(mean_weight * REL_W)\n\n # Calculate number of connections\n num_connections = get_scaled_num_connections(src_layer, src_pop,\n trg_layer, trg_pop,\n args.neuron_scale, args.connectivity_scale)\n\n if num_connections > 0:\n print(\"\\tConnection between '%s' and '%s': numConnections=%u, meanWeight=%f, weightSD=%f, meanDelay=%f, delaySD=%f\"\n % (src_name, trg_name, num_connections, mean_weight, weight_sd, MEAN_DELAY[src_pop], DELAY_SD[src_pop]))\n\n # Build parameters for fixed number total connector\n connect_params = {\"num\": num_connections}\n\n # Build distribution for delay parameters\n d_dist = {\"mean\": MEAN_DELAY[src_pop], \"sd\": DELAY_SD[src_pop], \"min\": 0.0, \"max\": max_delay[src_pop]}\n\n total_synapses += num_connections\n\n # Build unique synapse name\n synapse_name = src_name + \"_\" + trg_name\n\n matrix_type = \"PROCEDURAL\" if args.procedural_connectivity else \"SPARSE\"\n\n # Excitatory\n if src_pop == \"E\":\n # Build distribution for weight parameters\n # **HACK** np.float32 doesn't seem to automatically cast\n w_dist = {\"mean\": mean_weight, \"sd\": weight_sd, \"min\": 0.0, \"max\": float(np.finfo(np.float32).max)}\n\n # Create weight parameters\n static_synapse_init = init_weight_update(\"StaticPulseDendriticDelay\", {},\n {\"g\": init_var(\"NormalClipped\", w_dist),\n \"d\": init_var(\"NormalClippedDelay\", d_dist)})\n # Add synapse population\n syn_pop = model.add_synapse_population(synapse_name, matrix_type,\n neuron_populations[src_name], neuron_populations[trg_name],\n static_synapse_init, exp_curr_init,\n init_sparse_connectivity(\"FixedNumberTotalWithReplacement\", connect_params))\n\n # Set max dendritic delay and span type\n syn_pop.max_dendritic_delay_timesteps = max_dendritic_delay_slots\n\n if args.procedural_connectivity:\n syn_pop.num_threads_per_spike = NUM_THREADS_PER_SPIKE\n # Inhibitory\n else:\n # Build distribution for weight parameters\n # **HACK** np.float32 doesn't seem to automatically cast\n w_dist = {\"mean\": mean_weight, \"sd\": weight_sd, \"min\": float(-np.finfo(np.float32).max), \"max\": 0.0}\n\n # Create weight parameters\n static_synapse_init = init_weight_update(\"StaticPulseDendriticDelay\", {},\n {\"g\": init_var(\"NormalClipped\", w_dist),\n \"d\": init_var(\"NormalClippedDelay\", d_dist)})\n # Add synapse population\n syn_pop = model.add_synapse_population(synapse_name, matrix_type,\n neuron_populations[src_name], neuron_populations[trg_name],\n static_synapse_init, exp_curr_init,\n init_sparse_connectivity(\"FixedNumberTotalWithReplacement\", connect_params))\n\n # Set max dendritic delay and span type\n syn_pop.max_dendritic_delay_timesteps = max_dendritic_delay_slots\n\n if args.procedural_connectivity:\n syn_pop.num_threads_per_spike = NUM_THREADS_PER_SPIKE\n print(\"Total neurons=%u, total synapses=%u\" % (total_neurons, total_synapses))\n\n print(\"Building Model\")\n model.build()\n\n print(\"Loading Model\")\n duration_timesteps = int(round(args.duration / DT_MS))\n ten_percent_timestep = duration_timesteps // 10\n model.load(num_recording_timesteps=duration_timesteps)\n\n print(\"Simulating\")\n\n # Loop through timesteps\n sim_start_time = perf_counter()\n while model.t < args.duration:\n # Advance simulation\n model.step_time()\n\n # Indicate every 10%\n if (model.timestep % ten_percent_timestep) == 0:\n print(\"%u%%\" % (model.timestep / 100))\n\n sim_end_time = perf_counter()\n\n\n # Download recording data\n model.pull_recording_buffers_from_device()\n\n print(\"Timing:\")\n print(\"\\tSimulation:%f\" % ((sim_end_time - sim_start_time) * 1000.0))\n\n if args.kernel_profiling:\n print(\"\\tInit:%f\" % (1000.0 * model.init_time))\n print(\"\\tSparse init:%f\" % (1000.0 * model.init_sparse_time))\n print(\"\\tNeuron simulation:%f\" % (1000.0 * model.neuron_update_time))\n print(\"\\tSynapse simulation:%f\" % (1000.0 * model.presynaptic_update_time))\n\n if args.save_data:\n # Loop through populations and write spike data to CSV\n for n, pop in neuron_populations.items():\n np.savetxt(f\"{n}_spikes.csv\", np.column_stack(pop.spike_recording_data[0]),\n delimiter=\",\", fmt=(\"%f\", \"%d\"), header=\"Times [ms], Neuron ID\")\n\n else:\n import matplotlib.pyplot as plt\n\n # Create plot\n figure, axes = plt.subplots(1, 2)\n\n # **YUCK** re-order neuron populations for plotting\n ordered_neuron_populations = list(reversed(list(neuron_populations.values())))\n\n start_id = 0\n bar_y = 0.0\n for pop in ordered_neuron_populations:\n # Get recording data\n spike_times, spike_ids = pop.spike_recording_data[0]\n\n # Plot spikes\n actor = axes[0].scatter(spike_times, spike_ids + start_id, s=2, edgecolors=\"none\")\n\n # Plot bar showing rate in matching colour\n axes[1].barh(bar_y, len(spike_times) / (float(pop.num_neurons) * args.duration / 1000.0),\n align=\"center\", color=actor.get_facecolor(), ecolor=\"black\")\n\n # Update offset\n start_id += pop.num_neurons\n\n # Update bar pos\n bar_y += 1.0\n\n axes[0].set_xlabel(\"Time [ms]\")\n axes[0].set_ylabel(\"Neuron number\")\n\n axes[1].set_xlabel(\"Mean firingrate [Hz]\")\n axes[1].set_yticks(np.arange(0.0, len(neuron_populations), 1.0))\n axes[1].set_yticklabels([n.name for n in ordered_neuron_populations])\n\n # Show plot\n plt.show()"
]
}
],
"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.9.5"
}
},
"nbformat": 4,
"nbformat_minor": 0
}