```
In [4]:
```# This file is to analyze the spencer data and generate figures by calling
# source codes
import numpy as np
import scipy.io as sio

```
In [ ]:
```def get_response_mat(spiketrain, imgPara, stimType, goodCells, plotRaster):
numNeuron = len(goodCells) # the number of neurons
print numNeuron
# imgPara.stim_time = 32s, imgPara.dt = 0.075103,
# numFramesPerStim is the number of the frames within 32s movie stimulus
numFramesPerStim = int(round(imgPara['stim_time'] / imgPara['dt']))
# print numFramesPerStim
spikeMat = []
## generate the spike timing for all the neurons through all trials
for rep in range(imgPara['stimrep']):
spikesCurrentTrial = np.zeros((numNeuron, numFramesPerStim))
spikesRaster = []
cellI = 0
for i in range(len(goodCells)):
# spikesI: spiking timing of a specific neuron at a specific trial
# print i
spikesI = spiketrain[0,i][0][rep,stimType]
# print spikesI
spikesI = np.round(spikesI[np.nonzero(spikesI<=numFramesPerStim)])
#print spikesI
spikesI = spikesI[np.nonzero(spikesI>0)];
spikesI = spikesI.astype(int)
spikesI = spikesI - 1
# print spikesI
# along the 426 frames, spike timings was labeled
spikesCurrentTrial[cellI,spikesI] = 1
cellI = cellI +1;
spikesRaster.append(spikesI*imgPara['dt'] -1)
# return spikeMat as the spiking time for all neurons
spikeMat.append(spikesCurrentTrial)
# change spikeMat to be numpy array
spikeMat = np.array(spikeMat)
print spikeMat.shape
return spikeMat

```
In [ ]:
```def calculate_sparsity_over_trials(spikeMat, imgPara):
print spikeMat.shape
stimrep = int(imgPara['stimrep'])
sparsity = np.zeros((stimrep,1))
numFramesPerStim = int(round(imgPara['stim_time']/imgPara['dt']))
# the returned sparsity is everaged for all the neurons'
# firings within each specific trial
for rep in range(stimrep):
numFramesArray = np.array(range(numFramesPerStim))
spikesCurrentTrial = spikeMat[rep]
sparsity[rep] = np.mean(spikesCurrentTrial) # average by columns and by rows
return sparsity

```
In [3]:
```# this function code has to be debugged further with matrix dimension matching
from scipy import signal
def calculate_subset_index(spikeMat, numFramesPerStim):
numNeuron = spikeMat.shape[1]
numFrames = spikeMat.shape[2]
numRpts = spikeMat.shape[0]
spikeEarly = np.zeros(numNeuron, numFramesPerStim)
# spikeEarly only count the spikes in the first trial
spikeEarly = spikeMat[0]
timewindow = 1
sharedNeuronsAll = np.zeros(4)
numNeuronsLateAll = np.zeros(4)
tt = 0
for rep in range(17,21): # compute trials from 17 to 20
spikeLate = spikeMat[rep]
spikeLate = signal.convolve2d(spikeLate, np.ones((1, timewindow)),'same')
spikesLate[np.nonzero(spikesLate>0)] = 1
# find all the cells both fire both at the 1st trial and the last trial
sharedNeurons= numpy.multiply(spikeEarly, spikeLate).sum(axis=0)
numNeuronLate = spikeLate.sum(axis=0)
sharedNeuronsAll[tt] = sharedNeurons.sum()
numNeuronsLateAll[tt] = numNeuronsLate.sum()
tt = tt+1
# calculate the subsetindex for trials from 17 to 20, the last 4 trials
subsetIndex = sharedNeuronsAll.sum()/numNeuronsLateAll.sum()
return subsetIndex

```
In [6]:
```# this function code has to be debugged further with matrix dimension matching
def generate_poisson_spikes(spikeMat, imgPara):
meanFiringRate = np.mean(spikeMat, 1)
numNeuron = spikeMat.shape[1]
NT = spikeMat.shape[2]
# NT is number of frames through 20 trials
poissSpikes = np.zeros(numNeuron, NT)
for i in range(numNeuron):
poissSpikes[i,:] = np.random.poisson(meanFiringRate[i], (1, NT))
return poissSpikes

