PPopulation-level Task-evoked Functional Connectivity
Kun Meng and Ani Eloyan Department of Biostatistics, Brown University School of Public Health, Providence, RI 02903, USA * Corresponding Author, e-mail: kun [email protected] , Address: 121 S. Main St, Providence, RI 02903, USA.
March 2, 2021
Abstract
Functional magnetic resonance imaging (fMRI) is a non-invasive, in-vivo imaging technique essential formeasuring brain activity. Functional connectivity is used to study associations between brain regionseither at rest or while study participants perform tasks. This paper proposes a rigorous definition oftask-evoked functional connectivity at the population level (ptFC). Importantly, our proposed ptFC isinterpretable in the context of task-fMRI studies. Two algorithms for estimating ptFC are provided.We present the performance of the proposed algorithms compared to existing functional connectivityestimation approaches using simulations. Lastly, we apply the proposed framework to estimate functionalconnectivity in a motor-task study from the Human Connectome Project.
Keywords:
AMUSE algorithm, Human Connectome Project, motor-task, weakly stationary with mean zero
Functional magnetic resonance imaging (fMRI) is a non-invasive brain imaging technique used to estimateboth brain regional activity and interactions between brain regions due to its ability to detect neuronalactivity throughout the entire brain simultaneously. In a typical fMRI study, functional neuronal activitytime courses are measured on three-dimensional volume elements (called voxels ) during a certain period oftime. Each time course is observed at discrete time points. For example, if the fMRI study lasts for 20minutes and data are collected at every 2 seconds, i.e., the time between successive brain scans referred to
Abbreviations:
BOLD, blood-oxygenation level-dependent; FC, functional connectivity; fMRI, functional magnetic res-onance imaging; GLM, general linear models; HCP, Human Connectome Project; HRF, hemodynamic response function; iid,independent and identically distributed; ptFC, population-level task-evoked functional connectivity; ptFCE, population-leveltask-evoked functional connectivity estimation; ROI, region of interest; sd, standard deviation; TR, repetition time. a r X i v : . [ s t a t . A P ] F e b s repetition time (TR) is 2 seconds, then we obtain observations at 900 time points. In this paper, we areinterested in investigating neuronal activity. However, since neuronal activity occurs in milliseconds, it isnot possible to directly observe it using fMRI technology. It is implicitly captured by the blood-oxygenationlevel-dependent (BOLD) signals. Neuronal activity is followed by localized changes in metabolism. Whenthe localized neuronal activity in a brain area occurs, the corresponding local oxygen consumption increases,and then oxygen-rich blood flows to this area. This process results in an increase in oxyhemoglobin anda decrease in deoxyhemoglobin. The BOLD signal value at each time point is the difference between theoxyhemoglobin and deoxyhemoglobin levels. The time course consisting of the BOLD signal values acrossall time points measures the localized metabolic activity influenced by the local vasculature of the brain,and it indirectly measures the localized neuronal activity. FMRI captures BOLD signals while the studyparticipants are either resting (resting-state fMRI) or performing a task (task-fMRI).We are interested in modeling brain neuronal activity and BOLD signals of participants ω , where ω ∈ Ωand Ω is the collection of participants in the study. At each time point t , we denote by Y k ( ω ; t ) the BOLDsignal value at the k th node of participant ω ’s brain. The word “node” throughout this paper refers to eithera voxel or a regions of interest (ROI). Aggregated BOLD signals at a macro-area level are often of interest.That is, for each participant, BOLD signals are spatially averaged within pre-selected ROI and the resultingROI-specific signals are analyzed. The BOLD signals for participant ω are represented by a K -vector-valuedfunction { YYY ( ω ; t ) = ( Y ( ω ; t ) , Y ( ω ; t ) , · · · , Y K ( ω ; t )) T | t ∈ T } on T , where T is the collection of time indicesand K is the number of nodes. Furthermore, we assume that T is a compact subset of R . In this paper, wemodel three types of signals in the context of task-fMRI studies as illustrated in Figure 1: (i) the stimulussignals denoted by N ( t ) representing the experimental design of the task of interest, (ii) the task-evokedneuronal activity signals at the k th node Φ k [ ω ; N ( t )] stemming from the task stimulus N ( t ), where the“stimulus-to-activity” map Φ k [ ω ; • ] : N ( t ) (cid:55)→ Φ k [ ω ; N ( t )] depends on participants ω , and (iii) the observedBOLD signals Y k ( ω ; t ) = Ψ k { ω ; Φ k [ ω ; N ( t )] } that are associated with the task-evoked neuronal activityΦ k [ ω ; N ( t )], where the “activity-to-BOLD” map Ψ k { ω ; •} : Φ k [ ω ; N ( t )] (cid:55)→ Y k ( ω ; t ) is participant-dependent.In (ii) and (iii), the map Φ k [ ω ; • ] characterizes how neurons at the k th node of ω react to stimulus N ( t ), andthe map Ψ k { ω ; •} describes how neuronal activity Φ k [ ω ; N ( t )] induces BOLD signal Y k ( ω ; t ). Therefore, weare interested in modeling signals Φ k [ ω ; N ( t )]. However, only BOLD signals Y k ( ω ; t ) are observable. Thegoal of many fMRI studies, including our proposed method, is to recover the information in Φ k [ ω ; • ] byanalyzing Y k ( ω ; t ). In Section 2, we provide a nonparametric model for Φ k [ ω ; • ] and Ψ k { ω ; •} .Theoretically, time t is a continuous variable and T = [0 , t ∗ ], where 0 < t ∗ < ∞ denotes the end of the2 k [ ω ; N ( t )] Ψ k { ω ; Φ k [ ω ; N ( t )] } N ( t )ExperimentsFigure 1: Stimulus signals N ( t ) are known and determined by experimental designs of interest. For ex-ample, N ( t ) are finite linear combinations of boxcar functions in block designs, and N ( t ) are finite linearcombinations of point masses in event-related designs. The task-evoked neuronal activity signals Φ k [ ω ; N ( t )]responding to N ( t ) are participant-dependent functions of N ( t ), i.e., Φ k [ ω ; • ] : N ( t ) (cid:55)→ Φ k [ ω ; N ( t )]. Addi-tionally, we assume Φ k [ ω ; N ( t )] ≡ N ( t ) ≡
0, i.e., there is no task-evoked neuronal activity at the k th node when there is no task stimulus. BOLD signals stem from Φ k [ ω ; N ( t )] and are characterized by theparticipant-dependent maps Φ k [ ω ; N ( t )] (cid:55)→ Ψ k { ω ; Φ k [ ω ; N ( t )] } = Y k ( ω ; t ).experiment. In applications, we obtain data only at discrete and finite time points in T = { τ × T R | τ =0 , , · · · , T } , where T R is the predetermined repetition time and T indicates that BOLD signals are observedat T + 1 time points. A BOLD signal, e.g., the one at the k th node of participant ω , in task-fMRI consists ofthree components: (i) P k ( ω ; t ) denotes one evoked by the experimental task N ( t ), (ii) Q k ( ω ; t ) denotes theone stemming from spontaneous brain activity and the neuronal activity responding to stimuli independentof N ( t ), and (iii) the random error (cid:15) k ( ω ; t ). Throughout this paper, we assume that these three componentshave an additive structure, and the observed BOLD signals Y k ( ω ; t ) are of the following form. Y k ( ω ; t ) = P k ( ω ; t ) + Q k ( ω ; t ) + (cid:15) k ( ω ; t ) , k = 1 , , · · · , K, (1.1)where P k ( ω ; t ) are called task-evoked terms .We implement the methods proposed in this paper to estimate functional connectivity in the context ofunderstanding the human motor system. To investigate motor functional connectivity, we use a cohort ofparticipants from a task-evoked functional MRI study publicly available at the Human Connectome Project(HCP, https://protocols.humanconnectome.org/HCP/3T/task-fMRI-protocol-details.html ). Theblock-design motor task used in this study is adapted from experiments by Buckner et al. (2011) andYeo et al. (2011), while details on the HCP implementation are given by Barch et al. (2013). During theexperiment, the participants are asked to perform five different tasks when presented with a cue: tap theirleft/right fingers, squeeze their left/right toes, and move their tongues. BOLD signals Y k ( ω j ; t ) are collectedfrom 308 participants, i.e. the set of participants is defined as Ω = { ω j } j =1 , while each BOLD signal is ob-tained with T R = 0 .
72 (seconds). Each experiment lasted for about 205 seconds, i.e., T = 283. We focus onthe estimation of brain functional connectivity for the task of squeezing right toes. The stimulus signal hastwo 12-second task blocks, such that the onsets of these task blocks vary across participants ω j . Neverthe-less, the corresponding onsets of any two participants differ in less than 0.1 seconds. Therefore, we assume3hat all participants share the same stimulus signal defined by N ( t ) = N rt ( t ) = [86 . , . ( t ) + [162 , ( t ).In statistical analysis of task-fMRI data { YYY ( ω ; t ) } t ∈T defined by (1.1), two general topics are primarilyof interest: (i) identification of nodes presenting task-evoked neuronal activity, i.e., the indices k suchthat P k ( ω ; t ) (cid:54) = 0, and (ii) identification of task-evoked interactions between brain nodes. General linearmodels (GLM, Friston et al. (1994)) are commonly implemented to detect task-evoked nodes (Lindquistet al. (2008)). GLM is conducted at an individual node level and is not informative for investigating theinteraction between nodes. Interactions between neuronal activity at pairs of nodes captured by BOLDsignals are referred to as functional connectivity (FC) (Friston et al. (1993), Friston (1994), Friston et al.(1994), and Friston (2011)). FC observed during a task experiment tends to be different from that observedin a resting-state experiment (see Hampson et al. (2002) and Lowe et al. (2000)). The difference betweentask-fMRI and resting-state fMRI is presented by the non-zero task-evoked terms P k ( ω ; t ) (cid:54) = 0. In this paper,we investigate the interactions between task-evoked terms { P k ( ω ; t ) } Kk =1 . These interactions describe the FCstructure stemming solely from the task stimulus N ( t ) and reveal how information of N ( t ) is transmittedacross the brain. Specifically, rather than FC between BOLD signals { Y k ( ω ; t ) } Kk =1 , we investigate the FCbetween task-evoked terms { P k ( ω ; t ) } Kk =1 .A considerable amount of work has been done to define and estimate FC between brain nodes. Forexample, Friston et al. (1993) defines FC as the temporal Pearson correlation between a pair of BOLDsignals across time t . Following the Pearson correlation approach, one may measure the task-evoked FCbetween P k ( ω ; t ) and P l ( ω ; t ) by the absolute value of the Pearson correlation as follows (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:82) T (cid:104) P k ( ω ; t ) − µ ( T ) (cid:82) T P k ( ω ; s ) µ ( ds ) (cid:105) × (cid:104) P l ( ω ; t ) − µ ( T ) (cid:82) T P l ( ω ; s ) µ ( ds ) (cid:105) µ ( dt ) (cid:114)(cid:82) T (cid:12)(cid:12)(cid:12) P k ( ω ; t ) − µ ( T ) (cid:82) T P k ( ω ; s ) µ ( ds ) (cid:12)(cid:12)(cid:12) µ ( dt ) × (cid:82) T (cid:12)(cid:12)(cid:12) P l ( ω ; t ) − µ ( T ) (cid:82) T P l ( ω ; s ) µ ( ds ) (cid:12)(cid:12)(cid:12) µ ( dt ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (1.2)where, if T = [0 , T ∗ ], µ ( dt ) is the Lebesgue measure dt ; if T = { τ × T R | τ = 0 , , · · · , T } , µ ( dt ) is thecounting measure (cid:80) τ ∈ Z δ τ × T R ( dt ), where δ τ × T R is the point mass at time point τ × T R . There are manyother approaches to defining and estimating FC developed in the last three decades. For example, coherenceanalysis (M¨uller et al. (2001), Sun et al. (2004)) is an analog of Pearson correlation in the Fourier frequencyspace and is a widely used method. Additionally, beta-series regression (Rissman et al. (2004)) is widelyused as well for FC studies. These approaches have several limitations discussed in Section 2.2.The paper is organized as follows. In Section 2, we propose a model for task-evoked BOLD signals; basedon this model, we propose a rigorous and interpretable definition of task-evoked functional connectivity at4 population level (ptFC). In Section 3, we present two algorithms for estimating ptFC. We compare theperformance of our proposed algorithms with existing approaches using simulations in Section 4. In Section5, we apply the proposed approach to estimate the ptFC during a motor task using the publicly availableHCP data set. Section 6 concludes the paper and provides further discussions.
We first propose a model for BOLD signals at individual nodes. Using this model, we provide the definitionof ptFC. Lastly, we list advantages of our proposed ptFC compared to existing approaches. To describe thestatistical properties of the participant population Ω, we assume throughout this paper that an underlyingprobability space (Ω , P ) is defined on Ω, where P is a probability measure (see Chapter 1.1 of Durrett (2019)for the definition of probability spaces and related materials). Let h k ( t ) be the hemodynamic response function (HRF, Lindquist et al. (2008)) at the k th node correspondingto the task stimulus N ( t ) and shared by all participants ω ∈ Ω. For each participant ω , the GLM approachessentially models the task-evoked terms in (1.1) corresponding to N ( t ) by P k ( ω ; t ) = β k ( ω ) × ( N ∗ h k )( t − t ,k ) , t ∈ T , ω ∈ Ω , k = 1 , , · · · , K. (2.1)where β k ( ω ) is a GLM regression coefficient depending on ω , t ,k is a reaction delay time as the reaction ofthe k th node to N ( t ) at time t is not instantaneous, and the asterisk ∗ in ( N ∗ h k ) denotes the convolutionoperation. Specifically, this is a GLM with response { P k ( ω ; t ) } t ∈T and independent variable { ( N ∗ h k )( t − t ,k ) } t ∈T , and the error term of this GLM is absorbed by the error { (cid:15) k ( ω ; t ) } t ∈T in (1.1). The task-evokedterms P k ( ω ; t ) in (2.1) contain the information of the neuronal activity evoked by the stimulus N ( t ), andthe quantity β k ( ω ) measures the magnitude of the task-evoked neuronal activity. Visualizations of (2.1) arepresented in Figure 2. Motivated by (1.1), we model the BOLD signals Y k ( ω ; t ) as follows. Y k ( ω ; t ) = β k ( ω ) × ( N ∗ h k ) ( t − t ,k ) + R k ( ω ; t ) , t ∈ T , ω ∈ Ω , k = 1 , , · · · , K, (2.2)where the terms R k ( ω ; t ) = Q k ( ω ; t ) + (cid:15) k ( ω ; t ), with Q k ( ω ; t ) and (cid:15) k ( ω ; t ) defined in (1.1), are referred to as reference terms and contain the information of (i) the neuronal activity from stimuli independent of N ( t ),5ii) spontaneous neuronal activity, and (iii) random error. The model (2.2) is a random effects model andformally equivalent to the equation (1) in Warnick et al. (2018) and the equation (10) in Joel et al. (2011).As described above, the task-evoked terms P k ( ω ; t ) = β k ( ω ) × ( N ∗ h k )( t − t ,k ) = { [ β k ( ω ) × N ( · − t ,k )] ∗ h k } ( t ) contain the information of the neuronal activity evoked by N ( t ). In P k ( ω ; t ), HRF h k representmetabolism and vasculature and do not characterize neuronal activity. Therefore, we model neuronal activitysignals Φ k [ ω ; N ( t )] responding to N ( t ) (see Figure 1) as β k ( ω ) × N ( t − t ,k ). Using model (2.2), the “stimulus-to-activity” map Φ k [ ω ; • ] and the “activity-to-BOLD” map Ψ k { ω ; •} in Figure 1 are expressed as follows.Φ k [ ω ; • ] : N ( t ) (cid:55)→ β k ( ω ) × N ( t − t ,k ) , Ψ k { ω ; •} : Φ k [ ω ; N ( t )] = β k ( ω ) × N ( t − t ,k ) (cid:55)→ { Φ k [ ω ; N ( · )] ∗ h k } ( t ) + R k ( ω ; t ) = Y k ( ω ; t ) . (2.3)The reference terms R k ( t ) in (2.2) are general. Here, we provide an explicit example of R k ( t ) suitable formost applications. Let ˜ N γ ( t ), for γ = 1 , , · · · , Γ, define the stimulus signals other than the experimental taskof interest N ( t ). Let ˜ h k,γ ( t ) be the HRF at the k th node corresponding to stimulus ˜ N γ ( t ), for γ = 1 , , · · · , Γ.Then the explicit example of R k ( t ), for k = 1 , , · · · , K , is presented as follows. R k ( ω ; t ) = Q k ( ω ; t ) + (cid:15) k ( ω ; t ) , where Q k ( ω ; t ) = Γ (cid:88) γ =1 ˜ β k,γ ( ω ) × (cid:16) ˜ N γ ∗ ˜ h k,γ (cid:17) ( t − t ,k,γ ) + S k ( ω ; t ) , (2.4)where t ,k,γ is the reaction delay of the k th node corresponding to stimulus ˜ N γ ( t ), ˜ β γ ( ω ) depends on par-ticipants ω , S k ( ω ; t ) is the BOLD signal component arising from spontaneous neuronal activity at the k th node of ω , and (cid:15) k ( ω ; t ) is the random noise at the k th node of participant ω . In the HCP motor task study,Γ = 4 stimuli other than N rt ( t ) of interest are tapping left and right fingers, squeezing left toes, and movingthe tongue. These four stimuli are independent of N rt ( t ) and are not of interest when estimating the FCevoked by the task of interest; applying the notations in (2.4) to this example, stimulus signals of these fourstimuli are represented as ˜ N ( t ) = [71 . , . ( t ) + [177 . , . ( t ), ˜ N ( t ) = [11 , ( t ) + [116 . , . ( t ) , ˜ N ( t ) = [26 . , . ( t ) + [146 . , . ( t ) , and ˜ N ( t ) = [56 . , . ( t ) + [131 . , . ( t ) . FC is expected to characterize the mechanism of neuronal activity rather than metabolism or vasculature.Therefore, the task-evoked FC is defined using the task-evoked neuronal activity signals Φ k [ ω ; N ( t )]. Inmodel (2.3) of Φ k [ ω ; N ( t )], stimulus signals N ( t ) are fully determined by experimental design. Additionally,the reaction delay t ,k affects HRF h k since the task-evoked terms can be expressed as P k ( ω ; t ) = { [ β k ( ω ) × N ] ∗ [ h k ( · − t ,k )] } ( t ). Therefore, the task-evoked FC is expected to be determined by coefficients β k ( ω ).6ince β k ( ω ), for k = 1 , , · · · , K , measure the magnitude of the neuronal activity evoked by task stimulus N ( t ), the k th and l th nodes are functionally connected at the population level during the task if one of thefollowing scenarios holds: (i) for participants with strong reaction to N ( t ) at the k th node, their l th nodes’reaction to N ( t ) is strong as well and vice versa, i.e., they are positively correlated; (ii) for participants withstrong reaction to N ( t ) at the k th nodes, their l th nodes’ reaction to N ( t ) is weak and vice versa, i.e., they arenegatively correlated. Furthermore, we make the following two modeling assumptions on the distribution of β k ( ω ) across all participants ω ∈ Ω: (1) If there exists a connection evoked by task N ( t ) between k th and l th nodes, the corresponding random coefficients β k ( ω ) and β l ( ω ) are approximately linearly associated. (2) Eachvariance V β k = E (cid:0) β k (cid:1) − ( E β k ) , for = 1 , , · · · , K , is strictly positive, i.e., random variable β k : ω → β k ( ω )is not deterministic † ; this assumption is realistic as it allows the variability of task-evoked neuronal activity β k ( ω ) × N ( t − t ,k ) across population Ω. Since the correlation corr ( β k , β l ) = E ( β k β l ) − ( E β k )( E β l ) √ V β k × V β l measures thelinear statistical correlation between β k and β l , we define ptFC as follows. Definition 2.1.
For each ω in a population of interest Ω , the BOLD signal Y k ( ω ; t ) at the k th node, for k = 1 , , · · · , K , is of the form (2.2). Then the ptFC between the k th and l th nodes is defined as | corr ( β k , β l ) | . An advantage of Definition 2.1 is its scale-invariance. Since task-evoked term β k ( ω ) × ( N ∗ h k ) can berepresented as β k ( ω ) c × c × [( c N ) ∗ ( c h k )], the scale of β k ( ω ) changes if we change the scale of HRF h k or N ( t )in our model. However, corr ( β k , β l ) is invariant to scaling transforms β k (cid:48) (cid:55)→ cβ k (cid:48) , for k (cid:48) ∈ { k, l } , so is ptFC | corr ( β k , β l ) | . This scale-invariance solves the potential issue of identifiability and interpretability of β k . In this subsection, we present advantages of our proposed ptFC in Definition 2.1 compared to existingapproaches. We first discuss the limitations of the Pearson correlation approach (1.2) using model (2.1).Plugging (2.1) into (1.2), we obtain the following measurement of association between P k ( ω ; t ) and P l ( ω ; t ). | corr ( P k , P l ) | := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:82) T (cid:104) φ k ( t ) − µ ( T ) (cid:82) T φ k ( s ) µ ( ds ) (cid:105) × (cid:104) φ l ( t ) − µ ( T ) (cid:82) T φ l ( s ) µ ( ds ) (cid:105) µ ( dt ) (cid:114)(cid:82) T (cid:12)(cid:12)(cid:12) φ k ( t ) − µ ( T ) (cid:82) T φ k ( s ) µ ( ds ) (cid:12)(cid:12)(cid:12) µ ( dt ) × (cid:82) T (cid:12)(cid:12)(cid:12) φ l ( t ) − µ ( T ) (cid:82) T φ l ( s ) µ ( ds ) (cid:12)(cid:12)(cid:12) µ ( dt ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (2.5)where φ k (cid:48) ( t ) = N ∗ h k (cid:48) ( t − t ,k (cid:48) ) for k (cid:48) ∈ { k, l } , and (2.5) reveals the following limitations of Pearsoncorrelation approach (1.2). † Expectation E is associated with probability measure P . Specifically, E ( β νk ) = (cid:82) Ω β k ( ω ) ν P ( dω ), for all ν ∈ N .
20 40 60 80 − . . . . . Time T a sk T e r m s k − th nodel − th node 0 5 10 15 20 25 30 − . . . . . . . Time
DefaultNon−default 0 20 40 60 80 − . . . . . Time T a sk T e r m s k − th nodel − th node 0 20 40 60 80 − . . . . . Time T a sk T e r m s k − th nodel − th node (a): | corr ( P k , P l ) | = 0 .
571 (b) (c): | corr ( P k , P l ) | = 0 .
445 (d): | corr ( P k , P l ) | = 0 . N ( t ) = (cid:80) m =1 (20 m − , m ] ( t ) be the stimulussignal of a block design, T = [0 , h k be the canonicalHRF function with default parameters in the R package neuRosim . h k is illustrated by the blue curve in (b). Panel (a) shows the influence of the variationin reaction delay on (2.5), where h l ( t ) = h k ( t + 3), i.e., t ,k = −
3; the task-evoked terms at the k th and l th nodes are presented by blue and red curves, respectively. In (b), h l is replaced by the canonicalHRF function with parameters a1 = 10 , a2 = 15 , b1 = b2 = 0 . , c = 0 .
35, and presented by the red (dashed)curve. Panel (c) shows the influence of variation in HRF on (2.5), and the task-evoked terms at the k th and l th nodes are presented by blue and red curves, respectively, where h l is defined in (b) as the red curve.Panel (d) is a noise-contaminated version of (c), i.e., curves of P k ( t ) + (cid:15) k ( t ) (blue) and P l ( t ) + (cid:15) l ( t ) (red)with (cid:15) k ( t ) , (cid:15) l ( t ) ∼ iid N (0 , .
1) for each t . Panel (d) shows that random noise can further influence (2.5). Inpanels (a,c,d), the presented | corr ( P k , P l ) | is defined in (1.2). Only nuisance parameters:
Based on the reasoning in Section 2.1, the information on the mechanismof the task-evoked neuronal activity is contained in coefficients β k ( ω ). Hence, it is counterintuitive that(2.5) does not depend on β k ( ω ) and depends only on the nuisance parameters N ( t ) and h k when consideringquantification of brain neuronal activity. For example, if the k th node does not react to N ( t ), then thetask-evoked term β k ( ω ) × ( N ∗ h k )( t − t ,k ) is expected to be approximately zero, i.e., β k ( ω ) ≈
0, and thereshould be no task-evoked interaction between k th and all other nodes. However, (2.5) can still be verylarge as the non-reaction information contained in β k ( ω ) ≈ { ( β k ( ω ) , β l ( ω )) | ω ∈ Ω } , rather than nuisance parameters. Variation in reaction delay:
Reaction delay times t ,k may vary across nodes (for example, see Miezinet al. (2000) for the investigation of the left and right visual and motor cortices and the correspondingdifference between response onsets). If t ,k (cid:54) = t ,l , for example t ,k < t ,l , then it is likely that the k th nodereacts to the task-stimulus N ( t ) first, and then the neuronal activity at the k th node causes that at the l th node. Because of this potential causality represented by the difference between t ,k and t ,l , it is natural toexpect that the k th and l th nodes are likely functionally connected. However, the measurement (2.5) can bevery small if t ,k (cid:54) = t ,l and may not reveal the true interaction between neuronal activity of two nodes. An8xample of this issue is illustrated in Figure 2 (a). Since our proposed definition of ptFC does not involvereaction delay times, the variation in reaction delay does not influence the ptFC in Definition 2.1. Variation in HRF:
HRF can heavily vary across brain nodes (Buckner et al. (1998), Lee et al. (2001),Miezin et al. (2000), Saad et al. (2001), Schacter et al. (1997), Rosen et al. (1998)). Since task-evoked FC isnot expected to depend on HRF, variation of HRF in different brain regions should not influence task-evokedFC. However, the quantity defined in (2.5) can be small if h k (cid:54) = h l as illustrated by the example in Figure2 (c). Given that our proposed definition of ptFC does not involve the information of HRF, it is invariantto the variation in HRF across brain nodes.The beta-series regression approach is widely accepted. However, its estimation procedure is based onGLM and involves modeling HRF across all nodes. Since the HRF may vary across nodes, node-wise HRFmodeling is cumbersome, and the inaccuracy in modeling reduces the efficacy of beta-series regression. Ourproposed approach for defining ptFC and the forthcoming algorithms for estimating ptFC do not involvemodeling HRF and hence are not influenced by node dependent variability in HRF.Coherence analysis based estimation of task-evoked FC is not influenced by any of the issues discussedabove. In this approach, FC between BOLD signals Y k ( ω ; t ) and Y l ( ω ; t ) is measured by the coherence eval-uating the extent of the linear time-invariant relationship between these two signals via Fourier frequencies.The linear time-invariant relationship is equivalent to that there exists a (generalized) function κ such that Y k ( t ) = κ ∗ Y l ( t ) or Y l ( t ) = κ ∗ Y k ( t ) (Theorem 1.2 of H¨ormander (1960)), implying that the correspondingcoherence is equal to 1 at all frequencies. However, there is no guarantee that the relationship between twoBOLD signals is linear time-invariant if the corresponding two nodes are functionally connected. Addition-ally, the mathematical formulation implemented in M¨uller et al. (2001) and Sun et al. (2004) has limitedinterpretability. Equation (2) in M¨uller et al. (2001) and equation (3) in Sun et al. (2004) apply the auto-covariance defined as the expected value E [ Y k ( τ + t ) Y l ( t )]. However, the underlying probability space anddistribution for conducting the expectation E operation is not clearly interpreted therein. In the proposeddefinition of ptFC, we clearly specify the underlying probability space for expectation E . In this section, we propose two algorithms to estimate ptFC. Our proposed estimation algorithms do notdepend on any parametric assumptions and are based on Fourier transforms. Since BOLD signals in appli-cations are observed only at discrete and finite time points, in the sequel, we set the collection of time indices T to be { τ × T R | τ = 0 , , · · · , T } and t = τ × T R for some index τ ∈ { , , · · · , T } . Before introducing the9stimation algorithms, we provide notations for (discrete) Fourier transforms. Further materials on Fouriertransforms can be found in Bloomfield (2004). For any function f : T → R , we extend f as follows to be aperiodic function on { τ (cid:48) × T R | τ (cid:48) ∈ Z } f ( τ (cid:48) × T R ) = f ( τ × T R ) , τ (cid:48) ≡ τ (mod T + 1) for τ = 0 , , · · · , T and all τ (cid:48) ∈ Z . (3.1)In this paper, all functions defined on T are automatically extended using (3.1) to be periodic functions on { τ (cid:48) × T R | τ (cid:48) ∈ Z } . The (discrete) convolution ( N ∗ h k )( t ) is defined by ( N ∗ h k ) ( τ × T R ) = T +1 (cid:80) Tτ (cid:48) =0 N ( τ (cid:48) × T R ) h k (( τ − τ (cid:48) ) × T R ) for all t = τ × T R with τ = 0 , , · · · , T . We define the Fourier transform of f as (cid:98) f ( ξ ) := 1 T + 1 T (cid:88) τ =0 f ( τ × T R ) e − πiξ ( τ × T R ) for all ξ ∈ R . (3.2)It is straightforward that (cid:98) f ( ξ ) is a periodic function with period 1 /T R .Additionally, a technical assumption for our proposed algorithms is E β k = E R k ( t ) = 0, for k =1 , , · · · , K and t ∈ T motivated by the following centralization. † Y k ( ω ; t ) − E Y k ( t ) = [ β k ( ω ) − E β k ] ( N ∗ h k )( t − t ,k ) + [ R k ( ω ; t ) − E R k ( t )] . (3.3)Using [ β k ( ω ) − E β k ] as β k ( ω ) and [ R k ( ω ; t ) − E R k ( t )] as R k ( ω ; t ), we model the demeaned outcome. If thereference terms R k ( ω ; t ) satisfy (2.4), centralization (3.3) allows us to assume E ˜ β k,γ = E S k ( t ) = E (cid:15) k ( t ) = 0for all k = 1 , , · · · , K, γ = 1 , , · · · , Γ, and t ∈ T . In this subsection, we propose the ptFCE algorithm for the estimation of ptFC in studies where both task-evoked signals Y k ( ω ; t ) and reference signals R k ( ω ; t ) are observed. In the next subsection, we propose analgorithm based on ptFCE applicable in experiments where only task-evoked signals Y k ( ω ; t ) are available.In Section 2.1, our proposed definition of ptFC, given by | corr ( β k , β l ) | , does not depend on reactiondelay times t ,k or reference terms R k ( ω ; t ). To construct an estimator of | corr ( β k , β l ) | using the observedBOLD signals, we note that applying the periodic extension (3.1) to function ( N ∗ h k ), it is easy to verifythat the two stochastic processes ( N ∗ h k )( t − t ,k − U ( ω )) and ( N ∗ h k )( t − U ( ω )), where the randomvariable U : Ω → T , ω (cid:55)→ U ( ω ) is uniformly distributed on T , are identically distributed. In the meantime, † Throughout this paper, for any real-valued stochastic process { Z ( ω ; t ) | ω ∈ Ω } t ∈T , the expected value E Z ( t ) denotes theaverage of Z ( ω ; t ) across all ω ∈ Ω, i.e., (cid:82) Ω Z ( ω ; t ) P ( dω ), for each fixed t ∈ T . N ∗ h k )( t − U ( ω )), and hence the distribution of ( N ∗ h k )( t − t ,k − U ( ω )), does not depend on the nuisanceparameter t ,k . Therefore, we investigate the following time-shifted signals Y k ( ω ; t − U ( ω )) = β k ( ω )( N ∗ h k ) ( t − t ,k − U ( ω )) + R k ( ω ; t − U ( ω )) , for k = 1 , , · · · , K, (3.4)Based on the discussion above, the distributions of Y k ( ω ; t − U ( ω )) do not depend on the nuisance parameter t ,k . Then, we consider the autocovariance E [ Y k ( t − U ) Y l ( t + s − U )] when estimating the ptFC betweennodes k and l , where s = s × T R with indices s ∈ Z . Additionally, ptFC | corr ( β k , β l ) | only depends ontask-evoked terms [ Y k ( ω ; t − U ( ω )) − R k ( ω ; t − U ( ω ))] (see (3.4)). Hence, we investigate the autocovariancedifferences E [ Y k ( t − U ) Y l ( t + s − U )] − E [ R k ( t − U ) R l ( t + s − U )] obtained by filtering out the referenceterms from the observed BOLD signals. If U ( ω ), { β k ( ω ) } Kk =1 , and { R k ( ω ; t ) | t ∈ T } Kk =1 are independent,Theorem 8.1 in the Appendix implies that this autocovariance difference depends only on s and does notdepend on t . Then we denote this difference of interest as follows. A kl ( s ) = E [ Y k ( t − U ) Y l ( t + s − U )] − E [ R k ( t − U ) R l ( t + s − U )] , s = s × T R for all s ∈ Z , (3.5)for k, l = 1 , , · · · , K . Finally, we propose the following estimator for ptFC | corr ( β k , β l ) | by incorporating anormalization multiplier in (3.5). C kl ( ξ ) := (cid:12)(cid:12)(cid:12) (cid:98) A kl ( ξ ) (cid:12)(cid:12)(cid:12)(cid:114)(cid:12)(cid:12)(cid:12) (cid:98) A kk ( ξ ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) (cid:98) A ll ( ξ ) (cid:12)(cid:12)(cid:12) , for all ξ ∈ R , (3.6)where (cid:98) A kl ( ξ ) denotes the Fourier transform of A kl ( s ) defined in (3.2). Periodicity property of the Fouriertransform implies that C kl ( ξ ) is a periodic function of ξ with period 1 /T R . Lemma 8.1 in the Appendiximplies that the ratio in (3.6) is well-defined for all ξ ∈ [0 , /T R ] except for at most finitely many points ξ suchthat (cid:98) A kk ( ξ ) (cid:98) A ll ( ξ ) = 0. Importantly, based on Theorem 8.1, we obtain C kl ( ξ ) = | corr ( β k , β l ) | , for all ξ ∈ R ,i.e., C kl ( ξ ) is an unbiased estimator of ptFC for all ξ under the model (2.2). An algorithm for computing C kl ( ξ )in (3.6) will be provided in this subsection. In the implementation of the algorithm, the random variable U ( ω ) is artificially generated and automatically independent of { β k ( ω ) , R k ( ω ; t ) | t ∈ T , k = 1 , , · · · , K } .The model for BOLD signals (2.2) is based on the assumption that the reference terms R k ( ω ; t ) = Q k ( ω ; t ) + (cid:15) k ( ω ; t ) contain the total random noise (cid:15) k ( ω ; t ). However, in applications it is more realistic toassume that the reference terms contain only a part of the random noise (cid:15) k ( ω ; t ). Specifically, we assume that11he random noise (cid:15) k ( ω ; t ) can be additively decomposed into two parts (cid:15) k ( ω ; t ) = V k ( ω ; t ) + W k ( ω ; t ), wherethe reference terms can be used to model only a part of the random noise, i.e., R k ( ω ; t ) = Q k ( ω ; t ) + V k ( ω ; t ),and W k ( ω ; t ) = W k ( ω ; t ) is the “residual noise” not modeled by R k ( ω ; t ). Therefore, it is more realistic tomodel task-evoked BOLD signals Y k ( t ) by the following “noise-contaminated” version of (2.2). Y k ( ω ; t ) := β k ( ω ) × N ∗ h k ( t − t ,k ) + W k ( ω ; t ) + R k ( ω ; t ) , t ∈ T , ω ∈ Ω , k = 1 , , · · · , K. (3.7)Before modeling residual noise W k ( ω ; t ), we introduce a concept: a vector-valued stochastic process { G ( ω ; t ) =( G ( ω ; t ) , G ( ω ; t ) , · · · , G K ( ω ; t )) T } t ∈T is called weakly stationary with mean zero if E G k ( t ) = 0, for all t ∈ T , and E [ G k ( t ) G l ( t + s )] depends only on s , rather than t , for all k, l . We assume { W ( ω ; t ) :=( W ( ω ; t ) , W ( ω ; t ) , · · · , W K ( ω ; t )) T } t ∈T satisfies (i) W ( ω ; t ) is independent of W ( ω ; t ) if t (cid:54) = t , (ii)Σ kl = E [ W k ( t ) W l ( t )] for all t ∈ T , and (iii) W ( ω ; t ) is weakly stationary with mean zero. If { W ( ω ; t ) } t ∈T , U ( ω ), { β k ( ω ) } Kk =1 , and { R k ( ω ; t ) | t ∈ T } Kk =1 are independent, Theorem 8.2 in the Appendix implies C kl ( ξ ) = (cid:12)(cid:12)(cid:12)(cid:12) E ( β k β l ) (cid:98) h k ( ξ ) (cid:98) h l ( ξ ) | (cid:98) h k ( ξ ) (cid:98) h l ( ξ ) | e πiξ ( t ,k − t ,l ) + Σ kl ( T +1) | (cid:98) N ( ξ ) | | (cid:98) h k ( ξ ) (cid:98) h l ( ξ ) | (cid:12)(cid:12)(cid:12)(cid:12)(cid:115)(cid:18) E ( β k ) + Σ kk ( T +1) | (cid:98) N ( ξ ) (cid:98) h k ( ξ ) | (cid:19) (cid:18) E ( β l ) + Σ ll ( T +1) | (cid:98) N ( ξ ) (cid:98) h l ( ξ ) | (cid:19) , for all ξ ∈ R , (3.8)where C kl ( ξ ) is defined in (3.6). For a fixed T , (3.8) implies C kl ( ξ ) ≈ | corr ( β k , β l ) | at frequencies ξ suchthat | (cid:98) h k ( ξ ) | and | (cid:98) h l ( ξ ) | are sufficiently large. For most models of HRF h k , the absolute value of Fouriertransform | (cid:98) h k ( ξ ) | takes large values when ξ ∈ (0 , . | corr ( β k , β l ) | by the median of C kl ( ξ ) across ξ ∈ (0 , . C kl ( ξ ) in (3.6) andthe discussion above on the choice of ξ , we propose the population-level task-evoked functional connectivityestimation (ptFCE) algorithm (Algorithm 1) for the estimation of ptFC.To calculate the C kl ( ξ ), we need to estimate two components using data from the experiment of interest.The component E [ Y k ( t − U ) Y l ( t + s − U )] is estimated using the task-evoked BOLD signals Y k ( ω ; t ), and E [ R k ( t − U ) R l ( t + s − U )] is estimated using the reference terms R k ( ω ; t ). Therefore, the experiment maybe designed as follows: during the first half of the experiment, a resting-state scan is obtained and referencesignals R k ( ω ; t ) are collected; in the second half, the scans stimulus N ( t ) is implemented, and task-evokedBOLD signals Y k ( ω ; t ) are collected. However, not all experiments are designed to include a resting-statesegment. If the reference terms R k ( ω ; t ) are not available, a compromised approach proposed in the nextsubsection may be used for estimation. 12 lgorithm 1 ptFCE Algorithm: Input: (i) Task-evoked signals { YYY ( ω j ; τ × T R ) = ( Y k ( ω j ; τ × T R ) , Y l ( ω j ; τ × T R )) T } Tτ =0 of participants { ω j } Mj =1 in population Ω; (ii) reference signals { RRR ( ω j ; τ × T R ) = ( R k ( ω j ; τ × T R ) , R k ( ω j ; τ × T R )) T } Tτ =0 of participants { ω j } Mj =1 ; (iii) repetition time T R . Output:
An estimation of the ptFC | corr ( β k , β l ) | between the k th and l th nodes. Zero-mean:
For each τ = 0 , , · · · T and j = 1 , , · · · , M , do YYY ( ω j ; τ × T R ) ← YYY ( ω j ; τ × T R ) − M M (cid:88) j (cid:48) =1 YYY ( ω j (cid:48) ; τ × T R ) ,RRR ( ω j ; τ × T R ) ← RRR ( ω j ; τ × T R ) − M M (cid:88) j (cid:48) =1 RRR ( ω j (cid:48) ; τ × T R ) . For each j = 1 , , · · · , M , periodically extend YYY ( ω j , τ × T R ) and
RRR ( ω j , τ × T R ) using (3.1). Generate M integers { u j } Mj =1 from the uniform distribution on { , , , · · · , T } . For each s = 0 , , , · · · , T and ( k (cid:48) , l (cid:48) ) ∈ { ( k, l ) , ( k, k ) , ( l, l ) } , compute A k (cid:48) l (cid:48) ( s × T R ) ← T + 1 T (cid:88) τ =0 {E y ( τ ; s ) − E r ( τ ; s ) } , where E y ( τ ; s ) ← M M (cid:88) j =1 (cid:110) Y k (cid:48) ( ω j ; ( τ − u j ) × T R ) × Y l (cid:48) ( ω j ; ( t + s − u j ) × T R ) (cid:111) , E r ( τ ; s ) ← M M (cid:88) j =1 (cid:110) R k (cid:48) ( ω j ; ( τ − u j ) × T R ) × R l (cid:48) ( ω j ; ( t + s − u j ) × T R ) (cid:111) . For all ξ ∈ (0 , / (2 × T R )), compute the Fourier transforms (cid:100) A kl ( ξ ), (cid:100) A kk ( ξ ), and (cid:99) A ll ( ξ ) using (3.2). Compute C kl ( ξ ) ← | (cid:100) A kl ( ξ ) | / [ | (cid:100) A kk ( ξ ) (cid:99) A ll ( ξ ) | ] / . Estimate the ptFC | corr ( β k , β l ) | by the median of C kl ( ξ ) across all ξ ∈ (0 , . In many applications, the reference signals R k ( ω ; t ) are not available. To estimate ptFC when referencesignals are not available, we apply the AMUSE algorithm (Tong et al. (1991)) to estimate reference signalsand then ptFC are estimated using the ptFCE algorithm provided in Section 3.1. As a result, we proposethe AMUSE-ptFCE algorithm (Algorithm 2).First, we note that E [( N ∗ h k )( t − t ,k − U )] is a deterministic constant not depending on t , and we denotethis constant by C k . Define the stochastic processes J k ( ω ; t ) = β k ( ω ) [( N ∗ h k )( t − t ,k − U ( ω )) − C k ], for k = 1 , , · · · , K . For each fixed k , using (8.1) in the Appendix, one can verify that the scalar-valued stochasticprocess J k ( ω ; t ) is weakly stationary with mean zero conditioning on β k , then E (cid:8) J k ( t ) × J k ( t + s ) (cid:12)(cid:12) β k (cid:9) depends only on s . The following theorem gives the foundation for the AMUSE-ptFC algorithm. Theorem 3.1.
For each fixed k ∈ { , , · · · , K } , suppose task-evoked BOLD signals Y k ( ω ; t ) , for ω ∈ − . . . . . . . Time HR F s Canonical HRFNon−canonical HRF 0.0 0.2 0.4 0.6 0.8 1.0 . . . . . . Frequencies F ou r i e r t r an s f o r m s o f HR F s Canonical HRFNon−canonical HRFCanonical HRF with xi<0.1Non−canonical HRF with xi<0.1 0.00 0.05 0.10 0.15 0.20 0.25 . . . . . . Frequencies F ou r i e r t r an s f o r m s o f HR F s Figure 3: The canonical HRF (blue solid curve in the left panel) is plotted using the canonicalHRF functionin the R package neuRosim with default parameters. An example non-canonical HRF (red dahsed curve inthe left panel) is constructed by changing the function parameters to a1 =10, a2 =15, b1 = b2 =0.9, c =0.35.The middle panel shows the absolute values of the Fourier transforms of these two HRF when T R = 0 . T R = 2. In boththe middle and right panels, we plot the curves only in the corresponding semiperiod intervals (0 , / (2 T R )).Ω , are of the form in (2.2), the scalar-valued stochastic process { R k ( ω ; t ) } t ∈T is weakly stationary withmean zero, and the random variable U : Ω → T is uniformly distributed. Additionally, for each fixed k , β k ( ω ) , { R k ( ω ; t ) } t ∈T , and U ( ω ) are independent, and there exists t ∗ k ∈ T such that E (cid:110) J k ( t ) J k ( t − t ∗ k ) (cid:12)(cid:12) β k (cid:111) E (cid:110) J k ( t ) (cid:12)(cid:12) β k (cid:111) (cid:54) = E { R k ( t ) R k ( t − t ∗ k ) } E { R k ( t ) } . If the matrix AAA = ( a ij ) ≤ i,j ≤ and the stochastic process { sss ( ω ; t ) = ( s ( ω ; t ) , s ( ω ; t )) T } t ∈T satisfy (i) given β k , sss ( ω ; t ) is weakly stationary with mean zero, (ii) E [ s ( t ) s ( t + s ) | β k ] = 0 for all t and s ,(iii) E { s ( t ) s ( t − t ∗ k ) | β k } E { s ( t ) | β k } (cid:54) = E { s ( t ) s ( t − t ∗ k ) | β k } E { s ( t ) | β k } , and AAA s ( ω ; t ) s ( ω ; t ) = β k ( ω ) J k ( ω ; t ) R k ( ω ; t − U ( ω )) = Y k ( ω ; t − U ( ω )) − β k ( ω ) C k ( N ∗ h k )( t − t ,k − U ( ω )) − C k , (3.9) for all t ∈ T , then there exist a diagonal matrix ΛΛΛ and a permutation matrix
PPP , such that β k ( ω ) = AAA
ΛΛΛ − PPP − and J k ( ω ; t ) R k ( ω ; t − U ( ω )) = PPP
ΛΛΛ sss ( ω ; t ) , for all t ∈ T . (3.10)The proof of Theorem 3.1 is provided in the Appendix. For each k and ω , the pair ( AAA, sss ( ω ; t )) in (3.9) canbe derived by the AMUSE algorithm with the following 2-vector-valued signal as its input. (cid:110) ( Y k ( ω ; t − U ( ω )) − β k ( ω ) C k , ( N ∗ h k )( t − t ,k − U ( ω )) − C k ) T (cid:12)(cid:12) t ∈ T (cid:111) . (3.11)14n our proposed AMUSE-ptFCE algorithm, we assume that h k ( t ) and t ,k are known. For example, we mayset t ,k = 0 and h k ( t ) to be the R function canonicalHRF with default parameters. Since h k ( t ), t ,k , N ( t ),and the distribution of U ( ω ) are known, the constants C k are known. However, the coefficients β k ( ω ) areunknown, and then the input (3.11) is not available in applications. Since E β k = 0, in applications, weignore the term β k ( ω ) C k in (3.11) and apply (cid:110) ( Y k ( ω ; t − U ( ω )) , ( N ∗ h k )( t − t ,k − U ( ω )) − C k ) T (cid:12)(cid:12) t ∈ T (cid:111) as the input of the AMUSE algorithm for deriving ( AAA, sss ( ω ; t )) instead of (3.11). To approximately recoverreference terms R k ( ω ; t ) from the pair ( AAA, sss ( ω ; t )) in (3.9), we show the following result. Theorem 3.2.
For each k ∈ { , , · · · , K } and ω ∈ Ω , suppose the pair ( AAA, sss ( ω ; t )) satisfies (3.10). Thenthere exists i (cid:48) ∈ { , } such that a i (cid:48) s i (cid:48) ( ω ; t ) = J k ( ω ; t ) . The proof of Theorem 3.2 is in the Appendix. In applications, the index i (cid:48) in Theorem 3.2 is computedas i (cid:48) = argmax i (cid:8) corr ( s i ( ω ; t ) , ( N ∗ h k )( t − t ,k − U ( ω j )) − C k ) , across t (cid:12)(cid:12) i = 1 , (cid:9) . The AMUSE algorithmand Theorem 3.2 approximately recover J k ( ω ; t ) = β k ( ω )( N ∗ h k )( t − t ,k − U ( ω )) − β k ( ω ) C k . Again, weignore the term β k ( ω ) C k as E β k = 0, i.e., J k ( ω ; t ) ≈ β k ( ω )( N ∗ h k )( ω ; t − t ,k − U ( ω )). Then we obtain R k ( ω ; t ) ≈ Y k ( ω ; t ) − J k ( ω ; t ) , for all k = 1 , , · · · , K and ω ∈ Ω . (3.12)Given the estimated reference terms in (3.12), we may apply Algorithm 1 to estimate ptFC. The estimationprocedure provided above is summarized in the AMUSE-ptFC algorithm (Algorithm 2) for the estimationof ptFC in cases where reference terms are not available.Although the AMUSE-ptFCE algorithm is based only on the task-evoked signals and omits β k ( ω ) C k terms, its performance is comparable with that of ptFCE algorithm, which applies both task-evoked andreference signals, if random noise is large as illustrated in Section 4. While it may seem counterintuitivethat large noise helps improve the estimation accuracy of AMUSE-ptFCE, a possible explanation of thisphenomenon is that estimation bias from omitting β k ( ω ) C k is overwhelmed by large random noise. The R code for performing estimation using the ptFCE and AMUSE-ptFCE algorithms is available at https://github.com/KMengBrown/Population-level-Task-evoked-Functional-Connectivity.git .15 lgorithm 2 AMUSE-ptFCE Algorithm
Input: (i) Task-evoked signals { YYY ( ω j ; τ × T R ) = ( Y k ( ω j ; τ × T R ) , Y l ( ω j ; τ × T R )) T } Tτ =0 of participants { ω j } Mj =1 in population Ω; (ii) repetition time T R ;(iii) the stimulus signal { N ( τ × T R ) } Tτ =0 corresponding to the task of interest;(iv) reaction delay times t ,k (cid:48) = τ ,k (cid:48) × T R with some τ ,k (cid:48) ∈ Z for k (cid:48) ∈ { k, l } , and their default valuesare t ,k (cid:48) = 0 for k (cid:48) ∈ { k, l } ;(v) HRF h k (cid:48) for k (cid:48) ∈ { k, l } , and the default HRF are h k = h l = the R function canonicalHRF with itsdefault parameters. Output:
An estimation of the ptFC | corr ( β k , β l ) | between the k th and l th nodes. Generate M integers { u j } Mj =1 from the uniform distribution on { , , · · · , T } . for k (cid:48) ∈ { k, l } do , Compute C k (cid:48) ← M (cid:80) Mj =1 N ∗ h k (cid:48) (( τ − τ ,k (cid:48) − u j ) × T R ). for j = 1 , , · · · , M do Apply the AMUSE Algorithm using the following 2-dimensional signal as an input (cid:110)(cid:0) Y k (cid:48) ( ω j ; ( τ − u j ) × T R ) , N ∗ h k (cid:48) (cid:0) ( τ − τ ,k (cid:48) − u j ) × RT (cid:1) − C k (cid:48) (cid:1) T (cid:12)(cid:12) τ = 0 , , · · · , T (cid:111) . The AMUSE algorithm returns a 2-dimensional signal { ( s ( τ × T R ) , s ( τ × T R )) T | τ = 0 , , · · · , T } andan unmixing matrix QQQ . Let
MMM = ( m µν ) ≤ µ,ν ≤ be the (Moore-Penrose generalized) inverse of QQQ . for i = 1 , do Compute the Pearson correlation between { s i ( τ × T R ) } Tτ =0 and { ( N ∗ h k (cid:48) )(( τ − τ ,k (cid:48) − u j ) × T R ) − C k (cid:48) } Tτ =0 and denote the absolute value of this correlation as r i . end for Compute i ∗ = argmax l { r i | i = 1 , } . Update Y k (cid:48) ( ω j ; τ × T R ) ← m i ∗ × s i ∗ ( τ × T R ) and R k (cid:48) ( ω j ; τ × T R ) ← τ = 0 , , · · · , T . end for end for Apply Algorithm 1 using updated task-evoked signals { ( Y k ( ω j ; τ × T R ) , Y l ( ω j ; τ × T R )) T } τ,j , newlydefined reference signals { ( R k ( ω j ; τ × T R ) , R k ( ω j ; τ × T R )) T } τ,j , and the repetition time T R . Theestimate of ptFC | corr ( β k , β l ) | is returned. In this section, we first describe the data generating mechanism for simulations, which mimics the trends weobserved in the HCP motor-task dataset. Then, using the generated data, we illustrate the performance ofthe ptFCE and AMUSE-ptFCE algorithms in terms of estimating ptFC as compared to existing functionalconnectivity estimation approaches.
We simulate 4-dimensional synthetic BOLD signals { ( Y k ( ω j ; τ × T R ) , Y l ( ω j ; τ × T R ) , R k ( ω j ; τ × T R ) , R l ( ω j ; τ × T R )) } Tτ =0 , for j = 1 , , · · · , T R = 0 .
72 (seconds), T = 283; (ii) the stimulus signals in (2.2) and (2.4) applied in our simulationsare those used in HCP: N ( t ) = [86 . , . ( t ) + [162 , ( t ), ˜ N ( t ) = [71 . , . ( t ) + [177 . , . ( t ),˜ N ( t ) = [11 , ( t ) + [116 . , . ( t ), ˜ N ( t ) = [26 . , . ( t ) + [146 . , . ( t ), ˜ N ( t ) = [56 . , . ( t ) + [131 . , . ( t ); (iii) HRF h k and ˜ h k,γ , for γ ∈ { , , , } , are the double-gamma variate functions imple-mented in the R function canonicalHRF with default parameters; (iv) HRF h l and ˜ h l,γ , for γ ∈ { , , , } ,are the same function with parameters a1 = 10, a2 = 15, b1 = 0 . b2 = 0 .
9, and c = 0 .
35. We generate (cid:110) ( β k ( ω j ) , β l ( ω j )) T (cid:111) j =1 ∼ iid M V N ( µµµ β , VVV β ) and (cid:26)(cid:16) ˜ β k,γ ( ω j ) , ˜ β k,γ ( ω j ) (cid:17) T (cid:27) j =1 ∼ iid M V N (cid:16) µµµ ( γ ) , VVV ( γ ) (cid:17) , for γ ∈ { , , , } , where the collections of random vectors { ( ˜ β k,γ ( ω j ) , ˜ β k,γ ( ω j )) } j =1 , for γ ∈ { , , , } , aremutually independent and are independent of the collection { ( β k ( ω j ) , β l ( ω j )) T } j =1 , and µµµ β = µµµ ( γ ) = (0 , T , VVV ( γ ) = . × √ . × √ , for γ = 1 , , , , VVV β = √ ρ √ ρ , for ρ = 0 . , . , . . (4.1)The values ρ in (4.1) are the true values of ptFC in this simulation study. For each j ∈ { , , · · · , } and k (cid:48) ∈ { k, l } , the simulated random noise { (cid:15) k (cid:48) ( ω j ; τ × T R ) } Tτ =0 in (2.2) are iid generated from M V N ( µµµ n , VVV n ),where the random noise signals for different ω j are independent of each other and µµµ n = , VVV n = . . . (4.2)The BOLD signals in this simulation study are generated by the following two steps. Step 1 , for each j ∈ { , , · · · , } and k (cid:48) ∈ { k, l } , we generate { ˜ β k (cid:48) ,γ ( ω j ) | γ = 1 , , , } and { (cid:15) k (cid:48) ( ω j ; τ × T R ) } Tτ =0 using themechanism presented above, while the simulated reference signal { R k (cid:48) ( ω j ; τ × T R ) } Tτ =0 is defined by R k (cid:48) ( ω j ; τ × T R ) = 9000 + (cid:88) γ =1 ˜ β k (cid:48) ,γ ( ω j ) × (cid:16) ˜ N γ ∗ ˜ h k,γ (cid:17) ( τ × T R ) + (cid:15) k (cid:48) ( ω j ; τ × T R ) . (4.3) Step 2 , we generate another pair of realizations of { ˜ β k (cid:48) ,γ ( ω j ) | γ = 1 , , , } and { (cid:15) k (cid:48) ( ω j ; τ × T R ) } Tτ =0 andcompute another realization of the reference signal { R k (cid:48) ( ω j ; τ × T R ) } Tτ =0 using (4.3), then we generate17 ( β k ( ω j ) , β l ( ω j )) T } j =1 and compute the simulated task-evoked signals under the stimulus N ( t ) as follows. Y k (cid:48) ( ω j ; τ × T R ) = β k (cid:48) ( ω j ) × ( N ∗ h k (cid:48) ) ( τ × T R ) + R k (cid:48) ( ω j ; τ × T R ) for τ = 0 , , · · · , T and k (cid:48) ∈ { k, l } . These two simulation steps can be interpreted as follows: the experiment for each participant ω j has twoparts; in the first, participant ω j is at rest and reference signals R k (cid:48) ( ω j ; t ) are observed; in the second part,the participant perform stimulus N ( t ), and task-evoked signals Y k (cid:48) ( ω j ; t ) are observed. Although the signalsobtained during the second part of the experiment contain reference signals, these reference signals areindependent of the reference signals observed in the first part of the experiment. The upper panel of Figure4 presents a pair of the synthetic task-evoked signals { Y k ( ω j ; t ) , Y l ( ω j ; t ) } t ∈T for a participant ω j . In thesesimulated signals, we view N ( t ) as the task stimuli of interest and ˜ N γ ( t ), for γ = 1 , , ,
4, as stimuli not ofinterest, again following the set-up of the HCP motor-task study dataset analyzed in this paper.
Time B O L D Time B O L D B O L D Figure 4: In all three panels, the blue and red curves represent the synthetic task-evoked signals at the k th and l th nodes, respectively. Upper : A pair of synthetic task-evoked signals simulated using the datagenerating mechanism in Section 4.1.
Middle : Synthetic task-evoked signals simulated with λ = 2 in (4.4). Lower : Synthetic task-evoked signals simulated with λ = 5 in (4.4).18 .2 Performance of the ptFCE and AMUSE-ptFCE Algorithms Using the simulation mechanism in Section 4.1, we simulate { ( β k ( ω j ) , β l ( ω j )) } j =1 and the correspondingtask-evoked and reference signals. Then we compute the sample correlation between { β k ( ω j ) } j =1 and { β l ( ω j ) } j =1 across all ω j . The absolute value of this correlation is denoted as ˆ ρ true . Lastly, we apply theptFCE algorithm (Algorithm 1) to estimate ptFC using these task-evoked and reference signals and applythe AMUSE-ptFCE algorithm (Algorithm 2) to estimate ptFC from the task-evoked signals. The estimatesfrom Algorithms 1 and 2 are denoted by ˆ ρ and ˆ ρ , respectively. Repeating the simulation procedure 100times, we obtain 100 triplets ( ˆ ρ true , ˆ ρ , ˆ ρ ). The differences ( ˆ ρ − ˆ ρ true ) and ( ˆ ρ − ˆ ρ true ) measure the accuracyof Algorithms 1 and 2. The summary of these differences is presented in Table 1.Table 1: A summary of ( ˆ ρ − ˆ ρ true ) and (ˆ ρ − ˆ ρ true ) ρ = 0 . ρ = 0 . ρ = 0 . { ˆ ρ − ˆ ρ true } { ˆ ρ − ˆ ρ true } { ˆ ρ − ˆ ρ true } -0.02745 -0.08420 -0.1461sd of { ˆ ρ − ˆ ρ true } { (cid:15) k (cid:48) ( ω j ; τ × T R ) } Tτ =0 is determined by the variance VVV n in(4.2). To illustrate the influence of VVV n on the accuracy of the ptFCE and AMUSE-ptFCE algorithms, wereplace the VVV β and VVV n in (4.1, 4.2) with the following matrices. VVV β = √ × . √ × . , VVV n = λ . . , where λ ∈ (0 . , . (4.4)Examples of the task-evoked signals simulated using this mechanism with λ = 2 and λ = 5 are presented inthe middle and lower panels of Figure 4, respectively. For each fixed λ ∈ (0 . , ρ true , ˆ ρ , ˆ ρ ).Then the mean and sd of the differences ( ˆ ρ j − ˆ ρ true ), for j = 1 ,
2, across all 100 simulated populationsare computed. The relationship between ( ˆ ρ − ˆ ρ true ), ( ˆ ρ − ˆ ρ true ), and the noise variance scale λ in thissimulation study is presented by Figure 5.From the analyses above, we conclude that when both task-evoked and reference signals are available, theptFCE algorithm is applicable and performs well in terms of estimating ptFC. Although it tends to slightlyoverestimate ptFC, the corresponding bias is small. Additionally, Figure 5 shows that, as the random noisevariance increases, the bias and variance of ptFCE estimates increase as well, though moderately. When the19 − . . . Lambda a l g1_p t F C − t r ue_p t F C − . − . . Lambda a l g2_p t F C − t r ue_p t F C Figure 5: In each of the two panels, the blue curve presents the mean of 100 differences ( ˆ ρ j − ˆ ρ true ) with j ∈ { , } , the red dashed curves present the 95% Wald confidence interval, i.e., mean ± . × sd , and theblack dotted line indicates the zero value on the vertical axis. Upper : The relationship between ( ˆ ρ − ˆ ρ true )and λ presenting the influence of random noise on the ptFCE algorithm. Lower : The relationship between( ˆ ρ − ˆ ρ true ) and λ presenting the influence of random noise on the AMUSE-ptFCE algorithm.random noise variance is large, the AMUSE-ptFCE algorithm performs well in estimating ptFC. Althoughit tends to underestimate ptFC, the corresponding bias is small when the random noise variance is large.Next, we compare our proposed algorithms in terms of estimating pfFC with the following existing ap-proaches: beta-series regression, naive Pearson correlation , task Pearson correlation , and coherence analysis.We apply the simulation mechanism in Section 4.1 with VVV β defined by (4.1) except ρ ∈ { , . , . , . , } and VVV n defined by (4.2). For each ρ , the ptFC estimated using the ptFCE algorithm is denoted by ˆ ρ ptF CE ,while the estimated ptFC by the AMUSE-ptFCE algorithm is denoted by ˆ ρ ptF CE . In the implementationof the AMUSE-ptFCE algorithm, the input reaction delay t ,k (cid:48) and HRF h k (cid:48) are set to be correspondingdefault values in Algorithm 2. The existing methods are implemented as follows. Beta-series regression:
We apply the procedure described in Chapter 9.2 of Ashby (2019) to estimatetask-evoked FC using beta-series regression, except we ignore the “nuisance term” therein. For each ρ ∈{ , . , . , . , } and ω j , we apply beta-series regression to task-evoked signals { Y k ( ω j ; τ × T R ) , Y l ( ω j ; τ × T R ) | τ = 0 , , · · · , T } , estimate the FC evoked by N ( t ), and denote the estimated quantity by ˆ ρ betaSj . Thenwe compute the median of { ˆ ρ betaSj } j =1 across all j and denote it by ˆ ρ betaSM .20 aive Pearson correlation: For each ρ and ω j , we compute the Pearson correlation between { Y k ( ω j ; τ × T R ) } τ and { Y l ( ω j ; τ × T R ) } τ across all τ = 0 , , · · · , T and denote the absolute value of this correlation byˆ ρ naiveCorrj . Let ˆ ρ naiveCorrM denote the median of { ˆ ρ naiveCorrj } j =1 across all j . Task Pearson correlation:
For each ρ and ω j , we compute the Pearson correlation between { Y k ( ω j ; τ × T R ) | N ( τ × T R ) = 1 } and { Y l ( ω j ; τ × T R ) | N ( τ × T R ) = 1 } across all τ such that N ( τ × T R ) = 1 anddenote the absolute value of this correlation by ˆ ρ taskCorrj . We compute the median of { ˆ ρ taskCorrj } j =1 acrossall j and denote it by ˆ ρ taskCorrM . Coherence analysis:
For each ρ and ω j , we compute the coherence between task-evoked signals { Y k ( ω j ; τ × T R ) } τ and { Y l ( ω j ; τ × T R ) } τ by the R function coh in package seewave . The coherence is a function coh ρ,j ( ξ )of ξ ∈ (0 , / (2 × T R )). Since HRF acts as a band-pass filter (typically 0 − .
15 Hz, Aguirre et al. (1997)),we compute the median of { coh ρ,j ( ξ ) | < ξ < . } across all ξ ∈ (0 , .
15) and denote it by ˆ ρ Cohj . Finally,the median of { ˆ ρ Cohj } j =1 across all j is denoted by ˆ ρ CohM .For each ρ ∈ { , . , . , . , } , we repeat the procedures above for 100 simulated data sets and obtain100 estimated vectors ( ˆ ρ true , ˆ ρ ptF CE , ˆ ρ ptF CE , ˆ ρ betaSM , ˆ ρ naiveCorrM , ˆ ρ taskCorrM , ˆ ρ CohM ) . Instead of considering thespecific values of these measurements, we are interested in grading the strength of the pairwise connectivity,e.g., the FC between two nodes is either extremely strong , strong , median , weak , or extremely weak . Inaddition, these different approaches potentially use slightly different numerical scales for measuring FCstrength. Analogously, whether the Fahrenheit or Celsius scale system is used, we would conclude thatDeath Valley and Antarctica are the hottest and coldest places on earth, respectively, even though thesetwo systems are using different numerical measurements. Therefore, we compare the performance of theseapproaches by their “graded version” instead of considering the numerical values. Specifically, the fivevalues ρ ∈ { , . , . , . , } correspond to five grades: 0 (cid:55)→ grade 1 (extremely weak), 0 . (cid:55)→ grade 2(weak), 0 . (cid:55)→ grade 3 (median), 0 . (cid:55)→ grade 4 (strong), 1 (cid:55)→ grade 5 (extremely strong), i.e., ρ (cid:55)→ grade z = order ( ρ ), where order ( ρ ) denotes the z th smallest number ρ in { , . , . , . , } ; if an estimationapproach results in a correct estimate of grade order ( ρ ), we conclude that this approach is accurate, e.g.,the ptFCE algorithm is accurate if order ( ˆ ρ ptF CE ) = order ( ρ ) for all ρ . The accuracy of an approach ismeasured by its “correct grading rate,” i.e., the rate of order ( ˆ ρ ptF CE ) = order ( ρ ) across all 100 simulations.A summary of the accuracy of all estimation approaches described above is presented in Table 2.The ptFCE and AMUSE-ptFCE algorithms perform similarly well in grading ptFC, although the AMUSE-ptFCE algorithm is based on task-evoked signals only. The naive and task Pearson correlation approachesperform well in grading ptFC under the framework provided by Definition 2.1, although not as good as the21able 2: Performance in the estimation of order ( ρ ) ρ = 0 ρ = 0 . ρ = 0 . ρ = 0 . ρ = 1Rate of order ( ˆ ρ ptF CE ) = order ( ρ ) 94% 94% 100% 100% 100%Rate of order ( ˆ ρ ptF CE ) = order ( ρ ) 99% 99% 99% 98% 98%Rate of order ( ˆ ρ betaSM ) = order ( ρ ) 12% 18% 26% 21% 17%Rate of order ( ˆ ρ naiveCorrM ) = order ( ρ ) 86% 74% 76% 86% 92%Rate of order ( ˆ ρ taskCorrM ) = order ( ρ ) 96% 91% 92% 96% 99%Rate of order ( ˆ ρ CohM ) = order ( ρ ) 52% 34% 32% 36% 65%ptFCE and AMUSE-ptFCE algorithms. Under the framework of Definition 2.1, beta-series regression andcoherence analysis do not provide informative grading results.In all the simulation studies, using a PC with a processor and
32 GB2400 MHz DDR4 memory, all estimation approaches take less than 30 seconds to compute the correspondingestimate for each simulated data set.
In this section, we present results obtained by the estimation of ptFC using a data set in the task-evokedfunctional MRI component of HCP. For comparison, we apply the AMUSE-ptFCE algorithm and existingmethods to measure FC in the database of 308 participants from HCP performing motor tasks. Thedetails of the experimental design and corresponding parameters are available on the HCP website and inBarch et al. (2013). As described in Section 1, the participants in this experiment were asked to performfive tasks: tap their left/right fingers, squeeze their left/right toes, and move their tongue. We modelsqueezing right toes as the task of interest with two 12-second task blocks and the corresponding taskstimulus signal is N rt ( t ) = [86 . , . ( t ) + [162 , ( t ). The stimulus signals of the four nuisance stimuli arerepresented as ˜ N ( t ) = [71 . , . ( t ) + [177 . , . ( t ), ˜ N ( t ) = [11 , ( t ) + [116 . , . ( t ) , ˜ N ( t ) = [26 . , . ( t ) + [146 . , . ( t ) , and ˜ N ( t ) = [56 . , . ( t ) + [131 . , . ( t ) . The onsets of all task blocksvary across participants. Nevertheless, the corresponding onsets of any two participants differ in less than0.1 seconds. Therefore, we assume that all participants share the same stimulus signals.Before estimating FC, we compute region-specific time courses using the AAL atlas (Tzourio-Mazoyeret al. (2002)) that consists of 120 brain regions. For each region, we extract the voxel-specific time series inthat region and compute their spatial average for each time point. As a result, we obtain 120 time coursescorresponding to the 120 regions of interest. We select the left precentral gyrus (
PreCG.L ) as the seed regionsince it is located in the primary motor cortex and the motions of the right toe are associated with the22eft hemisphere of the brain. We measure the FC induced by N ( t ) at a population level between the seedregion and all other regions using the following five approaches: the AMUSE-ptFCE algorithm, naive andtask Pearson correlation methods, beta-series regression, and coherence analysis. † The detailed proceduresof applying these approaches are described in Section 4.2. . . . . . . . C onne c t i v i t y S F G . L M F G . L I F G ope r c . L I F G t r i ang . L I F G o r b . L R O L . L S M A . L O L F . L S F G m ed i a l . L P F C v en t m ed . L R E C . L O F C m ed . L O F C an t. L O F C po s t. L O F C l a t. L I N S . L A CC . L M CC . L P CC . L H I P . L P H G . L A M Y G . L C A L . L CUN . LL I N G . L S OG . L M OG . L I OG . L FF G . L P o C G . L SP G . L I P G . L S M G . L A N G . L P CUN . L P C L . L C A U . L P U T . L PA L . L T H A . L H ES . L S T G . L T P O s up . L M T G . L T P O m i d . L I T G . L C r u s . L C r u s . L C B . L C B . L C B . L C B . L C B . L V e r m i s . L V e r m i s . L V e r m i s . L V e r m i s . L P r e C G . R S F G . R M F G . R I F G ope r c . R I F G t r i ang . R I F G o r b . RR O L . R S M A . R O L F . R S F G m ed i a l . R P F C v en t m ed . RR E C . R O F C m ed . R O F C an t. R O F C po s t. R O F C l a t. R I N S . R A CC . R M CC . R P CC . RH I P . R P H G . R A M Y G . RC A L . RCUN . R L I N G . R S OG . R M OG . R I OG . R FF G . R P o C G . R SP G . R I P G . R S M G . R A N G . R P CUN . R P C L . RC A U . R P U T . R PA L . R T H A . RH ES . R S T G . R T P O s up . R M T G . R T P O m i d . R I T G . RC r u s . RC r u s . RC B . RC B . RC B . RC B . RC B . RC B C . R V e r m i s . R V e r m i s . R V e r m i s . R V e r m i s . R AMUSE−ptFCE algorithmNaive Pearson correlationTask Pearson correlationBeta−series regressionCoherence analysis
Figure 6: Illustration of estimation results from five FC estimation methods. The horizontal axis indicatesthe 116 regions. The abbreviations of region names are provided in the data set aal2.120 in the R package brainGraph . The vertical axis presents the standardized connectivity measurements between each regionand the seed region PreCG.R .Suppose X rj , for r = 1 , , · · · ,
116 and j = 1 , · · · ,
5, are the estimated FC between
PreCG.L and 116regions † computed using the AMUSE-ptFCE algorithm, naive Pearson correlation, task Pearson correlation,beta-series regression, and coherence analysis indexed by j . Since these estimation approaches use slightlydifferent scales, we standardize the estimates to enable comparisons as follows X ( st ) r,j = X r,j − min {X r (cid:48) ,j } r (cid:48) =1 max {X r (cid:48) ,j } r (cid:48) =1 − min {X r (cid:48) ,j } r (cid:48) =1 ,for r = 1 , , · · · ,
116 and j = 1 , · · · ,
5. The standardized versions of the estimates are presented in Figure6. As illustrated in Figure 6, the five methods result in similar estimated FC patterns. Specifically, theestimated FC tend to be relatively large/small simultaneously in most regions when comparing between thefive methods. Therefore, all five methods provide approximately the same population-level FC structure in-dicating that our proposed estimation approach has similar performance to widely used Pearson correlation,beta-series regression, and coherence analysis approaches. In the meantime, our proposed AMUSE-ptFCEAalgorithm is based on the model proposed in Section 2 providing a clear interpretation in applications andhas several advantages when compared to the competitor methods as described in Section 2.2. In addition,based on our proposed AMUSE-ptFCE algorithm estimates for FC, we identify high connectivity betweenseveral regions and the seed region that are missed by the competitor methods. Specifically, we obtain highFC between the seed region and left/right thalamus (passing motor signals to the cerebral cortex), left/rightparacentral lobule (motor nerve supply to the lower extremities), left superior temporal gyrus (containing † Because of numerical issues, we omit the regions
CB7.L , CB7.R , and
CB10.L and investigate the rest 117 regions.
PreCG.L
PreCG.RSFG.L SFG.RMFG.L MFG.RIFGoperc.L IFGoperc.RIFGtriang.L IFGtriang.RIFGorb .L IFGorb .RROL.L ROL.RSMA.L SMA.ROLF.L OLF.RSFGmedial.LSFGmedial.RPFCventmed.LPFCventmed.RREC.L REC.ROFCmed.L OFCmed.ROFCant.L OFCant.ROFCpost.L OFCpost.ROFClat.L OFClat.RINS.L INS.RACC.L ACC.RMCC.L MCC.RPCC.L PCC.RHIP.L HIP.RPHG.L PHG.RAMYG.L AMYG.RCAL.L CAL.RCUN.L CUN.RLING.L LING.RSOG.L SOG.RMOG.L MOG.RIOG.L IOG.RFFG.L FFG.RPoCG.L PoCG.RSPG.L SPG.RIPG.L IPG.RSMG.L SMG.RANG.L ANG.RPCUN.L PCUN.RPCL.L PCL.RCAU.L CAU.RPUT.L PUT.RPAL.L PAL.RTHA.L THA.RHES.L HES.RSTG.L STG.RTPOsup.L TPOsup.RMTG.L MTG.RTPOmid.L TPOmid.RITG.L ITG.RCrus1.L Crus1.RCrus2.L Crus2.RCB3.L CB3.RCB4.L CB4.RCB6.L CB6.RCB8.L CB8.RCB9.L CB9.R CBC10.RVermis1_2.LVermis3.RVermis4_5.LVermis6.RVermis7.LVermis8.RVermis9.LVermis10.R
PreCG.L
SFG.L MFG.L IFGoperc.L IFGtriang.L IFGorb .L ROL.L SMA.L SMA.R INS.L MCC.L MCC.R FFG.L PoCG.L SPG.L SPG.RIPG.L SMG.L SMG.R PCUN.L PCL.L PCL.RCAU.L PUT.L THA.L HES.L STG.L TPOsup.L MTG.L ITG.L
Figure 7: Each dot denotes a region from the AAL atlas, located using its corresponding MNI coordinates.The abbreviated region names are given next to each dot. We apply the AMUSE-ptFCE algorithm toestimate the ptFC between region
PreCG.L and each of the rest 116 regions. In the left panel, we usegrayscale coloring of the edges to indicate the magnitude of ptFC between the corresponding two vertices;specifically, the larger a ptFC, the darker the line segment connecting the corresponding region pair. In theright panel, the presented blue line segments indicate the first 30 largest pfFC estimated by the AMUSE-ptFCE algorithm among all 116 regions.To further visualize the ptFC induced by the task of interest, we apply the MNI space coordinates ofthe 117 regions from the AAL atlas, where the three-dimensional coordinates are obtained from the dataset aal2.120 in R package brainGraph . The regions are depicted by their MNI coordinates in Figure 7.The grayscale shade of the edges connecting each region to the seed region illustrates the magnitude of FCbetween the corresponding region and the seed region (Figure 7, left panel), while the edges in the rightpanel of Figure 7 present the highest FC values. These plots show that most of the large estimated ptFCvalues are in the left brain. This is expected, since it is known that behaviors of extremities are functionallyassociated with contralateral brain regions (Nieuwenhuys et al., 2014).24 Conclusions and Further Discussions
This paper proposes a model of the random effect form for task-evoked BOLD signals in fMRI studies. Basedon this model, we define a rigorous measurement of the task-evoked functional connectivity at a populationlevel. The newly defined ptFC framework has a clear interpretation in applications. We propose twoalgorithms, i.e., the ptFCE and AMUSE-ptFCE algorithms, for the estimation of ptFC. Simulation studiesshow that these two algorithms perform well in estimating ptFC. We apply the AMUSE-ptFCE algorithmto estimate FC in a motor-task dataset publicly available from the HCP. The data analysis results show thatthe AMUSE-ptFCE algorithm results in similar FC patterns to the widely used Pearson correlation, beta-series regression, and coherence analysis approaches. Additionally, the AMUSE-ptFCE algorithm indicatesthat the FC induced by squeezing right toes is mainly restricted in the left brain regions as expected.After estimating ptFC κ , an immediate question is to decide which pairs of nodes are functionallyconnected. In many applications, estimated FC values are compared to a pre-defined cut-off to identifyfunctionally connected nodes. This question can be presented by the hypothesis test H : κ ≥ κ ∗ vs. H : κ < κ ∗ , where κ ∗ is a predetermined threshold. If H is rejected, the corresponding pair of nodes arenot functionally connected; otherwise, they are functionally connected. We plan to develop a rigorous andcomputationally efficient approach for performing this hypothesis test. If a threshold κ ∗ is given, a boostrapapproach can be implemented to perform the hypothesis test. In most applications, the determination of athreshold κ ∗ for ptFC can be ambiguous. Cross-validation methods may be considered for determining thethreshold. Instead of choosing a single threshold κ ∗ and applying hypothesis testing, an alternative approachfor describing the FC structure using the estimated ptFC is to test all possible thresholds κ ∗ ∈ (0 , persistent homology for modeling all thresholds. Future researchmay consider combination of persistent homology theory and the proposed ptFC. The project herein was supported by Grant Number 5P20GM103645 from the National Institute of GeneralMedical Sciences. 25
Appendix
Lemma 8.1. If f : T → R is not constant zero, its Fourier transform (cid:98) f ( ξ ) has at most finitely many zeropoints - ξ ∈ R such that (cid:98) f ( ξ ) = 0 - in any compact subset of R .Proof. Define the complex function Φ f ( z ) := T +1 (cid:80) Tτ =0 f ( τ × T R ) e − πz ( τ × T R ) , for all z = η + iξ ∈ C . Sinceit is straightforward that Φ f ( z ) satisfies the Cauchy-Riemann equation ∂∂z Φ f ( z ) = 0, where ∂∂z = ( ∂∂η + i ∂∂ξ )is a Wirtinger derivative, the Looman–Menchoff theorem implies that Φ f ( z ) is a holomorphic function. Thenthe zero points of Φ f ( z ) are isolated, i.e., every zero point has a neighbourhood that does not contain anyother zero point. Therefore, (cid:98) f ( ξ ) = Φ f ( iξ ) implies the desired result. (cid:3) Theorem 8.1.
Suppose signals { Y k ( ω ; t ) | t ∈ T } Kk =1 , for ω ∈ Ω , are defined as in (2.2), E β k = E R k ( t ) = 0 ,and t ,k = τ ,k × T R with some τ ,k ∈ Z for all k = 1 , , · · · , K and t ∈ T . Let the random variable U : Ω → T be uniformly distributed on T . Furthermore, we assume that U ( ω ) , { β k ( ω ) } Kk =1 , and { R k ( ω ; t ) | t ∈ T } Kk =1 are independent. Then, the autocovariance differences in (3.5) depend only on s , the Fourier transforms of A ( s ) defined in (3.5) are (cid:98) A kl ( ξ ) = E ( β k β l ) × | (cid:98) N ( ξ ) | (cid:98) h k ( ξ ) (cid:98) h l ( ξ ) e πi ( t ,k − t ,l ) ξ , and C kl ( ξ ) = | corr ( β k , β l ) | forall ξ ∈ R , where C kl ( ξ ) is defined in (3.6).Proof. The independence between U ( ω ), { β k ( ω ) } Kk =1 , and { R k ( ω ; t ) | t ∈ T } Kk =1 implies E [ Y k ( t − U ) Y l ( t + s − U )] − E [ R k ( t − U ) R l ( t + s − U )]= E ( β k β l ) × E [( N ∗ h k )( t − t ,k − U ) × ( N ∗ h l )( t + s − t ,l − U )]= E ( β k β l ) × T + 1 T (cid:88) u =0 (cid:110) ( N ∗ h k ) (( τ − τ ,k − u ) × T R ) × ( N ∗ h l ) (( s + τ ,k − τ ,l ) × T R + ( τ − τ ,k − u ) × T R ) (cid:111) = E ( β k β l ) × T + 1 T − ( τ − τ ,k ) (cid:88) v = − ( τ − τ ,k ) (cid:110) ( N ∗ h k ) ( − v × T R ) × ( N ∗ h l ) (( s + τ ,k − τ ,l ) × T R − v × T R ) (cid:111) = E ( β k β l ) × [ N ∗ h k ( −· )] ∗ [ N ∗ h l ] ( s + ( t ,k − t ,l )) , (8.1)which only depends on s and does not depend on t , where the last equality follows from the periodic extension(3.1) and the definition of convolution. Then, we have the following Fourier transform (cid:100) A kl ( ξ ) = E ( β k β l ) × (cid:92) ( N ∗ h k )( ξ ) (cid:92) ( N ∗ h l )( ξ ) e πi ( t ,k − t ,l ) ξ = E ( β k β l ) × (cid:12)(cid:12)(cid:12) (cid:98) N ( ξ ) (cid:12)(cid:12)(cid:12) (cid:99) h k ( ξ ) (cid:98) h l ( ξ ) e πi ( t ,k − t ,l ) ξ , for all ξ ∈ R implying C kl ( ξ ) = | corr ( β k , β l ) | for all ξ ∈ R . (cid:3) emma 8.2. The vector-valued stochastic process
ZZZ ( ω ) = { ( Z ( ω ; t ) , Z ( ω ; t ) , · · · , Z K ( ω ; t )) T } t ∈T is definedon (Ω , P ) and weakly stationary with mean zero, U : Ω → T is a uniformly distributed, and U ( ω ) and ZZZ ( ω ) are independent. Then the vector-valued stochastic process { ( Z ( ω ; t − U ( ω )) , Z ( ω ; t − U ( ω )) , · · · , Z K ( ω ; t − U ( ω ))) T } t ∈T is weakly stationary with mean zero as well. Additionally, E [ Z k ( t − U ) Z l ( t + s − U )] = E [ Z k (0) Z l ( s )] , for all k, l = 1 , , · · · , K .Proof. Let P ZZZ = P ZZZ ( dz ) be the probability measure on the path space ( R T ) K associated with the K -vector-valued stochastic process ZZZ , i.e., P ZZZ = P ◦ ZZZ − . The probability measure on T associated with U is denotedas µ , i.e., µ = P ◦ U − . Since ZZZ and U are independent, the probability measure on ( R T ) K × T associatedwith the joint distribution ( ZZZ, U ) is the product measure P ZZZ ⊗ µ , i.e., P ZZZ ⊗ µ = P ◦ ( ZZZ, U ) − . Then the zeromean of Z k ( t ) implies E Z k ( t − U ) = (cid:82) T E [ Z k ( t − U ) | U = u ] µ ( du ) = (cid:82) T E Z k ( t − u ) µ ( du ) = 0. Let s ∈ T and k, l ∈ { , , · · · , K } be fixed, for each t ∗ ∈ T , we define the following map.Φ t ∗ : ( R T ) K → R , { ( z ( t ) , z ( t ) , · · · , z K ( t )) } t ∈T (cid:55)→ z k ( t ∗ ) z l ( t ∗ + s ) . Then Fubini’s theorem implies the following representation. E [ Z k ( t − U ) Z l ( t + s − U )] = E [Φ t − U ( ZZZ )] = (cid:90) ( R T ) K ×T Φ t − u ( z ) P ZZZ ⊗ µ ( dz, du ) = (cid:90) T (cid:34)(cid:90) ( R T ) K Φ t − u ( z ) P ZZZ ( dz ) (cid:35) µ ( du ) . Since
ZZZ is weakly stationary, (cid:82) ( R T ) K Φ t − u ( z ) P ZZZ ( dz ) = E [ Z k ( t − u ) Z l ( t + s − u )] = E [ Z k (0) Z l ( s )], for all u ∈ T . Then E [ Z k ( t − U ) Z l ( t + s − U )] = E [ Z k (0) Z l ( s )], which depends only on s . (cid:3) Theorem 8.2.
Suppose signals Y k ( ω ; t ) are defined in (3.7) for all k = 1 , , · · · , K , the random variable U :Ω → T is uniformly distributed, and the stochastic process { W ( ω ; t ) = ( W ( ω ; t ) , W ( ω ; t ) , · · · , W K ( ω ; t )) T } t ∈T satisfies (i) W ( ω ; t ) and W ( ω ; t ) are independent if t (cid:54) = t , (ii) Σ kl = E [ W k ( t ) W l ( t )] for all t ∈ T , and(iii) W ( ω ; t ) is weakly stationary with mean zero. Assuming { W k ( ω ; t ) | t ∈ T } Kk =1 , U ( ω ) , { β k ( ω ) } Kk =1 , and { R k ( ω ; t ) | t ∈ T } Kk =1 are independent, we have (3.8).Proof. That { W k ( ω ; t ) | t ∈ T } Kk =1 , U ( ω ), { β k ( ω ) } Kk =1 , and { R k ( ω ; t ) | t ∈ T } Kk =1 are independent implies A kl ( s ) = E ( β k β l ) × I + II , where I = T +1 (cid:80) Tu =0 N ∗ h k (( τ − τ ,k − u ) × T R ) N ∗ h l (( τ + s − τ ,l − u ) × T R )and II = E [ W k ( t − U ) W l ( t + s − U )]. The calculation in (8.1) implies that the Fourier transform of I is | (cid:98) N ( ξ ) | (cid:99) h k ( ξ ) (cid:98) h l ( ξ ) e πiξ ( t ,k − t ,l ) . Lemma 8.2 implies II = E [ W k (0) W l ( s )] = Σ kl × { } ( s ), where { } ( s ) inthe indicator function for the singleton { } , and the Fourier transform of II is the constant Σ kl / ( T + 1).Then (cid:98) A kl ( ξ ) = E ( β k β l ) (cid:12)(cid:12)(cid:12) (cid:98) N ( ξ ) (cid:12)(cid:12)(cid:12) (cid:99) h k ( ξ ) (cid:98) h l ( ξ ) e πiξ ( t ,k − t ,l ) + Σ kl T +1 , and (3.8) follows. (cid:3) roof of Theorem 3.1: Lemma 8.2 implies that, for each fixed k ∈ { , , · · · , K } , the scalar-valuedstochastic process { R k ( ω ; t − U ( ω )) } t ∈T is weakly stationary with mean zero. The independence between β k ( ω ) and { U ( ω ) , R k ( ω ; t ) | t ∈ T } implies that { R k ( ω ; t − U ( ω )) } t ∈T is weakly stationary with mean zeroconditioning on β k . Additionally, E { J k ( t ) R k ( t + s − U ) | β k } ( ω ) = β k ( ω ) E { [( N ∗ h k )( t − t ,k − U ) − C k ] E [ R k ( t + s − U ) | U ] } a.s.. Since E R k ( t ) = 0 for all t ∈ T and U ( ω ) is independent of { R k ( ω ; t ) } t ∈T , we have E [ R k ( t + s − U ) | U = u ] = E [ R k ( t + s − u )] = 0, for all u ∈ T . Then E { J k ( t ) R k ( t + s − U ) | β k } = 0, for all s and t , and the desiredresult follows from Theorem 2 of Tong et al. (1991). (cid:3) Proof of Theorem 3.2:
Since
PPP is a permutation matrix and ΛΛΛ is a diagonal matrix, the equation (3.10)implies β k ( ω ) = a i (cid:48) /a i (cid:48) for some i (cid:48) ∈ { , } and a j = 0 when j (cid:54) = i (cid:48) . Then (3.9) implies a i (cid:48) s i (cid:48) ( ω ; t ) =( N ∗ h k )( t − t ,k − U ( ω )) − C k . Therefore, a i (cid:48) s i (cid:48) ( ω ; t ) = a i (cid:48) a i (cid:48) × a i (cid:48) s i (cid:48) ( ω ; t ) = J k ( ω ; t ). (cid:3) References
G. K. Aguirre, E. Zarahn, and M. D’Esposito. Empirical analyses of bold fmri statistics.
NeuroImage , 5(3):199–212, 1997.F. G. Ashby.
Statistical analysis of fMRI data . MIT press, 2019.D. M. Barch, G. C. Burgess, M. P. Harms, S. E. Petersen, B. L. Schlaggar, M. Corbetta, M. F. Glasser,S. Curtiss, S. Dixit, C. Feldt, et al. Function in the human connectome: task-fmri and individual differencesin behavior.
Neuroimage , 80:169–189, 2013.P. Bloomfield.
Fourier analysis of time series: an introduction . John Wiley & Sons, 2004.R. L. Buckner, W. Koutstaal, D. L. Schacter, A. M. Dale, M. Rotte, and B. R. Rosen. Functional–anatomicstudy of episodic retrieval: Ii. selective averaging of event-related fmri trials to test the retrieval successhypothesis.
Neuroimage , 7(3):163–175, 1998.R. L. Buckner, F. M. Krienen, A. Castellanos, J. C. Diaz, and B. T. Yeo. The organization of the humancerebellum estimated by intrinsic functional connectivity.
Journal of neurophysiology , 106(5):2322–2345,2011.R. Durrett.
Probability: theory and examples , volume 49. Cambridge university press, 2019.28. Friston, C. Frith, P. Liddle, and R. Frackowiak. Functional connectivity: the principal-component analysisof large (pet) data sets.
Journal of Cerebral Blood Flow & Metabolism , 13(1):5–14, 1993.K. J. Friston. Functional and effective connectivity in neuroimaging: a synthesis.
Human brain mapping , 2(1-2):56–78, 1994.K. J. Friston. Functional and effective connectivity: a review.
Brain connectivity , 1(1):13–36, 2011.K. J. Friston, P. Jezzard, and R. Turner. Analysis of functional mri time-series.
Human brain mapping , 1(2):153–171, 1994.M. Hampson, B. S. Peterson, P. Skudlarski, J. C. Gatenby, and J. C. Gore. Detection of functional connec-tivity using temporal correlations in mr images.
Human brain mapping , 15(4):247–262, 2002.L. H¨ormander. Estimates for translation invariant operators inl p spaces.
Acta Mathematica , 104(1-2):93–140, 1960.S. E. Joel, B. S. Caffo, P. C. van Zijl, and J. J. Pekar. On the relationship between seed-based and ica-basedmeasures of functional connectivity.
Magnetic Resonance in Medicine , 66(3):644–657, 2011.H. Lee, M. K. Chung, H. Kang, B.-N. Kim, and D. S. Lee. Discriminative persistent homology of brainnetworks. In , pages841–844. IEEE, 2011.S.-P. Lee, T. Q. Duong, G. Yang, C. Iadecola, and S.-G. Kim. Relative changes of cerebral arterial and venousblood volumes during increased cerebral blood flow: implications for bold fmri.
Magnetic Resonance inMedicine: An Official Journal of the International Society for Magnetic Resonance in Medicine , 45(5):791–800, 2001.M. A. Lindquist et al. The statistical analysis of fmri data.
Statistical science , 23(4):439–464, 2008.M. J. Lowe, M. Dzemidzic, J. T. Lurito, V. P. Mathews, and M. D. Phillips. Correlations in low-frequencybold fluctuations reflect cortico-cortical connections.
Neuroimage , 12(5):582–587, 2000.F. M. Miezin, L. Maccotta, J. Ollinger, S. Petersen, and R. Buckner. Characterizing the hemodynamicresponse: effects of presentation rate, sampling procedure, and the possibility of ordering brain activitybased on relative timing.
Neuroimage , 11(6):735–759, 2000.29. M¨uller, G. Lohmann, V. Bosch, and D. Y. Von Cramon. On multivariate spectral analysis of fmri timeseries.
NeuroImage , 14(2):347–356, 2001.R. Nieuwenhuys, J. Hans, and C. Nicholson.
The central nervous system of vertebrates . Springer, 2014.J. Rissman, A. Gazzaley, and M. D’Esposito. Measuring functional connectivity during distinct stages of acognitive task.
Neuroimage , 23(2):752–763, 2004.B. R. Rosen, R. L. Buckner, and A. M. Dale. Event-related functional mri: past, present, and future.
Proceedings of the National Academy of Sciences , 95(3):773–780, 1998.Z. S. Saad, K. M. Ropella, R. W. Cox, and E. A. DeYoe. Analysis and use of fmri response delays.
Humanbrain mapping , 13(2):74–93, 2001.D. L. Schacter, R. L. Buckner, W. Koutstaal, A. M. Dale, and B. R. Rosen. Late onset of anterior prefrontalactivity during true and false recognition: an event-related fmri study.
Neuroimage , 6(4):259–269, 1997.F. T. Sun, L. M. Miller, and M. D’Esposito. Measuring interregional functional connectivity using coherenceand partial coherence analyses of fmri data.
Neuroimage , 21(2):647–658, 2004.L. Tong, R.-W. Liu, V. C. Soon, and Y.-F. Huang. Indeterminacy and identifiability of blind identification.
IEEE Transactions on circuits and systems , 38(5):499–509, 1991.N. Tzourio-Mazoyer, B. Landeau, D. Papathanassiou, F. Crivello, O. Etard, N. Delcroix, B. Mazoyer,and M. Joliot. Automated anatomical labeling of activations in spm using a macroscopic anatomicalparcellation of the mni mri single-subject brain.
Neuroimage , 15(1):273–289, 2002.R. Warnick, M. Guindani, E. Erhardt, E. Allen, V. Calhoun, and M. Vannucci. A bayesian approach forestimating dynamic functional network connectivity in fmri data.
Journal of the American StatisticalAssociation , 113(521):134–151, 2018.B. T. Yeo, F. M. Krienen, J. Sepulcre, M. R. Sabuncu, D. Lashkari, M. Hollinshead, J. L. Roffman, J. W.Smoller, L. Z¨ollei, J. R. Polimeni, et al. The organization of the human cerebral cortex estimated byintrinsic functional connectivity.