Repository for MATLAB code used for information theoretic analysis of neural signals. In particular, this software computes Mutual Information in the Frequency domain (MIF) (Malladi et al., 2018) between time series data. The code uses a new low-variance approach that is improved for correlating MIF between neural signals with relevent experimental measurements, such as the accuracy of a monkey during a task.
The key function that you can use for MIF estimation is helper_functions/miCompute.m, which takes in 6 arguments:
- x - input data for first signal in shape (time)X(trials/windows)
- y - input data for second signal in shape (time)X(trials/windows)
- Fs - sampling frequency in Hz
- Fm - maximum frequency to compute MIF for in Hz (e.g. Fm = 100 Hz means MIF will not be computed above 100 Hz)
- K - number of neighbors for k-nearest neighbors (k-nn) mutual information estimate (Kraskov et al., 2004), which should be half the number of samples for decreasing variance (Kraskov et al., 2004)
- Np - number of permutations for significance testing
Using an approach based on a prior work (Pantazis et al., 2005), significance testing relies on using a threshold for MIF values equal to the 95th percentile of a maximum distribution composed of maximum MIF samples across frequency for each permutation. In particular, trials/windows of x are permuted and then MIF is computed, from which the maximum MIF value across frequencies is contributed to the maximum distribution.
The output of the function is two variables:
- MI - vector of MIF values
- f - vector of frequencies corresponding to MIF values
Therefore, a full call would be: [ MI, f ] = miCompute( x, y, Fs, Fm, K, Np )
We provide an example file miExampleScript.m at the root of the repository that calls this function to compute the MIF between two simulated random processes which are composed of random sinusoids to reflect neural oscillations. This is inspired by prior work (Malladi et al., 2018). For this model, MIF is analytically known and the estimates of our implementation can be compared with those true MIF values. The script displays a scatter plot comparing these estimated and true MIF values, along with the correlation.
The three sinusoids we consider are as follows:
X(t) = Acos(2π ƒ0t+Θ)
W(t) = Bcos(2π ƒ0t+Φ)
Y(t) = X(t) + W(t)
where A ~ Rayleigh(σ_A = 1), B ~ Rayleigh(σ_B), Θ ~ uniform(0, 2π), and Φ ~ uniform(0, 2π).
Accordingly, the MIF MI_XY(ƒ0,ƒ0) is computed between X and Y for ƒ0, and is known to be (Malladi et al., 2018):
MI_XY(ƒ0,ƒ0) = log(1 + (S_X(ƒ0)/S_W(ƒ0))) = log(1 + (1/σ_B))
when σ_A = 1.
Therefore, the true MI values used in the example script are randomly selected over a range by varying σ_B. The resulting scatter plot shows the relationship between these true MI_XY(ƒ0,ƒ0) values and estimated values produced by our technique.
This software has been successfully tested in MATLAB R2020a on Windows 10, macOS Catalina, and Ubuntu 20.04. Pre-compiled mex files are included (compiled using modified code from TSTOOL (DPI Göttingen, 2009)), and the helper_functions/load_mex_directory.m script automatically adds the OS-relavant folder containing these mex files.
For future compatibility, if the pre-compiled mex files no longer work in a newer release of MATLAB, we have included a function tsool_files/makemex.m (modified from TSTOOL (DPI Göttingen, 2009)) that automatically will recompile the mex files from source and overwrite the existing mex files in the OS-appropriate folder. To utilize this function, change your directory to tstool_files/ and then in MATLAB call "makemex" with no input arguments.
The key idea of this augmentation to the prior MIF approach (Malladi et al., 2018) is that we capitalize on the mutlitaper method (Thomson 1982) to reduce estimation variance and sharpen the capability of MIF estimation in correlation analyses. In particular, we implement what we call the "post" approach, where MI_XY(ƒ0,ƒ0) is computed separately for each taper and then the average MIF value taken across tapers is treated as the estimated MI_XY(ƒ0,ƒ0). Each MI estimate is performed via k-nn MI estimation (Kraskov et al., 2004).
We thank DPI Göttingen for TSTOOL (DPI Göttingen, 2009), which is a software package that we have included some of the code from in order to perform k-nn estimation of MI (Kraskov et al., 2004). We have modified the code to produce up-to-date mex files, and both the modified code and mex files are included in the tstool_files/ directory. Since their code was licensed as GPLv2, we have also adopted the GPLv2 license for this toolbox (see LICENSE and GPLv2_note.txt files).
- Malladi et al. Mutual information in frequency and its application to measure cross-frequency coupling in epilepsy. IEEE Transactions on Signal Processing, 66, 2018.
- Kraskov et al. Estimating mutual information. Physical Review E, 69, 2004.
- Pantazis et al. A comparison of random field theory and permutation methods for the statistical analysis of MEG data. NeuroImage, 25, 2005.
- Thomson. Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70, 1982.
- DPI Göttingen. TSTOOL. University of Göttingen, Göttingen, Germany, 2009. http://www.dpi.physik.uni-goettingen.de/tstool/.