Single Unit Data and Spike Trains¶

Based on assignment 4
This Demo and the attached files from course PSYO/NESC 3505 Neural Data Science, at Dalhousie University. You should not disseminate, distribute or copy.
I am NOT post inappropriate or copyrighted content, advertise or promote outside products or organizations
© Copyright 2020.PSYO/NESC 3505 FAll 2020 https://dalpsychneuro.github.io/NESC_3505_textbook/intro.html
For demonstration purposes only

Background¶

Data Source¶

These data were provided courtesy of Dr. Nathan Crowder, Department of Psychology & Neuroscience, Dalhousie University. Data are from 23 mouse primary visual cortex neurons; for simplicity we are going to work with data firstly from only one neuron, then we’ll load data from 3 neurons to get experience working with a larger data set. Note that these are not data from a multielectrode array. Data were recorded from single electrodes that were embedded directly into single neurons, a technique called patch clamping. We have data from multiple neurons because Dr. Crowder’s team painstakingly recorded from several neurons in each of a few animals.

—–PSYO 3505 Assignment 4 cell1

Experimental Design¶

The experiment involved presenting drifting (oscillating) sine wave gratings to mice while recording from electrodes in the primary visual cortex. The gratings oscillated (drifted) at a speed of 2 Hz. The contrast of the stimuli were varied systematically over ten levels from 4 (virtually no contrast between the brightest and darkest parts of the sine wave grating, which would appear as uniform grey) to 100% (gratings varying between black and white, as in the example grating in the textbook). Contrast is therefore one of the experimental variables; the levels of contrast used were 4, 8, 12, 16, 24, 32, 48, 64, 84, and 100%.

A second variable is condition:the stimuli were presented under 2 experimental conditions: control (CTRL) and adaptation (ADAPT):

In the CTRL condition, each trial began with a 2000 ms blank grey screen, then the oscillating sine wave grating at the trial’s contrast level (see below) for 1000 ms, then 1000 ms blank grey screen. In the ADAPT condition, each trial started with the adaptation stimulus (a 50% contrast sine wave grating for 2000 ms), then the sine wave grating at the trial’s contrast level (see below), then 1000 ms blank grey screen. A total of 8 trials were presented at each contrast level, in each condition. Each trial was 4000 ms in duration, and data were sampled at 1000 Hz, meaning that we have a data point every 1 ms.

—–PSYO 3505 Assignment 4 cell2

Predictions¶

In the control condition, latency to first spike is expected to decrease with increasing contrast.

Some primary visual cortex cells are phase sensitive (so-called simple cells), so in response to a drifting sine wave grating they will oscillate between excitation and inhibition at the same temporal frequency of the stimulus (2 Hz in this case).

—–PSYO 3505 Assignment 4 cell3

Single neuron data¶

Read file and check it¶

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

First of all we need define the experiment parameters

# Conditions
cond_labels = ['CTRL', 'ADAPT'] 

# Labels for the levels of stimulus contrast
contr_labels = [4, 8, 12, 16, 24, 32, 48, 64, 84, 100]

# Labels for repetitions (trials)
rep_labels = list(np.arange(1, 8))
num_reps = len(rep_labels)

# Time point labels - each trial is 4000 ms long
time_labels = list(np.arange(4000))

# Stimulus timing
stim_on_time = 2000
stim_off_time = stim_on_time + 1000

adapt_on_time = 0
adapt_off_time = adapt_on_time + 2000

df = pd.read_csv('crowder_1_neuron.csv')

df.head()

	neuron	time	repetition	contrast	condition
0	m6_11	0	1	4	CTRL
1	m6_11	1	1	4	CTRL
2	m6_11	2	1	4	CTRL
3	m6_11	3	1	4	CTRL
4	m6_11	4	1	4	CTRL

df.tail()

	neuron	time	repetition	contrast	condition
639995	m6_11	3995	8	100	ADAPT
639996	m6_11	3996	8	100	ADAPT
639997	m6_11	3997	8	100	ADAPT
639998	m6_11	3998	8	100	ADAPT
639999	m6_11	3999	8	100	ADAPT

By using .unique() methods We can get idea about number of repertition and contrasts

print(df['repetition'].unique())
print(df['contrast'].unique())

