NeuroAnalysis - Neural Signal Analysis

Predefined ColorMaps

Default digital input channels

Time Unit set to millisecond by default.

Generates linearly interpolated colors between key colors inclusively.

return an Array of colors.

alpha weighted average of foreground and background colors

Masking alpha channel of an image, match the implementation in Experica shaders.

Try to locate cell layer.

1. y coordinate of cell postion
2. layer definition

Weighted average of banded masked matrix

2D cas function defined as $cas(x+y) = cos(x+y) + sin(x+y)$

• kx: Frequency in cycle/unit_x
• ky: Frequency in cycle/unit_y
• phase: Phase of a cycle in [0, 1] scale
• isnorm: scale cas in [-√2, √2] to [-1, 1]

cas function defined as $cas(x) = cos(x) + sin(x)$

• f: Frequency in cycle/unit_x
• phase: Phase of a cycle in [0, 1] scale
• isnorm: scale cas in [-√2, √2] to [-1, 1]

2D cas to 2D sin gratingf

cas phase to sin phase, phase is in [0, 1] scale

Check circuit consistency, remove conflicts and duplicates

Check Layer Boundaries

Logical mask for Neuropixels channels in probe shape 1 2 3 4 5 6 ...

probe shape channel positions

Map linear channel index to probe shape indices

Properties of Circular Tuning

- Prefered Angle with Peak Response
- Selectivity Index
    - version 1: (ResponsePeak - ResponseOpposing)/ResponsePeak
    - version 2: (ResponsePeak - ResponseOpposing)/(ResponsePeak + ResponseOpposing)
- Half Width at Half Peak-to-Trough

1. x: angles in radius
2. y: responses

• od: opposing angle distance to prefered angle, e.g. π for DSI, 0.5π for OSI
• fn: factor name

Confine ROI so that it is not out of the image

clamp value to min and max, and linearly map range [min, max] to [0, 1]

Multi-Taper channel pair-wise coherence spectrum estimation

Complex sum of angles and corresponding image responses

1. image responses [height, width, n]
2. angles in radius [n]

• centerfun: function to get center for sdfactor (median)
• presdfactor: center +- sdfactor * sd, for clamping before filtering
• filter: filter kernel
• sufsdfactor: center +- sdfactor * sd, for clamping after filtering
• responsesign: increasing/decreasing(+/-) response

return complex map, angle map, and magnitude map

Get factors of conditions

Mean image responses of each conditions

1. rs: responses of each epoch [height, width, nepoch]
2. ci: epoch indices of repeats for each condition

In Condition Tests, find each unique condition, number of repeats and its indices

Group repeats of Conditions, get Mean and SEM of responses

1. rs: responses of each trial
2. ci: trial indices of repeats for each condition

Print Condition in String

Get CondTest DataFrame

Get CondTestCond DataFrame

Normalized(coincidence/spike), condition/trial-averaged Cross-Correlogram of binary spike trains.

• Correction: shuffle=true

  (Bair, W., Zohary, E., and Newsome, W.T. (2001). Correlated Firing in Macaque Visual Area MT: Time Scales and Relationship to Behavior. J. Neurosci. 21, 1676–1697.)

          jitter

  ()

1D Current Source Density for voltage traces from a equidistant linear array of electrodes.

1. data: ch x sample or ch x sample x epoch(volts)

2. method:

• CSD: finite difference approximation of second spatial derivatives of potentials, C = - ∇⋅σ∇ϕ, σ is conductivity tensor, ϕ is potential. (Nicholson & Freeman, 1975, J Neurophysiol, 38(2): 356-68)

• iCSD: solve C -> ϕ forward model, then get Ĉ by the inverse transformation matrix, here C is modeled by δ-source, i.e. infinitely thin current source disc of radius r. (Petterson et al., 2006, J Neurosci Methods, 154(1-2):116-33)

• kCSD: kernal CSD ()

• h: channel spacing(micrometer)

• c: conductivity of extracellular medium(siemans/meter)

• r: radius(micrometer) of δ-source

• return: CSD(amps/meter^3)

