diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b5d19a6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,93 @@ +# MATLAB +*.asv +*.m~ +*.mex* +*.autosave + +# Results and Data directories (too large for git) +Results/ +Dataset/ +*.mat +*.nii +*.nii.gz + +# Brainstorm database +brainstorm_db/ +*brainstorm*/ +protocols/ + +# FreeSurfer outputs +subjects/ +fsaverage/ + +# Temporary files +*.tmp +*.temp +*.log +*.swp +*~ +.DS_Store +Thumbs.db + +# Python +__pycache__/ +*.py[cod] +*$py.class +*.so +.Python +env/ +venv/ +*.egg-info/ +dist/ +build/ + +# R +.Rhistory +.RData +.Rproj.user + +# IDE +.vscode/ +.idea/ +*.sublime-* + +# OS +.DS_Store +.DS_Store? +._* +.Spotlight-V100 +.Trashes +ehthumbs.db +Desktop.ini + +# Documentation builds +_build/ +site/ + +# Backup files +*.bak +*.backup + +# Large data files +*.fif +*.sqd +*.con +*.raw +*.eve +*.cov +*.trans + +# Specific to this project - keep these patterns +# But allow small example/test data if needed +!example_data.mat + +# Logs and reports +*.pdf +*.png +*.jpg +*.jpeg +*.svg +# But keep documentation images +!docs/*.png +!docs/*.jpg +!docs/*.svg diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..e5decd6 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,121 @@ +# Changelog + +All notable changes to the AD_PAC project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), +and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [Unreleased] + +### Added +- Comprehensive documentation suite + - Main README.md with project overview + - CITATION.md with preprint reference + - REQUIREMENTS.md with detailed dependency information + - QUICKSTART.md for rapid onboarding + - CONTRIBUTING.md with contribution guidelines + - Module-specific READMEs for each pipeline stage +- LICENSE file (MIT License) +- .gitignore for keeping repository clean +- CHANGELOG.md for tracking changes + +### Changed +- Initial documentation release + +### Deprecated +- None + +### Removed +- None + +### Fixed +- None + +### Security +- None + +## [1.0.0] - 2026-02-16 + +### Added +- Initial release of AD_PAC pipeline +- Head modeling pipeline (anatomyTransform, buildHeadModels, processBSFiles, removeElec) +- Preprocessing pipeline (preprocessing, interpolate_bad_channels, detect_bad_trials, plot_data) +- Source reconstruction pipeline (sourceReconstruction_p, lcmv_meg, and helper functions) +- PAC analysis pipeline (seedRoitoCortex, er_pac, er_pac_3, test_pac) +- Statistical analysis pipeline (prepareTable, powerAnalysis, clusterBasedModel, linearMixedEffectsModel) +- Python scripts for visualization (Fig_creator.py, cluster_based_perm.py) +- Supporting data files (roi_name.mat, dk_labels.mat, cm17.mat, mask.mat) + +## Release Notes + +### Version 1.0.0 (2026-02-16) + +First public release of the AD_PAC pipeline. This version includes: + +**Features:** +- Complete pipeline from raw MEG data to statistical analysis +- Support for CTF MEG systems (275 channels) +- Integration with FieldTrip, Brainstorm, and FreeSurfer +- Cluster-based permutation testing +- Linear mixed-effects modeling +- Comprehensive documentation + +**Data:** +- Compatible with data from OSF repository: https://osf.io/pd4h9/overview +- Supports Desikan-Killiany cortical parcellation +- Analyzes theta/alpha - low-gamma PAC + +**Requirements:** +- MATLAB R2018b or later +- FieldTrip toolbox +- Brainstorm toolbox +- FreeSurfer 6.0 or later + +**Known Limitations:** +- Currently optimized for CTF MEG systems +- Requires substantial computational resources +- FreeSurfer preprocessing can take 6-24 hours per subject + +**Citation:** +If you use this release, please cite: +- Preprint: https://www.medrxiv.org/content/10.64898/2026.02.06.26345635v1 + +--- + +## Future Plans + +### Planned for Version 1.1.0 +- [ ] Support for additional MEG systems (Neuromag, KIT) +- [ ] GPU acceleration for PAC computation +- [ ] Additional PAC metrics (MI, MVL, GLM) +- [ ] Automated quality control metrics +- [ ] Example dataset for testing + +### Planned for Version 1.2.0 +- [ ] Web-based visualization interface +- [ ] Batch processing utilities +- [ ] Docker container for easier deployment +- [ ] Integration with BIDS standard + +### Long-term Goals +- [ ] Real-time PAC analysis capabilities +- [ ] Machine learning classification features +- [ ] Support for combined MEG/EEG analysis +- [ ] Cloud computing integration + +--- + +## Contributing + +See [CONTRIBUTING.md](CONTRIBUTING.md) for guidelines on how to contribute to this project. + +## Versioning + +We use [Semantic Versioning](https://semver.org/): +- MAJOR version for incompatible API changes +- MINOR version for new functionality in a backwards compatible manner +- PATCH version for backwards compatible bug fixes + +## License + +This project is licensed under the MIT License - see [LICENSE](LICENSE) file for details. diff --git a/CITATION.md b/CITATION.md new file mode 100644 index 0000000..519509d --- /dev/null +++ b/CITATION.md @@ -0,0 +1,105 @@ +# Citation + +## Primary Citation + +If you use this code or methodology in your research, please cite our preprint: + +**Preprint Paper:** +- **URL**: [https://www.medrxiv.org/content/10.64898/2026.02.06.26345635v1](https://www.medrxiv.org/content/10.64898/2026.02.06.26345635v1) +- **DOI**: 10.64898/2026.02.06.26345635 + +### BibTeX Entry + +```bibtex +@article{AD_PAC_2026, + title={Phase-Amplitude Coupling Analysis in Alzheimer's Disease}, + author={[Authors to be filled from paper]}, + journal={medRxiv}, + year={2026}, + doi={10.64898/2026.02.06.26345635}, + url={https://www.medrxiv.org/content/10.64898/2026.02.06.26345635v1} +} +``` + +## Additional Citations + +### PAC Method + +The PAC estimation code is based on the method developed by: + +Pellegrini, F., Delgado Saa, J., & Haufe, S. (2023). Identifying good practices for detecting inter-frequency coupling in resting-state EEG: A simulation study. *NeuroImage*, 277, 120233. + +```bibtex +@article{pellegrini2023identifying, + title={Identifying good practices for detecting inter-frequency coupling in resting-state EEG: A simulation study}, + author={Pellegrini, Franziska and Delgado Saa, Josefina and Haufe, Stefan}, + journal={NeuroImage}, + volume={277}, + pages={120233}, + year={2023}, + publisher={Elsevier} +} +``` + +**Code Repository**: [https://github.com/fpellegrini/PAC](https://github.com/fpellegrini/PAC) + +### Required Toolboxes + +If you use the toolboxes employed in this pipeline, please also cite: + +#### FieldTrip + +Oostenveld, R., Fries, P., Maris, E., & Schoffelen, J. M. (2011). FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. *Computational Intelligence and Neuroscience*, 2011. + +```bibtex +@article{oostenveld2011fieldtrip, + title={FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data}, + author={Oostenveld, Robert and Fries, Pascal and Maris, Eric and Schoffelen, Jan-Mathijs}, + journal={Computational Intelligence and Neuroscience}, + volume={2011}, + year={2011}, + publisher={Hindawi} +} +``` + +#### Brainstorm + +Tadel, F., Baillet, S., Mosher, J. C., Pantazis, D., & Leahy, R. M. (2011). Brainstorm: a user-friendly application for MEG/EEG analysis. *Computational Intelligence and Neuroscience*, 2011. + +```bibtex +@article{tadel2011brainstorm, + title={Brainstorm: a user-friendly application for MEG/EEG analysis}, + author={Tadel, Fran{\c{c}}ois and Baillet, Sylvain and Mosher, John C and Pantazis, Dimitrios and Leahy, Richard M}, + journal={Computational Intelligence and Neuroscience}, + volume={2011}, + year={2011}, + publisher={Hindawi} +} +``` + +#### FreeSurfer + +Dale, A. M., Fischl, B., & Sereno, M. I. (1999). Cortical surface-based analysis: I. Segmentation and surface reconstruction. *NeuroImage*, 9(2), 179-194. + +```bibtex +@article{dale1999cortical, + title={Cortical surface-based analysis: I. Segmentation and surface reconstruction}, + author={Dale, Anders M and Fischl, Bruce and Sereno, Martin I}, + journal={NeuroImage}, + volume={9}, + number={2}, + pages={179--194}, + year={1999}, + publisher={Elsevier} +} +``` + +## Data Availability + +The MEG and MRI data used in this study are available at: +- **OSF Repository**: [https://osf.io/pd4h9/overview](https://osf.io/pd4h9/overview) + +## Software Availability + +The source code for this project is available at: +- **GitHub Repository**: [https://github.com/braindatalab/AD_PAC](https://github.com/braindatalab/AD_PAC) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..a25f399 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,237 @@ +# Contributing to AD_PAC + +Thank you for your interest in contributing to the AD_PAC project! This document provides guidelines for contributing to this repository. + +## How to Contribute + +### Reporting Issues + +If you encounter bugs or have suggestions for improvements: + +1. Check if the issue already exists in the [Issues](https://github.com/braindatalab/AD_PAC/issues) section +2. If not, create a new issue with: + - A clear, descriptive title + - Detailed description of the problem or suggestion + - Steps to reproduce (for bugs) + - Expected vs. actual behavior + - Your environment (MATLAB version, OS, toolbox versions) + - Relevant code snippets or error messages + +### Submitting Changes + +1. **Fork the repository** to your GitHub account + +2. **Clone your fork** locally: + ```bash + git clone https://github.com/YOUR_USERNAME/AD_PAC.git + cd AD_PAC + ``` + +3. **Create a new branch** for your changes: + ```bash + git checkout -b feature/your-feature-name + # or + git checkout -b fix/bug-description + ``` + +4. **Make your changes**: + - Follow the existing code style + - Add comments where necessary + - Update documentation if you change functionality + - Test your changes thoroughly + +5. **Commit your changes** with clear, descriptive messages: + ```bash + git add . + git commit -m "Add feature: brief description" + ``` + +6. **Push to your fork**: + ```bash + git push origin feature/your-feature-name + ``` + +7. **Submit a Pull Request** (PR): + - Go to the original repository on GitHub + - Click "New Pull Request" + - Select your fork and branch + - Provide a clear description of your changes + - Reference any related issues + +## Code Guidelines + +### MATLAB Code Style + +- Use clear, descriptive variable names +- Add comments for complex operations +- Follow MATLAB naming conventions: + - Functions: `camelCase` or `snake_case` + - Variables: `camelCase` or `snake_case` + - Constants: `UPPER_CASE` +- Include function documentation headers: + ```matlab + function output = myFunction(input1, input2) + % MYFUNCTION Brief description + % + % Detailed description of what the function does + % + % Input: + % input1 - Description of first input + % input2 - Description of second input + % + % Output: + % output - Description of output + % + % Example: + % result = myFunction(data, params); + + % Function implementation + end + ``` + +### Documentation + +- Update README files if you add new features +- Document new functions and scripts +- Include usage examples for new functionality +- Update CITATION.md if you add methods that should be cited + +### Testing + +- Test your changes with sample data +- Ensure backward compatibility when possible +- Document any breaking changes clearly +- Verify that existing functionality still works + +## Types of Contributions + +### Bug Fixes + +- Fix incorrect calculations +- Resolve errors or crashes +- Improve error handling +- Fix documentation errors + +### New Features + +- New analysis methods +- Additional visualization options +- Performance improvements +- Support for additional data formats + +### Documentation + +- Improve existing documentation +- Add tutorials or examples +- Clarify confusing sections +- Translate documentation + +### Code Quality + +- Refactor for clarity or efficiency +- Add unit tests +- Improve error messages +- Remove deprecated code + +## Development Setup + +### Prerequisites + +1. MATLAB (R2018b or later recommended) +2. Required toolboxes: + - FieldTrip + - Brainstorm + - FreeSurfer +3. Git for version control + +### Setting Up Development Environment + +1. Clone the repository +2. Install required toolboxes +3. Add toolboxes to MATLAB path +4. Download sample data (if available) +5. Run tests to ensure everything works + +### Testing Your Changes + +Before submitting a PR: + +1. Run the modified scripts with test data +2. Check for errors or warnings +3. Verify outputs are as expected +4. Test edge cases +5. Ensure documentation is accurate + +## Pull Request Process + +1. **PR Title**: Use a clear, descriptive title + - Good: "Fix memory leak in PAC computation" + - Bad: "Fixed bug" + +2. **PR Description**: Include: + - What changes were made + - Why the changes were needed + - How to test the changes + - Any breaking changes + - Related issue numbers (if applicable) + +3. **Review Process**: + - Maintainers will review your PR + - Address any requested changes + - Once approved, your PR will be merged + +4. **Merge**: + - PRs are typically merged via "Squash and merge" + - Your contribution will be acknowledged + +## Code of Conduct + +### Our Pledge + +We are committed to providing a welcoming and inclusive environment for all contributors, regardless of: +- Age, body size, disability +- Ethnicity, gender identity and expression +- Level of experience, education, socio-economic status +- Nationality, personal appearance, race, religion +- Sexual identity and orientation + +### Our Standards + +**Positive behaviors**: +- Using welcoming and inclusive language +- Being respectful of differing viewpoints +- Gracefully accepting constructive criticism +- Focusing on what is best for the community +- Showing empathy towards others + +**Unacceptable behaviors**: +- Trolling, insulting/derogatory comments, personal attacks +- Public or private harassment +- Publishing others' private information without permission +- Other conduct which could reasonably be considered inappropriate + +### Enforcement + +Instances of unacceptable behavior may be reported to the project maintainers. All complaints will be reviewed and investigated promptly and fairly. + +## Questions? + +If you have questions about contributing: +- Open an issue with the "question" label +- Contact the maintainers +- Check existing documentation + +## Recognition + +Contributors will be acknowledged in: +- The repository's contributor list +- Release notes (for significant contributions) +- Academic citations (for methodological contributions) + +## License + +By contributing to AD_PAC, you agree that your contributions will be licensed under the MIT License. + +--- + +Thank you for contributing to AD_PAC! Your efforts help advance research in Alzheimer's Disease and neuroscience. diff --git a/HeadModelling/README.md b/HeadModelling/README.md new file mode 100644 index 0000000..d9fbe5e --- /dev/null +++ b/HeadModelling/README.md @@ -0,0 +1,172 @@ +# Head Modelling Pipeline + +This directory contains scripts for constructing head models and leadfield matrices required for MEG source reconstruction. + +## Overview + +The head modeling pipeline processes anatomical MRI data and MEG sensor information to create forward models that describe how neural sources project to MEG sensors. These models are essential for accurate source localization. + +## Pipeline Steps + +### 1. anatomyTransform.m + +**Purpose**: Corrects MRI transformation matrices and converts MRI data to NIfTI format. + +**Input**: +- Raw MRI data in MAT format from the Dataset folder + +**Output**: +- NIfTI files with corrected transformation matrices +- Prepared for FreeSurfer processing + +**Usage**: +```matlab +anatomyTransform; +``` + +**Note**: This step must be run before FreeSurfer processing. + +### 2. removeElec.m + +**Purpose**: Removes the 'elec' field from MEG files to avoid conflicts in Brainstorm. + +**Input**: +- MEG data files from the Dataset folder + +**Output**: +- Cleaned MEG files compatible with Brainstorm + +**Usage**: +```matlab +removeElec; +``` + +**Why**: Some MEG files may contain electrode information that conflicts with Brainstorm's processing. This step ensures compatibility. + +### 3. buildHeadModels.m + +**Purpose**: Constructs realistic head models and computes leadfield matrices using Brainstorm and OpenMEEG. + +**Requirements**: +- Brainstorm must be installed and in the MATLAB path +- FreeSurfer must have been run on the MRI data +- OpenMEEG must be properly configured in Brainstorm + +**Input**: +- Processed MRI data (from FreeSurfer) +- MEG sensor information + +**Output**: +- Head models (skull, brain surfaces) +- Leadfield matrices for each subject +- Saved in Brainstorm database format + +**Usage**: +```matlab +buildHeadModels; +``` + +**Process**: +1. Imports MRI and MEG data into Brainstorm +2. Generates cortical surface from FreeSurfer output +3. Creates boundary element model (BEM) surfaces +4. Computes forward model using OpenMEEG +5. Generates leadfield matrices + +### 4. processBSFiles.m + +**Purpose**: Post-processes Brainstorm outputs for use in source reconstruction. + +**Input**: +- Brainstorm database with computed head models and leadfields + +**Output**: +- Leadfield matrices in MNI space +- High-resolution cortical surface extrapolation +- Saved in MAT format for easy loading + +**Usage**: +```matlab +processBSFiles; +``` + +**Process**: +1. Loads Brainstorm results +2. Transforms leadfields to MNI standard space +3. Extrapolates to high-resolution cortical mesh +4. Saves leadfields and related metadata + +## Helper Functions + +### eucl.m + +**Purpose**: Computes Euclidean distances between points. + +**Usage**: Called internally by other scripts for spatial calculations. + +## Dependencies + +### External Software + +1. **Brainstorm**: [https://neuroimage.usc.edu/brainstorm/](https://neuroimage.usc.edu/brainstorm/) + - Used for head model construction + - Provides OpenMEEG integration + +2. **FreeSurfer**: [https://surfer.nmr.mgh.harvard.edu/](https://surfer.nmr.mgh.harvard.edu/) + - Must be run independently on MRI data before this pipeline + - Generates cortical surfaces and parcellations + +3. **OpenMEEG**: Integrated within Brainstorm + - Boundary element method for forward modeling + - Provides accurate electromagnetic field computations + +### MATLAB Toolboxes + +- Image Processing Toolbox (for NIfTI handling) +- Signal Processing Toolbox + +## Expected Runtime + +The head modeling pipeline is computationally intensive: +- **anatomyTransform**: ~1-5 minutes per subject +- **removeElec**: <1 minute +- **buildHeadModels**: ~30-120 minutes per subject (depends on mesh resolution) +- **processBSFiles**: ~5-15 minutes per subject + +Total: ~1-2 hours per subject + +## Outputs Structure + +``` +Results/ +└── HeadModel/ + └── [SubjectID]/ + ├── leadfield_mni.mat # Leadfield in MNI space + ├── cortex_highres.mat # High-resolution cortical surface + └── head_model_info.mat # Metadata and transformation matrices +``` + +## Troubleshooting + +### Common Issues + +1. **FreeSurfer not run**: Ensure FreeSurfer's recon-all has completed successfully +2. **Brainstorm not found**: Add Brainstorm to MATLAB path before running +3. **OpenMEEG errors**: Check Brainstorm's OpenMEEG configuration +4. **Memory issues**: Head model computation requires significant RAM (16GB+ recommended) +5. **Transformation errors**: Verify MRI coordinate system and orientation + +### Quality Check + +After running the pipeline: +1. Visually inspect head models in Brainstorm +2. Verify sensor alignment with head surface +3. Check leadfield matrix dimensions +4. Ensure all subjects have complete outputs + +## Notes + +- This pipeline must be run before the source reconstruction steps +- Results are cached and don't need to be recomputed unless MRI/MEG data changes +- Different mesh resolutions can be configured in buildHeadModels.m +- The pipeline assumes CTF MEG system (275 channels); adapt for other systems if needed diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..a32d9d5 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Brain Data Lab + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Pac/README.md b/Pac/README.md new file mode 100644 index 0000000..bf3b29f --- /dev/null +++ b/Pac/README.md @@ -0,0 +1,374 @@ +# Phase-Amplitude Coupling (PAC) Analysis + +This directory contains scripts for computing phase-amplitude coupling between hippocampal and cortical regions. + +## Overview + +The PAC analysis pipeline investigates cross-frequency coupling between the phase of low-frequency oscillations (theta/alpha, 4-12 Hz) in the hippocampus and the amplitude of high-frequency oscillations (low-gamma, 30-55 Hz) in cortical regions. + +## Main Script + +### seedRoitoCortex.m + +**Purpose**: Computes phase-amplitude coupling between hippocampal seed regions and cortical ROIs. + +**Input**: +- Source-reconstructed data from `Results/source/[SubjectID]/source_rec_results.mat` +- Each subject's data contains: + - `source_roi_data`: Time series for ROIs (ROIs × PCs × time) + - `labels`: ROI names + - `regions_cortex`: Cortical region definitions + +**Output**: +- PAC values saved in `Results/SeedtoCortex/[SubjectID]/` +- Dimensions: (L/R hippocampus × ipsi/contra × 34 cortical ROIs × phase/amplitude × freq × pcs) + +**Processing Steps**: + +1. **Load source data**: Read source time series for each subject +2. **Select seed ROIs**: Extract left and right hippocampus activity +3. **Define frequency bands**: + - Phase frequency: 4-12 Hz (theta/alpha) + - Amplitude frequency: 30-55 Hz (low-gamma) +4. **Compute PAC**: For each hippocampus-cortex pair +5. **Save results**: Store PAC matrices for statistical analysis + +**Usage**: +```matlab +seedRoitoCortex; +``` + +## Core PAC Functions + +### er_pac.m + +**Purpose**: Computes phase-amplitude coupling using bispectral analysis (Event-Related PAC). + +**Method**: Bispectrum-based PAC estimation + +**Theory**: +Phase-amplitude coupling can be detected through the bispectrum, which measures nonlinear interactions between frequency components. If the phase of frequency f1 modulates the amplitude of frequency f2, a peak appears in the bispectrum at (f1, f2). + +**Algorithm**: +``` +1. Compute cross-spectrum between two signals X and Y: + CS(f) = FFT(X) * conj(FFT(Y)) + +2. Compute triple product (bispectrum): + B(f1, f2) = + where <> denotes averaging across trials/segments + +3. Normalize by power spectra to get bicoherence: + b(f1, f2) = |B(f1, f2)| / sqrt(<|CS(f1)|²> * <|CS(f2)|²>) + +4. Extract PAC strength at frequencies of interest +``` + +**Input**: +- `data`: Time series data (channels × time × trials) +- `segleng`: Segment length for FFT +- `segshift`: Shift between segments (for overlap) +- `epleng`: Total epoch length +- `freqpairs`: Frequency pairs to analyze + +**Output**: +- `erpac`: PAC values (frequency × frequency × channels) +- `freq`: Frequency vector +- Statistical significance (optional) + +**Parameters**: +```matlab +fs = 600; % Sampling frequency +segleng = 600; % 1-second segments +segshift = 120; % 0.2-second shift (80% overlap) +freq_phase = 4:12; % Phase frequencies (Hz) +freq_amp = 30:55; % Amplitude frequencies (Hz) +``` + +### er_pac_3.m + +**Purpose**: Extended version of er_pac with additional features. + +**Enhancements**: +- Three-way interactions +- Multiple segment lengths +- Advanced normalization options +- Significance testing + +**Usage**: Similar to er_pac.m but with extended capabilities + +### test_pac.m + +**Purpose**: Test and validation script for PAC computation. + +**Functions**: +- Validate PAC methods on synthetic data +- Test parameter sensitivity +- Compare different PAC metrics +- Quality control checks + +**Usage**: For development and validation purposes + +```matlab +% Generate synthetic PAC data +[data_pac, true_coupling] = generate_synthetic_pac(); + +% Test PAC detection +detected_pac = test_pac(data_pac); + +% Compare with ground truth +correlation = corr(true_coupling(:), detected_pac(:)); +``` + +## PAC Analysis Details + +### Hippocampus-to-Cortex Coupling + +The analysis examines several coupling patterns: + +1. **Left hippocampus phase → Left cortex amplitude** (ipsilateral) +2. **Left hippocampus phase → Right cortex amplitude** (contralateral) +3. **Right hippocampus phase → Right cortex amplitude** (ipsilateral) +4. **Right hippocampus phase → Left cortex amplitude** (contralateral) + +5. **Left cortex phase → Left hippocampus amplitude** (reverse coupling) +6. **Right cortex phase → Right hippocampus amplitude** (reverse coupling) + +### Frequency Bands + +**Phase Frequencies (4-12 Hz)**: +- Theta: 4-8 Hz +- Alpha: 8-12 Hz +- Combined theta/alpha: 4-12 Hz + +**Amplitude Frequencies (30-55 Hz)**: +- Low-gamma: 30-55 Hz +- Analysis at 1 Hz resolution + +### Cortical ROIs + +34 cortical regions per hemisphere (Desikan-Killiany parcellation): + +**Frontal**: +- Superior frontal, middle frontal, inferior frontal +- Precentral, postcentral, paracentral + +**Temporal**: +- Superior temporal, middle temporal, inferior temporal +- Fusiform, entorhinal, parahippocampal + +**Parietal**: +- Superior parietal, inferior parietal +- Precuneus, supramarginal + +**Occipital**: +- Lateral occipital, lingual, cuneus, pericalcarine + +**Other**: +- Cingulate regions, insula, etc. + +## Output Format + +### PAC Matrix Structure + +```matlab +% Dimensions: hcpac(L/R, ipsi/contra, cortical_ROI, phase/amp, freq_phase, freq_amp, pcs) +hcpac = zeros(2, 2, 34, 2, 9, 26, 9); + +% Indices: +% dim1: 1=Left hippocampus, 2=Right hippocampus +% dim2: 1=Ipsilateral cortex, 2=Contralateral cortex +% dim3: 1-34 = Cortical ROI index +% dim4: 1=Hippocampus phase, 2=Hippocampus amplitude +% dim5: Phase frequency index (4-12 Hz, 9 points) +% dim6: Amplitude frequency index (30-55 Hz, 26 points) +% dim7: Principal component index (1-9) +``` + +### Saved Files + +``` +Results/ +└── SeedtoCortex/ + └── [SubjectID]/ + ├── hcpac_results.mat # Full PAC matrix + ├── pac_summary.mat # Averaged PAC values + └── quality_metrics.mat # QC measures +``` + +## Expected Runtime + +PAC computation is computationally intensive: + +- **Data loading**: ~30 seconds per subject +- **PAC computation**: ~20-60 minutes per subject + - Depends on: number of segments, frequency resolution, number of ROI pairs +- **Saving results**: ~1 minute per subject + +Total: ~30-90 minutes per subject + +## Configuration + +Key parameters in `seedRoitoCortex.m`: + +```matlab +% Frequency bands +theta_alpha = 4:12; % Phase frequencies (Hz) +low_gamma = 30:55; % Amplitude frequencies (Hz) + +% Segmentation +fs = 600; % Sampling frequency (Hz) +segleng = 600; % Segment length (samples) = 1 second +segshift = 120; % Segment shift (samples) = 0.2 second + % (80% overlap for smooth estimates) + +% Hippocampus ROI indices +LH_id = 10; % Left hippocampus index in ROI list +RH_id = 11; % Right hippocampus index + +% Number of cortical ROIs +n_cortical_roi = 34; % Per hemisphere (Desikan-Killiany) +``` + +## Quality Control + +### Verification Steps + +1. **Check PAC values**: Should be in reasonable range (0 to ~0.1 for bicoherence) +2. **Frequency specificity**: PAC should show clear peaks at expected frequencies +3. **Spatial patterns**: Examine consistency across similar cortical regions +4. **Subject variability**: Compare PAC distributions across subjects + +### Validation Code + +```matlab +% Load PAC results +load('Results/SeedtoCortex/SubjectID/hcpac_results.mat'); + +% Check value ranges +min_pac = min(hcpac(:)); +max_pac = max(hcpac(:)); +mean_pac = mean(hcpac(:)); + +fprintf('PAC range: [%.4f, %.4f], mean: %.4f\n', min_pac, max_pac, mean_pac); + +% Plot frequency-specific PAC +figure; +imagesc(squeeze(mean(hcpac(1, 1, :, 1, :, :), [1,3,7]))); +xlabel('Amplitude Frequency (Hz)'); +ylabel('Phase Frequency (Hz)'); +title('Average Left Hippocampus PAC'); +colorbar; +``` + +## Troubleshooting + +### Common Issues + +1. **NaN values in PAC**: + - Check for zero-variance segments + - Verify sufficient data length + - Ensure proper normalization + +2. **Extremely low PAC values**: + - May indicate genuine lack of coupling + - Check source reconstruction quality + - Verify frequency bands are appropriate + +3. **Memory errors**: + - Process subjects individually + - Reduce frequency resolution + - Use fewer segments with longer averaging + +4. **Very long computation time**: + - Reduce segment overlap (increase segshift) + - Parallelize across subjects + - Optimize FFT computation + +### Optimization Tips + +```matlab +% Reduce frequency resolution +theta_alpha = 4:2:12; % Every 2 Hz instead of 1 Hz +low_gamma = 30:2:55; % Every 2 Hz + +% Reduce overlap +segshift = 300; % 50% overlap instead of 80% + +% Process specific ROIs only +roi_subset = [1, 5, 10]; % Subset of cortical ROIs +``` + +## Statistical Considerations + +### Multiple Comparisons + +PAC analysis involves many comparisons: +- 2 hippocampi × 2 hemispheres × 34 ROIs × 9 phase freq × 26 amp freq +- ≈ 32,000+ tests per subject + +**Correction methods** (applied in Stat/ pipeline): +- Cluster-based permutation testing +- False Discovery Rate (FDR) +- Bonferroni correction (conservative) + +### Effect Size + +Raw PAC values can be small; consider: +- Normalizing by baseline or control condition +- Computing z-scores across subjects +- Using percent change from mean + +## Advanced Usage + +### Custom Frequency Bands + +```matlab +% Define custom frequency pairs +freq_pairs = [4 30; % Theta-low gamma + 6 40; % Theta-gamma + 10 45; % Alpha-gamma + 8 60]; % Alpha-high gamma + +% Compute PAC for specific pairs +for pair = 1:size(freq_pairs, 1) + f_phase = freq_pairs(pair, 1); + f_amp = freq_pairs(pair, 2); + pac(pair) = compute_pac_pair(data, f_phase, f_amp); +end +``` + +### Directionality Analysis + +Test causal direction of PAC: + +```matlab +% Compare hippocampus→cortex vs cortex→hippocampus +pac_hc_to_ctx = hcpac(:, :, :, 1, :, :, :); % HC phase +pac_ctx_to_hc = hcpac(:, :, :, 2, :, :, :); % HC amplitude + +% Statistical test for asymmetry +[h, p] = ttest(pac_hc_to_ctx(:), pac_ctx_to_hc(:)); +``` + +## References + +### PAC Methods + +This implementation is based on: + +Pellegrini, F., Delgado Saa, J., & Haufe, S. (2023). Identifying good practices for detecting inter-frequency coupling in resting-state EEG: A simulation study. *NeuroImage*, 277, 120233. + +Original code: [https://github.com/fpellegrini/PAC](https://github.com/fpellegrini/PAC) + +### Theoretical Background + +- Tort, A. B., et al. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. *Journal of Neurophysiology*, 104(2), 1195-1210. +- Canolty, R. T., & Knight, R. T. (2010). The functional role of cross-frequency coupling. *Trends in Cognitive Sciences*, 14(11), 506-515. + +## Notes + +- PAC analysis is sensitive to signal quality and SNR +- Results should be interpreted in context of specific hypothesis +- Consider running validation analysis on known PAC patterns +- Save intermediate results for debugging and quality control diff --git a/PreProcessing/README.md b/PreProcessing/README.md new file mode 100644 index 0000000..9437415 --- /dev/null +++ b/PreProcessing/README.md @@ -0,0 +1,293 @@ +# Preprocessing Pipeline + +This directory contains scripts for preprocessing raw MEG data, including resampling, artifact rejection, and quality control. + +## Overview + +The preprocessing pipeline prepares raw MEG data for source reconstruction by removing artifacts, interpolating bad channels, and segmenting the continuous recordings into trials. + +## Main Script + +### preprocessing.m + +**Purpose**: Main preprocessing pipeline that coordinates all preprocessing steps. + +**Input**: +- Raw MEG data files (MAT format) from the Dataset folder +- Files should be named: `[SubjectID]_meg_rest_60sec.mat` + +**Output**: +- Preprocessed MEG data saved in `Dataset/Preprocessed/` +- Quality control plots in `Results/Preprocessing/[SubjectID]/` + +**Processing Steps**: +1. Load raw MEG data +2. Resample to 200 Hz +3. Detect and interpolate bad channels +4. Segment into 2-second trials +5. Detect and reject bad trials +6. Save preprocessed data + +**Configuration**: +```matlab +resamplefs = 200; % Target sampling frequency (Hz) +trial_length = 2; % Trial duration (seconds) +amplitude_threshold = 3e-12; % Artifact detection threshold (Tesla) +``` + +**Usage**: +```matlab +preprocessing; +``` + +## Helper Functions + +### interpolate_bad_channels.m + +**Purpose**: Detects and interpolates bad MEG channels using spatial interpolation. + +**Method**: +- Identifies channels with excessive variance or flat signals +- Uses neighboring channels for interpolation +- Based on FieldTrip's channel repair functions + +**Input**: +- FieldTrip data structure with MEG channels + +**Output**: +- Data with bad channels interpolated +- List of interpolated channel labels + +**Algorithm**: +```matlab +% Detect bad channels based on variance +variance_threshold = 3 * std(channel_variances); +bad_channels = channels exceeding threshold + +% Interpolate using neighbors +cfg = []; +cfg.method = 'weighted'; +cfg.badchannel = bad_channel_labels; +cfg.layout = 'CTF275_helmet.mat'; +data_interp = ft_channelrepair(cfg, data); +``` + +### detect_bad_trials.m + +**Purpose**: Identifies and marks trials containing artifacts. + +**Method**: +- Computes peak-to-peak amplitude for each trial and channel +- Flags trials exceeding amplitude threshold +- Uses robust statistical measures to avoid false positives + +**Input**: +- Segmented trial data (FieldTrip structure) +- Amplitude threshold (default: 3e-12 Tesla) + +**Output**: +- Vector of bad trial indices +- Quality metrics (percentage rejected, etc.) + +**Algorithm**: +```matlab +% For each trial +for trial = 1:ntrials + % Compute peak-to-peak amplitude across channels + ptp = max(data{trial}, [], 2) - min(data{trial}, [], 2); + + % Flag if any channel exceeds threshold + if any(ptp > threshold) + bad_trials(trial) = true; + end +end + +% Additional checks for gradient artifacts +% (sudden jumps between consecutive samples) +``` + +### plot_data.m + +**Purpose**: Creates quality control visualizations of MEG data. + +**Plots Generated**: +1. **Topographic maps**: Spatial distribution of signal power +2. **Time series**: Channel-wise or averaged signals over time +3. **Power spectrum**: Frequency content of the data +4. **Trial variance**: Distribution of signal variance across trials + +**Input**: +- FieldTrip data structure +- Layout information for sensor positions +- Output directory for saving figures + +**Output**: +- PNG/PDF files with visualizations +- Saved in subject-specific report directories + +**Usage**: +```matlab +plot_data(report_dir, data, layout, channel_names, plot_name); +``` + +## Dependencies + +### External Software + +1. **FieldTrip Toolbox**: [https://www.fieldtriptoolbox.org/](https://www.fieldtriptoolbox.org/) + - Core preprocessing functions (ft_preprocessing, ft_resampledata, etc.) + - Channel repair and artifact rejection + - Visualization tools + +### MATLAB Toolboxes + +- Signal Processing Toolbox (for filtering and resampling) +- Statistics and Machine Learning Toolbox (for outlier detection) + +## Configuration + +Edit `preprocessing.m` to adjust parameters: + +```matlab +% Resampling +resamplefs = 200; % Target frequency (Hz) + +% Bad channel detection +variance_multiplier = 3; % Standard deviations from mean + +% Trial segmentation +trial_length = 2; % Trial duration (seconds) +trial_overlap = 0; % Overlap between trials (seconds) + +% Artifact rejection +amplitude_threshold = 3e-12; % Peak-to-peak threshold (Tesla) +gradient_threshold = 1e-13; % Maximum allowed gradient +``` + +## Expected Runtime + +Preprocessing is relatively fast: +- **Load data**: <1 minute per subject +- **Resampling**: ~1-2 minutes per subject +- **Bad channel detection**: ~30 seconds per subject +- **Trial segmentation**: ~30 seconds per subject +- **Bad trial detection**: ~1-2 minutes per subject +- **Plotting**: ~1-3 minutes per subject + +Total: ~5-10 minutes per subject (60-second recordings) + +## Output Structure + +``` +Dataset/ +└── Preprocessed/ + └── [SubjectID]/ + └── preprocessed_data.mat # Cleaned, resampled, segmented data + +Results/ +└── Preprocessing/ + └── [SubjectID]/ + ├── raw_after_import_Grad.png # Raw data visualization + ├── after_resampling_Grad.png # After resampling + ├── after_interpolation_Grad.png # After bad channel repair + ├── after_segmentation_Grad.png # After trial segmentation + └── after_rejection_Grad.png # After artifact rejection +``` + +## Quality Control + +### Visual Inspection + +Review the generated plots for each subject: + +1. **Raw data**: Check for obvious artifacts, sensor malfunction +2. **After resampling**: Verify no aliasing artifacts introduced +3. **After interpolation**: Ensure bad channels properly recovered +4. **After segmentation**: Check trial length and overlap +5. **After rejection**: Verify artifact removal without excessive data loss + +### Metrics to Check + +```matlab +% Load preprocessed data +load('Dataset/Preprocessed/SubjectID/preprocessed_data.mat'); + +% Check data quality +nchannels_interpolated = length(bad_channels); +ntrials_rejected = length(bad_trials); +rejection_rate = ntrials_rejected / total_trials * 100; + +% Acceptable ranges: +% - Interpolated channels: 0-10% of total +% - Trial rejection: 5-30% typical +% - If >50% rejected: investigate data quality +``` + +## Troubleshooting + +### Common Issues + +1. **FieldTrip not in path**: Add FieldTrip to MATLAB path and run `ft_defaults` +2. **Layout file not found**: Download CTF275_helmet.mat or specify correct layout +3. **Too many trials rejected**: Lower amplitude threshold or check for systematic artifacts +4. **Memory errors**: Process subjects individually rather than in batch +5. **Resampling errors**: Check original sampling rate and ensure valid target rate + +### Data Quality Issues + +**High rejection rate (>50%)**: +- Check electrode impedances +- Review raw data for systematic artifacts +- Consider manual artifact removal before automated pipeline + +**Poor channel interpolation**: +- May indicate sensor malfunction +- Consider excluding channels if interpolation unsuccessful +- Verify sensor positions in layout file + +**Unexpected frequency content**: +- Check for line noise (50/60 Hz and harmonics) +- Apply notch filter if necessary +- Verify proper shielding during recording + +## Advanced Usage + +### Custom Filtering + +Add bandpass filtering to preprocessing.m: + +```matlab +% After resampling +cfg = []; +cfg.bpfilter = 'yes'; +cfg.bpfreq = [1 95]; % Bandpass 1-95 Hz +cfg.bpfiltord = 4; +data_filtered = ft_preprocessing(cfg, data_resampled); +``` + +### Alternative Artifact Rejection + +For more sophisticated artifact removal: + +```matlab +% Use ICA for artifact removal +cfg = []; +cfg.method = 'runica'; +comp = ft_componentanalysis(cfg, data); + +% Visualize components +ft_databrowser(cfg, comp); + +% Remove artifact components +cfg = []; +cfg.component = [1 3 7]; % Components to remove +data_clean = ft_rejectcomponent(cfg, comp, data); +``` + +## Notes + +- Preprocessing parameters may need adjustment for different datasets +- Visual inspection of quality control plots is essential +- Save preprocessing parameters for reproducibility +- Consider batch processing for large datasets +- The pipeline assumes CTF MEG system; adapt for other systems as needed diff --git a/QUICKSTART.md b/QUICKSTART.md new file mode 100644 index 0000000..04d2ad8 --- /dev/null +++ b/QUICKSTART.md @@ -0,0 +1,332 @@ +# Quick Start Guide + +This guide will help you get started with the AD_PAC pipeline quickly. + +## Prerequisites + +Before starting, ensure you have: +- ✅ MATLAB (R2018b or later) with required toolboxes +- ✅ FieldTrip, Brainstorm, and FreeSurfer installed +- ✅ MEG and MRI data downloaded from [OSF](https://osf.io/pd4h9/overview) +- ✅ FreeSurfer preprocessing completed for MRI data + +See [REQUIREMENTS.md](REQUIREMENTS.md) for detailed installation instructions. + +## Step 1: Clone the Repository + +```bash +git clone https://github.com/braindatalab/AD_PAC.git +cd AD_PAC +``` + +## Step 2: Set Up Data + +1. Create a `Dataset` folder in the repository: + ```bash + mkdir -p Dataset + ``` + +2. Place your data files in the Dataset folder: + ``` + Dataset/ + ├── Sub001_meg_rest_60sec.mat + ├── Sub002_meg_rest_60sec.mat + ├── ... + ├── Sub001_mri.nii + ├── Sub002_mri.nii + └── ... + ``` + +## Step 3: Configure Paths + +1. Open MATLAB +2. Navigate to the AD_PAC directory +3. Edit `main.m` to set the main path: + ```matlab + main_path = '/full/path/to/AD_PAC/'; + ``` + +4. Add required toolboxes to MATLAB path: + ```matlab + % Add FieldTrip + addpath('/path/to/fieldtrip/'); + ft_defaults; + + % Add Brainstorm + addpath('/path/to/brainstorm3/'); + + % Add AD_PAC and subdirectories + addpath(genpath(main_path)); + ``` + +## Step 4: Run FreeSurfer (if not already done) + +FreeSurfer must be run **before** the MATLAB pipeline: + +```bash +# Set up FreeSurfer environment +export FREESURFER_HOME=/path/to/freesurfer +source $FREESURFER_HOME/SetUpFreeSurfer.sh + +# Process each subject's MRI (takes 6-24 hours per subject) +recon-all -s Sub001 -i Dataset/Sub001_mri.nii -all +recon-all -s Sub002 -i Dataset/Sub002_mri.nii -all +# ... repeat for all subjects +``` + +## Step 5: Run the Pipeline + +### Option A: Run the Complete Pipeline + +In MATLAB, run the entire pipeline: + +```matlab +% Open main.m and run the entire script +main +``` + +This will execute all pipeline stages sequentially. + +### Option B: Run Individual Stages + +For more control, run each stage separately: + +#### 1. Head Modeling (do this first) + +```matlab +% Correct MRI transformations +anatomyTransform; + +% Remove electrode fields +removeElec; + +% Build head models (takes ~1-2 hours per subject) +buildHeadModels; + +% Process Brainstorm files +processBSFiles; +``` + +#### 2. Preprocessing (~5-10 minutes per subject) + +```matlab +preprocessing; +``` + +Check the quality control plots in `Results/Preprocessing/[SubjectID]/` + +#### 3. Source Reconstruction (~15-40 minutes per subject) + +```matlab +sourceReconstruction_p; +``` + +#### 4. PAC Analysis (~30-90 minutes per subject) + +```matlab +seedRoitoCortex; +``` + +#### 5. Statistical Analysis + +```matlab +% Prepare data table +prepareTable; + +% Analyze power spectra +powerAnalysis; + +% Cluster-based permutation testing +clusterBasedModel; + +% Linear mixed-effects modeling +linearMixedEffectsModel; +``` + +## Step 6: View Results + +Results are saved in the `Results/` directory: + +``` +Results/ +├── Preprocessing/ # QC plots for each subject +├── source/ # Source reconstruction results +├── SeedtoCortex/ # PAC analysis results +└── Stat/ # Statistical analysis outputs + ├── ClusterTest/ # Cluster-based permutation results + ├── LME/ # Mixed-effects model results + └── Power/ # Power analysis results +``` + +## Expected Timeline + +For a typical dataset (20-30 subjects): + +| Stage | Time per Subject | Total Time (25 subjects) | +|-------|-----------------|--------------------------| +| FreeSurfer | 6-24 hours | 150-600 hours (run in parallel) | +| Head Modeling | 1-2 hours | 25-50 hours | +| Preprocessing | 5-10 minutes | 2-4 hours | +| Source Reconstruction | 15-40 minutes | 6-17 hours | +| PAC Analysis | 30-90 minutes | 12-38 hours | +| Statistical Analysis | 1-3 hours | 1-3 hours (group-level) | + +**Total**: ~1-2 weeks of computation time (can be parallelized) + +## Tips for Faster Processing + +1. **Parallel FreeSurfer**: Run multiple subjects simultaneously + ```bash + # Example: Process 4 subjects in parallel + for subj in Sub001 Sub002 Sub003 Sub004; do + recon-all -s $subj -i Dataset/${subj}_mri.nii -all & + done + wait + ``` + +2. **MATLAB Parallel Processing**: Use parallel computing for preprocessing + ```matlab + parpool('local', 8); % Use 8 cores + parfor i = 1:length(subjects) + % Process subject i + end + ``` + +3. **Process in Batches**: Don't wait for all subjects; process in batches + - Process first batch through entire pipeline + - Start next batch while analyzing first batch + +## Common Issues and Solutions + +### Issue 1: FieldTrip Function Not Found +``` +Error: Undefined function 'ft_preprocessing' +``` +**Solution**: Add FieldTrip to path and run `ft_defaults` + +### Issue 2: Out of Memory Error +``` +Error: Out of memory +``` +**Solution**: +- Close other applications +- Process fewer subjects simultaneously +- Increase MATLAB Java heap size + +### Issue 3: Brainstorm Database Not Found +``` +Error: Cannot find Brainstorm database +``` +**Solution**: +- Check that `buildHeadModels.m` completed successfully +- Verify Brainstorm configuration + +### Issue 4: FreeSurfer Outputs Not Found +``` +Error: Cannot find FreeSurfer recon-all output +``` +**Solution**: +- Ensure FreeSurfer preprocessing completed without errors +- Check FreeSurfer SUBJECTS_DIR environment variable +- Verify subject directories exist + +## Verifying Results + +### After Preprocessing: +```matlab +% Load preprocessed data +load('Dataset/Preprocessed/Sub001/preprocessed_data.mat'); + +% Check dimensions +disp(size(data_preprocessed.trial{1})); % Should be [channels x time] + +% Visual inspection +plot(data_preprocessed.time{1}, data_preprocessed.trial{1}(1,:)); +xlabel('Time (s)'); ylabel('Amplitude (T)'); title('Channel 1'); +``` + +### After Source Reconstruction: +```matlab +% Load source results +load('Results/source/Sub001/source_rec_results.mat'); + +% Check ROI data +disp(size(source_roi_data)); % [n_rois x n_pcs x n_timepoints] +disp(labels); % ROI labels + +% Plot first ROI +figure; plot(squeeze(source_roi_data(1,1,:))); +title(['ROI: ' labels{1}]); +``` + +### After PAC Analysis: +```matlab +% Load PAC results +load('Results/SeedtoCortex/Sub001/hcpac_results.mat'); + +% Check PAC matrix +disp(size(hcpac)); % [2, 2, 34, 2, 9, 26, 9] + +% Plot average PAC for left hippocampus +pac_avg = squeeze(mean(hcpac(1,1,:,1,:,:,:), [3,7])); +imagesc(pac_avg); +xlabel('Amplitude Frequency (Hz)'); +ylabel('Phase Frequency (Hz)'); +title('Left Hippocampus PAC'); +colorbar; +``` + +## Next Steps + +After running the pipeline: + +1. **Review Quality Control**: Check all QC plots in `Results/Preprocessing/` +2. **Examine Statistical Results**: Review significant clusters and effects +3. **Visualize Findings**: Use plotting scripts to create publication figures +4. **Validate Results**: Compare with literature and expectations +5. **Customize Analysis**: Modify parameters for your specific research questions + +## Getting Help + +- **Documentation**: See README.md files in each subdirectory +- **Issues**: Open an issue on [GitHub](https://github.com/braindatalab/AD_PAC/issues) +- **Citation**: See [CITATION.md](CITATION.md) for references + +## Example: Minimal Working Example + +Here's a minimal example to process a single subject: + +```matlab +%% Setup +main_path = '/path/to/AD_PAC/'; +addpath(genpath(main_path)); +addpath('/path/to/fieldtrip/'); ft_defaults; +addpath('/path/to/brainstorm3/'); + +%% Process one subject +subject_id = 'Sub001'; + +% 1. Head model (assumes FreeSurfer completed) +% Run buildHeadModels.m for this subject + +% 2. Preprocessing +% Run preprocessing.m for this subject + +% 3. Source reconstruction +% Run sourceReconstruction_p.m for this subject + +% 4. PAC analysis +% Run seedRoitoCortex.m for this subject + +%% Check results +load(['Results/SeedtoCortex/' subject_id '/hcpac_results.mat']); +fprintf('PAC analysis complete for %s\n', subject_id); +fprintf('PAC values range: [%.4f, %.4f]\n', min(hcpac(:)), max(hcpac(:))); +``` + +--- + +For detailed documentation, see: +- [README.md](README.md) - Full documentation +- [REQUIREMENTS.md](REQUIREMENTS.md) - Detailed installation guide +- Module-specific READMEs in each subdirectory diff --git a/README.md b/README.md new file mode 100644 index 0000000..1f9e326 --- /dev/null +++ b/README.md @@ -0,0 +1,227 @@ +# AD_PAC: Phase-Amplitude Coupling Analysis in Alzheimer's Disease + +This repository contains MATLAB code for analyzing phase-amplitude coupling (PAC) in MEG data from patients with Alzheimer's Disease (AD) and healthy controls. + +## Citation + +If you use this code in your research, please cite our preprint: + +**Preprint:** [https://www.medrxiv.org/content/10.64898/2026.02.06.26345635v1](https://www.medrxiv.org/content/10.64898/2026.02.06.26345635v1) + +## Overview + +This pipeline processes MEG (Magnetoencephalography) data to investigate phase-amplitude coupling between hippocampal and cortical regions in Alzheimer's Disease. The analysis focuses on theta/alpha phase coupling with low-gamma amplitude, examining differences between AD patients and healthy controls. + +## Key Features + +- **Head Modeling**: Automated head model construction using Brainstorm and OpenMEEG +- **MEG Preprocessing**: Bad channel interpolation, resampling, and trial rejection +- **Source Reconstruction**: LCMV beamformer-based source localization +- **PAC Analysis**: Cross-frequency coupling between hippocampus and cortical ROIs +- **Statistical Analysis**: Cluster-based permutation tests and linear mixed-effects models + +## Requirements + +### Software Dependencies + +1. **MATLAB** (tested with recent versions) +2. **FieldTrip Toolbox** - For anatomy transformation and preprocessing + - Download: [https://www.fieldtriptoolbox.org/](https://www.fieldtriptoolbox.org/) +3. **Brainstorm Toolbox** - For head model construction + - Download: [https://neuroimage.usc.edu/brainstorm/](https://neuroimage.usc.edu/brainstorm/) +4. **FreeSurfer** - For cortical surface extraction (must be run on MRI files before the pipeline) + - Download: [https://surfer.nmr.mgh.harvard.edu/](https://surfer.nmr.mgh.harvard.edu/) + +### Data + +MEG and MRI data can be downloaded from: +- **OSF Repository**: [https://osf.io/pd4h9/overview](https://osf.io/pd4h9/overview) + +Data should be placed in a `Dataset` folder under the main project path. + +### PAC Code Attribution + +The PAC estimation code is based on the implementation by Franziska Pellegrini and Stefan Haufe: +- **GitHub**: [https://github.com/fpellegrini/PAC](https://github.com/fpellegrini/PAC) + +## Installation + +1. Clone this repository: + ```bash + git clone https://github.com/braindatalab/AD_PAC.git + cd AD_PAC + ``` + +2. Download and install required toolboxes (FieldTrip, Brainstorm, FreeSurfer) + +3. Download the MEG and MRI data from OSF and place in `Dataset/` folder + +4. Add required toolboxes to your MATLAB path in the `main.m` script + +## Usage + +The complete analysis pipeline is orchestrated through the `main.m` script. The pipeline consists of five main stages: + +### 1. Head Modeling Pipeline + +```matlab +anatomyTransform; % Correct MRI transformation matrices and save as NIfTI files +removeElec; % Remove 'elec' field from MEG files +buildHeadModels; % Build head models and leadfields using Brainstorm and OpenMEEG +processBSFiles; % Convert to MNI, extrapolate to high-res cortex, save leadfields +``` + +### 2. Preprocessing Pipeline + +```matlab +preprocessing; % Resample, interpolate bad channels, segment trials, reject bad trials +``` + +### 3. Source Reconstruction Pipeline + +```matlab +sourceReconstruction_p; % LCMV beamformer source reconstruction +``` + +### 4. PAC Estimation Pipeline + +```matlab +seedRoitoCortex; % Compute cross-site PAC between hippocampus and cortical ROIs + % (theta/alpha - low-gamma coupling) +``` + +### 5. Statistical Analysis Pipeline + +```matlab +prepareTable; % Prepare table with clinical/demographic data and PAC values +powerAnalysis; % Analyze PSD and region-wise power differences +clusterBasedModel; % Cluster-based statistical analysis on PAC differences +linearMixedEffectsModel; % Linear mixed-effects modeling +``` + +### Running the Complete Pipeline + +To run the full analysis: + +1. Open MATLAB and navigate to the repository directory +2. Edit `main.m` to set your `main_path` variable +3. Ensure all required toolboxes are in your MATLAB path +4. Run individual sections or the complete pipeline: + +```matlab +% Set the main path +main_path = './'; + +% Add all subdirectories to path +addpath(genpath(main_path)); + +% Run individual pipeline stages as needed +% (See main.m for complete pipeline) +``` + +## Directory Structure + +``` +AD_PAC/ +├── main.m # Main pipeline orchestration script +├── Dataset/ # Data directory (not included, download from OSF) +├── HeadModelling/ # Head model construction scripts +│ ├── anatomyTransform.m +│ ├── buildHeadModels.m +│ ├── processBSFiles.m +│ └── removeElec.m +├── PreProcessing/ # MEG data preprocessing scripts +│ ├── preprocessing.m +│ ├── interpolate_bad_channels.m +│ ├── detect_bad_trials.m +│ └── plot_data.m +├── SourceReconstruction/ # Source localization scripts +│ ├── sourceReconstruction_p.m +│ ├── lcmv_meg.m +│ └── ... +├── Pac/ # Phase-amplitude coupling analysis +│ ├── seedRoitoCortex.m +│ ├── er_pac.m +│ └── ... +└── Stat/ # Statistical analysis scripts + ├── prepareTable.m + ├── powerAnalysis.m + ├── clusterBasedModel.m + └── linearMixedEffectsModel.m +``` + +## Pipeline Details + +### Head Modeling + +The head modeling pipeline creates forward models for MEG source localization: +- Transforms anatomical MRI data to standard space +- Constructs realistic head models using OpenMEEG +- Generates leadfield matrices for source reconstruction +- Outputs are saved for each subject and used in subsequent analysis + +### Preprocessing + +MEG data preprocessing includes: +- Resampling to 200 Hz (configurable) +- Bad channel detection and interpolation +- Trial segmentation (2-second trials) +- Artifact rejection based on amplitude thresholds +- Quality control plots are generated for each step + +### Source Reconstruction + +Source-level analysis using LCMV beamformer: +- Computes spatial filters for each cortical location +- Projects sensor-level MEG data to source space +- Focuses on regions of interest (ROIs) including hippocampus and cortical areas +- Produces time series for each ROI + +### PAC Analysis + +Phase-amplitude coupling estimation: +- Computes coupling between hippocampal theta/alpha phase (4-12 Hz) +- And cortical low-gamma amplitude (30-55 Hz) +- Analyzes left/right hippocampus to ipsilateral/contralateral cortex +- Uses bispectral analysis methods +- Results saved per subject for statistical analysis + +### Statistical Analysis + +Group-level statistics: +- Prepares data tables with demographics and PAC metrics +- Power spectral density analysis +- Cluster-based permutation testing for PAC differences +- Linear mixed-effects models accounting for repeated measures +- Visualization of significant results and correlations + +## Output + +Results are organized in the `Results/` directory: +- `Results/Preprocessing/`: Preprocessing quality control plots +- `Results/source/`: Source reconstruction results per subject +- `Results/SeedtoCortex/`: PAC analysis results per subject +- Statistical analysis outputs including plots and model results + +## Troubleshooting + +### Common Issues + +1. **Missing toolboxes**: Ensure FieldTrip, Brainstorm, and FreeSurfer are properly installed and in your MATLAB path +2. **Data not found**: Verify that the Dataset folder exists and contains the downloaded MEG/MRI data +3. **Memory issues**: Large datasets may require substantial RAM; consider processing subjects in batches +4. **FreeSurfer preprocessing**: MRI files must be processed with FreeSurfer before running the head model pipeline + +## Support + +For questions or issues, please open an issue on the GitHub repository. + +## License + +See [LICENSE](LICENSE) file for details. + +## Acknowledgments + +- PAC analysis code adapted from Franziska Pellegrini and Stefan Haufe +- FieldTrip, Brainstorm, and FreeSurfer development teams +- Data contributors and participants diff --git a/REQUIREMENTS.md b/REQUIREMENTS.md new file mode 100644 index 0000000..1ed89ca --- /dev/null +++ b/REQUIREMENTS.md @@ -0,0 +1,384 @@ +# Requirements and Dependencies + +This document details the software requirements and dependencies needed to run the AD_PAC pipeline. + +## System Requirements + +### Hardware Requirements + +**Minimum**: +- CPU: Multi-core processor (4+ cores recommended) +- RAM: 16 GB +- Storage: 50 GB free space (for data and results) + +**Recommended**: +- CPU: 8+ cores for parallel processing +- RAM: 32 GB or more +- Storage: 100+ GB SSD for faster I/O +- GPU: Optional, for accelerated computations in some toolboxes + +### Operating System + +The pipeline has been tested on: +- **Linux**: Ubuntu 18.04+, CentOS 7+, or similar +- **macOS**: 10.14 (Mojave) or later +- **Windows**: Windows 10 or later (with some limitations for FreeSurfer) + +Note: FreeSurfer runs natively on Linux and macOS. Windows users may need to use WSL (Windows Subsystem for Linux) or a virtual machine. + +## Software Dependencies + +### 1. MATLAB + +**Version**: R2018b or later (R2020a+ recommended) + +**Required Toolboxes**: +- Signal Processing Toolbox +- Statistics and Machine Learning Toolbox +- Image Processing Toolbox (for NIfTI handling) +- Parallel Computing Toolbox (optional, for speedup) + +**Installation**: +- Download from [MathWorks](https://www.mathworks.com/products/matlab.html) +- Academic licenses typically include required toolboxes +- Verify installation: `ver` command in MATLAB + +### 2. FieldTrip Toolbox + +**Version**: Latest stable release (tested with 20240110 and later) + +**Purpose**: +- MEG/EEG preprocessing +- Data visualization +- Channel repair and artifact rejection + +**Installation**: +1. Download from [FieldTrip website](https://www.fieldtriptoolbox.org/download/) + ```bash + git clone https://github.com/fieldtrip/fieldtrip.git + ``` + +2. Add to MATLAB path in your script: + ```matlab + addpath('/path/to/fieldtrip/'); + ft_defaults; + ``` + +**Configuration**: +- No special configuration required +- Ensure all subdirectories are in the path +- Run `ft_defaults` once per MATLAB session + +**Documentation**: [https://www.fieldtriptoolbox.org/](https://www.fieldtriptoolbox.org/) + +### 3. Brainstorm Toolbox + +**Version**: Latest stable release (tested with version 3.230101 and later) + +**Purpose**: +- Head model construction +- Forward modeling with OpenMEEG +- MEG/MRI coregistration + +**Installation**: +1. Download from [Brainstorm website](https://neuroimage.usc.edu/brainstorm/Installation) + ```bash + # Or download directly in MATLAB + ``` + +2. Add to MATLAB path: + ```matlab + addpath('/path/to/brainstorm3/'); + ``` + +3. Start Brainstorm GUI once to complete setup: + ```matlab + brainstorm; + ``` + +**Configuration**: +- Set up OpenMEEG during first launch +- Configure paths in Brainstorm preferences +- Test with sample dataset + +**OpenMEEG**: +- Integrated within Brainstorm +- No separate installation needed +- Used for boundary element method (BEM) forward modeling + +**Documentation**: [https://neuroimage.usc.edu/brainstorm/](https://neuroimage.usc.edu/brainstorm/) + +### 4. FreeSurfer + +**Version**: 6.0 or later (7.0+ recommended) + +**Purpose**: +- Cortical surface extraction +- Brain segmentation +- Surface parcellation (Desikan-Killiany atlas) + +**Installation**: + +**Linux/macOS**: +1. Download from [FreeSurfer website](https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall) +2. Extract and set environment variables: + ```bash + export FREESURFER_HOME=/path/to/freesurfer + source $FREESURFER_HOME/SetUpFreeSurfer.sh + ``` +3. Add to your `.bashrc` or `.bash_profile` for persistent setup + +**Windows**: +- Use WSL (Windows Subsystem for Linux) +- Or use a Linux virtual machine +- Native Windows version has limitations + +**License**: +- Free for academic use +- Register for license at FreeSurfer website +- Place license.txt in $FREESURFER_HOME + +**Usage**: +FreeSurfer processing is a prerequisite, run before the MATLAB pipeline: +```bash +# Process anatomical MRI +recon-all -s SubjectID -i /path/to/MRI.nii -all + +# Typical processing time: 6-24 hours per subject +``` + +**Documentation**: [https://surfer.nmr.mgh.harvard.edu/](https://surfer.nmr.mgh.harvard.edu/) + +## Data Requirements + +### Input Data Format + +**MEG Data**: +- Format: MAT files (MATLAB) +- Structure: Should contain MEG data in FieldTrip-compatible format +- Naming: `[SubjectID]_meg_rest_60sec.mat` +- Channels: CTF MEG system (275 channels) +- Sampling rate: 600 Hz or similar + +**MRI Data**: +- Format: NIfTI (.nii, .nii.gz) or DICOM +- Type: T1-weighted anatomical MRI +- Resolution: 1mm isotropic (typical) +- Must be processed with FreeSurfer before pipeline + +**Clinical Data**: +- Format: CSV or Excel +- Required fields: + - Subject ID + - Group (AD / Control) + - Age + - Sex + - MMSE score (or other cognitive measures) +- Optional: Education, disease duration, medications, etc. + +### Data Organization + +``` +AD_PAC/ +├── Dataset/ +│ ├── Sub001_meg_rest_60sec.mat +│ ├── Sub002_meg_rest_60sec.mat +│ ├── ... +│ ├── Sub001_mri.nii +│ ├── Sub002_mri.nii +│ ├── ... +│ └── clinical_data.csv +``` + +### Data Availability + +Sample data can be downloaded from: +- **OSF Repository**: [https://osf.io/pd4h9/overview](https://osf.io/pd4h9/overview) + +## Installation Checklist + +- [ ] Install MATLAB (R2018b+) + - [ ] Verify required toolboxes: `ver` +- [ ] Install FieldTrip + - [ ] Add to MATLAB path + - [ ] Test: `ft_defaults` +- [ ] Install Brainstorm + - [ ] Add to MATLAB path + - [ ] Configure OpenMEEG + - [ ] Test with sample data +- [ ] Install FreeSurfer + - [ ] Set environment variables + - [ ] Obtain and install license + - [ ] Test: `recon-all -version` +- [ ] Download MEG/MRI data + - [ ] Place in `Dataset/` folder + - [ ] Verify file formats and naming +- [ ] Clone AD_PAC repository + - [ ] Set `main_path` in `main.m` + - [ ] Add repository to MATLAB path + +## Verification + +### Test MATLAB Installation + +```matlab +% Check MATLAB version +version + +% Check installed toolboxes +ver + +% Test basic operations +x = randn(100, 100); +y = fft2(x); +disp('MATLAB working correctly'); +``` + +### Test FieldTrip + +```matlab +% Add FieldTrip to path +addpath('/path/to/fieldtrip/'); +ft_defaults; + +% Check version +ft_version + +% Should display version info without errors +``` + +### Test Brainstorm + +```matlab +% Add Brainstorm to path +addpath('/path/to/brainstorm3/'); + +% Start Brainstorm (GUI will open) +brainstorm; + +% Should open without errors +``` + +### Test FreeSurfer + +```bash +# Check version +recon-all -version + +# Check license +cat $FREESURFER_HOME/license.txt + +# Test on sample data (optional) +recon-all -s bert -all +``` + +## Troubleshooting + +### Common Issues + +**MATLAB Toolbox Missing**: +``` +Error: Undefined function 'fft' for input arguments of type 'double'. +Solution: Install Signal Processing Toolbox +``` + +**FieldTrip Not Found**: +``` +Error: Undefined function 'ft_preprocessing'. +Solution: +1. Add FieldTrip to path: addpath('/path/to/fieldtrip/') +2. Run: ft_defaults +``` + +**Brainstorm OpenMEEG Error**: +``` +Error: OpenMEEG not configured +Solution: +1. Open Brainstorm GUI +2. Go to File > Edit preferences +3. Configure OpenMEEG path +``` + +**FreeSurfer License Missing**: +``` +Error: No valid license found +Solution: +1. Register at FreeSurfer website +2. Download license.txt +3. Place in $FREESURFER_HOME/ +``` + +### Getting Help + +- **MATLAB**: [MathWorks Support](https://www.mathworks.com/support.html) +- **FieldTrip**: [Discussion List](https://mailman.science.ru.nl/mailman/listinfo/fieldtrip) +- **Brainstorm**: [Forum](https://neuroimage.usc.edu/forums/) +- **FreeSurfer**: [Mailing List](https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSupport) + +## Version Compatibility + +| Software | Minimum Version | Recommended | Tested With | +|----------|----------------|-------------|-------------| +| MATLAB | R2018b | R2020a+ | R2020b, R2021a, R2022b | +| FieldTrip | 20200101 | Latest | 20240110 | +| Brainstorm | 3.200101 | Latest | 3.230101 | +| FreeSurfer | 6.0 | 7.0+ | 7.1, 7.2 | + +## Performance Optimization + +### MATLAB Optimization + +```matlab +% Use parallel processing +parpool('local', 8); % Use 8 cores + +% Increase Java heap size +java.lang.Runtime.getRuntime.maxMemory / 1024^3 % Check current +% Edit matlab.prf to increase if needed +``` + +### System Optimization + +- **Storage**: Use SSD for data and results +- **Memory**: Close unnecessary applications +- **CPU**: Disable power saving, enable high-performance mode + +## Updates and Maintenance + +### Keeping Software Updated + +**MATLAB**: +- Update through MATLAB interface or website +- Check for updates: Help > Check for Updates + +**FieldTrip**: +```bash +cd /path/to/fieldtrip/ +git pull origin master +``` + +**Brainstorm**: +- Update through Brainstorm GUI: Help > Update Brainstorm + +**FreeSurfer**: +- Download new version and reinstall +- Or use package manager if available + +### Checking for Breaking Changes + +Before updating: +1. Review release notes +2. Test with sample data +3. Keep backup of working versions + +## Additional Resources + +- [MATLAB Documentation](https://www.mathworks.com/help/matlab/) +- [FieldTrip Tutorial](https://www.fieldtriptoolbox.org/tutorial/) +- [Brainstorm Tutorials](https://neuroimage.usc.edu/brainstorm/Tutorials) +- [FreeSurfer Tutorial](https://surfer.nmr.mgh.harvard.edu/fswiki/Tutorials) +- [MEG Analysis Best Practices](https://www.sciencedirect.com/science/article/pii/S1053811918306827) + +## Contact + +For issues specific to this pipeline, please open an issue on the [GitHub repository](https://github.com/braindatalab/AD_PAC/issues). diff --git a/SourceReconstruction/README.md b/SourceReconstruction/README.md new file mode 100644 index 0000000..65360b6 --- /dev/null +++ b/SourceReconstruction/README.md @@ -0,0 +1,386 @@ +# Source Reconstruction Pipeline + +This directory contains scripts for reconstructing neural sources from sensor-level MEG data using beamforming techniques. + +## Overview + +The source reconstruction pipeline transforms preprocessed MEG sensor data into estimated neural activity on the cortical surface using the Linearly Constrained Minimum Variance (LCMV) beamformer. + +## Main Script + +### sourceReconstruction_p.m + +**Purpose**: Orchestrates the complete source reconstruction pipeline for all subjects. + +**Input**: +- Preprocessed MEG data from `Dataset/Preprocessed/` +- Head models and leadfields from `HeadModelling/` pipeline +- ROI definitions and cortical parcellations + +**Output**: +- Source-space time series for cortical ROIs +- Saved in `Results/source/[SubjectID]/source_rec_results.mat` + +**Processing Steps**: +1. Load preprocessed data and head model +2. Compute data covariance matrix +3. Calculate LCMV beamformer spatial filters +4. Project sensor data to source space +5. Extract time series for predefined ROIs +6. Save results for subsequent PAC analysis + +**Usage**: +```matlab +sourceReconstruction_p; +``` + +## Core Functions + +### lcmv_meg.m + +**Purpose**: Implements the LCMV beamformer algorithm for MEG source reconstruction. + +**Method**: Linearly Constrained Minimum Variance (LCMV) beamformer + +**Algorithm**: +``` +For each source location: +1. Extract leadfield vector L for that location +2. Compute spatial filter: W = (L' * C^-1 * L)^-1 * L' * C^-1 + where C is the data covariance matrix +3. Apply filter to sensor data: S = W * M + where M is the sensor data, S is the source activity +``` + +**Input**: +- Leadfield matrix (n_sources × n_sensors) +- MEG sensor data (n_sensors × n_timepoints) +- Data covariance matrix (n_sensors × n_sensors) + +**Output**: +- Source time series (n_sources × n_timepoints) +- Spatial filter weights (n_sources × n_sensors) + +**Parameters**: +```matlab +cfg.lambda = 0.05; % Regularization parameter (5%) +cfg.fixedori = 'yes'; % Use fixed dipole orientation +cfg.normalize = 'yes'; % Normalize spatial filters +``` + +**Key Features**: +- Regularized covariance inversion for numerical stability +- Handles rank-deficient covariance matrices +- Supports fixed or free dipole orientation +- Efficient vectorized computation + +### shrinkage.m + +**Purpose**: Performs covariance shrinkage regularization to stabilize matrix inversion. + +**Method**: Ledoit-Wolf shrinkage estimator + +**Why**: MEG covariance matrices are often ill-conditioned due to: +- Limited number of trials +- Correlated sensor signals +- Noise and artifacts + +**Algorithm**: +``` +C_shrunk = (1 - λ) * C_sample + λ * C_prior + +where: +- C_sample: Sample covariance matrix +- C_prior: Prior covariance (often diagonal) +- λ: Shrinkage parameter (automatically estimated) +``` + +**Input**: +- Sample covariance matrix +- Optional: Prior covariance structure + +**Output**: +- Regularized covariance matrix +- Optimal shrinkage parameter + +### fp_get_lf.m + +**Purpose**: Loads and prepares leadfield matrices for source reconstruction. + +**Input**: +- Subject ID +- Path to leadfield files + +**Output**: +- Leadfield matrix in appropriate format +- Source locations and orientations +- ROI definitions + +**Processing**: +1. Load leadfield from head model +2. Select channels matching MEG data +3. Apply any necessary transformations +4. Return in format compatible with beamformer + +### data2cs_event.m + +**Purpose**: Converts MEG trial data to cross-spectral density matrices. + +**Input**: +- Epoched MEG data (trials × channels × time) + +**Output**: +- Cross-spectral density matrix (channels × channels × frequency) + +**Method**: +- Fourier transform each trial +- Compute cross-spectrum: CS(f) = X(f) * X(f)' +- Average across trials + +**Usage**: For frequency-specific source reconstruction + +### data2spwctrgc.m + +**Purpose**: Computes spectral connectivity measures in source space. + +**Input**: +- Source time series +- Frequency bands of interest + +**Output**: +- Spectral power +- Coherence matrices +- Granger causality (optional) + +**Applications**: +- Frequency-specific source analysis +- Connectivity analysis between ROIs + +### cs2psd.m + +**Purpose**: Converts cross-spectral density to power spectral density. + +**Input**: +- Cross-spectral density matrix + +**Output**: +- Power spectral density for each channel/source + +**Method**: Extract diagonal elements (auto-spectra) from cross-spectral matrix + +## Helper Functions + +### colormap_interpol.m + +**Purpose**: Creates custom colormaps for visualization. + +**Usage**: Generates smooth color gradients for plotting brain activity + +## Dependencies + +### External Software + +1. **FieldTrip Toolbox**: [https://www.fieldtriptoolbox.org/](https://www.fieldtriptoolbox.org/) + - Source reconstruction functions + - Covariance estimation + +### MATLAB Toolboxes + +- Signal Processing Toolbox +- Statistics and Machine Learning Toolbox + +## Configuration + +Key parameters in `sourceReconstruction_p.m`: + +```matlab +% Beamformer settings +cfg.method = 'lcmv'; +cfg.lambda = 0.05; % Regularization (5%) +cfg.fixedori = 'yes'; % Fixed dipole orientation + +% Frequency bands +freq_bands = [1 4; % Delta + 4 8; % Theta + 8 12; % Alpha + 13 30; % Beta + 30 55]; % Low-Gamma + +% ROI selection +roi_atlas = 'Desikan-Killiany'; % Cortical parcellation +include_hippocampus = true; % Include subcortical ROIs +``` + +## Expected Runtime + +Source reconstruction is moderately intensive: +- **Load data and head model**: ~1 minute per subject +- **Compute covariance**: ~2-5 minutes per subject +- **Beamformer computation**: ~10-30 minutes per subject (depends on source space resolution) +- **ROI extraction**: ~1-2 minutes per subject + +Total: ~15-40 minutes per subject + +## Output Structure + +``` +Results/ +└── source/ + └── [SubjectID]/ + └── source_rec_results.mat + ├── source_roi_data # Time series for each ROI (ROIs × PCs × time) + ├── labels # ROI labels + ├── regions_cortex # Cortical region definitions + ├── spatial_filters # Beamformer weights + └── source_locations # 3D coordinates of sources +``` + +### Data Format + +**source_roi_data**: +- Dimensions: (n_rois × n_pcs × n_timepoints) +- n_rois: Number of cortical and subcortical ROIs +- n_pcs: Principal components per ROI (typically 3) +- n_timepoints: Time samples at 600 Hz sampling rate + +**labels**: Cell array of ROI names (e.g., 'L_hippocampus', 'R_superiorfrontal') + +**regions_cortex**: Structure with cortical parcellation information + +## ROI Extraction + +The pipeline extracts source activity for predefined regions: + +### Cortical ROIs (Desikan-Killiany Atlas) +34 regions per hemisphere, including: +- Frontal: Superior frontal, middle frontal, precentral, etc. +- Temporal: Superior temporal, middle temporal, fusiform, etc. +- Parietal: Superior parietal, inferior parietal, precuneus, etc. +- Occipital: Lateral occipital, cuneus, lingual, etc. + +### Subcortical ROIs +- Left hippocampus +- Right hippocampus + +### Dimensionality Reduction + +For each ROI: +1. Extract all source locations within the ROI +2. Compute principal component analysis (PCA) +3. Keep top 3 principal components (captures most variance) +4. Save PC time series for further analysis + +## Quality Control + +### Visual Inspection + +Check source reconstruction quality: + +```matlab +% Load results +load('Results/source/SubjectID/source_rec_results.mat'); + +% Plot source activity +figure; +for roi = 1:size(source_roi_data, 1) + subplot(7, 5, roi); + plot(squeeze(source_roi_data(roi, 1, :))); % First PC + title(labels{roi}); +end +``` + +### Validation Metrics + +1. **Spatial filter quality**: Check for dipolar patterns +2. **SNR**: Signal-to-noise ratio in source space +3. **Localization error**: Validate on known sources (if available) +4. **ROI coverage**: Ensure all ROIs have valid source estimates + +## Troubleshooting + +### Common Issues + +1. **Singular covariance matrix**: + - Increase regularization parameter (lambda) + - Use more trials for covariance estimation + - Apply covariance shrinkage + +2. **Memory errors**: + - Reduce source space resolution + - Process in batches + - Use sparse matrices where applicable + +3. **Poor localization**: + - Check head model alignment + - Verify leadfield quality + - Ensure sufficient SNR in data + +4. **Missing ROIs**: + - Verify parcellation atlas + - Check surface mesh quality + - Ensure FreeSurfer processing completed + +### Optimization + +For faster processing: + +```matlab +% Reduce source space (in head model creation) +cfg.grid.resolution = 10; % mm (larger = faster, less accurate) + +% Fewer PCs per ROI +n_components = 2; % Instead of 3 + +% Parallel processing +parfor sbj = 1:n_subjects + % Process subject +end +``` + +## Advanced Usage + +### Custom ROIs + +Define custom regions of interest: + +```matlab +% Load source space +source_space = load('HeadModel/SubjectID/source_space.mat'); + +% Define ROI by MNI coordinates +roi_center = [20, -30, 50]; % MNI coordinates +roi_radius = 15; % mm + +% Find sources within radius +distances = vecnorm(source_space.pos - roi_center, 2, 2); +roi_sources = find(distances < roi_radius); + +% Extract ROI time series +roi_data = source_timeseries(roi_sources, :); +``` + +### Alternative Source Methods + +The pipeline can be adapted for other source reconstruction methods: + +**Minimum Norm Estimate (MNE)**: +```matlab +cfg.method = 'mne'; +cfg.lambda = 1e-3; +source = ft_sourceanalysis(cfg, data); +``` + +**Dynamic Imaging of Coherent Sources (DICS)**: +```matlab +cfg.method = 'dics'; +cfg.frequency = 10; % Hz +source = ft_sourceanalysis(cfg, freq); +``` + +## Notes + +- LCMV beamformer is optimal for sources with distinct spatial patterns +- Requires good head model and sufficient SNR +- Results are sensitive to regularization parameter +- Multiple active sources can cause cancellation artifacts +- Consider source reconstruction validation before further analysis diff --git a/Stat/README.md b/Stat/README.md new file mode 100644 index 0000000..f9d276e --- /dev/null +++ b/Stat/README.md @@ -0,0 +1,471 @@ +# Statistical Analysis Pipeline + +This directory contains scripts for group-level statistical analysis of PAC and power data, comparing Alzheimer's Disease patients with healthy controls. + +## Overview + +The statistical pipeline performs group-level comparisons, correlations with clinical measures, and visualization of results. It includes both frequentist and mixed-effects modeling approaches. + +## Main Scripts + +### prepareTable.m + +**Purpose**: Prepares a comprehensive data table combining PAC values, power measures, and clinical/demographic information. + +**Input**: +- PAC results from `Results/SeedtoCortex/[SubjectID]/` +- Source power data from `Results/source/[SubjectID]/` +- Clinical data (demographics, MMSE scores, diagnosis) + +**Output**: +- Structured table with all variables for statistical analysis +- Saved as MAT file and CSV for use in MATLAB and external software (R, Python) + +**Data Table Structure**: +```matlab +% Columns: +% - Subject ID +% - Group (AD vs Control) +% - Age, Sex, Education +% - MMSE score (cognitive function) +% - PAC values (per ROI, frequency, condition) +% - Power values (per ROI, frequency band) +% - Derived metrics (asymmetry, ratios) +``` + +**Usage**: +```matlab +prepareTable; +``` + +### powerAnalysis.m + +**Purpose**: Analyzes power spectral density (PSD) in cortical and hippocampal regions. + +**Analyses Performed**: + +1. **PSD Computation**: + - Calculate power in standard frequency bands + - Cortical ROIs and hippocampus + - Delta (1-4 Hz), Theta (4-8 Hz), Alpha (8-12 Hz), Beta (13-30 Hz), Gamma (30-55 Hz) + +2. **Group Comparisons**: + - AD vs Control power differences + - Region-wise statistical tests + - Multiple comparison correction + +3. **Visualization**: + - Topographic maps of power differences + - Power spectrum plots per ROI + - Effect size distributions + +**Output**: +- Power values per subject, ROI, and frequency band +- Statistical test results (t-tests, effect sizes) +- Publication-quality figures + +**Usage**: +```matlab +powerAnalysis; +``` + +### clusterBasedModel.m + +**Purpose**: Performs cluster-based permutation testing on PAC differences between AD and control groups. + +**Method**: Cluster-based permutation test (non-parametric) + +**Why**: Controls family-wise error rate (FWER) in the presence of multiple comparisons while maintaining sensitivity to spatially/spectrally contiguous effects. + +**Algorithm**: +``` +1. Compute observed statistic (t-statistic) for each ROI and frequency +2. Identify clusters of contiguous significant effects +3. Permutation testing: + a. Randomly shuffle group labels + b. Recompute statistics + c. Find maximum cluster statistic + d. Repeat 1000+ times +4. Compare observed clusters to null distribution +5. Report significant clusters (p < 0.05) +``` + +**Analyses**: + +1. **PAC Group Differences**: + - AD vs Control PAC values + - Across cortical ROIs and frequencies + - Both hippocampi analyzed separately + +2. **Correlation with MMSE**: + - PAC vs cognitive function + - Within AD group + - Identifies regions related to cognitive decline + +3. **Visualization**: + - Significant ROIs plotted on cortical surface + - Frequency-by-frequency cluster maps + - Scatter plots for correlations + +**Input**: +- Data table from `prepareTable.m` +- PAC values per subject, ROI, frequency +- Group labels and MMSE scores + +**Output**: +- Significant clusters (ROIs, frequency bands) +- Permutation test statistics +- Plots of results +- Saved in `Stat/ClusterTest/` + +**Usage**: +```matlab +clusterBasedModel; +``` + +### linearMixedEffectsModel.m + +**Purpose**: Fits linear mixed-effects (LME) models to account for repeated measures and covariates. + +**Method**: Linear Mixed-Effects Modeling + +**Why**: +- Accounts for within-subject correlation (repeated measures across ROIs, frequencies) +- Includes random effects for subject-specific variability +- Can include covariates (age, sex, education) +- Provides estimates of fixed effects (group differences) and random effects + +**Model Structure**: +``` +PAC ~ Group + Age + Sex + MMSE + (1 | Subject) + (1 | ROI) + +Fixed Effects: +- Group: AD vs Control +- Age: Continuous +- Sex: Binary +- MMSE: Cognitive function + +Random Effects: +- Subject: Random intercept per subject +- ROI: Random intercept per cortical region +``` + +**Analyses**: + +1. **Main Effects**: + - Group effect on PAC + - Covariate effects + - Interaction terms if specified + +2. **Model Comparison**: + - Likelihood ratio tests + - AIC/BIC criteria + - R² (variance explained) + +3. **Post-hoc Tests**: + - Pairwise comparisons + - Simple effects analysis + +**Input**: +- Data table from `prepareTable.m` +- Model specification (formula) + +**Output**: +- Model fit statistics +- Fixed effect estimates with confidence intervals +- Random effect variances +- Diagnostic plots (residuals, Q-Q plots) +- Saved in `Stat/LME/` + +**Usage**: +```matlab +linearMixedEffectsModel; +``` + +## Helper Scripts + +### cluster_based_perm.py + +**Purpose**: Python implementation of cluster-based permutation testing (alternative to MATLAB). + +**Dependencies**: +- NumPy, SciPy +- MNE-Python (optional, for visualization) + +**Usage**: +```python +python cluster_based_perm.py --input pac_data.csv --output results/ +``` + +### Fig_creator.py + +**Purpose**: Creates publication-quality figures from statistical results. + +**Features**: +- Cortical surface plots with significant ROIs +- Frequency-by-frequency heatmaps +- Scatter plots with regression lines +- Customizable color schemes and layouts + +**Usage**: +```python +python Fig_creator.py --results cluster_results.mat --output figures/ +``` + +## Supporting Data Files + +### roi_name.mat + +**Content**: ROI names and labels for cortical regions. + +**Structure**: +```matlab +roi_names = { + 'Left-Hippocampus', + 'Right-Hippocampus', + 'Left-Superior-Frontal', + 'Right-Superior-Frontal', + ... +}; +``` + +### dk_labels.mat + +**Content**: Desikan-Killiany atlas label information. + +**Structure**: +- Label IDs +- ROI names +- Anatomical groupings +- MNI coordinates + +### cm17.mat + +**Content**: 17-color colormap for cortical parcellations. + +**Usage**: Visualization of ROI-specific results on brain surface. + +### mask.mat + +**Content**: Binary masks for including/excluding ROIs or frequency bins. + +**Usage**: Apply specific inclusion criteria in statistical tests. + +## Subdirectories + +### ClusterTest/ + +Contains results from cluster-based permutation testing: +- `significant_clusters.mat`: Cluster statistics +- `cluster_plots/`: Visualization of significant effects +- `permutation_distribution.mat`: Null distribution data + +### LME/ + +Contains results from linear mixed-effects modeling: +- `model_fits.mat`: Model parameters and statistics +- `diagnostic_plots/`: Residual and Q-Q plots +- `fixed_effects.csv`: Effect estimates with CI + +### Power/ + +Contains power analysis results: +- `power_spectra.mat`: PSD per subject and ROI +- `group_differences.mat`: Statistical test results +- `power_plots/`: Topographic and spectral plots + +## Expected Runtime + +Statistical analyses are generally fast: +- **prepareTable**: ~5-10 minutes (depending on number of subjects) +- **powerAnalysis**: ~10-20 minutes +- **clusterBasedModel**: ~30-120 minutes (depends on number of permutations) +- **linearMixedEffectsModel**: ~10-30 minutes (depends on model complexity) + +## Configuration + +### Cluster-Based Permutation + +```matlab +% Number of permutations +n_permutations = 5000; % More = better p-value resolution + +% Cluster-forming threshold +cluster_threshold = 0.05; % Initial p-value for cluster formation + +% Cluster statistic +cluster_stat = 'maxsum'; % 'maxsum', 'maxsize', or 'tfce' + +% Multiple comparison correction +correction_method = 'cluster'; % 'cluster', 'fdr', or 'bonferroni' +``` + +### Linear Mixed-Effects + +```matlab +% Model formula +formula = 'PAC ~ Group + Age + Sex + MMSE + (1|Subject) + (1|ROI)'; + +% Fitting method +fitmethod = 'REML'; % 'REML' or 'ML' + +% Optimizer +optimizer = 'fminunc'; % MATLAB optimization function +``` + +## Statistical Considerations + +### Multiple Comparisons + +**Problem**: Testing many hypotheses increases false positive rate. + +**Solutions Implemented**: + +1. **Cluster-based permutation**: Controls FWER while preserving sensitivity +2. **FDR correction**: Controls false discovery rate (less conservative) +3. **Bonferroni**: Most conservative, suitable for few planned comparisons + +### Effect Size + +Report effect sizes along with p-values: + +```matlab +% Cohen's d for group differences +cohens_d = (mean_AD - mean_Control) / pooled_std; + +% Interpretation: +% Small: d = 0.2 +% Medium: d = 0.5 +% Large: d = 0.8 +``` + +### Power Analysis + +Ensure sufficient sample size: + +```matlab +% Required sample size for 80% power, α=0.05, d=0.5 +n_required = sampsizepwr('t', [0, 1], 0.5, 0.8, [], 'Alpha', 0.05); +``` + +## Troubleshooting + +### Common Issues + +1. **Convergence failures in LME**: + - Simplify random effects structure + - Center and scale continuous predictors + - Check for multicollinearity + +2. **No significant clusters**: + - May indicate genuine null results + - Check effect sizes (may be small but real) + - Ensure sufficient power + +3. **Memory errors in permutation testing**: + - Reduce number of permutations + - Process frequency bands separately + - Use cluster computing if available + +4. **Unbalanced groups**: + - Use robust statistics + - Consider weighting or matching + - Report group sizes clearly + +## Validation + +### Sanity Checks + +```matlab +% Check data distribution +figure; histogram(pac_values); +title('PAC Distribution'); + +% Check for outliers +outliers = isoutlier(pac_values, 'median'); +fprintf('Detected %d outliers (%.1f%%)\n', sum(outliers), ... + 100*sum(outliers)/length(pac_values)); + +% Check group balance +group_counts = groupcounts(data.Group); +disp(group_counts); + +% Check for missing data +missing_data = sum(isnan(data_table), 2); +fprintf('Subjects with missing data: %d\n', sum(missing_data > 0)); +``` + +## Visualization Best Practices + +1. **Use consistent color schemes**: AD=red, Control=blue +2. **Include error bars**: SEM or 95% CI +3. **Report statistics on plots**: p-values, effect sizes +4. **Use appropriate scales**: Consider log scale for power +5. **Add legends and labels**: Clear axis labels and units + +## Output for Publication + +### Tables + +- Table 1: Demographics and clinical characteristics +- Table 2: PAC group differences (mean ± SD, statistics) +- Table 3: Significant clusters (location, frequency, statistics) +- Table 4: LME model results (fixed effects estimates) + +### Figures + +- Figure 1: Power spectra by group and ROI +- Figure 2: Significant PAC differences on brain surface +- Figure 3: Frequency-specific PAC clusters +- Figure 4: PAC vs MMSE correlations +- Figure 5: Effect sizes across ROIs + +## Advanced Usage + +### Custom Contrasts + +Test specific hypotheses: + +```matlab +% Test specific contrast +contrast = [1 -1 0 0]; % AD - Control, ignore covariates +[F, p] = coefTest(lme, contrast); +``` + +### Interaction Effects + +Include interactions in model: + +```matlab +% Group × MMSE interaction +formula = 'PAC ~ Group * MMSE + Age + Sex + (1|Subject)'; +lme = fitlme(data_table, formula); + +% Test interaction +[p, F] = coefTest(lme, [0 0 0 0 1]); % Interaction term +``` + +### Sensitivity Analysis + +Test robustness of results: + +```matlab +% Without outliers +data_clean = data_table(~outliers, :); +lme_clean = fitlme(data_clean, formula); + +% With different covariates +formula_alt = 'PAC ~ Group + Age + (1|Subject)'; +lme_alt = fitlme(data_table, formula_alt); + +% Compare results +compare(lme, lme_alt); +``` + +## Notes + +- Always visualize data before statistical testing +- Report effect sizes, not just p-values +- Consider biological significance, not just statistical significance +- Document all analysis choices for reproducibility +- Use consistent statistical thresholds across analyses