Browse Source

Update 'code/homotopy_denoising.m'

Elliot Layden 6 months ago
parent
commit
e10b8f7a33
1 changed files with 280 additions and 286 deletions
  1. 280 286
      code/homotopy_denoising.m

+ 280 - 286
code/homotopy_denoising.m

@@ -1,286 +1,280 @@
-% homotopy_denoising.m
-
-% Within-subjects:
-% 1. Linear trends
-% 2. Quadratic trends
-% 3. 6 rigid-body motion parameters
-% 4. 5 principal components of CSF signals
-% 5. Body temperature
-% 6. Isoflurane dose
-% 7. Global signal regression (optional)
-% 8. Bandpass filtering
-
-% Between-subjects:
-% 1. Age
-% 2. Mean Body Temperature
-% 3. Mean Isoflurane Dose
-% 4. ROI pair Euclidean distance
-% 5. ROI volume
-
-% Note: requires Tools for NIfTI and ANALYZE image
-% (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)
-
-% Use smoothed or unsmoothed data for ROI extraction?
-useSmoothedData = false;
-
-% Use homotopic-SC as well as homotopic-nSC?
-useSC = true;
-
-% Regress out global signal?
-useGlobalSignal = false;
-
-%% LOAD DATA
-
-% Collect file paths to normalized functional images:
-functPath = '';
-nscan = 19;
-timeSeries = cell(1,nscan); % 5 normals, 4 scans each (except normal #5, 3 scans)
-for scan = 1:nscan
-    timeSeries{scan} = cell(1,175);
-    for timepoint = 6:180 % skip first 5 volumes (T1 equilibration)
-        if useSmoothedData
-            timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
-                sprintf('s_w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
-        else
-            timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
-                sprintf('w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
-        end
-    end
-end
-
-% Load ROIs:
-ROI_path = '';
-SongSystem = fullfile(ROI_path,'SongSystemROIs.nii'); % n = 20
-mask = load_nii(fullfile(ROI_path,'new_temp_brain_mask.nii'));
-nROI = 20;
-
-% Get ROI Subscripts/Coordinates:
-rois = load_nii(SongSystem);
-roi_subscripts = cell(1,nROI);
-for i = 1:nROI
-    [x, y, z] = ind2sub(size(rois.img),find(rois.img(:)==i));
-    roi_subscripts{i} = [x,y,z];
-end
-
-% Get CSF Subscripts/Coordinates:
-CSF_mask = load_nii(fullfile(ROI_path, 'CSF_mask.nii'));
-[CSF_x,CSF_y,CSF_z] = ind2sub(size(CSF_mask.img), find(CSF_mask.img(:)>0));
-n_CSF = length(CSF_x); 
-
-% Extract ROI Data & CSF Signals:
-ROI_dat = cell(1,nscan); CSF = cell(1,nscan); globalSignal = zeros(175,nscan); % Initialize
-hWait = waitbar(0,'Extracting functional data...');
-for scan = 1:nscan
-    waitbar(scan/nscan, hWait);
-    
-    ROI_dat{scan} = zeros(175, nROI); CSF_dat = zeros(175, n_CSF); % Initialize
-    for j = 1:175 % iterate through time (j)
-        
-        img = load_nii(timeSeries{scan}{j});
-        globalSignal(j, scan) = mean(img.img(mask.img(:)>0)); % Average whole brain signal for time point j
-        
-        for k = 1:nROI % iterate through ROIs (k)
-            nvox = size(roi_subscripts{k}, 1); % # voxels in ROI
-            roi_holder = zeros(1, nvox);
-            for l = 1:nvox % iterate through voxels of ROI
-                roi_holder(l) = img.img(roi_subscripts{k}(l,1), roi_subscripts{k}(l,2), roi_subscripts{k}(l,3));
-            end
-            ROI_dat{scan}(j,k) = mean(roi_holder); % average voxels of ROI
-        end    
-        
-        for k = 1:n_CSF % iterate through CSF voxels (k)
-            CSF_dat(j,k) = img.img(CSF_x(k), CSF_y(k), CSF_z(k));
-        end  
-    end
-    
-    [~, PCs] = pca(CSF_dat); % pca across voxels of CSF mask
-    CSF{scan} = PCs(:,1:5); % save PC's 1-5
-    
-end; delete(hWait)
-
-% Load Motion Data:
-motionPath = '';
-motion_params = cell(1,nscan);
-for i = 1:nscan
-    load(fullfile(motionPath, sprintf('motion_scan_%02g.mat',i)))
-    motion_params{i} = motion;
-end
-
-% Load Physiology Data (Temperature (degrees Celsius), Isoflurane (%)):
-physio_path = '';
-physio_params = cell(1,nscan);
-physio_times = [60:60:540,576]; % 60 seconds to 576 seconds
-TR = 1.0671*3; scan_times = TR:TR:TR*175;
-for i = 1:nscan
-    tbl = readtable(fullfile(physio_path,'All_Normal_Condition_Physiology_Records.xlsx'),'Sheet',i);
-    pp_temp = spline(physio_times,tbl.Temperature); % piecewise polynomial of cubic spline
-    pred_temp = ppval(pp_temp,scan_times); % evaluate PP at x
-    pp_iso = spline(physio_times,tbl.Isoflurane); % piecewise polynomial of cubic spline
-    pred_iso = ppval(pp_iso,scan_times); % evaluate PP at x
-    physio_params{i} = [pred_temp',pred_iso'];
-end
-
-% Labels:
-SongSystem_labels = {'L Area X','R Area X','L LMAN','R LMAN','L RA','R RA',...
-    'L HVC','R HVC','L Aud. Forebrain','R Aud. Forebrain','L Field L',...
-    'R Field L','L Diencephalon','R Diencephalon','L Medulla','R Medulla',...
-    'L Uva','R Uva','L DM','R DM'};
-
-ages = [25, 45, 65, 89, 24, 45, 65, 90, 25, 46, 67, 89, 24, 46, 66, 90, 24, 45, 67];
-
-% Get ROI-level tSNR Data:
-tSNR_mat = zeros(nscan, nROI); 
-for i = 1:nscan
-    for j = 1:nROI
-        mean_sig = mean(ROI_dat{i}(:,j));
-        noise_SD = std(ROI_dat{i}(:,j));
-        tSNR_mat(i,j) = mean_sig ./ noise_SD;
-    end
-end
-
-% Get Mean Physiological Parameters:
-temps = zeros(19,1); iso = zeros(19,1);
-for i = 1:19
-    temps(i) = mean(physio_params{i}(:,1));
-    iso(i) = mean(physio_params{i}(:,2));
-end
-
-% Get Mean Volumes and Euclidean Distances:
-vox_dim = rois.hdr.dime.pixdim(2:4); 
-vox_volume = prod(vox_dim);
-m_volumes = zeros(nROI); euclideans = zeros(nROI);
-
-% Get ROI Centroids:
-centroids = zeros(nROI,3);
-for i = 1:nROI
-    [x,y,z] = ind2sub(size(rois.img),find(rois.img==i));
-    centroids(i,:) = [mean(x),mean(y),mean(z)].*vox_dim;
-end
-for i = 1:nROI
-    for j = 1:nROI
-        m_volumes(i,j) = mean([sum(rois.img(:)==i)*vox_volume,sum(rois.img(:)==j)*vox_volume]);
-        euclideans(i,j) = sqrt(sum((centroids(i,:)-centroids(j,:)).^2));
-    end
-end
-
-%% RUN DENOISING
-
-% Labels for FC:
-labels = {'L_Area_X','R_Area_X','L_LMAN','R_LMAN','L_RA','R_RA',...
-    'L_HVC','R_HVC','L_Aud_Forebrain','R_Aud_Forebrain','L_Field_L',...
-    'R_Field_L','L_Diencephalon','R_Diencephalon','L_Medulla','R_Medulla',...
-    'L_Uva','R_Uva','L_DM','R_DM'};
-
-useScans = 1:19;
-
-if useSC
-    nROI = 20;
-else
-    nROI = 12;
-end
-
-% Preprocess ROI Data / Within-Subjects Nuisance Regression:
-Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2; intercept = ones(175,1);
-hWait = waitbar(0,'Denoising...');
-for scan = useScans
-    if ishandle(hWait); waitbar(scan/nscan, hWait); end
-    
-    % Get nuisance regression covariates:
-    if all(physio_params{scan}(:,2)==mode(physio_params{scan}(:,2))) % iso is a constant, so don't need
-        if useGlobalSignal
-            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1), globalSignal(:,scan)];
-        else
-            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1)];
-        end
-    else
-        if useGlobalSignal
-            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}, globalSignal(:,scan)];
-        else
-            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}];
-        end
-    end
-    
-    % Perform nuisance regression, retain residuals:
-    for j = 1:nROI
-        [~, ~, ROI_dat{scan}(:,j)] = regress(ROI_dat{scan}(:,j), nuisance);
-    end
-    
-    % Perform bandpass filtering: 
-    [ROI_dat{scan},~,~] = bst_bandpass_filtfilt(ROI_dat{scan}', Fs, .008, .1, 0, 'iir');
-    ROI_dat{scan} = ROI_dat{scan}';
-
-end
-if ishandle(hWait); delete(hWait); end 
-
-% Calculate Correlation Matrix:
-zmats = zeros(nROI, nROI, nscan);
-for scan = useScans
-    zmats(:,:,scan) = fisherz(corr(ROI_dat{scan}(:,1:nROI)));
-end
-
-% Between-Subjects Nuisance Regression:
-nscan = length(useScans);
-numconn = nROI*(nROI-1)/2*nscan;
-birds = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5];
-scan = zeros(numconn,1); bird = zeros(numconn,1);  
-FC = zeros(numconn,1); age = zeros(numconn,1); 
-connection_type = zeros(numconn,1); euclidean = zeros(numconn,1); 
-m_volume = zeros(numconn,1); mean_temp = zeros(numconn,1); 
-mean_iso = zeros(numconn,1); conn_label = cell(numconn,1); count = 0; 
-for i = 1:nROI-1
-    for j = i+1:nROI
-        FC(count+1:count+nscan) = squeeze(zmats(i,j,:));
-        conn_label(count+1:count+nscan) = {[labels{i},'_',labels{j}]};
-        age(count+1:count+nscan) = ages(useScans);
-        bird(count+1:count+nscan) = birds(useScans);
-        scan(count+1:count+nscan) = useScans;
-        euclidean(count+1:count+nscan) = repmat(euclideans(i,j),nscan,1);
-        m_volume(count+1:count+nscan) = repmat(m_volumes(i,j),nscan,1);
-        if ((mod(i,2) && ~mod(j,2)) || (~mod(i,2) && mod(j,2))) && ~(mod(i,2) && j==(i+1)) % 4: heterotopic (contralateral)
-            connection_type(count+1:count+nscan) = 4*ones(nscan,1);
-        elseif mod(i,2) && j==(i+1)
-            if i>12 || j>12
-                connection_type(count+1:count+nscan) = 3*ones(nscan,1); % 3: homotopic-SC
-            else
-                connection_type(count+1:count+nscan) = 2*ones(nscan,1); % 2: homotopic-nSC
-            end
-        elseif (mod(i,2) && mod(j,2)) || (~mod(i,2) && ~mod(j,2)) % 1: ipsilateral
-            connection_type(count+1:count+nscan) = ones(nscan,1);
-        end
-        mean_temp(count+1:count+nscan) = temps(useScans);
-        mean_iso(count+1:count+nscan) = iso(useScans);
-        count = count+nscan;
-    end
-end
-bird = categorical(bird); scan = categorical(scan);
-if useSC
-    connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic-nSC','Homotopic-SC','Heterotopic'});
-else
-    connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic','Homotopic','Heterotopic'});
-end
-
-% Mixed Model:
-FC_data = table(FC, bird, scan, conn_label, connection_type, age, mean_temp, ...
-    mean_iso, euclidean, m_volume, 'VariableNames', {'FC', 'Bird', 'Scan', ...
-    'Connection', 'Type', 'Age', 'Temp', 'Iso', 'Euclidean', 'Volume'});
-
-% Select Between-subjects nuisance parameters for formula:
-FC_formula = 'FC ~ Type*Age + Type*Iso + Temp + Euclidean + Volume + (1|Bird)';
-
-% Run LME random intercept model:
-lme = fitlme(FC_data, FC_formula,  'StartMethod','random',  'verbose',true,  'CheckHessian',true); 
-lme.display
-
-% Get random effect coefficients:
-[B, Bnames, stats] = randomEffects(lme) %#ok
-
-% Check residuals:
-figure; lme.plotResiduals
-
-% Get marginal (fixed-effects only) vs. conditional (w/ random) R2
-lme.Rsquared
-
-resid_marginal = residuals(lme, 'Conditional',false); % fixed-effects only
-resid_conditional = residuals(lme, 'Conditional',true); % fixed-effects only
-R2_marginal = 1 - sum(resid_marginal.^2) / lme.SST
-R2_conditional = 1 - sum(resid_conditional.^2) / lme.SST
+% homotopy_denoising.m
+
+% Within-subjects:
+% 1. Linear trends
+% 2. Quadratic trends
+% 3. 6 rigid-body motion parameters
+% 4. 5 principal components of CSF signals
+% 5. Body temperature
+% 6. Isoflurane dose
+% 7. Global signal regression (optional)
+% 8. Bandpass filtering
+
+% Between-subjects:
+% 1. Age
+% 2. Mean Body Temperature
+% 3. Mean Isoflurane Dose
+% 4. ROI pair Euclidean distance
+% 5. ROI volume
+
+% Note: requires Tools for NIfTI and ANALYZE image
+% (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)
+
+% Use smoothed or unsmoothed data for ROI extraction?
+useSmoothedData = false;
+
+% Use homotopic-SC as well as homotopic-nSC?
+useSC = true;
+
+% Regress out global signal?
+useGlobalSignal = false;
+
+%% LOAD DATA
+
+% Collect file paths to normalized functional images:
+functPath = '';
+nscan = 19;
+timeSeries = cell(1,nscan); % 5 normals, 4 scans each (except normal #5, 3 scans)
+for scan = 1:nscan
+    timeSeries{scan} = cell(1,175);
+    for timepoint = 6:180 % skip first 5 volumes (T1 equilibration)
+        if useSmoothedData
+            timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
+                sprintf('s_w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
+        else
+            timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
+                sprintf('w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
+        end
+    end
+end
+
+% Load ROIs:
+ROI_path = '';
+SongSystem = fullfile(ROI_path,'SongSystemROIs.nii'); % n = 20
+mask = load_nii(fullfile(ROI_path,'new_temp_brain_mask.nii'));
+nROI = 20;
+
+% Get ROI Subscripts/Coordinates:
+rois = load_nii(SongSystem);
+roi_subscripts = cell(1,nROI);
+for i = 1:nROI
+    [x, y, z] = ind2sub(size(rois.img),find(rois.img(:)==i));
+    roi_subscripts{i} = [x,y,z];
+end
+
+% Get CSF Subscripts/Coordinates:
+CSF_mask = load_nii(fullfile(ROI_path, 'CSF_mask.nii'));
+[CSF_x,CSF_y,CSF_z] = ind2sub(size(CSF_mask.img), find(CSF_mask.img(:)>0));
+n_CSF = length(CSF_x); 
+
+% Extract ROI Data & CSF Signals:
+ROI_dat = cell(1,nscan); CSF = cell(1,nscan); globalSignal = zeros(175,nscan); % Initialize
+hWait = waitbar(0,'Extracting functional data...');
+for scan = 1:nscan
+    waitbar(scan/nscan, hWait);
+    
+    ROI_dat{scan} = zeros(175, nROI); CSF_dat = zeros(175, n_CSF); % Initialize
+    for j = 1:175 % iterate through time (j)
+        
+        img = load_nii(timeSeries{scan}{j});
+        globalSignal(j, scan) = mean(img.img(mask.img(:)>0)); % Average whole brain signal for time point j
+        
+        for k = 1:nROI % iterate through ROIs (k)
+            nvox = size(roi_subscripts{k}, 1); % # voxels in ROI
+            roi_holder = zeros(1, nvox);
+            for l = 1:nvox % iterate through voxels of ROI
+                roi_holder(l) = img.img(roi_subscripts{k}(l,1), roi_subscripts{k}(l,2), roi_subscripts{k}(l,3));
+            end
+            ROI_dat{scan}(j,k) = mean(roi_holder); % average voxels of ROI
+        end    
+        
+        for k = 1:n_CSF % iterate through CSF voxels (k)
+            CSF_dat(j,k) = img.img(CSF_x(k), CSF_y(k), CSF_z(k));
+        end  
+    end
+    
+    [~, PCs] = pca(CSF_dat); % pca across voxels of CSF mask
+    CSF{scan} = PCs(:,1:5); % save PC's 1-5
+    
+end; delete(hWait)
+
+% Load Motion Data:
+motionPath = '';
+motion_params = cell(1,nscan);
+for i = 1:nscan
+    load(fullfile(motionPath, sprintf('motion_scan_%02g.mat',i)))
+    motion_params{i} = motion;
+end
+
+% Load Physiology Data (Temperature (degrees Celsius), Isoflurane (%)):
+physio_path = '';
+physio_params = cell(1,nscan);
+physio_times = [60:60:540,576]; % 60 seconds to 576 seconds
+TR = 1.0671*3; scan_times = TR:TR:TR*175;
+for i = 1:nscan
+    tbl = readtable(fullfile(physio_path,'All_Normal_Condition_Physiology_Records.xlsx'),'Sheet',i);
+    pp_temp = spline(physio_times,tbl.Temperature); % piecewise polynomial of cubic spline
+    pred_temp = ppval(pp_temp,scan_times); % evaluate PP at x
+    pp_iso = spline(physio_times,tbl.Isoflurane); % piecewise polynomial of cubic spline
+    pred_iso = ppval(pp_iso,scan_times); % evaluate PP at x
+    physio_params{i} = [pred_temp',pred_iso'];
+end
+
+ages = [25, 45, 65, 89, 24, 45, 65, 90, 25, 46, 67, 89, 24, 46, 66, 90, 24, 45, 67];
+
+% Get ROI-level tSNR Data:
+tSNR_mat = zeros(nscan, nROI); 
+for i = 1:nscan
+    for j = 1:nROI
+        mean_sig = mean(ROI_dat{i}(:,j));
+        noise_SD = std(ROI_dat{i}(:,j));
+        tSNR_mat(i,j) = mean_sig ./ noise_SD;
+    end
+end
+
+% Get Mean Physiological Parameters:
+temps = zeros(19,1); iso = zeros(19,1);
+for i = 1:19
+    temps(i) = mean(physio_params{i}(:,1));
+    iso(i) = mean(physio_params{i}(:,2));
+end
+
+% Get Mean Volumes and Euclidean Distances:
+vox_dim = rois.hdr.dime.pixdim(2:4); 
+vox_volume = prod(vox_dim);
+m_volumes = zeros(nROI); euclideans = zeros(nROI);
+
+% Get ROI Centroids:
+centroids = zeros(nROI,3);
+for i = 1:nROI
+    [x,y,z] = ind2sub(size(rois.img),find(rois.img==i));
+    centroids(i,:) = [mean(x),mean(y),mean(z)].*vox_dim;
+end
+for i = 1:nROI
+    for j = 1:nROI
+        m_volumes(i,j) = mean([sum(rois.img(:)==i)*vox_volume,sum(rois.img(:)==j)*vox_volume]);
+        euclideans(i,j) = sqrt(sum((centroids(i,:)-centroids(j,:)).^2));
+    end
+end
+
+%% RUN DENOISING
+
+% Labels for FC:
+labels = {'R Area X','L Area X','R LMAN','L LMAN','R RA','L RA',...
+    'R HVC','L HVC','R Aud. Forebrain','L Aud. Forebrain','R Field L',...
+    'L Field L','R Diencephalon','L Diencephalon','R Medulla','L Medulla',...
+    'R Uva','L Uva','R DM','L DM'};
+
+useScans = 1:19;
+
+if useSC
+    nROI = 20;
+else
+    nROI = 12;
+end
+
+% Preprocess ROI Data / Within-Subjects Nuisance Regression:
+Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2; intercept = ones(175,1);
+hWait = waitbar(0,'Denoising...');
+for scan = useScans
+    if ishandle(hWait); waitbar(scan/nscan, hWait); end
+    
+    % Get nuisance regression covariates:
+    if all(physio_params{scan}(:,2)==mode(physio_params{scan}(:,2))) % iso is a constant, so don't need
+        if useGlobalSignal
+            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1), globalSignal(:,scan)];
+        else
+            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1)];
+        end
+    else
+        if useGlobalSignal
+            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}, globalSignal(:,scan)];
+        else
+            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}];
+        end
+    end
+    
+    % Perform nuisance regression, retain residuals:
+    for j = 1:nROI
+        [~, ~, ROI_dat{scan}(:,j)] = regress(ROI_dat{scan}(:,j), nuisance);
+    end
+    
+    % Perform bandpass filtering: 
+    [ROI_dat{scan},~,~] = bst_bandpass_filtfilt(ROI_dat{scan}', Fs, .008, .1, 0, 'iir');
+    ROI_dat{scan} = ROI_dat{scan}';
+
+end
+if ishandle(hWait); delete(hWait); end 
+
+% Calculate Correlation Matrix:
+zmats = zeros(nROI, nROI, nscan);
+for scan = useScans
+    zmats(:,:,scan) = fisherz(corr(ROI_dat{scan}(:,1:nROI)));
+end
+
+% Between-Subjects Nuisance Regression:
+nscan = length(useScans);
+numconn = nROI*(nROI-1)/2*nscan;
+birds = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5];
+scan = zeros(numconn,1); bird = zeros(numconn,1);  
+FC = zeros(numconn,1); age = zeros(numconn,1); 
+connection_type = zeros(numconn,1); euclidean = zeros(numconn,1); 
+m_volume = zeros(numconn,1); mean_temp = zeros(numconn,1); 
+mean_iso = zeros(numconn,1); conn_label = cell(numconn,1); count = 0; 
+for i = 1:nROI-1
+    for j = i+1:nROI
+        FC(count+1:count+nscan) = squeeze(zmats(i,j,:));
+        conn_label(count+1:count+nscan) = {[labels{i},'_',labels{j}]};
+        age(count+1:count+nscan) = ages(useScans);
+        bird(count+1:count+nscan) = birds(useScans);
+        scan(count+1:count+nscan) = useScans;
+        euclidean(count+1:count+nscan) = repmat(euclideans(i,j),nscan,1);
+        m_volume(count+1:count+nscan) = repmat(m_volumes(i,j),nscan,1);
+        if ((mod(i,2) && ~mod(j,2)) || (~mod(i,2) && mod(j,2))) && ~(mod(i,2) && j==(i+1)) % 4: heterotopic (contralateral)
+            connection_type(count+1:count+nscan) = 4*ones(nscan,1);
+        elseif mod(i,2) && j==(i+1)
+            if i>12 || j>12
+                connection_type(count+1:count+nscan) = 3*ones(nscan,1); % 3: homotopic-SC
+            else
+                connection_type(count+1:count+nscan) = 2*ones(nscan,1); % 2: homotopic-nSC
+            end
+        elseif (mod(i,2) && mod(j,2)) || (~mod(i,2) && ~mod(j,2)) % 1: ipsilateral
+            connection_type(count+1:count+nscan) = ones(nscan,1);
+        end
+        mean_temp(count+1:count+nscan) = temps(useScans);
+        mean_iso(count+1:count+nscan) = iso(useScans);
+        count = count+nscan;
+    end
+end
+bird = categorical(bird); scan = categorical(scan);
+if useSC
+    connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic-nSC','Homotopic-SC','Heterotopic'});
+else
+    connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic','Homotopic','Heterotopic'});
+end
+
+% Mixed Model:
+FC_data = table(FC, bird, scan, conn_label, connection_type, age, mean_temp, ...
+    mean_iso, euclidean, m_volume, 'VariableNames', {'FC', 'Bird', 'Scan', ...
+    'Connection', 'Type', 'Age', 'Temp', 'Iso', 'Euclidean', 'Volume'});
+
+% Select Between-subjects nuisance parameters for formula:
+FC_formula = 'FC ~ Type*Age + Type*Iso + Temp + Euclidean + Volume + (1|Bird)';
+
+% Run LME random intercept model:
+lme = fitlme(FC_data, FC_formula,  'StartMethod','random',  'verbose',true,  'CheckHessian',true); 
+lme.display
+
+% Get random effect coefficients:
+[B, Bnames, stats] = randomEffects(lme) %#ok
+
+% Check residuals:
+figure; lme.plotResiduals
+
+% Get marginal (fixed-effects only) vs. conditional (w/ random) R2
+lme.Rsquared
+
+resid_marginal = residuals(lme, 'Conditional',false); % fixed-effects only
+resid_conditional = residuals(lme, 'Conditional',true); % fixed-effects only
+R2_marginal = 1 - sum(resid_marginal.^2) / lme.SST
+R2_conditional = 1 - sum(resid_conditional.^2) / lme.SST