Get CondTest and CondTestCond DataFrame

Discrete Fourier Transform at frequencies

\[DFT[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{-2\Pi ink}{N}}, k = 0:N-1, N = length(x), fₛ sampling x, fₖ = fₛ/N sampling DTFT\]

1. signal
2. simpling frequency of signal
3. at which frequencies DFT are directly evaluated

Parse bit event time and value.

Get digital channel time and value.

Direction tuning significance using dot product with empirical orientation preference This function calculates the probability that the "true" direction tuning vector of a neuron has non-zero length. It performs this by empirically determing the unit orientation vector, and then computing the dot product of the direction vector for each trial onto the overall orientation vector, and then looking to see if the average is non-zero.

Inputs: ANGLES is a vector of direction angles at which the response has been measured. RATES is the response of the neuron in to each angle; each row should contain responses from a different trial. Output: P the probability that the "true" direction tuning vector is non-zero.

See: Mazurek, Kagan, Van Hooser 2014; Frontiers in Neural Circuits

concentric circular dog contour

concentric circular dog envelope

Binary mask for concentric circular dog envelope

Difference of Gaussians function

Drop array dims according to MATLAB convention.

• scalar: scalar instead of 1x1 array
• rvector: vector instead of 1xN array
• cvector: vector instead of Nx1 array

concentric orientated elliptical dog contour

concentric orientated elliptical dog envelope

Binary mask for concentric orientated elliptical dog envelope

Ellipse function

Convert epochs in time to sample index

Epochs of Channel Sample.

Epochs of Neuropixels data stream (Channels x Samples), optionally gain corrected(voltage), line noise(60,120,180Hz) removed, bandpass filtered and common average referenced

Epochs of a Spike Train in between binedges

Epoch of a Spike Train x, where min <= x[i] < max.

kwargs isminzero, ismaxzero and shift set epoch zero value to min+shift or max+shift.

return:

• y: epoch of the Spike Train
• n: number of elements in the epoch, or divided by duration(max-min) of the epoch, based on kwargs israte
• w: epoch window(min,max)
• i: indices of epoch elements in the original Spike Train, such that y = x[i]

See also: epochspiketrains, epochspiketrainresponse

Response of each epoch of a Spike Train, could be mean firing rate or number of spikes based on kwarg israte.

See also: epochspiketrainresponse_ono

Response of each epoch of a Spike Train, could be mean firing rate or number of spikes based on kwarg israte.

Epochs for each Spike Trains

Logical mask for Neuropixels excluded channels in probe shape

extrema value of max amplitude and it's delay index

Regular Expression to match Experica data file name

Estimate the F1 Ori and SpatialFreq of an image from its 2D powerspectrum.

1. x: 2D powerspectrum
2. freqd1: frequencies of dim 1
3. freqd2: frequencies of dim 2

return:

• ori: Orientation in radius[0,π], 0 is -, increase counter-clock wise
• sf: SpatialFreq along the line perpendicular to ori

Condition Response in Factor Space

Tuning properties of factor response

1. fl: factor levels

2. fr: factor responses for each level

   Angle, Orientation and Direction follow the same convention such that 0 is -/→, then increase counter-clock wise.

   For cases where Orientation and Direction are interlocked(drifting grating): - when Orientation is -(0), then Direction is ↑(90) - when Direction is →(0), then Orientation is |(-90)

fill data in the shape of mask, where masked channels are replaced with local average

Get Final factors of conditions

find closest distance and index in α to β, all in radius

Find Condition with Factor = level

Fit 1D model to data

Fit 2D model to image

Flat Spike Trains to SpikeTimes and Trials, optionally sort trial based on sv

Find levels for each factor and indices, repetition for each level

Frame Response

gabor contour

gabor envelope

Binary mask for gabor envelope

Gabor function

Convert raw Neuropixels data to gain-corrected voltages. The saved-channel id subset in data will be used to get corresponding gain.

The full conversion with gain is:

$dataVolts = dataInt * fi2v / gain$

Note that each channel may have its own gain.

