Unique brain systems architecture characterizes affective sequelae after traumatic brain injury (Science Transl Med, 2023)

%Variable definitions:
%TBI = history of TBI (1=yes, 0=no)
%MDD = presence vs absence of significant depressive symptoms
%DepSeverity = continuous measure of depression severity
%PTSD = presence vs absence of PTSD
%Dataset = dataset label

%DAN_SG, DAN_DMN, and DMN_SG are the relevant connectivity values.

%DISCOVERY COHORTS
mdl = fitlm([TBI MDD DepSeverity PTSD Dataset],DAN_SG,’y ~ x1:x2 + x1 + x2 + x3 + x4 + x5 + x1:x3′,’CategoricalVars’,[1 2 4 5]);
mdl = fitlm([TBI MDD DepSeverity PTSD Dataset],DAN_DMN,’y ~ x1:x2 + x1 + x2 + x3 + x4 + x5 + x1:x3′,’CategoricalVars’,[1 2 4 5]);
mdl = fitlm([TBI MDD DepSeverity PTSD Dataset],DMN_SG,’y ~ x1:x2 + x1 + x2 + x3 + x4 + x5 + x1:x3′,’CategoricalVars’,[1 2 4 5]);

%REPLICATION COHORT
mdl = fitlm([TBI MDD DepSeverity PTSD],DAN_SG,’y ~ x1:x2 + x1 + x2 + x3 + x4 + x1:x3′,’CategoricalVars’,[1 2 4]);
mdl = fitlm([TBI MDD DepSeverity PTSD],DAN_DMN,’y ~ x1:x2 + x1 + x2 + x3 + x4 + x1:x3′,’CategoricalVars’,[1 2 4]);
mdl = fitlm([TBI MDD DepSeverity PTSD],DMN_SG,’y ~ x1:x2 + x1 + x2 + x3 + x4 + x1:x3′,’CategoricalVars’,[1 2 4]);

%The main/interaction effects from this model were used to create Figure 2.
%The residuals from the model were used to test the combined effect of
%DAN_SG and DAN_DMN, as depicted in Figure 3.
%The connectivity values were computed using the corrmat_inputROI script,
%see below:


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




function zmat = corrmat(tc,ROIs,mask,within,Zscore)
%Generates correlation matrix from input timecourse and ROIs
%ROIs is an MxN matrix of ROIs, where M is the number of voxels and N is
%the number of ROIs.
%tc = input timecourse
%mask = a mask within which you are working (saves computational resources)
%within = should we compute within-ROI connectivity? (1=yes, 0=no) (this is resource-intensive)
%Zscore = should this be Zscored? (1=yes, 0=no) – if not, output will be Fisher z(r)


numrois = size(ROIs,2);
datapre = tc;
datapre = mask_4dfpimg(datapre,mask,’remove’);

d=size(datapre);

meanimg=mean(datapre,2);
meanimg=repmat(meanimg,[1 d(2)]);
datapre=datapre-meanimg;

for p=1:numrois
clear temproitc avgroitc seed;
seed = ROIs(:,p);
seed(find(isnan(seed)))=0;

%WEIGHTED ROI
[temproitc] = datapre.*seed;
[temproitc] = mask_4dfpimg(temproitc,seed,’remove’);
avgroitc=nanmean(temproitc,1);
tcmat(:,p)=avgroitc’;
end

pre(:,:)=(corrcoef(tcmat));
zmat = FisherTransform(pre);
if within==1
for p=1:numrois
clear temproitc avgroitc1 avgroitc2 seed;
seed = ROIs(:,p);
seed(find(isnan(seed)))=0;
[temproitc] = datapre.*seed;
[temproitc] = mask_4dfpimg(temproitc,seed,’remove’);
s = size(temproitc,1);
if s>10000
temproitc1 = temproitc(1:4:end,:);
temproitc2 = temproitc(2:4:end,:);
temproitc3 = temproitc(3:4:end,:);
temproitc4 = temproitc(4:4:end,:);
s = size(temproitc4,1);
temproitc = temproitc1(1:s,:) + temproitc2(1:s,:) + temproitc3(1:s,:) + temproitc4;
clear temproitc1 temroitc2 temproitc3 temproitc4;
end

tmpz = atanh(corr(temproitc’));
zmat(p,p) = nanmean(nanmean(tmpz(abs(tmpz)<3.6)));
end

clear tcmat;
end