wb_model_D1_sto_det.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Project introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project we will show how to simulate the whole-brain mouse dynamics, loosely following [1].\n",
    "Using the structural connectome obtained from diffusion-weighted imaing, we will build a network model of a mouse brain and simulate the resting state dynamics, which we then compare to the fMRI (functional magnetic resonance imaging) from the individual mouse. We are also going to assess the brain lateralization of the dataset in [1] by both an empirical investigation and an in-silico assessment with brain network modelling. The objective is to employ TVB to describe the structure-function interaction by mathematically interfering with the propertie of brain asymmetry and evaluate the variability of this interaction.\n",
    "\n",
    "For this project we will be using Reduced Wong Wang (RWW) model [2,3]. In this first part, we will explore its behavior in a single node, unconnected to the network.\n",
    "\n",
    "The model is described by a single differential equation for each node,\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    \\dfrac{\\text{d}S_{i}}{\\text{d}t} &=& \\dfrac{-S_{i}}{\\tau_{s}} + \\gamma \\ (1 - S_{i})\\ H(x_{i}) + \\sigma\\eta_{i}(t)\\\\\n",
    "    &\\\\\n",
    "    H(x_{i}) &=& \\dfrac{ax_{i} - b}{1 - \\exp(-d \\ (ax_{i} - b))}\\\\\n",
    "    &\\\\\n",
    "    x_{i} &=& wJ_{N}S_{i} + J_{N}G\\sum_{j}C_{ij}S_{j} + I_{0}\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $S_{i}$ denotes the average synaptic gating variable at the local area $i$, $H(x_{i})$ is a sigmoid function that converts the input synaptic activity $x_{i}$ into an output population firing rate. $a = 0.270$ (nA.ms<sup>-1</sup>), $b = 0.108$ (kHz), and $d = 154$ (ms) are the parameters of the input-output function $H$, and $w$ is the local excitatory recurrence that is equal to 1. The kinetic parameters are $\\gamma = 0.641$, and $\\tau=100$ ms. The synaptic couplings are $J_{N} = 0.2609$ (nA) and the overall effective external input is $I_0 = 0.3$ (nA).\n",
    "\n",
    "$C_{ij}$ are the entries of the anatomical structural connectivity matrix reweigthed by the global coupling parameter $G$. $\\eta_{i}(t)$ is a Gaussian white noise with a noise amplitude $\\sigma = 0.001.$\n",
    "\n",
    "\n",
    "**References**\n",
    "\n",
    "[1] Melozzi, Francesca, Eyal Bergmann, Julie A. Harris, Itamar Kahn, Viktor Jirsa, and Christophe Bernard. 2019. “Individual Structural Features Constrain the Mouse Functional Connectome.” Proceedings of the National Academy of Sciences 116 (52): 26961–69. https://doi.org/10.1073/pnas.1906694116.\n",
    "\n",
    "[2] Wong, Kong-Fatt, and Xiao-Jing Wang. 2006. “A Recurrent Network Mechanism of Time Integration in Perceptual Decisions.” Journal of Neuroscience 26 (4): 1314–28. https://doi.org/10.1523/JNEUROSCI.3733-05.2006.\n",
    "\n",
    "[3] Deco, Gustavo, Adrián Ponce-Alvarez, Dante Mantini, Gian Luca Romani, Patric Hagmann, and Maurizio Corbetta. 2013. “Resting-State Functional Connectivity Emerges from Structurally and Dynamically Shaped Slow Linear Fluctuations.” Journal of Neuroscience 33 (27): 11239–52. https://doi.org/10.1523/JNEUROSCI.1091-13.2013.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 1: Exploring single node dynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goals:\n",
    "- Plot the RHS function of the RWW model for $S \\in [0, 1]$ and parameter values: \n",
    "    - $I_o$ = 0.300 and $\\omega$ = 1.0  \n",
    "    - $I_o$ = 0.322 and $\\omega$ = 1.0\n",
    "    - $I_o$ = 0.350 and $\\omega$ = 1.0\n",
    "  \n",
    "  Keep all other parameters default. What does it tell you about the fixed points of the model?\n",
    "- In a 1-node network with Reduced Wong Wang model ($I_o$ = 0.322, $\\omega$ = 1.0), run deterministic simulations with initial conditions for $S$ ranging from 0 to 1 and plot the state over time. Compare the results with the previous point.\n",
    "  - Run 3 seconds long simulations, use `EulerDeterministic` scheme with `dt` = 0.1.\n",
    "  - Use `Monitors.Raw()`  or `Monitors.TemporalAverage(period=1)` to monitor the state space.\n",
    "  - Use the 1-node connectivity defined below.\n",
    "- In the same network, set the initial conditions to 0.1 and run stochastic simulations with additive noise and noise levels from 0.001 to 0.3. How does the behavior change with increasing noise?\n",
    "  - Run 60 seconds long simulations.\n",
    "    \n",
    "Bonus goals:\n",
    "- Find the fixed points of the models numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-09-03 15:39:21,322 - WARNING - tvb.simulator.common - psutil module not available: no warnings will be issued when a\n",
      "    simulation may require more memory than available\n",
      "   INFO  log level set to INFO\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mariolavanga/Documents/GitHub/eitn_school_2021_tvb_project/env/lib/python3.7/site-packages/tvb/datatypes/surfaces.py:61: UserWarning: Geodesic distance module is unavailable; some functionality for surfaces will be unavailable.\n",
      "  warnings.warn(msg)\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from tvb.simulator.lab import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "conn = connectivity.Connectivity(weights=np.array([[0]]), \n",
    "                                 tract_lengths=np.array([[0]]),\n",
    "                                 region_labels=np.array([\"Region\"]),\n",
    "                                 centres=np.array([[0, 0, 0]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RHS function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Deterministic simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulations:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "tags": []
   },
   "source": [
    "dt = 0.1\n",
    "\n",
    "results = []\n",
    "ics = np.linspace(0, 1, 11)\n",
    "\n",
    "for i, ic in enumerate(ics):\n",
    "    rww = models.ReducedWongWang(I_o=np.array(0.322), w=np.array(1.0))\n",
    "    ...\n",
    "    \n",
    "    results.append((t, raw[:, 0, 0, 0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the results:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stochastic simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulations:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "tags": []
   },
   "source": [
    "dt = 0.1\n",
    "\n",
    "nsigmas = [0.001, 0.007, 0.015, 0.03]\n",
    "results = []\n",
    "CPU_TOTAL = 0\n",
    "for i, nsigma in enumerate(nsigmas):\n",
    "    t0 = time.time()\n",
    "    print(['SIGMA IDX -->',str(i)])\n",
    "    rww = models.ReducedWongWang(I_o=np.array(0.322), w=np.array(1.0))    \n",
    "    ...\n",
    "    results.append((t, raw[:, 0, 0, 0]))\n",
    "    \n",
    "    CPU_TIME = time.time() - t0\n",
    "    print(f\"CPU TIME = {CPU_TIME:.5f}\")\n",
    "    CPU_TOTAL += CPU_TIME\n",
    "    \n",
    "print(f\"TOTAL CPU TIME = {CPU_TOTAL:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}