Current-Based-Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rate-Based Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remaining tasks:\n",
    "\n",
    "* refactor code for generalising simulations to tasks\n",
    "* ~~parallelising code for faster EA~~\n",
    "* putting code inside DEAP\n",
    "* testing R-BCM learning rule\n",
    "* attempting solving XOR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parameters import *\n",
    "\n",
    "import numba\n",
    "from numba import jit\n",
    "\n",
    "import numpy as np\n",
    "from scipy import special, integrate\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.23389231, -0.23389231, -0.23389231])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Jahr-Stevens formula functions\n",
    "@jit(nopython=True)\n",
    "def J(V):\n",
    "    result = np.empty_like(V)\n",
    "    for i in range(V.shape[0]):\n",
    "        result[i] = 1 + gamma_JahrStevens * np.exp(-beta_JahrStevens * V[i])\n",
    "    return result\n",
    "        \n",
    "@jit(nopython=True)\n",
    "def J_2(V):\n",
    "    result = np.empty_like(V)\n",
    "    J_V = J(V)\n",
    "    for i in range(V.shape[0]):\n",
    "        numerator = J_V[i] + beta_JahrStevens * (V[i] - V_E)*(J_V[i] - 1)\n",
    "        denominator = J_V[i]**2\n",
    "        result[i] = numerator / denominator\n",
    "    return result\n",
    "\n",
    "J_2(np.full(p, V_drive))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Functions to compute steady-state NMDA channels\n",
    "@jit(nopython=True)\n",
    "def my_factorial(n):\n",
    "    if n <= 1.:\n",
    "        return 1.\n",
    "    return my_factorial(n-1) * n\n",
    "\n",
    "@jit(nopython=True)\n",
    "def my_binomial(n, m):\n",
    "    return my_factorial(n) / (my_factorial(m) * my_factorial(n-m))\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def _get_Tn_summand(n, m, nu):\n",
    "    assert n >= 1 and m >= 0\n",
    "    binom_coeff = (-1)**m * my_binomial(n, m)\n",
    "    result = np.empty_like(nu)\n",
    "    for i in range(nu.shape[0]):\n",
    "        numerator = tau_NMDA_rise * (1 + nu[i] * tau_NMDA)\n",
    "        denominator = tau_NMDA_rise * (1 + nu[i] * tau_NMDA) + m * tau_NMDA_decay\n",
    "        result[i] = binom_coeff * numerator / denominator\n",
    "    return result\n",
    "\n",
    "def _get_Tn_summand_scipy(n, m, nu):\n",
    "    assert n >= 1 and m >= 0\n",
    "    binom_coeff = (-1)**m * special.binom(n, m)\n",
    "    numerator = tau_NMDA_rise * (1 + nu * tau_NMDA)\n",
    "    denominator = tau_NMDA_rise * (1 + nu * tau_NMDA) + m * tau_NMDA_decay\n",
    "    return binom_coeff * numerator / denominator\n",
    "\n",
    "@jit(nopython=True)\n",
    "def _get_Tn(n, nu):\n",
    "    assert n >= 1\n",
    "    Tn = np.zeros_like(nu)\n",
    "    for m in np.arange(0, n+1):\n",
    "        Tn += _get_Tn_summand(n, m, nu)\n",
    "    return Tn\n",
    "\n",
    "def _get_Tn_scipy(n, nu):\n",
    "    assert n >= 1\n",
    "    Tn = 0\n",
    "    for m in np.arange(0, n+1):\n",
    "        Tn += _get_Tn_summand_scipy(n, m, nu)\n",
    "    return Tn\n",
    "\n",
    "@jit(nopython=True)\n",
    "def psi(nu, n_truncate=5):\n",
    "    \"\"\"\n",
    "    Computes fraction of open NMDA channels given presynaptic firing rate nu \n",
    "    up until truncation round-off n_truncate.\n",
    "    \"\"\"\n",
    "    coeff = nu * tau_NMDA / (1 + nu * tau_NMDA)\n",
    "    summation = np.zeros_like(nu)\n",
    "    for n in np.arange(1, n_truncate+1):\n",
    "        Tn = _get_Tn(n, nu)\n",
    "        summand_coeff = ((-alpha * tau_NMDA_rise)**n) / my_factorial(n+1)\n",
    "        summation += summand_coeff * Tn\n",
    "    return coeff * (1 + summation / (1 + nu * tau_NMDA))\n",
    "\n",
    "def psi_scipy(nu, n_truncate=5):\n",
    "    \"\"\"\n",
    "    Computes fraction of open NMDA channels given presynaptic firing rate nu \n",
    "    up until truncation round-off n_truncate.\n",
    "    \"\"\"\n",
    "    coeff = nu * tau_NMDA / (1 + nu * tau_NMDA)\n",
    "    summation = 0\n",
    "    for n in np.arange(1, n_truncate+1):\n",
    "        summand_coeff = ((-alpha * tau_NMDA_rise)**n) / special.factorial(n+1)\n",
    "        summation += summand_coeff * _get_Tn_scipy(n, nu)\n",
    "    return coeff * (1 + summation / (1 + nu * tau_NMDA))\n",
    "\n",
    "# psi(rate_ext) == psi_scipy(rate_ext)\n",
    "test_rate = np.full(p, rate_ext)\n",
    "np.allclose(psi(test_rate), psi_scipy(test_rate))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# not used\n",
    "gamma = 0.8  # 0.641 in Wong&Wang\n",
    "@jit(nopython=True)\n",
    "def psi_approx(nu):\n",
    "    \"\"\"Following Wong&Wang 2006, we can approximate psi thusly\"\"\"\n",
    "    product = gamma * nu * tau_NMDA\n",
    "    return product/(1 + product)\n",
    "\n",
    "@jit(nopython=True)\n",
    "def psi_regular(nu):\n",
    "    \"\"\"Gives the steady-state value under regular firing\"\"\"\n",
    "    product = nu * tau_NMDA\n",
    "    return product/(1 + product)\n",
    "\n",
    "# psi_approx(np.full(p, rate_ext)), psi_regular(np.full(p, rate_ext))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Weights and mask\n",
    "pyramidal_mask = np.array([True] * (p+1) + [False])\n",
    "\n",
    "# (W * 100).astype(int).astype(float) / 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def _rate_upperbound_vectorised(\n",
    "    V_SS, sigma, tau_m, tau_rp,\n",
    "):\n",
    "    summand = (V_thr - V_SS)/sigma\n",
    "    summand *= 1 + 0.5 * tau_AMPA/tau_m\n",
    "    summand += 1.03 * np.sqrt(tau_AMPA/tau_m) - 0.5 * tau_AMPA/tau_m\n",
    "    return summand\n",
    "\n",
    "@jit(nopython=True)\n",
    "def _rate_lowerbound_vectorised(\n",
    "    V_SS, sigma,\n",
    "):\n",
    "    return (V_reset - V_SS) / sigma\n",
    "\n",
    "def _siegert_integrand(x):\n",
    "    return np.exp(x**2)*(1+special.erf(x))\n",
    "\n",
    "\n",
    "root_pi = np.sqrt(np.pi)\n",
    "def rate_vectorised(\n",
    "    V_SS, sigma, tau_m, tau_rp\n",
    "):\n",
    "    integration_results = np.empty(V_SS.shape)\n",
    "    UB = _rate_upperbound_vectorised(\n",
    "        V_SS, sigma,\n",
    "        tau_m=tau_m, tau_rp=tau_rp\n",
    "    )\n",
    "    LB = _rate_lowerbound_vectorised(\n",
    "        V_SS, sigma\n",
    "    )\n",
    "    # There would be a dynamic programming solution\n",
    "    # to computing multiple integrals by splitting \n",
    "    # this up into subintervals, integrating over \n",
    "    # them and adding them back together\n",
    "    for i, (lb, ub) in enumerate(zip(LB, UB)):\n",
    "        integral, error = integrate.quad(\n",
    "            _siegert_integrand, lb, ub\n",
    "        )\n",
    "        integration_results[i] = integral\n",
    "    return (tau_rp + tau_m * root_pi * integration_results)**-1\n",
    "\n",
    "# rate_vectorised(np.full(5, V_drive), np.full(5, 2e-3), tau_m_E, tau_rp_E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# @jit(nopython=True)  # faster without jit for now\n",
    "def ds_NMDA_dt(s_NMDA, nu):\n",
    "    psi_nu = psi(nu)\n",
    "    psi_nu[~pyramidal_mask] = 0.0\n",
    "    tau_NMDA_eff = tau_NMDA * (1 - psi_nu)\n",
    "    dsdt = -(s_NMDA - psi_nu) / tau_NMDA_eff\n",
    "    return dsdt\n",
    "\n",
    "# @jit(nopython=True)\n",
    "def dic_noise_dt(\n",
    "    ic_noise, sigma_noise=7e-12#*5/p\n",
    "):\n",
    "    eta = np.random.randn(*ic_noise.shape)\n",
    "    dicdt = (-ic_noise + eta * np.sqrt(tau_AMPA/defaultdt) * sigma_noise) / tau_AMPA\n",
    "    return dicdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Following W&W2006\n",
    "# 1. Treat all firing rates and synaptic variables as dynamic\n",
    "# 2. use tau_r=2ms for the firing rate dynamics\n",
    "# 3. phi(I_syn) = rates(V_SS_new = V_L - I_syn / g_m)\n",
    "\n",
    "# Can't jit rate (yet)\n",
    "# @jit(nopython=True)\n",
    "def phi(\n",
    "    I_syn, g_m, sigma, tau_m, tau_rp\n",
    "):\n",
    "    V_SS = V_L - I_syn/g_m\n",
    "    return rate_vectorised(V_SS, sigma, tau_m, tau_rp)\n",
    "\n",
    "# tau_r = 2e-3  # moved to parameters.py\n",
    "# @jit(nopython=True)\n",
    "def dnu_dt(\n",
    "    nu, I_syn, g_m, sigma, tau_m, tau_rp\n",
    "):\n",
    "    phi_Isyn = phi(\n",
    "        I_syn, g_m, sigma, tau_m, tau_rp\n",
    "    )\n",
    "    deriv = (-nu + phi_Isyn)/tau_r\n",
    "    return deriv\n",
    "\n",
    "@jit(nopython=True)\n",
    "def ds_AMPA_dt(s_AMPA, nu):\n",
    "    deriv = -s_AMPA/tau_AMPA + nu\n",
    "    deriv[~pyramidal_mask] = 0.0\n",
    "    return deriv\n",
    "\n",
    "@jit(nopython=True)\n",
    "def ds_GABA_dt(s_GABA, nu):\n",
    "    deriv = -s_GABA/tau_GABA + nu\n",
    "    deriv[pyramidal_mask] = 0.0\n",
    "    return deriv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BCM-rule from 1994 (Law and Cooper):\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "% y &= \\sigma \\left( \\sum_j w_j x_j \\right)\\\\\n",
    "\\frac{dw_{ij}}{dt} &= \\nu_i(\\nu_i - \\theta_{M,i}) \\nu_j / \\theta_{M,i}\\\\\n",
    "\\theta_M &= \\mathbb{E}[\\nu^2]\n",
    "\\end{split}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "plasticity_mask_source = pyramidal_mask.copy()\n",
    "plasticity_mask_target = np.full_like(plasticity_mask_source, True)\n",
    "# plasticity_mask = np.outer(plasticity_mask_target, plasticity_mask_source)\n",
    "\n",
    "tau_theta = 0.1\n",
    "@jit(nopython=True)\n",
    "def dtheta_BCM_dt(theta, nu):\n",
    "    \"\"\"\n",
    "    tau_theta * dtheta/dt = -theta + nu**2\n",
    "    \n",
    "    Units of theta are technically different to those of nu. \n",
    "    There should be a constant to fix this.\n",
    "    \"\"\"\n",
    "    dtheta_dt = (-theta + nu**2)/tau_theta\n",
    "    return dtheta_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau_W = 10.\n",
    "@jit(nopython=True)\n",
    "def dW_dt_BCM(W, nu, theta):\n",
    "    dW_dt = (np.outer(nu * (nu-theta), nu) / theta.reshape(-1, 1)) / tau_W\n",
    "    # BCM has a positive feedback loop on recurrent weights\n",
    "#     for i in range(dW_dt.shape[0]):\n",
    "#         dW_dt[i,i] = 0.0\n",
    "    for j, b in enumerate(plasticity_mask_source):\n",
    "        if not b:\n",
    "            dW_dt[:, j] = 0.0\n",
    "    for i, b in enumerate(plasticity_mask_target):\n",
    "        if not b:\n",
    "            dW_dt[i, :] = 0.0\n",
    "#     dW_dt[~plasticity_mask] = 0.0\n",
    "    return dW_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuSklEQVR4nO3deXgV5fn/8fdNQhJC2JewRRJkE1AQAsGlCooVqMi3WhUtKIsiClV/rVawVlutLa3aKlVBRGRT0AoUtAgKNbggEnbZAiRsYd8hIYQs9++PM9pjEsgJWWZycr+u61ycM/M8cz5zSM6dmWcWUVWMMcYYf1XcDmCMMcZ7rDgYY4wpwIqDMcaYAqw4GGOMKcCKgzHGmAJC3Q5QGurXr6+xsbEX3T8jI4Pq1auXXqBSYrmKx3IVj+UqnmDMtWrVqiOq2qDQmapa4R9dunTRkvj8889L1L+sWK7isVzFY7mKJxhzASv1PN+rtlvJGGNMAVYcjDHGFGDFwRhjTAFWHIwxxhRgxcEYY0wBVhyMMcYUYMXBGGNMAVYcjDHGFGDFwRhjTAFWHDwiKSmJpKQkt2MYYwwQJNdWCgZPPPEEAImJie4GMcYYrDh4xmuvveZ2BGOM+YEVB4/o0KGD2xGMMeYHNubgEcuWLWPZsmVuxzDGGMC2HDzjqaeeAmzMwRjjDVYcPOLNN990O4IxxvzAioNHtGnTxu0IxpgK5p2vd5BxJJceZbBsG3PwiKVLl7J06VK3YxhjKojPkw/x3Meb+Hpfdpks37YcPOLZZ58FbMzBGFO0HUcyeGTmGi5rVJPB7XPK5D2sOHjE5MmT3Y5gjKkA0rNyeGDaSkKrCG8O6kLK+hVl8j5WHDyiRYsWbkcwxnhcXp7y6/fXsuNIBtOHdiOmbiQpZfReNubgEYsXL2bx4sVuxzDGeNi4/27j000HearvZVzdsn6ZvpdtOXjEn/70JwB69erlchJjjBct3HCAVxZv4/bOzRh6TWyZv19AWw4i0ltEkkVku4iMLmS+iMg4Z/56EelcVF8Red5pu1ZEPhWRJn7zxjjtk0Xk5pKuZEUwffp0pk+f7nYMY4wHbTlwil9/sJaOMbV54ecdEJEyf88ii4OIhACvA32AdsDdItIuX7M+QCvnMRwYH0DfF1X1ClXtBHwMPOP0aQcMANoDvYE3nOUEtZiYGGJiYtyOYYzxmOMZ5xg+bRVR4aFMHNSFiKrl83UYyJZDN2C7qqaq6jlgFtA/X5v+wDT1WQ7UFpHGF+qrqqf8+lcH1G9Zs1Q1S1V3ANud5QS1hQsXsnDhQrdjGGM8JDs3j4ffXc2Bk2eZMKgL0TUjyu29AxlzaArs8XudBiQE0KZpUX1F5AXgXuAk0NNvWcsLWdaPiMhwfFspREdHl+j8gPT0dNfPLxg92rfHLSLif//5XshVGMtVPJareCzX/0zflMU3u3N44PIwTqWuIzG1/HIFUhwK27mlAba5YF9V/R3wOxEZA4wCng3w/VDVicBEgPj4eO3Ro0dh2QOSmJhISfqXhu+3Gho1avTDNC/kKozlKh7LVTyWy+fdb3exZPcGhl/Xgqf6XlbuuQLZrZQG+O8MbwbsC7BNIH0B3gNuL8b7BZ1GjRr9qDAYYyqvb1KO8uy8jVzfugFP9m7rSoZAikMS0EpE4kQkDN9g8fx8beYD9zpHLXUHTqrq/gv1FZFWfv1vBbb4LWuAiISLSBy+Qe6yOQXQQz766CM++ugjt2MYY1y2++gZHnp3Fc3rRTLu7isJqVL2RyYVpsjdSqqaIyKjgEVACDBZVTeKyAhn/gRgAdAX3+DxGWDIhfo6ix4rIm2APGAX8P3yNorIB8AmIAcYqaq5pbXCXvXyyy8D0K9fP5eTGGPccvpsNsOmJgHw9n1dqVWtqmtZAjoJTlUX4CsA/tMm+D1XYGSgfZ3ptxfS/Pt5LwAvBJItWHz44YduRzDGuCg3T3lk5hp2HMlg2rBuxNav7moeO0PaI+rXL9tT4Y0x3vaXBZv5PPkwf/q/Dlx9qfvfB3ZtJY+YM2cOc+bMcTuGMcYFM1fsZtJXOxh8dSwDuzd3Ow5gWw6eMW7cOABuu+02l5MYY8rTspQj/P7fG7iudQOe/tn5D1ktb1YcPGLevHluRzDGlLMdRzJ4aMZqYutX57V7riQ0xDs7c6w4eEStWrXcjmCMKUcnzpxj6JQkQqoIb98XT80I945MKox3ylQl9/777/P++++7HcMYUw7O5eTx4PRV7D2eycRBXWhez90jkwpjWw4eMX78eADuuusul5MYY8qSqvK7ud/x7Y5jvHJXJ+Jj67odqVBWHDxiwYICp4IYY4LQG4kp/GtVGo/e2Ir/u7LANUU9w4qDR0RGRrodwRhTxj5ev48XFyXTv1MTHuvVqugOLrIxB4+YMWMGM2bMcDuGMaaMrNp1nF9/sI6usXX46+1XlMvd3ErCthw8YtKkSQAMHDjQ5STGmNK2++gZHpi2kia1InhzUHy53c2tJKw4eMRnn33mdgRjTBk4nnGOwe+sIE+VyYO7Urd6mNuRAmLFwSOqVvXWMc7GmJI7m53L8OkrSTuRybv3J9CiQZTbkQJmYw4eMWXKFKZMmeJ2DGNMKcnLU574cD1JO4/z8h0d6erRQ1bPx4qDR1hxMCa4vPhpMh+t28eTvdvSr2MTt+MUm+1W8ggv3lDdGHNxZizfxfjEFO5JuIQR17dwO85FsS0HY4wpRYs3HeSZeRu4sW1Dnru1vecPWT0fKw4e8dZbb/HWW2+5HcMYUwLr9pzgVzPX0KFpLf7psausFlfFTR5k7MJ7xlRsO49kMHRKEvVrhPH2fV2JDKvYe+0DKg4i0ltEkkVku4iMLmS+iMg4Z/56EelcVF8ReVFEtjjt54pIbWd6rIhkisha5zEh//sFo8WLF7N48WK3YxhjLsKR9Czuc85lmDqkGw1qhLsdqcSKLA4iEgK8DvQB2gF3i0i7fM36AK2cx3BgfAB9PwM6qOoVwFZgjN/yUlS1k/MYcbErZ4wxZe3MuRyGTUni4KmzvD24a4U6l+FCAtly6AZsV9VUVT0HzAL652vTH5imPsuB2iLS+EJ9VfVTVc1x+i8HmpXC+lRYb7zxBm+88YbbMYwxxZCdm8fId1fz3d6T/PPuznS+pI7bkUpNIDvFmgJ7/F6nAQkBtGkaYF+AoYD/Dvc4EVkDnAKeVtUv83cQkeH4tlKIjo4u0aGg6enprh9KOnXqVADatfvfRpkXchXGchWP5SqeipJLVZn03Tm+3pfD4PZhVD20mcRDm13PVWpU9YIP4A5gkt/rQcA/87X5D3Ct3+slQJcA+/4OmAuI8zocqOc874KvuNS8UMYuXbpoSXz++ecl6l9WLFfxWK7isVzFkz/X2E82a/MnP9ZXPtvqTiBHST4vYKWe53s1kC2HNCDG73UzYF+AbcIu1FdE7gNuAW50gqKqWUCW83yViKQArYGVAWQ1xpgy987XO344ye2RG1u6HadMBDLmkAS0EpE4EQkDBgDz87WZD9zrHLXUHTipqvsv1FdEegNPAreq6pnvFyQiDZyBbESkBb5B7tQSrWUF8Oqrr/Lqq6+6HcMYU4R5a/fyx4828dN20Tzfv0OFPcmtKEVuOahqjoiMAhYBIcBkVd0oIiOc+ROABUBfYDtwBhhyob7Ool/DtwvpM+fDXa6+I5OuA54TkRwgFxihqsdKa4W9asmSJQA8+uijLicxxpzPF1sP8/i/1pEQV5dxd19JSJXgLAwQ4LWVVHUBvgLgP22C33MFRgba15le6LaYqs4GZgeSK5jMn59/Y8wY4yWpJ3J56b+raNmwBm/dVzFu2FMSdoa0McYUYfuh0/x91VnqRYUxdUhXakYE//1XrDh4xEsvvcRLL73kdgxjTD57T2Qy6O0VVBFh+tAEGtaMcDtSuajYF/8IIt98843bEYwx+RxNz2LQ29+SnpXDE/HhxNav7nakcmPFwSNmz650wyzGeNrps9kMfieJvcczmXF/Ahk717sdqVzZbiVjjMnnbHYuw6auZPP+U4wf2LnC3eKzNNiWg0eMHTsWgNGjC1z01hhTjrJz83j43dUk7TzGK3d14oa20W5HcoUVB49Yu3at2xGMqfTy8pTH/7WO/245xAs/70D/Tk3djuQaKw4eMWvWLLcjGFOpqSpPz9vAvLX7+G3vNvwyobnbkVxlYw7GmEpPVfnzgs289+1uHu5xKQ/3CM7rJRWHFQePeP7553n++efdjmFMpTRuyXbe+nIH917VnCdubuN2HE+w3UoekZyc7HYEYyqlSV+m8o/FW7mtc1P+0K990F5Ir7isOHjEjBkz3I5gTKUzY/ku/vSfzfzs8sb87fYrqBLEF9IrLtutZIyplGavSuPpf2/ghrYN+cddnQgNsa9Df/ZpeMQzzzzDM88843YMYyqFj9bt44kP13Fty/q88cvOhIXaV2F+tlvJI/bs2VN0I2NMiS3ccIDH3l9LfPO6TLy3S9BfevtiWXHwiHfeecftCMYEvf9uOcivZq7mima1mDykK5Fh9hV4PrYtZYypFL7YepgRM1bTtlFNpgzpRlS4FYYLseLgEWPGjGHMmDFuxzAmKH29/QgPTFvJpQ2imDa0G7WqBf/NekrKSqdHHD161O0IxgSlb1KOMmxqEnH1q/Pu/QnUqR7mdqQKwYqDR0ycONHtCMYEnW9TjzJ0ShIxdSKZcX8Cda0wBCyg3Uoi0ltEkkVku4gUuKa0+Ixz5q8Xkc5F9RWRF0Vki9N+rojU9ps3xmmfLCI3l3AdjTGV0IodxxgyJYkmtSN494EE6keFux2pQimyOIhICPA60AdoB9wtIu3yNesDtHIew4HxAfT9DOigqlcAW4ExTp92wACgPdAbeMNZTlB7/PHHefzxx92OYUxQWLHjGIPfWUHjWhHMHN6dhjUqx32fS1MgWw7dgO2qmqqq54BZQP98bfoD09RnOVBbRBpfqK+qfqqqOU7/5UAzv2XNUtUsVd0BbHeWE9QyMzPJzMx0O4YxFV7STl9haFQrgpkPWGG4WKKqF24g8gugt6re77weBCSo6ii/Nh8DY1X1K+f1EuBJILaovs70j4D3VXWGiLwGLFfVGc68t4FPVPXDfH2G49tKITo6uktJ7oeQnp5OVFTURfcvK5areCxX8QRjruRjufx91VnqRAiju0ZQO6L0DsgMxs+rZ8+eq1Q1vrB5gQxIF3YlqvwV5XxtiuwrIr8DcoB3i/F+qOpEYCJAfHy89ujRo5BugUlMTKQk/cuK5Soey1U8wZbrm5SjvLIkiZh6Ubz3QEKpbzEE2+dVlECKQxoQ4/e6GbAvwDZhF+orIvcBtwA36v82YQJ5v6Dz2GOPAfDKK6+4msOYimjZ9iMMneo7Kum9B7rToIYNPpdUINtcSUArEYkTkTB8g8Xz87WZD9zrHLXUHTipqvsv1FdEeuPb9XSrqp7Jt6wBIhIuInH4BrlXlGAdjTFBbOnWwwyZkkTzutWZOdwKQ2kpcstBVXNEZBSwCAgBJqvqRhEZ4cyfACwA+uIbPD4DDLlQX2fRrwHhwGfOzTWWq+oIZ9kfAJvw7W4aqaq5pbbGHmVbDMYU35LNB3loxmpaNoyy8xhKWUAnwanqAnwFwH/aBL/nCowMtK8z/bw3aVXVF4AXAslmjKmcFm08wKj3VnNZ45pMG9qN2pFWGEqTXVvJI0aOHMnIkYXWV2NMPvPX7ePhd1fToWktpg9LsMJQBuzyGR5RrVo1tyMYUyF8sHIPT85eT9fYukwe3NWurlpG7FP1iJdeesntCMZ43vRvdvL7eRv5Sav6TBwUT7WwoL94gmusOBhjKoQJS1MY+8kWel3WkNfu6Wx3cCtjVhw8Yvjw4YBdndWY/FSVlz/dymufb6dfxyb8/c6OVA2x4dKyZsXBI+rVq+d2BGM8Jy9Pee7jTUxZtpMBXWN44eeXE1KlsIsomNJmxcEj/vKXv7gdwRhPycnN48nZ3zF7dRrDro3j6Z9dhnNOlCkHVhyMMZ6TlZPLr95bw6ebDvLrm1rzqxtaWmEoZ1YcPGLIkCEAvPPOOy4nMcZdmTnK0ClJfL39KM/2a8eQa+LcjlQpWXHwiJiYmKIbGRPkjmWc429JZ9l9OpOX7ujIL7o0K7qTKRNWHDziueeeczuCMa7afzKTQW+vYM/pPCYMjOemdtFuR6rUrDgYY1y3/VA69779LafP5vB4fIQVBg+wg4U9YuDAgQwcONDtGMaUuzW7j3PHhGWcy1VmDu9O27p2cpsX2JaDR7Rp08btCMaUu8TkQzw0YzUNaoQzfVg3mterTuI2t1MZsOLgGb///e/djmBMuZq9Ko0nZ6+ndXQNpgztWuq39TQlY8XBGFOuVJXxS1P428JkrmlZjwkDu1AjoqrbsUw+Vhw8YsCAAQDMmjXL5STGlJ3cPOV553IYt3Zswkt3dCQs1IY+vciKg0d06tTJ7QjGlKmz2bk8OmsNizYe5P5r43iq72VUseskeZYVB48YPXq02xGMKTPHMs5x/9Qk1uw5we9vacewa+2sZ68LaHtORHqLSLKIbBeRAt9i4jPOmb9eRDoX1VdE7hCRjSKSJyLxftNjRSRTRNY6jwn5388YU3HsPnqGX4xfxoZ9p3jt7s5WGCqIIrccRCQEeB24CUgDkkRkvqpu8mvWB2jlPBKA8UBCEX03ALcBbxbytimq2umi16oCuv322wGYPXu2y0mMKT1rdh/n/qkryVVlxrAEusXVdTuSCVAgu5W6AdtVNRVARGYB/QH/4tAfmKaqCiwXkdoi0hiIPV9fVd3sTCutdanQrrrqKrcjGFOqFm08wKOz1tCwRgTvDOnKpQ2i3I5kikF83+cXaCDyC6C3qt7vvB4EJKjqKL82HwNjVfUr5/US4El8xaGovonA46q60nkdC2wEtgKngKdV9ctCcg0HhgNER0d3KclRPunp6URFee8H13IVj+UqnrLKpaos3JnDB8nniKtVhcc6R1AzPPA/Aivb51VSJcnVs2fPVaoaX9i8QLYcCvtfzV9RztcmkL757QcuUdWjItIF+LeItFfVUz9aiOpEYCJAfHy89ujRo4jFnl9iYiIl6V9WLFfxWK7iKYtcObl5PDt/I+8n76bv5Y14+Y5OVAsr3uUwKtPnVRrKKlcgxSEN8L+edDNgX4BtwgLo+yOqmgVkOc9XiUgK0BpYGUDWCuvWW28FYP78+S4nMebinDqbzaj31vDF1sM81ONSnvhpGztUtQILpDgkAa1EJA7YCwwA7snXZj4wyhlTSABOqup+ETkcQN8fEZEGwDFVzRWRFvgGuVOLs1IV0Y033uh2BGMu2p5jZxg6JYkdRzL46+2Xc1fXS9yOZEqoyOKgqjkiMgpYBIQAk1V1o4iMcOZPABYAfYHtwBlgyIX6AojIz4F/Ag2A/4jIWlW9GbgOeE5EcoBcYISqHivNlfaiRx991O0IxlyUpJ3HeHD6KnLzlGnDunH1pfXdjmRKQUAnwanqAnwFwH/aBL/nCowMtK8zfS4wt5DpswE7ntOYCuDDVWk8Nec7mtapxtv3xdPCjkgKGnaGtEf06dMHgE8++cTlJMYULTdP+dvCLbz5RSpXX1qPN37ZmdqRYW7HMqXIioNH9OvXz+0IxgTk9NlsHpu1liVbDjGoe3Oe6deOqiF28bxgY8XBIx5++GG3IxhTpB1HMnhg2kp2HMng+f7tGXRVrNuRTBmx4mCMCcgXWw8z6r3VhFQRptvAc9Cz4uARvXr1AmDx4sUuJzHmx1SVSV/u4C+fbKZVwxpMui+emLqRbscyZcyKg0fcddddbkcwpoDMc7mMnrOeeWv3cXP7aF6+sxNR4fa1URnY/7JHPPDAA25HMOZH0o6f4cHpq9i0/xSP/7Q1D/doaWc8VyJWHIwxBXy57TCPzFxDTp7y9n3x3NA22u1IppxZcfCI7y+clZiY6GoOU7mpKhOWpvLioi20aliDNwd1IbZ+dbdjGRdYcfCIwYMHux3BVHKnzmbzxL/WsWjjQW65ojF/vf0Kqtv4QqVl//MeYcXBuCn5wGlGzFjF7mNnePpnlzHs2ji7EVclZ8XBI7KzswGoWrWqy0lMZfPvNXsZM+c7oiJCmflAd7uVpwGsOHjGTTfdBNiYgyk/53KVp+Z+x3vf7qZbbF1eu+dKGtaMcDuW8QgrDh5x//33ux3BVCK7j57hhW/PsuvUbkZcfymP/7Q1oXZ9JOPHioNHDBw40O0IppL45Lv9/PbD9eTm5jHp3nh6tbPDVE1BVhw84syZMwBERtplCUzZyMrJ5c//2czUb3bRMaY2A+OyrDCY87LtSI/o27cvffv2dTuGCVI7jmRw+/hlTP1mF/dfG8e/HryKBpH262/Oz7YcPOKhhx5yO4IJUnPXpPH03A1UDa3CxEFd+Gn7Rm5HMhWAFQePsAvvmdKWkZXDM/M2Mnt1Gl1j6/DqgCtpUrua27FMBRHQdqWI9BaRZBHZLiKjC5kvIjLOmb9eRDoX1VdE7hCRjSKSJyLx+ZY3xmmfLCI3l2QFK4qTJ09y8uRJt2OYIPFd2klu+edXzF2TxiM3tGTmA92tMJhiKXLLQURCgNeBm4A0IElE5qvqJr9mfYBWziMBGA8kFNF3A3Ab8Ga+92sHDADaA02AxSLSWlVzS7SmHte/f3/AznMwJZOXp0z6KpUXFyVTPyqcmQ90J6FFPbdjmQookN1K3YDtqpoKICKzgP6Af3HoD0xTVQWWi0htEWkMxJ6vr6pudqblf7/+wCxVzQJ2iMh2J8M3F7eKFcMjjzzidgRTwe0/mclvPljHspSj9G7fiLG3X07tyDC3Y5kKKpDi0BTY4/c6Dd/WQVFtmgbYt7D3W17IsoLabbfd5nYEU4Et+G4/Y+Z8x7mcPP56++XcGR9j10YyJRJIcSjsJ0wDbBNI34t5P0RkODAcIDo6ukS7Y9LT013fnfP9eEOtWrV+mOaFXIWxXMVTlrnOZCvvbj7H1/tyiKtVhQfjw4nOSGXp0lRXc5WE5SqessoVSHFIA2L8XjcD9gXYJiyAvhfzfqjqRGAiQHx8vH5/P4SLkZiYSEn6l4bC7ufghVyFsVzFU1a5lqce5XcfrOPAqVweubEVv7qhJVWLcQmMyvZ5lVRlyxVIcUgCWolIHLAX32DxPfnazAdGOWMKCcBJVd0vIocD6JvffOA9Efk7vgHpVsCKQFeoovrNb37jdgRTQZzNzuXFRclM/noHsfWq8+GIq7jykjpuxzJBpsjioKo5IjIKWASEAJNVdaOIjHDmTwAWAH2B7cAZYMiF+gKIyM+BfwINgP+IyFpVvdlZ9gf4BrxzgJHBfqQSQL9+/dyOYCqAdXtO8OsP1pJyOIN7r2rO6D5tiQyz05VM6Qvop0pVF+ArAP7TJvg9V2BkoH2d6XOBuefp8wLwQiDZgsWBAwcAaNTIzl41BWXl5DJuyTYmLE2lYY1wZgxL4NpW9d2OZYKY/cnhEQMGDADsPAdT0Pq0Ezzxr/UkHzzNHV2a8fQt7ahVzW4KZcqWFQePGD26wInnppI7m+3bWnjzi1TqR4XxzuCu9Gzb0O1YppKw4uARvXv3djuC8ZBVu47z2w/XkXI4gzvjm/G7n9nWgilfVhw8Ys8e37mCMTExRbQ0wSwjK4cXFyUz9ZudNKlVjWlDu3Fd6wZuxzKVkBUHjxg0aBBgYw6V2efJh3h67gb2ncxkUPfm/LZ3W6LC7VfUuMN+8jzi6aefdjuCccnh01k8//Em5q/bR8uGUXw44iq6NK/rdixTyVlx8IhevXq5HcGUs7w85YOVe/jzgs2czc7jsV6teKjHpYSHhrgdzRgrDl6Rmuq7Fk6LFi1cTmLKw5YDp3h67gZW7jpOQlxdXvj55bRsGOV2LGN+YMXBI4YOHQrYmEOwy8jKYdySbUz6agc1I0L52+1XcEd8M7uCqvEcKw4e8cc//tHtCKYMqSqfbDjA8x9vYv/Js9wVH8PoPm2pU93ut2C8yYqDR1x//fVuRzBlZH96HvdOXsGX247QrnFNXrvnShtwNp5nxcEjkpOTAWjTpo3LSUxpOX02m3/+dztvf51JZHg2f7y1Pb9MuITQYlxW2xi3WHHwiAcffBCwMYdgkJenzF2zl7ELt3AkPYufNA3l74N7UD8q3O1oxgTMioNH/PnPf3Y7gikFq3Yd47mPNrEu7SSdYmoz6d54jqestcJgKhwrDh5x9dVXux3BlEDa8TP8bWEy89ftI7pmOP+4qyP9OzalShUhMcXtdMYUnxUHj9iwYQMAHTp0cDmJKY7TZ7MZn5jCpK92IMCvbmjJiOsvpbpd9sJUcPYT7BGjRo0CbMyhosjOzWPmit28ungbRzPOcduVTXn85jY0qV3N7WjGlAorDh7x4osvuh3BBEBVWbTxAH9bmEzqkQwS4uoyue9ldIyp7XY0Y0qVFQeP6Nq1q9sRTBGWpRzhrwuTWbfnBK0aRjF5cDw92zS0s5tNULLi4BFr164FoFOnTq7mMAV9l3aSFz9N5outh2lSK4K//eIKbu/cjJAqVhRM8AqoOIhIb+BVIASYpKpj880XZ35f4AwwWFVXX6iviNQF3gdigZ3Anap6XERigc1AsrP45ao64uJXsWJ47LHHABtz8JJtB0/z8qdbWbjxALUjq/JU37bce1UsEVXtqqkm+BVZHEQkBHgduAlIA5JEZL6qbvJr1gdo5TwSgPFAQhF9RwNLVHWsiIx2Xj/pLC9FVTuVxgpWFK+88orbEYwj5XA645ZsY/66fVQPC+WxXq0Ydm0cNSLsNp2m8ghky6EbsF1VUwFEZBbQH/AvDv2BaaqqwHIRqS0ijfFtFZyvb3+gh9N/KpDI/4pDpWO7k9yXejid1z7fzr/X7CU8NIQR11/K8J+0sIvjmUpJfN/nF2gg8gugt6re77weBCSo6ii/Nh8DY1X1K+f1Enxf9LHn6ysiJ1S1tt8yjqtqHWe30kZgK3AKeFpVvywk13BgOEB0dHSXWbNmXdwnAKSnpxMV5e619Lds2QJA27Ztf5jmhVyFCbZc+9Lz+CjlHMv351K1CtxwSSh948KoGV46YwrB9nmVNctVPCXJ1bNnz1WqGl/YvEC2HAr7DclfUc7XJpC++e0HLlHVoyLSBfi3iLRX1VM/WojqRGAiQHx8vPbo0aOIxZ5fYmIiJelfGv7whz/8kOV7XshVmGDJtXHfSd74PIUFG/YTERrC8OtacP9PWtCgRule6iJYPq/yYrmKp6xyBVIc0oAYv9fNgH0Btgm7QN+DItJYVfc7u6AOAahqFpDlPF8lIilAa2BlQGtUQb322mtuR6g0Vu48xvjEFJZsOURUeCgPXX8pw66No55d/8iYHwRSHJKAViISB+wFBgD35GszHxjljCkkACedL/3DF+g7H7gPGOv8Ow9ARBoAx1Q1V0Ra4BvkTi3BOlYIdtmMspWXpyRuPcT4xBSSdh6nTmRVfn1Ta+67KpZakTbQbEx+RRYHVc0RkVHAInyHo05W1Y0iMsKZPwFYgO8w1u34DmUdcqG+zqLHAh+IyDBgN3CHM/064DkRyQFygRGqeqxU1tbDli1bBtgF+EpbVk4u89bs460vU9l2KJ2mtavxbL923NU1hsgwO83HmPMJ6LdDVRfgKwD+0yb4PVdgZKB9nelHgRsLmT4bmB1IrmDy1FNPAXaeQ2k5lnGO977dxdRvdnH4dBaXNa7J3+/sSL+OTahqN9sxpkj2p5NHvPnmm25HCApbD57mna93MGf1XrJy8riudQP+cWcLrmlZzy5zYUwxWHHwCLs96MXLzVOWbD7IP1ZksnnhF4SHVuG2zk0Zek0craJruB3PmArJioNHLF26FIDrr7/e5SQVx5H0LN5P2sN73+5m74lM6kYIv+3dhgFdL6GunbhmTIlYcfCIZ599FrAxh6KoKstTj/Heit0s3LCf7Fzl6kvr8fTPLiPs8BZu7NHS7YjGBAUrDh4xefJktyN42pH0LOasTmPWij2kHsmgZkQov0xozsDuzWnZ0Hd2aGJichFLMcYEyoqDR7Ro0cLtCJ6Tk5vHF9sO837SHpZsPkROntKleR1e7tmSn13R2K6OakwZsuLgEYsXLwagV69eLidx35YDp5i9Ko25a/ZxJD2LetXDGHJNLHd1jaFlQxtgNqY8WHHwiD/96U9A5S0OB06eZf66vcxZvZctB04TWkW4oW1DftGlGT3aNCQs1M5NMKY8WXHwiOnTp7sdodwdyzjHJxv2M3/tPlbsPIYqdIqpzR9vbc8tVzS2ax0Z4yIrDh4RExNTdKMgcDzjHIs2HuA/3+1nWcpRcvOUFg2q8+iNrbi1YxNaNPDeJZGNqYysOHjEwoULAejdu7fLSUrfwVNn+XTjARZuPMDy1GPk5imX1I1k+HUtuOWKxrRrXNPOXjbGY6w4eMTYsb7bcgdDcVBVth5M57NNB/hs8yHW7TkBQIsG1RlxfQv6dGhM+yZWEIzxMisOHlGSO9l5Qea5XL5JPcJ/txzi8y2H2XsiE4COMbV5/Ket6d2hkR1pZEwFYsXBIxo1auR2hGJRVZIPnubLrUf4Ytthvt1xjHM5eUSGhXBNy/qMuqElN7RtSHTNCLejGmMughUHj/joo48A6Nevn8tJCqeq7DmWyTepR/j3urM8/tUSjqRnAdCqYRT3dm/Oda0bkNCiLuGhdnKaMRWdFQePePnllwHvFAdVJfVIBit2HCNpxzG+3XHsh11FNcOEnu3qcU3L+lzbsj5NaldzOa0xprRZcfCIDz/80NX3z8jK4bu9J1m9+zirdx1n1a7jHD+TDUD9qDC6xdXlwetbcFWLeqRtWknPnle6mtcYU7asOHhE/fr1y+29zmbnsnn/KTbsPcn6NN9j26HT5KlvfosG1el1WTRdmtehW1xd4upX/9GRRXs321FGxgQ7Kw4eMWfOHABuu+22UlumqrLv5Fm2HjhN8sHTbN5/ik37TpFyOP2HQlC3ehiXN61F7w6N6BRTm44xte1eCMaYwIqDiPQGXgVCgEmqOjbffHHm9wXOAINVdfWF+opIXeB9IBbYCdypqsedeWOAYUAu8IiqLirRWlYA48aNAy6uOGSey2XXsQx2HM5gx9EMth9KJ+VQOimHM0jPyvmhXeNaEbRvUpM+HRrRrkktLm9Wiya1Iux8A2NMAUUWBxEJAV4HbgLSgCQRma+qm/ya9QFaOY8EYDyQUETf0cASVR0rIqOd10+KSDtgANAeaAIsFpHWqppbOqvsTfPmzSt0uqpyKjOHA6fOcuDUWfadyGTv8Uz2nshkz7Ez7D52hkOns37UJ7pmOC0bRnF756a0jK5B20Y1aB1dg1rVqpbHqhhjgkAgWw7dgO2qmgogIrOA/oB/cegPTFNVBZaLSG0RaYxvq+B8ffsDPZz+U4FE4Eln+ixVzQJ2iMh2J8M3F7+ahTt8OovZq9M4sCebzO/2U6taVWpEVCUqIpTq4SFEhoVSrWoIIVVK/pe1qpKVk0fmuVzSs3JIz8rh9NkcTmZmczIzmxNnznH8zDmOZezmWEYWR9LPsfvQGU4vWcjZ7LwfLSukitCoZgQxdavRo00DYupEElu/OnH1q9O8XiQ1IqwIGGNKJpDi0BTY4/c6Dd/WQVFtmhbRN1pV9wOo6n4Raei3rOWFLOtHRGQ4MBwgOjr6om6vue14LmO/PQvAlI2rz9uuahXfI7SKULUKhFSBKgIhAv5lQ4E89T1yFXLyIDtPycmDc7m++RciKDXChKgwoVaY0Dwqj3rVq1InPJQ6EUKdCKFehFA7XJyCleU8jsNxOHIcjmwr9sdQbOnp6Z68nanlKh7LVTyVLVcgxaGwP5vzf8+dr00gfS/m/VDVicBEgPj4eO3Ro0cRiy3oelXu6ZvLwv9+QbtO8Zw4k03GD3/VZ5OZncuZc7lknsslKyePc7l5ZGXnkZuXR06ekpOrqF+0KiK+RxWhahUhLLSK7xFShciwECLCQqhWNYSo8FDfIyKUWtWqUqtaVQbeeRtVcrNY6vefnJiYyMWsV1mzXMVjuYrHchVPWeUKpDikAf7Xk24G7AuwTdgF+h4UkcbOVkNj4FAx3q9UiAjVw0OpV60KlzWuWRZvEbBFH8119f2NMcZfILfXSgJaiUiciIThGyyen6/NfOBe8ekOnHR2GV2o73zgPuf5fcA8v+kDRCRcROLwDXKvuMj1qzAiIyOJjIx0O4YxxgABbDmoao6IjAIW4TscdbKqbhSREc78CcACfIexbsd3KOuQC/V1Fj0W+EBEhgG7gTucPhtF5AN8g9Y5wMhgP1IJYMaMGQAMHDjQ5STGGBPgeQ6qugBfAfCfNsHvuQIjA+3rTD8K3HiePi8ALwSSLVhMmjQJsOJgjPEGO0PaIz777DO3IxhjzA+sOHhE1ap2boIxxjsCGZA25WDKlClMmTLF7RjGGANYcfAMKw7GGC+x3Uoe4cUzL40xlZdtORhjjCnAioMxxpgCrDgYY4wpwIqDMcaYAqw4GGOMKcCKgzHGmAKsOBhjjCnAioMxxpgCrDgYY4wpQHxX267YROQwsKsEi6gPHCmlOKXJchWP5Soey1U8wZiruao2KGxGUBSHkhKRlaoa73aO/CxX8Viu4rFcxVPZctluJWOMMQVYcTDGGFOAFQefiW4HOA/LVTyWq3gsV/FUqlw25mCMMaYA23IwxhhTgBUHY4wxBVTq4iAivUUkWUS2i8hoF3NMFpFDIrLBb1pdEflMRLY5/9ZxIVeMiHwuIptFZKOIPOqFbCISISIrRGSdk+uPXsjlly9ERNaIyMdeySUiO0XkOxFZKyIrPZSrtoh8KCJbnJ+zq9zOJSJtnM/p+8cpEXnM7VxOtv/n/MxvEJGZzu9CmeSqtMVBREKA14E+QDvgbhFp51KcKUDvfNNGA0tUtRWwxHld3nKA36jqZUB3YKTzGbmdLQu4QVU7Ap2A3iLS3QO5vvcosNnvtVdy9VTVTn7HxHsh16vAQlVtC3TE97m5mktVk53PqRPQBTgDzHU7l4g0BR4B4lW1AxACDCizXKpaKR/AVcAiv9djgDEu5okFNvi9TgYaO88bA8ke+MzmATd5KRsQCawGEryQC2jm/ILeAHzslf9LYCdQP980V3MBNYEdOAfGeCVXviw/Bb72Qi6gKbAHqAuEAh87+cokV6XdcuB/H/T30pxpXhGtqvsBnH8buhlGRGKBK4Fv8UA2Z9fNWuAQ8JmqeiIX8ArwWyDPb5oXcinwqYisEpHhHsnVAjgMvOPshpskItU9kMvfAGCm89zVXKq6F3gJ2A3sB06q6qdllasyFwcpZJod11sIEYkCZgOPqeopt/MAqGqu+jb7mwHdRKSDy5EQkVuAQ6q6yu0shbhGVTvj2406UkSuczsQvr9+OwPjVfVKIAP3drkVICJhwK3Av9zOAuCMJfQH4oAmQHURGVhW71eZi0MaEOP3uhmwz6UshTkoIo0BnH8PuRFCRKriKwzvquocL2UDUNUTQCK+MRu3c10D3CoiO4FZwA0iMsMDuVDVfc6/h/DtP+/mgVxpQJqz1QfwIb5i4Xau7/UBVqvqQee127l6ATtU9bCqZgNzgKvLKldlLg5JQCsRiXP+QhgAzHc5k7/5wH3O8/vw7e8vVyIiwNvAZlX9u1eyiUgDEantPK+G75dmi9u5VHWMqjZT1Vh8P0//VdWBbucSkeoiUuP75/j2U29wO5eqHgD2iEgbZ9KNwCa3c/m5m//tUgL3c+0GuotIpPO7eSO+AfyyyeXWQI8XHkBfYCuQAvzOxRwz8e1DzMb319QwoB6+gc1tzr91Xch1Lb5dbeuBtc6jr9vZgCuANU6uDcAzznTXPzO/jD3434C0259XC2Cd89j4/c+627mcDJ2Alc7/5b+BOh7JFQkcBWr5TfNCrj/i+0NoAzAdCC+rXHb5DGOMMQVU5t1KxhhjzsOKgzHGmAKsOBhjjCnAioMxxpgCrDgYY4wpwIqDMcaYAqw4GGOMKeD/A/0v3AEQ2T+8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta = 3**2\n",
    "y_hz = np.linspace(0, 80, 100)\n",
    "x_hz = 1\n",
    "plt.plot(y_hz, defaultdt * np.outer(y_hz * (y_hz-theta), x_hz) / (theta * tau_W))\n",
    "min_, max_ = plt.ylim()\n",
    "plt.vlines(theta, min_, max_, ls=':', color='k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate(\n",
    "    initialisation_steps=10,\n",
    "    lambda_=0.8,\n",
    "    total_time=runtime,\n",
    "    coherence=0.5,\n",
    "    plasticity=True,\n",
    "    W=None,\n",
    "    show_time=False\n",
    "):\n",
    "    ## Initialise arrays for computation\n",
    "    # For some reason this can't be numba.jit-ed\n",
    "    C_k = np.array([N_non] + [N_sub] * p + [N_I])\n",
    "    g_m = np.array([g_m_E] * (p+1) + [g_m_I])\n",
    "    C_m = np.array([C_m_E] * (p+1) + [C_m_I])\n",
    "    tau_m = np.array([tau_m_E] * (p+1) + [tau_m_I])\n",
    "    tau_rp = np.array([tau_rp_E] * (p+1) + [tau_rp_I])\n",
    "    nu = np.array([rate_pyramidal] * (p+1) + [rate_interneuron])\n",
    "    \n",
    "    # set weights\n",
    "    if W is None:\n",
    "        W = np.ones((p+2, p+2))  # from column to row\n",
    "        for i in range(0, p+1):\n",
    "            weights = np.full(p+2, w_minus)\n",
    "            if i > 0:\n",
    "                weights[i] = w_plus\n",
    "            else:\n",
    "                weights[i] = 1.0\n",
    "            weights[-1] = 1.0\n",
    "            weights[0] = 1.0\n",
    "            W[:,i] = weights\n",
    "    else:\n",
    "        assert W.shape[0] == p+2 \n",
    "        W = W.copy()\n",
    "\n",
    "    # AMPA\n",
    "    g_AMPA = np.array([g_AMPA_rec_E]* (p+1) + [g_AMPA_rec_I])\n",
    "    s_AMPA = tau_AMPA * nu\n",
    "    s_AMPA[~pyramidal_mask] = 0.  # inhibitory neurons won't feed AMPA-mediated synapses\n",
    "    ip_AMPA = (V_drive - V_E) * C_k * s_AMPA\n",
    "    ic_AMPA = g_AMPA * (W @ ip_AMPA)\n",
    "\n",
    "    # AMPA_ext\n",
    "    g_AMPA_ext = np.array([g_AMPA_ext_E]* (p+1) + [g_AMPA_ext_I])\n",
    "    s_AMPA_ext = np.full_like(\n",
    "        s_AMPA,\n",
    "        tau_AMPA * rate_ext\n",
    "    )  # array to allow for differing inputs\n",
    "    ip_AMPA_ext = (V_drive - V_E) * C_ext * s_AMPA_ext\n",
    "    ic_AMPA_ext = g_AMPA_ext * ip_AMPA_ext\n",
    "\n",
    "    # GABA\n",
    "    g_GABA = np.array([g_GABA_E]* (p+1) + [g_GABA_I])\n",
    "    s_GABA = tau_GABA * nu\n",
    "    s_GABA[pyramidal_mask] = 0.\n",
    "    ip_GABA = (V_drive - V_I) * C_k * s_GABA\n",
    "    ic_GABA = g_GABA * (W @  ip_GABA)\n",
    "\n",
    "    # NMDA - requires self-consistent calculation\n",
    "    g_NMDA = np.array([g_NMDA_E]* (p+1) + [g_NMDA_I])\n",
    "    # can I speed these up with @numba.jit?\n",
    "    g_NMDA_eff = lambda V: g_NMDA * J_2(V)\n",
    "    V_E_eff = lambda V: V - (1 / J_2(V)) * (V - V_E) / J(V)\n",
    "\n",
    "    s_NMDA = psi(nu)\n",
    "    s_NMDA[~pyramidal_mask] = 0.\n",
    "    V_avg_initial = -56e-3  # initial guess for determining V_avg\n",
    "    \n",
    "    # Initialise V_avg, V_SS\n",
    "    # can I speed this up with @numba.jit?\n",
    "    V_avg = np.full(p+2, V_avg_initial)\n",
    "    for k in range(initialisation_steps):\n",
    "        g_NMDA_eff_V = g_NMDA_eff(V_avg)\n",
    "        V_E_eff_V = V_E_eff(V_avg)\n",
    "        ip_NMDA = (V_drive - V_E_eff_V) * C_k * s_NMDA\n",
    "        ic_NMDA = g_NMDA_eff_V * (W @ ip_NMDA)\n",
    "        I_syn = ic_AMPA + ic_AMPA_ext + ic_NMDA + ic_GABA\n",
    "        V_SS = V_L - I_syn / g_m  # notice the minus, because current flows out?\n",
    "        V_avg = V_SS - (V_thr-V_reset)*nu*tau_m - (V_SS-V_reset)*nu*tau_rp\n",
    "    \n",
    "    # set sigma\n",
    "    sigma = np.sqrt(\n",
    "        g_AMPA_ext**2 * (V_drive - V_E)**2 * C_ext * rate_ext * tau_AMPA**2 / (g_m**2 * tau_m)\n",
    "    )\n",
    "    sigma[:-1] = lambda_ * (2e-3) + (1-lambda_) * sigma[:-1]\n",
    "\n",
    "    # initialise values for learning\n",
    "    theta_BCM = nu ** 2\n",
    "\n",
    "    ## Initialise arrays to track values, and for simulation\n",
    "    times = np.arange(0, total_time, defaultdt)\n",
    "#     nu_tracked = np.zeros((p+2, times.shape[0]))\n",
    "#     s_NMDA_tracked = np.zeros((p+2, times.shape[0]))\n",
    "#     s_AMPA_tracked = np.zeros((p+2, times.shape[0]))\n",
    "#     s_GABA_tracked = np.zeros((p+2, times.shape[0]))\n",
    "#     I_syn_tracked = np.zeros((p+2, times.shape[0]))\n",
    "    \n",
    "#     theta_BCM_tracked = np.zeros((p+2, times.shape[0]))\n",
    "#     W_tracked = np.zeros((p+2, p+2, times.shape[0]))\n",
    "    nu_tracked = np.full((p+2, times.shape[0]), np.nan)\n",
    "    s_NMDA_tracked = np.full((p+2, times.shape[0]), np.nan)\n",
    "    s_AMPA_tracked = np.full((p+2, times.shape[0]), np.nan)\n",
    "    s_GABA_tracked = np.full((p+2, times.shape[0]), np.nan)\n",
    "    I_syn_tracked = np.full((p+2, times.shape[0]), np.nan)\n",
    "    \n",
    "    theta_BCM_tracked = np.full((p+2, times.shape[0]), np.nan)\n",
    "    W_tracked = np.full((p+2, p+2, times.shape[0]), np.nan)\n",
    "    \n",
    "    if not plasticity:\n",
    "        W_tracked = np.full((p+2, p+2, times.shape[0]), W.reshape(p+2, p+2, 1))\n",
    "    \n",
    "    ic_noise = np.zeros_like(s_AMPA)\n",
    "    if show_time:\n",
    "        from tqdm import tqdm\n",
    "    else:\n",
    "        tqdm = lambda x: x\n",
    "    for itr, t in tqdm(enumerate(times)):\n",
    "        ip_AMPA = (V_drive - V_E) * C_k * s_AMPA\n",
    "        ic_AMPA = g_AMPA * (W @ ip_AMPA)\n",
    "\n",
    "        g_NMDA_eff_V = g_NMDA_eff(V_avg)\n",
    "        V_E_eff_V = V_E_eff(V_avg)\n",
    "        ip_NMDA = (V_drive - V_E_eff_V) * C_k * s_NMDA\n",
    "        ic_NMDA = g_NMDA_eff_V * (W @ ip_NMDA)\n",
    "\n",
    "        ip_GABA = (V_drive - V_I) * C_k * s_GABA\n",
    "        ic_GABA = g_GABA * (W @  ip_GABA)\n",
    "\n",
    "\n",
    "        ## LATEST CHANGE: added noise adapting to changing input by adapting s_AMPA_ext\n",
    "        s_AMPA_ext = np.full_like(\n",
    "            s_AMPA,\n",
    "            tau_AMPA * rate_ext\n",
    "        )\n",
    "        if t > 200e-3 and t < 400e-3:\n",
    "            # mu_0 ~= 0.05 in Wong&Wang2006, taken by comparison with mean external input\n",
    "            multiplier = np.ones(p+2)\n",
    "            multiplier[1] += 0.05 * (1+coherence)\n",
    "            multiplier[2] += 0.05 * (1-coherence)\n",
    "            s_AMPA_ext = s_AMPA_ext * multiplier\n",
    "\n",
    "            # increase input by roughly 1.5 times\n",
    "#             s_AMPA_ext[1] = s_AMPA_ext[1] * 1.5  #(25. * b2.Hz + rate_ext) / rate_ext\n",
    "        \n",
    "        ip_AMPA_ext = (V_drive - V_E) * C_ext * s_AMPA_ext\n",
    "        ic_AMPA_ext = g_AMPA_ext * ip_AMPA_ext\n",
    "        sigma = np.sqrt(\n",
    "            g_AMPA_ext**2 * (V_drive - V_E)**2 * C_ext * s_AMPA_ext * tau_AMPA / (g_m**2 * tau_m)\n",
    "        )\n",
    "        sigma[:-1] = lambda_ * (2e-3) + (1-lambda_) * sigma[:-1]\n",
    "        ## END OF LAST CHANGE\n",
    "\n",
    "        I_syn = ic_AMPA + ic_AMPA_ext + ic_NMDA + ic_GABA + ic_noise\n",
    "\n",
    "        dnu_dt_now = dnu_dt(\n",
    "            nu, I_syn,\n",
    "            sigma=sigma,\n",
    "            g_m=g_m,\n",
    "            tau_m=tau_m,\n",
    "            tau_rp=tau_rp\n",
    "        )\n",
    "        dic_noise_dt_now = dic_noise_dt(ic_noise)\n",
    "        ds_NMDA_dt_now = ds_NMDA_dt(s_NMDA, nu)\n",
    "        ds_AMPA_dt_now = ds_AMPA_dt(s_AMPA, nu)\n",
    "        ds_GABA_dt_now = ds_GABA_dt(s_GABA, nu)\n",
    "        \n",
    "        if plasticity:\n",
    "            dtheta_BCM_dt_now = dtheta_BCM_dt(theta_BCM, nu)\n",
    "            dW_dt_now = dW_dt_BCM(W, nu, theta_BCM)\n",
    "        \n",
    "\n",
    "        nu += dnu_dt_now * defaultdt\n",
    "        ic_noise += dic_noise_dt_now * defaultdt\n",
    "        s_NMDA += ds_NMDA_dt_now * defaultdt\n",
    "        s_AMPA += ds_AMPA_dt_now * defaultdt\n",
    "        s_GABA += ds_GABA_dt_now * defaultdt\n",
    "        \n",
    "        if plasticity:\n",
    "            theta_BCM += dtheta_BCM_dt_now * defaultdt\n",
    "            W += dW_dt_now * defaultdt\n",
    "            W = W.clip(0.0, np.inf)\n",
    "\n",
    "        nu_tracked[:, itr] = nu\n",
    "        s_NMDA_tracked[:, itr] = s_NMDA\n",
    "        s_AMPA_tracked[:, itr] = s_AMPA\n",
    "        s_GABA_tracked[:, itr] = s_GABA\n",
    "        I_syn_tracked[:, itr] = I_syn\n",
    "        \n",
    "        if plasticity:\n",
    "            theta_BCM_tracked[:, itr] = theta_BCM\n",
    "            W_tracked[:, :, itr] = W\n",
    "\n",
    "        # Not mentioned in W&W2006:\n",
    "        V_SS = V_L - I_syn / g_m\n",
    "        V_avg = V_SS - (V_thr-V_reset)*nu*tau_m - (V_SS-V_reset)*nu*tau_rp\n",
    "        \n",
    "        if np.any(np.isnan(nu)):\n",
    "            break\n",
    "\n",
    "    return  times,\\\n",
    "            nu_tracked,\\\n",
    "            s_NMDA_tracked,\\\n",
    "            s_AMPA_tracked,\\\n",
    "            s_GABA_tracked,\\\n",
    "            I_syn_tracked,\\\n",
    "            theta_BCM_tracked,\\\n",
    "            W_tracked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load_ext snakeviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %snakeviz simulate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "4516it [00:05, 1180.46it/s]<ipython-input-7-6bb2f7acddb1>:37: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
      "  the requested tolerance from being achieved.  The error may be \n",
      "  underestimated.\n",
      "  integral, error = integrate.quad(\n",
      "4635it [00:05, 1059.98it/s]<ipython-input-7-6bb2f7acddb1>:37: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
      "  If increasing the limit yields no improvement it is advised to analyze \n",
      "  the integrand in order to determine the difficulties.  If the position of a \n",
      "  local difficulty can be determined (singularity, discontinuity) one will \n",
      "  probably gain from splitting up the interval and calling the integrator \n",
      "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
      "  integral, error = integrate.quad(\n",
      "12000it [01:44, 114.96it/s]\n"
     ]
    }
   ],
   "source": [
    "# simulate(total_time=20e-3)\n",
    "times,\\\n",
    "nu_tracked,\\\n",
    "s_NMDA_tracked,\\\n",
    "s_AMPA_tracked,\\\n",
    "s_GABA_tracked,\\\n",
    "I_syn_tracked,\\\n",
    "theta_BCM_tracked,\\\n",
    "W_tracked = simulate(\n",
    "    total_time=0.6, coherence=0.7,\n",
    "    plasticity=False,\n",
    "    show_time=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "\n",
    "first_nan = np.any(np.isnan(W_tracked), axis=(0, 1)).argmax()\n",
    "plt.imshow(W_tracked[:,:,first_nan-1])\n",
    "plt.colorbar()\n",
    "for i in range(p+2):\n",
    "    for j in range(p+2):\n",
    "        plt.text(j, i, f\"{W_tracked[i,j,first_nan-1]:.2f}\",\n",
    "                       ha=\"center\", va=\"center\", color=\"w\")\n",
    "\n",
    "plt.figure()\n",
    "# plt.plot(times, np.max(W_tracked, axis=(0, 1)), label='max')\n",
    "# plt.plot(times, np.min(W_tracked, axis=(0, 1)), label='min')\n",
    "weight_lines = plt.plot(times, W_tracked[0, :, :].T)\n",
    "plt.legend(\n",
    "    weight_lines,\n",
    "    [f'{k}->non-sel' for k in ['non-sel']+list(range(1,p+1))+['inh'] ],\n",
    "#     title='neuron population',\n",
    ")\n",
    "\n",
    "sel_n = 1\n",
    "plt.figure()\n",
    "weight_lines = plt.plot(times, W_tracked[sel_n, :, :].T)\n",
    "plt.legend(\n",
    "    weight_lines,\n",
    "    [f'{k}->{sel_n}' for k in ['non-sel']+list(range(1,p+1))+['inh'] ],\n",
    "#     title='neuron population',\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 2, figsize=(14,11))\n",
    "\n",
    "# clear_output(wait=True)\n",
    "rates_lines = axes[0][0].plot(\n",
    "    times, nu_tracked.T,\n",
    "    alpha=0.8\n",
    ")\n",
    "axes[0][0].legend(\n",
    "    rates_lines,\n",
    "    ['non-sel'] + list(range(1, p+1)) +  ['inh'],\n",
    "    title='neuron population',\n",
    ")\n",
    "# _, ymax = axes[0][0].get_ylim()\n",
    "# axes[0][0].set_ylim([0., min(40., ymax)])\n",
    "\n",
    "axes[0][0].set_ylabel('Hz')\n",
    "axes[0][0].set_title('Firing Rates')\n",
    "axes[0][1].plot(times, s_NMDA_tracked[pyramidal_mask, :].T)\n",
    "axes[0][1].set_title('$s_{NMDA}$ (fraction of open NMDA-mediated channels)')\n",
    "axes[1][0].plot(times, s_AMPA_tracked[pyramidal_mask, :].T)\n",
    "axes[1][0].set_title('$s_{AMPA}$')\n",
    "axes[1][1].plot(times, s_GABA_tracked[~pyramidal_mask, :].T, color=rates_lines[-1].get_color())\n",
    "axes[1][1].set_title('$s_{GABA}$')\n",
    "axes[2][0].plot(times, I_syn_tracked.T * 1e12)\n",
    "axes[2][0].set_ylabel('pA')\n",
    "axes[2][0].set_title('$I_{syn}$ (outward current from neuron populations)')\n",
    "axes[2][1].plot(times, theta_BCM_tracked.T)\n",
    "axes[2][1].set_title(r'$\\theta_{BCM}$')\n",
    "axes[2][1].set_ylim([0., 40.])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show(fig)\n",
    "\n",
    "# add to discussion reasons (or lack thereof) for choosing times\n",
    "# for RL tasks\n",
    "\n",
    "# compare stable states in terms of s_NMDA levels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spiking Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parameters_spiking import *\n",
    "import brian2 as b2\n",
    "from brian2 import np, plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b2.start_scope()\n",
    "\n",
    "P_E = b2.NeuronGroup(\n",
    "    N_E, eqs_current_E, threshold='v > V_thr',\n",
    "    reset='v = V_reset', refractory='tau_rp_E', method='euler'\n",
    ")\n",
    "P_E.v = namespace['V_L']\n",
    "P_E.s_AMPA_ext = 4.8\n",
    "P_I = b2.NeuronGroup(\n",
    "    N_I, eqs_current_I, threshold='v > V_thr',\n",
    "    reset='v = V_reset', refractory='tau_rp_I', method='euler'\n",
    ")\n",
    "P_I.v = namespace['V_L']\n",
    "P_I.s_AMPA_ext = 4.8\n",
    "\n",
    "C_E_E = b2.Synapses(P_E, P_E, \n",
    "    model=eqs_glut, on_pre=eqs_pre_glut, on_post=eqs_post_glut,\n",
    "    method='euler', name='C_E_E')\n",
    "C_E_E.connect('i != j')\n",
    "C_E_E.w[:] = 1.0\n",
    "for pi in range(N_non, N_non + p * N_sub, N_sub):\n",
    "    # internal other subpopulation to current nonselective\n",
    "    C_E_E.w[C_E_E.indices[:, pi:pi + N_sub]] = w_minus\n",
    "    # internal current subpopulation to current subpopulation\n",
    "    C_E_E.w[C_E_E.indices[pi:pi + N_sub, pi:pi + N_sub]] = w_plus\n",
    "\n",
    "C_E_I = b2.Synapses(\n",
    "    P_E, P_I, model=eqs_glut, on_pre=eqs_pre_glut, method='euler',\n",
    "    name='C_E_I'\n",
    ")\n",
    "C_E_I.connect()\n",
    "C_E_I.w[:] = 1.0\n",
    "C_I_I = b2.Synapses(\n",
    "    P_I, P_I, on_pre=eqs_pre_gaba, method='euler',\n",
    "    name='C_I_I'\n",
    ")\n",
    "C_I_I.connect('i != j')\n",
    "# C_I_I.w[:] = 1.0\n",
    "C_I_E = b2.Synapses(P_I, P_E, on_pre=eqs_pre_gaba, method='euler')\n",
    "C_I_E.connect()\n",
    "# C_I_E.w[:] = 1.0\n",
    "\n",
    "C_P_E = b2.PoissonInput(\n",
    "    P_E, 's_AMPA_ext', C_ext, rate_ext*b2.Hz, '1',\n",
    "#         name='C_P_E'\n",
    ")\n",
    "C_P_I = b2.PoissonInput(\n",
    "    P_I, 's_AMPA_ext', C_ext, rate_ext*b2.Hz, '1',\n",
    "#         name='C_P_I'\n",
    ")\n",
    "\n",
    "\n",
    "\n",
    "C_selection = int(f * C_ext)\n",
    "stimtime = 1\n",
    "rate_selection = 25 * b2.Hz  # * (5/p)\n",
    "#     stimuli1 = b2.TimedArray(np.r_[np.zeros(40), np.ones(stimtime), np.zeros(100)], dt=25 * b2.ms)\n",
    "stimuli1 = b2.TimedArray(np.r_[np.zeros(8), np.ones(stimtime), np.zeros(100)], dt=25 * b2.ms)\n",
    "input1 = b2.PoissonInput(\n",
    "    P_E[N_non:N_non + N_sub], 's_AMPA_ext', C_selection, rate_selection, 'stimuli1(t)'\n",
    ")\n",
    "# stimuli2 = b2.TimedArray(np.r_[np.zeros(80), np.ones(stimtime), np.zeros(100)], dt=25 * b2.ms)\n",
    "# input2 = b2.PoissonInput(\n",
    "#     P_E[N_non + N_sub:N_non + 2 * N_sub], 's_AMPA_ext', C_selection, rate_selection, 'stimuli2(t)'\n",
    "# )\n",
    "# stimuli_reset = b2.TimedArray(np.r_[np.zeros(120), np.ones(2), np.zeros(100)], dt=25 * b2.ms)\n",
    "# input_reset_I = b2.PoissonInput(P_E, 's_AMPA_ext', C_ext, rate_selection, 'stimuli_reset(t)')\n",
    "# input_reset_E = b2.PoissonInput(P_I, 's_AMPA_ext', C_ext, rate_selection, 'stimuli_reset(t)')\n",
    "\n",
    "namespace['stimuli1'] = stimuli1\n",
    "# namespace['stimuli2'] = stimuli2\n",
    "# namespace['stimuli_reset'] = stimuli_reset\n",
    "\n",
    "N_activity_plot = 15\n",
    "sp_E_sels = [\n",
    "    b2.SpikeMonitor(P_E[pi:pi + N_activity_plot], name=f'sp_E_{int((pi-N_non)/N_sub) + 1}')\n",
    "    for pi in range(N_non, N_non + p * N_sub, N_sub)\n",
    "]\n",
    "sp_E = b2.SpikeMonitor(P_E[:N_activity_plot], name=f'sp_E_non')\n",
    "sp_I = b2.SpikeMonitor(P_I[:N_activity_plot], name=f'sp_I')\n",
    "\n",
    "r_E_sels = [\n",
    "    b2.PopulationRateMonitor(P_E[pi:pi + N_sub], name=f'r_E_{int((pi-N_non)/N_sub) + 1}')\n",
    "    for pi in range(N_non, N_non + p * N_sub, N_sub)]\n",
    "r_E = b2.PopulationRateMonitor(P_E[:N_non], name=f'r_E_non')\n",
    "r_I = b2.PopulationRateMonitor(P_I, name=f'r_I')\n",
    "\n",
    "st_E_sels = [\n",
    "    b2.StateMonitor(P_E[pi:pi + N_activity_plot], variables=True, record=True,\n",
    "                   name=f'st_E_{int((pi-N_non)/N_sub) + 1}')\n",
    "    for pi in range(N_non, N_non + p * N_sub, N_sub)]\n",
    "st_E = b2.StateMonitor(P_E[:N_activity_plot], variables=True, record=True, name=f'st_E_non')\n",
    "st_I = b2.StateMonitor(P_I[:N_activity_plot], variables=True, record=True, name=f'st_I')\n",
    "\n",
    "\n",
    "net = b2.Network(b2.collect())\n",
    "net.add(sp_E_sels)\n",
    "net.add(r_E_sels)\n",
    "net.add(st_E_sels)\n",
    "\n",
    "net.store('initialised')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "net.restore('initialised')\n",
    "\n",
    "net.run(\n",
    "    0.6*b2.second,\n",
    "    report='stdout',\n",
    "    namespace=namespace\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following plots can be used to compare the mean-field and true models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,6))\n",
    "\n",
    "st = net['st_E_non']\n",
    "mean_input = np.mean(st.I_syn / b2.pA, axis=0)\n",
    "std_input = np.std(st.I_syn / b2.pA, axis=0)\n",
    "b2.plot(st.t / b2.ms, mean_input, label='non-sel.')\n",
    "# b2.fill_between(st.t / b2.ms, mean_input-std_input, mean_input+std_input, alpha=0.5)\n",
    "\n",
    "for group in range(1,p+1):\n",
    "    st = net[f'st_E_{group}']\n",
    "    mean_input = np.mean(st.I_syn / b2.pA, axis=0)\n",
    "    std_input = np.std(st.I_syn / b2.pA, axis=0)\n",
    "    b2.plot(st.t / b2.ms, mean_input, label=f'sel. {group}')\n",
    "#     b2.fill_between(st.t / b2.ms, mean_input-std_input, mean_input+std_input, alpha=0.5)\n",
    "\n",
    "st = net[f'st_I']\n",
    "mean_input = np.mean(st.I_syn / b2.pA, axis=0)\n",
    "std_input = np.std(st.I_syn / b2.pA, axis=0)\n",
    "b2.plot(st.t / b2.ms, mean_input, label=f'inhib.')\n",
    "# b2.fill_between(st.t / b2.ms, mean_input-std_input, mean_input+std_input, alpha=0.5)\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "from matplotlib.patches import Rectangle\n",
    "\n",
    "\n",
    "def plot_average_membrane_potential(\n",
    "    axes, population='E_non', group='non-selective',\n",
    "    tmin=None, tmax=None, legend=True, title=True\n",
    "):\n",
    "    # using formula:\n",
    "    V_SS = V_L * b2.volt - net['st_'+population].I_syn.mean(axis=0) / (g_m_E * b2.siemens)\n",
    "    nu = net['r_'+population].smooth_rate('gaussian', 5*b2.ms)\n",
    "    V_avg = V_SS - (V_thr - V_reset)*b2.volt * nu * tau_m_E * b2.second\n",
    "    # using data:\n",
    "    v_mean = b2.mean(net['st_'+population].v, axis=0)\n",
    "    v_std = b2.std(net['st_'+population].v, axis=0)\n",
    "    \n",
    "    if tmin is None:\n",
    "        tmin = net['r_'+population].t[0]\n",
    "    if tmax is None:\n",
    "        tmax = net['r_'+population].t[-1]\n",
    "    mask = (net['r_I'].t >= tmin) & (net['r_I'].t <= tmax)\n",
    "    \n",
    "    true_line = axes.plot(\n",
    "        net['r_'+population].t[mask] / b2.ms,\n",
    "        v_mean[mask] / b2.mV,\n",
    "        label='true avg.',\n",
    "        ls='--'\n",
    "    )\n",
    "    colour = true_line[0].get_color()\n",
    "    axes.fill_between(\n",
    "        net['r_'+population].t[mask] / b2.ms,\n",
    "        (v_mean-v_std)[mask]/b2.mV,\n",
    "        (v_mean+v_std)[mask]/b2.mV,\n",
    "        alpha=0.6,\n",
    "        label=r'$\\pm 1$ std',\n",
    "        color=colour\n",
    "    )\n",
    "    axes.fill_between(\n",
    "        net['r_'+population].t[mask] / b2.ms,\n",
    "        (v_mean-2*v_std)[mask]/b2.mV,\n",
    "        (v_mean-v_std)[mask]/b2.mV,\n",
    "        alpha=0.3,\n",
    "        label=r'$\\pm 2$ std',\n",
    "        color=colour\n",
    "    )\n",
    "    axes.fill_between(\n",
    "        net['r_'+population].t[mask] / b2.ms,\n",
    "        (v_mean+2*v_std)[mask]/b2.mV,\n",
    "        (v_mean+v_std)[mask]/b2.mV,\n",
    "        alpha=0.3,\n",
    "        color=colour\n",
    "    )\n",
    "\n",
    "    axes.plot(\n",
    "        net['r_'+population].t[mask] / b2.ms,\n",
    "        V_avg[mask] / b2.mV,\n",
    "        label='computed',\n",
    "        color='darkorange',\n",
    "        lw=1\n",
    "    )\n",
    "    if legend:\n",
    "        axes.legend()\n",
    "    axes.grid(alpha=0.5)\n",
    "    if title:\n",
    "        axes.set_title(\n",
    "            r'Comparison of $\\langle V \\rangle$ formula to true values'\n",
    "            + f'\\nGroup: {group}'\n",
    "        )\n",
    "    return None\n",
    "\n",
    "        \n",
    "fig, axes = plt.subplots(1,2,figsize=(14, 5), sharey=True)\n",
    "\n",
    "axes[0].set_xlabel('time (ms)')\n",
    "axes[1].set_xlabel('time (ms)')\n",
    "axes[0].set_ylabel('membrane potential (mV)')\n",
    "\n",
    "plot_average_membrane_potential(\n",
    "    axes[0],\n",
    "    tmax=400*b2.ms\n",
    ")\n",
    "plot_average_membrane_potential(\n",
    "    axes[1], population='E_1', group='selective-1',\n",
    "    tmax=400*b2.ms,\n",
    "    legend=False\n",
    ")\n",
    "\n",
    "inset_ax = inset_axes(\n",
    "    axes[1],\n",
    "    width=\"40%\", # width = 30% of parent_bbox\n",
    "    height=1.2, # height : 1 inch\n",
    "    loc=4,\n",
    "    borderpad=2.0,\n",
    "    axes_kwargs={'alpha':0.5}\n",
    ")\n",
    "plot_average_membrane_potential(\n",
    "    inset_ax, population='E_1', group='selective-1',\n",
    "    legend=False, title=False,\n",
    "    tmin=190*b2.ms, tmax=240*b2.ms\n",
    ")\n",
    "inset_ax.set_xticklabels([])\n",
    "inset_ax.set_yticklabels([])\n",
    "\n",
    "left, right = inset_ax.get_xlim()\n",
    "bottom, top = inset_ax.get_ylim()\n",
    "\n",
    "ax_ticksy = axes[1].get_yticks()\n",
    "ax_ticksy = ax_ticksy[ax_ticksy <= top]\n",
    "ax_ticksy = ax_ticksy[ax_ticksy >= bottom]\n",
    "inset_ax.set_yticks(ax_ticksy)\n",
    "\n",
    "ax_ticksx = axes[1].get_xticks()\n",
    "ax_ticksx = ax_ticksx[ax_ticksx <= right]\n",
    "ax_ticksx = ax_ticksx[ax_ticksx >= left]\n",
    "inset_ax.set_xticks(ax_ticksx)\n",
    "\n",
    "rectangle_patch = Rectangle(\n",
    "    (left, bottom), right-left, top-bottom,\n",
    "    fill=False, clip_on=False, alpha=0.9, ls='--'\n",
    "    )\n",
    "\n",
    "axes[1].add_patch(rectangle_patch)\n",
    "\n",
    "fig.subplots_adjust(wspace=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b2.figure(figsize=(12, 6))\n",
    "\n",
    "window = 2*b2.ms\n",
    "\n",
    "nu_E = net['r_E_non']\n",
    "b2.plot(\n",
    "    nu_E.t / b2.ms,\n",
    "    nu_E.smooth_rate('gaussian', window),\n",
    "#     color='r',\n",
    "    label='non-sel.')\n",
    "\n",
    "for group in range(1,p+1):\n",
    "    nu_E_sel = net[f'r_E_{group}']\n",
    "    b2.plot(\n",
    "        nu_E_sel.t / b2.ms,\n",
    "        nu_E_sel.smooth_rate('gaussian', window),\n",
    "#         color='g',\n",
    "        label=f'sel. {group}')\n",
    "\n",
    "nu_I = net['r_I']\n",
    "b2.plot(\n",
    "    nu_I.t / b2.ms,\n",
    "    nu_I.smooth_rate('gaussian', window),\n",
    "#     color='b',\n",
    "    label='inhib.')\n",
    "b2.legend()\n",
    "b2.ylabel('Hz')\n",
    "b2.title('Firing Rates')\n",
    "b2.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# b2.title('Population rates')\n",
    "# b2.xlabel('ms')\n",
    "# b2.ylabel('Hz')\n",
    "\n",
    "# b2.plot(\n",
    "#     r_E.t / b2.ms, r_E.smooth_rate(width=5 * b2.ms) / b2.Hz,\n",
    "#     label='nonselective')\n",
    "# b2.plot(\n",
    "#     r_I.t / b2.ms, r_I.smooth_rate(width=5 * b2.ms) / b2.Hz,\n",
    "#     label='inhibitory')\n",
    "\n",
    "# for i, r_E_sel in enumerate(r_E_sels[::-1]):\n",
    "#     b2.plot(\n",
    "#         r_E_sel.t / b2.ms, r_E_sel.smooth_rate(width=5 * b2.ms) / b2.Hz,\n",
    "#         label='selective {}'.format(p - i))\n",
    "\n",
    "# b2.legend()\n",
    "# b2.figure()\n",
    "# b2.title('Population activities ({} neurons/pop)'.format(N_activity_plot))\n",
    "# b2.xlabel('ms')\n",
    "# b2.yticks([])\n",
    "\n",
    "# b2.plot(\n",
    "#     sp_E.t / b2.ms, sp_E.i + (p + 1) * N_activity_plot, '.',\n",
    "#     markersize=2, label='nonselective')\n",
    "# b2.plot(\n",
    "#     sp_I.t / b2.ms, sp_I.i + p * N_activity_plot, '.',\n",
    "#     markersize=2, label='inhibitory')\n",
    "\n",
    "# for i, sp_E_sel in enumerate(sp_E_sels[::-1]):\n",
    "#     b2.plot(\n",
    "#         sp_E_sel.t / b2.ms, sp_E_sel.i + (p - i - 1) * N_activity_plot, '.',\n",
    "#         markersize=2, label='selective {}'.format(p - i))\n",
    "\n",
    "# b2.legend()\n",
    "# b2.show()"
   ]
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}