{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Adjusting the low-gamma peak of microcircuit model\n\nHere we calculate the sensitivity measure of the :cite:t:`potjans2014`\nmicrocircuit model including modifications made in :cite:t:`bos2016` to study\nthe low-$\\gamma$ peak around $63\\,\\mathrm{Hz}$ in the power\nspectrum of the network. We subsequently adjust the network connectivity to\nmodify the peak.\n\nThis example reproduces parts of Fig. 5 and 8 in :cite:t:`bos2016`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nnmt\nimport numpy as np\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nfrom mpl_toolkits.axes_grid1 import make_axes_locatable\nfrom mpl_toolkits.axes_grid1.inset_locator import InsetPosition\nimport matplotlib.ticker\n\n\ndef colorbar(mappable, cax=None):\n    ax = mappable.axes\n    fig = ax.figure\n    divider = make_axes_locatable(ax)\n    if cax==None:\n        cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n    return fig.colorbar(mappable, cax=cax)\n\nplt.style.use('frontiers.mplstyle')\nmpl.rcParams.update({'legend.fontsize': 'medium',  # old: 5.0 was too small\n                     'axes.titlepad': 0.0})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we create an instance of the network model class `Microcircuit` using\nthe parameters from :cite:t:`bos2016`\n(:download:`Bos2016_network_params.yaml <../../../../tests/fixtures/integration/config/Bos2016_network_params.yaml>`,\n:download:`Bos2016_analysis_params.yaml <../../../../tests/fixtures/integration/config/Bos2016_analysis_params.yaml>`).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "microcircuit = nnmt.models.Microcircuit(\n    network_params=\n    '../../tests/fixtures/integration/config/Bos2016_network_params.yaml',\n    analysis_params=\n    '../../tests/fixtures/integration/config/Bos2016_analysis_params.yaml')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The frequency resolution used in the original publication was quite high.\nHere, we reduce the frequency resolution for faster execution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reduce_frequency_resolution = False\n\nif reduce_frequency_resolution:\n    microcircuit.change_parameters(changed_analysis_params={'df': 1},\n                                overwrite=True)\n    derived_analysis_params = (\n        microcircuit._calculate_dependent_analysis_parameters())\n    microcircuit.analysis_params.update(derived_analysis_params)\n\nfrequencies = microcircuit.analysis_params['omegas']/(2.*np.pi)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We calculate all required dynamical properties, the power spectra, and\nfinally the sensitivity measure for all eigenmodes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# calculate working point for exponentially shape post synaptic currents\nnnmt.lif.exp.working_point(microcircuit, method='taylor')\n# calculate the transfer function\nnnmt.lif.exp.transfer_function(microcircuit, method='taylor')\n# calculate the delay distribution matrix\nnnmt.network_properties.delay_dist_matrix(microcircuit)\n# calculate the effective connectivity matrix\nnnmt.lif.exp.effective_connectivity(microcircuit)\n\n# calculate the power spectra\npower_spectra = nnmt.lif.exp.power_spectra(microcircuit)\n\n# calculate sensitivity measure\nsensitivity_dict = nnmt.lif.exp.sensitivity_measure_all_eigenmodes(\n    microcircuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The sensitivity measure tells us that the low-$\\gamma$ peak in the\npower spectra can be modified by altering the connection from population 4I\nto itself. Here we change this connection by 5 percent. In order to keep the\nnetwork's working point constant when changing the connection, we modify the\nconnections from the external input to the target population as well (use Eq.\n2 in :cite:t:`bos2016` for this).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "percentage_of_change = 0.05\n\ng =  microcircuit.network_params['g']\npopulation_rate = microcircuit.results['lif.exp.firing_rates'][3]\nexternal_input_rate = microcircuit.network_params['nu_ext']\nconnections_4I_4I = microcircuit.network_params['K'][3,3]\nconnection_external_4I = microcircuit.network_params['K_ext'][3]\n\nexternal_connections_to_add = (\n    abs(g) * connections_4I_4I * percentage_of_change * population_rate\n    ) / external_input_rate\n\nK_5_percent = microcircuit.network_params['K'].copy()\nK_5_percent_ext = microcircuit.network_params['K_ext'].copy()\nK_5_percent[3,3] = K_5_percent[3,3]*(1+percentage_of_change)\nK_5_percent_ext[3] = K_5_percent_ext[3] + external_connections_to_add\n\nmicrocircuit_5_percent = microcircuit.change_parameters(\n    {'K': K_5_percent,\n     'K_ext': K_5_percent_ext})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After changing the network connectivity, we calculate the power spectra\nagain.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# calculate working point for exponentially shape post synaptic currents\nnnmt.lif.exp.working_point(microcircuit_5_percent, method='taylor')\n# calculate the transfer function\nnnmt.lif.exp.transfer_function(microcircuit_5_percent, method='taylor')\n# calculate the delay distribution matrix\nnnmt.network_properties.delay_dist_matrix(microcircuit_5_percent)\n# calculate the effective connectivity matrix\nnnmt.lif.exp.effective_connectivity(microcircuit_5_percent)\n# calculate the power spectra\npower_spectra_5_percent = nnmt.lif.exp.power_spectra(microcircuit_5_percent)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we change the connectivity by 10 percent, similar to above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "percentage_of_change = 0.10\n\nexternal_connections_to_add = (\n    abs(g) * connections_4I_4I * percentage_of_change * population_rate\n    ) / external_input_rate\n\nK_10_percent = microcircuit.network_params['K'].copy()\nK_10_percent_ext = microcircuit.network_params['K_ext'].copy()\nK_10_percent[3,3] = K_10_percent[3,3]*(1+percentage_of_change)\nK_10_percent_ext[3] = K_10_percent_ext[3] + external_connections_to_add\n\nmicrocircuit_10_percent = microcircuit.change_parameters(\n    {'K': K_10_percent,\n     'K_ext': K_10_percent_ext})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After changing the network connectivity, we calculate the power spectra\nagain.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# calculate working point for exponentially shape post synaptic currents\nnnmt.lif.exp.working_point(microcircuit_10_percent, method='taylor')\n# calculate the transfer function\nnnmt.lif.exp.transfer_function(microcircuit_10_percent, method='taylor')\n# calculate the delay distribution matrix\nnnmt.network_properties.delay_dist_matrix(microcircuit_10_percent)\n# calculate the effective connectivity matrix\nnnmt.lif.exp.effective_connectivity(microcircuit_10_percent)\n# calculate the power spectra\npower_spectra_10_percent = nnmt.lif.exp.power_spectra(microcircuit_10_percent)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we print the critical frequencies for each eigenmode\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Look at the critical frequencies per eigenmode\nfor k, v in sensitivity_dict.items():\n    print(k, v['critical_frequency'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "and we plot the results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "width = 180. / 25.4\nheight = 60. / 25.4\n\nfig = plt.figure(figsize=(width, height),\n                 constrained_layout=True)\n\nlabels = ['2/3E', '2/3I', '4E', '4I', '5E', '5I', '6E', '6I']\n\ncolormap = mpl.cm.get_cmap('coolwarm').copy()\n\n# set colorbar max and min\nz = 1\n\n# by default np.nan are set to black with full transparency\n# .eps can't handle transparency\ncolormap.set_bad('w',1.)\n\nfig = plt.figure(figsize=(width, height),\n                 constrained_layout=True)\ngrid_specification = gridspec.GridSpec(1, 2,\n                                       height_ratios=[1],\n                                       width_ratios=[2.2, 1], figure=fig)\n\ngs = gridspec.GridSpecFromSubplotSpec(1,3,\n                                        height_ratios=[1],\n                            width_ratios=[1, 1, 0.2],\n                            subplot_spec=grid_specification[0])\n\nev = '5'\nax = fig.add_subplot(gs[0])\nlabel_prms = dict(x=-0.3, y=1.27, fontsize=10, fontweight='bold',\n                  va='top', ha='right')\n\nax.text(s='(A)', transform=ax.transAxes, **label_prms)\n\nfrequency = sensitivity_dict[ev]['critical_frequency']\nprojection_of_sensitivity_measure = sensitivity_dict[ev][\n    'sensitivity_amp']\n\nrounded_frequency = str(int(np.round(frequency,0)))\n\nplot_title = r'$\\mathbf{Z}_{b=%s}^{\\mathrm{amp}}(' % ev + \\\n    f'{rounded_frequency}' + r'\\,\\mathrm{Hz})$'\nax.set_title(plot_title)\n\ndata = np.ma.masked_where(projection_of_sensitivity_measure == 0,\n                            projection_of_sensitivity_measure)\n\n# np.flipud needed to cope for different origin compared to imshow\nheatmap = ax.pcolormesh(np.flipud(data),\n                    vmin=-z,\n                    vmax=z,\n                    cmap=colormap,\n                    edgecolors='k',\n                    linewidth=0.6)\n\nax.set_aspect('equal')\n\nif labels is not None:\n    ax.set_xticks(np.arange(len(labels))+0.5)\n    ax.set_yticks(np.arange(len(labels))+0.5)\n    ax.set_xticklabels(labels, rotation=90)\n    ax.set_yticklabels(list(reversed(labels)))\n\nax.set_xlabel('sources')\nax.set_ylabel('targets')\n\n\n# sensitivity_measure_frequency\nax = fig.add_subplot(gs[1])\n\nfrequency = sensitivity_dict[ev]['critical_frequency']\nprojection_of_sensitivity_measure = sensitivity_dict[ev][\n    'sensitivity_freq']\n\nrounded_frequency = str(int(np.round(frequency,0)))\n\nplot_title = r'$\\mathbf{Z}_{b=%s}^{\\mathrm{freq}}(' % ev + \\\n    f'{rounded_frequency}' + r'\\,\\mathrm{Hz})$'\nax.set_title(plot_title)\n\ndata = np.ma.masked_where(projection_of_sensitivity_measure == 0,\n                            projection_of_sensitivity_measure)\n\nheatmap = ax.pcolormesh(np.flipud(data),\n                    vmin=-z,\n                    vmax=z,\n                    cmap=colormap,\n                    edgecolors='k',\n                    linewidth=0.6)\n\nax.set_aspect('equal')\n\nif labels is not None:\n    ax.set_xticks(np.arange(len(labels))+0.5)\n    ax.set_yticks(np.arange(len(labels))+0.5)\n    ax.set_xticklabels(labels, rotation=90)\n    ax.set_yticklabels([])\n\nax.set_xlabel('sources')\n\ncolorbar_ax = fig.add_subplot(gs[2])\n\ncolorbar_width = 0.1\nip = InsetPosition(ax, [1.05,0,colorbar_width,1])\ncolorbar_ax.set_axes_locator(ip)\ncolorbar(heatmap, cax=colorbar_ax)\n\nax = fig.add_subplot(grid_specification[1])\nlabel_prms = dict(x=-0.2, y=1.2, fontsize=10, fontweight='bold',\n                  va='top', ha='right')\nax.text(s='(B)', transform=ax.transAxes, **label_prms)\nj = 3\nax.plot(microcircuit.analysis_params['omegas']/(2*np.pi),\n        power_spectra[:, j],\n        color='#4055C8', zorder=2,\n        label='original'\n        )\n\nax.plot(microcircuit_5_percent.analysis_params['omegas']/(2*np.pi),\n        power_spectra_5_percent[:, j],\n        color='#6072D5', zorder=2, ls='dashed',\n        label='5% increase'\n        )\n\nax.plot(microcircuit_10_percent.analysis_params['omegas']/(2*np.pi),\n        power_spectra_10_percent[:, j],\n        color='#8A98E5', zorder=2, ls='dotted',\n        label='10% increase'\n        )\n\nax.set_xlim([0.0, 100.0])\nax.set_ylim([1e-6, 1e0])\nax.set_yscale('log')\nax.set_xticks([20, 40, 60, 80])\n\nax.set_xlabel(r'frequency $\\omega/2\\pi\\quad(1/\\mathrm{s})$')\n\ny_minor = matplotlib.ticker.LogLocator(\n    base = 10.0,\n    subs = np.arange(1.0, 10.0) * 0.1,\n    numticks = 10)\nax.yaxis.set_minor_locator(y_minor)\nax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())\nax.set_yticks([])\nax.set_yticks([1e-5,1e-3,1e-1])\nax.set_ylabel(r'power spectrum $P_{\\mathrm{4I}}(\\omega)\\quad(1/\\mathrm{s}^2)$')\nax.legend(loc='lower right')\n\nplt.savefig('figures/sensitivity_measure_low_gamma_Bos2016.eps')"
      ]
    }
  ],
  "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.9.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}