Perceived-wealth-and-difficulty-modulation/Code/neuralAnalyses.m

% Computes behavioral analyses on the dACC dataset.
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0

clear all; close all; clc;

savefolder='./Figures/';

dACC =load('../Data/dACC_all.mat');
% NOTE: CELLS 1-55 subject1, 56-129 subject2
dACC_subject(1:55)=1; dACC_subject(56:length(dACC.data))=2;
Ntrs_dACC=arrayfun(@(cc)  dACC.data{cc}.numTrials, 1:length(dACC.data));
Nccs_dACC=length(dACC.data);

%% Loading main data from vars data structures
Vtl_dACC=reshapeACCstruct2varsmat(dACC,5);
Vbl_dACC=reshapeACCstruct2varsmat(dACC,6);
Vtr_dACC=reshapeACCstruct2varsmat(dACC,7);
Vbr_dACC=reshapeACCstruct2varsmat(dACC,8);
Ptl_dACC=reshapeACCstruct2varsmat(dACC,12)./100;
Ptr_dACC=reshapeACCstruct2varsmat(dACC,13)./100;
Ord_dACC=reshapeACCstruct2varsmat(dACC,15);

start_trial_dACC=reshapeACCstruct2varsmat(dACC,1);
start_outcm_dACC=reshapeACCstruct2varsmat(dACC,18)-start_trial_dACC;
startchoice_dACC=reshapeACCstruct2varsmat(dACC,19)-start_trial_dACC;
valid_trs_dACC=double(startchoice_dACC<10 & start_outcm_dACC<10);
startchoice_dACC(valid_trs_dACC==0)=nan;
start_outcm_dACC(valid_trs_dACC==0)=nan;

%% Computing EV and Risk

EVl_dACC = Ptl_dACC.*Vtl_dACC + (1-Ptl_dACC).*Vbl_dACC;
EVr_dACC = Ptr_dACC.*Vtr_dACC + (1-Ptr_dACC).*Vbr_dACC;

EVf_dACC=nan(size(Ord_dACC));
EVs_dACC=nan(size(Ord_dACC));
EVf_dACC(Ord_dACC==1)=EVl_dACC(Ord_dACC==1);
EVf_dACC(Ord_dACC==2)=EVr_dACC(Ord_dACC==2);
EVs_dACC(Ord_dACC==1)=EVr_dACC(Ord_dACC==1);
EVs_dACC(Ord_dACC==2)=EVl_dACC(Ord_dACC==2);

%% Computing psth and for valid i-1 trials
% tb_m1_p3 is -1s, +3s around first offer onset (not choice/outcome event-locked)
% tb_tlock is -1s, +3s around first offer + time-locked choice, + time-locked outcome

timeAxis=(-5e3+20:20:10e3)/1e3; % in seconds
Ntps=length(timeAxis);

psth_dACC0=nan(Nccs_dACC,max(Ntrs_dACC),Ntps);
psth_dACC_offers=nan(Nccs_dACC,max(Ntrs_dACC),Ntps);
psth_dACC_chs=nan(Nccs_dACC,max(Ntrs_dACC),Ntps);
psth_dACC_ocs=nan(Nccs_dACC,max(Ntrs_dACC),Ntps);

tr_st_bin=find(timeAxis<0,1,'last');

for cc=1:Nccs_dACC
    psth_dACC0(cc,1:size(dACC.data{cc}.psth,1),:)=dACC.data{cc}.psth;
    for tr=1:Ntrs_dACC(cc)
        if valid_trs_dACC(cc,tr) %&& ~isnan(start_trial_dACC(cc,tr))
            ichs_cc_tr=find(timeAxis<=startchoice_dACC(cc,tr),1,'last');
            iche_cc_tr=find(timeAxis>=start_outcm_dACC(cc,tr),1,'first');
            psth_dACC_offers(cc,tr,1:(ichs_cc_tr-tr_st_bin+1000/20))=psth_dACC0(cc,tr,[(-(1000/20)+1:0)+tr_st_bin, tr_st_bin+1:ichs_cc_tr]);
            %psth_dACC_chs(cc,tr,1:(iche_cc_tr-ichs_cc_tr))=psth_dACC0(cc,tr,ichs_cc_tr+1:iche_cc_tr);
            iocs_cc_tr=find(timeAxis<=start_outcm_dACC(cc,tr),1,'last');
            psth_dACC_chs(cc,tr,1:(iocs_cc_tr-ichs_cc_tr))=psth_dACC0(cc,tr,ichs_cc_tr+1:iocs_cc_tr);
            psth_dACC_ocs(cc,tr,1:(Ntps-iocs_cc_tr+1))=psth_dACC0(cc,tr,iocs_cc_tr:end);
        end
    end
