{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Quickstart

\n", "To run the code below:\n", "
    \n", "
  1. Click on the cell to select it.
    2. \n", "
  2. Press SHIFT+ENTER on your keyboard or press the play button\n", " () in the toolbar above
    3. \n", "
\n", "Feel free to create new cells using the plus button\n", "(), or pressing SHIFT+ENTER while this cell\n", "is selected.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This article demonstrates how a control flow, where simulation parameters depend on the results of previous simulations, can be expressed by making use of standard control structures in Python. By having access to the full expressivity of a general purpose programming language, expressing such control flow is straight-forward; this would not be the case for a declarative model description.\n", "\n", "Our goal in this toy example is to find the threshold voltage of neuron as a function of the density of sodium channels.\n", "\n", "This example is from our eLife paper [(Stimberg et al. 2019)](https://elifesciences.org/articles/47314).\n", "\n", "\n", "We start with the basic setup:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from brian2 import *\n", "%matplotlib notebook\n", "\n", "defaultclock.dt = 0.01*ms # small time step for stiff equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our model of the neuron is based on the classical model of from Hodgkin and Huxley (1952). Note that this is not actually a model of a neuron, but rather of a (space-clamped) axon. However, to avoid confusion with spatially extended models, we simply use the term \"neuron\" here. In this model, the membrane potential is shifted, i.e. the resting potential is at 0mV:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "El = 10.613*mV\n", "ENa = 115*mV\n", "EK = -12*mV\n", "gl = 0.3*msiemens/cm**2\n", "gK = 36*msiemens/cm**2\n", "gNa_max = 100*msiemens/cm**2\n", "gNa_min = 15*msiemens/cm**2\n", "C = 1*uF/cm**2\n", "\n", "eqs = '''\n", "dv/dt = (gl * (El-v) + gNa * m**3 * h * (ENa-v) + gK * n**4 * (EK-v)) / C : volt\n", "gNa : siemens/meter**2\n", "dm/dt = alpham * (1-m) - betam * m : 1\n", "dn/dt = alphan * (1-n) - betan * n : 1\n", "dh/dt = alphah * (1-h) - betah * h : 1\n", "alpham = (0.1/mV) * (-v+25*mV) / (exp((-v+25*mV) / (10*mV)) - 1)/ms : Hz\n", "betam = 4 * exp(-v/(18*mV))/ms : Hz\n", "alphah = 0.07 * exp(-v/(20*mV))/ms : Hz\n", "betah = 1/(exp((-v+30*mV) / (10*mV)) + 1)/ms : Hz\n", "alphan = (0.01/mV) * (-v+10*mV) / (exp((-v+10*mV) / (10*mV)) - 1)/ms : Hz\n", "betan = 0.125*exp(-v/(80*mV))/ms : Hz\n", "'''" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We simulate 100 neurons at the same time, each of them having a density of sodium channels between 15 and 100 mS/cm²:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "neurons = NeuronGroup(100, eqs, method='rk4', threshold='v>50*mV')\n", "neurons.gNa = 'gNa_min + (gNa_max - gNa_min)*1.0*i/N'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We initialize the state variables to their resting state values, note that the values for $m$, $n$, $h$ depend on the values of $\\alpha_m$, $\\beta_m$, etc. which themselves depend on $v$. The order of the assignments ($v$ is initialized before $m$, $n$, and $h$) therefore matters, something that is naturally expressed by stating initial values as sequential assignments to the state variables. In a declarative approach, this would be potentially ambiguous." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "neurons.v = 0*mV\n", "neurons.m = '1/(1 + betam/alpham)'\n", "neurons.n = '1/(1 + betan/alphan)'\n", "neurons.h = '1/(1 + betah/alphah)'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We record the spiking activity of the neurons and store the current network state so that we can later restore it and run another iteration of our experiment:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "S = SpikeMonitor(neurons)\n", "store()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The algorithm we use here to find the voltage threshold is a simple bisection: we try to find the threshold voltage of a neuron by repeatedly testing values and increasing or decreasing these values depending on whether we observe a spike or not. By continously halving the size of the correction, we quickly converge to a precise estimate.\n", "\n", "We start with the same initial estimate for all segments, 25mV above the resting potential, and the same value for the size of the \"correction step\":" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "v0 = 25*mV*ones(len(neurons))\n", "step = 25*mV" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For later visualization of how the estimates converged towards their final values, we also store the intermediate values of the estimates:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "estimates = np.full((11, len(neurons)), np.nan)*mV\n", "estimates[0, :] = v0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now run 10 iterations of our algorithm:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "for i in range(10):\n", " # Reset to the initial state\n", " restore()\n", " # Set the membrane potential to our threshold estimate\n", " neurons.v = v0\n", " # Run the simulation for 20ms\n", " run(20*ms)\n", " # Decrease the estimates for neurons that spiked\n", " v0[S.count > 0] -= step\n", " # Increase the estimate for neurons that did not spike\n", " v0[S.count == 0] += step\n", " # Reduce step size and store current estimate\n", " step /= 2.0\n", " estimates[i + 1, :] = v0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After the 10 iteration steps, we plot the results:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('