plotting/trgc-MedOffOn.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Plot MED OFF vs. ON TRGC (Figure 3c & S5).\"\"\"\n",
    "\n",
    "import os\n",
    "import sys\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import copy\n",
    "import pandas as pd\n",
    "import pte_stats\n",
    "import pte_decode\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import sem\n",
    "\n",
    "cd_path = Path(os.getcwd()).absolute().parent\n",
    "sys.path.append(os.path.join(cd_path, \"coherence\"))\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.rc('xtick', labelsize=6)\n",
    "matplotlib.rc('ytick', labelsize=6)\n",
    "matplotlib.rc('legend', fontsize=6)\n",
    "matplotlib.rc(\"font\", size=6, family=\"Arial\")\n",
    "matplotlib.rc('axes', labelsize=7)\n",
    "matplotlib.rc('axes', titlesize=7)\n",
    "matplotlib.rcParams[\"pdf.fonttype\"] = 42\n",
    "matplotlib.rcParams[\"ps.fonttype\"] = 42\n",
    "\n",
    "prop_cycle = plt.rcParams[\"axes.prop_cycle\"]\n",
    "\n",
    "n_perm = 100_000\n",
    "two_sided = True\n",
    "alpha = 0.05\n",
    "\n",
    "FOLDERPATH_ANALYSIS = \"Path_to\\\\Project\\\\Analysis\"\n",
    "FOLDERPATH_FIGURES = os.path.join(os.path.dirname(os.getcwd()), \"figures\")\n",
    "\n",
    "lfreq = 4.0\n",
    "hfreq = 40.0\n",
    "\n",
    "fband_func = np.mean\n",
    "\n",
    "subregions_mapping = {\n",
    "    \"Motor\": \"motor\",\n",
    "    \"Sensory\": \"sensory\",\n",
    "}\n",
    "\n",
    "fbands = {\n",
    "    \"alpha\": (8, 12),\n",
    "    \"low_beta\": (12, 20),\n",
    "    \"high_beta\": (20, 30)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MOTOR CORTEX -> STN TRGC (Figure 3c)\n",
    "data = pd.read_pickle(os.path.join(\n",
    "    FOLDERPATH_ANALYSIS, \"task-Rest_acq-multi_run-multi_con_granger_regional-MedOffOn_multi_sub.pkl\")\n",
    ")\n",
    "freqs = copy.deepcopy(data[\"frequencies\"][0])\n",
    "freqs_delta = np.unique(np.subtract(freqs[1:], freqs[:-1]))[0]\n",
    "data[\"frequencies\"] = [freqs] * len(data[\"seed_types\"])\n",
    "data = pd.DataFrame.from_dict(data)\n",
    "\n",
    "trgc_off = np.array(data[\"connectivity-trgc\"][(data[\"med\"] == \"Off\") & (data[\"seed_subregions\"] == \"motor\")].to_list())\n",
    "trgc_on = np.array(data[\"connectivity-trgc\"][(data[\"med\"] == \"On\") & (data[\"seed_subregions\"] == \"motor\")].to_list())\n",
    "\n",
    "np.random.seed(44)\n",
    "\n",
    "lfreq_i = freqs.index(lfreq)\n",
    "hfreq_i = freqs.index(hfreq)\n",
    "plot_freqs = freqs[lfreq_i:hfreq_i+1]\n",
    "\n",
    "trgc_off = trgc_off[:, lfreq_i:hfreq_i+1]\n",
    "trgc_on = trgc_on[:, lfreq_i:hfreq_i+1]\n",
    "\n",
    "trgc_off_sem = sem(trgc_off)\n",
    "trgc_on_sem = sem(trgc_on)\n",
    "\n",
    "fig, axis = pte_decode.lineplot_prediction_compare(\n",
    "    x_1=trgc_off.T,\n",
    "    x_2=trgc_on.T,\n",
    "    times=np.arange(lfreq, hfreq+0.5, 0.5),\n",
    "    data_labels=[\"OFF therapy\", \"ON levodopa\"],\n",
    "    x_label=\"Frequency (Hz)\",\n",
    "    y_label=\"Time-reversed\\nGranger causality (A.U.)\",\n",
    "    two_tailed=two_sided,\n",
    "    paired_x1x2=True,\n",
    "    n_perm=n_perm,\n",
    "    correction_method=\"cluster_pvals\",\n",
    "    title=\"Motor cortex -> STN\",\n",
    "    colour=[\"#DF4A4A\", \"#376EB4\"]\n",
    ")\n",
    "\n",
    "for fband_name, fband_lims in fbands.items():\n",
    "    stat, pval = pte_stats.permutation_onesample(\n",
    "        data_a=fband_func(trgc_off[:, freqs.index(fband_lims[0]):freqs.index(fband_lims[1])+1], 1),\n",
    "        data_b=fband_func(trgc_on[:, freqs.index(fband_lims[0]):freqs.index(fband_lims[1])+1], 1),\n",
    "        n_perm=n_perm, two_tailed=two_sided\n",
    "    )\n",
    "    if pval < alpha:\n",
    "        axis.axvspan(fband_lims[0], fband_lims[1], color=[0.8, 0.8, 0.8], alpha=0.3, linewidth=0)\n",
    "        prefix = \"[SIGNIFICANT]\"\n",
    "    else:\n",
    "        prefix = \"\"\n",
    "    print(f\"{prefix} {fband_name}: stat={stat}; p-val={pval}\")\n",
    "\n",
    "axis.spines['top'].set_visible(False)\n",
    "axis.spines['right'].set_visible(False)\n",
    "axis.legend(handlelength=0.5, handletextpad=0.5, labelspacing=0, borderaxespad=0, loc=\"upper left\")\n",
    "leg = axis.get_legend()\n",
    "leg.legend_handles[0].set_linewidth(1.5)\n",
    "leg.legend_handles[1].set_linewidth(1.5)\n",
    "leg.set_frame_on(False)\n",
    "axis.set_xticks(np.arange(5, 45, 5))\n",
    "axis.tick_params(\"x\", pad=2, size=2)\n",
    "axis.tick_params(\"y\", pad=2, size=2)\n",
    "axis.xaxis.labelpad = 0\n",
    "axis.yaxis.labelpad = 1\n",
    "axis.set_yticks([0, 0.05, 0.1])\n",
    "\n",
    "fig.set_size_inches(1.8, 1.35)\n",
    "\n",
    "fig.savefig(os.path.join(FOLDERPATH_FIGURES, \"Manuscript_TRGC_Med_motor_cortex-STN.pdf\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SENSORY CORTEX -> STN TRGC (Figure S5)\n",
    "data = pd.read_pickle(os.path.join(\n",
    "    FOLDERPATH_ANALYSIS, \"task-Rest_acq-multi_run-multi_con_granger_regional-MedOffOn_multi_sub.pkl\")\n",
    ")\n",
    "freqs = copy.deepcopy(data[\"frequencies\"][0])\n",
    "freqs_delta = np.unique(np.subtract(freqs[1:], freqs[:-1]))[0]\n",
    "data[\"frequencies\"] = [freqs] * len(data[\"seed_types\"])\n",
    "data = pd.DataFrame.from_dict(data)\n",
    "\n",
    "trgc_off = np.array(data[\"connectivity-trgc\"][(data[\"med\"] == \"Off\") & (data[\"seed_subregions\"] == \"sensory\")].to_list())\n",
    "trgc_on = np.array(data[\"connectivity-trgc\"][(data[\"med\"] == \"On\") & (data[\"seed_subregions\"] == \"sensory\")].to_list())\n",
    "\n",
    "np.random.seed(44)\n",
    "\n",
    "lfreq_i = freqs.index(lfreq)\n",
    "hfreq_i = freqs.index(hfreq)\n",
    "plot_freqs = freqs[lfreq_i:hfreq_i+1]\n",
    "\n",
    "trgc_off = trgc_off[:, lfreq_i:hfreq_i+1]\n",
    "trgc_on = trgc_on[:, lfreq_i:hfreq_i+1]\n",
    "\n",
    "trgc_off_sem = sem(trgc_off)\n",
    "trgc_on_sem = sem(trgc_on)\n",
    "\n",
    "fig, axis = pte_decode.lineplot_prediction_compare(\n",
    "    x_1=trgc_off.T,\n",
    "    x_2=trgc_on.T,\n",
    "    times=np.arange(lfreq, hfreq+0.5, 0.5),\n",
    "    data_labels=[\"OFF therapy\", \"ON levodopa\"],\n",
    "    x_label=\"Frequency (Hz)\",\n",
    "    y_label=\"Time-reversed\\nGranger causality (A.U.)\",\n",
    "    two_tailed=two_sided,\n",
    "    paired_x1x2=True,\n",
    "    n_perm=n_perm,\n",
    "    correction_method=\"cluster_pvals\",\n",
    "    title=\"Sensory cortex -> STN\",\n",
    "    colour=[\"#DF4A4A\", \"#376EB4\"]\n",
    ")\n",
    "\n",
    "for fband_name, fband_lims in fbands.items():\n",
    "    stat, pval = pte_stats.permutation_onesample(\n",
    "        data_a=fband_func(trgc_off[:, freqs.index(fband_lims[0]):freqs.index(fband_lims[1])+1], 1),\n",
    "        data_b=fband_func(trgc_on[:, freqs.index(fband_lims[0]):freqs.index(fband_lims[1])+1], 1),\n",
    "        n_perm=n_perm, two_tailed=two_sided\n",
    "    )\n",
    "    if pval < alpha:\n",
    "        axis.axvspan(fband_lims[0], fband_lims[1], color=[0.8, 0.8, 0.8], alpha=0.3, linewidth=0)\n",
    "        prefix = \"[SIGNIFICANT]\"\n",
    "    else:\n",
    "        prefix = \"\"\n",
    "    print(f\"{prefix} {fband_name}: stat={stat}; p-val={pval}\")\n",
    "\n",
    "axis.spines['top'].set_visible(False)\n",
    "axis.spines['right'].set_visible(False)\n",
    "axis.legend(handlelength=0.5, handletextpad=0.5, labelspacing=0, borderaxespad=0, loc=\"upper left\")\n",
    "leg = axis.get_legend()\n",
    "leg.legend_handles[0].set_linewidth(1.5)\n",
    "leg.legend_handles[1].set_linewidth(1.5)\n",
    "leg.set_frame_on(False)\n",
    "axis.set_xticks(np.arange(5, 45, 5))\n",
    "axis.tick_params(\"x\", pad=2, size=2)\n",
    "axis.tick_params(\"y\", pad=2, size=2)\n",
    "axis.xaxis.labelpad = 0\n",
    "axis.yaxis.labelpad = 1\n",
    "axis.set_yticks([0, 0.04, 0.08])\n",
    "\n",
    "fig.set_size_inches(1.8, 1.35)\n",
    "\n",
    "fig.savefig(os.path.join(FOLDERPATH_FIGURES, \"Manuscript_TRGC_Med_sensory_cortex-STN.pdf\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find significant periods of connectivity (whole cortex -> STN; OFF)\n",
    "data = pd.read_pickle(os.path.join(\n",
    "    FOLDERPATH_ANALYSIS, \"task-Rest_acq-multi_run-multi_con_granger_whole-MedOffOn_multi_sub.pkl\")\n",
    ")\n",
    "freqs = copy.deepcopy(data[\"frequencies\"][0])\n",
    "freqs_delta = np.unique(np.subtract(freqs[1:], freqs[:-1]))[0]\n",
    "data[\"frequencies\"] = [freqs] * len(data[\"seed_types\"])\n",
    "data = pd.DataFrame.from_dict(data)\n",
    "\n",
    "lfreq_i = freqs.index(lfreq)\n",
    "hfreq_i = freqs.index(100)\n",
    "plot_freqs = freqs[lfreq_i:hfreq_i+1]\n",
    "\n",
    "trgc_off = np.array(data[\"connectivity-net_trgc\"][data[\"med\"] == \"Off\"].to_list())\n",
    "trgc_on = np.array(data[\"connectivity-net_trgc\"][data[\"med\"] == \"On\"].to_list())\n",
    "trgc = np.mean([trgc_off, trgc_on], axis=0)\n",
    "trgc = trgc[:, lfreq_i:hfreq_i+1]\n",
    "\n",
    "fig, axis = pte_decode.lineplot_prediction_compare(\n",
    "    x_1=trgc.T,\n",
    "    x_2=np.zeros_like(trgc).T,\n",
    "    times=np.array(plot_freqs),\n",
    "    data_labels=[\"TRGC\", \"Baseline\"],\n",
    "    x_label=\"Frequency (Hz)\",\n",
    "    y_label=\"TRGC (A.U.)\",\n",
    "    two_tailed=two_sided,\n",
    "    paired_x1x2=True,\n",
    "    n_perm=n_perm,\n",
    "    correction_method=\"cluster_pvals\",\n",
    "    title=\"Cortex -> STN (OFF)\",\n",
    "    colour=[\"#DF4A4A\", \"k\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find significant periods of connectivity (motor cortex -> STN; OFF)\n",
    "data = pd.read_pickle(os.path.join(\n",
    "    FOLDERPATH_ANALYSIS, \"task-Rest_acq-multi_run-multi_con_granger_regional-MedOffOn_multi_sub.pkl\")\n",
    ")\n",
    "freqs = copy.deepcopy(data[\"frequencies\"][0])\n",
    "freqs_delta = np.unique(np.subtract(freqs[1:], freqs[:-1]))[0]\n",
    "data[\"frequencies\"] = [freqs] * len(data[\"seed_types\"])\n",
    "data = pd.DataFrame.from_dict(data)\n",
    "\n",
    "lfreq_i = freqs.index(lfreq)\n",
    "hfreq_i = freqs.index(hfreq)\n",
    "plot_freqs = freqs[lfreq_i:hfreq_i+1]\n",
    "\n",
    "trgc_off = np.array(data[\"connectivity-trgc\"][(data[\"med\"] == \"Off\") & (data[\"seed_subregions\"] == \"motor\")].to_list())\n",
    "trgc_off = trgc_off[:, lfreq_i:hfreq_i+1]\n",
    "\n",
    "fig, axis = pte_decode.lineplot_prediction_compare(\n",
    "    x_1=trgc_off.T,\n",
    "    x_2=np.zeros_like(trgc_off).T,\n",
    "    times=np.array(plot_freqs),\n",
    "    data_labels=[\"TRGC\", \"Baseline\"],\n",
    "    x_label=\"Frequency (Hz)\",\n",
    "    y_label=\"TRGC (A.U.)\",\n",
    "    two_tailed=two_sided,\n",
    "    paired_x1x2=True,\n",
    "    n_perm=n_perm,\n",
    "    correction_method=\"cluster_pvals\",\n",
    "    title=\"Cortex -> STN (OFF)\",\n",
    "    colour=[\"#DF4A4A\", \"k\"]\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pte",
   "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.10.9"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}