{
 "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 2 (Smooth pursuit eye movements)\n",
    "This is an idealized model of the smooth pursuit reflex, including two ocular muscles, a moving visual stimulus and spiking neural control."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from brian2 import *\n",
    "seed(79620)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the two antagonistic ocular muscles is modelled as a spring of elasticity $k$ and some friction. To simplify, we consider that the eye moves laterally, rather than rotate. If $x$ is the position of the eye with 0 being the center, then the lengths of the springs are $L+x$ and $L-x$. The dynamics of the eye is then given by a second-order differential equation:\n",
    "$$ m\\frac{d^2x}{dt^2} = - k\\left(\\left(L+x\\right)-x_L\\right) + k\\left(\\left(L-x\\right)-x_R\\right) - f\\frac{dx}{dt}$$\n",
    "or:\n",
    "$$ m\\frac{d^2x}{dt^2} = k(x_L-x_R-2x) - f\\frac{dx}{dt}$$\n",
    "\n",
    "We see that $x_0 = \\frac{1}{2}(x_L-x_R)$ is the equilibrium position of the eye. Here we have assumed that spring elasticities are identical. We can rewrite this equation with just two parameters $\\alpha$ and $\\beta$:\n",
    "$$ \\frac{d^2x}{dt^2} = \\alpha(x_0 - x) - \\beta\\frac{dx}{dt}$$\n",
    "\n",
    "We will assume that eye position can move between -1 and 1.\n",
    "\n",
    "We consider that the resting length is the variable on which motoneurons act. Each spike from a motoneuron produces a waveform of contraction (\"twitch\"), i.e., a change in resting length. We consider that contractions add linearly, and resting lengths relax exponentially. By linearity it follows that we can simply express the action of motoneurons on the equilibrium position of the eye $x_0$, which relaxes exponentially to the center position 0.\n",
    "\n",
    "Finally, we consider a visual object that performs a random walk according to an Ornstein–Uhlenbeck with 0 as the central location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = (1/(50*ms))**2 # characteristic relaxation time is 50 ms\n",
    "beta = 1/(50*ms) # friction parameter\n",
    "tau_muscle = 20*ms # relaxation time of muscle contraction\n",
    "tau_object = 500*ms # time constant of object movement\n",
    "\n",
    "eqs_eye = '''\n",
    "dx/dt = velocity : 1\n",
    "dvelocity/dt = alpha*(x0-x)-beta*velocity : 1/second\n",
    "dx0/dt = -x0/tau_muscle : 1\n",
    "dx_object/dt = (noise - x_object)/tau_object:  1\n",
    "dnoise/dt = -noise/tau_object + tau_object**-0.5*xi : 1\n",
    "'''\n",
    "eye = NeuronGroup(1, model=eqs_eye, method='euler')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define two motoneurons, one for each muscle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taum = 20*ms\n",
    "motoneurons = NeuronGroup(2, model= 'dv/dt = -v/taum : 1', threshold = 'v>1',\n",
    "                          reset = 'v=0', refractory = 5*ms, method='exact')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The motoneurons project to the eye, and each spike produces a small contraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "motosynapses = Synapses(motoneurons, eye, model = 'w : 1', on_pre = 'x0+=w')\n",
    "motosynapses.connect() # connects all motoneurons to the eye\n",
    "motosynapses.w = [-0.5,0.5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now implement the sensory neurons, which we simplify by considering spiking neurons which directly respond to light, i.e they represent both photoreceptors and retinal ganglion cells. \n",
    "TODO: Talk about Gaussian dependence, et.c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20\n",
    "width = 2./N # width of receptive field\n",
    "gain = 4.\n",
    "eqs_retina = '''\n",
    "I = gain*exp(-((x_object-x_eye-x_neuron)/width)**2) : 1\n",
    "x_neuron : 1 (constant)\n",
    "x_object : 1 (linked) # position of the object\n",
    "x_eye : 1 (linked) # position of the eye\n",
    "dv/dt = (I-(1+gs)*v)/taum : 1\n",
    "gs : 1 # total synaptic conductance\n",
    "'''\n",
    "retina = NeuronGroup(N, model = eqs_retina, threshold = 'v>1', reset = 'v=0', method='exact')\n",
    "retina.v = 'rand()'\n",
    "retina.x_eye = linked_var(eye, 'x')\n",
    "retina.x_object = linked_var(eye, 'x_object')\n",
    "retina.x_neuron = '-1.0 + 2.0*i/(N-1)'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we connect sensory neurons to motoneurons. Sensory neurons on each hemifield connects to the corresponding motoneuron, with a strength that scales with eccentricity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensorimotor_synapses = Synapses(retina, motoneurons, model = 'w : 1 (constant)', on_pre = 'v+=w')\n",
    "sensorimotor_synapses.connect(j = 'int(x_neuron_pre > 0)')\n",
    "sensorimotor_synapses.w = '20*abs(x_neuron_pre)/N_pre'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We record the position of the eye, of the object, and spikes produced by the retina and motoneurons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = StateMonitor(eye, ('x', 'x0', 'x_object'), record = True)\n",
    "S_retina = SpikeMonitor(retina)\n",
    "S_motoneurons = SpikeMonitor(motoneurons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now can run the simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run(10*second, report='text')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from plotly import tools\n",
    "from plotly.offline import iplot, init_notebook_mode\n",
    "import plotly.graph_objs as go\n",
    "\n",
    "init_notebook_mode(connected=True)\n",
    "\n",
    "fig = tools.make_subplots(3, 1, specs=[[{'rowspan': 2}], [None], [{}]],\n",
    "                          shared_xaxes=True, print_grid=False)\n",
    "\n",
    "trace = go.Scatter(x=S_retina.t/second,\n",
    "                   y=S_retina.i,\n",
    "                   marker={'symbol': 'line-ns', 'line': {'width': 1, 'color':'black'}},\n",
    "                   mode='markers',\n",
    "                   name='retina',\n",
    "                   showlegend=False)\n",
    "fig.append_trace(trace, 1, 1)\n",
    "motoneuron_spikes = S_motoneurons.spike_trains()\n",
    "trace = go.Scatter(x=motoneuron_spikes[0]/second,\n",
    "                   y=np.ones(S_motoneurons.count[0])*N,\n",
    "                   marker={'symbol': 'line-ns', 'line': {'width': 1, 'color':'#1f77b4'},\n",
    "                          'color':'#1f77b4'},\n",
    "                   mode='markers',\n",
    "                   name='left motoneuron',\n",
    "                   showlegend=False)\n",
    "fig.append_trace(trace, 1, 1)\n",
    "trace = go.Scatter(x=motoneuron_spikes[1]/second,\n",
    "                   y=np.ones(S_motoneurons.count[1])*(N+1),\n",
    "                   marker={'symbol': 'line-ns', 'line': {'width': 1, 'color':'#ff7f03'},\n",
    "                           'color':'#ff7f03'},\n",
    "                   mode='markers',\n",
    "                   name='right motoneuron',\n",
    "                   showlegend=False)\n",
    "fig.append_trace(trace, 1, 1)\n",
    "\n",
    "trace = go.Scatter(x=M.t/second,\n",
    "                   y=M.x[0],\n",
    "                   mode='lines',\n",
    "                   line={'color': 'black'},\n",
    "                   name='eye')\n",
    "fig.append_trace(trace, 3, 1)\n",
    "trace = go.Scatter(x=M.t/second,\n",
    "                   y=M.x_object[0],\n",
    "                   mode='lines',\n",
    "                   line={'color': '#2ca02c'},\n",
    "                   name='object')\n",
    "fig.append_trace(trace, 3, 1)\n",
    "\n",
    "fig['layout'].update(xaxis1={'showline': False,\n",
    "                             'zeroline': False,\n",
    "                             'title': 'time (in s)'},\n",
    "                     yaxis1={'title': 'neuron index',\n",
    "                             'showticklabels': False,\n",
    "                             'showline': True},\n",
    "                     yaxis2={'tickmode': 'array',\n",
    "                             'ticktext': ['left', 'right'],\n",
    "                             'tickvals': [-1, 1],\n",
    "                             'range': [-1.05, 1.05],\n",
    "                             'zeroline': True,\n",
    "                             'showline': True}\n",
    "                     )\n",
    "iplot(fig)"
   ]
  }
 ],
 "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
}