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

Original code as archived. Formatting and dependencies have not been revalidated. The archive contains typographic quote characters.

%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