Networks in NeuroDyn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Networks in NeuroDyn\n",
    "\n",
    "## Overview\n",
    "\n",
    "This notebook provides a guide to specyfing neural interconnections and designing networks using Hodgkin-Huxley and NeuroDyn models in Python environment. We first start with an overview of two main types of neural interconnections: gap junctions and synapses, and how to implement them in the software environment. The environment allows defining arbitrary conductance-based models of synapses, and we give an example here using an excitatory AMPA synapse. Networks can be constructed using both the standard Hodgkin-Huxley models as well as NeuroDyn neurons, and the notebook shows how we can specify both types. \n",
    "\n",
    "Finally, we introduce two additional classes `NeuroDynBoard` and `NeuroCube` which specify the NeuroDyn hardware structure:\n",
    "- `NeuroDynBoard` is a single NeuroDyn board containing 4 neurons and 12 synapses.\n",
    "- `NeuroCube` is a parallel interconnection of 4 NeuroDyn boards, where the boards are connected through mutural resistive connections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the environment\n",
    "\n",
    "Let's first import the required Python modules and model definitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from cb_models import HHModel, NeuroDynModel, NeuronalNetwork, Synapse, NDSynapse\n",
    "from cb_models import AMPA\n",
    "from cb_models import NeuroDynBoard, NeuroCube\n",
    "from fitting_utilities import FitND\n",
    "\n",
    "# **Ignore overflow warnings**\n",
    "import numpy as np\n",
    "np.seterr(all=\"ignore\")\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview of neural interconnections\n",
    "\n",
    "There are two main types of neural interconnections that we can implement: **gap junctions** and **synaptic connections**. In the conductance-based modelling framework, these interconnections lead to additional currents in the membrane equation of neuron $i$ that are dependent on membrane voltages of other neurons $j$ connected to it:\n",
    "\n",
    "\\begin{equation}\n",
    "C \\dot{V}_i = - I_{int,i} - \\sum_{j} I_{gap, i, j} - \\sum_{j} I_{syn, i, j} + I_{app, i}\n",
    "\\end{equation}\n",
    "\n",
    "where $I_{int,i}$ is total internal current of the neuron, $I_{gap, i, j}$ is the current due to the gap junction between neurons $j$ and $i$, $I_{syn, i, j}$ is the synaptic current from neuron $j$, and $I_{app, i}$ is the external current applied to neuron $i$.\n",
    "\n",
    "### Gap junctions\n",
    "\n",
    "Gap junctions are modeled as simple resistive connections between neurons, so that a resistive connection between neurons $i$ and $j$ leads to a current in neuron $i$:\n",
    "\n",
    "\\begin{equation}\n",
    "I_{gap, i, j} = g_{gap, i, j} (V_i - V_j)\n",
    "\\end{equation}\n",
    "\n",
    "Gap junction interconnections are symmetric $g_{gap, i, j} = g_{gap, j, i}$, so that:\n",
    "\n",
    "\\begin{equation}\n",
    "I_{gap, j, i} = -I_{gap, i, j}\n",
    "\\end{equation}\n",
    "\n",
    "and neuron $j$ receives the same current of the opposite sign.\n",
    "\n",
    "We can represent this as a matrix $G_{gap}$ where each element $g_{gap, i, j}$ is the conductance of the connection between neurons $i$ and $j$. This matrix is necessarily **symmetric**, i.e. $G_{gap} = G^T_{gap}$.\n",
    "\n",
    "### Synaptic connections\n",
    "\n",
    "Unlike gap junctions, synaptic connections are directed connections and are not symmetric. If there is a synapse from neuron $j$ to neuron $i$, then the activity of the *presynaptic* neuron ($j$) modulates the activity of the *postsynaptic* neuron ($i$). Synaptic currents have a similar form to the internal ionic currents, with the difference that their activation depends on the activity of the presynaptic neuron. This is observed in the dependency of the opening rate $\\alpha_r$ on the presynaptic neuron voltage $V_j$:\n",
    "\n",
    "\\begin{equation}\n",
    "I_{syn, i, j} = \\bar{g}_{syn, i, j} r (V_i - E_{syn}) \\\\\n",
    "\\dot{r} = \\alpha_r(V_j) (1 - r) - \\beta_r(V_i) r\n",
    "\\end{equation}\n",
    "\n",
    "The opening of the channels $\\alpha_r$ depends on the presynaptic neuron's voltage $V_j$, while the closing depends on the postsynaptic neuron's voltage $V_i$ (sometimes this relationship is modeled as constant, so that $\\beta_r = const$).\n",
    "\n",
    "The reversal potential of the synapse $E_{syn}$ determines if the synapse is *excitatory* or *inhibitory*: excitatory synapses tend to increase the postsynaptic neuron's voltage, while the inhibitory synapses act to decrease the postynaptic voltage. On the other hand, dynamics of activation $r$ determines the timescale of the synaptic current, i.e. how fast it activates in relation to the changes in the presynaptic neuron."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining neural interconnections\n",
    "\n",
    "Let's now see how we can define these interconnections in the Python environment. We start with a simple example of two Hodgkin-Huxley neurons, which we can define as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neuron1 = HHModel()\n",
    "neuron2 = HHModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gap junctions\n",
    "\n",
    "Let's start with a simple network consisting of a single gap junction between these two neurons with $g_{gap} = 0.1$. We do this by defining the matrix `G_gap` and passing it as a keyword argument along with a list of the neurons to our `NeuronalNetwork` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neurons = [neuron1, neuron2]\n",
    "g = 0.1\n",
    "G_gap = [[0, g], [g, 0]] # needs to be symmetric\n",
    "\n",
    "network = NeuronalNetwork(neurons, gap = G_gap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate the network exactly as before by calling the method `simulate`. This time, the applied current function needs to return an array of external currents corresponding to each neuron in the `neurons` list. Additionally, we need to make sure that the initial condition array `x0` has the right number of states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulation time\n",
    "T = 200\n",
    "trange = (0, T)\n",
    "\n",
    "# External current\n",
    "I1 = 15 # current to neuron 1\n",
    "I2 = 0 # current to neuron 2\n",
    "Iapp = lambda t: [I1, I2]\n",
    "\n",
    "# Initial condition\n",
    "x0 = [0, 0, 0, 0] * 2\n",
    "\n",
    "# Simulate\n",
    "sol = network.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time [ms]\")\n",
    "plt.ylabel(\"Voltage [mV]\")\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can change the conductance of the gap junction and observe the different behavior. For very small values of the conductance, the neurons are effectively disconnected and act independently, while gradually increasing the conductance will bring the membrane voltages of the two neurons closer together.\n",
    "\n",
    "Making $g_{gap}$ large effectively acts as the short-circuit between the two neurons (but be careful about making it too large as it may lead to numerical issues)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Synaptic connections\n",
    "\n",
    "Let's now see how we can introduce a synapse from neuron 1 to neuron 2. We first need to create an object for the synapse of the `Synapse` class which takes three arguments: `gsyn` - maximal conductance, `Esyn` - reversal potential, `r` - gating variable object.\n",
    "\n",
    "The `r` gating variable should be derived from the class `HHKinetics` and needs to provide methods `alpha` and `beta` that define the functions $\\alpha_r(V)$ and $\\beta_r(V)$.\n",
    "\n",
    "#### AMPA synapse\n",
    "\n",
    "As an example, let's take a look at implementing an excitatory AMPA synapse. There is an `AMPA` gating kinetics class already provided that takes the values from Ermentrout et al. 2010 that we can use here. Let's define a synapse using the `AMPA` kinetics class, setting the value of the reversal potential to $E_{AMPA} = 65mV$ (the reversal potential value is shifted to accomodate for the shift in the Hodgkin-Huxley model that sets $V_{rest} = 0mV$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = AMPA()\n",
    "gsyn = 0.05\n",
    "Esyn = 65\n",
    "synapse = Synapse(gsyn, Esyn, r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having defined the model of the synapse that we would like to use, we now need to define the connectivity matrix describing all of the synaptic connections in the network. For this, we will construct a matrix `syns` where each element of the matrix $(i, j)$ contains a **list** of all synaptic connections from neuron $j$ to neuron $i$. If there are no synaptic connections from neuron $j$ to neuron $i$, then we set that matrix element to `None`.\n",
    "\n",
    "Since we would like to define a single excitatory synapse from neuron 1 to neuron 2, we will put a single `Synapse` object into a list at position (1, 0), filling the rest of the matrix with `None`. Let's create a new network with the single excitatory synapse and simulate it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syns = [[None, None], [[synapse], None]]\n",
    "\n",
    "network = NeuronalNetwork(neurons, syns = syns)\n",
    "\n",
    "# Simulation time\n",
    "T = 200\n",
    "trange = (0, T)\n",
    "\n",
    "# External current\n",
    "I1 = 15 # current to neuron 1\n",
    "I2 = 0 # current to neuron 2\n",
    "Iapp = lambda t: [I1, I2]\n",
    "\n",
    "# Initial condition\n",
    "x0 = [0, 0, 0, 0] * 2 + [0] # extra state due to the synapse!\n",
    "\n",
    "# Simulate\n",
    "sol = network.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time [ms]\")\n",
    "plt.ylabel(\"Voltage [mV]\")\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can observe that each time the first neuron spikes it excites the second neuron. If the strength of the synapse is large enough, the second neuron will start spiking in response to the activity of the second neuron. We can play around with the synaptic strength and observe the different behavior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "synapse.gsyn = 0 # change the synapse strength\n",
    "\n",
    "# Simulate\n",
    "sol = network.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time [ms]\")\n",
    "plt.ylabel(\"Voltage [mV]\")\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Defining different synapse models\n",
    "\n",
    "You can use the `Synapse` class to define arbitrary synaptic connections - the main requirement is implementing a gating variable class inherited from the `HHKinetics` class that needs to define implementations for the methods `alpha(V)` and `beta(V)`. You can take a look at how the `AMPA` class is defined in the `cb_models` file to get a better idea.\n",
    "\n",
    "**Note** that since our original Hodgkin-Huxley model has parameter values shifted so that the resting potential is at $V_{rest} = 0mV$, you might need to shift biophyiscal values accordingly in order to make the synapse compatible with the model.\n",
    "\n",
    "Another important thing to note is that when we try to fit activation dynamics to the NeuroDyn hardware, the values should be in **SI units**.\n",
    "\n",
    "Finally, if we would like to stretch the voltage range of the Hodgkin-Huxley model using the `scl_v` parameter, the range of the synapse would need to be stretched accordingly to accommodate for this. It is therefore best to implement an optional keyword `scl_v` to the gating variable class that would multiply the voltage-related parameters (you can check this for the `AMPA` class and the `HHModel` class in the file `cb_models`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NeuroDyn network\n",
    "\n",
    "Let's see now how we can do the same using NeuroDyn models for the individual neurons, as well as the synapses.\n",
    "\n",
    "### NeuroDyn synapse\n",
    "\n",
    "A synaptic connection in the NeuroDyn hardware is implemented in very much the same way as ionic conductance dynamics was for a single neuron (see \"NeuroDyn equations\" section in **NeuroDyn Python model** notebook). Therefore, a NeuroDyn synapse class `NDSynapse` will be defined using the following digital parameters:\n",
    "\n",
    "- `dg`: digital value for the maximal conductance\n",
    "- `dE`: digital value for the reversal potential\n",
    "- `dIb`: digital values for the alpha and beta sigmoid weights\n",
    "    - `dIb` is a list containing `[dIalpha, dIbeta]`\n",
    "    - each `dIalpha` and `dIbeta` is an array of 7 digital values for the sigmoid weights\n",
    "- `ND`: optional keyword argument containing a `NeuroDynModel` object\n",
    "    - the synapse will take all physical parameters such as `I_master` and `I_voltage` from this object, so it should be the 'parent' NeuroDyn neuron for the synapse\n",
    "    - if no `ND` object is passed, `NDSynapse` will be initialized with a default `NeuroDynModel` object\n",
    "    \n",
    "In order to find these digital values we can again try fitting the parameters so that they replicate a particular biophysical synapse. As an example, let's try fitting the parameters to the AMPA synapse of the previous section.\n",
    "\n",
    "### Fitting a synapse\n",
    "\n",
    "In order to fit a synapse model, we can again use the class `FitND` that we have previously used to fit a Hodgkin-Huxley model to a NeuroDyn neuron. We first fit an appropriate Hodgkin-Huxley model as before, after which we can pass any other gating variables that we would like to fit to. In this case, we want to fit the AMPA gating variable.\n",
    "\n",
    "First, let's repeat the procedure of fitting a NeuroDyn neuron as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit Hodgkin-Huxley model\n",
    "\n",
    "scl_v = 3 # voltage scaling\n",
    "\n",
    "# Hodgkin-Huxley model we want to fit to\n",
    "hh_neuron = HHModel(scl_v = scl_v, SI_units = True)\n",
    "\n",
    "# Create a fitting object and obtain the initial fit\n",
    "fit_obj = FitND(hh_neuron)\n",
    "\n",
    "# Get the sigmoid weights\n",
    "weights = fit_obj.fitHH()\n",
    "\n",
    "g0 = [120e-3,36e-3,0.3e-3] # maximal conductances\n",
    "E0 = [120e-3,-12e-3,10.6e-3] # reversal potentials\n",
    "\n",
    "# Get the digital parameters for NeuroDyn\n",
    "dIb = fit_obj.get_digital_Ib(weights)\n",
    "dg = fit_obj.get_digital_g(g0)\n",
    "dE = fit_obj.get_digital_E(E0)\n",
    "dIb[2][1] = dIb[2][1]*15\n",
    "\n",
    "I_master = fit_obj.I_master\n",
    "I_voltage = fit_obj.I_voltage\n",
    "\n",
    "V_ref = 0.9\n",
    "\n",
    "# Create a NeuroDyn object\n",
    "ND1 = NeuroDynModel(dg, dE, dIb, V_ref, I_voltage, I_master)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's follow the same procedure in order to get the parameters for an AMPA synapse. We first make an AMPA gating variable (remember that it needs to be scaled appropriately, which we do by passing the `scl_v` keyword).\n",
    "\n",
    "We can then pass the `AMPA` gating variable, along with any other gating variables we would like to fit to, to the `fit` method. Importantly, note that a **list** of gating variables should be passed as an argument, along with an optional list `labels` for the figure labels.\n",
    "\n",
    "Another important thing to note is that each time we fit new gating variables using an existing `FitND` object, the `scl_t` variable determining the time scaling will adjust so that the maximal coefficient encountered so far fits within the NeuroDyn constraints. If we would like to fit gating variables without readjusting the scale, we can pass an optional keyword argument `update_scale = False`. This means that our fitting process might return parameters that are outside the digital range. Here we will keep the default `update_scale` that is set to `True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit AMPA gating variable\n",
    "r = AMPA(scl_v = scl_v, SI_units = True)\n",
    "\n",
    "# Fit weights to AMPA\n",
    "weights_syn = fit_obj.fit([r], labels = ['r'], plot_alpha_beta = True) # note that gating variables are passed in a list\n",
    "\n",
    "weights_syn = weights_syn[0] # fitting method returns a list of gating variables, so get the only element"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having obtained the weights for the gating variable kinetics, we can then get the digital parameters for the maximal conductance and the reversal potential of the synapse, and use that to create an `NDSynapse` object. Apart from the maximal conductance, reversal potential, and the gating variable, the `NDSynapse` object also accepts a reference to a NeuroDyn object from which it takes physical constants and currents $I_{voltage}$ and $I_{master}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsyn = 0.3e-3 # max coductance of the synapse\n",
    "Esyn = 65e-3 # reversal potential of the synapse\n",
    "\n",
    "dIb_syn = fit_obj.get_digital_Ib(weights_syn)\n",
    "dgsyn = fit_obj.get_digital_g(gsyn)\n",
    "dEsyn = fit_obj.get_digital_E(Esyn)\n",
    "\n",
    "nd_syn = NDSynapse(dgsyn, dEsyn, dIb_syn, ND1)\n",
    "\n",
    "print(\"Alpha and beta weights:\\n\", dIb_syn)\n",
    "print(\"Digital value for maximal conductance: \", dgsyn)\n",
    "print(\"Digital value for the reversal potential: \", dEsyn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a network\n",
    "\n",
    "Having created an object for our synapse, we can now define a network in exactly the same way as before! Let's define a second neuron `ND2` with the same parameters, but let's decrease the applied current to this neuron in order to make it silent. For convenience, we can `deepcopy` our first neuron, so that the second neuron has the same parameters but does not share any variables with the first neuron. Let's set the appropriate current $I_{app}$ so that the neuron is silent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "\n",
    "ND2 = deepcopy(ND1) # makes an independent copy of the ND1 object\n",
    "\n",
    "# Set the applied current for ND2 to be silent\n",
    "I2 = -1e-6\n",
    "Iapp = lambda t : fit_obj.convert_I(I2)\n",
    "\n",
    "V_ref = 0.9\n",
    "\n",
    "T = 0.02\n",
    "trange = (0, T)\n",
    "\n",
    "sol = ND2.simulate(trange,[0.7,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t [s]')\n",
    "plt.ylabel('V [V]')\n",
    "plt.title('ND 2')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our two neurons `ND1` and `ND2`, let's set up the excitatory synapse we have created from `ND1` to `ND2` so that the first neuron spiking excites the resting second neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nd_neurons = [ND1, ND2]\n",
    "\n",
    "nd_syns = [[None, None], [[nd_syn], None]]\n",
    "\n",
    "nd_network = NeuronalNetwork(nd_neurons, syns = nd_syns)\n",
    "\n",
    "# Simulation time\n",
    "T = 0.04\n",
    "trange = (0, T)\n",
    "\n",
    "# External currents\n",
    "I1 = fit_obj.convert_I(0e-6) # current to neuron 1\n",
    "I2 = fit_obj.convert_I(-1e-6) # current to neuron 2\n",
    "Iapp = lambda t: [I1, I2]\n",
    "\n",
    "# Initial condition\n",
    "x0 = [0.7, 0, 0, 0] * 2 + [0] # extra state due to the synapse\n",
    "\n",
    "# Simulate\n",
    "sol = nd_network.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t [s]')\n",
    "plt.ylabel('V [V]')\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may notice that after the initial transient has died out, individual spikes in the first neuron induce excitations in the voltage of the second neuron. We can try changing the digital value of the maximal conductance of the synapse to observe its effect on the activity of the second neuron. You should observe that by increasing the synaptic strength the second neuron eventually starts spiking, and its spikes appear more closely to the spikes of the first neuron as it is further increased."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change synapse strength\n",
    "nd_syn.update_dg(0)\n",
    "\n",
    "# Simulate\n",
    "sol = nd_network.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t [s]')\n",
    "plt.ylabel('V [V]')\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NeuroDyn board classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we have seen how to define a network of Hodgkin-Huxley or NeuroDyn neurons. The arguments of the network are a full list of the individual neural objects, a conductance matrix defining the gap junction connections, and a synapse matrix defining the synapse connections of the circuit. Using the `NeuronalNetwork` class we can define arbitrary neural networks given by the structure of the gap junction and synapse matrices.\n",
    "\n",
    "The NeuroDyn hardware however has a fixed structure: a single **NeuroDyn board** consists of **4 neurons** with **12 synapses** defining a network of four fully-connected neurons. The class `NeuroDynBoard` defines a network with this particular architecture, allowing in addition to optionally connect neurons with resistive connections.\n",
    "\n",
    "By connecting four NeuroDyn boards in parallel through resistive connections, we get a **NeuroCube**. The NeuroCube consists of four NeuroDyn boards stacked on top of each other, so that Neuron 1 of each board can be connected to corresponding Neuron 1 of every other board, Neuron 2 of every board can be connected to every correspodning Neuron 2, and so on, giving in total 24 additional resistive connections that can be controlled. This network is represented by the `NeuroCube` class.\n",
    "\n",
    "Let's see how we can use these classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NeuroDynBoard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a single NeuroDyn board object by calling the class `NeuroDynBoard`, which will initialize all neurons and synapses with default parameters, as well as set conductance of all resistive connections to $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nd_board = NeuroDynBoard()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try simulating the board - since all synapses and gap junctions are initialized with $0$ conductance, we should observe four neurons operating individually. Note that we can reuse are previously defined fitting object in order to convert currents to the right scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulation time\n",
    "T = 0.02\n",
    "trange = (0, T)\n",
    "\n",
    "# External currents\n",
    "I1 = fit_obj.convert_I(-2e-6) # silent\n",
    "I2 = fit_obj.convert_I(-1e-6) # silent\n",
    "I3 = fit_obj.convert_I(0e-6) # spiking\n",
    "I4 = fit_obj.convert_I(2e-6) # spiking\n",
    "Iapp = lambda t: [I1, I2, I3, I4]\n",
    "\n",
    "# Initial conditions\n",
    "x0 = [0.7, 0, 0, 0] * 4 + [0]*12 # 4 neurons, 12 synapses\n",
    "\n",
    "# Simulate\n",
    "sol = nd_board.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "V3 = sol.y[8]\n",
    "V4 = sol.y[12]\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"NeuroDyn Board\")\n",
    "plt.xlabel('t [s]')\n",
    "plt.ylabel('V [V]')\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.plot(sol.t, V3, label = \"Neuron 3\")\n",
    "plt.plot(sol.t, V4, label = \"Neuron 4\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can change the parameters of the circuit by accessing individual neurons and synapses. We can do this by using the following methods:\n",
    "\n",
    "- `get_neuron(i)`: retrive neuron $i$ ($i$ goes from $0$ to $3$)\n",
    "- `get_syn(i, j)`: retrieve synapse from neuron $j$ to neuron $i$\n",
    "\n",
    "We can then use the neural and synaptic class methods to update the parameters. In addition we can set the gap junction conductances by calling:\n",
    "\n",
    "- `set_gap(g, i, j)`: set the maximal conductance of the $(i, j)$ gap junction to $g$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a test, let's see how we can create a network where neurons excite each other sequentially 1 -> 2 -> 3 -> 4. The default parameters are set to an excitatory synapse, so we only need to increase the maximal conductances of the appropriate connections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get synapses\n",
    "syn1 = nd_board.get_syn(1, 0) # synapse 1 -> 2\n",
    "syn2 = nd_board.get_syn(2, 1) # synapse 2 -> 3\n",
    "syn3 = nd_board.get_syn(3, 2) # synapse 3 -> 4\n",
    "\n",
    "# Update synapses\n",
    "syn1.update_dg(5)\n",
    "syn2.update_dg(5)\n",
    "syn3.update_dg(2)\n",
    "\n",
    "# External currents (set neurons 2, 3, 4 to silent)\n",
    "I1 = -0e-6\n",
    "I2 = -1e-6\n",
    "I3 = -1e-6\n",
    "I4 = -1e-6\n",
    "\n",
    "Iapp = lambda t: fit_obj.convert_I([I1, I2, I3, I4])\n",
    "\n",
    "# Initial conditions\n",
    "x0 = [0.7, 0, 0, 0] * 4 + [0]*12 # 4 neurons, 12 synapses\n",
    "\n",
    "# Simulate\n",
    "sol = nd_board.simulate(trange, x0, Iapp)\n",
    "V1 = sol.y[0]\n",
    "V2 = sol.y[4]\n",
    "V3 = sol.y[8]\n",
    "V4 = sol.y[12]\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"NeuroDyn Board with excitatory synapses\")\n",
    "plt.xlabel('t [s]')\n",
    "plt.ylabel('V [V]')\n",
    "plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "plt.plot(sol.t, V3, label = \"Neuron 3\")\n",
    "plt.plot(sol.t, V4, label = \"Neuron 4\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can change the maximal conductance parameters of the synapses to make sure the network behaves as expected (for example turning the synapses on or off)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NeuroCube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In exactly the same manner, we can define an object of the `NeuroCube` class which will initialize 4 disconnected NeuroDyn boards with the default parameters. Let's first set up a default neurocube object and try simulating it. We can for example set the applied current to silence all neurons apart from the neurons on the third board:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "neurocube  = NeuroCube()\n",
    "\n",
    "I_silent = [-1e-6, -1e-6, -1e-6, -1e-6] # all neurons silent\n",
    "I3 = [0e-6, 1e-6, 2e-6, 3e-6] # all neurons spiking\n",
    "Itotal = I_silent * 2 + I3 + I_silent # set the neurons of board 3 to spike\n",
    "\n",
    "Iapp = lambda t: fit_obj.convert_I(Itotal)\n",
    "\n",
    "x0 = [0.7, 0, 0, 0] * 16 + [0]*48 # 16 neurons, 48 synapses\n",
    "\n",
    "sol = neurocube.simulate(trange, x0, Iapp)\n",
    "\n",
    "for i in range(4):\n",
    "    V1 = sol.y[16 * i + 0]\n",
    "    V2 = sol.y[16 * i + 4]\n",
    "    V3 = sol.y[16 * i + 8]\n",
    "    V4 = sol.y[16 * i + 12]\n",
    "\n",
    "    plt.figure()\n",
    "    plt.title(f\"NeuroDyn Board {i + 1}\")\n",
    "    plt.xlabel('t [s]')\n",
    "    plt.ylabel('V [V]')\n",
    "    plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "    plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "    plt.plot(sol.t, V3, label = \"Neuron 3\")\n",
    "    plt.plot(sol.t, V4, label = \"Neuron 4\")\n",
    "    plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we would like to change the parameters of an individual board we can use the method `get_board(i)` which returns the board with the index $i$, after which we can use the methods described in **NeuroDynBoard** section. For example:\n",
    "\n",
    "`board1 = neurocube.get_board(0)`\n",
    "\n",
    "An additional method that the `NeuroCube` class provides is `connect_boards(board_i, board_j, neuron_no, g)` which sets the resistive connection between neuron numbered `neuron_no` on board `board_i` and neuron numbered `neuron_no` on board `board_j` with conductance `g`.\n",
    "\n",
    "Let's for example connect boards 3 and 1 through neuron 1, and boards 3 and 4 through neuron 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 1e-8\n",
    "neurocube.connect_boards(2, 0, 0, g) # boards 3 and 1 connected through neuron 1\n",
    "neurocube.connect_boards(2, 3, 1, g) # boards 3 and 4 connected through neuron 2\n",
    "\n",
    "sol = neurocube.simulate(trange, x0, Iapp)\n",
    "\n",
    "for i in range(4):\n",
    "    V1 = sol.y[16 * i + 0]\n",
    "    V2 = sol.y[16 * i + 4]\n",
    "    V3 = sol.y[16 * i + 8]\n",
    "    V4 = sol.y[16 * i + 12]\n",
    "\n",
    "    plt.figure()\n",
    "    plt.title(f\"NeuroDyn board {i + 1}\")\n",
    "    plt.xlabel('t [s]')\n",
    "    plt.ylabel('V [V]')\n",
    "    plt.plot(sol.t, V1, label = \"Neuron 1\")\n",
    "    plt.plot(sol.t, V2, label = \"Neuron 2\")\n",
    "    plt.plot(sol.t, V3, label = \"Neuron 3\")\n",
    "    plt.plot(sol.t, V4, label = \"Neuron 4\")\n",
    "    plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}