{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Quadratic Integrate-and-Fire (QIF) Neuron Mean-Field Model\n\nHere, we will introduce the QIF population mean-field model, which has been derived from a population of all-to-all\ncoupled QIF neurons in [1]_. The model equations are given by:\n\n\\begin{align}\\tau \\dot r &= \\frac{\\Delta}{\\pi\\tau} + 2 r v, \n\n \\tau \\dot v &= v^2 +\\bar\\eta + I(t) + J r \\tau - (\\pi r \\tau)^2,\\end{align}\n\nwhere $r$ is the average firing rate and $v$ is the average membrane potential of the QIF population [1]_.\nIt is governed by 4 parameters:\n - $\\tau$ --> the population time constant\n - $\\bar \\eta$ --> the mean of a Lorenzian distribution over the neural excitability in the population\n - $\\Delta$ --> the half-width at half maximum of the Lorenzian distribution over the neural excitability\n - $J$ --> the strength of the recurrent coupling inside the population\nThis mean-field model is an exact representation of the macroscopic firing rate and membrane potential dynamics of a\nspiking neural network consisting of QIF neurons with Lorentzian distributed background excitabilities.\nWhile the mean-field derivation is mathematically only valid for all-to-all coupled populations of infinite size,\nit has been shown that there is a close correspondence between the mean-field model and neural populations with\nsparse coupling and population sizes of a few thousand neurons [2]_. In the same work, it has been demonstrated how to\nextend the model by adding synaptic dynamics or additional adaptation currents to the single cell network, that can be\ncarried through the mean-field derivation performed in [1]_. For example, a QIF population with\nspike-frequency adaptation would be given by the following 3D system:\n\n\\begin{align}\\tau \\dot r &= \\frac{\\Delta}{\\pi\\tau} + 2 r v, \n\n \\tau \\dot v &= v^2 +\\bar\\eta + I(t) + J r \\tau - a - (\\pi r \\tau)^2, \n\n \\tau_a \\dot a &= -a + \\alpha \\tau_a r,\\end{align}\n\nwhere the evolution equation for $a$ expresses a convolution of $r$ with a mono-exponential kernel, with\nadaptation strength $\\alpha$ and time constant $\\tau_a$.\n\nIn the sections below, we will demonstrate for each model how to load the model template into pyrates, perform\nsimulations with it and visualize the results.\n\n## References\n\n.. [1] E. Montbri\u00f3, D. Paz\u00f3, A. Roxin (2015) *Macroscopic description for networks of spiking neurons.* Physical\n Review X, 5:021028, https://doi.org/10.1103/PhysRevX.5.021028.\n\n.. [2] R. Gast, H. Schmidt, T.R. Kn\u00f6sche (2020) *A Mean-Field Description of Bursting Dynamics in Spiking Neural\n Networks with Short-Term Adaptation.* Neural Computation 32 (9): 1615-1634.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Basic QIF Model\n\nWe will start out by a tutorial of how to use the QIF model without adaptation. To this end, we will use the\nhigh-level interfaces from PyRates which allow to perform model loading and simulation in a single step.\nFor a step-to-step tutorial of the more low-level interfaces, have a look at the Jansen-Rit example.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 1: Numerical simulation of a the model behavior in time\n\nHere, we use the :code:`integrate` function imported from PyRates. As a first argument to this function, either a path\nto a YAML-based model definition or a :code:`CircuitTemplate` instance can be provided. The function will then compile\nthe model and solve the initial value problem of the above defined differential equations for a time interval from\n0 to the given simulation time. This solution will be calculated numerically by a differential equation solver in\nthe backend, starting with a defined step-size. Here, we use the default backend and solver. Furthermore,\nwe provide a step-function extrinsic input that excites all QIF neurons in a time window from :code:`start` to\n:code:`stop`. This input is defined on a time vector with fixed time steps of size :code:`step_size`.\nCheck out the arguments of the code:`CircuitTemplate.run()` method for a detailed explanation of the\narguments that you can use to adjust this numerical procedure.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from pyrates import integrate\nimport numpy as np\n\n# define simulation time and input start and stop\nT = 100.0\nstep_size = 5e-4\nstart = 20.0\nstop = 80.0\n\n# extrinsic input definition\nsteps = int(np.round(T/step_size))\nI_ext = np.zeros((steps,))\nI_ext[int(start/step_size):int(stop/step_size)] = 3.0\n\n# perform simulation\nresults = integrate(\"model_templates.neural_mass_models.qif.qif\", step_size=step_size, simulation_time=T,\n outputs={'r': 'p/qif_op/r'}, inputs={'p/qif_op/I_ext': I_ext}, clear=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 2: Visualization of the solution\n\nThe output of the :code:`simulate()` function is a :code:`pandas.Dataframe`, which allows for direct plotting of the\ntimeseries it contains.\nThis timeseries represents the numerical solution of the initial value problem solved in step 1 with respect to the\nstate variable $r$ of the model. You can observe that the model expresses hysteresis and converges to a\ndifferent steady-state after the input is turned off, than the steady-state it started from when the input was turned\non. This behavior has been described in detail in [1]_.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\nfig, ax = plt.subplots()\nax.plot(results)\nax.set_xlabel('time')\nax.set_ylabel('r')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# QIF SFA Model\n\nNow, lets have a look at the QIF model with spike-frequency adaptation. We will follow the same steps as outlined\nabove.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"results2 = integrate(\"model_templates.neural_mass_models.qif.qif_sfa\", simulation_time=T, step_size=step_size,\n outputs={'r': 'p/qif_sfa_op/r'}, inputs={'p/qif_sfa_op/I_ext': I_ext}, clear=True)\n\nfig2, ax2 = plt.subplots()\nax2.plot(results2)\nax2.set_xlabel('time')\nax2.set_ylabel('r')\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"you can see that, by adding the adaptation variable to the model, we introduced synchronized bursting behavior to\nthe model. Check out [2]_ if you would like to test out some different parameter regimes and would like to know what\nkind of model behavior to expect if you make changes to the adaptation parameters. To change the parameters, you need\nto derive a new operator template from the given operator template in a yaml file and simply set the parameter you\nwould like to change. For a detailed introduction on how to handle model definitions via YAML files, have a look at\nthe model definition gallery.\n\n"
]
}
],
"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.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}