Surveying Neural Spiking in Visual Cortex

Import Libraries

Import the modules required for this notebook

import pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
import pingouin as pg

Download Data

Download the data required for this notebook

import requests
url = "https://uni-bonn.sciebo.de/public.php/webdav/YTgrGyB5MekRP2n"
print("Downloading Data ...")
response = requests.get(("/").join(url.split("/")[:-1]), auth=(url.split("/")[-1], "ibots"))
response.raise_for_status()
with open("flash_spikes.parquet", "wb") as file:
    file.write(response.content)
print("Done!")

Downloading Data ...
Done!

In this notebook, we’ll work with data from a recent study conducted at the Allen Institute for Brain Science. The data we’ll use contains the spikes recorded from over 300 neurons recorded across 6 different visual brain areas. In this section, we are going to read the data into a pandas data frame and explore its contents.

Exercise: Load the data stored in flash_spikes.parquet into a data frame and print the first 5 rows.

df = pd.read_parquet("flash_spikes.parquet")
df.head(5)

Exercise: How many spikes are there in the data frame (Hint: each row contains one spike)

df.shape[0]

457706

Exercise: At what time was the first and last spike recorded?

df.spike_time.min(), df.spike_time.max()

(np.float64(1276.2973389650292), np.float64(1572.793837848813))

Exercise: From which (unique) brain areas were the spikes recorded?

df.brain_area.unique()

array(['LM', 'PM', 'AL', 'AM', 'V1', 'RL'], dtype=object)

Exercise: What is the total number of units recorded?

df.unit_id.nunique()

387

Exercise: What is the number of spikes recorded from every unit?

df.groupby("unit_id").spike_time.count()

unit_id
951015628      84
951015643      75
951015704      31
951015717     598
951015731      30
             ... 
951039462     753
951039560     168
951039854      76
951039863     114
951039870    3403
Name: spike_time, Length: 387, dtype: int64

Exercise: What is the number of units recorded in each brain area (Hint: use df.groupby("brain_area")["unit_id"])?

df.groupby("brain_area").unit_id.nunique()

brain_area
AL    82
AM    77
LM    62
PM    68
RL    13
V1    85
Name: unit_id, dtype: int64

Exercise: What is the smallest and largest number of spikes recorded from any single unit?

print(df.groupby("unit_id").spike_time.count().min())
print(df.groupby("unit_id").spike_time.count().max())

2
15498

Even though the metrics we computed in the previous section are useful for describing and understanding the data, they only provide a very coarse and incomplete picture. In this section, we are going to visualize the spike trains to get a more detailed view of what is going on. For this, we will use rasterplots, one of the most common visualizations of neural spiking activity. We can make rasterplots for single neurons or groups of neurons to show how the firing of neurons evolves over time.

Exercise: Get all spikes recorded from unit 951021564

df[df.unit_id == 951021564]

Exercise: Create a rasterplot of the spikes recorded from unit 951021564 (Hint: use plt.eventplot).

Solution

plt.eventplot(df[df.unit_id == 951021564].spike_time)

Exercise: Create a rasterplot of the spikes recorded from unit 951021564 and reduce the linewidth to 0.5.

Solution

plt.eventplot(df[df.unit_id == 951021564].spike_time, linewidth=0.5)

Exercise: Create a rasterplot for the units 951021564 and 951021835 (Hint: plt.eventplot accepts a lists as input).

Solution

spike_times = [
    df[df.unit_id == 951021564].spike_time,
    df[df.unit_id == 951021835].spike_time,
]
plt.eventplot(spike_times, linewidth=0.5)

Example: Extract the spikes recorded from every unit in the primary visual cortex "V1".

mask = df.brain_area == "V1"
spike_times = df[mask].groupby("unit_id").spike_time.unique().tolist()

Exercise: Create a rasterplot displaying the spike_times from each unit between 1300 and 1350 seconds.

Solution

plt.eventplot(spike_times, linewidth=0.5, linelengths=0.5, lineoffsets=1.5)
plt.xlim(1300, 1350)

(1300.0, 1350.0)

Exercise: Create a rasterplot displaying the spikes from all units recorded in the posteriomedial area "PM" between 1300 and 1350 seconds.

Solution

mask = df.brain_area == "PM"
spike_times = df[mask].groupby("unit_id").spike_time.unique().tolist()
plt.eventplot(spike_times, linewidth=0.5, linelengths=0.5, lineoffsets=1.5)
plt.xlim(1300, 1350)

(1300.0, 1350.0)

Exercise: Create a rasterplot displaying the spikes from all units recorded in the rostrolateral area "RL" between 1400 and 1500 seconds.

Solution

mask = df.brain_area == "RL"
spike_times = df[mask].groupby("unit_id").spike_time.unique().tolist()
plt.eventplot(spike_times, linewidth=0.5, linelengths=0.5, lineoffsets=1.5)
plt.xlim(1400, 1500)

(1400.0, 1500.0)

Raster plots are great for detailed visualization of neural firing but with all this detail, it is difficult to make out larger trends and compare the activity between larger groups of neurons. In this section, we are going to use the seaborn library to create plots that aggregate information across spikes and units which allows us to compare the different areas in the visual cortex.