[1 2 3 4 5 6 7 8]
[  4   8  12  16  24  32  48  64  84 100]

We want to know number of spikes in DF

len(df[df['spike'] == 1])

Because the data is stored in the data frames, ‘Advanced indexing’ technique can apply here. like if I want to know the repetion 1 Contrast 100 CTRL and spike data I can write code like this :

df[(df['repetition']==1)&(df['contrast']==100)&(df['condition']=='CTRL')&(df['spike']==1)]

	neuron	time	repetition	contrast	condition	spike
288797	m6_11	797	1	100	CTRL	1
290153	m6_11	2153	1	100	CTRL	1
290362	m6_11	2362	1	100	CTRL	1
290420	m6_11	2420	1	100	CTRL	1
290439	m6_11	2439	1	100	CTRL	1
290473	m6_11	2473	1	100	CTRL	1
290500	m6_11	2500	1	100	CTRL	1
290557	m6_11	2557	1	100	CTRL	1
290633	m6_11	2633	1	100	CTRL	1
290657	m6_11	2657	1	100	CTRL	1
290715	m6_11	2715	1	100	CTRL	1
290759	m6_11	2759	1	100	CTRL	1
290897	m6_11	2897	1	100	CTRL	1
291050	m6_11	3050	1	100	CTRL	1

Visualization¶

This is a complex experiment so we’ll need lots of subplots to see all trials for all contrast levels, and both conditions. Here we can using Raster plots show the experitement

fig = plt.figure(figsize=[12,18])

subplot_counter = 1 

for contr in contr_labels:
    for cond in cond_labels:

        ax = fig.add_subplot(len(contr_labels), len(cond_labels), subplot_counter)

        # Show adaptation grating at 50% contrast
        if cond == 'ADAPT':
            plt.axvspan(0, stim_on_time-1, alpha=0.5, color='grey')

        # shading indicates stimulus on period
        plt.axvspan(stim_on_time, stim_off_time, alpha=contr/100, color='grey')

        for rep in rep_labels:
            spike_times = df[(df['repetition']==rep)&(df['contrast']==contr)&(df['condition']==cond)&(df['spike']==1)]['time']
            plt.vlines(spike_times, rep-1, rep);

        # Pretty formatting
        plt.xlim([0, max(time_labels)+1])
        plt.ylim([0, len(rep_labels)])
        plt.title(cond + ' ' + str(contr) + '% contrast')
        plt.xlabel('Time (ms)')
        plt.ylabel('Repetition',)
        plt.yticks([x + 0.5 for x in range(num_reps)], [str(x + 1) for x in range(num_reps)], size=8) 
        plt.tight_layout()
        subplot_counter += 1
        
plt.show()

Also, we can draw peri-stimulus time histogram (PTSH) plot rather than rasters. By summing the number of spikes in short time bins, PSTHs allow us to visualize when the neuron is most likely to spike.

hist_bin_width = 50 
time_bins = np.arange(0, max(time_labels), hist_bin_width)

fig = plt.figure(figsize=[12,18])

subplot_counter = 1 

for contr in contr_labels:
    for cond in cond_labels:

        ax = fig.add_subplot(len(contr_labels), len(cond_labels), subplot_counter)

        # Show adaptation grating at 50% contrast
        if cond == 'ADAPT':
            plt.axvspan(0, stim_on_time-1, alpha=0.5, color='grey')
            
        # shading indicates stimulus on period
        plt.axvspan(stim_on_time, stim_off_time, alpha=contr/100, color='grey')
        spike_times = df[(df['contrast']==contr)&(df['condition']==cond)&(df['spike']==1)]['time']
        spike_count, bins = np.histogram(spike_times, bins=time_bins)
        spike_count = spike_count/(num_reps * hist_bin_width)
        plt.bar(bins[:-1], spike_count,width=hist_bin_width)
        
        # Pretty formatting
        plt.xlim([0, max(time_labels)+1])
        plt.ylim([0, 0.05])
        plt.title(cond + ' ' + str(contr) + '% contrast')
        plt.xlabel('Time (ms)')
        plt.ylabel('Spike Probability')
        
        plt.tight_layout()
        subplot_counter += 1
        
plt.show()

Next we try to plot heat map

