1 INTRODUCTION

Computational modeling can deepen our understanding of how the brain processes information and produces overt behavior. In psychology, computational modeling has a long history of describing and explaining behavioral concepts. For example, reinforcement learning algorithms have been used to explain classical conditioning (Rescorla & Wagner, 1972), drift-diffusion models have been used to model reaction times in decision-making tasks (Ratcliff, 1978), and race models of reaction times have been used as theoretical formulations of visual–spatial attention (Bundesen, 1990). Similarly, different computational modeling approaches have been employed in neuroscience and neuroimaging. For example, generative graphical models of brain connectivity describing blood oxygenation level-dependent (BOLD) amplitudes in response to experimental inputs can be estimated using dynamic causal modeling (DCM) (Friston et al., 2017; Friston, Harrison, & Penny, 2003), and multivariate temporal response functions have been used to model ongoing sensory stimulation, like speech, in electrophysiological recordings (Crosse, Di Liberto, Bednar, & Lalor, 2016).

Although computational models are very prominent in the two fields, behavioral and neural responses are mostly treated separately (Turner, Forstmann, Love, Palmeri, & van Maanen, 2017). However, a combined modeling framework could provide deeper insights into the neural processes and the emergence of behavior. Different approaches have been proposed here: one possibility is to correlate the parameters of neural and behavioral models to describe how the different measures are related across different participants (Vossel, Weidner, Moos, & Fink, 2016). Alternatively, in model-based fMRI, the behavioral computational model's outputs (or hidden states) are used as a factor in a classical GLM analysis. One such factor could be a participant's perceived cue validity in a probabilistic spatial cueing task, recovered from reaction times (e.g., Dombert, Kuhns, Mengotti, Fink, & Vossel, 2016). Leveraging the theory-driven outputs of cognitive models allowed to determine more specific brain activation patterns of cognitive processes than by using nonspecific measures such as reaction times (Turner et al., 2017). A third option is a joint modeling approach (Turner et al., 2017). Here, an overarching set of parameters is used to describe both brain activity and behavior. An example is a study by Nunez (2015), where the drift-diffusion model parameters were constrained with task-based brain activity, incorporating the covariation between reaction times and neural activity on a trial-by-trial basis.

Although these approaches are tremendously useful, none of them employs an integrative model describing the generation of brain activity and behavior, which would allow us to directly investigate the hidden processes behind the two measurements. Rigoux and Daunizeau (2015) provided such a framework, where an additional output function extends DCM to describe behavioral responses (behavioral DCM, bDCM). This simultaneous modeling has high descriptive power and allows thorough diagnostics of the model. For example, by disabling specific nodes in the network (i.e., artificial lesions), conclusions can be drawn about the contribution or necessity of different brain regions to the emergence of behavioral patterns. So far—to our knowledge—bDCM has been applied to a larger dataset in one study only, which modeled binary choices in an economic decision-making task (Shaw et al., 2019).

The current study shows that bDCM can be extended to continuous measures (i.e., reaction times). Furthermore, we provide a direct comparison between bDCM and classical DCM, and between bDCM and an adjusted version of the Rescorla–Wagner model (Rescorla & Wagner, 1972; Vossel et al., 2014). We employ Bayesian model comparison based on the free energy of competing models and classical metrics of accuracy (mean absolute error and R2-score).

As a testing ground, we modeled the effects of attentional reorientation along the horizontal and vertical meridians in a spatial cueing paradigm, where participants had to report the orientation of a pre-cued Gabor patch. In trials in which invalid cues indicated an incorrect location of the target Gabor patch (20% of the trials), participants had to reorient their attention to the opposite location (Posner, 1980). This paradigm has been found to elicit reliable reaction time differences between invalid and valid trials, both on the individual and the group level (Hedge, Powell, & Sumner, 2017). Additionally, it has been shown that the internal representation of cue validity can be modeled using the Rescorla–Wagner model as a generative model of reaction times (Mengotti, Dombert, Fink, & Vossel, 2017; Vossel, Mathys, et al., 2014).

