How to use the code: Please use this code for ethical pursuits only, see lisence.
What it does: This code takes a set of timeseries (e.g. parcellated fMRI data) and clusters them into phase consistent communities at each timepoint based on the instantaneous phase of the mode decomposed timeseries. Each cluster label references a limited phase range across time and subjects. See below for details.
If you use this code, please cite "Time-resolved brain network community detection based on instantaneous phase of fMRI data", doi: https://doi.org/10.64898/2026.03.17.712372
1.Band-pass filtering and Variational Mode Decomposition: BandPassAndVMD.m
Data is band-pass filtered based on the choice of lowpass and highpass thresholds. This step is optional but reduces noise in the data which in turn reduces mode mixing in the next step. In the reference manuscript, a band width of 0.01–0.15 Hz was used (TR = 0.72). Five time points will be cropped at the start and end of the time series to remove potential edge effects from bandpass filtering.
The band-pass filtered data is decomposed into modes using the VMD algorithm. The number of modes must be defined prior to the decomposition. In the manuscript, two modes were used (see the manuscript for details).
Run as: BandPassAndVMD(data, TR, lowpass, highpass, M) Requires VMD.m (included in MATLAB 2025b, else you can find it here: https://se.mathworks.com/matlabcentral/fileexchange/44765-variational-mode-decomposition)
data = pre-processed parcellated data TR = the time repetition of your data lowpass = frequency for the lowpass filter highpass = frequency for the highpass filter M = number of modes for the VMD decomposition
2.Phase scrambled the mode timeseries and preform the same analysis as for the real data to establish statistical thresholds for significant results. iaaft_batch_scramble.m (requires iaaft_1d.m). Run as: iaaft_batch_scramble(VMD_out) VMD_out = mode decomposed timeseries, a.k.a intrinsic mode functions (imfs)
3.For both datasets (real and phase scrambled), for each mode run the code that generate the communities: Run_Phase_Cluster_Analysis.m. Five additional timepoints are cropped at the start and end of the mode timeseries to remove edge effects created by the VMD. Thus, a total of 10 time points at the beginning and end of the timeseries will be lost. If the original length of the raw timeseries was T, it will now be T-20. Depending on your data and choices of parameters, you may be able to crop less.
Run as: Run_Phase_Cluster_Analysis (VMD, TR, mode, thr1, seg, perm_nr, parcellation) Requires : AddCom_permute.m AddComSingle_permute.m FindCom_NoOpp_permute.m FindCom_permute.m PhaseCom_NoOpp_permute.m PhaseCom_permute.m relabel_matrix.m TimeResolvedNets_permute.m
VMD = mode-decomposed data (S, N, T, M), where S = number of subjects, N = number of parcels, T= time (after the two cropping steps above), and M = number of modes (i.e. intrinsic mode functions a.k.a. IMFs).
TR = the time repetition of your data
Mode = which mode you are analysing.
Perm_nr = the number of permutations if you want to improve the optimization for the largest among several potential clusters.
Thr1 = the integration threshold where thr1=cos〖(θ_1)〗. Abs(θ_1) is the maximum difference in instantaneous phase between any two parcels in any community. The manuscript focused on thr1 = 0.5 which corresponds to θ_1= π/3 =60°.
Seg = the segregation threshold. If you want the standard use (i. e. thr2= -thr1), set seg = 50. It can also be omitted. If you want to omit it set thr2 =1000. Else set thr2 to your choice in the range [-1,1] while thr2 < thr1.
Parcellation = a vector of network numbers for all the parcels.
Main group level outputs as described in the paper:
Com1.mat, Com2.mat …. Etc (Phase distribution of parcels across time and subject for each cluster)
AccumSize.mat
SizesCom.mat (containes “MeanSizeCom”, “SizeComAll”, “ComSortedAll”). “MeanSizeCom” is the average (across subjects) community size in terms of fraction of total parcels. Used for whole-brain time-resolved integration/segregation analysis.
Time-resolved Integration/Segregation for individual networks are found in: RelativeComSizeActivation.mat MeanComNet Requires NetName.mat which is a table where the first column is the number of the network and the third column is the size in number of parcels.
CrossSubjectSynchronization.mat contains “Cross_subj_sync”, “Parcel_Com_All” and “miniall”. Cross_subj_sync is the synchronization timeseries for all parcels in the order of the parcellation. Parcel_Com_All contains the subjects specific parcel community allegiances across time. Miniall is only for plotting purposes in python.
SortedSYNC.mat contains SortedSync and MaxiTime. SortedSync are the cross-subject-synchronization timeseries for all parcels sorted by networks. For the Schaefer parcellation, left hemisphere parcels are followed by right hemisphere parcels within each network. If your parcellation from the start is sorted by network primarily and not hemispheres, then SortedSync is the same as Cross_subj_sync. MaxiTime is the maximum cross-subject-synchronization per parcel and time across communities. Sorted in the same way as SortedSync.
For questions contact: marika.strindberg@ki.se