diff --git a/examples/tsodyksmarkramstp/tmevaluator.py b/examples/tsodyksmarkramstp/tmevaluator.py new file mode 100644 index 0000000..75f2a0e --- /dev/null +++ b/examples/tsodyksmarkramstp/tmevaluator.py @@ -0,0 +1,78 @@ +"""Tsodyks-Markram model Evaluator. +This module contains an evaluator for the Tsodyks-Markram model. + +@author: Giuseppe Chindemi +@remark: Copyright (c) 2017, EPFL/Blue Brain Project + This file is part of BluePyOpt + This library is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License version 3.0 as + published by the Free Software Foundation. + This library is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + See the GNU Lesser General Public License for more details. + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +""" +import bluepyopt as bpop +import numpy as np +import tmodeint + + +class TsodyksMarkramEvaluator(bpop.evaluators.Evaluator): + def __init__(self, t, v, tstim, params): + """ + Parameters + ---------- + t : numpy.ndarray + Time vector. + v : numpy.ndarray + Voltage vector, must have same dimension of t. + tstim : numpy.ndarray + Time of the stimuli. + params : list + List of parameters to fit. Every entry must be a tuple + (name, lower bound, upper bound). + """ + super(TsodyksMarkramEvaluator, self).__init__() + self.v = v + self.t = t + self.stimidx = np.searchsorted(t, tstim) + self.dx = t[1] - t[0] + self.nsamples = len(v) + self.params = params + # Find voltage baseline + bs_stop = np.searchsorted(t, tstim[0]) + self.vrest = np.mean(v[:bs_stop]) + # Compute time windows where to compare model and data + offset = 0.005 # s + window = 0.04 # s + window_samples = int(np.round(window/self.dx)) + psp_start = np.searchsorted(t, tstim+offset) + psp_stop = psp_start + window_samples + psp_stop[-1] += 2*window_samples # Extend last psp window (RTR case) + self.split_idx = zip(psp_start, psp_stop) + # Parameters to be optimized + self.params = [bpop.parameters.Parameter(name, bounds=(minval, maxval)) + for name, minval, maxval in self.params] + # Objectives + self.objectives = [bpop.objectives.Objective('interval_%d' % (i,)) + for i in xrange(len(self.split_idx))] + + def generate_model(self, individual): + v, _ = tmodeint.integrate(self.stimidx, self.nsamples, self.dx, + self.vrest, *individual) + return v + + def evaluate_with_lists(self, individual): + candidate_v = self.generate_model(individual) + errors = [np.linalg.norm(self.v[t0:t1] - candidate_v[t0:t1]) + for t0, t1 in self.split_idx] + return errors + + def evaluate_with_lists_aggregate(self, param_values): + errors = self.evaluate_with_lists(param_values) + return np.linalg.norm(errors) + diff --git a/examples/tsodyksmarkramstp/tmodeint.py b/examples/tsodyksmarkramstp/tmodeint.py new file mode 100644 index 0000000..a61cc8a --- /dev/null +++ b/examples/tsodyksmarkramstp/tmodeint.py @@ -0,0 +1,90 @@ +"""Tsodyks-Markram model ODEINT. +This module contains functions to numerically integrate the Tsodyks-Markram +model. The current implementation is a port of an Igor Pro (WaveMetrics) +routine, written by Rodrigo Perin with the help of Raphael Holzer. Many thanks +to Misha Tsodyks and Henry Markram for their support during all development +stages. + +@author: Giuseppe Chindemi, Rodrigo Perin +@remark: Copyright (c) 2017, EPFL/Blue Brain Project + This file is part of BluePyOpt + This library is free software; you can redistribute it and/or modify it + under the terms of the GNU Lesser General Public License version 3.0 as + published by the Free Software Foundation. + This library is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + See the GNU Lesser General Public License for more details. + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +""" +import numpy as np + + +def integrate(sampstim, nsamples, dt, vRest, Trec, Tfac, ASE, USE, Rinput, Tmem, + Tinac, latency): + """Integrate Tsodyks-Markram model and produce corresponding voltage trace. + + Parameters + ---------- + sampstim : numpy.ndarray + Array of sample indices where stimulation occurs. + nsamples : int + Number of samples in the trace. + dt : float + Time step of the trace. + vRest : float + Membrane resting potential (V). + Trec : float + Recovery time constant (ms). + Tfac : float + Facilitation time constant (ms). + ASE : float + Absolute Synaptic Efficacy. + USE : float + Utilization of Synaptic Efficacy. Has to be in the interval [0, 1]. + Rinput : float + Input resistance (MOhm). + Tmem : float + Membrane time constant (ms). + Tinac : float + Inactivation time constant (ms). + latency : float + PSP latency (ms). + + Returns + ------- + vtrace : numpy.ndarray + Voltage trace corresponding to the input parameters. + tm_statevar : dictionary + Integral of all state variables. + """ + AP = np.zeros(nsamples) + I = np.zeros(nsamples) + R = np.zeros(nsamples) + E = np.zeros(nsamples) + U = np.zeros(nsamples) + P = np.zeros(nsamples) + # Set stimuli in AP vector + psp_offset = int(np.round(1e-3*latency/dt)) + AP[sampstim+psp_offset] = 1 + # Initialize state vectors + R[0] = 1 + E[0] = 0 + P[0] = 0 + U[0] = USE + # Integrate TM model ODE + for i in xrange(1, nsamples): + R[i] = R[i-1] + dt*(1-R[i-1]-E[i-1])*1e3/Trec - U[i-1]*R[i-1]*AP[i-1] + E[i] = E[i-1] - dt*E[i-1]*1e3/Tinac + U[i-1]*R[i-1]*AP[i-1] + U[i] = U[i-1] - dt*(U[i-1]-USE)*1e3/Tfac + USE*(1-U[i-1])*AP[i-1] + P[i] = P[i-1] + dt*(Rinput*ASE/10**6*E[i-1]-P[i-1])*1e3/Tmem + # Update state + P = P + vRest + E = E*ASE/10**12 + I = 1-R-E + tm_statevar = {'recovered': R, 'effective': E, 'used': U, 'inactive': I} + return P, tm_statevar + diff --git a/examples/tsodyksmarkramstp/trace.pkl b/examples/tsodyksmarkramstp/trace.pkl new file mode 100644 index 0000000..0c23664 Binary files /dev/null and b/examples/tsodyksmarkramstp/trace.pkl differ diff --git a/examples/tsodyksmarkramstp/tsodyksmarkramstp.ipynb b/examples/tsodyksmarkramstp/tsodyksmarkramstp.ipynb new file mode 100644 index 0000000..5845e04 --- /dev/null +++ b/examples/tsodyksmarkramstp/tsodyksmarkramstp.ipynb @@ -0,0 +1,1737 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Tsodyks-Markram model of short-term synaptic plasticity\n", + "In this notebook we demonstrate how to fit the parameters of the Tsodyks-Markram model to a given in vitro somatic recording. The in vitro trace used here shows a typical L5TTPC-L5TTPC depressing connection, kindly provided by Rodrigo Perin (EPFL)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\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", + " 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", + " this.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 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);\n", + " canvas.attr('height', height);\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 = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Get v trace for best model\n", + "vmodel = evaluator.generate_model(best)\n", + "\n", + "fig, ax = plt.subplots()\n", + "ax.plot(trace['t'], trace['v'], label='in vitro')\n", + "ax.plot(trace['t'], vmodel, label='model')\n", + "ax.legend(loc=0)\n", + "ax.set_xlabel('time (s)')\n", + "ax.set_ylabel('soma voltage (V)')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}