end
psth_dACC1=cat(3,rmnans(psth_dACC_offers,3),rmnans(psth_dACC_chs,3),rmnans(psth_dACC_ocs,3));
offs_dACC_Ntps=size(rmnans(psth_dACC_offers,3),3);
chcs_dACC_Ntps=size(rmnans(psth_dACC_chs,3),3);
otcs_dACC_Ntps=size(rmnans(psth_dACC_ocs,3),3);

psth_dACC_ntrav=(sum(~isnan(reshape(psth_dACC1,[Nccs_dACC*max(Ntrs_dACC) size(psth_dACC1,3)])))./sum(valid_trs_dACC(:)));

Ntps_psth_dACC_offs_used=sum(psth_dACC_ntrav(1:offs_dACC_Ntps)>.9);
Ntps_psth_dACC_chcs_used=sum(psth_dACC_ntrav(offs_dACC_Ntps+(1:chcs_dACC_Ntps))>.9);
Ntps_psth_dACC_otcs_used=sum(psth_dACC_ntrav(offs_dACC_Ntps+chcs_dACC_Ntps+(1:otcs_dACC_Ntps))>.9);

psth_dACC=cat(3,psth_dACC_offers(:,:,1:Ntps_psth_dACC_offs_used),psth_dACC_chs(:,:,1:Ntps_psth_dACC_chcs_used),psth_dACC_ocs(:,:,1:Ntps_psth_dACC_otcs_used));

timeAxis_psth_tl=[-1000:20:0 20:20:(Ntps_psth_dACC_offs_used+Ntps_psth_dACC_chcs_used+Ntps_psth_dACC_otcs_used-50-1)*20 ]/1000;
timeAxis_psth_li=cumsum([0 Ntps_psth_dACC_offs_used-50-1 Ntps_psth_dACC_chcs_used Ntps_psth_dACC_otcs_used])*20/1000;

psth_dACC_rn=renormtoprestim(psth_dACC,psth_dACC(:,:,1:(1000/20)),3);

