Download
Download this notebook: raw_history_feature.ipynb!
Fit GLMs For Neural Coupling#
GLMs can be used to capture pairwise interactions (couplings) between simultaneously recorded neurons. Let’s see how to do this in NeMoS.
Learn More
You can learn more about GLMs for pairwise interactions by following our head direction tutorial.
Setting up a Fully Coupled GLM#
Raw Spike History as a Feature#
We will start by learning how to model the pairwise interaction using the raw spike history. Let’s first generate some population counts.
import numpy as np
import nemos as nmo
import matplotlib.pyplot as plt
np.random.seed(42)
n_samples = 100
n_neurons = 4
# generate population counts
counts = np.random.poisson(size=(n_samples, n_neurons))
Let’s use the HistoryConv basis to construct a feature matrix capturing the effect of 10 samples of spiking history.
basis = nmo.basis.HistoryConv(window_size=10, label="spike-history")
X = basis.compute_features(counts)
print("Shape of the feature matrix: ", X.shape)
# visualize the output (skip the first window_size of nans)
plt.title("Features", fontsize=16)
plt.imshow(X[10:], cmap="Greys")
plt.xticks([5, 15, 25, 35], [0, 1, 2, 3])
plt.xlabel("Neurons", fontsize=12)
plt.ylabel("Samples", fontsize=12)
plt.show()
Shape of the feature matrix: (100, 40)
This is all you need to fit a fully coupled GLM.
model = nmo.glm.PopulationGLM().fit(X, counts)
model.coef_
Array([[ 0.14045489, -0.18120752, -0.23780271, 0.1358288 ],
[-0.30346444, -0.18104576, -0.06980755, 0.0662384 ],
[-0.34999514, 0.0483201 , 0.28963235, 0.14443153],
[-0.03209282, 0.17395245, 0.16040781, -0.02493234],
[-0.27149802, -0.01220886, -0.11103661, -0.05205673],
[-0.024914 , -0.14321466, 0.15962641, -0.09260146],
[-0.09070015, -0.01457935, -0.11154344, -0.06969942],
[-0.01675784, 0.02318468, -0.07783294, -0.3467661 ],
[-0.03421415, 0.22455376, 0.10233752, -0.09798998],
[-0.29132608, 0.23937738, -0.20096514, -0.01809505],
[ 0.04558323, -0.06524562, -0.22478367, 0.00161728],
[ 0.11102811, -0.04553459, -0.09167013, -0.07718034],
[-0.10426072, 0.19783951, 0.16141443, -0.12835504],
[-0.24666363, 0.06206059, -0.17483383, 0.0792407 ],
[-0.15952398, 0.17811559, 0.00090788, 0.2149349 ],
[ 0.16248585, -0.10513984, -0.05813771, 0.04830596],
[-0.37141153, -0.062241 , -0.0402254 , -0.11875752],
[ 0.16814691, -0.23893166, -0.2647168 , 0.09999399],
[-0.22498278, 0.05915617, -0.10762664, -0.15632086],
[-0.03732371, -0.17941539, 0.10566666, -0.10561956],
[ 0.09111549, 0.14330034, -0.17656595, -0.01338421],
[ 0.14174922, 0.18194321, 0.2282579 , -0.02694508],
[-0.11461368, 0.09806205, -0.00868411, 0.11179899],
[-0.05810085, -0.10380704, -0.3964075 , 0.09961612],
[ 0.3342301 , -0.39886993, -0.01799975, 0.01241543],
[ 0.262705 , -0.23467347, 0.18236944, 0.03809161],
[ 0.08066791, 0.07005222, 0.01067582, 0.11961341],
[ 0.00258666, 0.08246824, -0.23653123, 0.17044774],
[-0.35167852, 0.08866626, 0.08480833, -0.1845658 ],
[ 0.4113944 , -0.16316186, 0.19129303, 0.19340724],
[ 0.14134929, -0.02621229, -0.11538692, -0.18995126],
[ 0.13749726, -0.08522836, -0.01878163, -0.0390583 ],
[ 0.0498317 , 0.06143894, 0.02854554, 0.03721639],
[ 0.07477795, 0.1670006 , -0.00234499, 0.18589541],
[-0.15235828, 0.24513938, 0.28770435, 0.10981987],
[-0.02563755, 0.00379717, -0.27599865, 0.30415523],
[ 0.08002874, -0.21151477, 0.206062 , 0.10593122],
[ 0.33365312, 0.07702804, 0.06801722, -0.00739034],
[-0.02813852, -0.07401617, 0.00837274, -0.16002993],
[ 0.13233314, 0.17294073, 0.13590576, -0.1433093 ]], dtype=float32)
Reducing Dimensionality With Basis#
Alternatively, one can use basis to reduce the number of parameters (useful when you have a history filter the spans hundreds of samples). Reducing the dimensionality of the features helps to prevent overfitting and speeds up computation, particularly for long history windows.
raised_cos = nmo.basis.RaisedCosineLogConv(3, window_size=3)
X2 = raised_cos.compute_features(counts)
# 12 Features: 3 basis funcs x 4 neurons
print("Shape of the feature matrix: ", X2.shape)
plt.title("Basis Feature Matrix", fontsize=16)
plt.imshow(X2[10:], cmap="Greys")
plt.xticks([1, 4, 7, 10], [0, 1, 2, 3])
plt.xlabel("Neurons", fontsize=12)
plt.ylabel("Samples", fontsize=12)
plt.show()
Shape of the feature matrix: (100, 12)
This model can be fit with the usual model_reduced = nmo.glm.PopulationGLM().fit(X2, counts).
Interpreting the coefficients#
The learned model coefficients are stored in a 2D array of shape (n_features, n_neurons). This array concatenates the coefficients representing pairwise couplings along the first dimension. Using the split_by_feature method of the basis object, the coefficients can be reshaped into a 3D array of shape (n_neurons, n_basis_funcs, n_neurons), where the first dimension corresponds to sender neurons, the second dimension contains the basis function coefficients, and the third dimension corresponds to receiver neurons.
# get the dictionary (1 key per basis component, here there is single basis component)
split_coef = basis.split_by_feature(model.coef_, axis=0)
print(split_coef["spike-history"].shape)
(4, 10, 4)
This makes it easy to retrieve the coefficients capturing how the spike history of the neuron i affects the firing of neuron j,
sender_neuron_i = 1
receiver_neuron_j = 2
coeff_ij = split_coef["spike-history"][sender_neuron_i, :, receiver_neuron_j]
coeff_ij
Array([-0.22478367, -0.09167013, 0.16141443, -0.17483383, 0.00090788,
-0.05813771, -0.0402254 , -0.2647168 , -0.10762664, 0.10566666], dtype=float32)
Selecting The Connectivity Map#
It is also possible select which couplings we want to learn by specifying a boolean mask (a mask of 0s and 1s) of shape (n_features, n_neurons). The mask will be elementwise multiplied to the coefficients, selecting which will be included in the GLM.
For instance, if mask[i, j] == 0, the i-th feature won’t be included in the GLM fit of the j-th neuron.
For example, we can exclude the bi-directional coupling between neuron 0 and 1, and the directional coupling between neuron 2 (sender) and neuron 3 (receiver).
# initialize a neuron x neuron mask
mask = np.ones((n_neurons, n_neurons))
# remove the bi-directional coupling between 0 and 1
mask[0, 1] = 0
mask[1, 0] = 0
# remove the directional coupling from 2 to 3
mask[2, 3] = 0
# repeat over the rows by the number of basis functions
mask = np.repeat(mask, basis.n_basis_funcs, axis=0)
We are now ready to fit the GLM masking the coefficients.
model_with_mask = nmo.glm.PopulationGLM(feature_mask=mask).fit(X, counts)
print(model_with_mask.coef_)
[[ 0.04630114 0. -0.23782578 0.06813474]
[-0.27391973 0. -0.07002713 0.14300758]
[-0.2698096 0. 0.28980285 0.20106782]
[ 0.0633865 0. 0.16020632 -0.06749498]
[-0.06823029 0. -0.11079037 0.08256442]
[-0.08302427 0. 0.15963702 -0.1068478 ]
[-0.15966156 0. -0.11132955 -0.00046913]
[ 0.03471004 0. -0.07757225 -0.3169972 ]
[-0.09083457 0. 0.10234725 -0.10740117]
[-0.18620567 0. -0.2009712 0.00472539]
[ 0. 0.05497848 -0.22489585 0.03454009]
[ 0. -0.05261067 -0.09183352 -0.08446524]
[ 0. 0.15176393 0.16126904 -0.09870093]
[ 0. 0.01657126 -0.17475788 0.08635148]
[ 0. 0.22074749 0.00067707 0.17821233]
[ 0. -0.08273246 -0.05826272 0.05252425]
[ 0. -0.03634679 -0.0403231 -0.10383549]
[ 0. -0.21660514 -0.2648237 0.02635365]
[ 0. 0.06628797 -0.10806276 -0.12586506]
[ 0. -0.16933341 0.10577829 -0.09974751]
[ 0.12467519 0.23134491 -0.17698134 0. ]
[ 0.09236328 0.07699226 0.22820951 0. ]
[ 0.13449037 -0.07053307 -0.00875364 0. ]
[-0.13304639 -0.06465386 -0.39669162 0. ]
[ 0.29952595 -0.15363216 -0.01847102 0. ]
[ 0.16719241 -0.1149743 0.18232046 0. ]
[ 0.13905984 -0.03036673 0.01050939 0. ]
[-0.06460074 -0.003047 -0.23702075 0. ]
[-0.24841776 0.23920667 0.08470093 0. ]
[ 0.24449126 -0.01037522 0.19145684 0. ]
[ 0.06033821 0.01672796 -0.11526815 -0.17216173]
[ 0.10661148 -0.08722092 -0.01869801 -0.04801634]
[ 0.00245442 -0.00621818 0.02852625 0.07815327]
[ 0.08082574 0.08473439 -0.00226621 0.2032633 ]
[-0.12612696 0.193724 0.28753555 0.09120496]
[-0.01689474 -0.00863638 -0.2759687 0.26601073]
[ 0.22023158 -0.10151734 0.20594643 0.13836515]
[ 0.11793084 0.09283365 0.06813051 0.02484599]
[ 0.09320085 -0.18433885 0.00833217 -0.22409207]
[ 0.04329374 0.07402241 0.1360465 -0.13390823]]
Check that the coefficient for the connection we removed, are actually null.
split_coef = basis.split_by_feature(model_with_mask.coef_, axis=0)["spike-history"]
# check sender=2 receiver=3
print("Coupling 2 -> 3: ", split_coef[2, :, 3], "\n")
# check that we kept the sender=3, receiver=2
print(f"Coupling 3 -> 2: ", split_coef[3, :, 2])
Coupling 2 -> 3: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Coupling 3 -> 2: [-0.11526815 -0.01869801 0.02852625 -0.00226621 0.28753555 -0.2759687
0.20594643 0.06813051 0.00833217 0.1360465 ]