Adjusting the low-gamma peak of microcircuit model

Here we calculate the sensitivity measure of the Potjans and Diesmann [2014] microcircuit model including modifications made in Bos et al. [2016] to study the low- peak around in the power spectrum of the network. We subsequently adjust the network connectivity to modify the peak.

This example reproduces parts of Fig. 5 and 8 in Bos et al. [2016].

import nnmt
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
import matplotlib.ticker


def colorbar(mappable, cax=None):
    ax = mappable.axes
    fig = ax.figure
    divider = make_axes_locatable(ax)
    if cax==None:
        cax = divider.append_axes("right", size="5%", pad=0.05)
    return fig.colorbar(mappable, cax=cax)

plt.style.use('frontiers.mplstyle')
mpl.rcParams.update({'legend.fontsize': 'medium',  # old: 5.0 was too small
                     'axes.titlepad': 0.0})

First, we create an instance of the network model class Microcircuit using the parameters from Bos et al. [2016] (Bos2016_network_params.yaml, Bos2016_analysis_params.yaml).

microcircuit = nnmt.models.Microcircuit(
    network_params=
    '../../tests/fixtures/integration/config/Bos2016_network_params.yaml',
    analysis_params=
    '../../tests/fixtures/integration/config/Bos2016_analysis_params.yaml')

The frequency resolution used in the original publication was quite high. Here, we reduce the frequency resolution for faster execution.

reduce_frequency_resolution = False

if reduce_frequency_resolution:
    microcircuit.change_parameters(changed_analysis_params={'df': 1},
                                overwrite=True)
    derived_analysis_params = (
        microcircuit._calculate_dependent_analysis_parameters())
    microcircuit.analysis_params.update(derived_analysis_params)

frequencies = microcircuit.analysis_params['omegas']/(2.*np.pi)

We calculate all required dynamical properties, the power spectra, and finally the sensitivity measure for all eigenmodes.

# calculate working point for exponentially shape post synaptic currents
nnmt.lif.exp.working_point(microcircuit, method='taylor')
# calculate the transfer function
nnmt.lif.exp.transfer_function(microcircuit, method='taylor')
# calculate the delay distribution matrix
nnmt.network_properties.delay_dist_matrix(microcircuit)
# calculate the effective connectivity matrix
nnmt.lif.exp.effective_connectivity(microcircuit)

# calculate the power spectra
power_spectra = nnmt.lif.exp.power_spectra(microcircuit)

# calculate sensitivity measure
sensitivity_dict = nnmt.lif.exp.sensitivity_measure_all_eigenmodes(
    microcircuit)

The sensitivity measure tells us that the low- peak in the power spectra can be modified by altering the connection from population 4I to itself. Here we change this connection by 5 percent. In order to keep the network’s working point constant when changing the connection, we modify the connections from the external input to the target population as well (use Eq. 2 in Bos et al. [2016] for this).

percentage_of_change = 0.05

g =  microcircuit.network_params['g']
population_rate = microcircuit.results['lif.exp.firing_rates'][3]
external_input_rate = microcircuit.network_params['nu_ext']
connections_4I_4I = microcircuit.network_params['K'][3,3]
connection_external_4I = microcircuit.network_params['K_ext'][3]

external_connections_to_add = (
    abs(g) * connections_4I_4I * percentage_of_change * population_rate
    ) / external_input_rate

K_5_percent = microcircuit.network_params['K'].copy()
K_5_percent_ext = microcircuit.network_params['K_ext'].copy()
K_5_percent[3,3] = K_5_percent[3,3]*(1+percentage_of_change)
K_5_percent_ext[3] = K_5_percent_ext[3] + external_connections_to_add

microcircuit_5_percent = microcircuit.change_parameters(
    {'K': K_5_percent,
     'K_ext': K_5_percent_ext})

After changing the network connectivity, we calculate the power spectra again.

# calculate working point for exponentially shape post synaptic currents
nnmt.lif.exp.working_point(microcircuit_5_percent, method='taylor')
# calculate the transfer function
nnmt.lif.exp.transfer_function(microcircuit_5_percent, method='taylor')
# calculate the delay distribution matrix
nnmt.network_properties.delay_dist_matrix(microcircuit_5_percent)
# calculate the effective connectivity matrix
nnmt.lif.exp.effective_connectivity(microcircuit_5_percent)
# calculate the power spectra
power_spectra_5_percent = nnmt.lif.exp.power_spectra(microcircuit_5_percent)

Here we change the connectivity by 10 percent, similar to above.

percentage_of_change = 0.10

external_connections_to_add = (
    abs(g) * connections_4I_4I * percentage_of_change * population_rate
    ) / external_input_rate

K_10_percent = microcircuit.network_params['K'].copy()
K_10_percent_ext = microcircuit.network_params['K_ext'].copy()
K_10_percent[3,3] = K_10_percent[3,3]*(1+percentage_of_change)
K_10_percent_ext[3] = K_10_percent_ext[3] + external_connections_to_add

microcircuit_10_percent = microcircuit.change_parameters(
    {'K': K_10_percent,
     'K_ext': K_10_percent_ext})

After changing the network connectivity, we calculate the power spectra again.

# calculate working point for exponentially shape post synaptic currents
nnmt.lif.exp.working_point(microcircuit_10_percent, method='taylor')
# calculate the transfer function
nnmt.lif.exp.transfer_function(microcircuit_10_percent, method='taylor')
# calculate the delay distribution matrix
nnmt.network_properties.delay_dist_matrix(microcircuit_10_percent)
# calculate the effective connectivity matrix
nnmt.lif.exp.effective_connectivity(microcircuit_10_percent)
# calculate the power spectra
power_spectra_10_percent = nnmt.lif.exp.power_spectra(microcircuit_10_percent)