% Safety-check: plot histogram of trials available along time axis
%plot((sum(~isnan(reshape(psth_dACC,[Nccs_dACC*max(Ntrs_dACC) size(psth_dACC,3)]))./sum(valid_trs_dACC(:))));

clearvars psth_dACC1 psth_sgACC1

%% REGRESSING EVs JOINTLY, EVENT-LOCKED

Nshf=100; %This will impact comptational time, consider setting >=100
Dsfactor=1;
EVfs_dACC=cat(3,EVf_dACC,EVs_dACC);
clearvars pvstruct*

if 1% joint regression for EV1 and EV2
    [~, ~, ~, ~, pv12sig, pv12sig_sh]=fitlmspksvars(psth_dACC_rn,EVfs_dACC,Nshf,Dsfactor);
    pvsig_struct.pvdACCsig=pv12sig;
    pvsig_struct.pvdACCsig_sh=pv12sig_sh;
    pvsig_struct.Ntrs=sum(all(isnan(EVfs_dACC),3),2);
    save('./pvsigEV_all_struct_v1.mat','pvsig_struct');
else
    load pvsigEV_all_struct_v1.mat
end

hff=figure('Units','Normalized','Position',[0 0 .5 .33]);
plot(timeAxis_psth_tl,squeeze(mean(pvsig_struct.pvdACCsig(:,:,2)))*100,'r'); hold on;
plot(timeAxis_psth_tl,squeeze(mean(pvsig_struct.pvdACCsig(:,:,3)))*100,'b');
plot(timeAxis_psth_tl,prctile(squeeze(mean(pvsig_struct.pvdACCsig_sh(:,:,2,:)))*100,95,2),'r:');
plot(timeAxis_psth_tl,prctile(squeeze(mean(pvsig_struct.pvdACCsig_sh(:,:,3,:)))*100,95,2),'b:');
ylim([0 40]); xlim(timeAxis_psth_tl([1 end])); yl=get(gca,'YLim'); 
plot(repmat([0 .6 .75 1.35 1.5 1.62 1.74 8.34],2,1),repmat([0 40]',1,8),'k:');
text(0.05,yl(end)*.85,'offer1','rotation',90); text(0.65,yl(end)*.85,'delay1','rotation',90);
text(0.80,yl(end)*.85,'offer2','rotation',90); text(1.40,yl(end)*.85,'delay2','rotation',90); text(8.40,yl(end)*.85,'outcome','rotation',90);
title('Encoding of offer EV_1, EV_2, all trials'); legend({'EV_1','EV_2'},'box','off','Location','NorthEastOutside');
xlabel('time (s)'); ylabel('% significant cells');
hff.Theme='light';
saveas(hff,[savefolder '/Fig2_EVfs_all_dACC_eventlocked.png']);
saveas(hff,[savefolder '/Fig2_EVfs_all_dACC_eventlocked.svg']);


% ---------------------------------SUPPORT----FUNCTIONS----------------------------
%% Removes nans along specified dimension
function [rx,ri]=rmnans(x,dim)
% rx are the non-nan elements
% ri are the non-nan indices

if nargin<2
    dim=1;
end

if isvector(x)
    rx=x;
    ri=1:length(x);
    ri(isnan(x))=[];
    rx(isnan(x))=[];
elseif ismatrix(x)
    rx=x;
    if dim==1
        ri=1:size(x,1);
        ri(all(isnan(x),2))=[];
        rx(all(isnan(x),2),:)=[];
    elseif dim==2
        ri=1:size(x,2);
        ri(all(isnan(x),1))=[];
        rx(:,all(isnan(x),1))=[];
    end
elseif length(size(x))==3
    ri=1:size(x,dim);
    ri=ri(~all(isnan(x),setdiff(1:3,dim)));
    if dim==1
        rx=x(ri,:,:);
    elseif dim==2
        rx=x(:,ri,:);
    elseif dim==3
        rx=x(:,:,ri);
    end
else
    error('more than 3D is not supported');
end
end

%% Baseline-normalization for psth
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0
function Xnorm = renormtoprestim(X,Xprestim,dim)
% !important time has in X and Xprestim to lie along dim dimension

Ndims=numel(size(X));
Ndim=size(X,dim);

%if dim==3 -> [1 1 Ndim]
%if dim==2 -> [1 Ndim 1]

NX=size(X);
NP=size(Xprestim);
elsedims=setdiff(1:Ndims,dim);

if ~all(NX(elsedims)==NP(elsedims))
    error('X and Xprestim need to match size in all dimensions except dim.');
end

repvec=ones(1,Ndims);
repvec(dim)=Ndim;

Xnorm=(X-repmat(nanmean(Xprestim,dim),repvec)); % De-mean

end

%% Compute Linear model fits
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0
function [pvalues, betas, pvalues_shf, betas_shf, pvalues_sig, pvalues_sig_shf, r2values, r2values_shf]=fitlmspksvars(Xspks,Yvars,Nshf,Dsfactor,pvsig_th)

if nargin<5
    if nargin<4
        if nargin<3
            Nshf=0;
        end
        Dsfactor=3;
    end
    pvsig_th=.05;
end
% X is [cells, trials,  time ] the neural signal
% Y is [cells, trials~, nvars] the trial variable

[Nccx, Ntrx, Ntp] = size(Xspks);
[Nccy, Ntry, Nvars] = size(Yvars); %for Nvars=1 only computes beta0 and beta1

assert(Nccx==Nccy, 'The two first variables have to match size along first dimension.');
assert(Ntrx==Ntry, 'The two first variables have to match size along second dimension.');

Ntp=Ntp/Dsfactor;
betas=nan(Nccx,Ntp,Nvars+1);
pvalues=nan(Nccx,Ntp,Nvars+1);
r2values=nan(Nccx,Ntp);
pvalues_shf=nan(Nccx,Ntp,Nvars+1,Nshf);
betas_shf=nan(Nccx,Ntp,Nvars+1,Nshf);
r2values_shf=nan(Nccx,Ntp,Nshf);

warning('off','all')

for cc=1:Nccx
    Xccmm=movmean(squeeze(Xspks(cc,:,:))',[0 10],"omitnan")'; % works well only for DSfactor==1
    Ycc=squeeze(Yvars(cc,:,:));
    if Nvars>1
        nans2rem=all(isnan(Ycc),2);
        Ycc=Ycc(~nans2rem,:);
        Nycc=size(Ycc,1);
    else
        nans2rem=isnan(Ycc);
        Ycc=Ycc(~nans2rem);
        Nycc=length(Ycc);
    end
    Xccmm=Xccmm(~nans2rem,:);
    Nxcm=size(Xccmm);
    parfor tt=1:Nxcm(2)%Ntp
        lastwarn('');
        gfitccttmdl=fitlm(Ycc,Xccmm(:,tt));
        % faster withouth function call
        if isempty(lastwarn())
            pvalues(cc,tt,:)=gfitccttmdl.Coefficients.pValue;
            betas(cc,tt,:)=gfitccttmdl.Coefficients.Estimate;
            r2values(cc,tt)=gfitccttmdl.Rsquared.Ordinary;
        end
        if Nshf>0
            shpvs=nan(Nvars+1,Nshf);
            shbvs=nan(Nvars+1,Nshf);
            shrvs=nan(1,Nshf);
            currx=Xccmm(:,tt);
            for sh=1:Nshf
                lastwarn('');
                if Nvars>1
                    shmdl=fitlm(Ycc(randperm(Nycc),:),currx(randperm(Nxcm(1))));
                else
                    shmdl=fitlm(Ycc(randperm(Nycc)),currx(randperm(Nxcm(1))));
                end
                if isempty(lastwarn())
                    shpvs(:,sh)=shmdl.Coefficients.pValue;
                    shbvs(:,sh)=shmdl.Coefficients.Estimate;
                    shrvs(sh)=shmdl.Rsquared.Ordinary;
                end
            end
            pvalues_shf(cc,tt,:,:)=shpvs;
            betas_shf(cc,tt,:,:)=shbvs;
            r2values_shf(cc,tt,:)=shrvs;
        end
        %gfitccttmdl=[];
    end
    disp([sprintf('%s',datetime('now','Format','HH:mm')) ' - done: ' num2str(cc) '/' num2str(Nccx) ', ' sprintf('%2.2f',cc/Nccx.*100) '%']);
end

%Safety checks: plot(sum(squeeze(pvalues(:,:,2))<0.05,1))
pvalues_sig=uint8(nanlt(pvalues,pvsig_th)); %uint8 reduces storage overhead
pvalues_sig_shf=uint8(nanlt(pvalues_shf,pvsig_th));
warning('on','all')
end

%% Compare elements in A to be lower than B but keeping NaNs in place
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0
function AltB=nanlt(A,B)
if(isscalar(B) && isnan(B))
    error('nan is not allowed as input 2');
end

if  isscalar(B)
    nanmask=zeros(size(A));
    nanmask(isnan(A))=nan;
elseif all(size(A)==size(B))
    nanmask=zeros(size(A));
    nanmask(isnan(A))=nan;
else
    error('check inputs size, input 2 scalar is allowed.');
end
AltB=sum(cat(length(size(A))+1,double(A<B),nanmask),length(size(A))+1);
end


%% Get task variables
function ACCmat=struct2varsmat(ACCstruct, vars_index, start_offset,end_offset)
% The function takes as first input ACCstruct loaded from "../Data/dACC_all.mat" and provides task variables,
% respectively indicized by the second input vars_index, in trial range start_offset:end_offset
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0
if nargin<4
    if nargin<3
        start_offset=0;
    end
    end_offset=0;
end
Nccs=length(ACCstruct.data);
Ntrs=arrayfun(@(cc) ACCstruct.data{cc}.numTrials, 1:Nccs);

ACCmat=nan(Nccs,max(Ntrs));

for cc=1:Nccs
    Ntrcc=size(ACCstruct.data{cc}.vars,1);
    ACCmat(cc,:)=cat(1,ACCstruct.data{cc}.vars(:,vars_index),nan(max(Ntrs)-Ntrcc,1));
    ACCmat(cc,1:start_offset)=nan;
    ACCmat(cc,Ntrcc-end_offset+1:Ntrcc)=nan;
end
end