Besides the reliable behavioral effects, the cortical networks involved in this task have been characterized by multiple studies. We have previously analyzed the present dataset using classical DCM (Steinkamp, Vossel, Fink, & Weidner, 2020), which has also been used in similar cueing paradigms (c.f., Vossel, Weidner, Driver, Friston, & Fink, 2012). Moreover, studies in patients with stroke-induced lesions have revealed brain regions critically involved in spatial cueing tasks (Corbetta & Shulman, 2011; Malherbe et al., 2018; Posner, Walker, Friedrich, & Rafal, 1984). It is well established that the orientation of visual–spatial attention is mediated by a dorsal fronto-parietal attention network consisting of the intraparietal sulci (IPS) and the frontal eye fields (FEF). This network interacts with a ventral fronto-parietal attention network of ventral frontal cortex and the temporoparietal junction (TPJ) when a sudden reorientation of attention is necessary (Corbetta, Kincade, Lewis, Snyder, & Sapir, 2005; Corbetta & Shulman, 2011). In patients with spatial neglect, damage to ventral parietal regions such as TPJ causes a deficit in reorienting to contralesional targets. Moreover, it leads to dysfunctions in structurally intact dorsal regions such as the IPS (Corbetta et al., 2005), and direct lesions to the IPS have also been associated with impaired reorienting (Gillebert et al., 2011).

As IPS, FEF, and TPJ may differentially contribute to the behavioral outcome (RT), we used Bayesian model comparison to identify which regions convey information about the behavioral dynamics after accounting for the complexity of the network model.

2 METHODS 2.1 Participants

Data were collected from 29 participants (15 female, 21–39 years old, M = 25, SD = 3) with normal or corrected-to-normal vision [all right-handed, Edinburgh handedness Inventory (Oldfield, 1971), M = 0.86, SD = 0.14], who provided written informed consent to participate in the study. Participants had to be older than 18 and younger than 40 years old and had to be right-handed. Participants with neurological or psychiatric disorders were excluded from the study. Due to the fMRI protocol, we also excluded participants with metal implants and tattoos. One participant had to be excluded subsequently because of noncompliance. Another participant was excluded due to excessive head movement (predefined criteria translation >3 mm, rotation >3°). Furthermore, we could not extract the time series for the left-TPJ volume of interes (VOI) in one participant. Therefore, the final sample included 26 participants. The ethics board of the German Psychological Association had approved the study. Volunteers were paid 15€ per hour for their participation. The dataset has been used in a previous study (see Steinkamp et al., 2020).

2.2 Task

Participants performed a spatial cueing task while lying in a 3 T Trio (Siemens, Erlangen) MRI scanner. Stimuli were displayed on a screen behind the scanner bore, which could be seen via a mirror (mirror to display distance: 245 cm) mounted on a 32-channel head coil. The participants' task was to report the orientation (horizontal/vertical) of a target Gabor patch (size 1° visual angle) by button presses of either the left or the right index finger while continuously fixating a diamond in the center of the screen (0.5° visual angle). A brightening of the central diamond (500 ms) indicated the beginning of a trial and was followed by a spatial cue after 1,000 ms (brightening of one of the diamond's edges for 200 ms) that indicated the location of the following target stimulus with 80% probability. Participants were explicitly informed about the percentage of cue validity. The possible target locations were indicated by empty boxes (1° width) located to the fixation diamond's left, right, top, and bottom (4° visual angle). After 400 or 600 ms, the target stimulus appeared for 250 ms at the cued location or in the box opposite to it. Distractor stimuli (constructed from two overlapping Gabor patches that were rotated by −45° and 45°, respectively) appeared simultaneously in the remaining three locations. Participants performed two runs of the spatial cueing paradigm. In one run, targets and cues occurred along the vertical axis, in another along the horizontal axis (see Figure 1).

Illustration of the spatial cueing paradigm. In the upper row, a valid trial of the horizontal run is shown. The lower row depicts an example of an invalid trial in the vertical run. Reused from Steinkamp et al. (2020), licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/)

Each run consisted of five blocks of 40 trials (32 valid, 8 invalid). All possible combinations of target location, target orientation, and interstimulus interval were presented in random order within each block. The time between the trials was drawn from the set of 2.0, 2.7, 3.2, 3.9, or 4.5 s with equal probability. Between the blocks, there was a break of 10–13 s.

Run order (vertical or horizontal first) and the response mapping (left index finger for vertical orientations/right index finger for horizontal orientations or vice versa) were counterbalanced across participants. Before the experiment, participants performed a rapid detection task to train the mapping of stimulus–response associations. Here, targets appeared rapidly in the middle of the screen, and participants had to press the corresponding button as fast as possible. Immediate feedback and a running score of their accuracy were given. Additionally, there were 20 practice trials with feedback before each run of the main experiment.

Stimulus presentation and response collection were controlled using PsychoPy (version 1.85.3, Peirce, 2007, 2008; Peirce et al., 2019).

2.3 Behavioral analysis

The mean reaction times were calculated for each participant, cueing condition, and target location. Before calculating the mean reaction times, we preprocessed the data for each participant separately. First, incorrect, missed, and outlier trials were removed. Outliers were defined as trials with reaction times below 0.2 s or greater than the 75th percentile + 3 × interquartile range (IQR). The higher threshold for outlier exclusion was chosen to retain as many trials as possible in the analysis (removed trials, including errors, in the horizontal run: invalid M = 2.54, SD = 2.63; valid M = 6.62, SD = 5.91; in the vertical run: invalid M = 3.12, SD = 1.8; valid M = 6.0, SD = 3.94).

For the analysis of the “validity effect” (i.e., the slowing of reaction times in invalid as compared with valid trials), the data were pooled across the two runs (horizontal/vertical). The mean reaction times of the 2 × 4 (cueing × target location) factorial design were then analyzed in a repeated-measures ANOVA. The analysis was conducted in Python 3.7 using pingouin (version 0.3.3; Vallat, 2018).

2.4 fMRI analyses

For each participant and each run, we collected 557 T2*-weighted images using an echo-planar imaging (EPI) sequence [time of repetition (TR) 2.2 s; echo time (TE) 30 ms; flip angle 90°]. Each recorded volume consisted of 36 transverse slices with a slice thickness of 3 mm and a field of view of 200 mm. The voxel size was 3.1 × 3.1 × 3.3 mm. The first five images were discarded to account for T1 equilibrium artifacts. Next to functional images, we also obtained an anatomical T1-weighted image for each participant, which was used in the preprocessing.