Like a PSTH, the heat map has time on the x axis, and plots the probability of spikes at each time point. However, the heat maps a more condensed, information-rich, and efficient version of the PSTHs: rather than needing 10 plots to show the histograms for all 10 levels of intensity.

—Neural Data Science Heat Maps, Interpolation, and Colour Map Choice

Using nested list comprehension to create two columns: condition and contrast

condition_labels = pd.DataFrame([[cond, contr] for cond in cond_labels for contr in contr_labels ],columns=['condition', 'contrast'])

condition_labels.head()

	condition	contrast
0	CTRL	4
1	CTRL	8
2	CTRL	12
3	CTRL	16
4	CTRL	24

psth_df_list=[]
for cond in cond_labels:
    for contr in contr_labels:
        spike_times = df[(df['contrast']==contr)&(df['condition']==cond)&(df['spike']==1)]['time']
        spike_count, bins = np.histogram(spike_times, bins=time_bins)
        psth_df_list.append(spike_count) 

# Combine labels and data
psth_df = pd.concat([condition_labels, pd.DataFrame(psth_df_list)], axis=1)     

psth_df = psth_df.set_index(['condition','contrast'])

Before the final plot, we need check the dataframe

psth_df.head()

		0	1	2	3	4	5	6	7	8	9	...	69	70	71	72	73	74	75	76	77	78
condition	contrast
CTRL	4	0	0	0	0	0	0	0	0	0	0	...	0	0	0	0	0	0	0	0	0	0
	8	0	0	0	0	1	0	1	0	0	0	...	0	0	0	0	0	0	0	2	0	0
	12	0	1	0	2	0	0	0	0	1	0	...	0	0	0	1	0	0	0	0	0	0
	16	0	0	0	1	0	2	0	0	0	0	...	0	0	0	0	0	1	0	0	0	0
	24	0	0	1	0	1	0	0	0	0	0	...	0	0	0	0	0	1	0	0	0	0

5 rows × 79 columns

fig = plt.figure(figsize=[15,30])
subplot_counter = 1

for cond in cond_labels:
    
    ax = fig.add_subplot(1, 2, subplot_counter)
    plt.imshow(psth_df.loc[cond],origin='lower',cmap='hot',interpolation='bilinear',aspect =5)
    
    plt.title(cond)
    xticks = range(0, len(time_bins), 10)
    plt.xticks(xticks,
               time_bins[xticks])
    plt.xlabel('Time (ms)')
    plt.yticks([x for x,y in enumerate(contr_labels)], 
               [x for x in contr_labels])
    if cond == cond_labels[0]:
        plt.ylabel('Stimulus Contrast')

    subplot_counter += 1
    
plt.show()

Additionally, the colour map parameter can be changed. Consider the colour-blindness and have good perceptually, like ‘viridis’below

fig = plt.figure(figsize=[15,30])
subplot_counter = 1 # used to track subplots

for cond in cond_labels:
    
    ax = fig.add_subplot(1, 2, subplot_counter)
    plt.imshow(psth_df.loc[cond],origin='lower',cmap='viridis',interpolation='bilinear',aspect =5)
    
    plt.title(cond)
    xticks = range(0, len(time_bins), 10)
    plt.xticks(xticks,
               time_bins[xticks])
    plt.xlabel('Time (ms)')
    plt.yticks([x for x,y in enumerate(contr_labels)], 
               [x for x in contr_labels])
    if cond == cond_labels[0]:
        plt.ylabel('Stimulus Contrast')

    subplot_counter += 1
    
plt.show()

Working with multiple neurons¶

Read file and check it¶

df3 = pd.read_csv('crowder_3_neurons.csv')

df3.head()

	neuron	time	repetition	contrast	condition
0	m1_6	0	1	4	CTRL
1	m1_6	1	1	4	CTRL
2	m1_6	2	1	4	CTRL
3	m1_6	3	1	4	CTRL
4	m1_6	4	1	4	CTRL

df3.tail()

	neuron	time	repetition	contrast	condition
1919995	m6_3a2	3995	8	100	ADAPT
1919996	m6_3a2	3996	8	100	ADAPT
1919997	m6_3a2	3997	8	100	ADAPT
1919998	m6_3a2	3998	8	100	ADAPT
1919999	m6_3a2	3999	8	100	ADAPT

In this case, we want to know how many neurons in the data