Then we print the critical frequencies for each eigenmode

# Look at the critical frequencies per eigenmode
for k, v in sensitivity_dict.items():
    print(k, v['critical_frequency'])

and we plot the results.

width = 180. / 25.4
height = 60. / 25.4

fig = plt.figure(figsize=(width, height),
                 constrained_layout=True)

labels = ['2/3E', '2/3I', '4E', '4I', '5E', '5I', '6E', '6I']

colormap = mpl.cm.get_cmap('coolwarm').copy()

# set colorbar max and min
z = 1

# by default np.nan are set to black with full transparency
# .eps can't handle transparency
colormap.set_bad('w',1.)

fig = plt.figure(figsize=(width, height),
                 constrained_layout=True)
grid_specification = gridspec.GridSpec(1, 2,
                                       height_ratios=[1],
                                       width_ratios=[2.2, 1], figure=fig)

gs = gridspec.GridSpecFromSubplotSpec(1,3,
                                        height_ratios=[1],
                            width_ratios=[1, 1, 0.2],
                            subplot_spec=grid_specification[0])

ev = '5'
ax = fig.add_subplot(gs[0])
label_prms = dict(x=-0.3, y=1.27, fontsize=10, fontweight='bold',
                  va='top', ha='right')

ax.text(s='(A)', transform=ax.transAxes, **label_prms)

frequency = sensitivity_dict[ev]['critical_frequency']
projection_of_sensitivity_measure = sensitivity_dict[ev][
    'sensitivity_amp']

rounded_frequency = str(int(np.round(frequency,0)))

plot_title = r'$\mathbf{Z}_{b=%s}^{\mathrm{amp}}(' % ev + \
    f'{rounded_frequency}' + r'\,\mathrm{Hz})$'
ax.set_title(plot_title)

data = np.ma.masked_where(projection_of_sensitivity_measure == 0,
                            projection_of_sensitivity_measure)

# np.flipud needed to cope for different origin compared to imshow
heatmap = ax.pcolormesh(np.flipud(data),
                    vmin=-z,
                    vmax=z,
                    cmap=colormap,
                    edgecolors='k',
                    linewidth=0.6)

ax.set_aspect('equal')

if labels is not None:
    ax.set_xticks(np.arange(len(labels))+0.5)
    ax.set_yticks(np.arange(len(labels))+0.5)
    ax.set_xticklabels(labels, rotation=90)
    ax.set_yticklabels(list(reversed(labels)))

ax.set_xlabel('sources')
ax.set_ylabel('targets')


# sensitivity_measure_frequency
ax = fig.add_subplot(gs[1])

frequency = sensitivity_dict[ev]['critical_frequency']
projection_of_sensitivity_measure = sensitivity_dict[ev][
    'sensitivity_freq']

rounded_frequency = str(int(np.round(frequency,0)))

plot_title = r'$\mathbf{Z}_{b=%s}^{\mathrm{freq}}(' % ev + \
    f'{rounded_frequency}' + r'\,\mathrm{Hz})$'
ax.set_title(plot_title)

data = np.ma.masked_where(projection_of_sensitivity_measure == 0,
                            projection_of_sensitivity_measure)

heatmap = ax.pcolormesh(np.flipud(data),
                    vmin=-z,
                    vmax=z,
                    cmap=colormap,
                    edgecolors='k',
                    linewidth=0.6)

ax.set_aspect('equal')

if labels is not None:
    ax.set_xticks(np.arange(len(labels))+0.5)
    ax.set_yticks(np.arange(len(labels))+0.5)
    ax.set_xticklabels(labels, rotation=90)
    ax.set_yticklabels([])

ax.set_xlabel('sources')

colorbar_ax = fig.add_subplot(gs[2])

colorbar_width = 0.1
ip = InsetPosition(ax, [1.05,0,colorbar_width,1])
colorbar_ax.set_axes_locator(ip)
colorbar(heatmap, cax=colorbar_ax)

ax = fig.add_subplot(grid_specification[1])
label_prms = dict(x=-0.2, y=1.2, fontsize=10, fontweight='bold',
                  va='top', ha='right')
ax.text(s='(B)', transform=ax.transAxes, **label_prms)
j = 3
ax.plot(microcircuit.analysis_params['omegas']/(2*np.pi),
        power_spectra[:, j],
        color='#4055C8', zorder=2,
        label='original'
        )

ax.plot(microcircuit_5_percent.analysis_params['omegas']/(2*np.pi),
        power_spectra_5_percent[:, j],
        color='#6072D5', zorder=2, ls='dashed',
        label='5% increase'
        )

ax.plot(microcircuit_10_percent.analysis_params['omegas']/(2*np.pi),
        power_spectra_10_percent[:, j],
        color='#8A98E5', zorder=2, ls='dotted',
        label='10% increase'
        )

ax.set_xlim([0.0, 100.0])
ax.set_ylim([1e-6, 1e0])
ax.set_yscale('log')
ax.set_xticks([20, 40, 60, 80])

ax.set_xlabel(r'frequency $\omega/2\pi\quad(1/\mathrm{s})$')

y_minor = matplotlib.ticker.LogLocator(
    base = 10.0,
    subs = np.arange(1.0, 10.0) * 0.1,
    numticks = 10)
ax.yaxis.set_minor_locator(y_minor)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
ax.set_yticks([])
ax.set_yticks([1e-5,1e-3,1e-1])
ax.set_ylabel(r'power spectrum $P_{\mathrm{4I}}(\omega)\quad(1/\mathrm{s}^2)$')
ax.legend(loc='lower right')

plt.savefig('figures/sensitivity_measure_low_gamma_Bos2016.eps')