1. y: the saved data stream, excluding the last Sync channel.
2. meta: corresponding meta for y

Gaussian contour

Gaussian envelope

Gaussian envelope mask

Gaussian function

Get matched file names in path

Get matched Experica files in directory

Get Optical Imaging VDAQ matched files in directory

Goodness of Fit Metrics:

• r: Pearson Correlation Coefficient
• mae: Mean Absolute Error
• rmse: Root Mean Squared Error
• rae: Relative Absolute Error
• rse: Relative Squared Error
• r2: R Squared
• adjr2: Adjusted-R²
• s: Residual Standard Error
• aic: Akaike Information Criterion
• bic: Bayesian Information Criterion

1. y: responses
2. ŷ: model predictions

• n: sample size
• e: errors(y - ŷ)
• k: number of predictors
• df: degree of freedom(n - k - 1)

Generate Grating Image, match the implementation in Experica grating shader.

• θ: Orientation (radius), 0 is -, increase counter-clock wise
• sf: SpatialFreq (cycle/deg)
• tf: TemporalFreq (cycle/sec)
• t: Time (second)
• phase: Phase of a cycle in [0, 1] scale
• size: Tuple of image size in degree
• ppd: pixel per degree
• isnorm: return image in [0, 1] or [-1, 1]

2D sin grating function

• μ₁: x offset
• μ₂: y offset
• θ: Orientation in radius, 0 is -, increase counter-clock wise
• f: Frequency in cycle/unit_x/y
• phase: Phase of a cycle in [0, 1] scale

sin grating function

• μ: x offset
• f: Frequency in cycle/unit_x
• phase: Phase of a cycle in [0, 1] scale

Generalized von Mises function ^[1]

\[f(α) = βe^{κ₁(cos(α - μ₁) - 1) + κ₂(cos2(α - μ₂) - 1)}\]

^[1]

Gatto, R., and Jammalamadaka, S.R. (2007). The generalized von Mises distribution. Statistical Methodology 4, 341–353.

left and right half width of v relative to y[start], -Inf/Inf when no v is found.

• start: index of y which is the center of the width
• v: value on which width is cutoff
• circ: whether y is defined on circular domain and width can wrap around
• x: domain of y, return width when x provided, otherwise return cutoff indices

Generate gratings in Hartley space (PL)

• kx, ky: wavenumber along x, y axis
• bw: black and white (phase) flip
• stisize: stimulus size in visual angle (degree); note that sz[1]=sz[2]

if norm=true,Return image in [0,1], otherwise return image in [-1,1]

Generate Hartley Subspace, where k is frequency in cycle/unit_x/y. ^[1]

• kbegin: k >= kbegin
• kend: k <= kend
• dk: Δk, step on k axis
• phase: phase in [0, 1] scale, default 0.
• shape: :square or :circle shape subspace
• addhalfcycle: add half cycle shifted hartleys
• blank: the element of hartley as blank, default uniform gray.
• nblank: number of blanks to add

^[1]

Ringach, D.L., Sapiro, G., and Shapley, R. (1997). A subspace reverse-correlation technique for the study of visual neurons. Vision Research 37, 2455–2464.

High pass and low pass filtering

Hotelling's T-Squared test for one multivariate sample.

Performs Hotelling's T^2 test on multivariate samples X to determine if the data have mean MU. X should be a NxP matrix with N observations of P-dimensional data, and the mean MU to be tested should be 1xP. the significance level α is 0.05 by default.

H is 1 if the null hypothesis (that the mean of X is equal to MU) can be rejected at significance level α. P is the actual P value.

The code is based on HotellingT2.m by A. Trujillo-Ortiz and R. Hernandez-Walls.

Check if any group in response is significently different from at least one other group in response

Check if any sub group of response is significently different from baseline by Wilcoxon Signed Rank Test

Check if response is significently different from baseline by Wilcoxon Signed Rank Test

Check if std of a spatial-temporal kernal within response time window significently higher than that of the baseline time window

Load images in directory

Weighted average of local channel pair coherences for each channel

local RMS contrast of each image, highlighting local structure regions