Exercise: Create a sns.countplot that displays the number of spikes recorded from each "brain_area".

Solution

sns.countplot(df, x="brain_area")

Exercise: Create a sns.countplot that displays the number of spikes recorded from each "unit" (Hint: you can remove the x-ticks with plt.xticks([])).

Solution

sns.countplot(df, x="unit_id")
plt.xticks([])

Exercise: Create a sns.countplot that displays the number of spikes recorded from each "unit" and add a hue to encode "brain_area".

Solution

sns.countplot(data=df, x="unit_id", hue="brain_area")
plt.xticks([])

Exercise: Create the same plot again but with horizontal bars

Solution

sns.countplot(data=df, y="unit_id", hue="brain_area")
plt.yticks([])

Exercise: Create a sns.kdeplot to visualize the distribution of "spike_time"

Solution

sns.kdeplot(df, x="spike_time")

Exercise: Create a sns.kdeplot to visualize the distribution of "spike_time" and add a hue for "brain_area".

Solution

sns.kdeplot(df, x="spike_time", hue="brain_area")

Exercise: Compute the number of units recorded in each brain area and plot the result as a sns.barplot (Hint: the function only requires a single argument).

Solution

grouping = df.groupby("brain_area").unit_id.nunique()
sns.barplot(grouping)

Exercise: Get the number of spikes recorded from each brain area and divide it by the number of units recorded in each area to get the average number of spikes per unit in each area. Then, plot the result as a sns.barplot.

Solution

grouping = (
    df.groupby("brain_area").spike_time.count()
    / df.groupby("brain_area").unit_id.nunique()
)

sns.barplot(grouping)

The visualizations from the previous session seem to suggest differences in firing rates between visual areas. In this section we are going to quantify these differences. We are going to compute the firing rates for every unit and visualize their distribution across the different brain areas. Then, we are going to use the Mann-Whitney-U (MWU) test to check whether the observed differences are statistically significant. The MWU test is a non-parametric statistical test that estimates how likely it is that two sets of measurements (i.e. firing rates in two areas) come from the same distribution. Based on this probability, captured by the p-value, we can assess whether the difference is above what we would expect by chance.

Exercise: Compute the number of spike times for each "brain_area" and "unit_id".

df.groupby(["brain_area", "unit_id"]).spike_time.count()

brain_area  unit_id  
AL          951034946      53
            951034962       7
            951034978     168
            951034993     402
            951035003    5060
                         ... 
V1          951028461       2
            951028472      35
            951028481      97
            951028618      70
            951028751     269
Name: spike_time, Length: 387, dtype: int64

Exercise: Divide the result from @exr-count by the duration of the recording to get the firing rate of every unit (Hint: you can approximate the duration as the difference between the first and last "spike_time"). Assign the result to a new variable fr.

dur = df.spike_time.max() - df.spike_time.min()
fr = df.groupby(["brain_area", "unit_id"]).spike_time.count() / dur

Exercise: Use fr.reset_index() to create a new data frame where each row represents one unit (Hint: the method does not modify fr so you must expliclity overwrite the original value). Then, rename the columns to label the new column with "firing_rate" and print the renamed data frame.

fr = fr.reset_index()
fr.columns = ["brain_area", "unit_id", "firing_rate"]
fr

Example: Create a sns.barplot for the mean "firing_rate" of every "brain_area".

sns.barplot(fr, x="brain_area", y="firing_rate", estimator="mean")

Exercise: Create a sns.barplot for the standard deviation "std" of the "firing_rate" of every "brain_area".

Solution

sns.barplot(fr, x="brain_area", y="firing_rate", estimator="std")

Exercise: Create a sns.kdeplot to visualize the distribution of "firing_rate" and add a hue to encode the "brain_area".

Solution

sns.kdeplot(fr, x="firing_rate", hue="brain_area")

Exercise: Re-create the sns.kdeplot but use common_norm=False to account for the difference in the number of recorded units per area.

Solution

sns.kdeplot(fr, x="firing_rate", hue="brain_area", common_norm=False)

Example: Perform a Mann-Whitney-U-Test to test if there is a significant difference in firing rate between the primary visual cortex "V1" and the anterolateral area "AL"

x = fr[fr.brain_area == "V1"].firing_rate
y = fr[fr.brain_area == "AL"].firing_rate
pg.mwu(x, y)

Exercise: Look at the graph created in @exm-mean that showed the mean firing rates for each brain area. Based on this plot, which differences would you expect to be significant? Test you expectation by performing Mann-Whitney-U-Test on the respective pair(s).

x = fr[fr.brain_area == "PM"].firing_rate
y = fr[fr.brain_area == "AM"].firing_rate
pg.mwu(x,y)

Example: Apply the Bonferroni correction to adjust the p-values for multiple comparisons and compare them to the given threshold alpha.

pvals = [0.01, 0.04]
alpha = 0.05
pg.multicomp(pvals, alpha, method="bonferroni")

(array([ True, False]), array([0.02, 0.08]))

Exercise: Apply the Bonferroni correction to adjust the p-values returned by the statistical tests you performed in the exercises before. Are the results still significant after correction?