Spikes Module

Using the bmtool spike analysis module¶

By Gregory Glickert

In [1]:

from bmtool.analysis.spikes import load_spikes_to_df, get_population_spike_rate, compute_firing_rate_stats
from bmtool.bmplot import raster, plot_firing_rate_pop_stats, plot_firing_rate_distribution
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

from bmtool.analysis.spikes import load_spikes_to_df, get_population_spike_rate, compute_firing_rate_stats from bmtool.bmplot import raster, plot_firing_rate_pop_stats, plot_firing_rate_distribution import matplotlib.pyplot as plt import numpy as np %matplotlib inline

Warning: no DISPLAY environment variable.
--No graphics will be displayed.

load_spikes_df¶

We can use the load_spikes_df to load our network output into a dataframe. It can also be useful to provide the config for the network so we can see which population the cells come from.

In [2]:

config_path = '/home/gjgpb9/cortex_modeling/V1_Layer5/Model-Parameters/simulation_config_baseline.json'
output_path = '/home/gjgpb9/cortex_modeling/V1_Layer5/Run-Storage/final_result_2/block1/long/spikes.h5'

df = load_spikes_to_df(spike_file=output_path,network_name='cortex',config=config_path)
df

config_path = '/home/gjgpb9/cortex_modeling/V1_Layer5/Model-Parameters/simulation_config_baseline.json' output_path = '/home/gjgpb9/cortex_modeling/V1_Layer5/Run-Storage/final_result_2/block1/long/spikes.h5' df = load_spikes_to_df(spike_file=output_path,network_name='cortex',config=config_path) df

Out[2]:

	node_ids	timestamps	pop_name
0	8622	7.1	FSI
1	9477	7.3	LTS
2	9903	7.7	LTS
3	9300	8.4	FSI
4	3791	9.0	CP
...	...	...	...
464167	9864	14499.6	LTS
464168	9681	14499.6	LTS
464169	1890	14499.6	CP
464170	6613	14499.7	CS
464171	1647	14499.9	CP

464172 rows × 3 columns

You can actually have the function label the spikes with any attribute that is there when you build the nodes. By default it is pop_name but you could do anything. So if your build looks like this.

net.add_nodes(N=n_I, pop_name='PV',     # N = number of inhibitory cells
        model_type='biophysical',
        model_template='hoc:WBInhCell',  # WBInhCell hoc definition
        morphology='blank.swc',
        positions = positions_cuboid(N=n_I,center=[0.,0.,0.],xside_length=10.,yside_length=10.,height=10.,min_dist=2))

You could use, N,pop_name,model_type,model_template,morphology,and position to label the cells. For this example N and positions would be poor choices to label due to being unique for every cell.

In [3]:

df = load_spikes_to_df(output_path,network_name='cortex',config=config_path,groupby=['model_template','pop_name'])
df

df = load_spikes_to_df(output_path,network_name='cortex',config=config_path,groupby=['model_template','pop_name']) df

Out[3]:

	node_ids	timestamps	model_template	pop_name
0	8622	7.1	hoc:FSI_Cell	FSI
1	9477	7.3	hoc:LTS_Cell	LTS
2	9903	7.7	hoc:LTS_Cell	LTS
3	9300	8.4	hoc:FSI_Cell	FSI
4	3791	9.0	hoc:CP_Cell	CP
...	...	...	...	...
464167	9864	14499.6	hoc:LTS_Cell	LTS
464168	9681	14499.6	hoc:LTS_Cell	LTS
464169	1890	14499.6	hoc:CP_Cell	CP
464170	6613	14499.7	hoc:CS_Cell	CS
464171	1647	14499.9	hoc:CP_Cell	CP

464172 rows × 4 columns

This can be helpful if you have nodes with the same pop_name, but different layers so you may want more than one label to tell which cells are spiking.

raster plot¶

You may want a way to look at the spikes to see when they are firing. One way to do this is by a raster plot

In [4]:

raster(spikes_df=df)
plt.show()

raster(spikes_df=df) plt.show()

You can also play with some settings in the raster function if you want to view the data better or in a different way.

In [5]:

fig, ax = plt.subplots(1, 1, figsize=(8, 5))

color_map = {
    'CS': 'red',
    'CP': 'red', 
    'FSI': 'blue', 
    'LTS': 'green'
}

raster(spikes_df=df,ax=ax,
       groupby='pop_name',
       color_map=color_map,
       tstart=0,tstop=1000)
plt.show()

fig, ax = plt.subplots(1, 1, figsize=(8, 5)) color_map = { 'CS': 'red', 'CP': 'red', 'FSI': 'blue', 'LTS': 'green' } raster(spikes_df=df,ax=ax, groupby='pop_name', color_map=color_map, tstart=0,tstop=1000) plt.show()

Spike Rate¶

There is also a population spike rate function. This can turn your dataframe into a time series, which can then be futher analyzed.

In [6]:

df = load_spikes_to_df(output_path,network_name='cortex',config=config_path)
spike_rate = get_population_spike_rate(df,config=config_path)

fig, ax = plt.subplots(4, 1, figsize=(8, 5))
ax = np.ravel(ax)
for i, (name, rate) in enumerate(spike_rate.items()):
    ax[i].plot(rate)
    ax[i].set_title(name)
    
plt.tight_layout()
plt.show()

df = load_spikes_to_df(output_path,network_name='cortex',config=config_path) spike_rate = get_population_spike_rate(df,config=config_path) fig, ax = plt.subplots(4, 1, figsize=(8, 5)) ax = np.ravel(ax) for i, (name, rate) in enumerate(spike_rate.items()): ax[i].plot(rate) ax[i].set_title(name) plt.tight_layout() plt.show()