local channel pairs and gaussian weights within a circle centered for each channel

Map time to data timeline, optionally add latency.

Get matched file names in directory, optionally join directory to get file path

convert Matlab struct of array to DataFrame

Normalize values of an array to be between -1 and 1

original sign of the array values is maintained.

Regular Expression to match Optical Imaging VDAQ block file name

Get digital rising and falling edges in a analog stream return di: edges index dv: edges active/inactive state

ROI encompass peak value and its delay index

Plot Analog Signals

Plot Mean and SEM of responses for each condition

plot hartley subspace

scatter plot of spike trains

Generate a Poisson Spike Train

Multi-Taper power spectrum estimation

2D powerspectrum of an image.

2D powerspectrums of images in same size.

Read and Prepare Dataset in MATLAB MAT File

Prepare Experica Experiment Data

Prepare Optical Imaging Block Data

Prepare Ripple Data

excitatory and inhibitory projections between two spiking units from cross-correlogram

PSTH of Spike Trains

Read variables of a MATLAB MAT file into a Dictionary

1. f: MAT file path
2. vars: variable names in the varset to read

• varset: variable names to select in vars
• scalar: scalar instead of 1x1 matrix
• rvector: vector instead of 1xN matrix
• cvector: vector instead of Nx1 matrix

Read DataExport Metadata MAT File

Read list of files encoding raw mono 8bit image

Map time between ref and target, according to the timing diff between target and ref of the same sync signal

Logical mask for Neuropixels reference channels in probe shape

Remove line noise and its harmonics by notch filter

This function calculates the ROC curve, which represents the 1-specificity and sensitivity of two classes of data, (i.e., class1 and class2).

The function also returns all the needed quantitative parameters: threshold position, distance to the optimum point, sensitivity, specificity, accuracy, area under curve (AROC), positive and negative predicted values (PPV, NPV), false negative and positive rates (FNR, FPR), false discovery rate (FDR), false omission rate (FOR), F1 score, Matthews correlation coefficient (MCC), Informedness (BM) and Markedness; as well as the number of true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN).

Rewrite from Víctor Martínez-Cagigal (2020). ROC Curve, MATLAB Central File Exchange. Retrieved February 11, 2020.

Get ROI from region indices

sample index to time

Find the extreme value

Peichao:

Join the results from differernt Hartleys after hartleyFourier analysis

peak ROI region and its delay

Properties of Spatial Frequency Tuning

- Prefered Spatial Frequency with Peak Response
- Half Width at Half Peak-to-Trough
- Freq Passing Type {A:All Pass, H:High Pass, L:Low Pass, B:Band Pass}
- Bandwidth ``log2(H_cut/L_cut)``
- Passwidth at Half Peak-to-Trough constrained by low freq lim and high freq lim

1. x: sf in cycle/degree
2. y: responses

2D sin gratingf to 2D cas

sin phase to cas phase, phase is in [0, 1] scale

Spike Triggered Average of Images

1. x: Matrix where each row is one image
2. y: Vector of image response

• norm: normalization factor, default no normalization. it could be $sum(y)$ if y is number of spike or spike rate, then STA would be spiking probability.
• whiten: whiten factor, default no whiten. it could be $(xᵀx)⁻¹$ or inverse of covariance matrix, that decorrelate STA.

Extract state time of condtest

Spatial and Temporal filter of data stream

time to sample index

Organize each spiking unit from Kilosort result.

Organize each spiking unit from Kilosort3 result.

Regular expression to match VLab data file names

Vertical Mean and SEM of a matrix

von Mises function ^[1]

\[f(α) = βe^{κ(cos(n(α - μ)) - 1)}\]

^[1]

Swindale, N.V. (1998). Orientation tuning curves: empirical description and estimation of parameters. Biol Cybern 78, 45–56.

• β: amplitude at μ
• μ: angle of peak
• κ: width parameter
• n: frequency parameter

Vertical stack same length vectors to matrix

Generates linearly interpolated ColorGradient

Generates linearly interpolated ColorGradient between key colors.