neuron_labels = df3['neuron'].unique()
print(neuron_labels)

['m1_6' 'm3_11' 'm6_3a2']

Visualization¶

Again, we create labels dataframe first

condition_labels = pd.DataFrame([[neur,cond, contr] for neur in neuron_labels for cond in cond_labels for contr in contr_labels ],columns=['neuron','condition','contrast',])
condition_labels.head(10)

	neuron	condition	contrast
0	m1_6	CTRL	4
1	m1_6	CTRL	8
2	m1_6	CTRL	12
3	m1_6	CTRL	16
4	m1_6	CTRL	24
5	m1_6	CTRL	32
6	m1_6	CTRL	48
7	m1_6	CTRL	64
8	m1_6	CTRL	84
9	m1_6	CTRL	100

Due to we have three neurons and each neruons has two conditions and 10 contrasts, we need douvle check it and make sure no error in here

condition_labels

	neuron	condition	contrast
0	m1_6	CTRL	4
1	m1_6	CTRL	8
2	m1_6	CTRL	12
3	m1_6	CTRL	16
4	m1_6	CTRL	24
5	m1_6	CTRL	32
6	m1_6	CTRL	48
7	m1_6	CTRL	64
8	m1_6	CTRL	84
9	m1_6	CTRL	100
10	m1_6	ADAPT	4
11	m1_6	ADAPT	8
12	m1_6	ADAPT	12
13	m1_6	ADAPT	16
14	m1_6	ADAPT	24
15	m1_6	ADAPT	32
16	m1_6	ADAPT	48
17	m1_6	ADAPT	64
18	m1_6	ADAPT	84
19	m1_6	ADAPT	100
20	m3_11	CTRL	4
21	m3_11	CTRL	8
22	m3_11	CTRL	12
23	m3_11	CTRL	16
24	m3_11	CTRL	24
25	m3_11	CTRL	32
26	m3_11	CTRL	48
27	m3_11	CTRL	64
28	m3_11	CTRL	84
29	m3_11	CTRL	100
30	m3_11	ADAPT	4
31	m3_11	ADAPT	8
32	m3_11	ADAPT	12
33	m3_11	ADAPT	16
34	m3_11	ADAPT	24
35	m3_11	ADAPT	32
36	m3_11	ADAPT	48
37	m3_11	ADAPT	64
38	m3_11	ADAPT	84
39	m3_11	ADAPT	100
40	m6_3a2	CTRL	4
41	m6_3a2	CTRL	8
42	m6_3a2	CTRL	12
43	m6_3a2	CTRL	16
44	m6_3a2	CTRL	24
45	m6_3a2	CTRL	32
46	m6_3a2	CTRL	48
47	m6_3a2	CTRL	64
48	m6_3a2	CTRL	84
49	m6_3a2	CTRL	100
50	m6_3a2	ADAPT	4
51	m6_3a2	ADAPT	8
52	m6_3a2	ADAPT	12
53	m6_3a2	ADAPT	16
54	m6_3a2	ADAPT	24
55	m6_3a2	ADAPT	32
56	m6_3a2	ADAPT	48
57	m6_3a2	ADAPT	64
58	m6_3a2	ADAPT	84
59	m6_3a2	ADAPT	100

This time we no longer plot Raster or PTSH (actually we have 3 neurons * 2 conditions * 10 contrasts need 60 histograms !)instead Heat maps directly

Reformat the data from raw data set

psth_df3_list = []
for neuron in neuron_labels:
     for cond in cond_labels:
        for contr in contr_labels:
            spike_times = df3[(df3['neuron']==neuron)&(df3['contrast']==contr)&(df3['condition']==cond)&(df3['spike']==1)]['time']
            spike_count, bins = np.histogram(spike_times,bins=time_bins)
            psth_df3_list.append(spike_count)

#Combine labels and data together 
psth_df3 = pd.concat([condition_labels,pd.DataFrame(psth_df3_list)],axis =1)
psth_df3 = psth_df3.set_index(['neuron','condition','contrast'])
psth_df3.head(20)

			0	1	2	3	4	5	6	7	8	9	...	69	70	71	72	73	74	75	76	77	78
