Dynamics of Intersubject Brain Networks during Anxious Anticipation

15 Pages·2017·4.16 MB·English

by  Luiz Pessoa

ORIGINALRESEARCH
published:21November2017
doi:10.3389/fnhum.2017.00552
Dynamics of Intersubject Brain
Networks during Anxious
Anticipation
MahshidNajafi,JoshuaKinnisonandLuizPessoa*
DepartmentofPsychologyandMarylandNeuroimagingCenter,UniversityofMaryland,CollegePark,CollegePark,MD,
UnitedStates
How do large-scale brain networks reorganize during the waxing and waning of
anxious anticipation? Here, threat was dynamically modulated during human functional
MRI as two circles slowly meandered on the screen; if they touched, an unpleasant
shock was delivered. We employed intersubject correlation analysis, which allowed
the investigation of network-level functional connectivity across brains, and sought to
determinehownetworkconnectivitychangedduringperiodsofapproach(circlesmoving
closer) and periods of retreat (circles moving apart). Analysis of positive connection
weights revealed that dynamic threat altered connectivity within and between the
salience, executive, and task-negative networks. For example, dynamic functional
connectivity increased within the salience network during approach and decreased
during retreat. The opposite pattern was found for the functional connectivity between
the salience and task-negative networks: decreases during approach and increases
during approach.Functionalconnectionsbetweensubcorticalregions andthesalience
network also changed dynamically during approach and retreat periods. Subcortical
regions exhibiting such changes included the putative periaqueductal gray, putative
Editedby:
habenula,andputativebednucleusofthestriaterminalis.Additionalanalysisofnegative
LucinaQ.Uddin,
UniversityofMiami,UnitedStates functional connections revealed dynamic changes, too. For example, negative weights
Reviewedby: within the salience network decreased during approach and increased during retreat,
ChristinaB.Young,
oppositewhatwasfoundforpositiveweights.Together,ourfindingsunraveleddynamic
NorthwesternUniversity,UnitedStates
XinDi, featuresoffunctionalconnectivityoflarge-scalenetworksandsubcorticalregionsacross
NewJerseyInstituteofTechnology, participants while threat levels varied continuously, and demonstrate the potential of
UnitedStates
characterizingemotionalprocessingatthelevelofdynamicnetworks.
AaronShainHeller,
UniversityofMiami,UnitedStates
Keywords: emotion, networks, intersubject correlation, dynamics, threat, amygdala, bed nucleus of the stria
*Correspondence:
terminalis(BST),periaqueductalgray(PAG)
LuizPessoa
[email protected]
INTRODUCTION
Received:28July2017
Accepted:31October2017
Imagineyourselfrecliningonadentist’schair.Mostofuswaitanxiouslyasthedentistgradually
Published:21November2017
movesthedrilltowardourmouth.Atthesametime,ifthedrillismovedaway(perhapsthedentist
Citation: neededanadditionaladjustment),anxiousapprehensionlikelywillsubside.Agrowingliterature
NajafiM,KinnisonJandPessoaL
ofbothnon-humanandhumanresearchindicatesthatanticipatoryprocessingofnegativeevents
(2017)DynamicsofIntersubjectBrain
engages multiple brain regions (Davis et al., 2010; Grupe and Nitschke, 2013; Tovote et al.,
NetworksduringAnxiousAnticipation.
Front.Hum.Neurosci.11:552. 2015),includingmedialprefrontalcortex,insula,andorbitofrontalcortex,cortically.Subcortically,
doi:10.3389/fnhum.2017.00552 implicatedregionsincludetheamygdala,periaqueductalgray(PAG),andthebednucleusofthe
FrontiersinHumanNeuroscience|www.frontiersin.org 1 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
stria terminalis (BST) (see Davis et al., 2010; Fox et al., 2015).
Anticipatory negative processing allows participants to prepare
andpossiblyminimizetheimpactofharmfulstimuli.However,
aberrant responding to uncertain future negative events is
believedtobecentraltoanxietydisorders(GrupeandNitschke,
2013; Fox and Kalin, 2014). Thus, further elucidation of the
mechanisms of anticipatory processing is important from both
basicandclinicalperspectives.
Anxiousanticipationofnegativeeventshaswidespreadeffects
on brain function (Thomason et al., 2011; Kang et al., 2016;
Raz et al., 2016; Young et al., 2017). However, understanding
how the organization of large-scale brain networks is affected
during anxious apprehension is poorly understood. At least
three networks are altered by processing threat (Pessoa, 2013):
a salience network that responds to motivationally salient
stimuli (Seeley et al., 2007; Menon and Uddin, 2010); a task-
negative (or “default mode”) network that is engaged when
attention is directed internally and during some forms of
emotional processing (Gusnard et al., 2001; Greicius et al.,
2003, 2009); and an executive control network that is engaged
when cognitively demanding tasks require attention (Seeley
et al., 2007; Vanhaudenhuyse et al., 2011). In the context
FIGURE1|Experimentalparadigm.Tovarythreatleveldynamically,two
of anxious anticipation, one study described greater salience- circlesmovedaroundonthescreenwithsomedegreeofrandomness,
network connectivity while participants watched an aversive sometimesapproachingeachother,othertimesretreatingfromeachother.
movie(Hermansetal.,2011).Inapreviousstudy,weinvestigated Whentheycollidedwitheachother,anunpleasantmildelectricshockwas
delivered.
networkinteractionswhenparticipantswereineitherprolonged
threat(unpredictablemildshockscouldbeadministered)orsafe
(noshockspossible)conditions,andcharacterizedtransientand
sustained changes to the three networks above (McMenamin used to compute correlations between the same region (across
etal.,2014). brains) as well as correlations between different regions (again
Anxious anticipation is inherently temporal. Although acrossbrains)(Figure2B)(NajafiandPessoa,2016;Simonyetal.,
previous studies have investigated how brain responses are 2016). Simony et al. (2016) showed that this method increased
sensitivetothreatproximity(Mobbsetal.,2010;Somervilleetal., the signal-to-noise ratio in detecting functional correlations
2010; Grupe et al., 2013), little is known about how patterns (comparedtocomputingthemwithinbrains)likelyfromfiltering
of brain co-activation (thus networks) change during dynamic outprocessingunrelatedtoongoingstimulusprocessing,aswell
manipulations of threat. To address this gap in the literature, as non-neuronal artifacts (for example, respiratory rate, head
here we modulated threat dynamically during functional MRI motion) that can influence correlation patterns within a brain
scanning.Twocirclesmovedonthescreen,sometimesmoving butaretypicallynotcorrelatedacrossbrains.Anotherimportant
closerandsometimesmovingapart(Figure1).Iftheytouched, propertyofintersubjectnetworkanalysisisthatitcanconsider
an unpleasant shock was delivered to the participant. We thecorrelationofaregionwithitself;inintersubjectanalysisthis
soughttodeterminehowfunctionalconnectivitychangedduring correlationismeaningfulbecausethetimeseriesdatacomefrom
periodsofapproach(circlesmovingcloser)andperiodsofretreat differentbrains.
(circles moving apart). As in our previous study (McMenamin
et al., 2014), we studied a set of regions spanning the salience, MATERIALS AND METHODS
executive, and task-negative networks, given their involvement
incognitiveandemotionalprocessing(Yeoetal.,2011;Pessoa, Participants
2013).Inaddition,weinvestigatedsubcorticalregionsimportant Ninety-three participants with normal or corrected-to-normal
foremotionalprocessing. visionandnoreportedneurologicalorpsychiatricdiseasewere
Overall, our approach allowed us to test several questions recruited from the University of Maryland community. Data
aboutthebrainbasisofanxiousanticipation.Howdofunctional from84participants(44malesand40females,ages18–40years;
connectivity properties of large-scale networks evolve during average:22.62,STD:4.85)wereemployedfordataanalysis(ofthe
periodsofthreatapproachandretreat?Duringdynamicthreat, original sample of 93, data from 7 subjects were discarded due
whatistherelationshipbetweencorticalandsubcorticalregions to technical issues during data transfer [specifically, field maps
important for threat processing? We investigated functional were lost], 1 subject was removed because of poor structural-
connectivity based on intersubject correlation analysis (Hasson functionalalignment,and1subject’sdatawerelost).Theproject
et al., 2004), where time series data from voxels/ROIs are was approved by the University of Maryland College Park
correlatedacrossparticipants(Figure2A).Thisapproachcanbe InstitutionalReviewBoardandallparticipantsprovidedwritten
FrontiersinHumanNeuroscience|www.frontiersin.org 2 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
FIGURE2|Intersubjecttimeseriesanalysis.(A1)Instandardintersubjectanalysis,correlationbetweenthesameregionacrossdifferentparticipants’brainsis
calculated.(A2)Tocalculateintersubjectcorrelation,foreachvoxelorregionofinterest(ROI),thetimeseriesforonesubject(say,S1)iscorrelatedwiththeaverage
timeseriesofallothersubjects(S2,…SN),forthesamevoxel/ROI.Thisprocessistheniteratedandaveragedtodeterminegroup-levelcorrelations,namelyavector
thatcontainsthecorrelationofeveryvoxel/ROIwithitself(acrossparticipants)(A3).(B1)Themethodcanbegeneralizedtostudymultipleregionsbycomputingthe
correlationsbetweenallpairsofvoxels/ROIsacrossdifferentbrains.(B2)Foreachvoxel/ROI,thetimeseriesofaonesubject(say,S1)iscorrelatedwiththeaverage
timeseriesacrossothersubjects(S2,…,SN).Thisprocessistheniteratedandaveragedtodeterminegroup-levelcorrelations,namelyamatrixthatcontainsthe
correlationofeveryvoxel/ROIwitheveryothervoxel/ROI(acrossparticipants)(B3).NotethatthevectorinA3correspondstothediagonalofthematrixinB3,
illustratingthatintersubjectnetworksprovideamoregeneralcharacterizationoftimeseriesrelationships.
informed consent before participation. Results are reported for wherec=0.2andN(0,1)indicatesthenormaldistributionwith
a group of 49 participants (23 males and 26 females, ages 0 mean and standard deviation 1. The position and amount of
18–40years;average:22.78,STD:5.40).Aseparategroupof35 displacementofeachcirclewereupdatedindependently.
participants was used as an exploratory dataset to fix specific Visual stimuli were presented using PsychoPy (http://www.
processing choices and to define subcortical regions of interest psychopy.org/) and viewed on a projection screen via a mirror
(ROIs). mounted to the scanner’s head coil. The total experiment
included6runs,eachofwhichhad6blocks.Ineachblock,the
circlesappearedonthescreenandmovedaroundfor60s;blocks
Procedure and Stimuli
wereprecededbya15-sblankscreen.Eachrunendedwith7sof
Two circles with different colors moved around on the screen
ablankscreen.
randomly, and when they collided an unpleasant mild electric
In each of the 6 runs the circles collided a total of 8 times
shockwasdelivered.Overall,theproximityandrelativevelocity
in 4 out of the 6 blocks (3 shocks maximum per block); each
ofthecircleswereusedtoinfluencethreatlevel.Thepositionof
collisionresultedinthedeliveryofanelectricshock.The500-ms
eachcircle(ontheplane),x,wasdefinedbasedonitsprevious
t electricshockwasdeliveredbyanelectricstimulator(Coulbourn
position, xt−1, plus a random displacement, xt. The magnitude Instruments,PA,USA)tothefourthandfifthfingersofthenon-
anddirectionofthedisplacementwascalculatedbycombining
dominantlefthandviaMRI-compatibleelectrodes.Tocalibrate
anormalrandomdistributionwithamomentumtermtoensure
theintensityoftheshock,eachparticipantwasaskedtochoose
motionsmoothness,whileatthesametimeremaining(relatively)
his/her own stimulation level immediately prior to functional
unpredictabletoparticipants.Specifically,thedisplacementwas
imaging, such that the stimulus would be “highly unpleasant
updatedevery50msasfollows:
but not painful.” After each run, participants were asked about
theunpleasantnessofthestimulusinordertore-calibrateshock
1xt =(1−c)1xt−1+ cN(0,1) strength,ifneeded.
FrontiersinHumanNeuroscience|www.frontiersin.org 3 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
MRI Data Acquisition in the brain-to-skull segmentation can be easily identified.
Functional and structural MRI data were acquired using a 3T Accordingly, to skull strip the T1 high-resolution anatomical
Siemens TRIO scanner with a 32-channel head coil. First, a image (which was rotated to match the oblique plane of
high-resolution T2-weighted anatomical scan using Siemens’s the functional data with AFNI’s 3dWarp), we employed six
SPACEsequence(0.8mmisotropic)wascollected.Subsequently, differentpackages:ANTs(Avantsetal.,2011),AFNI(Cox,1996;
we collected 457 functional EPI volumes using a multiband http://afni.nimh.nih.gov/),ROBEX(Iglesiasetal.,2011;https://
scanning sequence (Feinberg et al., 2010) with TR = 1.0s, www.nitrc.org/projects/robex),FSL(http://fsl.fmrib.ox.ac.uk/fsl/
TE = 39ms, FOV = 210mm, and multiband factor = 6. Each fslwiki),SPM(http://www.fil.ion.ucl.ac.uk/spm),andBrainsuite
volume contained 66 non-overlapping oblique slices oriented (ShattuckandLeahy,2002).Weemployeda“votingscheme”as
30◦Clockwise relative to the AC-PC axis (2.2mm isotropic). In follows:basedonT1data,avoxelwasconsideredtobepartofthe
addition, a high-resolution T1-weighted MPRAGE anatomical brainif4/6packagesestimatedittobeabrainvoxel;otherwise
scan(0.8mmisotropic)wascollected.Finally,double-echofield the voxel was not considered to be brain tissue (for 6 subjects
maps (TE1 = 4.92ms, TE2 = 7.38ms) were acquired with whoseT1datawerelostduetoissuesduringdatatransfer,theT2
acquisitionparametersmatchedtothefunctionaldata. imagewasusedinsteadandonlytheANTspackagewasusedfor
skull-stripping).
Exploratory and Test Data Sets Subsequently,FSLwasusedtoprocessfieldmapimagesand
Data from 84 participants were employed in this study. One create a phase-distortion map for each participant (with bet
of our goals was to attempt to enhance reproducibility in a andfsl_prepare_fieldmap).FSL’sepi_regwasthenusedtoapply
researchareathatfacesthe“curseofflexibility.”Forexample,a boundary-based co-registration to align the unwarped mean
recent review enumerated 69,120 different workflows for basic volumeregisteredEPIimageswiththeskull-strippedanatomical
functional MRI analysis alone (Poldrack et al., 2017). We thus image (T1 or T2) along with simultaneous EPI distortion-
employed an exploratory data set to fix specific processing correction (Greve and Fischl, 2009). Next, ANTS was used to
choices, and to define subcortical ROIs, as described in this determine a non-linear transformation that mapped the skull-
section. At a point during data collection when approximately stripped anatomical image (T1 or T2) to the MNI152 template
40 participants had been studied, we labeled the data as (interpolatedto1-mmisotropicvoxels).Finally,ANTscombined
“exploratory” and to be used for data exploration. A total of N the non-linear transformations from co-registration/unwarping
=35usableparticipantswereusedintheexploratoryset.With (from mapping mean functional EPI images to the anatomical
the entire processing pipeline and ROIs fixed, statistical testing T1 or T2) and normalization (from mapping T1 or T2 to the
wasthenappliedtoaseparatedataset(N =49;theoriginalgoal MNItemplate)intoasingletransformationthatwasappliedto
beingtohaveapproximately50participants). mapregisteredfunctionalvolumesoffunctionaldatatostandard
space (interpolated to 2-mm isotropic voxels). In this process,
Stimulus Conditions ANTsalsoutilizedthefieldmapstosimultaneouslyminimizeEPI
The exploratory set was used to define two conditions, distortion.
“approach” and “retreat,” based on whether the circles were Additional preprocessing steps included the following. The
movingtowardorawayfromeachother.Timepointswereonly first three volumes of each functional run were discarded to
considered for analysis if the Euclidian distance between the account for equilibration effects. Slice-timing correction (with
circleswasatleast75%ofthemaximumdistancethatthecircles AFNI’s 3dTshift) used Fourier interpolation to align the onset
exhibitedduringthewholeexperiment;otherwise,thedatawere times of every slice in a volume to the first acquisition slice,
notemployedintheanalysis.Therationalebehindthiswasthat, andthenasix-parameterrigidbodytransformation(withAFNI’s
whenthecircleswerefarfromeachother,participantsreported 3dvolreg) corrected head motion within and between runs by
thattheydidnotreallypayasmuchattentiontothem.Therefore, spatiallyregisteringeachvolumetothefirstvolume.
we reasoned that the analysis should focus on the time points
during which the circles were in (relative) closer proximity to
eachother.Investigationoftheexploratorysetrevealedthatthe Subcortical ROIs
particularcutoffwasnotcriticalfortheeffectsinvestigatedand We used the exploratory data set to investigate subcortical
thatvaluesatleastbetween65and85%wereadequate;wechose ROIs focusing on functional to anatomical co-registration.
a cutoff value of 75%. After applying the cutoff, the approach The subcortical ROIs selected for evaluation in the test set
condition included 481 data points (i.e., TRs) and the retreat included the amygdala, PAG, habenula, and BST. For the
condition290datapoints. amygdala,weconsideredtwosubregionsdefinedintheNacewicz
et al. (2014) parcellation, specifically: lateral amygdala and
Functional MRI Preprocessing central/medial amygdala. For the PAG, we modified the mask
We employed the exploratory data set to define the following by Roy et al. (2014), which was dilated by 1 voxel; in addition,
preprocessingsteps.Skullstrippingdetermineswhichvoxelsare we manually removed voxels overlapping cerebrospinal fluid.
to be considered part of the brain and, although conceptually The habenula has been implicated in emotional/motivational
simple, plays a very important role in successful subsequent processing (Hikosaka, 2010; Mizumori and Baker, 2017), and
co-registration and normalization steps. Currently, available hereweemployedamaskdefinedaccordingtotheMorelatlas,
packages perform sub-optimally in specific cases, and mistakes as defined in Krauth et al. (2010). For the BST, we employed a
FrontiersinHumanNeuroscience|www.frontiersin.org 4 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
recentlydevelopedmaskbasedon7Tesladatabutdefinedhaving TABLE1|ListofcorticalandsubcorticalRegionsofInterest(ROIs).
inmind3Teslaacquisition(Theissetal.,2017).
ROInames Coordinates(MNI)
x y z
ROIs for Large-Scale Networks
Weinvestigatedthreenetworkswidelystudiedintheliterature: SALIENCE
salience, executive, and task negative. From these networks, we 1) FrontoinsularcortexL −34 18 4
employedcorticalROIs(definedas5-mmradiusspheres)based 2) FrontoinsularcortexR 34 22 4
onthecentercoordinatesprovidedbypreviousstudies(Table1): 3) Dorsalanteriorcingulatecortex 2 10 40
salience network (Hermans et al., 2011) (13 regions), executive
4) Temporo-parietaljunctionL −62 −26 36
network (Seeley et al., 2007) (12 regions), and task-negative
5) Temporo-parietaljunctionR 62 −26 36
network (Fox et al., 2005) (12 regions). If two ROIs abutted
6) InferotemporalcortexL −54 −62 −4
eachother,eachmaskwaserodedby1voxelfromthetouching
7) InferotemporalcortexR 54 −54 −8
boundarytominimizeanypotentialdata“spillover.”
8) PrecentralL −26 −6 64
9) PrecentralR 26 −2 64
Intersubject Correlation Analysis 10) DorsolateralprefrontalcortexL −38 42 24
We investigated functional connections based on intersubject 11) DorsolateralprefrontalcortexR 34 46 28
correlationanalysis(Hassonetal.,2004).Thiswasthefirststep 12) InferiorfrontalgyrusL −54 6 20
fromwhichtemporaldynamicswasevaluated(seebelow). 13) InferiorfrontalgyrusR 54 10 12
Toperformintersubjectcorrelationanalysis,timeseriesdata EXECUTIVE
from ROIs were correlated across participants, that is, across 14) OrbitalfrontoinsulaL −36 24 −10
different brains (Figure2A). A simple, yet powerful extension 15) DorsolateralprefrontalcortexR 46 46 14
to intersubject correlation analysis is to consider intersubject 16) DorsolateralprefrontalcortexL −34 46 6
correlationsacrossallpairsofROIs,whichallowstheapplication 17) VentrolateralprefrontalcortexR 34 56 −6
of the technique to networks (Figure2B) (Najafi and Pessoa, 18) VentrolateralprefrontalcortexL −32 54 −4
2016; Simony et al., 2016). The procedure to generate an
19) FrontaloperculumR 56 14 14
intersubject network is as follows (Figure3). For a given ROI,
20) Dorsolateralprefrontalcortex/frontal 30 12 60
onesubject’ssdataareheldout(ys),andtherestofthesubjects’ eyefieldR
time series is averaged (y−s). Then, the Pearson correlation 21) Dorsolateralprefrontalcortex/frontal −32 18 50
betweenthesetwodatavectorsiscomputed:corr(y,y ).This eyefieldL
s −s
basic operation is repeated for all pairs of ROIs to compute an 22) Dorsomedialprefrontalcortex 0 36 46
intersubject matrix for the held-out subject (ISMS). Thus, the 23) LateralparietalR 38 −56 44
ISMSisanNxNmatrix,whereNisthenumberofROIs,andthe 24) LateralparietalL −48 −48 48
ij-thmatrixelementcontainsthecorrelationcoefficientbetween 25) InferiortemporalR 58 −54 −16
thei-thROItimeseriesoftheheld-outsubjectandthej-thROI TASKNEGATIVE
timeseriesaveragedacrosstheremainingsubjects. 26) Posteriorcingulateprecuneus(PCC) −3 −38 38
This procedure is repeated for all subjects. A group matrix 27) Retro-splenial 2 −52 9
(ISMG) is then obtained by averaging across all subjects; the 28) Lateralparietalcortex(LP)L −50 −64 38
distribution of intersubject correlations was well approximated 29) Lateralparietalcortex(LP)R 50 −64 38
byaGaussian(mean:−0.0022;STD:0.0181)soaveragingthem 30) Medialprefrontalcortex(MPF)L −4 42 −9
isreasonable.Notethattheresultingintersubjectnetworkisnot 31) Medialprefrontalcortex(MPF)R 0 59 16
necessarily symmetric, because, for each ROI, the time series 32) SuperiorfrontalL −16 44 51
in the held-out subject (y) is not necessarily equal to average
s 33) SuperiorfrontalR 17 43 51
of all other subjects’ time series (y−s) (symmetry is desirable 34) InferiortemporalL −62 −33 −20
because the intersubject correlation between regions i and j
35) InferiortemporalR 66 −18 −20
should, in principle, be equal to the intersubject correlation
36) ParahippocampalgyrusL −22 −26 −20
betweenregionsjandi).WhiletheISM matrixinpracticewill
G 37) ParahippocampalgyrusR 25 −26 −18
benear-symmetric,asimpleandintuitivewaytomathematically
SUBCORTICALREGIONS
accomplish symmetry is to average the group-level intersubject
38) AmygdalaCentral-MedialL
networkwithitstranspose(whererowandcolumnindexesare
39) AmygdalaCentral-MedialR
flipped),leadingtoafinalsymmetricgroupmatrix.Finally,the
40) AmygdalaLateralL
procedures above were performed, separately, for the approach
41) AmygdalaLateralR
and retreat conditions (generating one matrix per condition).
42) Periaqueductalgray(PAG)L
Note that in our method the matrix diagonal is also computed
43) Periaqueductalgray(PAG)R
because data at a given diagonal entry [i, i] is computed across
differentbrains.
(Continued)
FrontiersinHumanNeuroscience|www.frontiersin.org 5 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
TABLE1|Continued pairs each time t (Figure4C), thus computing the matrices
ISMt=0,ISMt=1,etc.(Figure4D).Weconsideredintersubject
ROInames Coordinates(MNI) matrices for t = 0,...,6 s for the approach condition and
x y z t = 0,...,5 s for the retreat condition. For both conditions,
atleast20datapoints(Figure4B)wereavailablepercondition
44) HabenulaL
and time slice (but note that less data were available for the
45) HabenulaR
retreat compared to the approach condition, as some points
46) Bednucleusofthestriaterminalis
werediscardedfollowingshockevents).Whereas,longerperiods
(BST)L
of approach/retreat occurred, they did not occur at least 20
47) Bednucleusofthestriaterminalis
times, which was the minimum number of repetitions that we
(BST)R
judged needed for stable assessment of the correlations. Given
CorticalROIsweredefinedvia5-mmradiusspherescenteredontheMNIcoordinates thisconstraint,therewere443datapointsforapproachand255
providedbelow.SubcorticalROIsweredefinedanatomically(seetextfordetails).
forretreat.
Time Series Data Within- and Between-Network Weight
Because intersubject analysis seeks to investigate correlations To measure the strength of connections within and between
across brains of different individuals, the contributions of the networks,weutilizedwithin-andbetween-networkconnectivity
paradigm were not regressed out as done in some within- weights, respectively. The weight between network N and
1
subject analyses. However, because the moving circles were networkN ,W ,wasdefinedasfollows:
2 N1,N2
presented in 60-s blocks that alternated with 15-s rest period,
we regressed out (with AFNI’s 3dDeconvolve) the block effect 1
(whichwasconvolvedwiththestandardhemodynamicresponse; WN1,N2 = n n X cij (1)
(see Cohen, 1997). Other regressors in the model included 1 2 i∈N1,j∈N2
6 motion parameters (3 linear displacements and 3 angular
where n is the number of ROIs in network N. The value c
i i ij
rotations),andtheirdiscretetemporalderivatives.Additionally,
is the ISM value when the i-th ROI belongs to network N
G 1
to model baseline and drifts of the MR signal, regressors were
and the j-th ROI belongs to network N . When N and N
2 1 2
included corresponding to polynomial terms up to 4th order
are the same network, the formula calculates within-network
(for each run separately). To minimize the shock effect, data
weight. For most of the analyses reported here, we considered
points in a 15-s window after shock delivery were discarded
only positive weights, as commonly done in network research
from all analyses. In addition, to minimize the impact of
(butseebelowfornegativeweights).Mostnetworkmeasuresdo
potential block onsets (which can elicit strong transients), for
nothandleself-connections(Newman,2010),whichinstandard
each block, the first 15 time points (15s) were discarded.
analysis(whencorrelationsaredeterminedwithinparticipants)
Finally, the residual time series for each run was z-scored are not informative (c = 1). Here, we considered functional
ii
separately. connections between the same region (across participants),
The residual time series as defined was used for the which could be incorporated in within-network weight by
intersubjectcorrelationanalysis.Asapproachandretreatvaried consideringthetermsc inEquation(1).
ii
dynamically throughout the blocks, we employed a windowing To study dynamic changes to network cohesion, ISM was
G
procedure to extract data segments corresponding to approach computed as outlined above at each time t and within- and
and retreat periods. Intuitively, the windowing allowed us to between-networkweightscomputed.Linearregressionwasused
selectsegmentsofthetimeseriesassociatedwitheachcondition to evaluate changes in network weight during approach and
andconcatenatethemacrossallruns,asoutlinednext. retreat conditions. To evaluate functional connections between
subcortical regions and the salience network, we computed a
Dynamic Intersubject Correlation weightindexthatsummedallfunctionalconnectionsbetweena
Toinvestigatedynamicaspectsoffunctionalconnectivity,wefirst specificsubcorticalregionandallnodesofthesaliencenetwork.
computed intersubject correlations as indicated previously (see This was performed for the approach and retreat conditions,
Figure3) at each time t. The method considered all approach separately.
and retreat segments separately, and computed one vector for Because intersubject correlations are not independent (each
eachtimet,ROI,andcondition(Figure4).Specifically,foreach participant is present in more than one participant pair),
segment type (approach and retreat, separately), we considered statistical tests were performed via non-parametric resampling
all time series data at t = 0 time points, then t = 1 time tests based on the approach proposed by Kauppi et al.
points, and so on, separately (Figure4A). The goal was to (2014). Briefly, we computed a null resampling distribution by
generate a vector of data at t = 0 by concatenating all of circularly shifting each participants’ time-series by a random
the t = 0 data across approach/retreat segments. To do so, amount so that they were no longer aligned in time across
we concatenated the t = k points (for a fixed k) across participants, and then calculated intersubject correlations.
segments (Figure4B). To account for the hemodynamic delay, The null distribution was estimated by computing 100,000
we discarded the first 5s of each segment. The overall process realizations (fewer iterations produced very similar results
allowed us to compute correlations across subjects for ROI showing that the estimation was stable). Statistically, we were
FrontiersinHumanNeuroscience|www.frontiersin.org 6 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
FIGURE3|Computationofgroup-levelintersubjectcorrelationmatrix(ISM).TheoperatorcorrcorrespondstoPearsoncorrelation.
FIGURE4|Evaluationofintersubjectfunctionalconnectivitydynamics.Intersubjectcorrelationmatriceswerecomputedateachtimepointofapproachandretreat
segments,separately.(A)Timewasusedtoindexdataduringapproachandretreat(notshown)segments.Toaccountforthehemodynamicdelay,wediscardedthe
first5sofeachsegment(graypart).(B)Atimeseriesvectoratt=0wasgeneratedbyconcatenatingallofthet=0samples(acrosssame-conditionsegments
acrossblocksandruns);likewise,fort=1,andsoon.(C)Atagiventimet,timeseriesvectorswereusedtocomputetheintersubjectcorrelationmatrix(D).This
processisshownfort=0andrepeatedfortimest=1,…,6forapproachandt=1,…5forretreat.
interested in evaluating changes to network weight during values (which were due to the unequal number of time
approach and retreat, separately, and comparing the slopes of points used for segments of different duration; for example,
thelinearfitsduringapproachvs.retreat.Thus,observedslopes there were fewer segments of duration 5s than 4s). Thus,
(based on data) were compared to values of the corresponding the mean values expected by chance (determined via the
null distribution and p-values determined; for the contrast of null distribution) were subtracted from the observed values
approach vs. retreat, the observed difference (based on data) (Figures5,7).
was compared to values of the corresponding difference null To explore dynamic changes in negative weights, we
distribution. considered equation (1) by excluding all positive weights and
We observed that the mean of the null distributions followed the same procedures described above for positive
for approach and retreat were not zero but small positive weights.Thenulldistributionwasbasedon100,000realizations.
FrontiersinHumanNeuroscience|www.frontiersin.org 7 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
FIGURE5|Temporalevolutionofintersubjectnetworkweightduringapproachandretreatsegments(positiveweights).Within-andbetween-networkweightsduring
approachandretreatforthesalience,executive,andtask-negativenetworks.Forexample,asthecirclesapproachedeachother,theweightswithinthesalience
networkincreased;whenthecirclesretreated,within-networkweightsdecreased.Slightlylongerapproachperiodswereavailabletocomputedynamicsforthis
condition.Theorangelineshowsapproachweights(90%confidenceband);thecyanlineshowsretreatweights(90%confidenceband).Thered/bluelineshowthe
least-squareslinearfittoweightvaluesduringapproach/retreat;solidlinesindicatefitssuchthatp<0.05.Timeisinseconds;they-axisshowswithin/between
weight.Theasteriskindicatesthattheslopedifferencewassuchthatp<0.05.
FrontiersinHumanNeuroscience|www.frontiersin.org 8 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
FIGURE6|Snapshotsofintersubjectcorrelationmatrices.Resultsareshownforapproachandretreatconditions,aswellasapproachminusretreat.Thetoprow
showsintersubjectcorrelationsatt=1s,andthebottomrowshowsintersubjectcorrelationsatt=5s(inbothcases,ahemodynamiclagof5swasemployed).
RESULTS correlation matrix of the three large-scale networks (negative
valueswereremovedfromeachparticipant’smatrix).
We performed intersubject correlation analysis during threat We also investigated the evolution of the functional
approach and retreat. We focused on regions of the salience, interactionsbetweentargetedsubcorticalregionsandthesalience
executive, and task-negative networks, as well as a targeted set network. In most cases, their correlations with the salience
of subcortical ROIs. As threat level was varied dynamically, networkweredynamic(Figure7).ThePAG,habenula,andBST
weinvestigatedhowtheintersubjectcorrelationmatrixevolved exhibited similar patterns of dynamic intersubject correlations
temporally. Figure5 shows the temporal evolution of within- with the salience network, which increased during approach
and between-network weight during approach and retreat for anddecreasedduringretreat(fortherightBST,changesduring
the three networks. Statistically, we evaluated changes during retreat were not detected). Figure8 illustrates changesbetween
approach and retreat via linear regression, as well as the therightPAGandtheentiresetofsaliencenetworkROIs.Forthe
differencebetweenthetwo;foreaseofreference,statisticalvalues amygdala regions, when changes were detected, only decreases
areprovidedinthefigure(slopesatp<0.05areindicatedbysolid in correlation were observed during both approach and retreat
lines,andp-valuesforthedifferenceinslopesareprovided). segments,andtheslopeswerecomparableinbothconditions.
Evaluation of within-network weight revealed that all
networks exhibited functional connectivity changes during Negative Weights
approach and/or retreat periods, with the strongest change The results above did not consider negative weights. Whereas,
observed for the salience network. During retreat, weight thestepofthresholdingnegativeweightsiscommonlyadopted
decreased within all networks. Changes in intersubject in analysis of brain networks, investigation of changes in
correlation were also found between networks; specifically, negative weights is likely to be informative. For example,
between the salience and executive networks, and between the negative correlations between two systems might indicate that
salience and task-negative networks. The weight between the they are complementary or that they work in some form of
salience and executive networks increased during approach opposition (such as appetitive and aversive systems). More
anddecreasedduringretreat.However,theoppositewasfound generally, effective ways of handling negative ways remain an
betweenthesalienceandtask-negativenetworks.Inthiscase,the openquestion(RubinovandSporns,2011;Fornitoetal.,2016),
correlation actually decreased during approach. For illustration aspositiveandnegativeweightsmayindexqualitativelydifferent
purposes only, Figure6 shows changes of the group-averaged processes.Possibleapproachesincludeconsideringtheabsolute
FrontiersinHumanNeuroscience|www.frontiersin.org 9 November2017|Volume11|Article552
Najafietal. DynamicsofIntersubjectBrainNetworks
FIGURE7|Temporalevolutionofintersubjectweightsbetweensubcorticalregionsandthesaliencenetwork(positiveweights).Dynamicchangesinintersubject
correlationbetweensubcorticalregionsandthesaliencenetworkweredetectedforbothapproachandretreat.Forexample,fortherightPAG,intersubjectcorrelation
increasedduringapproachanddecreasedduringretreat.ConventionsasinFigure5.
value of connections (thus disregarding polarity), adding a (Figure10). Salience-network negative weights increased
constant positive value to all weights (effectively making all during retreat and decreased during approach. Negative
weights positive and negative weights smaller than positive weights between the salience and task-negative networks
ones),orinthepresentcontextsummingpositiveandnegative increased during approach and decreased during retreat. These
connections(whichcouldleadtopositiveandnegativeweights results were thus opposite to what was observed for positive
“cancelingout”). weights.
Inspectionoftheintersubjectcorrelationmatrixwithoutany
thresholding (that is, considering positive and negative weights
in equation 1) revealed a predominance of negative weights DISCUSSION
betweenthesalienceandtask-negativenetworks(Figure9).To
explorepotentialinformationofnegativeweights,weconsidered In the present study, we employed intersubject functional
equation(1)butonlyfornegativeweights(thus,positiveweights correlation analysis to investigate large-scale networks
were excluded). We focused on the weights within the salience during threat approach and retreat. We found that functional
network and between the salience and task-negative networks connectivitywithinandbetweennetworkschangeddynamically
FrontiersinHumanNeuroscience|www.frontiersin.org 10 November2017|Volume11|Article552