```
In [9]:
```# this function code has to be debugged further with matrix dimension matching
def shuffle_spikes(spikeMat, imgPara):
[numNeuron, NT]= spikeMat.shape
fakeSpikes = np.zeros(numNeuron, NT)
numFramesPerStim = int(round(imgPara['stim_time'] / imgPara['dt']))
for i in range(numNeuron):
spikeR = spikeTrains[i,:].reshape((numFramesPerStim, imgPara['stimrep']))
shuffleSpikes = np.zeros(numFramesPerStim, imgPara['stimrep'])
for t in range(numFramesPerStim):
# shuffle the spikeMat among 20 trials
shuffleSpikes[t, :] = spikeR[t, np.random.permutation(imgPara['stimrep'])]
# generate the fake-shuffled spikes for each specific neuron
fakeSpikes[i, :] = shuffleSpikes.reshape(( 1, NT))
return fakeSpikes

```
In [8]:
```# this function code has to be debugged further with matrix dimension matching
from itertools import combinations
def synchrony_analysis_efficient(spikeMat, numCoactive):
minRpts = 1
#numSteps is the total number of frames through all 20 trials
[numNeurons, numSteps] = spikeMat.shape
popRate = spikeMat.sum(axis=0) # sum all the spike numbers across rows
binaryWords = []
frequency = []
for t in range(numSteps): # numSteps = 20*426 = 8520
if popRate[t] >= numCoactive: # match requirement of coactive number of neurons
activeCells = np.nonzeros(spikeMat[:,t]>0) # return the indices of the neurons firing at time t
# enumerate all permutations
c = combinations(range(activeCells), numCoactive) # return all the possible combinations of fired cell assembly
if numCoactive > 1:
binaryWords.append(c)
# frequency returns the counts of firings of cell assembly together through 20 trials
frequency = np.zeros((len(binaryWords),1)) # of cell assembly by 1
for i in range(len(binaryWords)):
frequency[i] = np.nonzeros(spikeMat[binaryWords[i,:],:].sum()>=numCoactive).sum()
## merge
idx= np.nonzeros(frequency==1)
binaryWordsNew = binaryWords[idx, :]
frequencyNew = frequency[idx]
for numRpts in range(2, max(frequency)):
idx= np.nonzeros(frequency==numRpts)
spikePatterns = binaryWords[idx,:]
uniqueSpikes = np.zeros((len(idx),numCoactive))
k=1 # the next two for loops is to remove duplicates from spikePattern
for i in range(len(spikePatterns)):
presence = 0
for j in range(k):
if max(abs(spikePatterns[i,:]-uniqueSpikes[j,:]))==0:
presence = 1
break
if presence ==0:
uniqueSpikes[k,:] = spikePatterns[i,:]
k = k + 1
uniqueSpikes = uniqueSpikes[1:k-1,:]
binaryWordsNew.append(uniqueSpikes)
frequencyNew.append(np.ones((len(uniqueSpikes),1))*numRpts)
# returned frequency is the firing frequency for each cell assembly
# returned binaryWords is the corresponding cell indices of cell assembly
frequency = frequencyNew
binaryWords = binaryWordsNew
idx = np.nonzeros(frequency>=minRpts)
frequency = frequency[idx]
binaryWords = binaryWords[idx,:]
return frequency, binaryWords