neuron	condition	contrast
m1_6	CTRL	4	0	0	3	2	0	1	4	0	1	4	...	0	2	4	1	1	2	2	4	2	0
		8	1	0	0	1	3	5	2	1	5	2	...	1	3	1	0	1	0	1	1	2	2
		12	1	2	2	0	3	2	4	3	5	1	...	2	4	1	2	3	2	5	0	4	4
		16	5	1	4	5	3	1	2	2	2	1	...	3	3	1	3	2	2	0	5	0	2
		24	1	1	2	1	1	4	2	12	5	1	...	2	0	1	3	1	0	1	0	1	4
		32	1	3	3	2	2	1	1	0	1	0	...	1	3	0	0	0	1	1	1	0	1
		48	2	5	5	0	2	1	3	2	1	2	...	2	1	1	4	0	0	1	0	1	0
		64	0	3	4	3	2	1	0	0	4	0	...	1	2	6	0	3	2	1	0	1	0
		84	1	5	1	4	2	2	2	3	3	3	...	1	0	4	3	1	3	1	3	1	1
		100	2	2	2	3	3	2	1	2	0	1	...	3	0	0	0	0	2	1	0	0	1
	ADAPT	4	4	0	13	4	17	10	14	16	10	9	...	0	4	0	0	3	2	2	0	1	0
		8	0	5	16	10	14	8	12	11	9	12	...	4	4	3	2	3	6	4	4	2	4
		12	2	1	9	9	13	5	10	7	8	7	...	1	1	0	3	0	2	0	0	2	3
		16	3	6	14	11	16	11	12	9	10	11	...	2	1	0	0	0	4	0	4	2	2
		24	3	6	13	10	17	12	9	13	10	6	...	0	2	3	1	1	0	3	1	0	0
		32	1	2	12	7	11	11	9	11	8	7	...	2	1	2	1	0	3	0	1	2	1
		48	1	4	12	7	15	14	11	14	11	9	...	1	2	0	0	2	1	1	2	0	1
		64	1	1	14	7	14	13	10	9	8	8	...	0	3	1	1	0	2	2	3	0	0
		84	4	2	18	10	14	8	6	10	11	7	...	0	0	0	0	0	1	0	0	1	0
		100	1	5	11	9	9	8	11	9	11	8	...	0	0	0	0	1	0	1	3	0	1

20 rows × 79 columns

psth_df3.tail(20)

			0	1	2	3	4	5	6	7	8	9	...	69	70	71	72	73	74	75	76	77	78
neuron	condition	contrast
m6_3a2	CTRL	4	4	4	0	0	1	3	1	1	3	2	...	2	1	1	4	2	7	6	3	2	4
		8	5	1	2	2	1	5	3	1	2	10	...	4	7	3	3	0	0	3	4	3	1
		12	4	1	3	4	2	2	2	3	1	3	...	1	3	2	6	5	5	3	2	1	1
		16	3	4	1	0	1	1	3	4	8	0	...	2	1	0	2	2	3	3	3	5	3
		24	3	4	0	5	2	3	1	6	3	4	...	3	0	2	3	9	1	2	3	4	2
		32	2	2	1	0	3	0	1	0	4	4	...	2	3	1	5	5	1	3	3	2	1
		48	2	2	4	0	7	3	2	2	2	2	...	1	3	3	0	6	5	4	6	0	4
		64	4	4	6	3	4	0	7	3	2	5	...	2	3	0	3	2	0	0	0	1	1
		84	4	0	2	4	6	4	7	2	0	1	...	3	4	4	2	0	1	2	0	1	1
		100	7	4	0	3	4	2	9	4	2	4	...	1	0	0	0	0	0	1	2	1	5
	ADAPT	4	0	0	11	11	13	15	21	12	12	12	...	3	0	0	5	3	1	4	0	1	1
		8	3	3	18	12	12	14	18	13	7	10	...	4	1	2	6	9	3	2	3	2	3
		12	7	5	16	5	12	21	18	15	12	9	...	2	2	0	0	0	0	0	0	0	0
		16	4	5	17	6	21	14	13	19	10	13	...	0	0	3	0	4	1	0	1	5	2
		24	6	7	13	7	10	21	10	19	9	8	...	3	3	3	1	5	2	1	0	3	1
		32	1	0	10	15	10	13	19	12	12	6	...	2	0	3	1	3	2	1	2	2	1
		48	2	2	18	9	15	12	21	22	12	8	...	5	2	2	4	3	4	2	2	3	0
		64	2	6	8	10	19	16	18	15	9	18	...	0	2	0	2	6	0	1	2	2	4
		84	1	1	11	7	9	16	12	11	12	7	...	4	4	2	1	0	4	1	0	3	3
		100	1	2	10	11	16	20	19	11	15	15	...	1	0	0	3	2	2	3	3	1	0

