diff --git a/notebooks/experiment_example.ipynb b/notebooks/experiment_example.ipynb new file mode 100644 index 0000000..a50e0d9 --- /dev/null +++ b/notebooks/experiment_example.ipynb @@ -0,0 +1,345 @@ +{ + "nbformat": 4, + "nbformat_minor": 0, + "metadata": { + "colab": { + "provenance": [], + "toc_visible": true, + "authorship_tag": "ABX9TyMvW7lhxBXG2gi+iXDgU0PT", + "include_colab_link": true + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + }, + "language_info": { + "name": "python" + } + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "view-in-github", + "colab_type": "text" + }, + "source": [ + "\"Open" + ] + }, + { + "cell_type": "markdown", + "source": [ + "#Install CompNeuroPy and ANNarchy" + ], + "metadata": { + "id": "MQwRmyFdmUuD" + } + }, + { + "cell_type": "code", + "source": [ + "!pip install CompNeuroPy\n", + "!git clone https://github.com/ANNarchy/ANNarchy && cd ANNarchy && git checkout develop && pip install .\n", + "!rm -rf ANNarchy" + ], + "metadata": { + "collapsed": true, + "id": "VGC0ujzTm-z3" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "#Imports and setup ANNarchy timestep" + ], + "metadata": { + "id": "Mcoz2oq76flU" + } + }, + { + "cell_type": "code", + "source": [ + "from IPython.display import Image, display\n", + "from CompNeuroPy import (\n", + " CompNeuroExp,\n", + " CompNeuroMonitors,\n", + " CompNeuroModel,\n", + " current_step,\n", + " current_ramp,\n", + " PlotRecordings,\n", + ")\n", + "from CompNeuroPy.full_models import HHmodelBischop\n", + "from ANNarchy import dt, setup, get_population\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "setup(dt=0.01)" + ], + "metadata": { + "id": "bMhApHy-6iJ3" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Define the Experiment\n", + "A single run resets the model, sets the leakage potential of population 0 and runs a current ramp and current step." + ], + "metadata": { + "id": "pyf2gSTyN80I" + } + }, + { + "cell_type": "code", + "source": [ + "class MyExp(CompNeuroExp):\n", + "\n", + " def run(self, model: CompNeuroModel, E_L: float):\n", + " # PREPARE RUN\n", + " self.reset()\n", + " self.monitors.start()\n", + " # SET E_L PARAMETER\n", + " get_population(model.populations[0]).E_L = E_L\n", + " # SIMULATION\n", + " ret_current_ramp = current_ramp(pop=model.populations[0], a0=0, a1=100, dur=1000, n=50)\n", + " self.reset(parameters=False)\n", + " ret_current_step = current_step(pop=model.populations[0], t1=500, t2=500, a1=0, a2=50)\n", + " # OPTIONAL DATA OF RUN\n", + " self.data[\"population_name\"] = model.populations[0]\n", + " self.data[\"time_step\"] = dt()\n", + " self.data[\"current_arr\"] = np.concatenate([ret_current_ramp[\"current_arr\"], ret_current_step[\"current_arr\"]])\n", + " # RETURN RESULTS\n", + " return self.results()" + ], + "metadata": { + "id": "8coZOl3oOc1B" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Create/Compile Model\n", + "The model is a single population (consisting of 1 neuron) of a Hodgkin & Huxley neuron. The `HHmodelBischop` class is a child of the `CompNeuroModel` class with a predefined model creation function." + ], + "metadata": { + "id": "XCbqszREO0U2" + } + }, + { + "cell_type": "code", + "source": [ + "model = HHmodelBischop()\n", + "model.populations" + ], + "metadata": { + "id": "gHmfSnNlO2WP" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Initialize the experiment\n", + "Recording the membrane potential of the models first population." + ], + "metadata": { + "id": "EQL2vn40PA_O" + } + }, + { + "cell_type": "code", + "source": [ + "my_exp = MyExp(monitors=CompNeuroMonitors({model.populations[0]: [\"v\"]}))" + ], + "metadata": { + "id": "I6SphV4qPHTU" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Set the \"reset-state\" for the experiment\n", + "Set the membrane potential of the model to -90 mV." + ], + "metadata": { + "id": "bJn_adfcPRz5" + } + }, + { + "cell_type": "code", + "source": [ + "print(f\"Compilation state v = {get_population(model.populations[0]).v}\")\n", + "get_population(model.populations[0]).v = -90.0\n", + "print(f\"Changed state v = {get_population(model.populations[0]).v}\")\n", + "my_exp.store_model_state(compartment_list=model.populations)" + ], + "metadata": { + "id": "wyE8p0SmPZ0c" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Run the experiment twice with different leakage potentials" + ], + "metadata": { + "id": "QOuRAfGhPeyT" + } + }, + { + "cell_type": "code", + "source": [ + "results_run1: CompNeuroExp._ResultsCl = my_exp.run(model=model, E_L=-68.0)\n", + "results_run2: CompNeuroExp._ResultsCl = my_exp.run(model=model, E_L=-90.0)" + ], + "metadata": { + "id": "PKqJScoJPf7q" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# PlotRecordings\n", + "This allows to easily get overview plots of the recordings of a single recording chunk." + ], + "metadata": { + "id": "Xg_52O9uPkgd" + } + }, + { + "cell_type": "code", + "source": [ + "for chunk in [0,1]:\n", + " PlotRecordings(\n", + " figname=f\"example_experiment_chunk_{chunk}.png\",\n", + " recordings=results_run1.recordings,\n", + " recording_times=results_run1.recording_times,\n", + " chunk=chunk,\n", + " shape=(1, 1),\n", + " plan={\n", + " \"position\": [1],\n", + " \"compartment\": [results_run1.data[\"population_name\"]],\n", + " \"variable\": [\"v\"],\n", + " \"format\": [\"line\"],\n", + " },\n", + " )" + ], + "metadata": { + "id": "ftBPfNZgPmok" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "source": [ + "display(Image(filename='example_experiment_chunk_0.png', height=700))" + ], + "metadata": { + "id": "o8Lvap-8QFIL" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "source": [ + "display(Image(filename='example_experiment_chunk_1.png', height=700))" + ], + "metadata": { + "id": "p4-rWESeQJkm" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Get Data and Time Arrays\n", + "Each experiment run created 2 recording chunks. They all start at time 0 (because of resetting the model, see above). The function combine_chunks() can be used to combine the chunks into a single recording time and value array." + ], + "metadata": { + "id": "OkEO6pyePzmT" + } + }, + { + "cell_type": "code", + "source": [ + "time_arr1, data_arr1 = results_run1.recording_times.combine_chunks(\n", + " recordings=results_run1.recordings,\n", + " recording_data_str=f\"{results_run1.data['population_name']};v\",\n", + " mode=\"consecutive\",\n", + ")\n", + "time_arr2, data_arr2 = results_run2.recording_times.combine_chunks(\n", + " recordings=results_run2.recordings,\n", + " recording_data_str=f\"{results_run2.data['population_name']};v\",\n", + " mode=\"consecutive\",\n", + ")\n", + "current_arr = results_run1.data[\"current_arr\"]" + ], + "metadata": { + "id": "SoiA-7cMP2bj" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "source": [ + "# Plot Data of Both Runs" + ], + "metadata": { + "id": "mUso4pycQTLf" + } + }, + { + "cell_type": "code", + "source": [ + "plt.figure()\n", + "plt.subplot(211)\n", + "plt.plot(time_arr1, data_arr1, label=\"E_L = -68.0\")\n", + "plt.plot(time_arr2, data_arr2, label=\"E_L = -90.0\")\n", + "plt.plot(\n", + " [time_arr1[0], time_arr1[-1]], [-90, -90], ls=\"dotted\", label=\"initial v = -90.0\"\n", + ")\n", + "plt.legend()\n", + "plt.ylabel(\"Membrane potential [mV]\")\n", + "plt.subplot(212)\n", + "plt.plot(time_arr1, current_arr, \"k--\")\n", + "plt.ylabel(\"Input current\")\n", + "plt.xlabel(\"Time [ms]\")\n", + "plt.tight_layout()\n", + "plt.savefig(\"example_experiment_combined.png\")" + ], + "metadata": { + "id": "61hXVuUDQYJ2" + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "source": [ + "display(Image(filename='example_experiment_chunk_1.png', height=700))" + ], + "metadata": { + "id": "aMVQLAi-QcFG" + }, + "execution_count": null, + "outputs": [] + } + ] +} \ No newline at end of file diff --git a/src/CompNeuroPy/simulation_functions.py b/src/CompNeuroPy/simulation_functions.py index a674150..96c0546 100644 --- a/src/CompNeuroPy/simulation_functions.py +++ b/src/CompNeuroPy/simulation_functions.py @@ -24,6 +24,11 @@ def attr_sim(pop: str, attr_dict, t=500): dictionary containing the attributes and their values t (int): duration in ms + + Returns: + attr_list_dict (dict): + dictionary containing the attribute values for each time step, keys are the + attribute names """ ### save prev attr @@ -42,6 +47,10 @@ def attr_sim(pop: str, attr_dict, t=500): for attr, v in v_prev_dict.items(): setattr(get_population(pop), attr, v) + ### return the values for the attribute for each time step + attr_list_dict = {attr: [v] * int(round(t / dt())) for attr, v in attr_dict.items()} + return attr_list_dict + def attribute_step(pop: str, attr, t1=500, t2=500, v1=0, v2=100): """ @@ -66,16 +75,18 @@ def attribute_step(pop: str, attr, t1=500, t2=500, v1=0, v2=100): dictionary containing: - duration (int): duration of the simulation + - v_arr (np.array): array of attribute values for each time step """ - + v_list = [] ### first/pre step simulation - attr_sim(pop, {attr: v1}, t=t1) - + attr_list_dict = attr_sim(pop, {attr: v1}, t=t1) + v_list.extend(attr_list_dict[attr]) ### second/post step simulation - attr_sim(pop, {attr: v2}, t=t2) + attr_list_dict = attr_sim(pop, {attr: v2}, t=t2) + v_list.extend(attr_list_dict[attr]) - ### return duration of the simulation - return {"duration": t1 + t2} + ### return duration of the simulation and the attribute values + return {"duration": t1 + t2, "v_arr": np.array(v_list)} def attr_ramp(pop: str, attr, v0, v1, dur, n): @@ -107,6 +118,7 @@ def attr_ramp(pop: str, attr, v0, v1, dur, n): - dv (int): step size of attribute - dur_stim (int): duration of single steps + - v_arr (np.array): array of attribute values for each time step Raises: ValueError: if resulting duration of one stimulation is not divisible by the @@ -125,11 +137,13 @@ def attr_ramp(pop: str, attr, v0, v1, dur, n): dv = (v1 - v0) / (n - 1) # for n stimulations only n-1 steps occur dur_stim = dur / n v = v0 + v_list = [] for _ in range(n): - attr_sim(pop, attr_dict={attr: v}, t=dur_stim) + attr_list_dict = attr_sim(pop, attr_dict={attr: v}, t=dur_stim) + v_list.extend(attr_list_dict[attr]) v = v + dv - return {"dv": dv, "dur_stim": dur_stim} + return {"dv": dv, "dur_stim": dur_stim, "v_arr": np.array(v_list)} def increasing_attr(pop: str, attr, v0, dv, nr_steps, dur_step): @@ -155,15 +169,18 @@ def increasing_attr(pop: str, attr, v0, dv, nr_steps, dur_step): dictionary containing: - attr_list (list): list of attribute values for each step simulation + - v_arr (np.array): array of attribute values for each time step """ attr_list = [] v = v0 + v_list = [] for _ in range(nr_steps): attr_list.append(v) - attr_sim(pop, {attr: v}, t=dur_step) + attr_list_dict = attr_sim(pop, {attr: v}, t=dur_step) + v_list.extend(attr_list_dict[attr]) v += dv - return {"attr_list": attr_list} + return {"attr_list": attr_list, "v_arr": np.array(v_list)} def current_step(pop: str, t1=500, t2=500, a1=0, a2=100): @@ -188,8 +205,13 @@ def current_step(pop: str, t1=500, t2=500, a1=0, a2=100): dictionary containing: - duration (int): duration of the simulation + - current_arr (np.array): array of current values for each time step """ - return attribute_step(pop, "I_app", t1=t1, t2=t2, v1=a1, v2=a2) + attribute_step_ret = attribute_step(pop, "I_app", t1=t1, t2=t2, v1=a1, v2=a2) + return { + "duration": attribute_step_ret["duration"], + "current_arr": attribute_step_ret["v_arr"], + } def current_stim(pop: str, t=500, a=100): @@ -206,8 +228,13 @@ def current_stim(pop: str, t=500, a=100): duration in ms a (int): current amplitude + + Returns: + current_arr (np.array): + array of current values for each time step """ - attr_sim(pop, {"I_app": a}, t=t) + attr_list_dict = attr_sim(pop, {"I_app": a}, t=t) + return np.array(attr_list_dict["I_app"]) def current_ramp(pop: str, a0, a1, dur, n): @@ -239,13 +266,18 @@ def current_ramp(pop: str, a0, a1, dur, n): - da (int): current step size - dur_stim (int): duration of one stimulation + - current_arr (np.array): array of current values for each time step Raises: ValueError: if resulting duration of one stimulation is not divisible by the simulation time step without remainder """ attr_ramp_return = attr_ramp(pop, "I_app", a0, a1, dur, n) - return {"da": attr_ramp_return["dv"], "dur_stim": attr_ramp_return["dur_stim"]} + return { + "da": attr_ramp_return["dv"], + "dur_stim": attr_ramp_return["dur_stim"], + "current_arr": attr_ramp_return["v_arr"], + } def increasing_current(pop: str, a0, da, nr_steps, dur_step): @@ -272,9 +304,13 @@ def increasing_current(pop: str, a0, da, nr_steps, dur_step): dictionary containing: - current_list (list): list of current amplitudes for each stimulation + - current_arr (np.array): array of current values for each time step """ increasing_attr_return = increasing_attr(pop, "I_app", a0, da, nr_steps, dur_step) - return {"current_list": increasing_attr_return["attr_list"]} + return { + "current_list": increasing_attr_return["attr_list"], + "current_arr": increasing_attr_return["v_arr"], + } class SimulationEvents: @@ -285,6 +321,11 @@ class SimulationEvents: SimulationEvents. Do never simulate within the effect functions of the events. The simulation is done between the events. + !!! warning + The onset of events and trigger times should be given in simulation steps (not + in ms). The 'end' event has to be triggered to end the simulation (otherwise + it will be triggered right at the beginning). + Example: ```python from CompNeuroPy import SimulationEvents @@ -433,9 +474,15 @@ def run(self): Run the simulation. The simulation runs until the end event is triggered. The simulation can be run multiple times by calling this function multiple times. """ + ### for all events with given onset, change onset to current step + onset + ### (otherwise run would need to be called at time 0) + for event in self.event_list: + if event.onset is not None: + event.onset = get_current_step() + event.onset + ### check if there are events which have no onset and are not triggered by other ### events and have no model_trigger --> they would never start - ### --> set their onset to current step --> they ar run directly after calling run + ### --> set their onset to current step --> they are run directly after calling run triggered_events = [] for event in self.event_list: if event.trigger is not None: @@ -523,9 +570,8 @@ def _run_model_trigger_events(self): """ ### loop to check if model trigger got active for model_trigger in self.model_trigger_list: - if ( - int(get_population(model_trigger).decision[0]) == -1 - ): ### TODO this is not generalized yet, only works if the model_trigger populations have the variable decision which is set to -1 if the model trigger is active + ### TODO this is not generalized yet, only works if the model_trigger populations have the variable decision which is set to -1 if the model trigger is active + if int(get_population(model_trigger).decision[0]) == -1: ### -1 means got active ### find the events triggerd by the model_trigger and run them for event in self.event_list: