Detecting Task-Dependent Functional Connectivity in Group Iterative Multiple Model Estimation with Person-Specific Hemodynamic Response Functions

Abstract

Introduction:

Group iterative multiple model estimation (GIMME) has proven to be a reliable data-driven method to arrive at functional connectivity maps that represent associations between brain regions across time in groups and individuals. However, to date, GIMME has not been able to model time-varying task-related effects. This article introduces HRF-GIMME, an extension of GIMME that enables the modeling of the direct and modulatory effects of a task on functional magnetic resonance imaging data collected by using event-related designs. Critically, hemodynamic response function (HRF)-GIMME incorporates person-specific modeling of the HRF to accommodate known variability in onset delay and shape.

Methods:

After an introduction of the technical aspects of HRF-GIMME, the performance of HRF-GIMME is evaluated via both a simulation study and application to empirical data. The simulation study assesses the sensitivity and specificity of HRF-GIMME by using data simulated from one slow and two rapid event-related designs, and HRF-GIMME is then applied to two empirical data sets from similar designs to evaluate performance in recovering known neural circuitry.

Results:

HRF-GIMME showed high sensitivity and specificity across all simulated conditions, and it performed well in the recovery of expected relations between convolved task vectors and brain regions in both simulated and empirical data, particularly for the slow event-related design.

Conclusion:

Results from simulated and empirical data indicate that HRF-GIMME is a powerful new tool for obtaining directed functional connectivity maps of intrinsic and task-related connections that is able to uncover what is common across the sample as well as crucial individual-level path connections and estimates.

Impact statement

Group iterative multiple model estimation (GIMME) is a reliable method for creating functional connectivity maps of the connections between brain regions across time, and it is able to detect what is common across the sample and what is shared between subsets of participants, as well as individual-level path estimates. However, historically, GIMME does not model task-related effects. The novel HRF-GIMME algorithm enables the modeling of direct and modulatory task effects through individual-level estimation of the hemodynamic response function (HRF), presenting a powerful new tool for assessing task effects on functional connectivity networks in functional magnetic resonance imaging data.

Introduction

Neuroscientists have increasingly employed functional connectivity analyses to arrive at graphical maps that represent time-dependent associations among brain regions. Group Iterative Multiple Model Estimation (GIMME) (Gates and Molenaar, 2012) has emerged as a reliable data-driven method for arriving at functional connectivity maps (Mumford and Ramsey, 2014). GIMME performs well on benchmark simulated data sets (Gates and Molenaar, 2012; Smith et al., 2011) and across numerous types of data modalities (Gates et al., 2017; Lane et al., 2019). Though useful for data obtained from block designs and resting-state functional magnetic resonance imaging (fMRI), until now it was not possible to model task-based effects within the GIMME framework.

Connectivity patterns vary meaningfully across different experimental tasks (Cohen and D'Esposito, 2016; Gonzalez-Castillo et al., 2017) and task conditions (Friston et al., 1997). For instance, the functional coupling between two brain regions might increase when presented with task stimuli. Hence, making inferences regarding brain processes elicited by particular stimuli or task demands depends on modeling task-based connectivity across time. To accomplish this goal, we extend GIMME to enable modeling of the direct influence of tasks on brain region activity as well as how tasks modulate, or change, the relations among brain regions. Termed HRF-GIMME, this novel approach incorporates person-specific modeling of the hemodynamic response function (HRF) into the existing framework to obtain sparse, directed functional connectivity maps in data obtained from event-related task designs.

The known heterogeneity across individuals that exists in connectivity patterns (Finn et al., 2015; Laumann et al., 2015) necessitates balancing the discovery of generalizable findings (i.e., connections that likely exist for all individuals) and individual differences. From a statistical standpoint, the traditional approach in fMRI of aggregating across all individuals to arrive at one composite connectivity map may identify spurious connections and generate connectivity maps that fail to describe any individual in the sample (Beltz et al., 2018; Molenaar, 2004). The GIMME avoids this pitfall by allowing for individual nuances in connectivity patterns and not aggregating by averaging across individuals. Similar to multitask approaches in machine learning, the algorithm first arrives at a group-level pattern that represents the generalizable aspects of the connectivity map. For each individual, the model search then builds on this baseline group-level model to uncover the individual-specific connectivity patterns. In addition to arriving at valid generalizable connections, this approach greatly reduces false positives (i.e., spurious connections) and overfitting to the data when arriving at individual-level connections (Gates and Molenaar, 2012).

Hitherto, GIMME has only been evaluated and developed for use in arriving at intrinsic connectivity patterns, or connectivity patterns that represent the direct relations among brain regions believed to be constant across time (e.g., resting state or task block designs). HRF-GIMME additionally enables the discovery of task-based effects in the connectivity maps, which can be direct or modulating. Direct effects of a task on region of interest (ROI) activity can be understood as the degree to which a task influences variability in blood-oxygen-level-dependent (BOLD) activity. This is often what is investigated in statistical parametric maps to identify which regions have significant change in activity during a task or stimulus presentation. Modulating effects may also arise when the strength or direction of relations between two ROIs is influenced by the task stimuli. For instance, perhaps in the absence of a task there is no functional connection between two regions yet a connection surfaces during a task. Or, perhaps an intrinsic functional connection exists but becomes attenuated in the presence of a task. This article presents innovations that enable these direct and modulating task effects to be identified.

The multiplication of the HRF-convolved task vector with the brain region activity across time allows for investigation into how the connections between brain regions vary in the presence of the task, much like is done in psychophysiological interaction (PPI) analysis (Friston et al., 1997). The difference in HRF-GIMME is that intrinsic connections are also included, meaning that the relations among brain regions emerge that may exist regardless of the presence or absence of a task. A benefit of examining these relations in the context of functional connectivity is that one can identify indirect effects, or the possibility that the relation between task and BOLD activity for a specific region is mediated by another region. By modeling intrinsic, direct task, and modulatory effects simultaneously the inferences made can be similar to those in dynamic causal modeling (DCM) (Friston et al., 2003; Gates et al., 2011). These modulatory or bilinear effects in combination with the intrinsic connections among brain regions enable the use of the HRF-GIMME algorithm with task-based data, in which connectivity patterns and strengths likely change across time.

Modeling the task-related connections among brain regions typically requires the use of an HRF to arrive at the expected BOLD response generated by task stimuli. The HRF describes the expected time course of the changes in deoxyhemoglobin that occur after neural activity. Previous research has shown that the parameters governing the shape of the HRF vary between people more so than within people (Handwerker et al., 2004). As such, to account for HRF variability, best practice involves deriving person-specific HRF parameters (Friston et al., 1997). The HRF-GIMME uses the smoothed finite impulse response (sFIR) (Goutte et al., 2000) to estimate HRF parameters separately for each individual. The FIR makes no assumptions about the shape of the HRF (Lindquist and Wager, 2007) and has been found to correctly recover the true HRF shape even when it varies across individuals (Lindquist et al., 2009). By using the sFIR to model individual-level HRF, HRF-GIMME can flexibly control for potential direct effects of the task on each brain region's activity while also investigating connectivity among the brain regions themselves (Price et al., 2020).

The article is structured as follows. First, we introduce the technical details regarding the GIMME algorithm. Next, we describe the generation of person-specific HRF-convolved task vectors and their use in HRF-GIMME. We then evaluate the performance of HRF-GIMME for data generated to be from slow and fast event-related task designs. Given that modulatory effects have the potential to be low powered (McClelland and Judd, 1993), we also investigate the benefits and pitfalls of having a more lenient threshold for including connections in models. We then present results from empirical data examples. Finally, we conclude with summaries of the findings and suggestions for users.

Technical Description of HRF-GIMME

Original GIMME

The original GIMME algorithm (Gates and Molenaar, 2012) provides the basis of the current extension. GIMME obtains reliable group- and individual-level patterns of dynamic effects, with all effects being estimated uniquely for each individual. The primary heuristic behind GIMME is to search for relations that are consistent across individuals and to use those relations as a foundation for searching for relations that may be unique to individuals. GIMME is based on a unified structural equation modeling (uSEM) framework (Gates et al., 2010; Kim et al., 2007). It detects both contemporaneous (i.e., instantaneous) as well as lagged relations, both of which occur in fMRI data due to the rate of measurement relative to the underlying neuronal activity of interest (Logothetis, 2008). Simulation studies suggest that contemporaneous relations best reflect the brain processes when using fMRI data (Ramsey et al., 2011; Smith et al., 2011).

The uSEM implemented in GIMME is: $y_{i, t} = (A_{i}^{i} + A_{i}^{g}) y_{i, t} + (Φ_{i}^{i} + Φ_{i}^{g}) y_{i, t - 1} + 1_{t},$ (1)

where $y_{i, t}$ is a vector of the p-variate ROI time series for an individual i at time t, each A matrix is a p × p matrix of contemporaneous effect estimates among the p ROIs (with a zero diagonal), each $Φ$ is a p × p matrix of lagged effect estimates with autoregressive effects on the diagonal, and $1_{t}$ is the p-variate vector of the errors that are assumed to be uncorrelated white noise processes. The superscript g for the parameter matrices indicates that the matrix has the same structure across the group (here, the entire sample). The superscript i indicates patterns of connections that are unique for the individual. The subscript i indicates individual-level estimates for those connections identified to be estimated. Note that every parameter included in the individual- and group-level patterns is estimated uniquely for each individual, even when the pattern is the same, as seen in the g-superscript matrices.

The search for connectivity patterns begins by estimating a null model defined here as having no connections among the brain regions; they are considered independent. Typically, autoregressive effects, or how a brain region's activity predicts itself at the next time point, are estimated to start. Thereafter, a score function, or modification index (Sörbom, 1989), is calculated for each potential connection among brain regions for each individual. The GIMME algorithm first selects the connection that is significant for the most individuals, as indicated by the modification indices after a strict Bonferroni correction. If this connection is significant for at least the majority of individuals, then it is added. The threshold for what constitutes the majority can be user-defined, and we explore the impact of that choice in this article. The current default is 75% of the individuals in the sample, which reflects the expected missed-signal due to noise (Smith et al., 2011). The sequence continues in a forward-selection fashion until no modification indices meet this criterion. GIMME continues by estimating for each individual the connectivity pattern of effects obtained at the group-level search, which can include both contemporaneous and lagged relations. From this starting point it then searches for connections that are unique to individuals, again using modification indices to guide the addition of connections. Full details can be found in Gates and Molenaar (2012) and Lane and Gates (2017).

HRF-GIMME

HRF-GIMME extends GIMME to model the direct effects of tasks as well as how tasks modulate (or change the weight of) connections among brain regions. Task effects have previously been integrated into the individual-level uSEM in the form of the extended unified SEM (euSEM) (Gates et al., 2011), which for one task onset vector is: $\begin{matrix} y_{i, t} = (A_{i}^{i} + A_{i}^{g}) y_{i, t} + (Φ_{i}^{i} + Φ_{i}^{g}) y_{i, t - 1} + (γ_{i}^{i} + γ_{i}^{g}) u_{i, t} \\ + (τ_{i}^{i} + τ_{i}^{g}) u_{i, t} y_{i, t} + 1_{t}, \end{matrix}$ (2)

where subscript and superscript conventions are as explained earlier. The new $γ$ parameter matrices provide the estimated coefficients for the direct relation between the convolved task vector ( u , described below) and the timeseries for each brain region. The $τ$ parameter matrices provide the estimated bilinear coefficients, or the modulatory effects for how the relation between two brain regions varies based on values of the convolved task vector. This can be extended to include additional task vectors if the experiment involves multiple trial types or phases. When task-onset vectors and specific multiplied variables are included, the direct and modulatory task effects are selected by using the same GIMME algorithm described earlier.

Obtaining person-specific convolved task vectors

The first step in arriving at the person-specific convolved task vectors is to estimate the expected shape of the BOLD response given the task-onset vector. Of all ROIs included in the data passed into HRF-GIMME, HRF-GIMME estimates the HRF for each ROIs using sFIR, and it selects the estimated HRF that explains the greatest proportion of variance. Although it may be ideal to generate an HRF for each ROI within an individual, this is not feasible within the present estimation framework. Given this limitation and prior work (Handwerker et al., 2004), it is imperative that HRF-GIMME estimates HRFs separately for individuals, but not within each individual.

The FIR contains one free parameter for each time point in a window of time after the onset of a stimulus (Glover, 1999). Each time point after the onset of the stimulus is modeled as a separate weighted basis function. In this way, the response function is estimated for the provided stimulus onset throughout the provided length of the response window (e.g., 16 sec). The smooth FIR is simply an FIR model with a Gaussian prior placed on adjacent beta estimates to constrain the fit to be smoother (Goutte et al., 2000).

Should missing values exist for some time points, such as those that have been scrubbed due to motion, HRF-GIMME utilizes the Kalman filter to impute values before the FIR analysis by using the R package imputeTS (Moritz and Bartz-Beielstein, 2017). This imputation is only used to generate the convolved task vector.^* Once sFIR is used to obtain unique parameter estimates of the HRF for each individual, the onset task vector is convolved with this basis set. This results in a unique HRF-convolved task vector for each individual that matches the expected change in the BOLD signal after a task stimulus presentation (Fig. 1).

FIG. 1.

Schematic representation of convolution of task onset vector with HRF. (a) Demonstrates the process of estimating the person-specific HRF after the onset of a stimulus (represented as the red line) using the sFIR, done internally in HRF-GIMME. (b) Represents the stimuli in an event-related design, with each red line representing a stimulus. (c) Demonstrates the final convolved event-onset vector, where the person-specific HRF model estimated in (a) is overlaid with the events in (b). GIMME, group iterative multiple model estimation; HRF, hemodynamic response function; sFIR, smoothed finite impulse response. Color images are available online.

Simulation Study

Simulation conditions

We conducted a simulation study to assess whether HRF-GIMME could recover a known data-generating connectivity pattern under conditions likely to be encountered by fMRI researchers. The simulation study sought to evaluate two sets of conditions: (1) which group cutoff to use for defining “the majority” when arriving at the pattern of connections common across individuals, and (2) how well the algorithm performs on data obtained from slow versus rapid onset designs. We opted not to investigate the impact of altering the HRF parameters (e.g., slow time to peak), as the sFIR has already extensively been evaluated for this purpose (Lindquist et al., 2009). Provided that the HRF shape is estimated appropriately for each individual, which can be expected (Linquist et al., 2009), there is no reason to believe variations in the HRF such as delay to onset will impact results here. One might argue that magnitude differences may be important. By standardizing each resulting HRF vector as well as all ROI time series and multiplied variables, we allay this concern.

The default for what constitutes the majority in GIMME has been informed by a large-scale simulation study conducted by Smith and colleagues (2011) on data generated to emulate resting-state fMRI data. However, Smith and colleagues (2011) did not consider modulatory effects, which can be difficult to detect due to low power (McClelland and Judd, 1993). For this reason, we evaluated two criteria for the majority, 75% and 51%, to determine their performance in recovering simulated connectivity patterns.

The second simulation factor aimed at testing the ability of HRF-GIMME to detect changes in connectivity that may occur quickly. Given that most investigations into HRF performance (Lindquist et al., 2009) have evaluated the HRF estimation in slow event-related designs, we expected that HRF-GIMME would perform well when the inter-trial interval (ITI) was approximately as long as the expected length of the HRF (e.g., around 16 sec to recovery). This allows for the elicitation of task-related BOLD activity and changes in the connectivity patterns. For this simulation parameter, we generated data with an ITI of 30 sec. To model a more rapid event-related design, two conditions were used, one with a mean ITI of 2.42 sec, to emulate the design of the rapid event-related empirical data used, and the other with a mean ITI of 3.92 sec, to emulate designs suggested in the literature. Exact specification of the event designs is described in the subsequent section. It is suspected that HRF-GIMME may not perform as well in the rapid event-related conditions as the slow event-related condition, since the expected task-evoked response may not have enough variability to capture changes in connectivity between relations, as illustrated in Figure 2. We expect that task-modulated relations (i.e., connections that are strong enough to reach significance only during task) will surface as intrinsic connections in this case. In addition, having such short intervals between tasks may lead to violations of the HRF linearity assumption (Miezin et al., 2000).

FIG. 2.

Exemplar data of rapid (left) and slow (right) event-related onset vectors convolved with an HRF for 50 time points. Black lines indicate onset of task events.

Data generation procedure

Generating task vectors

For the rapid event onset condition of the simulations, binary event onset vectors were created by using OptSeq (Greve, 2002), a tool developed to automatically generate onsets for rapid event-related designs with optimal “jittering” of null events to lessen the overlap of the HRF estimates inherent in rapid stimulus presentation and optimize the power of the design. Two rapid conditions with an ITI of 2 sec were used: The first had 20% of trials as null events, with a mean ITI of 2.40 sec (min 2 sec, max 6 sec), to better mimic the rapid condition of the empirical data we used (Mennes et al., 2013), and the second had 50% of trials as null events, with a mean ITI of 3.92 sec (min 2 sec, max 18 sec), to better mimic what is recommended more often in the literature (Liu, 2004). For the slow event-onset condition of the simulations, binary event onset vectors were generated in line with a previous simulation study on modeling the HRF (Lindquist et al., 2009). Events were 1 sec in length, with an ITI of 30 sec. Since the response curve of the HRF is ∼16 sec long, “jittering” with null events was not necessary.

Generating BOLD time series

All onset vectors were then convolved with a canonical HRF (difference of two gamma functions) by using the SimTB toolbox (Erhardt et al., 2012) to arrive at vectors of anticipated BOLD response following task stimuli ( u _i). These convolved onset vectors were then used in data simulation. The data-generating equation is found via algebraic manipulation of Equation (2): $\begin{matrix} y_{i, t, k} = {(I - A_{i} - τ_{i, t} u_{i, t})}^{- 1} (Φ_{i}) y_{i, t - 1} \\ + {(I - A_{i} - (τ_{i}) u_{i, t})}^{- 1} (γ_{i}) u_{i, t} \\ + {(I - A_{i} - τ_{i} u_{i, t})}^{- 1} 1_{i, t} . \end{matrix}$ (3)

Here, the individual and group pattern connectivity matrices are combined into one (e.g., $A_{i} = A_{i}^{i} + A_{i}^{g}$ ) for succinct presentation. Data were simulated for 150 participants. For all conditions, the total number of time points T = 500. Time series were started with a burn-in of 1000 time points. The $1_{i}$ noise vectors were generated to have a variance of 2 to ensure an average signal-to-noise ratio (SNR) of 2.65. The SNR here was defined as the ratio of the standard deviation of the convolved onsets to the standard deviation of the noise in the data following the use of SNR in DCM studies (Welvaert and Rosseel, 2013). The pattern of connectivity in the data was simulated as represented in Figure 3, and it includes both modulatory (on the relationship between ROIs 3 and 4) and direct task effects. All code for simulations, as well as output, is available at: https://osf.io/d6wbx/.

FIG. 3.

Schematic representation of connectivity for generation of simulated data. “ROI” = brain region of interest; “Task” = person-specific HRF-convolved task vector; “by task” = ROI is multiplied by convolved task vector. All paths between other variables represent intrinsic connectivity. Solid lines indicate a contemporaneous connection, and dashed lines indicate a lagged connection of the first order. Intrinsic connections are seen among ROIs 1–3; the direct task effect exists for ROI2 and ROI4; the task modulates the relation between ROI3 and ROI4 such that no relation exists in the absence of the task but is present during the task. All variables have (first-order) autoregressive effects, indicated by looped dashed line, except the convolved task vector and the modulatory effect. The numbers above the lines indicate the connection weight used in simulating the data.

Simulation results

Simulation results (Table 1) showed that for the slow event-related design, there was a 100% recovery rate for all paths, regardless of the group cutoff value (0.75 vs. 0.51). Both rapid event-related designs, however, had an ∼80% recovery rate of all paths, for group cutoff values of 0.75 and 0.51. This was due to a very low (near 0%) recovery of the Task-ROI4 path, which was consistently recovered in the slow event onset condition. The Task-ROI4 path has a lower connection weight (0.2) than the other direct effect (0.6), as well as the bilinear connection (0.6), so it may be more difficult for HRF-GIMME to capture these weaker effects in the rapid event-related designs where there is less variability within an individual's HRF.

Table 1.

Recovery Rate of All True Paths As Well As Overall Sensitivity and Specificity Across Each Condition of the Simulation

Condition	ROI1→ROI2	ROI2→ROI3	Task→ROI2	ROI3_by_Task→ROI4	Task→ROI4	Sensitivity, %	Specificity, %
0.75 cutoff, slow design	1.00	1.00	1.00	1.00	1.00	100.00	91.46
0.51 cutoff, slow design	1.00	1.00	1.00	1.00	1.00	100.00	91.23
0.75 cutoff, rapid design 20% null	1.00	1.00	1.00	1.00	0.0007	80.01	99.97
0.51 cutoff, rapid design 20% null	1.00	1.00	1.00	1.00	0.0004	80.00	99.99
0.75 cutoff, rapid design 50% null	1.00	1.00	1.00	1.00	0.00	80.00	92.31
0.51 cutoff, rapid design 50% null	1.00	1.00	1.00	1.00	0.00	80.00	92.31

ROI, region of interest.

False positive rates were relatively low across all conditions. The slow event onset condition had a specificity of ∼92%, and there was a negligible (0.2%) increase in false positive rates between the 0.75 and 0.51 group cutoff values. This was similar to the rapid onset condition with 50% of null events, which had a 92.23% specificity at both a group cutoff value of 0.51 and a cutoff value of 0.75. The rapid event onset condition with only 20% null events had a lower rate of false positives, with 99.97% specificity in the 0.75 group cutoff condition and a 99.99% specificity in the 0.51 group cutoff condition.

Empirical Data Examples

To investigate the use of HRF-GIMME with empirical data, we obtained publicly available data from the OpenfMRI project (PI: Poldrack; NSF Grant: OCI-1131441). The first dataset tested (Accession No. ds101) was a rapid event-related Simon task, similar to the rapid event onset simulation condition with 20% null events, with a sample size of N = 21. The second dataset tested (Accession No. ds102) (Kelly et al., 2008) was a slow event-related Flanker task, similar to the slow event onset simulation condition, with N = 25. Both studies were approved by the institutional review boards of the data collection sites. Details regarding the designs and processing of data can be found in the Supplementary Data, along with additional analyses sub-sampling the data to ensure consistency of results due to smaller sample sizes consistency of results due to smaller sample sizes (Supplementary Figs. S1–S4).

Eight regions were selected separately for each task by using centers of activation (Tables 2 and 3) from a meta-analysis of multiple tasks, including the Flanker and Simon (Nee et al., 2007) as well as a meta-analysis of the frontal eye fields (Paus, 1996). A 6 mm sphere was used to construct the ROIs, and ROI time series were extracted. The ROI time series and binarized event onset vectors for the two runs were concatenated for each subject, and each dataset was separately analyzed by using HRF-GIMME. Modulating effects were then explored by multiplying the task vector with three ROIs determined to be likely related to the task from previous literature (Kerns, 2006; Mennes et al., 2010; Nee et al., 2007). In the Flanker task, this was the left precuneus (lPrec), medial frontal/anterior cingulate cortex (ACC), and the right dorsolateral prefrontal cortex (rdlPFC). In the Simon task, this was the thalamus, ACC, and right inferior parietal. All variables were mean centered before multiplying, and the data were standardized. Both datasets were analyzed with a group cutoff value of both 0.51 and 0.75.

Table 2.

Coordinates in Montreal Neurological Institute Space for Centers of Constructed Region of Interests for Flanker Task, Along with Their Corresponding Brodmann's Area, Anatomical Region, and Abbreviation Denoted in Output Figures

x	y	z	BA	Region
40	14	28	9	rdlPFC
36	16	20	13	rIns
0	20	48	6/8/32	Medial frontal/ACC
−36	16	4	13	lifIns
−18	−72	42	7	lPrec
2	16	46	8	mfACC
−32	−5	50	6	lFEF
31	−5	50	6	rFEF

ACC, anterior cingulate cortex; BA, Brodmann's area; lFEF, left frontal eye fields; lifIns, left inferior frontal/insula; lPrec, left precuneus; mfACC, medial frontal/anterior cingulate cortex; rdlPFC, right dorsolateral prefrontal cortex; rFEF, right frontal eye fields; rIns, right insula.

Table 3.

Coordinates in Montreal Neurological Institute Space for Centers of Constructed Region of Interests for Simon Task, Along with Their Corresponding Brodmann's Area, Anatomical Region, and Abbreviation Denoted in Output Figures

x	y	z	BA	Region
40	14	28	9	rdlPFC
−22	−64	46	7	lPrec
0	20	48	6/8/32	Medial frontal/ACC
50	−44	32	9	rIP
−42	16	28	9/6/46/8/13	ldlPFC
0	−26	6	50	Thal
−32	−5	50	6	lFEF
31	−5	50	6	rFEF

ldlPFC, left dorsolateral prefrontal cortex; rIP, right inferior parietal lobule; Thal, thalamus.

Figure 4 depicts the results from the Flanker design data. At a group cutoff value of 0.51 (Fig. 4, right panel), a group-level direct effect of the task on the lPrec was detected, as well as group-level connections between the left frontal eye fields (lFEF) and lPrec, lFEF and the ACC, the rdlPFC-right FEF, as well as the rdlPFC-right insula (rIns) and the medial-frontal ACC-ACC connections (the only two group-level paths detected at the 0.75 group cutoff [Fig. 4, left panel]). At the individual level, direct effects of the task on all ROIs and modulatory effects of the selected ROIs and the task were detected, and the total number of modulatory paths did not differ between the cutoff values.

FIG. 4.

HRF-GIMME summary paths plot for the slow event-related Flanker task, at group cutoff values of 0.75 (left) and 0.51 (right). Note that variables indicating “X_by_X” represent modulating effects, and the variable “task” indicates the direct effect of the task on the other variables. All paths between other variables represent intrinsic connectivity. Solid lines indicate a contemporaneous connection, and dashed lines indicate a lagged connection of the first order. Black lines indicate a group-level connection, and gray lines indicate an individual-level connection. The width of any given line corresponds to the count of individuals for whom the path was generated.

The results aligned with known circuitry. The FEF, insula, ACC, and dlPFC show strong activation during the Flanker task (Mennes et al., 2010; Nee et al., 2007), and the precuneus, dlPFC, and ACC are expected to activate together, as all have been shown to have the highest responses to certain Flanker task conditions (van Veen et al., 2001). Further, the consistent group-level rdlPFC-rIns connection may represent the linkage between the two distinct task-positive cognitive control networks proposed by Dosenbach and colleagues (2007), the frontoparietal network (which includes the FEF, dlPFC, and precuneus) and cinguloopercular network (which includes the insula and ACC). Thus, the presence of group-level connection between these regions and related to the task suggests valid detection by HRF-GIMME, particularly at the 0.51 group cutoff.

For the rapid event-related Simon task, results were identical for both the 0.75 and 0.51 group cutoff values (Fig. 5). Only one group-level path was detected, the lFEF-ACC. Far fewer task-related effects were detected than in the slow event-related design, where all tested ROI-task modulations had significant connections at the individual level. Hence similar to the simulation results, HRF-GIMME detected modulatory effects in the slow event-related design more so than the rapid event-related design, perhaps reflecting minimized variability in the HRF for rapid designs. Again, results from the rapid event-related Simon task align with known circuitry, particularly at the individual level. The dlPFC, ACC, and parietal cortex have all been shown to activate strongly during the Simon task (Nee et al., 2007), and once again, the precuneus, dlPFC, and ACC are expected to activate together (van Veen et al., 2001), and all selected ROIs are generally part of cognitive control networks (Dosenbach et al., 2007). Thus, direct effects of the task on nearly all ROIs indicate valid detection of connections by HRF-GIMME.

FIG. 5.

HRF-GIMME summary paths plot for the rapid event-related Simon task, at a group cutoff value of 0.75 (left) and 0.51 (right). Note that variables indicating “X_by_X” represent modulating effects, and the variable “task” indicates the direct effect of the task on the other variables. All paths between other variables represent intrinsic connectivity. Solid lines indicate a contemporaneous connection, and dashed lines indicate a lagged connection of the first order. Black lines indicate a group-level connection, and gray lines indicate an individual-level connection. The width of any given line corresponds to the count of individuals for whom the path was generated.

Although group-level paths may be less likely to appear in rapid event-related designs, it is important to underscore that a benefit of GIMME is its ability to create individual-level networks. As previously mentioned, HRF-GIMME is able to uniquely estimate each individual's HRF, as well as obtain individual-level path estimates. This allows the user to understand the heterogeneity that exists between individuals in terms of functional connectivity. Figure 6 highlights this by comparing the estimated HRF and resultant connectivity plots of two individuals from the Simon task. The estimated HRFs of the two individuals differ, with subject 8's (Fig. 6a) being quite canonical, and subject 13's (Fig. 6b) being comparatively delayed. The connectivity plots also differ, as both subjects have different direct effects of the task, and subject 8 has both positive and negative direct task effects, highlighting two potentially distinct circuitries.

FIG. 6.

From the rapid event-related Simon task data analyzed. (a) The estimated HRF (left) using the sFIR function of HRF-GIMME and the individual path connectivity plot (right) for subject 8. (b) The estimated HRF (left) using the sFIR function of HRF-GIMME and the individual path connectivity plot (right) for subject 13. For the connectivity plots, a red line indicates a positive connection, a blue line indicates a negative connection, a solid line indicates a contemporaneous connection, and a dashed line indicates a lagged connection of the first order. The weight of the line corresponds to the magnitude of the connection, such that a thicker line indicates a stronger magnitude of connectivity. Color images are available online.

Discussion

HRF-GIMME enables the modeling of task effects directly within a robust functional connectivity framework. A number of benefits occur through this novel integration, allowing for the HRF-convolved task vector to be included in the connectivity maps derived from the original GIMME algorithm. For one, both intrinsic and task-modulated connectivity patterns are obtained. In this way, researchers can identify which connections are impacted by task conditions and demands. Two, any influence the task has directly on specific brain regions will be accounted for. In this way, HRF-GIMME controls for the potential direct task effects when estimating and identifying the connectivity among regions. This also enables investigation into mediation effects, thus clarifying results seen in univariate spatial parametric tests by allowing for the possibility that other brain regions mediate the relationship between task onset and BOLD response. Three, the power to detect connectivity can be increased by altering the threshold of what is considered the “majority” in the search procedure, thus enabling greater detection even if the signal is weak across some individuals.

GIMME attends to individual differences in two ways: allowing for person-specific connections and estimating each connection weight separately for each individual. The algorithm allows for qualitative (i.e., presence or absence of connections) and quantitative (i.e., numeric) differences to surface in the values for connections. As demonstrated elsewhere, GIMME also has the capacity to subgroup individuals based on all paths (Gates et al., 2017), including intrinsic, extrinsic, and modulatory. Taken together, HRF-GIMME attends to known heterogeneity in brain processes while also assessing what is common in brain processes across the sample. Using both simulated and empirical data, the HRF-GIMME performed well in terms of its ability to recover expected relations between convolved task vectors and brain regions.

When applied to three sets of simulated data in which one condition emulated a slow event-related design and two utilized rapid event-related designs, the sensitivity, or ability to recover true connections, was greater in the slow onset condition (100%) compared with the rapid onset conditions (80%). In the rapid event-related conditions, the consistent finding was that the weaker of the two simulated direct effects of the task on BOLD values was not detected. That the modulatory effects were detected was expected, given that the modeling approach within used in HRF-GIMME has been found in prior studies to have more power to detect person-specific true effects in simulated data than DCM (Gates et al., 2011). Empirical fMRI data yielded similar results, with the slow event-related design detecting many more task-related effects overall than the rapid design (although it is possible that fewer connections consistently exist across individuals in rapid event-related designs). In the simulated data, varying the cutoff of the proportion of individuals that determines a group-level path from 0.75 to 0.51 did not affect the sensitivity of any of the results, although this might be due to ceiling effects. In the empirical data, varying the group cutoff value did not impact the results of the rapid event-related design, but it did yield a group-level direct effect of the task when decreasing the cutoff from 0.75 to 0.51 in the slow event-related design. Importantly, both sets of results conformed to expectations of connectivity patterns given prior literature, indicating that the results may be providing accurate inferences.

The specificity in the simulated data was high across both the slow event-related onset (∼91%) and rapid event-related onset (∼99.9% for 20% null, and 92% for 50% null) designs, meaning that very few false positives were recovered. Notably, varying the cutoff of the proportion of individuals that determines a group-level path did very little to affect the specificity of the results. Decreasing the group cutoff from 0.75 to 0.51 only decreased the specificity by 0.2% in the slow event-related condition, and it did not decrease the specificity in either rapid event-related condition. This suggests that relaxing the criterion for what constitutes a group-level path did not introduce false positives into the search procedure.

Considering these results together, due to the greater detection of task-related relations at a group cutoff of 0.51 in the empirical data, as well as the lack of increased false positives in this condition for the simulated data, we recommend that the user set the cutoff value for group-level paths to 0.51 for HRF-GIMME given the lower power to detect effects from modulating variables. Results from the slow event-related empirical data indicate that the high levels of heterogeneity between individuals may make detecting task-related effects, and particularly the nuanced modulating effects, difficult. Reducing the group cutoff value should improve the chances of group-level detection of task effects, with minimal risk for increasing false positives.

Researchers selecting the optimal analytic approach should consider design features and research questions when selecting a method. First and foremost, HRF-GIMME necessitates event-related designs. For designs where there is no time-varying task vector (e.g., resting state designs), original GIMME, vector autoregression, or partial correlation approaches would be more suitable as they do not require a time-varying task vector. Another design feature concerns the number of individuals and number of brain regions of interest. There is no upper limit to the number of individuals. The lowest number that has been tested has been 10 participants, in which case GIMME performed well (Gates and Molenaar, 2012). There is a limit to the number of brain regions to consider, with a maximum of 25 suggested when counting both the brain regions and task-by-brain region interactions tested. In terms of research questions, HRF-GIMME is ideal for applications where the researchers are interested in detecting the direct and mediating relations among brain regions as well as the modulatory and direct task effects. This is similar to DCM and contrasts PPI, for instance, where only the modulatory task effects are modeled. However, prior work suggests that the modeling approach underlying HRF-GIMME (euSEM) (Gates et al., 2011) has more power to detect modulatory effects than DCM, suggesting a benefit over this approach. Finally, HRF-GIMME is also well suited for data where some heterogeneity in the functional connectivity maps is expected across the individuals. Given that connectivity tends to vary across individuals (Laumman et al., 2015), the assumption of heterogeneity across individuals is likely required when conducting analyses.

Several limitations should be considered. As previously mentioned, although each individual has a uniquely estimated HRF, one drawback of HRF-GIMME is that each brain region within an individual is modeled using the same HRF, yet the HRF can vary between brain regions (Miezin et al., 2000). Although it is known that the variation between brain regions is less than the variation between individuals (Handwerker et al., 2004), enabling ROI-specific HRF estimates could easily be done within this analytic framework provided computationally efficient approaches are used. This represents an area for future work.

In addition, due to the low detection of task-related paths in both the simulated and empirical data, irrespective of group cutoff value, we caution that some weaker task-related effects may be missed when using HRF-GIMME to model rapid event-related designs. As shown in Figure 2, due to the nature of the HRF, rapid event-related designs often result in less variability in the expected HRF, which makes detection of the direct and modulating effects of the task difficult. The user may want to consider designing a slow event-related task for evaluating task-modulated connections. Future extensions of HRF-GIMME will aim at improving the power to detect task-related effects, particularly by including adaptations for rapid event-related designs.

Conclusion

HRF-GIMME is a novel and powerful new method for obtaining directed functional connectivity maps of event-related designs. Individual differences in connectivity patterns are obtained along with connections that are found to be consistent across individuals, allowing for generalizable inferences. Direct, mediating, and modulating influences of event-related tasks on brain connectivity can be captured by using HRF-GIMME. Importantly, person-specific estimation of the HRF is utilized within the modeling approach. HRF-GIMME is freely and publicly available as an option in the gimme R package (Lane et al., 2019). Exemplar simulated data and code are provided with the package as well as the OSF page associated with this article.

Footnotes

Authors' Contributions

All co-authors of this article contributed significantly and meaningfully to its creation, including writing and revision. The following details each author's particular contributions:

K.A.D.: Software, Data curation, Formal analysis, Visualization, Writing-Original Draft, Writing-Review and editing. Z.F.F.: Software, Formal analysis, Methodology, Writing-Original Draft, Writing-Review and editing. C.A.A.: Software, Methodology, Formal analysis, Writing-Original Draft, Writing-Review and editing. P.C.M.M.: Software, Methodology, Writing-Review and editing. J.H.: Conceptualization, Writing-Review and editing. J.R.C.: Conceptualization, Supervision, Writing-Review and editing. A.M.B.: Methodology, Conceptualization, Writing-Review and editing. M.A.L.: Methodology, Writing-Review and editing. M.N.H.: Conceptualization, Writing-Review and editing. K.M.G.: Software, Formal analysis, Methodology, Supervision, Project administration, Funding acquisition, Writing-Original Draft, Writing-Review and editing.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

This project was supported by the National Institutes of Health [Grant No. T32DA050560] and the National Institute of Biomedical Imaging and Bioengineering [Grant No. R01EB021299]. Zachary F. Fisher acknowledges the support of the Population Science Training Program (T32 HD007168) and the general support of the Carolina Population Center (P2C HD050924).

Supplementary Material

Supplementary Data

Supplementary Figure S1

Supplementary Figure S2

Supplementary Figure S3

Supplementary Figure S4

References

Beltz

, Moser

, Zhu

, et al. 2018. Using person-specific neural networks to characterize heterogeneity in eating disorders: illustrative links between emotional eating and ovarian hormones. Int J Eat Disord, 51:730–740.

Cohen

, D'Esposito

. 2016. The segregation and integration of distinct brain networks and their relationship to cognition. J Neurosci, 36:12083–12094.

Dosenbach

NUF

, Fair

, Miezin

, et al. 2007. Distinct brain networks for adaptive and stable task control in humans. Proc Natl Acad Sci U S A, 104:11073–11078.

Erhardt

, Allen

, Wei

, et al. 2012. SimTB, a simulation toolbox for fMRI data under a model of spatiotemporal separability. Neuroimage, 59:4160–4167.

Finn

, Shen

, Scheinost

, et al. 2015. Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. Nat Neurosci, 18:1664–1671.

Friston

, Buechel

, Fink

, et al. 1997. Psychophysiological and modulatory interactions in neuroimaging. Neuroimage, 6:218–229.

Friston

, Harrison

, Penny

. 2003. Dynamic causal modelling. Neuroimage, 19:1273–1302.

Gates

, Lane

, Varangis

, et al. 2017. Unsupervised classification during time-series model building. Multivariate Behav Res, 52:129–148.

Gates

, Molenaar

PCM

. 2012. Group search algorithm recovers effective connectivity maps for individuals in homogeneous and heterogeneous samples. Neuroimage, 63:310–319.

10.

Gates

, Molenaar

PCM

, Hillary

, et al. 2010. Automatic search for fMRI connectivity mapping: an alternative to Granger causality testing using formal equivalences among SEM path modeling, VAR, and unified SEM. Neuroimage, 50:1118–1125.

11.

Gates

, Molenaar

PCM

, Hillary

, et al. 2011. Extended unified SEM approach for modeling event-related fMRI data. Neuroimage, 54:1151–1158.

12.

Glover

GH.

1999. Deconvolution of impulse response in event-related bold fMRI. Neuroimage, 9:416–429.

13.

Gonzalez-Castillo

, Chen

, Nichols

, et al. 2017. Variance decomposition for single-subject task-based fMRI activity estimates across many sessions. Neuroimage, 154:206–218.

14.

Goutte

, Nielsen

, Hansen

. 2000. Modeling the hemodynamic response in fMRI using smooth FIR filters. IEEE Trans Med Imaging, 19:1188–1201.

15.

Greve

. 2019. Optseq2 [Computer software]. https://surfer.nmr.mgh.harvard.edu/optseq/ Last accessed June 8, 2019 .

16.

Handwerker

, Ollinger

, D'Esposito

. 2004. Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses. Neuroimage, 21:1639–1651.

17.

Kelly

AMC

, Uddin

, Biswal

, Castellanos

, Milham

. 2008. Competition between functional brain networks mediates behavioral variability. NeuroImage, 39:527–537.

18.

Kerns

JG.

2006. Anterior cingulate and prefrontal cortex activity in an fMRI study of trial-to-trial adjustments on the Simon task. Neuroimage, 33:399–405.

19.

Kim

, Zhu

, Chang

, et al. 2007. Unified structural equation modeling approach for the analysis of multisubject, multivariate functional MRI data. Hum Brain Mapp, 28:85–93.

20.

Lane

, Gates

. 2017. Automated selection of robust individual-level structural equation models for time series data. Struct Equ Modeling, 24:768–782.

21.

Lane

, Gates

, Pike

, et al. 2019. Uncovering general, shared, and unique temporal patterns in ambulatory assessment data. Psychol Methods, 24:54–69.

22.

Laumann

, Gordon

, Adeyemo

, et al. 2015. Functional system and areal organization of a highly sampled individual human brain. Neuron, 87:657–670.

23.

Lindquist

, Meng Loh

, Atlas

, et al. 2009. Modeling the hemodynamic response function in fMRI: efficiency, bias and mis-modeling. Neuroimage, 45:S187–S198.

24.

Lindquist

, Wager

. 2007. Validity and power in hemodynamic response modeling: a comparison study and a new approach. Hum Brain Mapp, 28:764–784.

25.

Liu

TT.

2004. Efficiency, power, and entropy in event-related fMRI with multiple trial types. Neuroimage, 21:401–413.

26.

Logothetis

NK.

2008. What we can do and what we cannot do with fMRI. Nature, 453:869–878.

27.

McClelland

, Judd

. 1993. Statistical difficulties of detecting interactions and moderator effects. Psychol Bull, 114:376–390.

28.

Mennes

, Kelly

, Colcombe

, et al. 2013. The extrinsic and intrinsic functional architectures of the human brain are not equivalent. Cereb Cortex, 23:223–229.

29.

Mennes

, Kelly

, Zuo

X.-N.

, et al. 2010. Inter-individual differences in resting-state functional connectivity predict task-induced BOLD activity. Neuroimage, 50:1690–1701.

30.

Miezin

, Maccotta

, Ollinger

, et al. 2000. Characterizing the hemodynamic response: effects of presentation rate, sampling procedure, and the possibility of ordering brain activity based on relative timing. Neuroimage, 11:735–759.

31.

Molenaar

PCM.

2004. A manifesto on psychology as idiographic science: bringing the person back into scientific psychology, this time forever. Measurement (Mahwah N J), 2:201–218.

32.

Moritz

, Bartz-Beielstein

. 2017. Imputets: time series missing value imputation in R. R J, 9:207.

33.

Mumford

, Ramsey

. 2014. Bayesian networks for fMRI: a primer. Neuroimage, 86:573–582.

34.

Nee

, Wager

, Jonides

. 2007. Interference resolution: insights from a meta-analysis of neuroimaging tasks. Cogn Affect Behav Neurosci, 7:1–17.

35.

Paus

1996. Location and function of the human frontal eye-field: a selective review. Neuropsychologia, 34:475–483.

36.

Price

, Beltz

, Woody

, et al. 2020. Neural connectivity subtypes predict discrete attentional-bias profiles among heterogeneous anxiety patients. Clin Psychol Sci, 8:491–505.

37.

Ramsey

, Spirtes

, Glymour

. 2011. On meta-analyses of imaging data and the mixture of records. Neuroimage, 57:323–330.

38.

Smith

, Miller

, Salimi-Khorshidi

, et al. 2011. Network modelling methods for FMRI. Neuroimage, 54:875–891.

39.

Sörbom

1989. Model modification. Psychometrika, 54:371–384.

40.

van Veen

, Cohen

, Botvinick

, et al. 2001. Anterior cingulate cortex, conflict monitoring, and levels of processing. Neuroimage, 14:1302–1308.

41.

Welvaert

, Rosseel

. 2013. On the definition of signal-to-noise ratio and contrast-to-noise ratio for fMRI data. PLoS One, 8:e77089.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.98 MB