20 rows × 79 columns

fig = plt.figure(figsize=[15,15])
subplot_counter = 1 

for neuron in neuron_labels:
      for cond in cond_labels:
        ax = fig.add_subplot(3, 2, subplot_counter)

        data = psth_df3.loc[neuron]
        
        plt.imshow(data.loc[cond],origin='lower',cmap='viridis',interpolation='bilinear',aspect = 5)

        
        # Nice formatting
        plt.title('Neuron ' + neuron + ': ' + cond)

        xticks = range(0, len(time_bins), 10)
        plt.xticks(xticks,time_bins[xticks])
        
        if neuron == neuron_labels[-1]:
            plt.xlabel('Time (ms)')

        plt.yticks([x for x,y in enumerate(contr_labels)], 
                   [x for x in contr_labels])

        if cond == cond_labels[0]:
            plt.ylabel('Stimulus Contrast')
        subplot_counter += 1

plt.show()

Interpretation¶

According to above heat maps, there are relationship between spike rate (increase) and stimulus onset (regardless of contrast level). Specifically, in m1-6 neuron CTRL condition the ‘hot spot’ is between 2000-2500 (approximately) ms. In m3-11 neuron CTRL condition the ‘hot spot’ is between 2000-2600 (approximately) ms. In m6-3a2 neuron CTRL condition the ‘hot spot is between 2000-3000 (approximately) ms. And this experiment stimulus onset was at 2000ms.

The PTSH graph suggests that 2100 ms at 24% contrast, 2300 at 32%, 2100 at 64%, 2050 at 84%, and 2100 at 100% in the control condition, these number not support the prediction:

in the control condition, latency to first spike is expected to decrease with increasing contrast

It is remarkable that In the adaptation condition the neurons do not show a similar or systematic response to the onset of the “main” stimulus (the one that starts at 2000 ms) as in the control condition. The reason is neuronal ‘habituation (or desensitization)’ probably, namely, The neuron will not stimulus again under similar contrasts.

Additionally, With the contrast increase i.e., from 32% to 100%. the spike rate ‘density’ has significant increase. For example, like above heat maps, neuron m3-11 CTRL, from 48% to 100% the ‘hot spot’ area is larger than other stimulus contrast. Same thing we can found in other two neurons.

Neuron m6_3a2 behavior like simple cells, (response to a drifting sine wave grating they will oscillate between excitation and inhibition at the same temporal frequency of the stimulus), We can observe this phenomenon through heatmaps Neuron m6_3a2: ADAPT . About 100ms excitation then inhibition after 500ms, next 800ms excitation 1000ms inhibition, 1300ms excitation 1400 ms inhibition and so on. Therefore these number support the prediction:

Some primary visual cortex cells are phase sensitive (so-called simple cells), so in response to a drifting sine wave grating they will oscillate between excitation and inhibition at the same temporal frequency of the stimulus (2 Hz in this case,2 Hz means two peaks per second).

References¶

[1] NESC 3505 Neural Data Science, at Dalhousie University. Textbook
[2] NESC 3505 Neural Data Science, at Dalhousie University. Assignment 5

Based on assignment 4
This Demo and the attached files from course PSYO/NESC 3505 Neural Data Science, at Dalhousie University. You should not disseminate, distribute or copy.
I am NOT post inappropriate or copyrighted content, advertise or promote outside products or organizations
© Copyright 2020.PSYO/NESC 3505 FAll 2020 https://dalpsychneuro.github.io/NESC_3505_textbook/intro.html
For demonstration purposes only

Art of Thinking

Single Unit Data and Spike Trains¶

Background¶

Data Source¶

Experimental Design¶

Predictions¶

Single neuron data¶

Read file and check it¶

Visualization¶

Working with multiple neurons¶

Read file and check it¶

Visualization¶

Interpretation¶

References¶