We preprocessed the fMRI data using fmriprep (version 1.1.1; Esteban et al., 2019), a robust and standardized pipeline, which applies slice-time correction, realignment, and normalization to MNI space. A detailed preprocessing report can be created automatically (see http://fmriprep.readthedocs.io/en/1.1.1/workflows.html) and is included in the Supporting Information.

Data was further spatially smoothed using an 8 × 8 × 8 mm FWHM Gaussian kernel. This step was done in Matlab 2018b (The MathWorks, Inc., Natick, Massachusetts), using SPM12 (version 7.771; Friston, 2007).

2.5 fMRI—GLM

A classical GLM analysis was performed to identify activation peaks during attentional orientation and reorientation, later used to extract BOLD time-series data for the DCM analysis. The GLM analysis was conducted using SPM12. First-level models were created with four regressors of interest for each run, representing invalidly cued targets on the left (iL) and on the right (iR), validly cued targets on the left (vL) and the right (vR) for the horizontal run, and invalidly and validly cued targets in the lower (iD, vD) and the upper (iU, vU) part of the screen in the vertical run.

To account for other physiological noise in the BOLD signal, we added the three rotation and three translation estimates of the rigid body transform, the average white matter signal, and the average cerebral spinal fluid (CSF) signal as nuisance regressors. We further included the squared time-series of the eight regressors, the time-shifted time-series (t − 1), as well as the square of the shifted time-series, resulting in a total of 32 nuisance regressors (Friston, Williams, Howard, Frackowiak, & Turner, 1996). We also applied a high pass filter at 128 s. For each run, four first-level contrasts were calculated: T-contrasts of valid and invalid trials versus baseline, an F-contrast of target onset versus baseline, which was used in the VOI analysis, and a differential contrast of invalid trials greater than valid trials. The latter contrast isolates brain regions involved in attentional reorientation.

At the group (second)-level, we investigated the differential contrast of invalid > valid trials using two planned one sample permutation t-tests against 0 using SnPM 13 (Nichols & Holmes, 2002), with default settings, 10,000 permutations, and no additional variance smoothing, using the initial set of 27 participants. The cluster forming threshold was estimated during the processes with a predefined voxel-level cutoff of p < .001.

2.6 Modeling analysis

In the following, we will describe the modeling approaches used in our analysis, followed by a description of our model assessments and further analyses.

2.6.1 Rescorla–Wagner model

We employed a variant of the Rescorla–Wagner model used previously (Mengotti et al., 2017). While this study was interested in the parameter (the learning rate that describes how quickly participants adjust their internal assessment of the cue-validity), we applied this modeling approach to simulate reaction times in a trial-by-trial fashion. For parameter estimation, we defined new functions for the variational Bayesian analysis (VBA) toolbox (clone from master, in January 2020, Daunizeau, Adam, & Rigoux, 2014).

We used the following reinforcement learning formula as the evolution function, describing the hidden process governing the generation of reaction times: where describes the prediction error at trial . The external input describes whether the cue at time was either valid (0) or invalid (1), is the learning rate, and is the participant's perceived cue invalidity (i.e., the probability, that the cue will be invalid) after observation of trial . The observation function (i.e., the mapping from perceived cue invalidity to reaction times) was defined as: According to this formulation, the perceived cue invalidity of the previous trial governs the responses, with different bias parameters for valid and invalid trials and a general scaling parameter of the predictions.

To keep the behavioral dynamics as close as possible to the observed data, we set the reaction time of missing and outlier trials to 0 but ignored these trials during model inversion. The mean and standard deviation over participants of the posterior estimates can be found in the Section S6 in Supporting Information.

Table 1 depicts the Gaussian priors used in our estimation.

TABLE 1. Overview of parameters and prior values of the Rescorla–Wagner model Parameters 0.5 0.5 To ensure , was logit and inverse logit transformed during parameter updating 0 1 0 1 0 1 0.5 1 Initial state of 2.6.2 Behavioral DCM In the following, we will provide a short overview of key concepts of DCM. For a full derivation and detailed description of DCM (see, Friston et al., 2003; Rigoux & Daunizeau, 2015; Stephan et al., 2008). DCM is a full Bayesian approach to create a generative model of brain dynamics and infer effective connectivity between selected brain regions. In principle, DCM describes how experimental variations (described by the input ) drive the neural activity (, the hidden states) in brain regions of interest in a dynamical system. The evolution function describes the temporal dynamics of the hidden states () and how they are influenced by external inputs (). In DCM for fMRI, the evolution function is typically described as: corresponds to the number of inputs. The neural evolution parameters in correspond to the entries in (fixed connectivity between brain regions), (modulation of connection strength by input ), and (direct effects of inputs). Hemodynamic states z (dependent on the neural states x) are then gated through an observation function: This function captures BOLD signal variations based on the hemodynamic states () and the hidden neural activity (), with hemodynamic parameters . This mapping allows to observe and infer the hidden neural dynamics via the BOLD signal. BDCM augments the described formulation of DCM by adding new hidden states () for observed responses. The dynamic of is defined by , with as an additional evolution function, and as the observation function to map the hidden neural dynamics to behavioral responses. The evolution function of the new “behavioral” state follows the same rationale as the function in the DCM formulation: Here the parameter vector describes the linear () components of the behavioral state and the direct () and modulatory () influences of experimental manipulations. is an analogy of the weight vector in a regression model. In the original article, the neural states were mapped to binary behavioral observations (button press absent or present) via a sigmoidal function: Here, is an unknown bias term, and is the response or decision state. In our study, we slightly adjusted the sigmoid mapping by changing the scale on which it operates. As we are not expecting reaction times slower than 3 s, we used this as an upper bound: Regions

As in our previous study (Steinkamp et al., 2020), we included bilateral IPS and FEF in our DCM model, which correspond to the central nodes of the dorsal fronto-parietal attention network (Vossel, Geng, & Fink, 2014). Additionally, as part of the ventral attention network, we included the TPJ bilaterally. As additional inclusions (e.g., the inferior/middle frontal gyrus) would have increased model complexity and computational resources (and time), we did not include other brain regions, which may also play a role in attention reorienting.

Based on our assumptions about the dorsal and ventral attention network's interplay, we created three automatic meta-analyses using NeuroSynth (https://www.neurosynth.org/; Yarkoni, Poldrack, Nichols, van Essen, & Wager, 2011) to define the seed coordinates for the subsequent VOI analysis (see Table 2). Our regions of interest were bilateral IPS (search term: “intraparietal sulcus”), bilateral FEF (search term: “frontal eye”), and bilateral TPJ (search term: “tpj”). We downloaded the corresponding association maps (associations, p < .01 FDR corrected) and identified the seed location as the peak voxel in the cluster of interest, using the Anatomy toolbox (v2, Eickhoff et al., 2005). In all three maps, the two largest clusters encompassed our regions of interest in either the left or right hemispheres.

TABLE 2. Regions and search-terms for automated NeuroSynth meta-analyses Region NeuroSynth (accessed 10.10.19) Z-statistic X Y Z IPS—Left “Intraparietal sulcus” 14.6 −30 −50 42 IPS—Right “Intraparietal sulcus” 13.5 40 −38 44 FEF—Left “Frontal eye” 13.9 −30 −4 52 FEF—Right “Frontal eye” 14.6 32 −6 52 TPJ—Left “Tpj” 8.56 −60 −54 20 TPJ—Right “Tpj” 11.4 58 −50 14

We used the participant level t-maps (thresholded at p < .1 uncorrected) in each run to search for individual local maxima in a 12 mm sphere around the seed coordinates. The first principal component of BOLD time-courses in a 9 mm VOI around the participant's maximum was extracted and adjusted based on the F-contrast defined in the first-level analysis. Task-related activity for the IPS and FEF VOIs was defined by the contrast of valid trials against baseline and for TPJ by the contrast of invalid trials against the baseline.

Preprocessing

We preprocessed the BOLD signal by detrending each VOI signal (spm_detrend) and scaling the BOLD amplitude across VOIs to a maximum value of 4 (see spm_dcm_estimate). Behavioral data were extracted from the event data, and as in the previous analyses, error trials, trials with missed responses, and RTs fulfilling the outlier criterion (RT < 0.2 s and RT > 3 × IQR + UQ) were excluded.

BOLD data were resampled from a TR of 2.2 s to a sampling rate of 1.1 s (by interspersing “NaN” values). The behavioral observations were set to occur at the corresponding target onset, which was also downsampled to a resolution of 1.1 s. No resampling of BOLD data was performed for the classical DCM analysis. As the Rescorla–Wagner model represents trial-by-trial dynamics, the corresponding preprocessed reaction times were used, excluding error, missed, and outlier trials.

For our modeling, we assumed homogenous HRF dynamics across the six regions, fixing the initial states of the model to 0, and estimating the shape of the observation noise hyper-prior distributions. For this, we assumed that we would be able to explain 10–90% of the variance in both the BOLD and the reaction time data. The prior distributions over the other parameters were set to the defaults of the VBA toolbox. We used the same hyperpriors for the explained variance of the BOLD signal in the classical DCM analysis and the Rescorla–Wagner model's behavioral responses.

To define the inputs into the DCM models, we created separate SPM-design matrices that were only used to define the input streams. Stream one (u1) was defined as the driving input to all six regions, containing an impulse every time a target stimulus appeared (irrespective of the cueing condition or target location). The second stream (u2) was used purely for the modulatory effects, containing an impulse only in invalidly cued targets. The input streams were extracted from the SPM design matrix and were centered before entering the model inversion (spm_detrend). As mentioned above, the Rescorla–Wagner model is modeling trial-by-trial variations (rather than continuous time), so the input to this model was a vector consisting of ones and zeros, indicating whether the current trial is invalid or valid.

Model definition and comparison Before bDCM was compared against models of single modalities (classical DCM and Rescorla–Wagner model, respectively), we conducted Bayesian model selection to identify the most plausible configuration of output connections (i.e., the region in which neural activity is linked to the behavioral output). We used IPS, FEF, and TPJ as our brain regions of interest as described above. The fixed connectivity structures of our model (i.e., the A-matrix) had complete connections in each hemisphere and between homologous regions (Figure 2). As we did not include visual areas in our modeling approach, all six regions received driving input (C-matrix). For bidirectional intrahemispheric and interhemispheric modulatory connections (B-matrix), we considered the IPS and TPJ. Connections in both hemispheres to the FEF were unidirectional, assuming that there were no feedback modulations from FEF to the other brain regions. We then investigated how neural dynamics in the included brain regions are related to behavioral dynamics. More specifically, we tested the following alternative hypotheses regarding bDCM's matrix: Neural activity in IPS drives behavior. IPS is a major hub region in the dorsal attention network. It is thought to initiate top-down modulation of visual areas when attention is oriented in space. Moreover, it is often co-activated together with ventral frontoparietal regions during attentional reorienting in invalid trials, and isolated IPS lesions in stroke patients can lead to reorienting impairments (Gillebert et al., 2011; Vossel, Geng, & Fink, 2014). Neural activity in FEF drives behavior. FEF is also part of the dorsal attention system and is crucially involved in covert and overt attentional orienting (Corbetta, Kincade, & Shulman, 2002; Rizzolatti, Riggio, Dascola, & Umiltá, 1987). Neural activity in TPJ drives behavior. TPJ is the ventral region in our network model and has critically been related to detecting unattended behaviorally relevant events such as invalidly cued targets and mismatches between predicted and observed inputs (Corbetta, Patel, & Shulman, 2008; Mengotti, Käsbauer, Fink, & Vossel, 2020). All three regions drive behavior. This model was included to test whether a mix of all brain regions most plausibly describes behavior, despite the additional complexity. Bayesian model comparison was conducted with random effects (RFX) in the VBA-toolbox (VBA_groupBMCbtwConds) and fixed effects (FFX) to select the most likely output region (IPS, TPJ, FEF, or all).

Top row, basic structure of the DCM and bDCM models. Regions were fully interconnected in each hemisphere, and homologous regions were connected. All regions received driving input. We assumed that all regions' connections were modulated by invalid trials, except for feedback and interhemispheric connections from FEF. In the lower part, we conducted Bayesian model selection to select the most plausible output region(s), indicated by green arrows

After selecting the most plausible output region, we inverted a classical DCM model by setting the prior mean and variance of the parameter set to 0, essentially disabling the additional parameters necessary to fit behavioral dynamics. This included both parameters in the observation and evolution function. Moreover, we inverted models that included predictions from the behavioral Rescorla–Wagner model.

The following models were compared: bDCM (with output region selected as described above) DCM: The bDCM model above, but with the prior mean and variance of the evolution (i.e., ) and observation parameters set to 0, so that they are not considered in the model inversion. This model was chosen to compare the single modality model of fMRI and test whether bDCM merits its additional complexity. bDCM + Rescorla–Wagner (bDCM + RW): In this model, we added an additional input stream (), including the RT predictions derived by the Rescorla–Wagner model. The input-stream, however, was not included in the DCM part of the bDCM (i.e., A, B, C, and D) but was directly gated to the output function via the matrix. This model was included to assess if bDCM's predictions have additional value compared with the Rescorla–Wagner model's predictions. bDCM infused with the Rescorla–Wagner model (bDCM × RW): In this model, we replaced the input coding for invalid trials () with the prediction of the cue invalidity () of the Rescorla–Wagner model—note that this input was also centered. This tested whether the cognitive processes modeled by the Rescorla–Wagner model provide information over and above bDCM.