Grabbing first network; specify a network name to ensure correct node population is selected.

Firing rate stats¶

You can easily get some firing rate stats and plot them using the below functions

In [7]:

df = load_spikes_to_df(output_path, network_name="cortex", config=config_path,groupby=['model_template','pop_name'])
pop_stats,indivdual_stats = compute_firing_rate_stats(df, groupby=['model_template','pop_name'])
display(indivdual_stats)
display(pop_stats)

df = load_spikes_to_df(output_path, network_name="cortex", config=config_path,groupby=['model_template','pop_name']) pop_stats,indivdual_stats = compute_firing_rate_stats(df, groupby=['model_template','pop_name']) display(indivdual_stats) display(pop_stats)

	node_ids	spike_count	model_template	pop_name	firing_rate
0	9848	852	hoc:LTS_Cell	LTS	58.787812
1	8541	707	hoc:FSI_Cell	FSI	48.782844
2	9430	617	hoc:LTS_Cell	LTS	42.572864
3	8527	617	hoc:FSI_Cell	FSI	42.572864
4	9214	558	hoc:FSI_Cell	FSI	38.501877
...	...	...	...	...	...
9986	7224	1	hoc:CS_Cell	CS	0.069000
9987	8315	1	hoc:CS_Cell	CS	0.069000
9988	8980	1	hoc:FSI_Cell	FSI	0.069000
9989	9341	1	hoc:FSI_Cell	FSI	0.069000
9990	9218	1	hoc:FSI_Cell	FSI	0.069000

9991 rows × 5 columns

	model_template	pop_name	firing_rate_mean	firing_rate_std
0	hoc:CP_Cell	CP	2.210251	1.516544
1	hoc:CS_Cell	CS	2.173785	1.428758
2	hoc:FSI_Cell	FSI	4.494283	5.867914
3	hoc:LTS_Cell	LTS	14.776037	7.553584

In [8]:

plot_firing_rate_distribution(indivdual_stats,groupby='pop_name',plot_type='violin') # type can be box,swarm,violin or a combo 
plt.show()

plot_firing_rate_distribution(indivdual_stats,groupby='pop_name',plot_type='violin') # type can be box,swarm,violin or a combo plt.show()

In [9]:

plot_firing_rate_pop_stats(pop_stats, groupby="pop_name")
plt.show()

plot_firing_rate_pop_stats(pop_stats, groupby="pop_name") plt.show()

Like the raster function you can also make some changes to the appearance of the plot if you want. Note both plotting functions have the ability only one is shown

In [10]:

df = load_spikes_to_df(output_path, network_name="cortex", config=config_path,groupby='pop_name')
pop_stats,indivdual_stats = compute_firing_rate_stats(df, groupby="pop_name")

fig, ax = plt.subplots(1, 1, figsize=(8, 5))
plot_firing_rate_pop_stats(pop_stats, groupby='pop_name',color_map=color_map,ax=ax)

plt.show()

df = load_spikes_to_df(output_path, network_name="cortex", config=config_path,groupby='pop_name') pop_stats,indivdual_stats = compute_firing_rate_stats(df, groupby="pop_name") fig, ax = plt.subplots(1, 1, figsize=(8, 5)) plot_firing_rate_pop_stats(pop_stats, groupby='pop_name',color_map=color_map,ax=ax) plt.show()

compare_firing_over_times¶

You may want to compare the firing rates during different times during your simulation and compute if the firing is different. This can be useful if you are giving a changing input and want to see if you cells are responding differently.

In [11]:

from bmtool.analysis.spikes import compare_firing_over_times

spike_df = load_spikes_to_df("/home/gjgpb9/cortex_modeling/V1_Layer5/Run-Storage/final_result_2/block1/long/spikes.h5", "cortex",config='/home/gjgpb9/cortex_modeling/V1_Layer5/Model-Parameters/simulation_config_baseline.json')
compare_firing_over_times(spike_df,group_by='pop_name',time_window_1=(0,1000),time_window_2=(1000,2000))

from bmtool.analysis.spikes import compare_firing_over_times spike_df = load_spikes_to_df("/home/gjgpb9/cortex_modeling/V1_Layer5/Run-Storage/final_result_2/block1/long/spikes.h5", "cortex",config='/home/gjgpb9/cortex_modeling/V1_Layer5/Model-Parameters/simulation_config_baseline.json') compare_firing_over_times(spike_df,group_by='pop_name',time_window_1=(0,1000),time_window_2=(1000,2000))

Population: FSI
    Average firing rate in window 1: 4.60 Hz
    Average firing rate in window 2: 4.69 Hz
    U-statistic: 366210.50
    p-value: 0.40101877211121495
    Significant difference (p<0.05): No
Population: LTS
    Average firing rate in window 1: 2.14 Hz
    Average firing rate in window 2: 20.67 Hz
    U-statistic: 22456.00
    p-value: 6.613772768519174e-176
    Significant difference (p<0.05): Yes
Population: CP
    Average firing rate in window 1: 1.57 Hz
    Average firing rate in window 2: 2.68 Hz
    U-statistic: 9692980.50
    p-value: 1.629101314611588e-10
    Significant difference (p<0.05): Yes
Population: CS
    Average firing rate in window 1: 1.57 Hz
    Average firing rate in window 2: 2.44 Hz
    U-statistic: 9945849.00
    p-value: 1.0708612377729946e-18
    Significant difference (p<0.05): Yes