example_1_pyloric_network_complex.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quickstart\n",
    "To run the code below:\n",
    "\n",
    "1. Click on the cell to select it.\n",
    "2. Press `SHIFT+ENTER` on your keyboard or press the play button\n",
    "   (<button class='fa fa-play icon-play btn btn-xs btn-default'></button>) in the toolbar above.\n",
    "\n",
    "Feel free to create new cells using the plus button\n",
    "(<button class='fa fa-plus icon-plus btn btn-xs btn-default'></button>), or pressing `SHIFT+ENTER` while this cell\n",
    "is selected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1 (Pyloric network of the crustacean stomatogastric ganglion) – biophysically detailed model\n",
    "\n",
    "This example is equivalent to the [simplified model](example_1_pyloric_network.ipynb), but uses a more biophysically detailed neuron model. For a detailed explanation of the modelling approach and the conversion of equations into Brian syntax, see the comments in that file.\n",
    "\n",
    "For details about the neuron model, see:\n",
    "\n",
    "Golowasch, J., Casey, M., Abbott, L. F., & Marder, E. (1999).  \n",
    "Network Stability from Activity-Dependent Regulation of Neuronal Conductances.  \n",
    "Neural Computation, 11(5), 1079-1096.  \n",
    "https://doi.org/10.1162/089976699300016359"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from brian2.only import *\n",
    "import numpy as np\n",
    "\n",
    "set_device('cpp_standalone', build_on_run=False)\n",
    "\n",
    "defaultclock.dt = 0.01*ms\n",
    "\n",
    "init_time = 0.5*second\n",
    "observe_time = 3*second\n",
    "adapt_time = 24 * second\n",
    "\n",
    "### Class-independent constants\n",
    "# Reversal potentials\n",
    "E_L = -68*mV\n",
    "E_Na = 20*mV\n",
    "E_K = -80*mV\n",
    "E_Ca = 120*mV\n",
    "E_proc = -10*mV\n",
    "\n",
    "# Capacitance\n",
    "C_s = 0.2*nF\n",
    "C_a = 0.02*nF\n",
    "\n",
    "# maximal conductances\n",
    "g_E = 10*nS\n",
    "g_La = 7.5*nS\n",
    "g_Na = 300*nS\n",
    "g_Kd = 4*uS\n",
    "G_Ca = 0.2*uS\n",
    "G_K = 16*uS\n",
    "\n",
    "# time constants (independent of V)\n",
    "tau_h_Ca = 150*ms\n",
    "tau_m_A = 0.1*ms\n",
    "tau_h_A = 50*ms\n",
    "tau_m_proc = 6*ms\n",
    "tau_m_Na = 0.025*ms\n",
    "tau_z = 5*second\n",
    "\n",
    "# Synapses\n",
    "s_fast = 0.2/mV\n",
    "V_fast = -50*mV\n",
    "s_slow = 1/mV\n",
    "V_slow = -55*mV\n",
    "E_syn = -75*mV\n",
    "k_1 = 1/ms\n",
    "\n",
    "### Neuronal equations\n",
    "eqs = '''\n",
    "# somatic compartment\n",
    "dV_s/dt = (-I_syn - I_L - I_Ca - I_K - I_A - I_proc - g_E*(V_s - V_a))/C_s : volt\n",
    "\n",
    "I_L = g_Ls*(V_s - E_L) : amp\n",
    "I_K = g_K*m_K**4*(V_s - E_K) : amp\n",
    "I_A = g_A*m_A**3*h_A*(V_s - E_K) : amp\n",
    "I_proc = g_proc*m_proc*(V_s - E_proc) : amp\n",
    "\n",
    "I_syn = I_fast + I_slow: amp\n",
    "I_fast : amp\n",
    "I_slow : amp\n",
    "\n",
    "I_Ca = g_Ca*m_Ca**3*h_Ca*(V_s - E_Ca)            : amp\n",
    "dm_Ca/dt = (m_Ca_inf - m_Ca)/tau_m_Ca            : 1\n",
    "m_Ca_inf = 1/(1 + exp(0.205/mV*(-61.2*mV - V_s))): 1\n",
    "tau_m_Ca = 30*ms -5*ms/(1 + exp(0.2/mV*(-65*mV - V_s))) : second\n",
    "dh_Ca/dt = (h_Ca_inf - h_Ca)/tau_h_Ca            : 1\n",
    "h_Ca_inf = 1/(1 + exp(-0.15/mV*(-75*mV - V_s)))  : 1\n",
    "\n",
    "dm_K/dt = (m_K_inf - m_K)/tau_m_K                : 1\n",
    "m_K_inf = 1/(1 + exp(0.1/mV*(-35*mV - V_s)))     : 1\n",
    "tau_m_K = 2*ms + 55*ms/(1 + exp(-0.125/mV*(-54*mV - V_s))) : second\n",
    "\n",
    "dm_A/dt = (m_A_inf - m_A)/tau_m_A                : 1\n",
    "m_A_inf = 1/(1 + exp(0.2/mV*(-60*mV - V_s)))     : 1\n",
    "dh_A/dt = (h_A_inf - h_A)/tau_h_A                : 1\n",
    "h_A_inf = 1/(1 + exp(-0.18/mV*(-68*mV - V_s)))   : 1\n",
    "\n",
    "dm_proc/dt = (m_proc_inf - m_proc)/tau_m_proc    : 1\n",
    "m_proc_inf = 1/(1 + exp(0.2/mV*(-55*mV - V_s)))  : 1\n",
    "\n",
    "# axonal compartment\n",
    "dV_a/dt = (-g_La*(V_a - E_L) - g_Na*m_Na**3*h_Na*(V_a - E_Na)\n",
    "           -g_Kd*m_Kd**4*(V_a - E_K) - g_E*(V_a - V_s))/C_a : volt\n",
    "\n",
    "dm_Na/dt = (m_Na_inf - m_Na)/tau_m_Na            : 1\n",
    "m_Na_inf = 1/(1 + exp(0.1/mV*(-42.5*mV - V_a)))  : 1\n",
    "dh_Na/dt = (h_Na_inf - h_Na)/tau_h_Na            : 1\n",
    "h_Na_inf = 1/(1 + exp(-0.13/mV*(-50*mV - V_a)))  : 1\n",
    "tau_h_Na = 10*ms/(1 + exp(0.12/mV*(-77*mV - V_a))) : second\n",
    "\n",
    "dm_Kd/dt = (m_Kd_inf - m_Kd)/tau_m_Kd            : 1\n",
    "m_Kd_inf = 1/(1 + exp(0.2/mV*(-41*mV - V_a)))    : 1\n",
    "tau_m_Kd = 12.2*ms + 10.5*ms/(1 + exp(-0.05/mV*(58*mV - V_a))) : second\n",
    "\n",
    "# class-specific fixed maximal conductances\n",
    "g_Ls   : siemens (constant)\n",
    "g_A    : siemens (constant)\n",
    "g_proc : siemens (constant)\n",
    "\n",
    "# Adaptive conductances\n",
    "g_Ca = G_Ca/2*(1 + tanh(z)) : siemens\n",
    "g_K = G_K/2*(1 - tanh(z))   : siemens\n",
    "I_diff = (I_target + I_Ca) : amp\n",
    "dz/dt = tanh(I_diff/nA)/tau_z : 1\n",
    "I_target : amp (constant)\n",
    "\n",
    "# neuron class\n",
    "label : integer (constant)\n",
    "'''\n",
    "ABPD, LP, PY = 0, 1, 2\n",
    "\n",
    "circuit = NeuronGroup(3, eqs, method='rk2',\n",
    "                      threshold='m_Na > 0.5', refractory='m_Na > 0.5')\n",
    "\n",
    "# class-specific constants\n",
    "circuit.label = [ABPD, LP, PY]\n",
    "circuit.I_target = [0.4, 0.3, 0.5]*nA\n",
    "circuit.g_Ls = [30, 25, 15]*nS\n",
    "circuit.g_A = [450, 100, 250]*nS\n",
    "circuit.g_proc = [6, 8, 0]*nS\n",
    "\n",
    "# Initial conditions\n",
    "circuit.V_s = E_L\n",
    "circuit.V_a = E_L\n",
    "circuit.m_Ca = 'm_Ca_inf'\n",
    "circuit.h_Ca = 'h_Ca_inf'\n",
    "circuit.m_K = 'm_K_inf'\n",
    "circuit.m_A = 'm_A_inf'\n",
    "circuit.h_A = 'h_A_inf'\n",
    "circuit.m_proc = 'm_proc_inf'\n",
    "circuit.m_Na = 'm_Na_inf'\n",
    "circuit.h_Na = 'h_Na_inf'\n",
    "circuit.m_Kd = 'm_Kd_inf'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The definition of the synapses and the simulation protocol are identical to the simplified model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eqs_fast = '''\n",
    "g_fast : siemens (constant)\n",
    "I_fast_post = g_fast*(V_s_post - E_syn)/(1+exp(s_fast*(V_fast-V_s_pre))) : amp (summed)\n",
    "'''\n",
    "fast_synapses = Synapses(circuit, circuit, model=eqs_fast)\n",
    "fast_synapses.connect('label_pre != label_post and not (label_pre == PY and label_post == ABPD)')\n",
    "fast_synapses.g_fast['label_pre == ABPD and label_post == LP'] = 0.015*uS\n",
    "fast_synapses.g_fast['label_pre == ABPD and label_post == PY'] = 0.005*uS\n",
    "fast_synapses.g_fast['label_pre == LP and label_post == ABPD'] = 0.01*uS\n",
    "fast_synapses.g_fast['label_pre == LP and label_post == PY']   = 0.02*uS\n",
    "fast_synapses.g_fast['label_pre == PY and label_post == LP']   = 0.005*uS\n",
    "\n",
    "eqs_slow = '''\n",
    "k_2 : 1/second (constant)\n",
    "g_slow : siemens (constant)\n",
    "I_slow_post = g_slow*m_slow*(V_s_post-E_syn) : amp (summed)\n",
    "dm_slow/dt = k_1*(1-m_slow)/(1+exp(s_slow*(V_slow-V_s_pre))) - k_2*m_slow : 1 (clock-driven)\n",
    "'''\n",
    "slow_synapses = Synapses(circuit, circuit, model=eqs_slow, method='exact')\n",
    "slow_synapses.connect('label_pre == ABPD and label_post != ABPD')\n",
    "slow_synapses.g_slow['label_post == LP'] = 0.025*uS\n",
    "slow_synapses.k_2['label_post == LP']    = 0.03/ms\n",
    "slow_synapses.g_slow['label_post == PY'] = 0.015*uS\n",
    "slow_synapses.k_2['label_post == PY']    = 0.008/ms\n",
    "\n",
    "\n",
    "M = StateMonitor(circuit, 'V_s', record=True)\n",
    "spikes = SpikeMonitor(circuit)\n",
    "\n",
    "M.active = False\n",
    "run(init_time, report='text')\n",
    "M.active = True\n",
    "run(observe_time, report='text')\n",
    "M.active = False\n",
    "run(adapt_time, report='text')\n",
    "M.active = True\n",
    "run(observe_time, report='text')\n",
    "device.build(directory='example_1_complex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the same analysis/plotting procedure as in the simplified model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spike_trains = spikes.spike_trains()\n",
    "\n",
    "# Plot\n",
    "from example_1_pyloric_net_plot import do_pyloric_net_plot\n",
    "do_pyloric_net_plot(spike_trains, M.t, M.V_s, 'V_s',\n",
    "                    init_time, observe_time, adapt_time);"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}