```
In [10]:
```# this function code has to be debugged further with matrix dimension matching
def calculate_occurance_time(frequency, numFramesPerStim, spikeMat, binaryWords, numCoactive):
# check occurence of high-order patterns
idx = np.nonzeros(frequency >=2) # the frequency of spikings
timeInMovieSTD = np.zeros((idx.shape[0],1))
binEdges = range(-20,21)
timeJitterDistribution = np.zeros((1, len(binEdges)))
trialDist = []
interval = []
timeAll = []
timeInMovieDist = np.zeros(numFramesPerStim, 1)
occurance_time = np.zeros(len(frequency), 1)
for i in range(len(idx)): # search through the neurons that fired at least twice
# find out the time points when there are at least numCoactive neurons firing together
times = np.nonzeros(spikeMat[binaryWords[i,:],:].sum()>=numCoactive)
# kl: find out which the exact time frame the cell assembly respond to
timeInMovie = times % numFramesPerStim
# kl: the response at 0 is assumed to be responding to the last movie frame
timeInMovie[timeInMovie==0] = numFramesPerStim
# distribution of the time points when coactive firing
timeInMovieDist[timeInMovie] = timeInMovieDist[timeInMovie]+1
# what trials includes the coactive firings
trial = (times-timeInMovie)/numFramesPerStim + 1
# average the time points of coactive firings for all the frequency
occurance_time[idx[i]] = np.mean(timeInMovie)
# histogram of the coactive firings around average time points within binEdges
counts = np.histogram(timeInMovie-np.mean(timeInMovie), range=binEdges)
# normalize counts
if counts.sum()>1 :
counts = counts/counts.sum()
# the histogram of distribution of firing time points of coactive assembly around average time points
timeJitterDistribution = timeJitterDistribution + counts
timeInMovieSTD[i] = np.std(timeInMovie)
# the histogram of distribution of firing time points of coactive assembly around average time points
timeJitterDistribution = timeJitterDistribution/len(idx)

```
In [ ]:
```# this is the general code for draw pictures for spencer's data
# Combo3_V1andAL.mat includes all the 5 recordings of calcium imaging data
# within V1 and AL recorded simultaneously.
data = sio.loadmat("./data/2016-10-26_1/Combo3_V1.mat")
data_c = data['data']
# analyze the 3rd data, and spiketrain recorded from V1
data_c = data['data']
# imgPara
imgPara = data_c[0][0][1]
print imgPara
stim_type = imgPara['stim_type']
stim_time = imgPara['stim_time']
updatefr = imgPara['updatefr']
intertime = imgPara['intertime']
stimrep = imgPara['stimrep']
dt = imgPara['dt']
F = imgPara['F']
F_gray = imgPara['F_gray']
# spiketrain
spiketrain = data_c[0][0][2]
# number of neurons
numNeuron = spiketrain.shape[1]
numFramesPerStim = int(round(stim_time / dt))
spikesPerNeuron = np.zeros((numNeuron,3))
for i in range(numNeuron):
for j in range(stim_type):
numSpike = 0
for rep in range(stimrep):
spikesI = spiketrain[0,i][0][rep,j]
numSpike = numSpike + len(spikesI[0]);
spikesPerNeuron[i, j] = numSpike;
# print spikesPerNeuron
print "Number of cells: %d" % numNeuron
# Population Response to Natural Stimuli
# goodCells = range(numNeuron); # choose all the cells to be goodCells
goodCells = np.nonzero(spikesPerNeuron[:,stimType]>3)
print goodCells[0].shape
spikeMat = get_response_mat(spiketrain, imgPara, stimType,goodCells[0], 0);
# show sparsity over trials
sparsity = calculate_sparsity_over_trials(spikeMat, imgPara)
# sparsity plot figure(1)
plt.plot(sparsity,'k')
plt.xlabel('trial #')
plt.ylabel('sparseness')
plt.show()
# subset analysis, an array of numbers trial 17 to 20
subsetIndex = calculate_subset_index(spikeMat, numFramesPerStim);
# generate Poisson spike trains
subsetIndexShuffleList = []
subsetIndexPoissonList=[]
for rep in range(20):
# generate Poisson spike trains
poissSpikes = generate_poisson_spikes(spikeMat, imgPara)
# generate shuffled spike trains
shuffledSpikes = shuffle_spikes(spikeMat, imgPara);
subsetIndexShuffle = calculate_subset_index(shuffledSpikes, numFramesPerStim)
subsetIndexPoisson = calculate_subset_index(poissSpikes, numFramesPerStim)
subsetIndexPoissonList.append(subsetIndexPoisson)
subsetIndexShuffleList.append(subsetIndexShuffle)
# from here, begin to analyze the cell assembly result
numCoactive = 3 # 3 cell assembly
[frequency, binaryWords] = synchrony_analysis_efficient(spikeMat, numCoactive)
[binEdges, timeJitterDistribution3, timeInMovieDist3] = calculate_occurance_time(...
frequency, numFramesPerStim, spikeMat, binaryWords, numCoactive);